Table Of ContentFluctuations and correlations in population models with age structure
Martin Howard1 and R.K.P. Zia2
1 Department of Physics, Simon Fraser University, Burnaby, British Columbia, Canada V5A 1S6
2 Physics Department, Virginia Polytechnic Institute and State University, Blacksburg, Virginia 24061-0435
(February 1, 2008)
0
0
Westudythepopulationprofileinasimplediscretetimemodelofpopulationdynamics. Ourmodel,
0
which is closely related to certain “bit–string” models of evolution, incorporates competition for
2
resourcesviaapopulationdependentdeathprobability,aswellasavariablereproductionprobability
n foreachindividualasafunctionofage. Wefirstsolveforthesteady–stateofthemodelinmeanfield
a theory, before developing analytic techniques to compute Gaussian fluctuation corrections around
J
themeanfieldfixedpoint. Ourcomputationsarefoundtobeingood agreementwithMonte–Carlo
2 simulations. Finally we discuss how similar methods may be applied to fluctuations in continuous
1
time population models.
]
h PACS numbers: 87.23.-n, 02.50.-r, 05.40.-a, 87.23.Cc
c
e The problem of population dynamics has attracted ently non–linear. For example, n–dependent Verhulst
m enormous interest over many years (for some introduc- factors are often used to mimic competition for finite re-
tions and recent applications see Refs. [1–8]). Beginning sources.
-
t with simple logistic growthmodels [1], a tremendous va- However one deficiency of the models discussed hith-
a
t riety of systems have been studied displaying a diverse erto is that they are deterministic. Real population sys-
s
. range of behavior, varying from stable fixed points to tems areofcourseaffectedby randomfluctuations,com-
t
a strange attractors. Broadly speaking these models split ing from the environment and/or from the intrinsic dy-
m naturallyintotwocategories: thoseusingcontinuousand namics of the birth/death processes. Such stochastic
- those using discrete time. The simplest discrete time Leslie models have also been investigated[3,13,14],how-
d
models describe species where there is no overlap be- ever only for cases where the birth/death probabilities
n
tween successive generations, leading to difference equa- were independent of the population vector n. To date
o
c tions of the form N(t + 1) = g[N(t)], where N(t) is noinformationhasbeenavailableregardingthemorere-
[ the total population at time t. However these models alistic case of fluctuations in stochastic age structured
may easily be generalized to species with multiple dis- models with population dependent birth/death probabili-
1
v crete agegenerations (for example: eggs,larvae,adults), ties[14]. Thisisthesituationwewillstudyinthisletter.
2 where one or more generations may be present simul- Population models have also been intensively stud-
7 taneously. Instead of a single variable N, information ied by physicists in recent years in the context of so–
1
about the age distribution is now carried in a “vector” called“bit–string”modelsofevolution[6]. Thesemodels
1
n(t) n (t),n (t),...,n (t) , where n (t) is the num- are based on the mutation accumulation hypothesis [5],
0 0 1 D a
≡{ }
0 ber of individuals of age a at time t. Note that D is the which assumes that during the aging process each in-
0 maximum age, while n stands for the number of “new- dividual accumulates exclusively late–acting deleterious
0
/ borns”. Beginning with the pioneering work of Leslie [9] genetic mutations. In “bit–string” models the genome
t
a
models of this type have been extensively analyzed. The of a particular species is encoded as a series of ‘0’s and
m
simplestLesliemodelislinearinn,sothattheevolution ‘1’s (deleterious mutations), and as an individual ages
- equation is just n(t+1) = An(t). Here, A is the Leslie the bits (genes) are activated one by one. When the
d
n matrix accumulated sum of bad genes reaches a certain num-
o ber the individual dies (although death may also occur
c f0 f1 ... fD atyoungeragesdue to Verhulstcompetition). Note that
: v 0 ... 0
v 0 the“bit–strings”oftheoffspringmaydifferfromthoseof
i A= 0 v1 0 , (1) the parent due to additional beneficial (‘1’ ‘0’) or dele-
X ... ... ... terious(‘0’ ‘1’)mutations. Thistypeofm→odelisclearly
r →
a 0 ... 0 vD−1 0 well suited to efficient computer simulation [6,7]. In its
simplestcase,whereindividualsdieafterthefirstdelete-
wheretheelementsf arethefecundities (number ofoff- riousmutation,“bit–string”modelssimplycorrespondto
a
spring produced) of individuals of age a, and the v are multiple genome population models with age structure,
a
Verhulstfactors (the fractionofindividuals of agea who where the different genomes can be distinguished by dif-
surviveto become agea+1). Thoughthe originalLeslie ferent maximum ages. Deterministic versions of some of
model had all the f and v as constants, general- these models have already been treated analytically [7].
a a
izations have since{be}en ma{de }to n–dependent factors However,ananalysisoftheimportantroleplayedbyfluc-
[2,10–12],sothatthe evolutiondynamicsbecomesinher- tuations has so far been lacking. Our calculations form
1
the first step towards filling this gap. from old age. Let us emphasize that this equation is ac-
We begin our analysis by defining our discrete time tually quite complex, since V is a function of the total
population model. For simplicity, we consider only a population. Multiplying (2) by n and summing over all
c
single species reproducing asexually, without mutations. the other indices, we obtain
Thus, our system can be described by a single vector
n. At each time step we compute the Verhulst factor hnait+1 =hV(N) na−1it, (a>0), (3)
V(n) and let eachindividual survive with probabilityV. n = F r V(N) n , (4)
h 0it+1 d dh dit
After this “pruning”, each of the remaining individuals
d
X
of age a may give birth to Fa offspring with probability where denotes the average of over P(n,t). These
r . At this point the remaining population is aged by h•it •
a equations are exact. However, due to the presence of N
one time step, with the exception of the new offspring
through V, all moments of P may be coupled together.
who make up n . Individuals who exceed the maximum
0 Themeanfield(MF)approximationconsistsofreplacing
age D then die immediately and are removed from the
thehigherordermomentsbyappropriateproductsofthe
system. Since our model allows for a variable reproduc-
first moment. Hence we find
tion probability r as a function of age, features such as
a
puberty and menopause can be naturally incorporated.
However, we do assume that reproductive individuals of hnaiMt+F1 ="V( hnciMt F)#hna−1iMt F, (a>0), (5)
c
the same age produce an identical number (F ) of off- X
a
springs. Note that we are not restricting ourselves to n MF = V( n MF) F r n MF, (6)
specific forms for V, ra, or Fa beyond some general fea- h 0it+1 " c h cit # d d dh dit
tures, so that this analysis may easily be applied to the X X
where, to be clear, we have written the explicit expres-
various “bit–string” models [6,7]. These generalfeatures
sion for N. These non-linear equations are known to
include V,r [0,1], since they represent probabilities.
a
∈ contain a rich variety of behavior, depending on the de-
Furthermore, we assume that V depends only on the to-
tails of F ,r , and V. For example, if F r < 1, the
tal population N ≡ ana, via the ratio N/N0, where reproductcivecrates are too low and the pcopculcation will
N represents a characteristic population size that the
0 P
P eventually die out. On the other hand, if the reproduc-
resources can support. Also, to be reasonable, V is as-
tive rates are large enough, the population will display
sumedto be monotonicallydecreasingwith N. Forcom-
period doubling bifurcations and chaos [15]. Let us fo-
parisons with simulations, we use V = s exp( N/N ),
0 0
− cus on the “moderate” range, so that a simple non–zero
wheres isanotherconstant,aformfrequentlychosenin
0
steady–stateexists. Inthat casethese equationsareeas-
the biology literature [1]. In contrast,the algebraicform
ily solved to give
V =1 (N/N )isthefavoriteintherecentphysicslitera-
0
−
ture[6,7]. Weprefertheexponentialform,sinceabsolute za(1 z)
n MF =N(z) − , (7)
cut-offs seem unrealistic in a real population system. h ai (1 zD+1)
As mentioned above, deterministic versions of this −
where z is the unique, positive, real root of the equation
model have been studied [10–12]; in particular the pop-
F r zc+1 =1,andN(z)isthesteadystatetotalpop-
ulation dynamics of a semivoltine species studied in c c c
ulation, given by the value that satisfies V(N)=z. Our
Ref.[12]isquitesimilartoaD =1versionofourmodel.
P
objectiveisto gobeyondsuchwell–knownmeanfieldso-
Our goal is to go beyond these deterministic treatments
lutionsandinvestigatefluctuationsandcorrelations,i.e.,
andanalyze the fluctuations andcorrelationsin this sys-
tem. Therefore, we need to consider P(n,t), the prob- the second moments of P. Thus, we multiply Eq. (2) by
n n andsumoveralltheotherindices. Fora>0,b>0,
ability of finding the population with a particular dis- a b
tribution n at time t. Its evolution obeys the master we have
equation hnanbit+1 = V2na−1nb−1 t+δabhV(1−V)na−1it, (8)
P(n,t+1)= P(m,t)× hnan0it+1 =(cid:10) Fcrc V2n(cid:11)a−1nc t+
m0,...,XmD,nD+1 Xc (cid:10) (cid:11)
D+1 ma−1 Vna[1 V]ma−1−na δ n D F b +Fa−1ra−1hV(1−V)na−1it, (9)
"aY=1(cid:18) na (cid:19) − #b0X,...,bD" 0−Xc=0 c c# hn0n0it+1 = FcrcFdrd V2ncnd t+
c,d
D+1 X (cid:10) (cid:11)
×"aY=1(cid:18)bna−a1(cid:19)raba−−11[1−ra−1]na−ba−1#. (2) +Xc (cid:2)Fc2(rchVncit−rc2(cid:10)V2nc(cid:11)t)(cid:3). (10)
Notethatthen isjusta“temporary”variable,which Note that, like Eqs.(3)and (4), these areexact. Assum-
D+1
keeps track of the number of m ’s who survive the Ver- ingN 1,andthatthesystemiswellawayfrom“criti-
D 0
≫
hulst “pruning” so that they can give birth before dying cal”points(e.g.,bifurcationsandthesurvival/extinction
2
transition), it is reasonable to postulate a Gaussian dis- ∂[Vnc] 1 ∂2[Vnc]
S = F r , U = F r N G ,
tribution for P∗(n) with width of O(√N ). Rewriting 0b c c ∂n 0 2 c c∂n ∂n 0 de
∗ 0 c b c,d,e d e
Eqs. (3,4) and (8-10) for , we then have a closed set X X
h•i
of equations for the first and second moments. This ap- where I is the unit matrix and S is the stability matrix
proach should form the first step in a systematic expan- associated with the mean field (zeroth order) stationary
sion of all quantities in decreasing powers of N0. Fur- solution. Note that both S and U need to be evaluated
thermore, in the same spirit, we will let n/N0 assume atzerothorderonly. WithourassumptionsaboutFa,ra,
continuous values. As a check, we will compare the re- and V(N), the eigenvalues of the stability matrix S usu-
sults from this approachwith those from a Monte–Carlo ally lie within the unit circle, implying that our mean
simulation, for a simple case. field solution is stable. However, for sufficiently high
Proceeding, let us write reproductive rates, perturbations with δn n (i.e.,
a a
∝
populations with the same relative age distribution, but
1 (D+1)/2 1
P∗(n)= (11) withwithdifferenttotalnumbersofindividuals)canhave
(cid:18)2πN0(cid:19) √detG × eigenvalues of less than 1. This is the signal of a pe-
−
riod doubling bifurcation, leading to the breakdown of
exp 1 (n n¯ )G−1(n n¯ ) , our Gaussian perturbation expansion (see also below).
× −2N c− c cd d− d
0 Applyingthesameanalysistothesecondmoments,we
c,d
X
obtain, after some lengthy algebra,
where we expect the unknown (to be determined) pa-
rametersn¯ andGto be ofO(N0)andO(1), respectively. Gab SacSbdGcd =Kab , a,b, (15)
Notethatwewillintegratenfrom ratherthan − ∀
−∞→∞ Xc,d
from 0 , a simplification which will introduce only
→ ∞
negligible errors of O(exp[ N ]). Averages can now be where
0
−
computed using n
a
K =δ (1 V) , a,b>0, (16)
ab ab
∗ 1 ∂2f − N0
f(n ) =f(n¯ )+ N G +... . (12) n
h a i a 2Xc,d ∂n¯c∂n¯d 0 cd Ka0 =K0a =Fa−1ra−1(1−V)Na0, a>0, (17)
1
Note that the second term in Eq. (12) is expected to be K = VF2r n (1 Vr ). (18)
00 N b b b − b
suppressed by a factor of O(1/N ) compared to the first 0
0 b
X
termf(n¯ ). HencetherighthandsideofEq.(12)isactu-
a
Again, all quantities need to be evaluated only at the
allyanexpansioninpowersof1/N . Thisorderingallows
0
zeroth order, so that, e.g., V is just z. In compact form
us to setup a systematicperturbationtheory, whichcan
this equation can be written as G SGST = K, which
be pushed to higher orders if desired.
−
may be solved by series
From now on we drop the bars for clarity (n¯
a
′ →
n ). Defining ξ n /N and V dV/dξ, we have
∂Va/∂n = V′(ξ)/≡N acndc∂20V/∂n ∂n≡= V′′(ξ)/N2, in- G=K+SKST +S2K ST 2+...= SnK ST n.
dependecnt of c or d.P0Applying Eq.c(12d) to Eq. (3) g0ives (cid:0) (cid:1) Xn (cid:0) (cid:1)
(19)
′ na−1 ′′
na =Vna−1+V Gca−1+ V Gcd+... .
2N0 SincetheeigenvaluesandeigenvectorsofSareknown,let
c c,d
X X us write S = MEM−1, where E is in Jordan form, with
(13)
the eigenvalues on the diagonal, and M is the matrix
An equation for n can be similarly derived. With the (with its columns) composed of the corresponding right
0
assumptions n O(N ) and G O(1), the latter two eigenvectors. Note that M is not necessarily orthogonal
terms in Eq. (1∼3) repre0sent O(1/∼N ) corrections to the or unitary. If we define G˜ = M−1G MT −1 and K˜ =
0
mean field results, while the remaining (lowest order) M−1K M−1 T, then it is straightforward to show that
(cid:0) (cid:1)
pieces make up the mean field equation (5). Writing a G˜ = EnK˜ En. For simplicity, let us focus on the case
perturbative expansion: na = n(a0) + n(a1) + ... (with where E(cid:0)nis d(cid:1)iagonal. Then the sum is easily performed,
n(k) O(N1−k)), we see that n(0) is given by Eq. (7), so thPat
a ∼ 0 a
while the first order result is
G˜ =K˜ /(1 e e ) (no sum), (20)
n(1) = [I S]−1 U , with (14) ab ab − a b
∂[Vnaa−1] Xc (cid:16) 1− ∂(cid:17)2[aVcnac−1] wwehecraenthdeire{cet}lyaorebttahienetihgeenmvaalturiexs.GSianncdewGit=h iMt aG˜llMthTe,
S = , U = N G , (a>0),
ab ∂n a 2 ∂n ∂n 0 cd information about the Gaussian probability distribution
b c d
c,d
X
3
(11). The explicit formula for computing G is our prin- TABLEI. Comparison of results for the2 age model.
cipalresult. Given a particularformof V(N) andrepro- Gaussian Approximation Simulations
ductive parameters ra,Fa, we can compute G and find hn0i 156.8 156.4
the fluctuations in, as well as the correlations between, hn1i 118.9 118.6
the populations of various ages. hn20i−hn0i2 122.0 130.7
The result (20) contains a further appealing feature: hn0n1i−hn0ihn1i 63.9 67.8
the signal of bifurcation. From stability analysis, hn21i−hn1i2 69.9 71.7
we know that period doubling emerges when the eigen-
value associated with δna na reaches 1. Examining theleadingnoisetermsenteradditively, sothatapertur-
∝ −
Eq. (20), we see that it is precisely this feature which bation theory analogous to the above approach can be
signals the breakdown of the Gaussian approximation. set up.
Furthermore, in many studies of, e.g., the Penna “bit– In conclusion, we have developed analytic techniques
string” model [6], menopause sets in before death, so for dealing with fluctuation effects in a general class of
that E is not diagonal. Then the final expression for population models with age structure. The results we
G˜ will be slightly more complicated, although the above havepresentedalsoformafirststeptowardsanimproved
conclusions will remain qualitatively unchanged. analyticunderstandingofthe“bit–string”modelsofevo-
Tochecktheaboveanalysis,westudythesimplestpos- lution. Finally,itwouldbeinterestingtoperformfurther
siblecase: a2agesystem(i.e. D =1),withra =Fa =1. investigations near the bifurcation point, since interest-
Choosing the exponential form for V with s0 = 1 and ing collective behavior can be expected in that region.
N0 = 100, mean field theory yields n(00) = 157.4 and WethankR.Desai,B.Schmittmann,andU.C.T¨auber
n(0) = 119.3. Performing our analysis, we arrive at the for illuminating discussions, and also R. J. Astalos for
1
first order corrections to n and n , the fluctuations in sharing simulation results prior to publication. This re-
0 1
the populations of eachage, andthe correlationbetween search has been supported in part by the NSERC of
the populations of the two ages. The results are listed CanadaandtheUSNationalScienceFoundationthrough
inTableI,alongsidethosefromMonte–Carlosimulations the Division of Materials Research.
[16]. Theagreementisexcellent,validatingourapproach.
[1] J.D. Murray, Mathematical Biology (Springer-Verlag,
Notethatthecorrectionsn(a1)/n(a0) arelessthan1%,vin- Berlin, 1989).
dicating our assertion that they should be O(1/N0). [2] E.C. Pielou, Mathematical Ecology (John Wiley, New
Uptothispointwehavebeenconsideringmodelswith York, 1977).
discrete time steps. However it is perfectly possible, and [3] J.H.Pollard,ModelsforHumanPopulations(Cambridge
sometimesmoreappropriatebiologically,toanalyzemod- University Press, Cambridge, 1973).
els in continuous time [1,4]. Let us conclude with some [4] B. Charlesworth, Evolution in Age–Structured Popula-
tions (Cambridge University Press, Cambridge, 1994).
brief remarks about fluctuations in this context. A suit-
[5] M. Rose,Evolutionary Biology of Aging(OxfordUniver-
able equation for the mean field population dynamics is
sity Press, New York,1991).
∂n(x,t) ∂n(x,t) D [6] T.J.P. Penna, J. Stat. Phys. 78, 1629 (1995); T.J.P.
= λn(x,t) n(x′,t)dx′, (21) Penna and D. Stauffer, Int. J. Mod. Phys. C 6, 233
∂t − ∂x −
Z0 (1995).
[7] A.F.R. de Toledo Piza, Physica A 242, 195 (1997);
with boundary conditions for birth at x=0 and certain
R.M.C. de Almeida, S. Moss de Oliveira, and T.J.P.
death at x = D. However the birth/death/aging pro-
Penna, Physica A 253, 366 (1998).
cessesgivingrisetoEq.(21)canalsobewrittenasabal-
[8] W. Hwang, P.L. Krapivsky, and S. Redner, Phys. Rev.
listicreactionmodelonadiscretespatiallatticebutwith Lett. 83, 1251 (1999); adap-org/9912004.
continuoustime. Asiswellknown[17],startingfromthe [9] P.H. Leslie, Biometrika 33, 183 (1945); 35, 213 (1948).
correspondingmicroscopic lattice master equation, tech- [10] L.LiuandJ.E.Cohen,J.Math.Biol.25,73(1987);R.A.
niques now exist to map this model onto a field theory Desharnais and L. Liu, J. Anim. Ecol. 56, 885 (1987).
in continuous space–time. The ensuing action can be [11] A. Wikan,Theor. Pop. Biol. 53, 85 (1998).
recast as a Langevin equation, with the result being Eq. [12] R.M. Nisbet and L.C. Onyiah, J. Math. Biol. 32, 329
(21),butwithextramultiplicative noiseterms. Theform (1994).
[13] S. Tuljapurkar, Theor. Pop. Biol. 35, 227 (1989).
of these noise terms would then be completely specified,
[14] S. Engen and B.-E. Sæther,Oikos 83, 345 (1998).
without any ad–hoc guesses. Unfortunately the field–
[15] J. Guckenheimer, G. Oster, and A. Ipaktchi, J. Math.
theoreticactionisratherawkward,duetothepresenceof
Biol. 4, 101 (1977).
non–localinteractionsandnon–local,multiplicativenoise
[16] R.J. Astalos, privatecommunication.
(from fluctuations in the birth process at x = 0). How- [17] L. Peliti, J. Phys. (Paris) 46, 1469 (1984); J. Cardy in
ever,simplificationsoccurifweareinterestedonlyinthe Proceedings of the mathematical beauty of physics, ed.
simple, non–zero steady state, where expansions about J.-B. Zuber, Adv. Ser. in Math. Phys. 24, 113 (1997);
the mean field solution should be adequate. In this case cond-mat/9607163.
4
