Oscillations in two-species models: tying the stochastic and deterministic approaches Sebasti´an Risau-Gusman1 and Guillermo Abramson2 1Centro At´omico Bariloche and CONICET, 8400 S. C. de Bariloche, Argentina 2Centro At´omico Bariloche, CONICET and Instituto Balseiro, 8400 S. C. de Bariloche, Argentina (Dated: February 7, 2008) Weanalyzegeneraltwo-speciesstochasticmodels,ofthekindgenerallyusedforthestudyofpop- ulation dynamics. We show that the conditions for the stochastic (microscopic) model to display approximatesustainedoscillatorybehavioraregovernedbytheparametersofthecorrespondingde- terministic(macroscopic)model. Weprovideaquantitativecriterionforthequalityofthestochastic oscillation, using a dimensionless parameter that depends only on the deterministic model. When 7 this parameter is small, the oscillations are clear, and the frequencies of the stochastic and deter- 0 0 ministicoscillations areclose, forallstochasticmodels compatiblewiththesamedeterministicone. 2 On theother hand,when it is large, theoscillations cannot be distinguished from a noise. n PACSnumbers: 87.23.Cc,02.50.Ey,05.40.-a a J 6 It is well known that the dynamics of many systems foralargeclassofstochastictwo-speciessystemsandwe 2 of two species can display an oscillatory behavior in the show,bymeansofavanKampenexpansionofthemaster populations of both agents. This happens in predator- equation, that the information given by the determinis- ] E prey systems [1], in models of measles epidemics [2], in tic system—embodied in a deterministic parameter—is P chemicalsystemssuchasthoseexemplifiedbytheBrusse- enough to provide good bounds on this quality. In other . lator[3],etc. Thesesystemsareusuallymodelledbyaset words, we show that the quality of the oscillations is o i of two coupled ordinary differential equations, which are onlyweaklydependentonthedetailsofthedemographic b assumed to represent a macroscopic level of description noise. Moreover, it is shown that oscillations become - q of the system. Oscillations can appear in these mod- clear if and only if the deterministic parameter vanishes. [ els as limit-cycle solutions to the equations. However, We also suggest a heuristic value for the quality below 1 it frequently happens that the macroscopic model only which one can be almost certain that the evolution of v has damped oscillatory solutions, even though the mod- both populations does “look” oscillatory. 4 elledsystemdisplayssustainedoscillationsinits popula- We consider systems of two populations, A and B, de- 4 tions in the same region of parameter values. Examples scribed by stochastic variables m(t) and n(t). The state 0 of this are not uncommon in population dynamics (see, ofthesystemisdefinedbythejointprobabilityP(m,n;t) 1 0 e.g., the discussion in [4] with regard to predator-prey 7 and measles problems). It has often been noted that the 0 stochastic counterpartof these models—assumedto rep- 0.15 b=b =0.1 o/ resent a more microscopic level description of the same dAA=0.B1A3 sinufseccetpetdible i system— usually do display a kind of sustained oscilla- d=0.49 b B B - torybehavior,withafrequencyverysimilartotheoneof f, 0.10 p=3.02 q thedampedsolutionsofthedifferentialequations[5,6,7] A f v: (seeFig.1foranexamplebasedonasusceptible-infected , Xi epidemic model). These oscillations are said to be gen- Wn/ 0.05 erated by the demographic, or intrinsic, noise [8]. The , r a problem is that stochasticity precludes a clear-cut defi- Wm/ nition of “oscillations” for such systems. Therefore, the 0.00 comparisonbetweenthe results ofthe stochasticand de- 0 200 400 600 terministicapproachesisoftenmadeonaqualitativeba- t sis. In this Letter we address the problem of sustained os- FIG. 1: Deterministic dynamics (smooth lines) and one cillations in stochastic models, in an attempt to charac- stochastic realization (fluctuating lines) of an SI epidemic terizetheoscillatoryregime. Weshowthattheconditions model (susceptible and infected, respectively A and B). The for well defined oscillations are given by the parameters dynamics includes birth and death processes in both popu- of the corresponding macroscopic (deterministic) model lations, and contagion. A self-limiting intraspecific compe- only,disregardingthedetailsofthemicroscopic(stochas- tition mechanism is implemented as in [8, 10], with a to- tic) one. In other words, this general result proves that tal system size Ω = 5×105. Transition rates are: T10 = 2(bAφA+bBAφB)(1−φA−φB),T−10 =dAφA,T0−1 =dBφB, theconditionsarethesameforany stochasticmodelthat T−11 =2pφAφB (see Eq. (1)). corresponds to the same macroscopic one. To this end, we define a criterion of “quality” of oscillatory behavior 2 thatthesystemhasmindividualsofspeciesA,andnin- macroscopic considerations of the population dynamics, dividuals of species B. The transition from a state with disregarding its individual level origin. To analyze the (m,n)individualstoastatewith(m+i,n+j)individuals differences between the stochastic (individual level) and takes place at a rate: the deterministic (population level) approaches one usu- m n allychoosesastochasticmodelthatgivestherightdeter- T(m+i,n+j m,n)=f(Ω)Tij( , ), (1) ministic equations. In the limit of infinite size (Ω ), | Ω Ω →∞ Eqs. (2) are also satisfied by the average populations. where k < i < k and k < j < k. Ω is a scale pa- − − The deterministic equilibria are obtained by solving rameter that governs the fluctuations of the stochastic the system C (φ ,φ ) = C (φ ,φ ) = 0, and their A A B B A B evolution. Its precise definition depends on the system, stability is studied by means of a linear stability analy- but one chooses it in such a way that for large Ω the sis. When the system is close to a stable equilibrium, its fluctuations are small. It usually represents the volume evolution can be approximated by that of a damped os- containing the reactants in chemical systems [9], or the cillator. Intheunderdampedregime,thedampingfactor available resources in biological ones [8]. The constant and the frequency of oscillation are, respectively: k gives the maximal number of elements that can ap- pear, or disappear, from a given population at each step γ = ∆ǫ/2, of the dynamics. The most common choice are one-step ω2 = ∆(1 ǫ2/4), (4) d − processes, with k =1. with The evolution of the probability P(m,n;t) is given by the master equation [9]: ǫ= T /√∆, (5) | | ∂P(m,n;t) m i n j where ∆ is the Jacobian of C~, and T its trace: = P(m i,n j;t)T ( − , − ) ij ∂t − − Ω Ω Xij ∆ = C C C C , A,A B,B A,B B,A m n − −P(m,n;t) Tij(Ω,Ω), (2) T = CA,A+CB,B, (6) Xij where C = ∂Ci. The underdamped regime is there- where,asintherestofthisLetter,thesummationindices i,j ∂φj fore given by the condition ǫ < 2. We show below that run from k to k. thisparameter,whichdependsonlyontheparametersof − Except for a few simple cases, this equation is ex- the macroscopicEq. (3), plays a fundamental role in the tremely difficult to solve exactly. For this reason many characterization of the oscillations of stochastic origin. methodshavebeendevisedtolookforapproximatesolu- Notice also that the number of oscillations observed in tions. Perhaps the best known, and most applied, is the the characteristic time γ−1 depends only on ǫ (for small van Kampen expansion [9]. In the following we sketch ǫ, it is just 2/ǫ). the main steps leading to the series solution (a detailed ThefollowingorderinthevanKampenexpansiongives account can be found in van Kampen’s book [9]). the evolution of Π(ξ ,ξ ,t), the joint probability func- A B If one assumes that, at time zero, the system is in a tion of the fluctuations, in the form of a Fokker-Planck state where both populations have well defined macro- equation. To look for oscillations in the fluctuations it scopicalvalues, P(m,n)=δ(m m0)δ(n n0), with the iseasiertoworkwiththeequivalentLangevinequations, − − initialvaluesoforderO(Ω),itisreasonabletoexpectthat as shown by McKane [8]: atlatertimesP(m,n)willhaveasharppeakatsomepo- sition of order O(Ω) (in both populations), and a width ξ˙A = CA,AξA CA,BξB +LA(t) − − of order O(Ω1/2). That is, the fluctuating populations ξ˙ = C ξ C ξ +L (t) (7) B B,A A B,B B B will satisfy m = Ωφ +√Ωξ and n = Ωφ +√Ωξ , − − A A B B where L (t) and L (t) are delta-correlated Gaussian where the variables φ represent the “macroscopic” evo- A B ′ noises of zero mean, satisfying L (t)L (t) =D δ(t lution, while the stochastic variables ξ represent fluctu- A A A ′ ′ h′ i ′ − t), L (t)L (t) = D δ(t t), and L (t)L (t) = ations around them. Replacing this in Eq. (2), equating B B B A B h ′ i − h i D δ(t t). The noise intensities are given by: terms of the same order in Ω and adequately rescaling AB − the time, one obtains, for the leading order: 2 D (φ ,φ ) = i T (φ ,φ ), A A B ij A B φ˙ = iT (φ ,φ ) C (φ ,φ ), Xi,j A ij A B A A B Xij ≡ DB(φA,φB) = j2Tij(φA,φB), (8) φ˙ = jT (φ ,φ ) C (φ ,φ ). (3) Xi,j B ij A B B A B ≡ Xij DAB(φA,φB) = ijTij(φA,φB). Theseequations,calleddeterministic ormacroscopic,are Xi,j usuallythestartingpointofmanymodelsofchemicaland By Fourier transforming Eqs. (7) it is straightforward biologicalsystems. Theyaregenerallywrittendownfrom to obtain the power spectrum of the fluctuations around 3 For small ǫ, ωˆ2 tends to 1, which means that not only A 1.0 the frequencies of possible oscillations for both popula- Fig. 3(b) tions become close, but also that they become close to 0.8 ω , the frequency of the damped oscillations of the de- d terministic model (which also tends to 1 as ǫ 0). It 2 B0.6 → 2^ ^w , A 0.4 Fig. 3(a) idsyinnamthiicssrcehgiamraecttehraitstticheopfosptouclahtaisotnics sohsocwillatthieoncso.heTrehnist w mityotoiofnthweisllpbecetrfuumrthpeerackh.aFraigcuterreiz2edshboewloswωˆ2byanthdeωˆq2uaals- 0.2 A B susceptible functions of ǫ for the SI model presented in Fig. 1, for infected 0.0 a wide range of system parameters. The bounds given 0 1 2 3 4 by Eq. (12) are shown by continuous lines. Each point e represents the normalized squared frequency for one set of parameters, for both populations. The deterministic FIG. 2: ωˆA2 and ωˆB2 for the same SI model as in Fig. 1. The frequency, ω , is also shown, to emphasize the difference pointscorrespondtoauniformscanningofaportionofphase d between the three frequencies present in the system. space: bA = bBA = 0.1, dA ∈ (0,0.2), dB ∈ (0,0.5), p ∈ (2dB,3.2). The full lines show the bounds of Eq. (12), while When√2<ǫ<2therecanbesomestochasticmodels the dashed one corresponds to ωd. The arrows point to the for which no peak is present in S(ωA) or S(ωB). And, values corresponding to theparameters used in Fig. 3. for some values of F or F , it can happen that the A B power spectrum of any population has a maximum even ifǫ>2,i.e.evenwhenthedeterministicsystemdoesnot display damped oscillations (see Fig. 2: all the points to thedeterministicequilibrium[8]. Inthefollowingwecon- the right of ǫ=2 correspond to systems with a peak in centrateonpopulationA.Thecorrespondingexpressions the spectrum of the susceptible (A) population, no peak forpopulationB areobtainedbyexchangingA andB in in the infected (B) one, and no damped oscillations in all the subindices. The average power spectrum of ξ is A the deterministic model). These two features show that F +ωˆ2 the peaks of the stochastic power spectrum on the one A hSA(ω)i= (1 ωˆ2)2+ωˆ2ǫ2, (9) hand, and the deterministic damped oscillations on the − other, are not necessarily closely related. where Theabovediscussionestablishestheconditionsforthe existenceofapeakinthe powerspectrumofoneorboth 2 2 ωˆ =ω /∆, populations. That is, for the existence of a preferred C2 D +C2 D 2C C D frequency in their dynamics. But, should all peaks in FA = A,B B B,B A− A,B B,B AB. (10) the power spectrum be regarded as “oscillations”? The ∆D A answer to this question is certainly negative, and leads WestressthatF (andcorrespondinglyF ),throughits one to look for a criterion to quantify how close a time A B dependenceonC ,C ,D andD ,dependsultimately seriesistoanoscillatorymovement. Thiscanbedoneby A B A B on the transition probabilities that define the model. defining the “quality” of the oscillation as a measure of Itisstraightforwardtoseethat S (ω) iseithermono- thesharpnessofthepeak. Weproposeonesuchmeasure A h i tonically decreasing or it has a single maximum at in the following. Given a power spectrum of the form (9) we define the ωˆA2 =−FA+ (FA+1)2−ǫ2FA. (11) quality of a peak at ωpeak as p Theconditionofpositivityfortheargumentofthesquare ωpeak SA(ωpeak) Q (ω)= h i. (13) rootgivestheregioninphasespacewherethepowerspec- A S (ω) dω A h i trum has a singlemaximum. Notice that for ǫ<√2 this R This quantity is dimensionless and scale invariant. It conditionisfulfilledregardless of the exact dependence of is related to Fisher’s kappa, which measures the non- F on the parameters of themodel. Itcanalsobeproved A stationarity of a signal, given its periodogram [11]. For 2 that ωˆ satisfies [12]: A functions with only one peak, Q increases as the peaks A 2 2 grows. Forpowerspectraoftheform(9),Q canberead- 1 ǫ /2<ωˆ <1 (12) A − A ily calculated (using that S (ω) dω = ξ2 , see [9]): h A i h Ai (these bounds seem to be tight). In particular, this im- R 2 ωˆ ǫ F +ωˆ plies that in all the possible stochastic models that lead Q (ω )= A A A . (14) A A (ωˆ2 1)2+ωˆ2ǫ2 (cid:18) F +1 (cid:19) to the same deterministic equations (same C’s, different A− A A D’s) the position of the maximum can only vary within ThequalityQ divergesasǫvanishes,regardlessofthe A a finite range, that shrinks with ǫ. exactdependence ofF onthe parametersofthe model. A 4 0.02 Q = 2.11, e = 0.30 10 susceptible i 8 infected b=b =0.1 S( f )6 A BA d=0.13, d=0.49 4 1 A B 2 p=3.02 00.00 0.03 0.06 0.09 2 f ep W/ Q n 2 2 1-e 2/2 0.1 ep1 +e 2/4 (a) 0.01 1500 1600 1700t 1800 1900 2000 0 1 e 2 3 4 0.095 Qi = 0.25, e = 1.4 4 FIG. 4: QA and QB as functions of ǫ for the SI model of b =b =0.1 S( f )3 Fig. 1. The points correspond to the same portion of phase A BA 2 d =0.05, d=0.25 spaceasinFig.2. Thelinesshowtheupperandlowerbounds A B 1 0.090 p=2.2 0 of Eq. (15). 0.00 0.03 0.06 0.09 f W/ n series will indeed “look” oscillatory? As it is to be ex- 0.085 pected,thecontinuousnatureofQprecludesaconclusive (b) answer. From exhaustive observations of different mod- 500 600 700 800 900 1000 els we find that when Q(ωpeak) > 1, the oscillations are t welldefined and notably different froma noisy evolution (see Fig. 4). Insummary,wehaveshownthat,bydefiningaquality FIG.3: TwostochasticrealizationsoftheSImodel,withdif- measure,onecanquantifythe“oscillatorylook”ofatime ferentqualities. Theinsetsshowthecorrespondinganalytical averagepowerspectra. Onlytheinfectedpopulationisshown. series. Interestingly, we find that oscillations are present The arrows in Fig. 2 point to the corresponding frequencies: only when ǫ is small. This means that, given a deter- ωA ∼ ωB ∼ ωd for the good quality case shown in (a), and ministic model, one can know, using the bounds (15), ωA≫ωB for thebad quality one shown in (b). Ω=105. whether the time series given by any stochastic counter- partofthemodelwilllookoscillatoryornot. Inaddition, we have shown that, when oscillations are clear, the cor- Therefore, one can assure that the corresponding time responding frequencies of both populations will be close serieswilllookoscillatorywhenǫissufficientlysmall(see to each other and to the frequency of the damped oscil- Fig. 3 for an example of this). In such a case, we have lation of the deterministic system. already shown that the frequencies of both populations Given that our conclusions are based on the analysis areveryclose,andalsoveryclosetothefrequencyofthe of the first two terms of the systematic van Kampen’s deterministic damped oscillations. expansionof the master equation, they are exact only in Could it also happen that, for large values of ǫ, when the limit Ω . These analytical results, nevertheless, →∞ the frequencies of populations A and B can be rather compare well with the numerical observations made on different, one gets very sharp peaks? It can be shown finite systems. More details about the validity of the that this is not the case by giving bounds of Q that expansion for finite systems will be given elsewhere [12]. A depend solely on ǫ [12]: WearegratefultoE.Andr´es,I.Peixoto,A.Aguirre,L. 2 1 ǫ2/2 2 On˜a and H. Solarifor valuable discussions. We acknowl- − <Q (ω)< . (15) edge financial support from ANPCyT (PICT-R 2002- πǫ(cid:18)1+ǫ2/4(cid:19) A πǫ 87/2),CONICET(PIP5414)andUNCUYO (06/C209). The upper bound shows that, when ǫ is not small, the peak cannot be arbitrarily sharp. On the other hand, the lower bound shows that, when ǫ is small the peak is sharpforallthestochasticcounterpartsofadeterministic model. In Fig. 4 we illustrate this by showing several [1] M. Begon, C. R. Townsend and J. L. Harper, Ecology: From Individuals to Ecosystems (Blackwell, 2006). values of Q and Q for the SI model, along with the A B [2] E. B. Wilson and O. M. Lombard, Pathology 31, 367 corresponding bounds. (1945). One practical question remains: what is the critical [3] I. Prigogine and R. Lefever, J. Chem. Phys. 48, 1695 quality value above which one can be sure that the time (1968). 5 [4] E. Renshaw, Modelling Biological Populations in Space 218102 (2005). andTime,(Cambridge,1991).SeeSections6.2and10.4. [9] N. G. van Kampen, Stochastic Processes in Physics and [5] M. S. Bartlett, J. R.Stat. Soc. A 120, 48 (1957). Chemistry (Elsevier Sience B.V., Amsterdam, 2003). [6] H.W.HethcoteandS.A.Levin,inAppliedMathematical [10] A. J. McKane and T. J. Newman, Phys. Rev. E 70, Ecology, L. Gross, T. G. Hallam and S.A. Levin, (eds.), 041902 (2004). pp.193-211 (Springer,Berlin, 1989). [11] R. A.Fisher, Proc. R.Soc. Lon. A 125, 54 (1929). [7] J. P. Aparicio and H. G. Solari, Math. Biosciences 169, [12] S. Risau-Gusman and G. Abramson, in preparation. 15 (2001). [8] A. J. McKane and T. J. Newman, Phys. Rev. Lett. 94,