Influenza Example (original) (raw)
Seasonal human Influenza: An extended SIR model
In this example we model a single season of Influenza A virus (IAV) H3N2 and apply it to 102 HA-1 sequences collected between 2004 and 2005 in New York state. Our model builds on a simple susceptible-infected-recovered (SIR) model which accounts for importations of lineages from a global reservoir of IAV.
A typical SIR model has the following form:

where I(t) is the number of infected and infectious hosts, β is the per-capita transmission rate, S(t) is the number of hosts susceptible to infection, R(t) the number of hosts that have been infected and are now immune to this particular seasonal variant and N(t) denotes the total population size.
Our extended SIR model has two demes: I, corresponding to the IAV lineages circulating in New York, and _N_r, denoting the effective population size of the global IAV reservoir. The global reservoir is modelled as a constant-size coalescent process whereas transmission (and death) in I follow the standard SIR model: new infections within New York state are generated at the rate βI(t)S(t)/N and deaths from the pool of infected are proportional to γ.
Births within the global reservoir do not vary through time and depend on the effective population size in that deme. Births balance deaths in the reservoir population, both with rates γ _N(t)_r. Migration between the reservoir and New York is symmetric, with rate η I(t).
The resulting set of matrix equations for the model are:

The reproduction number _R_0 for this model is defined by β / γ . We introduce _R_0 and re-write the equations accordingly. The final equations look as follows:

Note that the introduction of the global reservoir do not effect the dynamics of I(t) over time: migrations between the reservoir and New York are balanced.
The population model described above can be described as a PopModelODE PhyDyn/BEAST2 object:
N = S + I + R foi = if (t > 2004.75) then R0gammaS/N else gamma foi * I gammaI if (t > 2004.75) then -foiI else 0.0 if (t > 2004.75) then gamma * I else 0.0 gamma * src gamma * src eta * I eta * I
Note that we have adapted our formulae to reflect the 2014-2015 flu season dates, beginning on September 2014 (t=2014.75). We assume a constant size coalescent process before the 2014-15 flu season, with the extended SIR model kicking off on t = 2014.75.
The PhyDyn population model is completed with the following:
We sample the migration rate eta, recovery rate gamma and reproduction number R0 , as well as the initial values for the number of infected and susceptible, and the total population size of the global reservoir srcSize. Note that an informative prior is used for the recovery rate, which increases identifiability of other parameters. The priors used are:
- Migration rate η (event per year): lognormal (log mean=1.38, log sd=1).
- Recovery rate γ (event per year): lognormal (log mean=4.8, log sd=0.25).
- Reproduction number _R_0: lognormal (log mean=0, log sd=1).
- Reservoir size _N_r: lognormal (log mean=9.2, log sd=1).
- Initial number of infected in September 2004: lognormal (log mean=0, log sd=1).
- Initial number of susceptible in September 2004: lognormal (log mean=9.2, log sd=1).
A complete BEAST2 xml implementing this example can be found at PhyDyn/examples/influenza/nyFluSeason.xml .