mrc-ide / PhyDyn

PhyDyn: Epidemiological modelling in BEAST
GNU General Public License v3.0
22 stars 7 forks source link

Setting initial conditions as functions of model parameters #14

Open francescopinotti92 opened 3 weeks ago

francescopinotti92 commented 3 weeks ago

Hello,

is there any way to set initial values of deme and nondeme compartments depending on model parameters? Let me consider for example a SEIR model at endemic equilibrium; I would like $S(t=0) = 1/R_0$ where $R_0$ is the basic reproductive number and is a function of transmission and recovery rates.

Thanks in advance!

emvolz commented 3 weeks ago

As an alternative to starting the system in equilibrium, you could start the system with a long period of time before sampling so that it has time to equilibrate in the period you are interested in. Starting the simulation with an initially small number of infections will also insure that lineages coalesce before reaching time zero; in practice this gives more sensible results, since it reverts to an unstructured constant Ne coalescent of lineages are not yet coalesced at time zero.

If you really want to proceed with S(0)=1/R0, I think you can achieve something like this, but I would set S to be on a natural scale, not proportional scale, because the ratio of lineages to size of each deme is important. Firstly I would define a parameter for the initial condition of S & population size N. Note that N may be different (generally smaller) than census population size, unless your model is accounting for all important factors that can change the ratio of Ne to infections; it is easier to estimate an "effective" number of infections that is proportional to I, and in this case you should estimate N. Secondly, define a variable for R0 in terms of N and S(0), and you can reference R0 in your other model equations. Here is an example of a definition: https://github.com/mrc-ide/PhyDyn/blob/fc81194b13a420714cca9bc244067623b9924bbf/examples/ebola/ebolaSEIRhet.xml#L432
I think your variable for initial S would need to appear in both the "rates" and "trajparams" blocks of the model spec, but I'm not 100% sure this would work. As always, for every parameter that is estimated, you should make sure there is an initial state, prior, operator, and logger.

francescopinotti92 commented 3 weeks ago

Thank you so much for the detailed answer, that's really helpful. I realised only later that I could start simulations long before sampling. Thanks also for the additional insights.