Fluctuations and oscillations in a simple epidemic model
aa r X i v : . [ q - b i o . P E ] J un Fluctuations and oscillations in a simple epidemic model
G. Rozhnova and A. Nunes Centro de F´ısica Te´orica e Computacional and Departamento de F´ısica,Faculdade de Ciˆencias da Universidade de Lisboa, P-1649-003 Lisboa Codex, Portugal
We show that the simplest stochastic epidemiological models with spatial correlations exhibit twotypes of oscillatory behaviour in the endemic phase. In a large parameter range, the oscillationsare due to resonant amplification of stochastic fluctuations, a general mechanism first reported forpredator-prey dynamics. In a narrow range of parameters that includes many infectious diseaseswhich confer long lasting immunity the oscillations persist for infinite populations. This effect isapparent in simulations of the stochastic process in systems of variable size, and can be understoodfrom the phase diagram of the deterministic pair approximation equations. The two mechanismscombined play a central role in explaining the ubiquity of oscillatory behaviour in real data and insimulation results of epidemic and other related models.
PACS numbers: 87.10.Mn; 87.19.ln; 05.10.Gg
Cycles are a very striking behaviour of prey-predatorsystems also seen in a variety of other host-enemy sys-tems — a case in point is the pattern of recurrent epi-demics of many endemic infectious diseases [1]. The con-troversy in the literature over the driving mechanismsof the pervasive noisy oscillations observed in these sys-tems has been going on for long [2], because the sim-plest deterministic models predict damped, instead ofsustained, oscillations. One of the aspects of this contro-versy is whether these mechanisms are mainly external orintrinsic, and the effects of seasonal forcing terms [3], [4]and of higher order non-linear interaction terms [5] havebeen explored in the framework of a purely determinis-tic description of well-mixed, infinite populations. Thesemore elaborate models exhibit oscillatory steady statesin certain parameter ranges, and have led to successfulmodelling when external periodic forcing is of paramountimportance [6], but they fail to explain the widespreadnon-seasonal recurrent outbreaks found, for instance, inchildhood infectious diseases [7].During the last decade, important contributions havecome from studies that highlight the inherently stochasticnature of population dynamics and the interaction pat-terns of the population as important endogenous factorsof recurrence or periodicity [8]. A general mechanism ofresonant amplification of demographic stochasticity hasbeen proposed to describe the cycling behaviour of prey-predator systems [9] and applied recently to recurrentepidemics of childhood infectious diseases [10]. The roleof demographic stochasticity modelled as additive Gaus-sian white noise of arbitrary amplitude in sustaining in-cidence oscillations had long been acknowledged in theliterature [11]. The novelty in [9] and [10] was that of pro-viding an analytical description of demographic stochas-ticity as an internal noise term whose amplitude is de-termined by the parameters and the size of the systemusing a method originally proposed by van Kampen [12].Our goal is to extend this approach by relaxing thehomogeneous mixing assumption to include an implicit representation of spatial dependence. We show that theinclusion of correlations at the level of pairs leads to dif-ferent quantitative and qualitative behaviours in a re-gion of parameters that corresponds to infectious diseaseswhich confer long lasting immunity. Our motivation wastwofold. On one hand, the homogeneous mixing assump-tion is known to give poor results for lattice or networkstructured population [13], [14]. On the other hand, sys-tematic simulations of infection on small-world networkshave shown that the resonant amplification of stochasticfluctuations is significantly enhanced in the presence ofspatial correlations [15]. Therefore, apart from stochas-ticity, the correlations due to the contact structure areanother key ingredient to understand the patterns of re-current epidemics. One of the main difficulties in in-cluding this ingredient is that the relevant network ofcontacts for disease propagation is not well known [14].In this paper we shall consider a stochastic Susceptible-Infective-Recovered-Susceptible (SIRS) epidemic modelthat leads to the ordinary pair approximation (PA) equa-tions of [13] in the thermodynamic limit as the simplestrepresentation of the spatial correlations on an arbitrarynetwork of fixed coordination number k . The power spec-trum of the fluctuations around the steady state can becomputed following the approach of [9] and [12]. Thecombined effect of stochasticity and spatial correlationshas been much studied through simulations, but this isan analytical treatment of a model that includes boththese ingredients.Consider then a closed population of size N at a giventime t , consisting of n individuals of type S , n indi-viduals of type I , and ( N − n − n ) individuals of type R , modelled as network of fixed coordination number k .Recovered individuals lose immunity at rate γ , infectedindividuals recover at rate δ , and infection of the suscep-tible node in a susceptible-infected link occurs at rate λ .Let n (respectively, n and n ) denote the number oflinks between nodes of type S and I (respectively, S and R and R and I ). In the infinite population limit, withthe assumptions of spatial homogeneity and uncorrelatedpairs, the system is described by the deterministic equa-tions of the standard or uncorrelated PA as follows [13]:˙ p = γ (1 − p − p ) − kλp , (1)˙ p = kλp − δp , ˙ p = γp − ( λ + δ ) p + ( k − λp p ( p − p − p ) , ˙ p = δp + γ (1 − p − p − p − p ) − ( k − λp p p , ˙ p = δ ( p − p ) − ( γ + 2 δ ) p + ( k − λp p p . In the above equations the variables stand for the limitvalues of the node and pair densities p = n /N , p = n /N , p = n / ( kN ), p = n / ( kN ), and p = n / ( kN )as N → ∞ . As expected, neglecting the pair correlationsand setting p = p p in the first two equations leads tothe classic equations of the randomly mixed or mean-fieldapproximation (MFA) SIRS model,˙ p = γ (1 − p − p ) − kλp p , (2)˙ p = kλp p − δp . The phase diagrams of the two models are plotted in Fig.1 [solid lines for Eq. (1) and dashed line for Eq. (2), bothwith k = 4]. We have set the time scale so that δ = 1.The critical line separating a susceptible-absorbing phasefrom an active phase where a stable steady state existswith nonzero infective density is given by λ MFA c = (dashed black line) for the MFA, and by λ PA c ( γ ) = γ +13 γ +2 [solid orange (gray) line] for the PA. In addition, in theactive phase of the PA we find for small values of γ a newphase boundary [solid blue (black) line] that correspondsto a Hopf bifurcation and seems to have been missed inprevious studies of this model [13]. This boundary sep-arates the active phase with constant densities (regionI) from an active phase with oscillatory behavior (regionII). The maximum of the curve is situated at λ ≈ . γ ≈ .
03, which means that the PA model predicts sus-tained oscillations in the thermodynamic limit when lossof immunity is much slower than recovery from infection.According to published data for childhood infections inthe pre-vaccination period [4], taking the average immu-nity waning time to be of the order of the length of theelementary education cycles at that time (10 years forthe data points in Fig. 1) many of the estimated param-eter values for these diseases fall into oscillatory regionII, and the others are in region I close to the boundary.Different data points for the same disease correspond toestimates for λ based on different data records. Althoughsmall enough to be missed in a coarse grained numeri-cal study, the oscillatory phase is large in the admissibleparameter region of an important class of diseases. Asystematic study of the dependence of this oscillatory λ γ ΙΙ Ι
FIG. 1: (Color online) Phase diagram in the ( λ, γ ) plane forthe MFA and the PA deterministic models and parameter val-ues for measles ( △ ), chicken pox ( ◦ ), rubella ( (cid:3) ) and pertussis( ⋄ ) from data sources for the pre-vaccination period. The bluestars are the parameter values used in Fig. 3. phase on the parameter k and of its relevance to under-stand the behaviour of simulations on networks will bereported elsewhere [16]. Preliminary results indicate thatthe oscillatory phase persists in the range 2 < k .
6, andthat it gets thinner as k increases. There are indicationsthat this oscillatory phase is robust also with respect tovariations of the underlying model [17].Let us now study the combined effect of correlationsand demographic stochasticity in region I by taking N large but finite. In the stochastic version of the MFASIRS model the state of the system is defined by n and n which change according to the transition rates as T n +1 ,n = γ ( N − n − n ) , (3) T n ,n − = δ n ,T n − ,n +1 = kλ n N n , associated to the processes of immunity waning, recov-ery and infection. Here T n + k ,n + k denotes the transi-tion rate from state ( n , n ) to state ( n + k , n + k ), k i ∈ {− , , } , where i = 1 ,
2. As in [9], the powerspectrum of the normalized fluctuations (PSNF) aroundthe active steady state of system (2) can be computedapproximately from the next-to-leading-order terms ofvan Kampen’s system size expansion of the correspond-ing master equation. Setting n ( t ) = N p ( t ) + √ N x ( t )and n ( t ) = N p ( t ) + √ N x ( t ), the equations of motionfor the average densities (2) are given by the leading-order terms of the expansion. The next-to-leading-orderterms yield a linear Fokker-Planck equation for the prob-ability distribution function Π( x ( t ) , x ( t ) , t ). The equiv-alent Langevin equation for the normalized fluctuationsis ˙ x i ( t ) = P j =1 J ij x j ( t ) + L i ( t ), where J is the Jacobianof Eq. (2) at the endemic equilibrium and L i ( t ) are Gaus-sian white noise terms whose amplitudes are given by theexpansion. Taking the Fourier transform we obtain forthe PSNF, P i ( ω ) ≡ h| ˜ x i ( ω ) | i = X j,k M − ik ( ω ) B kj M − ji ( − ω ) , (4)where M ik ( ω ) = i ωδ ik − J ik and h ˜ L i ( ω ) ˜ L j ( ω ′ ) i = B ij δ ( ω + ω ′ ). For k = 4 and δ = 1 this expression be-comes P MFAS = B ( J + ω )( D − ω ) + T ω , (5) P MFAI = B ( J + J J + J + ω )( D − ω ) + T ω , (6)where D and T are the determinant and the trace of J and B = B = − B = − B = γ (4 λ − λ ( γ +1) , for thesusceptible and the infected PSNFs, respectively.In a stochastic version of the PA SIRS model thestate of the system is defined by the integers n i , where i = 1 , ...,
5, and recovery, loss of immunity, and infec-tion induce different transitions according to the pairsor triplets involved in the process. The simplest set oftransitions and transition rates compatible with Eq. (1)is T n +1 ,n ,n + k,n ,n − k = γk n , (7) T n +1 ,n ,n ,n − k,n = γk n ,T n +1 ,n ,n ,n + k,n = γk ( k ( N − n − n ) − n − n ) ,T n ,n − ,n − k,n + k,n = δk n ,T n ,n − ,n ,n ,n + k = δk ( kn − n − n ) ,T n ,n − ,n ,n ,n − k = δk n ,T n − ,n +1 ,n − k,n ,n = λk n n n ,T n − ,n +1 ,n − ,n − ( k − ,n +( k − = λk n n n ,T n − ,n +1 ,n +( k − ,n ,n = λk n n ( kn − n − n ) . This is a coarse grained description where the effect ofthe change in state of a given node on the k pairs that itforms is averaged over each pair type. For instance, theevent of loss of immunity occurs at a rate γn R , where n R is the number of recovered nodes, and changes the k pairs formed by the node that switches from recoveredto susceptible. On average, each pair type will changeby k units at a rate proportional to its density, accord-ing to the equation γn R = γn R ( n kn R + n kn R + n RR kn R ),where n RR is the number of pairs of recovered nearestneighbours. Taking this level of description and using kn R = n + n + 2 n RR and n + n + n R = N , we obtainthe first three equations of Eq. (7) for the rates of the three different pair events associated with loss of immu-nity. A full microscopic description would require con-sidering separately all possible five-node configurationsfor the central node that switches from R to S and itsfour nearest neighbours. We have checked that the de-tailed stochastic model involving 40 different transitionsfor k = 4 gives essentially the same results [16] as thecoarse grained model (7) that we consider here.For the fluctuations of the pair densities we set n ( t ) = N kp ( t ) + √ N kx ( t ), n ( t ) = N kp ( t ) + √ N kx ( t ), and n ( t ) = N kp ( t ) + √ N kx ( t ). The leading order termsof van Kampen’s system size expansion of the masterequation associated to (7) yield the deterministic PA Eqs.(1). An approximate analytical expression for the PSNFcan be obtained as before from the next-to-leading-orderterms. Formula (4) is still valid taking now J as theJacobian of Eq. (1) at the endemic equilibrium and thenoise cross correlation matrix B computed directly fromthe expansion. P I ( ω ) ω P I ( ω ) ω −2 P I ( ω ) ω λ γ A N A N a)c)e) b)d)f)III FIG. 2: (Color online) (a) Analytical and numerical PSNFsof the infectives for model (3) with γ = 0 . λ = 2 .
5; (b)the same for model (7); (c) a similar plot for model (7) with γ = 0 .
034 and λ = 2 .
5, notice the lin-log scale; (d) location ofthe parameter values chosen for (a) and (b) (circle) and for (c)(square); and (e) and (f) plots of the peak amplitude of thePSNF of the PA model as a function of N for the parametervalues chosen for (b) and (c). In Fig. 2 the approximate PSNFs given by Eq. (4)are plotted (black lines) and compared with the numer-ical power spectra of stochastic simulations for N = 10 [green (gray) lines] for models (3) and (7) [Figs. 2(a)and 2(b)]. For this system size, there is almost per-fect agreement between the analytical approximate ex-pressions and the results of the simulations across thewhole region I. The plots show that for the same param-eter values the fluctuations are larger and more coherentfor model (7), in agreement with the results of simula-tions reported in [15] for small-world networks on a lat-tice and variable small-world parameter. This effect be-comes much more pronounced as the boundary betweenregions I and II, where the analytical PSNF of model(7) diverges, is approached. Close to this boundary [seeFig. 2(c)], there is a significant discrepancy between theanalytical (black line) and the numerical [green (gray)line] PSNFs associated with the appearance of secondarypeaks at multiples of the main peak frequency. This isa precursor of the oscillatory phase, and the breakdownof van Kampen’s approximation for this system size maybe understood as an effect of the loss of stability of theendemic equilibrium close to the boundary. Relaxationtowards equilibrium becomes slow compared with the pe-riod of the damped oscillations, and a significant part ofthe power spectrum energy shows up in the secondaryharmonics. For these parameter values, van Kampen’sexpansion becomes a good approximation only for largersystem sizes. Also shown in Fig. 2(e) [respectively, Fig.2(f)] is the scaling with system size of the peak ampli-tude of the infectives PSNF of the PA model [pink (black)dots] for the parameter values considered in (b) [respec-tively, (c)] and the peak amplitude (dashed blue line) ofthe approximated PSNF given by Eq. (4). Away fromthe phase boundary of the oscillatory phase we find thatthe simulations exhibit the amplitude and scaling pre-dicted by Eq. (4) down to system sizes of ∼ . By con-trast, close to the phase boundary the match is reachedonly for system sizes larger than 5 × . t p i t p i a) b) FIG. 3: (Color online) The steady-state infective densitygiven by the PA deterministic model (dashed blue lines), andby simulations of the PA and the MFA stochastic models [solidblack and green (gray) lines, respectively] in regions II and Ifor N = 10 . Parameters are (a) λ = 7 . γ = 0 .
01 and (b) λ = 9, γ = 0 . Examples of typical time series predicted by the PAmodel in the parameter region of childhood infectiousdiseases are shown in Fig. 3 [dashed blue lines for the de-terministic model (1) and solid black lines for simulationsof the stochastic model (7)]. The results of simulationsof the MFA stochastic model (3) for the same parametervalues are also shown for comparison [solid green (gray) lines]. Fig. 3(a) illustrates the regular high-amplitudeoscillations of region II. All over this region, simulationsof the stochastic model (7) reproduce the behaviour ofthe solutions of Eq. (1) with added amplitude fluctua-tions. The only limitation to observe these oscillationsin finite systems is that N has to be taken large enoughfor the deep interepidemic troughs to be spanned with-out stochastic extinction of the disease. In region I [Fig.3(b)] there are no oscillations in the thermodynamic limitbut, in contrast to the stochastic MFA model, the reso-nant fluctuations in the PA model are large and coherentenough to provide a distinct cycling pattern, which ispartially described by van Kampen’s expansion (4).In conclusion, we have considered a stochastic versionof the basic model of infection dynamics including a rep-resentation of the spatial correlations of an interactionnetwork through the standard PA. We have shown thatin general the resonant amplification and the coherenceof stochastic fluctuations are much enhanced with re-spect to the MFA model. This quantitative differencebecomes qualitative in a region of parameter space thatcorresponds to diseases for which immunity waning ismuch slower than recovery. In this region the nonlineari-ties of the model and demographic stochasticity give riseeither to oscillations that persist in the thermodynamiclimit or to high amplitude, coherent resonant fluctua-tions, providing realistic patterns of recurrent epidemics.These results are relevant for other population dynam-ics models in the slow driving regime that corresponds tosmall γ in our model, suggesting that in systems of mod-erate size intrinsic stochasticity together with the sim-plest representation of spatial correlations may be enoughto produce distinct oscillatory patterns. This favours theview that, for a large class of systems, noisy oscillations inpopulation dynamics data may be intrinsic, rather thanexternally driven.Financial support from the Foundation of the Uni-versity of Lisbon and the Portuguese Foundation forScience and Technology (FCT) under Contracts No.POCI/FIS/55592/2004 and No. POCTI/ISFL/2/618is gratefully acknowledged. The first author (G.R.)was also supported by FCT under Grant No.SFRH/BD/32164/2006 and by Calouste GulbenkianFoundation under its Program ”Stimulus for Research”. REFERENCES [1] P. Rohani, D. J. D. Earn, and B. T. Grenfell, Science ,968 (1999); N. C. Grassly, C. Fraser, and G. P. Garnett,Nature (London) , 417 (2005).[2] O. N. Bjørnstad and B. T. Grenfell, Science , 638(2001).[3] M. J. Keeling, P. Rohani, and B. T. Grenfell, Physica D , 317 (2001). [4] C. T. Bauch and D. J. D. Earn, Proc. R. Soc. London,Ser. B , 1573 (2003).[5] H. W. Hethcote and P. van den Driessche, J. Math. Biol. , 271 (1991); W. Wang, Math. Biosci. Eng. , 267(2006).[6] D. J. D. Earn, P. Rohani, B. M. Bolker, andB. T. Grenfell, Science , 667 (2000); B. T. Grenfell,O. N. Bjørnstad, and J. Kappey, Nature (London) ,716 (2001).[7] C. T. Bausch, in Mathematical Epidemiology , edited byF. Brauer, P. van den Driessche, and J. Wu (Springer,Berlin, 2008), p. 297.[8] J. E. Satulovsky and T. Tom´e, Phys. Rev. E , 5073(1994); A. Lipowski, Phys. Rev. E , 5179 (1999);M. Kuperman and G. Abramson, Phys. Rev. Lett. ,2909 (2001); T. Gross, Carlos J. Dommar DLima, andB. Blasius, Phys. Rev. Lett. , 208701 (2006).[9] A. J. McKane and T. J. Newman, Phys. Rev. Lett. ,218102 (2005). [10] D. Alonso, A. J. McKane, and M. Pascual, J. R. Soc.,Interface , 575 (2007).[11] M. S. Bartlett, Stochastic Population Models in Ecologyand Epidemiology (Methuen, London, 1960).[12] N. G. van Kampen,
Stochastic Processes in Physics andChemistry (Elsevier, Amsterdam, 1981).[13] J. Joo and J. L. Lebowitz, Phys. Rev. E , 036114(2004).[14] M. J. Keeling and K. T. D. Eames, J. R. Soc., Interface , 295 (2005).[15] M. Sim˜oes, M. M. Telo da Gama, and A. Nunes, J. R.Soc., Interface , 555 (2008).[16] G. Rozhnova and A. Nunes, e-print arXiv:0812.1812.[17] D. A. Rand, in Advanced Ecological Theory: Principlesand Applications , edited by J. McGlade (Blackwell Sci-ence, Oxford, 1999), p. 100; J. Benoit, A. Nunes, andM. M. Telo da Gama, Eur. Phys. J. B50