How bad is to be slow-reacting ? On the effect of the delay in response to a changing environment on a population’s survival Ioana Bena∗ and Michel Droz† University of Geneva, Theoretical Physics Department, Quai E. Ansermet no. 24, 1211 Geneva 4, Switzerland Janusz Szwabin´ski‡ University of Geneva, Theoretical Physics Department, Quai E. Ansermet no. 24, 1211 Geneva 4, Switzerland and Institute of Theoretical Physics, University of Wrocl aw, Pl. M. Borna 9, 50-204 Wrocl aw, Poland 8 0 Andrzej Pe¸kalski§ 0 Institute of Theoretical Physics, University of Wrocl aw, Pl. M. Borna 9, 50-204 Wrocl aw, Poland 2 (Dated: February 2, 2008) n a Weconsiderasimple-modelpopulation,whoseindividualsreactwithacertaindelaytotemporal J variations of their habitat. We investigate the impact of such a delayed-answer on the survival 7 chances of the population, both in a periodically changing environment, and in the case of an 1 abruptchangeofit. Itisfoundthatforpopulationwithlowdegreeofmutation-inducedvariability, being“slow-reacting” decreasestheextinctionriskfacetoenvironmentalchanges. Onthecontrary, ] for populations with high mutation amplitude, thedelayed reaction reduces thesurvivalchances. E P PACSnumbers: 87.10.+e,0.2.70.Lq . o bi I. INTRODUCTION was clarified. However, the intrinsic stochasticity, the - dynamically-built correlations between the individuals, q and the role of the mutation-induced variety in popula- [ In the nowadays context of global warming and habi- tion’sevolutioncannotbeappropriatelyaccountedforat tat destruction, there is an enhanced general interest in 1 a mean-field level. the impactof environmentalchangesonbiologicalpopu- v A more refined level of description, which is an indivi- 8 lations evolution. Despite this, little has been done from dual-basedone, wastherefore alsoconsidered. The main 0 a theoretical point of view, in the frame of evolution- conclusions were that the inherent fluctuations do not 7 ary dynamics modeling, towards a systematic approach destroythephasetransition“extinct-alive”,andthemu- 2 of the role of various elements involvedin these complex 1. circumstances on the population dynamics. In a recent tation amplitude strongly influences the value of the critical selection pressure, giving rise, in particular, to 0 paper[1]weinvestigatedsystematicallytheroleofthese- a diversity-induced resonance phenomenon [5, 6]. The 8 lection pressure and mutation amplitude, as well as the 0 impactofthequalityandquantityofthehabitatchanges phase diagram in the plane of the selection and muta- : tion parameters was discussed as a function of the en- v on the behavior of a single-species population. vironmental variation characteristics. In particular, an i For simplicity and in order to extract the generic fea- X important aspect well-known to experimental biologists, tures, we considered the case of a periodically chang- r see e.g. [7], was emerging naturally, namely that a small a ing environment, as in Refs. [2, 3, 4]. The case of an amount of randomness, due to mutations, is beneficial abrupt change in the environment was also addressed. for population’s survival in the changing environment, The mean-field level of description of the chosen model while a too large amount definitely is detrimental to it. allowed us to put the finger, for the first time, on the The differences between a smooth variation of the envi- very origin of the emerging complex behavior of this ronment and an abrupt, catastrophic change were also highly-nonlinearsystem,thatisthedelicateinterplay be- clarified, pointing to the beneficial role of the mutation tween the different time-scale processes. The role of the in ensuring species survival after a catastrophe. amplitude and period of the environmental changes on In this short paper we shall address another aspect of thecriticalvalueoftheselectionpressure(corresponding this survival problem, namely the role of the delay in to a phase-transition “extinct-alive” of the population) the “reactions” of the individuals. The lagged response to environmental changes is a phenomenon widespread in nature [8, 9, 10, 11, 12, 13]. However, an extensive theoretical analysis of its impact on population dynam- ∗Electronicaddress: [email protected] ics is still lacking. The role and effects of time-delay in †Electronicaddress: [email protected] ‡Electronicaddress: [email protected] biologicallysystemshasbeenaddressedpreviouslyinthe §Electronicaddress: [email protected] contextofLotka-Volterratypeofdynamicsofinteracting 2 species [14], where the “delay” was included at the level following expression: of the coupling between the species. Here we are consid- S eringadifferentproblem,namelythedelayed-responseof p (t) =1− exp − , (3) i the individuals ofa single-speciespopulationto a chang- (cid:20) fi(t)(cid:21) ing environment. Using a simple model, we shall try to where S is a parameter which models the selection pres- clarifythedegreeandlimitsofvalidityofthecommonly- sure of the environment and constitutes a main control spred belief that “a population of fast-reacting individ- parameter of the system. During its life-time, and indi- uals has better survival chances face to changes in their vidual oscillates cyclically from being perfectly-adapted, environment”. when z = ϕ(t−τ ), i.e., from a minimum possible ex- i i tinction rate p (t) = exp(−S), to a worse adaptation, i II. MODEL which corresponds to zi 6= ϕ(t−τi) and to a larger in- stantaneous extinction probability, and finally to a total lack of adaptation, when p (t)=1. The pool of adapted We consider the same type of model as in Ref. [1], i individuals changes thus at each time step. namely a population of hermaphrodite individuals (i.e., The choice (3) we made of the extinction probability which,althoughbisexual,needmatingforreproduction), and the implicit definition of the selection pressure pa- living on a two-dimensional square lattice of size L×L. rameter S are frequently encountered in the biological We assumethatthe individuals cannotcrossthe borders literature, see e.g. [15]. Other choices and thus other of the lattice. Moreover, the lattice has a finite carry- ways of measuring the “selection pressure” are of course ing capacity, whichcomesfromtheexclusionassumption possible. However, most of them can be mapped one that there is at most one individual in eachlattice node. onto the other and/or account for equivalent qualitative The dynamics of the population takes place at dis- aspects of the interaction between the individuals and crete time-steps and is the result of: natural selection their environment. (interaction with the environment), individual motion, The individual delayed-response time τ is fixed once mating and reproduction, as described below. i and for all at one individual’s birth. We consider here a simple case when the τ ’s are random variables drawn A.Naturalselection. Individualtrait,time-dependentop- i from an uniform distribution within an interval [0, T ]. timum, fitness, delayed-response, selection pressure, ex- d The upper limit of this interval T represents another tinction probability. d control parameter of the model. It is obvious that for Each individual i is characterized by its trait or pheno- a periodic variation of the environment only the values type, which for simplicity is represented here through a realnumber z ∈ [0,1]. The traitis fixed once and for all Td <2T0 are relevant. i Note also that an equal delay-time for all individuals at the birth of the individual. amounts simply to a change in the time-origin. As The population lives in an environment whose influ- such, a mean-field level description of the population ence on the individuals is encoded in the value of the dynamics (which is already known as inappropriate for so-called optimum, ϕ ∈ [0,1], which we suppose to be describing mutation, see [1]) will not be able to account homogeneous in space, but periodically variable in time for the effects of the individual delay-times on the ϕ = ϕ(t). Moreover, we consider here the simplest pos- global evolution of the population. We shall therefore sibility, focus exclusively on the individual-based numerical t−tinit simulations. ϕ(t)=0.5+Asin 2π Θ(t−tinit). (1) (cid:18) T0 (cid:19) B. Individual motion. Here A denotes the amplitude of the environmental os- An individual can move to its surroundings, and the cillation, with 0 < A 6 0.5, T0 is its period, and tinit is simplest possibility that we shall adopt hereafter is a the moment of onset of the optimum perturbation; Θ is random-walk. Namely, in one time step the individ- the Heaviside step function. ual jumps on the lattice, from its initial location to The case of an abrupt change in the environment, for a randomly chosen nearest-neighbor one (i.e., a site which the optimum jumps at t = tinit from ϕ = 0.5 to within the von Neumann neighborhood of the initial ϕ=0.5+A was also considered. node), provided that the chosen site is empty, and that Anindividuali“reacts”withacertainspecificdelayτi it lies within the boundaries of the system. If none to the changes in the environment. This means that its of the four first-neighbor nodes is empty, then the in- instantaneousfitness(or“adequacytotheenvironment”, dividualcannotmove,andthuscannotmate(seebelow). see below) at time t is determined by the value of the optimum at a previous time (t−τi), C. Mating and reproduction. Heredity and mutation. Suppose an individual i reaches a destination node. If f (t)=1−|z −ϕ(t−τ )|. (2) i i i there are other individuals (“neighbors”) in the nearest- The fitness determines the instantaneous individual ex- neighborhoodofthisdestinationsite,thentheindividual tinction probability per time step p (t), according to the “i” choses at random one of these neighbors, call it “j”, i 3 for mating. The pair of individuals i and j may then IfatthetimetthereareN(t)individualsinthesystem, give birth to offsprings, which are placed at random on then the above steps A–C are repeated N(t) times; this the empty nodes of the joint nearest-neighborhoods of constitutes one MCS, the unit-time of the simulations. thetwoparents(thatcounts6sites);therefore,themax- Afterwards, the time is advanced by one step, t→t+1, imum number of offsprings of one pair of parents equals and the above algorithm is repeated. 6. If there is no room in this neighborhood for putting As a last remark on the model, it is known on general an offspring, then this one is not born. backgrounds[16]thatthesystemsizeisplayingacertain Thetraitofaprogenykcomingfromparentsiandj is role on the location of the phase transitin point, as well determined by the parents’ traits (heredity), but it can as on its “sharpness”. We used for all our Monte Carlo also present some “variations” due to different random simulations a system of 100×100 lattice sites, for which factors, such as recombination, mutations, etc. We shall wehadshownpreviously,seeRef.[1],thatthequalitative assume that featuresofthephasediagramarepracticallynotaffected by finite-size effects. 1 z = (z +z )+m , (4) k i j k 2 III. RESULTS where m represents these variations. It brings diversi- k fication into the phenotypic pool of the population and we call it conventionally mutation. For simplicity, we Foraperiodicoscillationoftheoptimumweinvestigate shall admit that m is a random number, uniformly dis- the temporal evolution of a population starting from a k tributed in the interval [−M, M], where 0 < M < 1 is given initial concentration c(0). Depending on the char- called hereafter the mutation amplitude and is a control acteristic parameters, the population can evolve, on the parameter of the system [17]. Moreover, if Eq. (4) leads average, either to an “alive phase”, for which its con- to z > 1 or z < 0, then one “renormalizes” it by re- centrationisactuallyoscillatingperiodically,withperiod k k setting zk to 1, respectively 0, which means simply that T0/2, around a nonzero mean value, or can get extinct the trait of the individuals cannot overcome some fixed afteratransientperiodoftime. Inourpreviouspaper[1] limits. Thischoice(4)forthetraitofanoffspringisoften we investigated in detail the phase diagram “extinct– made in the biological literature [15]. alive”ofthepopulationintheplaneofthecontrolparam- The population dynamics is thus driven by two main eters S and M, for different values of the characteristics “forces” that are acting, to some extent, in opposite A and T0 of the optimum oscillation. The same type of directions: selection and mutation, characterized, re- phase diagram was also constructed for the case of an spectively, through the values of the control parameters abrupt jump of the the optimum. S and M. Selection, combined with heredity, tries to Our principal concern in this paper is to determine bring the average trait close to the optimum, while mu- howthe delay inthe individualresponse to the changing tation introduces diversity in the individual traits, and environment– i.e., the value ofthe controlparameterT d thus is broadening the distribution of the population’s – affects the phase diagram extinct-alive of the system. traits. We performed extensive simulations for various range of parametersandthe mainresultsareillustratedinFig.1. The Monte-Carlo simulation algorithm considers the One notices several interesting features exhibited by individuals distributed on the lattice nodes, the initial these figures: condition being represented by their positions, the pre- scribed values of the individual traits and delay-times. a) Consider first the “intermediate” values of T0 for The initial N(0) individuals are randomly-distributed which, as described in Ref. [1] for the no-delay case, one with a mean concentration c(0) = N(0)/L2, and their encountersthediversity-inducedresonancephenomenon, individual traits are randomly assigned from an uniform i.e., the “peak” in the phase-diagram illustrated in the distribution between 0 and 1. upper and middle panels of Fig. 1. Then: The individuals are evolving, at discrete Monte-Carlo (i) For small values of the mutation amplitude M, the time steps (MCS, defined hereafter), according to the existence of a delay in the response of the individuals stages A–C of the dynamics as described above, namely: to environmental changes (i.e., T 6= 0) is increasing the d A. Ata giventime t anindividuali is pickedatrandom, survival chances ofthe population. The diversityrelated and its extinction probability p (t), corresponding to to the randomness in the response of the individuals i one MCS, is determined according to Eqs. (2) and (3). can contribute to the appearance of a larger pool of Then a random number r is extracted from an uniform well-adapted individuals and is thus formally equivalent distribution in the range [0, 1]; if r < p , the individual to an increase in the “effective” mutation amplitude, i dies, otherwise it survives. which is beneficial for the survival [1]. B. If it survives, the individual i jumps at random to (ii) For large values of M, however, adding the ran- one of the empty nearest-neighbor nodes on the lattice. domness ofthe delayed-responseto the mutation-related C. Then it possibly mates and produces offsprings. one is leading to an even higher “effective” mutation amplitude. Assuch,theextinctionriskofthepopulation 4 0.22 is increased: as seen in the plots, the phase diagram T = 0 for the populations with delayed-response (Td 6= 0) lies 0.21 Td = 0.25 T always below the one of the instantaneously-reacting d 0 Td = 0.5 T0 population (Td =0). 0.2 T = 0.75 T (iii) Finally, the peak related to the mutation-induced d 0 diversity is generally still present for the systems with S0.19 time delay. However, in this case the randomness in the delayed-response can turn a part of the pool of 0.18 well-adapted individuals into less-adapted ones, and thus the height of the peak is reduced as compared to 0.17 the case of an instantaneously-adapting population. For T0 = 5000 large delays (like Td = 0.75T0 in the figures) this peak 0.160 0.1 0.2 0.3 can be even suppressed. M b) One concludes therefore that the role of the delay-induced diversity is an increase in the “effective” 0.22 mutation amplitude. As such, it can be easily predicted T = 0 that for small values of T0 (rapid oscillations of the 0.21 Td = 0.25 T environment) the dynamics of the system will be only d 0 T = 0.5 T slightly affected by the delay, since it is already only d 0 0.2 T = 0.75 T slightly sensitive to changes in M. This is illustrated d 0 in the lower panel of Fig. 1. No diversity-induced S0.19 peak, i.e., no optimal “effective” mutation amplitude is encountered in these cases, any mutation and any 0.18 delay in response being harmful for the surviving of the population. 0.17 T0 = 1000 A way to get a better insight into the reasons for this 0.160 0.1 0.2 0.3 behavior is the monitoring of the temporal evolution of M the pool of fittest individuals (i.e., the individuals with f (t) = 1). Figure 2 illustrates this point, for a fixed i value of T0 and three values of Td 6= 0, coresponding to 0.22 themiddlepanelofFig.1. TheupperpanelofFig.2per- tains to the region of small mutation amplitudes in the T = 0 0.21 d phase diagram, for which a delayed-response enhances T = 0.25 T d 0 the survival chances. The lower panel refers to the re- T = 0.5 T d 0 0.2 T = 0.75 T gion of the peak in the phase diagram, for which delay d 0 increases the extinction risk. S0.19 One can see that for the surviving populations the number n(f =1) of the instantaneously fittest individu- 0.18 als is oscillating periodically in time (but never reaching zero), while it decays (with oscillations) to zero for the 0.17 populations that will get extinct. The pool of the fittest T = 50 individuals isenhancedby the delay-induceddiversityin 0 0.16 systems with small mutation amplitude (upper panel of 0 0.1 0.2 0.3 M Fig.2)and,onthecontrary,itisdepletedbythedelayed- responseinpopulations with intermediate andlargemu- tation amplitude (lower panel of Fig. 2). FIG.1: Phase diagram extinct(abovethecurve)–alive(be- low the curve) in the plane of the selection pressure S and Finally, we addressed also the effects of a delayed an- mutation amplitudeM,for differentvaluesofthedelay time swer in the case of a catastrophic, abrupt change in the Td and of the optimum oscillation period T0. From the up- environment. As illustrated in Fig. 3, one encounters perto thelower panel, T0=5000,1000, and 50, respectively; the same type of phenomena as in the case of a smooth the values of the other parameters are L = 100, c(0) = 0.7, variationof the optimum, namely the fact that for small A = 0.3, and tinit = 1000. The average was taken over 10 mutation rate the largest the delay parameter Td, the realisations of thestochastic dynamics and the estimated er- bigger the survival chances of the population. rorsinthevalueofthecriticalselectionpressurearelessthan Inorder to understandthe mechanismunderlying this ±0.002. behavior of the populations with small mutation ampli- tude M, it is useful to follow the temporal evolution of 5 300 0.22 T = 0 T = 0 d d T = 0.25 T 0.21 T = 50 d 0 d T = 0.50 T T = 1000 d 0 d 200 Td = 0.75 T0 0.2 Td = 5000 ) 1 = (f S0.19 n 100 0.18 0.17 abrupt jump 0 tinit 2000 3000 4000 5000 0.160 0.1 0.2 0.3 t M FIG. 3: The phase diagram extinct (above thecurve) – alive 600 (belowthecurve)intheplaneoftheselectionpressureS and mutation amplitude M, for an abrupt jump in the value of 500 the optimum, from ϕ = 0.5 to ϕ = 0.8, for different values of thedelay time Td. The values of the other parameters are 400 L=100,c(0)=0.7,tinit=1000. Theaveragewastakenover ) 10 realisations, and the estimated errors in the value of the 1 =300 critical selection pressure are less than ±0.002. (f n 200 sure this persistence of the high-fitted individuals pool 100 for a long enough time, and the population dies, since the adaptation through mutations is not rapid enough. 0 t 2000 3000 4000 5000 As seen in Fig. 3, on the contrary, for large mutation init t amplitudesthelargerthedelayT ,thehighertheextinc- d tion risk, since, as in the case of a periodically-varying environment, in this case the delay-induced stochastic- FIG. 2: Time evolution of thenumbern(f =1) of thefittest ity adds up to the mutation, leading to an even higher individuals of a population for different values of the delay effective mutation amplitude, which is harmful for the time Td. Upper panel: M = 0.025 (small mutation ampli- system. tude), lower panel: M = 0.1 (intermediate mutation am- plitude). The legend in the upper panel also applies to the lower panel.The values of the other parameters are L=100, c(0)=0.7, A=0.3, T0 =1000, and tinit=1000, correspond- IV. CONCLUSIONS ing to themiddle panel of Fig. 1. We considered a simple model of single-species popu- lationdynamicsinachangingenvironmentandweinves- the fitness histogram“numberofindividualsn(f)versus tigated the role of a delayedanswerof the individuals to fitness f”. This is done in Fig. 4 for two populations these habitat changes. In the case of a smooth variation thatdifferonlythroughthevalueofthedelayparameter of the environment, it was found that, in general, for T , such that one of them gets extinct, while the other populations with small mutation amplitudes it is more d one survives after the catastrophe. Before the catastro- beneficial, in terms of the survival chance, to be slow- phe, the histogram had an important peak at f = 1, reactingthantoanswerinstantaneouslytothevariations and a tail (due to mutations) to low-fittnesses. After of the environment. However,for intermediate and large the catastrophe,a new peak of low-fittedindividuals ap- mutation amplitudes, faster reactions are preferable to pears, such that the histogram becomes bimodal. One slower ones. In case of a very-rapidly oscillating envi- notices that the existence of a larger delay time ensures ronment, the rapidity of reactioninfluences only slightly the persistence of a sufficiently large pool of high-fitted the survivalchances. The same type of statements holds individuals evenafter the carastrophe,andthis poolwill true for the case of a catastrophic, abrupt jump in the ensure the survival of the species till the new-born in- optimum. dividuals get adapted slowly, through small mutations, Assuch,onehastobe rathercautiouswith“common- to the new environment. For a surviving population the sense” statements of the kind “a population of fast- histogrambecomespeakedagain,inthe longrun,around reacting individuals has better survival chances face to f = 1. A shorter delay time T , however, cannot en- changesin their environment”. Of course,more complex d 6 200 and realistic models than the one we presented here T = 1000 are needed in order to make more detailed quantitative d T = 50 statements and reliable predictions for real biological d 150 t = 1000 systems, and to investigate further aspects of the intricate problem of a population evolving in a changing environment. n(f)100 50 0 0.6 0.7 0.8 0.9 1 f 200 T = 1000 d T = 50 d 150 t = 1100 n(f)100 50 0 0.6 0.7 0.8 0.9 1 f 200 T = 1000 d T = 50 d 150 t = 2100 n(f)100 50 0 0.6 0.7 0.8 0.9 1 f FIG. 4: Fitness histograms: number of individuals n(f) ver- susfitnessf fortwopopulationswithlowmutationamplitude M=0.05anddifferentdelay-responseparametersTd =1000 (continuousline)andTd =50(dottedline)incaseofacatas- trophic event. The optimum jumps at tinit = 1000, from ϕ = 0.5 to ϕ = 0.8. The upper panel corresponds to time t = 1000, just before the optimum jump. The middle panel corresponds to t = 1100, and the lower panel to t = 2100 M.D. and I.B. acknowledge partial support from the (when only the population with Td = 1000 survived). The Swiss National Science Foundation. M. D. and J. S. ac- other parameters are L = 100 and c(0) = 0.7. Single runs knowledge the COST10-SER-No.C06.0027 program for were considered. support. 7 [1] I. Bena, M. Droz, J. Szwabin´ski, and A. P¸ekalski, Phys. [11] P. J. Weatherhead,Oecologia 144, (2005) 168. Rev.E 76, (2007) 011908. [12] Ch.Both,S.Bouwhuis, C.M.Lessells andM.E.Visser, [2] C.M.Pease,R.Lande,andJ.J.Bull,Ecology70,(1989) Nature 441, (2006) 81. 1657. [13] I. Hanski, Nature396, (1998) 41. [3] R.Lande, Am.Nat. 142, (1993) 911 . [14] One of the first paper on this subject (that counts by [4] K.Sznajd-WeronandA.P¸ekalski,PhysicaA269,(1999) now ahugebibliography) is: P. J. Wangerskyand W.J. 527. Cunningham, Ecology 38, (1957) 136. [5] C. Tessone, C. Mirasso, R. Toral, and J. D. Gunton , [15] R. Bu¨rger and M. Lynch,Evolution 49, (1995) 201. Phys.Rev.Lett. 97, (2006) 193101. [16] Forapedagogicalintroductiontofinite-sizeeffectsinthe [6] L.Gaimmaitoni, P.H¨anggi,P.Jung,andF.Marchesoni, contextofdifferentnonequilibriumphasetransitions,see Rev.Mod. Phys. 70, (1996) 223. R.Toral andC.J. Tessone, Commun.Comput.Phys. 2, [7] M. Zawierta, P. Biecek, W. Waga, and S. Cebrat, Th. (2007) 177. Biosc. 125, (2007) 123. [17] In the biological literature parameters analogous to M [8] P. M. Thompson and J. C. Ollason, Nature 413, (2001) are often referred to as mutation rate. However, because 417. of the physical aspect M designates, the term mutation [9] Ch. Both and M. E. Visser, Nature 411, (2001) 296. amplitude looks more appropriate tous. [10] A.D.AndersandE.Post,J.Anim.Ecol.75,(2006)221.