A Simple Stochastic Model with Environmental Transmission Explains Multi-Year Periodicity in Outbreaks of Avian Flu

Avian influenza virus reveals persistent and recurrent outbreaks in North American wild waterfowl, and exhibits major outbreaks at 2–8 years intervals in duck populations. The standard susceptible-infected- recovered (SIR) framework, which includes seasonal migration and reproduction, but lacks environmental transmission, is unable to reproduce the multi-periodic patterns of avian influenza epidemics. In this paper, we argue that a fully stochastic theory based on environmental transmission provides a simple, plausible explanation for the phenomenon of multi-year periodic outbreaks of avian flu. Our theory predicts complex fluctuations with a dominant period of 2 to 8 years which essentially depends on the intensity of environmental transmission. A wavelet analysis of the observed data supports this prediction. Furthermore, using master equations and van Kampen system-size expansion techniques, we provide an analytical expression for the spectrum of stochastic fluctuations, revealing how the outbreak period varies with the environmental transmission.

So far there is no approximation made in the derivation of Eqs. (A2)-(A4), however, we can now substitute those transition rates of Eqs. (1)-(3) into these equations and take the mean-field limit N, N V → ∞, which allows us to take the replacement si = s i . After introduce the fractional variables we can thus write down the deterministic equations as where let constant κ = lim

A1.2 The system-size expansion and analysis of the fluctuations
Now we are in a position to implement the application of the van Kampen [25] system-size expansion to study the fluctuations. First, we introduce new continuous variables x 1 , x 2 , x 3 , in place of the previously used discrete variables s, i, v and make following replacement in the transition probabilities which appear in the master equation Next, defining a new probability distribution function Π by P (s, i, v; t) = Π(x 1 , x 2 , x 3 ; t), which implies that When using this formalism it is useful to rewrite the master equation (A1) in another way that is by introducing the step operators [21, 26] the Eq. (A1) with transition rates of Eq. (1)-(3) can be rewritten as Expanding the step operators E ±1 s and E ±1 i in a power series in and then substituting these expanding equations into Eq. (A8), drawing a comparison of which with Eq. (A7) order by order yields the so-called macroscopic equations to leading order, where which are indeed the equations (A6). Therefore, an alternative and systematic way of obtaining the mean-field equations is just the leading order terms in the system-size expansion. The next-to-leading order gives rise to a Fokker-Planck equation for the fluctuation variables Since we are interested in fluctuations about the fixed point E * = (φ * 1 , φ * 2 , ψ * ) defined in Eq. (6) of the determinist model, both matrix A = (A kl ) 3×3 and B = (B kl ) 3×3 are evaluated at this fixed point, whose explicit forms are found to be and with B 11 = βφ 1 φ 2 + ρφ 1 ψ + µ(1 − φ 1 ),

A1.3 Power spectral calculation and its peak
To calculate the power spectra of the fluctuations around the stationary state, we have to make a Fourier analysis, so it is first essential to formulate a set of Langevin equations of the stochastic variables x k (t), (k = 1, 2, 3). The Langevin equations corresponding to Eq. (A11) are which are three differential equations describing the stochastic behavior of the model at large, but finite N . The variables x l (l = 1, 2, 3) are stochastic corrections to the deterministic variables s, i and v, and ξ k (t)(k = 1, 2, 3) are Gaussian white noises with zero mean and a correlation function given by Taking the temporal Fourier transform of (A14) gives Actually, this Fourier transform is a system with three coupled linear algebraic equations which can be used to obtain a closed form expression for the power spectra. Therefore, solving Eq. (A15), we obtaiñ where and the denominator D is given by Averaging the squared moduli ofx k , (k = 1, 2, 3) gives the power-spectra of variables S, I and V : where By using these methods, we can analytically predict the epidemic outbreaks and fade-outs on a certain disease, as done by Alonso etal [21] for several childhood diseases.
A2 Stability analysis of fixed points E 0 and E * of system (A6) First, we give the stability analysis on equilibria E 0 , the Jacobian matrix of which given by and its characteristic polynomial can be written as Hence, it is easy know that the trivial fixed point E 0 is stable if and only if the conditions (η − δ) + (γ + µ − β) > 0 and (η − δ)(γ + µ − β) − κρτ > 0 are satisfied, which equal to the condition of basic reproduction number, i.e., R 0 < 1, and R 0 = β µ+γ + κρτ (η−δ)(µ+γ) . Second, the stability analysis on equilibria E * . The Jacobian matrix at E * is where where the stars indicate the values at equilibrium E * , which are φ * 1 = 1 R0 , φ * 2 = µ µ+γ (1 − 1 R0 ), and ψ * = κµτ (η−δ)(µ+γ) (1 − 1 R0 ). The characteristic polynomial of Jacobian matrix J 2 is written as a cubic polynomial about λ denoted as L(λ): all of the coefficients of which are both of the above inequations (A19) and (A20) are obtained in the condition: Then we can obtain According to the criterion of Routh-Hurwitz, the conditions a > 0, ab − c > 0, and c > 0 are satisfied, thus the equilibrium E * is stable.

A3 Environmental transmission
The environmental transmission rate is defined as the per capita rate at which susceptible individuals get infected through the environmental transmission route. This rate will be higher the higher the concentration of virion particles in the environment. This is so because individual birds are continuously exposed to the disease via ingestion of lake water while feeding. In general, whenever dose-response curves have been empirically measured, it has been observed that the probability of infection indeed depends on the "dose", or quantity of infection agents ingested by individuals. From these curves, a ID 50 is defined as the virion concentration at which half of the tested population gets infected.
In order to approximate how the environmental transmission rate depends on virion concentration in the water we make the following argument.
If N birds are exposed to contaminated lake water with a virion concentration V in the water, and monitored for a time t, the number, n, of birds that will develop the infection will be distributed following a binomial distribution: where the probability p can be written as: and α is the volume of water ingested per unit time per individual bird. Notice that p is a number (total ingested virion particles), lacking physical dimensions, as corresponds to a probability. However, this model only works for times short enough to make p < 1. The limit of large bird populations (N ), and low total virion ingestion (αV t), leads to the familiar Poison distribution, under which the probability of observing n individual infections during a finite time t is: Notice that the probability of observing at least one infection in a time interval t is 1 − P (0). This can be seen as a the pre-factor reducing a faster asymptotic per capita rate, ρ a , at which susceptible individuals would get infected through the environmental transmission route if the virion concentration in the water were very big, effectively infinite. So, we can write: where we have defined α t ≡ α t to make contact with the expression used by Breban et al's (2009).
can be regarded at the probability distribution function describing the time until an individual gets infected when ingesting water at rate α, contaminated with viruses at a concentration V . Therefore, this waiting time is, in fact, a random variable governed by the exponential distribution. As a consequence, the instantaneous ingestion rate α can be related to the ID 50 , as defined above: where ID 50 is virion concentration at 50% of infection and the time T is the exposure time dose-response tests take. This equation allows to approximate α through typical dose-response experiments. By exposing individuals to higher and higher viral concentrations, the same doseresponse experiments can also inform about the value of the asymptotic value ρ a .
Furthermore, Eq. (A24) can be approximated by: as long as α V t is small, in accordance with our initial assumption. In addition, what is clear from Eq. (A26) is that, within the range of this approximation, the environmental transmission rate,ρ(V ), can be assumed to be proportional to the concentration of virion particles in the the water, V . In fact, throughout our paper, we have used this proportionality, expressed as: where N V is a scaling factor for the virion concentration in the environment. As you see from Eq. (A26), the proportionality factor is, in the context of our approximation, the product (α t ρ a ).
The proportionality factor in Eq. (A27) is actually one of our model parameters, ρ (see our Eq (5) in the main text). The following argument and data from the literature allow us to obtain an estimation of some of our model parameters.
The volume of water ingested per unit time per individual bird has been estimated to be of about 2500 to 25000 liters/year (Bennett and Hughes, 2003). Experimental studies show that ID 50 take values of 10 1.8 to 10 5.7 viral particles per ml (Lu et al 2004 andIto et al 1995). With this data in mind, and using the relation between α and ID 50 above, we can write: On one hand, we can identify N V ≈ ID 50 , and, on the other, t and T are both short times that can be considered of the same order, t/T ≈ 1. So we can finally write the environmental transmission rate as:ρ This shows that the model parameter ρ corresponds, in the context of our approximation, to the asymptotic value ρ a times log e (2). Moreover, this argument supports our choice of roughly N V = 10 5 throughout our paper (see Table 2), and points to ways of estimating our model parameter ρ by measuring the asymptotic rate at which individuals get infected when they are exposed at higher and higher concentrations of virion particles in the water they drink. In fact, in this context (see Eq. (A29)), the parameter ρ represents approximately the rate at which individual ducks get sick at viral concentrations in the water about the ID 50 . Given the quantity of water birds drink daily (7 to 70 l per day), at moderate concentrations of viral particles in the water, they ingest during a day several times the dose that would infect them if they took that dose as a shot in one go. However, probably the same quantity of viral particles have a very mild effect if they are diluted during a day compared to concentrated in a single shot. Experimental tests to estimate ID 50 are usually done by providing a single shot of highly concentrated infected water. This makes it difficult to translate lab results to the field.
Finally, notice that the exponentially saturating functional form for the environmental transmission (see Eq. A24), suggested by Breban et al (2009), relies, as these authors already noted, on the Markovian assumption that the probability that a duck escapes infection when exposed to a virial concentration V does not depend on how much that duck has been exposed to the disease in the past, so it can be modeled by the exponential distribution. However, the derivation presented here makes clear the different approximations that such a representation involves. In addition, this simple analysis has also allowed us to give a rough estimate for the typical values of some of our model parameters.

A4 Sensitivity analysis
To better understand the effects of the uncertainty of parameter values on fluctuating periodicity, we performed a sensitivity analysis of the periodicity with respect to two model parameters: the environmental transmission rate ρ, and the host recovery rate, γ. In Fig. 1, we show the variation of the dominant period for different values of ρ and γ. The dominant period was calculated by averaging wavelet power spectra (see section wavelet power spectrum in the main text) over 20 independent simulations. On the left panel, we can see that decreased environmental transmission rate can significantly lengthen the period of the disease outbreak (P < 0.001). By comparison, ρ should increase to higher values to see a significant decrease in the main periodicity of time-series. Since both rates are given in the same units, we can also safely say that periodicity is less sensitive to changes in the the recovery rate γ than it is to changes in the environmental transmission rate ρ.
We also checked that when the parameter δ, representing the reproduction of the infectious agents in alternative hosts different from the focal host, tend to zero, our main conclusions hold. For instance, when we plot the theoretical PSDs given by Eqs (7) by taking δ = 0, the reported shift to longer dominant periods as environmental transmission decreases is also seen (comparare Fig. 2 Table 2. Error bars denote ±1 mean standard errors. The capital letters A and B denote whether or not there are significant differences in the average dominant period across the different insilico experiments based on Tukey's honest significant difference.

Figure S 2. (Online version in colour)
The qualitative behavior of the PSD when the parameter δ tends to zero is the same as reported in the main text. The frequency values at which these curves peak do not differ significantly from those reported in Fig. 3b