— physics pp. 1047–1053
Assessing the quality of stochastic oscillations
GUILLERMO ABRAMSON∗ and SEBASTI ´AN RISAU-GUSMAN
Statistical and Interdisciplinary Physics Group, Centro At´omico Bariloche, CONICET and Instituto Balseiro, R8402AGP Bariloche, Argentina
∗Corresponding author. E-mail: g.abrambon@gmail.com
Abstract. We analyze the relationship between the macroscopic and microscopic de- scriptions of two-state systems, in particular the regime in which the microscopic one shows sustained ‘stochastic oscillations’ while the macroscopic tends to a fixed point. We propose a quantification of the oscillatory appearance of the fluctuating populations, and show that good stochastic oscillations are present if a parameter of the macroscopic model is small, and that no microscopic model will show oscillations if that parameter is large.
The transition between these two regimes is smooth. In other words, given a macroscopic deterministic model, one can know whether any microscopic stochastic model that has it as a limit, will display good sustained stochastic oscillations.
Keywords. Population dynamics; stochastic oscillations.
PACS Nos 87.23.Cc; 02.50.Ey; 05.40.-a
1. Bridging the gap between the mean field and individual-based models The subject of this contribution is a spin-off of an interdisciplinary research effort on the epidemiology of Hantavirus in wild rodent populations. This project, de- voted to providing a theoretical framework for ecological and epidemiological field findings, has explained and predicted a number of phenomena in relation with this widespread epizootic (see, e.g. [1–3]). Most of that work is based on the analysis of mean field-like models, continuous models of the kind generally used for large populations, chemical kinetics, etc. They constitute a macroscopic,population level description of the system under study. While powerful and useful, they are unsatis- factory in many regards. For fundamentally discrete and not very large populations, a complementary microscopic,individual level description is called for. Moreover, the stochastic nature of the processes and interactions that take place at the indi- vidual level are responsible for a fluctuating dynamics that can only be captured by a stochastic model. One expects that the continuous macroscopic models, in certain sense and regimes in which one-individual changes are small, constitute a sensible limit of the microscopic ones. The determination of those regimes becomes a pre-eminent question [4,5]. In addition, there exists the possibility of some phe- nomena that, due to fluctuations, are exclusive to the microscopic models. The
0 200 400 600 0.00
0.05 0.10 0.15
m/Ω, n/Ω, φ A, φ Β
t
susceptible infected bA=bBA=0.1
dA=0.13 dB=0.49 p=3.02
Figure 1. Persistence of oscillations in a microscopic SI model. The fluctuat- ing populations of susceptibles and infected in a single realization are shown, together with the smooth lines corresponding to both the deterministic model and the average of 10,000 stochastic realizations, which are indistinguishable.
stabilization of stochastic oscillations, which is the subject of this work, is one of these.
This interesting phenomenon is displayed by certain stochastic microscopic mod- els, whose macroscopic counterpartsdo not oscillate. The oscillations arise purely from fluctuations, by the collective synchronization of individual transitions through their interactions. They may explain the non-seasonal oscillations observed in many epidemiological and ecological systems (such as in measles or chicken pox [6,7]). Mi- croscopic fluctuations are responsible for other interesting phenomena in chemical or chemical-like systems, such as the transition between a chaotic attractor and a fixed point suppression of chaotic oscillations [8–10].
Figure 1 shows the persistence of oscillations exemplified in a susceptible-infected (SI) model. The smooth analytic curves show just damped oscillations. So does an average of 10,000 realizations of a microscopic model. A single realization, however, clearly shows persistent oscillations. In this work we characterize the regime in which the stochastic oscillations are clearly visible against the random fluctuations, as exemplified in figure 2.
2. Formulation of the stochastic model
The formulation of a stochastic model proceeds by identifying the ‘microscopic’
processes first, and writing down a master equation next, i.e. an evolution equation for the joint probability distribution of the system.
Let us consider a population, composed of m individuals of one kind, and nof another kind (for example, susceptible and infected animals in an epidemic sys- tem), that dwell in a region of volume (area) Ω. (The generalization to systems with more than two species is in principle possible, while much of the analysis still needs to be carried out in detail.) The state of the system is defined by the joint probability distribution P(m, n, t). The animals are born, reproduce, and die ac- cording to certain rules, which are a simplification of the corresponding processes in a natural population. They are also converted from one class into the other, such as in the contagion of a susceptible individual after contact with an infected
1500 1600 1700 1800 1900 2000 0.01
0.02
Q = 6.6
n /
Ωt
0.00 0.03 0.06 0.0902468 10
S(f) f
500 600 700 800 900 1000
0.085 0.090 0.095
Q = 0.79
n /
Ωt
0.00 0.03 0.06 0.090 1
2 3 4
S(f) f
Figure 2. Infected population of a microscopic SI model, showing two regimes in which the power spectrum has a peak, but with very different oscillating characteristics. Insets: Power spectra of the time series.
one. These processes correspond to transitions in the configuration space of the system, and we suppose that they occur stochastically, according to certain tran- sition probabilitiesW(m0, n0|m, n) (the rate at which the system goes from state (m, n) to state (m0, n0)).
The master equation is the gain–lose equation that describes the evolution ofP: dP(m, n)
dt = X
m0,n0
W(m, n|m0, n0)P(m0, n0)
− X
m0,n0
W(m0, n0|m, n)P(m, n), (1)
where the sums are taken over all the states that have allowed transitions to, or from, the state (m, n). A common situation is that of one-step processes, for which we write the transition probabilities asWij, withi, j∈ {−1,0,1}.
3. Expansion of the master equation
Even though it is always a linear equation, the master equation can rarely be solved exactly, and either an approximate method, or a numerical scheme, or a combination of both, needs to be used. As shown by McKane and Newman [11], the systematic expansion of the master equation devised by van Kampen [12] is particularly adequate in this case. The expansion relies on the following ansatz, that assumes the existence of a size parameter Ω, and decomposes each stochastic variable into a macroscopic, analytical one, plus a noise of smaller amplitude:
m(t) = Ωφ1(t) +√
Ωξ1(t), (2)
and correspondingly for n(t). The use of the ansatz (2) into eq. (1) leads to an expansion in powers of Ω−1/2. Equating successive orders of the series provides
corresponding orders of dynamical description: the leading order consists of equa- tions for φ1(2) only, that one identifies as the macroscopic (mean field) law; the following order consists of a Fokker–Planck equation for the probability distribu- tion function of the fluctuationsξ1(2), which is the main component of the dynamics of the fluctuations; higher orders provide additional levels of detail in the fluctuating dynamics.
The leading order in the expansion is that of√
Ω. This needs to be made null, for the expansion to behave regularly when Ω→ ∞. One arrives, at this order, at the following set of two ordinary differential equations for the macroscopic variables:
dφ1
dt =X
ij
iWij(φ1, φ2)≡C1(φ1, φ2), (3a) dφ2
dt =X
ij
jWij(φ1, φ2)≡C2(φ1, φ2). (3b) These are the equations of the mean field, macroscopic level of description. They give the dynamics of the smooth population variables that can be observed, for example, averaging over realizations in a numerical simulation of the system (see figure 1). Observe, also, that the summation that appears in eqs (3) implies that many microscopic systems(corresponding to different choices of the transition prob- abilities)map onto the same macroscopic model.
A linear stability analysis shows that a stable fixed point exists in eqs (3) if
²= T
√∆
½<2, stable spiral,
>2, stable node, (4)
where T and ∆ are, respectively, the trace and the determinant of the matrix Cij =∂Ci/∂φj. The case ²= 2 corresponds to a center which, being structurally unstable, is of lesser interest in applications to real problems. The case² <2 is of particular interest, since it corresponds todamped oscillations of the macroscopic model.
The following order in the expansion, found by equating the terms of order Ω0, represents the main behavior of the deviations from the macroscopic variables. It takes the form of a linear Fokker–Planck equation for the probability distribution function of the fluctuationsξ1 andξ2. It has the form
∂Π(ξ1, ξ2, t)
∂t =X
ij
aij ∂
∂ξi
ξjΠ(ξ1, ξ2, t) +X
ij
bij ∂2
∂ξi∂ξj
Π(ξ1, ξ2, t), (5) with coefficientsaijandbijthat depend on the transition probabilities and, through φ1(t) and φ2(t), on time. Further orders in the expansion modify this Gaussian nature of the fluctuations, but one need not consider them in general. The next order, for example, is of order Ω−1 relative to the macroscopic values, and so of the order of a single animal.
For the purpose of analyzing the oscillations in the fluctuating solutions of the master equation, it is more convenient to pose eq. (5) in the form of an equivalent Langevin equation:
0 1 2 3 4 0.0
0.2 0.4 0.6 0.8 1.0
^ ^ω2 Α ,ω2 B
ε
susceptible infected
0 1 2 3 4
0.1
1 πε2
Q A, Q B
ε
susceptible infected
f(ε)
0.0 0.2 0.4 0.6 0.8 1.0 10-4
10-3 10-2 10-1
δQ
ε
Figure 3. Bounds in the frequency (left) and in the quality (right) of sto- chastic oscillations. Many microscopic systems are shown in each plot, with circles and crosses corresponding to each population.
ξ˙1=C11ξ1+C12ξ2+η1(t), (6a)
ξ˙2=C21ξ1+C22ξ2+η2(t), (6b)
where the noisesη1(2) are Gaussian (since eq. (5) is linear): hηi(t)ηj(t0)i=Dijδ(t− t0), withD11=P
iji2Wij, D12=P
ijijWij, andD22=P
ijj2Wij.
Equations (6) can be Fourier-transformed and an analytic expression for the power spectra can be found. For the first population, for example, it reads
S1(ω) =
· F1+ ˆω2 (1−ωˆ2)2+ ˆω2²2
¸D11
∆ , (7)
where the normalized frequency ˆω2 = ω2/∆ and F1 = (C122 D22 +C22D11 − 2C12C22D12)/(∆D11), are parameters defined by the system. It can be shown easily that S1(ω) either has a single maximum or is monotonically decreasing (in which case it has a single maximum atω= 0). The maximum lies at
ˆ
ω21=−F1+p
(F1+ 1)2−²2F1, (8)
which can be naturally associated with the ‘frequency’ of the stochastic oscillation.
There are corresponding expressions for the second population, and one sees that the two frequencies do not coincide in general. Indeed, it is possible to show that, for any given microscopic model, both frequencies are bounded by
1−²2/2<ωˆ1(2)2 <1. (9)
Figure 3 (left) shows these bounds as a function of ²as continuous lines. Between them lie a multitude of points, which correspond to many microscopic models.
The frequency of the damped oscillations of the single deterministic model that constitute their macroscopic limit is shown as a dashed line. It is seen that, in the microscopic description, the two populations oscillate with different frequencies, which collapse to a single one as²→0.
Now, the following legitimate question arises: does a peak in the power spectrum identifies, in any sensible way, a ‘good stochastic oscillation’ ? Figure 2 shows that, naturally, this is not the case. Both panels show the fluctuating infected population of a simulated SI model. The insets show the power spectra calculated from them, and the form of the peaks confirm what could be said directly from the field observations: that the one on the left ‘oscillates’, while the one on the right has apparently random fluctuations. Both spectra, however, have a peak at a finite frequency. The difference between them lies on the form of this peak, and to this feature we turn now our attention.
We can show that the following definition of thequalityof the oscillation, based on a property of the power spectrum, enables a good characterization of the oscillating regime:
Q1(2)= ω1(2)hS1(2)(ω1(2))i
RhS1(2)(ω)idω , (10)
and a corresponding expression for population 2. This definition of the quality is arbitrary, in the sense that many other definitions of the same concept are possible based on the form of the power spectrum. It is justifieda posteriori because very good bounds can be derived for it, showing its general behavior for a whole range of microscopic models. It can be shown that both qualities (10) obey
f(²)< Q1(2)< 2
π, (11)
wheref(²) has a cumbersome analytical expression that can nevertheless be found without difficulty from Q1(2) considered as a function of F1(2) for fixed ². These bounds are shown in figure 3, together with the qualities calculated for many micro- scopic models which tend to the same deterministic model in the macroscopic limit.
As we see, the existence of the lower bound tells that the quality of the stochastic oscillation can be large only when ²is small. And, because of the existence of an upper bound that so tightly follows the lower one, when it is large, it is so for all microscopic models.
In summary, the definition of an adequate measure of the quality of stochastic oscillations, and the analytical bounds found for it, show that the parameters that define a deterministic model are enough to characterize how oscillatory the time series of any of its stochastic counterparts looks like. We have found that good stochastic oscillations are clearly present only when ² is small (². 0.5, from our survey of known cases). Our bounds provide aquantitative and accurate measure of this effect. For increasing values of², there is a smooth transition to a behavior that is indistinguishable from noise.
Acknowledgements
This work has been partially supported by ANPCyT (PICT-R 02-87/2 and 04-943) and from CONICET (PIP 5114). GA presented the talk at the PNLD conference and his participation in the PNLD conference was possible, thanks to the hospitality of the Abdus Salam ICTP.
References
[1] G Abramsonet al,Bull. Math. Biol.65, 519 (2003) [2] I D Peixoto and G Abramson,Ecology 87, 873 (2006) [3] L Giuggioliet al,Bull. Math. Biol.67, 1135 (2005)
[4] G Abramson, Mathematical modelling of Hantavirus: From the mean field to the individual level, in: Progress in mathematical biology research edited by James T Kelly (Nova Science Publishers, 2008), pp. 219–245
[5] S Risau-Gusm´an and G Abramson,Eur. Phys. J. B60, 515 (2007) [6] E B Wilson and O M Lombard,Pathology31, 367 (1945)
[7] M S Bartlett,J. R. Stat. Soc.A120, 48 (1957)
[8] J G¨u´emez and M A Mat´ıas,Phys. Rev.E48, R2351 (1993) [9] M A Mat´ıas and J G¨u´emez,J. Chem. Phys.103, 1597 (1995) [10] H Wang and Q Li,J. Chem. Phys.108, 7555 (1998)
[11] A J McKane and T J Newman,Phys. Rev.E70, 041902 (2004)
[12] N G van Kampen, Stochastic processes in physics and chemistry (Elsevier Science B.V., Amsterdam, 2003)
