Oscillations in two-species models: tying the stochastic and deterministic approaches

0.35 MB·English

by  Sebastián Risau-Gusman

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,