Time-delayed Reaction-diffusion Fronts

A time-delayed second-order approximation for the front speed in reaction-dispersion systems was obtained by Fort and Méndez ͓Phys. Rev. Lett. 82, 867 ͑1999͔͒. Here we show that taking proper care of the effect of the time delay on the reactive process yields a different evolution equation and, therefore, an alternate equation for the front speed. We apply the new equation to the Neolithic transition. For this application the new equation yields speeds about 10% slower than the previous one.


I. INTRODUCTION
Reaction-diffusion systems have been applied to many complex biological and physical systems such as population dispersals ͓1͔, viral infections ͓2͔, chemical reaction processes ͓3͔, combustion flames ͓4͔, etc.In Ref. ͓5͔ a timedelayed model for the front speed was presented including terms up to second order.However, here we will show that there was an error in the mathematical derivation, and we will derive and analyze the behavior of the correct timedelayed equation for the front speed.
In biological systems, variations in the population number density, p, are due to two processes: population growth ͑reproduction minus deaths͒ and migration ͑dispersion͒.The variation due to population growth can be expressed as a Taylor series, ͓p͑x,y,t + T͒ − p͑x,y,t͔͒ g = Tͯ ‫ץ‬p ‫ץ‬t where the subindex g denotes growth, we have introduced the growth function as F͑p͒ = ‫ץ‬p ‫ץ‬t ͉ g , and T is the time delay ͑one generation in most applications ͓5͔͒.As usual, we assume that F͑p͒ Ͼ 0.
On the other hand, for the migration ͑dispersion͒ we will define the dispersion kernel ͑⌬ x , ⌬ y ͒ which gives the probability per unit area that an individual initially placed at ͑x + ⌬ x , y + ⌬ y ͒ has moved to ͑x , y͒ after a time interval T. Thus, the variation in population number density due to migration can be expressed as ͓5͔ ͓p͑x,y,t + T͒ − p͑x,y,t͔͒ m = ͵͵ p͑x + ⌬ x ,y + ⌬ y ,t͒͑⌬ x ,⌬ y ͒d⌬ x d⌬ y − p͑x,y,t͒.͑2͒ In a system involving the two processes ͑population growth and migration͒, the total variation in population density during a time interval T can be expressed as the sum of both contributions, p͑x,y,t + T͒ − p͑x,y,t͒ = ͵͵ p͑x + ⌬ x ,y + ⌬ y ,t͒͑⌬ x ,⌬ y ͒d⌬ x d⌬ y − p͑x,y,t͒ We assume that the kernel is isotropic, i.e., ͑⌬ x , ⌬ y ͒ = ͑⌬͒, with ⌬ = ͱ ⌬ x 2 + ⌬ y 2 , and we Taylor expand Eq. ͑3͒ up to second order in time and space, thus, obtaining the following reaction-diffusion equation: where 2T .Since F͑p͒ depends only on the population density p, then the last term in Eq. ͑4͒ can be written as In addition, as the density at the leading edge of the front is low, p Ϸ 0, we have that F͑p͒Ϸ pFЈ͑0͒ and FЈ͑p͒ϷFЈ͑0͒.Therefore, for p Ϸ 0 Eq.͑4͒ may be rewritten as

͑6͒
We now assume that for t → ϱ and r → ϱ the front can be considered locally planar.Thus, choosing the x axis parallel to the local speed of the front, c ϵ͉c x ͉, we can look for constant-shape solutions with the form p = p ¯exp͓͑x − ct͔͒.Applying this ansatz to Eq. ͑6͒ we see that the value of can be obtained from As has to be real, we obtain a lower bound for the front speed However, the result obtained in Ref. ͓5͔ was the so-called HRD speed, namely, which is different from Eq. ͑8͒.
The reason for this difference is the following.Here we have used Eq.͑5͒, which allows us to rewrite Eq. ͑4͒ as In contrast, in Ref. ͓5͔ the following equation was used: We can see that the last term is different.The reason is that in Ref. ͓5͔ the subindex g was omitted in the last term in Eq. ͑4͒.Therefore, in Ref. ͓5͔, the last term in Eq. ͑4͒ was not written as in Eq. ͑5͒ but as follows: and thus leading to Eq. ͑11͒ instead of Eq. ͑10͒.This is why in Ref. ͓5͔, speed ͑9͒ was obtained instead of Eq. ͑8͒.However, the derivation above clearly shows that Eq. ͑8͒ is the right result.In this Brief Report, we will apply variational analysis and show that Eq. ͑8͒ is not only a lower bound but the exact speed ͑Sec.II͒.We will also analyze the difference between the new Eq.͑8͒ and the HRD speed ͑9͒ by applying both equations to the Neolithic transition ͑Sec.III͒.In Sec.
IV we present our conclusions.

II. VARIATIONAL ANALYSIS: UPPER BOUND
Equation ͑8͒ is just a lower bound for the speed of front solutions to the new differential equation ͑4͒ ͓or Eq. ͑10͔͒.In order to find an upper bound, we apply variational analysis ͓6͔ to Eq. ͑10͒.As mentioned above, we assume that the fronts have a profile p͑z͒ = p͑x − ct͒ traveling with a speed c Ͼ 0, so all of the derivatives in Eq. ͑10͒ can be expressed in terms of z.We also assume that the population number density p Ͼ 0 cannot attain values above some value p max , the so-called saturation density.Then, defining n͑p͒ =−p z and assuming that n͑0͒ = n͑p max ͒ = 0 and n Ͼ 0 in ͑0, p max ͒, differential equation ͑10͒ can be rewritten as Now, introducing an arbitrary function g͑p͒ such that g͑p͒ Ͼ 0 and h͑p͒ =−gЈ͑p͒ Ͼ 0, we multiply Eq. ͑13͒ by g͑p͒ / n͑p͒.Integrating the resulting expression by parts, we obtain

͑14͒
Now, we can eliminate n͑p͒ from Eq. ͑14͒ applying that for any positive numbers r and s, it follows from ͑r − s͒ 2 Ն 0 that ͑r + s͒ Ն 2 ͱ rs.Let us assume that the condition holds for all p ͑0, p max ͒.As g͑p͒, h͑p͒, n͑p͒, F͑p͒, and ͑D − T 2 c 2 ͒ are positive ͓7͔, we may choose r ϵ͑D − T 2 c 2 ͒hn and s ϵ g n F͑1+ T 2 FЈ͒ into ͑r + s͒ Ն 2 ͱ rs and use Eq.͑14͒ to get the following restriction: Following the method in Ref. ͓8͔, Sec.3.3, it is easy to show that there is a function g for which the equality holds.Then,

͑17͒
In order to obtain the upper bound for the front speed we will use Jensen's inequality ͓9͔ ͐ 0 where ͑p͒ Ͼ 0 and ␣͑p͒ Ն 0. We define ͑p͒ϵg͑p͒ and ␣͑p͒ϵ͕h͑p͒F͑p͓͒1+ T 2 FЈ͑p͔͖͒ / g͑p͒.Using these functions into Jensen's inequality ͑18͒, and applying the result to Eq. ͑17͒, we obtain that

͑19͒
We want an upper bound independent of g͑p͒, so we will first find an expression in which h͑p͒ =−gЈ͑p͒ no longer appears by integrating by parts the numerator in the right-hand side of Eq. ͑19͒,

͑20͒
where we have assumed that F͑0͒ = F͑p max ͒ =0 ͑this holds, for example, for the logistic growth considered in Sec.III͒.Moreover, from Eq. ͑20͒ we obviously have

͑22͒
Let us assume that the population growth function F͑p͒ is a continuous function with FЉ͑p͒ Յ 0 and F͑0͒ =0 ͑again these assumptions are true for the logistic growth considered in Sec.III͒.Then FЈ͑p͒ is a decreasing function for increasing values of p.Its maximum value is reached for p =0.Thus, using the value p = 0 in Eq. ͑22͒ we obtain that the upper bound for the front speed is As the lower bound given by Eq. ͑8͒ is the same as upper bound ͑23͒, we can predict the speed of front solutions to Eq. ͑10͒ without any uncertainty, In contrast, for the HRD Eq. ͑11͒, the exact speed was previously shown to be ͓5͔ . ͑25͒

III. APPLICATION TO THE NEOLITHIC TRANSITION
In order to compare the predictions from Eqs. ͑24͒ and ͑25͒, we will apply them to the spread of the Neolithic transition in Europe, because this is the case to which Eq. ͑25͒ was initially applied ͓5͔.The Neolithic transition is the change from hunter-gatherer to farming economics.In Europe, it took place as an invasion of agricultural populations from the Southeast, which spread across Europe from 13 000 to 5000 years before present ͓10͔.
In order to make quantitative predictions we will use the logistic growth function, which has been widely applied to human populations ͓5,11͔: where a is called the initial growth rate and p max is the saturation density.Using logistic function ͑26͒, Eq. ͑24͒ can be rewritten as whereas the HRD speed ͑25͒, used in Ref. ͓5͔, is Both equations for the front speed depend on three parameters: the initial growth rate, a, the diffusion coefficient, D = ͗⌬ 2 ͘ 4T , and the generation time, T. We will use the ranges a = 0.028Ϯ 0.005 yr −1 ͓12͔, ͗⌬ 2 ͘ = 900-2200 km 2 ͓10͔, and the characteristic value T =32 yr ͓13͔, which have been measured for preindustrial farming populations.For these ranges, condition ͑15͒ is fulfilled, so Eq.͑27͒ gives the speed of fronts.
Figure 1 shows the front speeds obtained from Eqs. ͑27͒ and ͑28͒ for a characteristic mobility value ͗⌬ 2 ͘ = 1531 km 2 .We can see that, for the range of values for the initial growth rate a appropriate to this application, the new Eq.͑27͒ yields slower speeds than Eq.͑28͒ ͑about 8% slower͒.However, this is not the case for all values of a, as can be seen from the inset graph in Fig. 1.
In order to check the validity of Eq. ͑27͒, we have also numerically integrated Eq. ͑10͒, with F͑p͒ given by Eq. ͑26͒, and initially p = p max in a finite region and p = 0 elsewhere.FIG. 1. Comparative plot between the front speed for Eq.͑27͒ ͑solid line͒ and Eq.͑28͒ ͑dashed line͒.The symbols correspond to the speed obtained from numerically integrating Eq. ͑10͒, with F͑p͒ given by Eq. ͑26͒.All results have been calculated for a characteristic mobility value ͗⌬ 2 ͘ = 1531 km 2 .
The speed obtained from the numerical integrations corresponds to the circles in Fig. 1.They agree with the new Eq.͑27͒ within less than 0.8%.
The range of speeds for the Neolithic transition front obtained from archeological data is 0.6-1.3km/yr ͓10͔.We can see in Fig. 1 that the results from Eq. ͑27͒ lie within this range.
To what extent does our result depend on the uncertainty in the value of the mobility?In Fig. 2, we consider the front speed values 0.6, 0.95, and 1.3 km/yr, corresponding to the range obtained from archeological data, for Eq.͑27͒ ͑full lines͒ and Eq.͑28͒ ͑dashed lines͒ ͓14͔.It is seen that the predictions of the new model ͑full lines͒ are consistent with the observed front speed for most of the values of the mobility appropriate to this system.

IV. CONCLUDING REMARKS
In this Brief Report we have improved the derivation of the HRD speed in Ref. ͓5͔.We have obtained the correct evolution Eq. ͑10͒ and the new Eq.͑24͒ for the front speed.
We have applied the new Eq.͑24͒ to the Neolithic transition.Using realistic parameters the front speeds are consistent with the observed range for the Neolithic transition in Europe ͑0.6-1.3 km/yr ͓10͔͒.Comparing these results with those from the HRD speed, we see that for the Neolithic transition our new equation leads to slower speeds.
In this case, the correction obtained is only about 10%, but it could be higher in other systems where generalizations of Eq. ͑24͒ can be useful.For example, our framework could be applied in order to improve the predicted speeds of viral infection fronts ͓2͔.