Reaction-diffusion wave fronts : Multigeneration biological species under climate change

A generalization of reaction-diffusion models to multigeneration biological species is presented. It is based on more complex random walks than those in previous approaches. The new model is developed analytically up to infinite order. Our predictions for the speed agree to experimental data for several butterfly species better than existing models. The predicted dependence for the speed on the number of generations per year allows us to explain the change in speed observed for a specific invasion.


I. INTRODUCTION
Many attempts have been made in order to describe biological migrations and colonizations by physical methods ͓1-4͔.A possible approach to these problems is based on reaction-diffusion equations ͓4͔.Extension of such equations by considering time-delayed processes has focused the attention of many physicists in recent years ͓5-7͔.In this paper, we expand this approach to allow the study of multigeneration species, i.e., to explain range expansions for species that have several generations per year ͑gen/yr͒, separated by different delay ͑or resting͒ times.The problem we want to solve is fundamentally different from that of a waiting time distribution function considered by other authors ͓8,9͔.In their models, particles or individuals may ''jump'' after a rest time 1 with probability p 1 , after a rest 2 with probability p 2 , etc., and this happens at any instant of time.In contrast, in the case we shall introduce below, there is a seasonal, nonoverlapping succession of resting times.
Delayed diffusion-reaction models can be derived from random-walk movements ͓5͔.Every particle or individual is supposed to move at successive steps ͑with time of travel t), separated by a time of rest ͓Fig.1͑a͔͒.Here we allow for more complex situations, by introducing the possibility of secondary steps with different travel and rest times t 1 , 1 ,t 2 , 2 . . .͓Fig.1͑b͔͒.
In Sec.II we derive our model.In Sec.III, we apply it using typical dispersion and reproduction data for British butterflies, some of which present several gen/yr ͓10͔.Such species have been observed to expand their ranges northwards in the past years, and biologists have pointed out climate change as one main reason ͓11-14͔.We use our equation to predict the typical rates of spread, and compare them to experimental data and previous models.There is good agreement between theory and observations.We argue that this shows ͑i͒ the convenience of analytical models such as the one presented and ͑ii͒ that climatic change, on its own, does not explain the observed speed.

II. MODEL
Let n(r ជ ,t) be the density of particles or individuals in a two-dimensional space.If we assume that the particles jump in the way depicted in Fig. 1͑b͒, we may write for the change in particle number in a differential of area dA during a secondary step of duration

͑1͒
where the first term in the right is due to diffusion, and the second one to net reproduction ͑if dealing with a biological species͒.Following Einstein's approach ͓15͔, we write the diffusive term as the number of particles reaching the area dA minus those leaving it during T i , ͓n͑ r ជ ,tϩT i ͒Ϫn͑ r ជ ,t ͔͒ D dAϭ⌽Ϫn͑r ជ ,t ͒dA, ͑2͒ where with i (⌬x,⌬y)dA the fraction of particles which have jumped from an area differential centered at (xϩ⌬x,y ϩ⌬y) at time t into another area differential centered at (x,y) in tϩT i .We also assume, as in ͓15͔, that all dispersion kernels are isotropic, i (⌬x,⌬y)ϭ i (⌬), where ᭝ ϵͱ⌬x 2 ϩ⌬y 2 .Equations ͑1͒ and ͑2͒ can be approximated by Taylor series if the experimental data on the range expansion span along large enough times (tӷ ͚ iϭ1 N T i ) and distances (xӷ⌬x, yӷ⌬y).Following the same approach as in Ref. ͓16͔ we arrive at the expression up to kth order *Email address: daniel.campos@uab.es† Email address: joaquim.fort@udg.es‡ Email address: enric.llebot@uab.es

͑4͒
where a represents the initial growth rate ͑i.e., ͓‫ץ‬n/‫ץ‬t͔ R Ӎan for nӍ0 ͓16͔͒ and D i is the diffusion coefficient of the ith substep, As usual in this kind of analysis ͓4,16͔, we assume the existence of wave front solutions, by using in Eq. ͑4͒ solutions with the form nϷexp"(xϪvt)…, with Ͻ0.Next, to analyze the process during a time interval encompassing N secondary steps ͓i.e., for a time interval we just have to write the expression resulting from Eq. ͑4͒ for the time interval (t,tϩT 1 ), for the interval (t ϩT 1 ,tϩT 1 ϩT 2 ), etc., and add up all these equations.This yields which in the limit k→ϱ acquires the form

͑7͒
The wave front speed v can be found numerically from this equation in the usual way: for given parameter values, the speed is the minimum value of v such that a solution Ͻ0 exists ͓16͔.
In order to show the generality of Eq. ͑7͒, let us take the limit T i →0 ͑weakly delayed systems͒.From Eq. ͑7͒ one finds, up to first and second order in T i , respectively, D* 2 ϩv (1) ϩaϭ0, ͑8͒ In these special cases, using the fact that must be real, we reach the expressions for speed which correspond to the wave front speed for the wellknown parabolic ͑i.e., nondelayed͒ ͓17͔ and hyperbolic ͑i.e., weakly delayed͒ ͓4,18͔ approximations, respectively, except that here T* and D* appear ͑instead of T and D in Ref. ͓4͔͒.This is due to our assumption of secondary steps.The classical, single-substep case ͓4͔ is recovered for Nϭ1.Equation ͑7͒ is more general than Eq.͑12͒, because it applies even to strongly delayed diffusion.Moreover, the new diffusive behavior due to i secondary steps ͓Fig.1͑b͔͒ is considered in the general expression ͑7͒ and taken into account even for the parabolic and hyperbolic approximations by means of the new parameters T* and D*.

III. APPLICATIONS
(1) Invasion speeds.One of the most direct applications of our new model, and indeed its original motivation, is the spread of biological species with N generations per year.Let us see whether our model is able to predict the magnitude of the observed speeds.
Many butterfly species have up to 4 ͑or sometimes even more͒ summer generations, whereas they do not appear in colder seasons.The periods during which individuals are able to disperse, known as flight times, are typically cept for the last annual generation.For this Nth generation, because of the annual behavior (Tϵ ͚ iϭ1 N T i ϭ1 yr), we have N ϭTϪt 1 Ϫ 1 Ϫt 2 Ϫ•••.Thus, if Nϭ4 gen/yr, 4 Ӎ215 days ͑choosing similar values would not change considerably our final results͒.
The diffusion coefficient D i ͑5͒ can be rewritten, according to a well-known result of basic diffusion theory ͓19͔, as where the subscript obs means that these values correspond to direct observations.The values of ͗⌬ 2 ͘ obs and t obs are determined in ecology from mark-recapture experiments, which yield 2ϽD i Ͻ31 km 2 /yr for butterflies ͓20,21͔, depending on the species considered and the ecological conditions.Such a range may seem rather wide, but in fact, in biological systems the range of diffusion coefficients observed spans many orders of magnitude ͑from 10 Ϫ5 to 10 4 km 2 /yr ͓22͔͒.Finally, typical values of a between 0.6 and 1 yr Ϫ1 follow from the population field data as a function of time in Ref. ͓23͔.
We may now compare the predictions experimental observations.Nonmigratory butterfly expansions have been measured for the last thirty years in Great Britain and the US ͓23͔.Their observed spread rates are in the range 1.3-8.6 km/yr ͓24͔.Introducing the values above for t i ,T i ,a,D i into our new Eq.͑7͒, we obtain the prediction ͑for 4 gen/yr͒ 1.9 ϽvϽ8.4 km/yr, which agrees quite well with the observed range above.
The usefulness of assuming several generations can be seen by comparing with the results of a model ͓16͔ based in the same Eq.͑7͒, but with Nϭ1 gen/yr.In that model we take tϭ͚ iϭ1 4 t i Ӎ120 days, so the annual flight time is the same in both cases.The range predicted by the previously known model is 3.3ϽvϽ14.2,which is less consistent with the experimental data ͓24͔.This is also observed in Fig. 2, where our new model ͑solid lines͒ is more conistent with the observed values of a and D i than the model with N ϭ1 gen/yr ͑dotted lines͒.
The inset in Fig. 2 shows the convenience of assuming different waiting times from variations in v/v o (v o is the speed value for 1 gen/yr͒ as a function of N. Triangles correspond to the case with identical waiting times T 1 ϭT 2 ϭ•••ϭT i ͑that is, T*ϭT i and D*ϭD i ).Squares correspond to experimental values for t i , i above.The differences between both models are qualitatively important, especially for low N, simply because in our case ͑squares͒ generations are concentrated in a specific season ͑summer͒.When N is high enough, generations become uniformly distributed throughout the year, and both models lead to similar results, as they should.Of course, for other biological species, differences between both models could be important up to much higher values of N.
(2) Changes of speed.One of the very few butterfly species whose spreads have been measured in detail is Pieris Rapae.Andow et al. ͓25͔ analyzed its southward invasion across the US in the 19th century.They mentioned that the number of gen/yr for this species is higher in Missouri ͑6-7 gen/yr͒ than in more Northern locations ͑typically 3 gen/yr͒, all generations with similar flight times (t 1 Ӎt 2 Ӎ•••).They even argued that this could explain the increase in the front speed observed during this expansion.This gives us a unique oportunity to check whether the tendency of the speed on the number of generations agrees with our model or not.We checked the original data from Scudder ͓26͔ before P. Rapae reached Missouri (N BM Ӎ3 gen/yr) and after passing Missouri (N AM Ӎ6 -7 gen/yr) and found that the observed ratio of both speeds was v AM /v BM ϭ1.5Ϯ0.2.We will compute the predicted speed ratio, instead of both speeds, because P. Rapae is an extremely migrating species so its diffusion coefficient is very difficult to measure and the value of a for this species is also unknown.Equation ͑7͒ does not allow us to perform speed ratios, and we use the aproximation ͑11͒.If we recall the annual behavior (͚ iϭ1 N T i ϵTϭ1 yr) and assume, as mentioned, the same flight times t i for all generations i, Eqs.͑10͒ and ͑13͒ yield D*ϭNϽ⌬ obs 2 Ͼt i / (4t obs •1 yr).Then, Finally, we predict from Eq. ͑14͒ v AM /v BM ϭ1.6, which is very similar to the experimental value above.Additional agreement and more detailed comparisons to observations are expected to be achieved in the future, as new data for specific species will be obtained.

IV. DISCUSSION
In this work we have extended reaction-diffusion models to consider, up to infinite order, a more general diffusive pattern than previous authors.The new random walk is characterized by a sequential, periodic succession of move and rest times ͓Fig.1͑b͔͒.The need for this more general behavior comes from biological applications, namely from the fact that some species present more than one generation per year.
The situation considered can be regarded as a rapidly time-varying value for the rest time.In principle, this could lead to a rapidly changing wave front speed.But the experimental data are not accurate enough to show such an effect.This is why we have presented an averaged description, leading to a constant speed.
We stress that we have presented a specific application for illustration purposes, and because of the interest that phenological response have as an indicator of climate change.Specifically, nonmigratory butterflies have been intensively studied with the aim to predict the effects of climate, weather and habitat variations ͓11,14,27,28͔.Numerical methods have recently been used to describe the polewards expansions of butterflies and other biological species ͓29,13͔.These approaches are able to predict whether the species can be expected to be present in a given locality, but such a detailed, microscopic method assumes in fact that butterfly range shifts can be as high as allowed by the speed of climate change, which has been criticized ͓30͔.A very recent paper ͓14͔ summarizes a list of species which are expanding their ranges in response to climate change.The speeds listed there have very different values.Thus, it is unlikely that these expansions can be explained just by the speed of climatic isotherms ͓11,31͔.
This argument makes us think that analytical approaches to these processes, as proposed here, could be useful.The results obtained by our model are quite satisfactory: predicted ͑1.9-8.4 km/yr͒ and observed ͑1.3-8.6 km/yr͒ ranges for butterflies are very similar ͑see Fig. 2͒.Hence, we consider that our model may be applied to the study of other multigeneration species, as well as to new applications which are in accordance with the diffusive behavior in Fig. 1͑b͒.

FIG. 1 .
FIG. 1. Two kinds of random walks for the trajectory of a particle ͑or individual͒.͑a͒ Classical diffusion, where the travel time ͑t͒ alternates with the rest time () ͓4͔.͑b͒ The more general case in which there is a periodic distribution of travel (t i ) and rest ( i ) times.

FIG. 2 .
FIG. 2. Hatched regions define the range of independently observed values for the parameters a and D i .The solid lines correspond to the new multigeneration model ͑4 gen/yr͒ with vϭ1.3 and 8.6 km/yr, which is the experimental range, T 1 ϭT 2 ϭT 3 ϭ40 days, and T 4 ϭ245 days.Dotted lines are for the model with N ϭ1 gen/yr ͓16͔ and flight time tϭ ͚ iϭ1 4 t i .Inset: Speed ratio v/v 0 vs number of generations N for the multigeneration model ͑squares, T 1 ϭT 2 ϭ•••ϭT NϪ1 ϭ40 days, T N ϭ1 yrϪ ͚ iϭ1 NϪ1 T i days͒ and assuming all T i are equal ͑triangles, T 1 ϭT 2 ϭ•••ϭ1 yr/N).