Reaction-diffusion waves of advance in the transition to agricultural economics

In a previous paper @J. Fort and V. Me ́ndez, Phys. Rev. Lett. 82, 867 ~1999!#, the possible importance of higher-order terms in a human population wave of advance has been studied. However, only a few such terms were considered. Here we develop a theory including all higher-order terms. Results are in good agreement with the experimental evidence involving the expansion of agriculture in Europe. @S1063-651X~99!19110-4#


I. INTRODUCTION
Allowance of a time delay between cause and effect yields equations that are more reasonable from a conceptual perspective.For example, the Fourier heat conduction equation predicts that a temperature gradient ٌ ជ T causes the instanta- neous appearance of a heat flux q ជ ( is the thermal conduc- tivity, x ជ is the position vector, and t is the time͒.This physi- cally unpleasant property was noted long ago.Authors such as Cattaneo ͓1͔ and Vernotte ͓2͔ proposed to avoid it by letting the heat flux be retarded with respect to the temperature gradient, i.e., using a relationship of the form q ជ ͑ x ជ ,tϩ⌼ ͒ϭϪٌ ជ T͑x ជ ,t ͒, ͑2͒ where ⌼ plays the role of a delay or relaxation time.Such a simple modification leads to generalized heat conduction equations that have been used in the description of second sound in crystals ͓3͔.Similarly, time-delayed equations for viscous flow ͓4͔, diffusion ͓5͔, and heat radiation ͓6͔ have been considered, as well as for electrical ͓7͔ and chemical ͓8͔ systems.Applications include shear waves ͓9͔, ultrasound propagation ͓10͔, shock waves ͓11͔, pores in biological membranes ͓12͔, rheology ͓13͔, etc.It is worth stressing, however, that simple theories based on replacing, e.g., the left-hand side in Eq. ͑2͒ by its first-order Taylor expansion ͓1͔, usually provide only a qualitatively valid description ͓14-16͔.Such approaches lead to the so-called telegrapher equation ͑see Sec.II͒, which has the appealing property that it predicts a finite speed for the propagation of signals ͓14,15,17͔.A special case of time-delayed transport is relaxational diffusion, which has been applied to turbulence ͓5͔, propagation of light in turbid media ͓18͔, diffusion in glassy polymers ͓19͔, photon emission from stellar atmospheres ͓17͔, Taylor dispersion ͓20͔, etc. Again, these approaches are based on a linearization of the type of Eq. ͑3͒ for the diffusion flux J ជ instead of the heat flux q ជ , which ignores additional terms in the expansion.For this reason, it is important to develop more general models.This is one of the purposes of this paper, which has nevertheless been inspired by a specific application that we shall now summarize within its proper context.
In the last few years, a lot of interest has been focused on the application of time-delayed models to systems in which diffusion and reaction processes coexist.Applications include chemically reacting systems ͓21,22͔ as well as many biological applications such as epidemics ͓23͔, forest fire models ͓24͔, and population growth ͓25͔.Most authors have presented formalisms based on simplifications which are essentially of the type of Eq. ͑4͒.This leaves doubt as to the possible importance of the additional, neglected terms.In particular, application of such a model to the expansion of human populations has very recently led ͓26͔ to corrections higher than 40% with respect to the usual, nondelayed model.Since this modification is very large, there is no reason a priori to expect that keeping only a first-order correction in the series ͑4͒ will give quantitatively trustworthy results.It is thus necessary to analyze carefully the role of all higher-order terms, and this is our main purpose here.We will focus our attention on a specific application of the model, namely, the population expansion in the European Neolithic transition, in order to determine whether or not the conclusions in Ref. ͓26͔ remain valid or not when additional terms are included.However, we would like to remark that the formalism we will present here is valid in general, and should be useful in a variety of systems, specially those dealing with time-delayed approaches to reaction-diffusion ͓21-25͔.
The plan of the paper is as follows.In Sec.II, we derive a time-delayed reaction-diffusion equation including terms of up to an arbritrarily high order.Its wave-front solutions are analyzed in Sec.III.This generalizes the theory presented in Ref. ͓26͔.In Sec.IV, we explain why such an equation is a reasonable approach to the modeling of human expansions ͑with special emphasis on the transition to agri-cultural economics͒, we briefly survey how the parameters in the equations are determined experimentally, and find good agreement between the predictions of the new equations and the rate of spread of agricultural communities as determined experimentally from the archeological record.Sec.V is devoted to concluding remarks.

II. GENERALIZED REACTION-DIFFUSION EQUATION
The usual approach to reaction-diffusion is based on the so-called Fisher equation ͓27͔, in fact already derived by Luther in 1906 ͓28͔.This equation can be obtained from Fick's law of diffusion, namely, where D is the diffusion coefficient and p is the particle concentration.This is the diffusion analog to Fourier's heat conduction equation ͑1͒.
where ‫ץ‬ 0 F/‫ץ‬t 0 ϵF.It remains to calculate the migration term.We follow Einstein's approach to diffusion ͓34͔ by letting ⌬x and ⌬y stand for the changes in the position coordinates of a given family during the time interval , and writing the migration term as We replace the right-hand side in Eq. ͑11͒ by its Nth-order Taylor expansion and make use of Eq. ͑13͒, d⌬xd⌬y.͑14͒ In Ref. ͓26͔ the approximation Nϭ2 was analyzed, i.e., only terms up to ‫ץ‬ 2 p ‫ץ‬x‫ץ‬y ⌬x⌬yϩ ‫ץ‬ 2 p ‫ץ‬y 2 ⌬y 2 were considered.Here we include an arbitrary number of such terms, which can be written analogously.After inserting Eq. ͑12͒ into ͑14͒, Eq. ͑10͒ becomes where is the mean square displacement in the x direction during the time interval , etc.One may in principle introduce an infinite set of generalized diffusion coefficients and use them in the terms containing ͗⌬x 4 ͘, ͗⌬x 2 ⌬y 2 ͘, etc. in Eq. ͑15͒.However, this would require the estimation of many parameters, which would complicate or even preclude the comparison of theory to experiment.A much simpler model can be built by assuming that all families move approximately the same distance Ϯ⌬x in the x and y directions during the time interval .Then ͗⌬x k ͘ϭ⌬x k ϭ͗⌬y k ͘ for kϭ2,4, etc.Such lattice models are widely used in biological applications ͓35,29͔, although they have not been previously applied to time-delayed reactiondiffusion.Then, Eq. ͑15͒ becomes Equation ͑16͒ has been derived from the series expansions in Eqs.͑10͒ and ͑14͒.A possible approximation is to include only terms of up to second order, i.e., to neglect time and space derivatives of third and higher order.Then we recover from Eq. ͑16͒ the hyperbolic reaction-diffusion equation de- which is Eq.͑8͒ with a diffusion coefficient and relaxation time given by respectively.We have introduced ⌬ϵͱ⌬x 2 ϩ⌬y 2 .According to Eq. ͑19͒, the relaxation time appearing in the phenomenological equation ͑7͒ is half the mean time between two subsequent migrations.For the reasons explained in Sec.I, this microscopic interpretation is necessary in order to compare theory to experiment in the application considered.The former derivation is valid for an arbitrary system: one needs only to consider the mean time between collisions instead of between migrations.This can also be of interest in the analysis of chemically reacting systems ͓21,22͔.However, such applications are not within the scope of the present paper.

B. Higher-order equations
From Eqs. ͑16͒ and ͑18͒ we find Equation ͑20͒ is the fundamental equation we have been looking for: it generalizes the time-delayed reactiondiffusion equation considered in Refs.͓32,33͔ and ͓23-26͔ by including terms of up to an arbitrary order N.This equation can be used in order to find better solutions than those following from Eq. ͑17͒.In Eq. ͑17͒, only time and spatial derivatives of up to second order were retained.Less approximate results will be obtained by application of Eq. ͑20͒ including spatial and temporal derivatives of up to order N Ͼ2.

III. WAVE-FRONT SOLUTIONS
Wave fronts can be defined as traveling waves with constant shape and speed of propagation ͓29͔.It is observed both numerically and experimentally that, although a continuous range of wave-front speeds is consistent with the stability requirements, the system rapidly evolves toward the minimum possible speed ͓29͔.In the application considered here, propagation of such a wave across a given geographical area describes the immigration and establishment of farming communities.Simple calculations are possible for the generalized reaction-diffusion equation ͑20͒ if we assume that when a sufficiently long time has elapsed from the onset of agriculture, the farmers' wave of advance is approximately planar at scales much larger than that of individual migrations.We may then choose the x axis parallel to the local velocity of the wave.Let vϭ͉v x ͉ stand for its speed (v y ϭ0).We introduce the variable zϭxϪvt and look for constant-shape solutions, i.e., solutions such that p depends only on z.In general, we have F(p)ϭapϩb 2 p 2 ϩ•••, but the migration waves of advance under consideration travel into areas where farming communities were previously absent, so that pϷ0 and thus F(p)Ϸap.It is now easy to rewrite Eq. ͑20͒ as a differential equation involving only derivatives of f with respect to the variable z, where, as explained below Eq. ͑20͒, for a given order of approximation N one should keep temporal and spatial derivatives of up to order N. We are interested in determining the speed of propagation v.A usual method is based on reducing the reactiondiffusion equation to a system of first-order differential equations and finding its eigenvalues.For Nϭ2 the problem thus reduces to a second-order equation ͓25͔.However, this will not hold for NϾ2.Thus we will use a different method: the existence and stability of wave fronts can be studied by considering small perturbations of the form pϭexp͓z͔ about the state pϭ0 ͓29͔.We can require that R, since otherwise we would have an oscillatory behavior with pϽ0 for some values of z, which is a meaningless result.The system evolves toward a stable state provided that Ͻ0 ͓29͔.This method is applied below to four increasingly complicated cases.

A. Fisher's model
Fisher's equation ͑6͒ is recovered from the approximation Nϭ2 ͓Eq.͑17͔͒ in the limit of vanishingly small delay time, →0.We denote the corresponding speed v by v →0 .In this case, use of pϭexp͓z͔ yields the dispersion relation and the requirement R gives Fisher's well known minimal speed v →0 ϭ2ͱaD.͑23͒

B. Simplest time-delayed model "N‫…2؍‬
We denote the corresponding speed by v (2) .Substitution of pϭexp͓z͔ into Eq.͑17͒ leads to and after finding the solutions for we find that stable wave fronts can only exist with speeds equal to or higher than the critical value with aϽ2, which is the expression used in Ref. ͓26͔.It generalizes Fisher's classical result and is in agreement with recent results from the linearization ͓25͔ and path-integral ͓36͔ methods.Here we are interested in determining whether this result can be trusted for application to human population wave fronts.This is the reason we have developed a more general approach that will now be applied.

C. Higher-order solutions
Use of the same form for p as above into Eq.͑21͒ yields which is a polynomial equation.It can be solved analytically, e.g., in the third-order approximation (Nϭ3).However, the results are rather lengthy.One may solve Eq. ͑25͒ numerically for increasing values of N and study their convergence with increasing N. Nevertheless, we prefer to study the exact solution.Its validity will be checked in the next section, by means of the result ͑25͒.

D. Exact solution
In the limit N→ϱ, Eq. ͑25͒ can be written as It could at first sight seem possible to follow an alternative approach based on the diffusive analogue to Eq. ͑1͒, i.e.J ជ (x ជ ,tϩ⌼)ϭϪDٌ ជ p(x ជ ,t) and the mass balance equation.Such an approach would lead to a phenomenological, macroscopic, time-delayed Fisher equation, but not to the microscopic results ͑18͒ and ͑19͒; these results are necessary because D and ⌼ are not directly measurable: in the application we are interested in, what has been derived from experimental observations are the values of ͗⌬ 2 ͘/ and ͑see, e.g., Ref. ͓26͔͒.Moreover, Eq. ͑18͒ cannot be simply borrowed from Fickian diffusion, which holds near equilibrium; here we have shown that Eqs.͑18͒ and ͑19͒ hold arbitrarily far away from equilibrium.In our application ͑see Sec.IV͒, ͗⌬ 2 ͘ is the mean square displacement per generation and is the generation time ͑in chemically reactive systems, they correspond to the mean square displacement and the mean free time between reactive collisions, respectively͒.Thus a microscopic approach, such as that presented in Sec.II, is necessary in order to compare the theory to experiment ͑Sec.IV͒.Consistency with the classical results can be checked by noting that Eq. ͑17͒ has been obtained by dividing Eq. ͑16͒-in the second-order approximation-by .Accordingly, if we divide Eq. ͑26͒ by and consider the limit →0, we recover Fisher's dispersion Eq. ͑22͒, as it should be.
Although Eq. ͑26͒ cannot be solved analytically, we can show that it leads to wave fronts with a finite minimum velocity of propagation.In order to see this, in Fig. 1 we plot the functions on the left-and right-hand sides of Eq. ͑26͒.
For given values of , D, and , the RHS has the shape shown in the figure.Then, to each possible value of the reaction ͑or growth͒ parameter a there will correspond a minimum possible value of the velocity v, since the requirement Ͻ0 implies that the LHS of Eq. ͑26͒ increases with increasing v: the functions on the LHS and RHS will certainly not cross for low enough values of v. Thus for given values of the parameters, a real solution to Eq. ͑26͒ will exist only for speeds above minimum value, in complete analogy to the Fisher and hyperbolic results obtained above by the same method.

IV. APPLICATION TO THE NEOLITHIC TRANSITION IN EUROPE
Ammerman and Cavalli-Sforza were the first to present a scientific, testable model of one of the most important processes in human history ͓37͔: the change from huntergatherer to agricultural economics ͑i.e., the Neolithic transi-tion͒.This transition triggered the acceleration of human population growth ͓38͔.The motivation to build a mathematical model of this process was the discovery that, according to archeological data, agriculture did not arise independently in different European regions.Instead, it spread gradually ͓39͔.Ammerman and Cavalli-Sforza proposed that this was not purely a process of cultural diffusion ͑imitation͒ but of physical ͑or demic͒ diffusion, i.e., of movement of farming communities.It has been pointed out that this hypothesis is backed by the experimentally measured genetic gradients in human populations ͓40,41͔, as well as by the common origin of Indo-European languages ͓42,43͔.Based on the hypothesis of physical diffusion, the wave-of-advance model was proposed ͓44͔ by making use of Fisher's equation.This is a very reasonable choice to find approximate results because of the simplicity of Fisher's approach.Recently it has been shown that agreement with the archeological data is improved by taking into account that this process took place in two dimensions ͓26͔.This certainly contradicts those criticisms of the wave-of-advance model based on the claim that it predicts a much higher velocity for the spread of agriculture as compared to that determined experimentally ͓45,46͔.The role of second-order terms has also been discussed recently ͓26͔.In this context, it is very important to explain how the values of the parameters, used in Ref. ͓26͔, were derived from anthropological observations.We think this point requires a very brief discussion here, so that the reader can judge the scientific character of the application considered.Such values will then be applied to the equations derived above, in order to present a more rigorous analysis than that in Ref. ͓26͔.

A. Determination of the values of a, D, and from field data
As explained in Sec.III, when the first farmers arrive in a geographical area, the population density is very small, p FIG. 1. Plot of the left-and right-hand sides of Eq. ͑26͒.It is seen that a solution to this equation exists provided that v is higher than a certain minimum value which will depend on the values of the diffusive and reactive parameters D, , and a. Ϸ0, and the growth function is FϷap.It means that, in the absence of migrations, the initial evolution of the population density will follow the approximate law dp/dtϭap ͓see, e.g., Eq. ͑6͔͒, which corresponds to an exponential growth.The value of a cannot be derived from archeological data because they do not reach such a level of detail.However, plausible values for a can be inferred from observations of populations that settled in empty areas.Birdsell was able to collect such data from the 18th century on the island of Pitcairn ͑East of Chile͒ and also from the 19th century on the islands in Bass Strait ͑between Australia and Tasmania͒.
What is impressive about these data is that a plot of the population size versus the generations of elapsed time gives almost exactly the same curve in both cases ͓47͔ .Combining this result with a mean generation time of 25 yr one finds a value of aϭ0.032Ϯ0.003yr Ϫ1 , with 80% confidence level.An estimation for the diffusion coefficient D ϭ͗⌬ 2 ͘/4 ͓see Eq. ͑18͔͒ can also be made from anthropo- logical data of the mobility of Ethiopian shifting agriculturalists and Australian aborigines.The corresponding values ͓44,48͔ yield a mean square displacement per generation of ͗⌬ 2 ͘/ϭ1544Ϯ368 km 2 /generation, with 80% confidence level.The parameter ͓which is twice the phenomenological delay time; see Eq. ͑19͔͒ is more difficult to measure.As in Refs.͓49͔ and ͓26͔, we assume that it can be approximated to the mean generation time.This corresponds to one migration per generation, although we would like to stress that future archeological observations could be very useful in determining to what extent this is a realistic estimate: the use of manure and crop rotation avoids having the land becoming exhausted ͓50͔, but in case early agriculturalists did not use these techniques, then the value of could certainly be smaller than the mean generation time.In the present paper, we will use a mean generation time of ϭ25 yr ͓44͔.This completes our brief discussion on the parameter values used in Ref. ͓26͔.We apply them below to the theory developed in the preceding sections.

B. Comparison to observations
The archeological data on the earliest recorded farming settlements in Europe yield a rate of advance for the expansion of the farming communities of 1.0Ϯ0.2km/yr ͓39,44͔ ͑see, e.g., Fig. 1 and the discussion in Ref. ͓26͔͒.By contrast, use of the mean values for a, D, and given above in Fisher's velocity ͑23͒ gives v →0 ϭ1.41 km/yr.Use of the same values in the hyperbolic equation ͑24͒, which was used in Ref. ͓26͔, yields v (2) ϭ1.00 km/yr.This value is completely within the experimental range, but we note that the difference with respect to Fisher's result is higher than 40%.Such a large difference makes it necessary to use the results in the present paper in order to determine whether or not the second-order approximation v (2) can be trusted.The point is thus to analyze the role of the terms involving derivatives of order higher than the second in Eq. ͑21͒.In order to do so, we use the same values of a, D, and as those given above in the exact solution ͑26͒.After solving this equation numerically in the manner explained in Sec.III and Fig. 1, we find vϭ0.98 km/yr.Thus we see that the second-order result v (2) is very similar to the exact one, and both of them yield a wave-front velocity that lies completely within the experimental range (1.0Ϯ0.2 km/yr͒.This shows that the hyperbolic velocity ͑24͒ is quite a reasonable approximation, in contrast to the classical velocity ͑23͒, which neglects the role of the time delay, for the application and parameter values considered here.For the sake of completeness, we mention that in fact one can also obtain the result v ϭ0.98 km/yr by solving Eq. ͑25͒ numerically for increasing values of N: the solutions are seen to converge rapidly ͑for example, the approximation Nϭ4, which corresponds to including derivatives of up to fourth order, already yields v (4) ϭ0.98 km/yr͒.However, we think that the procedure based on the application of the ϱ-order Eq. ͑26͒ should be more practical in general, since it will certainly require less computation time for possible applications in which many terms could be necessary.
Because the calculations above have been performed for the mean values of the reaction and diffusion parameters, a and D respectively, it is important to analyze the results for other possible values.These results are shown in Fig. 2. We see that there is good agreement between theory and experiment, and that this figure is very similar to Fig. 3 of Ref.
͓26͔.This result is important because Fig. 2 in the present paper has been obtained by including an infinite number of terms, whereas Fig. 3 of Ref. ͓26͔ corresponded to the second-order approximation.Comparing both figures, we also note that use of all higher-order terms yields a slightly FIG. 2. This figure shows the predictions taking into account infinite higher-order terms in the mathematical model ͓see Eq. ͑26͔͒.By contrast, Fig. 2 in Ref. ͓26͔ was obtained by including only terms of up to second order.The hatched rectangle gives the values for the reaction and diffusion parameters (a and D, respec-tively͒ implied by independent observations.The curves give the values of a and D for wavefront velocities of 0.8 km/yr, 1.0 km/yr, and 1.2 km/yr, according to the model derived in the text.Since the velocity inferred from archaeology is 1.0Ϯ0.2km/yr, there is good agreement between theory and experiment.lower wave-front velocity than that corresponding to the second-order equation ͑24͒.This velocity is, in turn, lower than Fisher's result ͑23͒, which neglects any possible effect of the time delay on the evolution of the system.
We have seen that the approach reported shows that the hyperbolic model presented in Ref. ͓26͔ is a reasonable approximation.This can be relevant from two points of view: ͑i͒ it shows that the hyperbolic wave-front velocity ͑24͒ is rather accurate in the application here considered; ͑ii͒ it opens the way to a more general approach to chemically reactive systems, a case in which it is very important to understand the propagation speed of wave fronts ͑see, e.g., Refs.͓21,22͔͒.Finally, it is worth mentioning that we have here derived and analyzed linear equations, i.e., F͑ p,p t ,p tt , . . .,p x , p y , p xy , p xx , p yy , p xy , . . .͒ϭ0, where F is a linear function in p, p t , etc. ͓51͔, but the method in Sec.III also applies to some nonlinear equations.For example, consider the time-delayed equation which, neglecting the delay time (⌼Ϸ0), has been proposed in the analysis of the profile evolution of a growing interface in solidification and crystallization processes ͓52͔.The second term on the right-hand side is nonlinear, and is a parameter related to nonlocal effects; without this nonlinear term we recover the hyperbolic equation ͑8͒ in one dimension.It is simple to follow exactly the same steps as in Sec.
III, but in this case the exponentials do not cancel out.Thus the linearization method does not yield a lower bound for the speed of the front.It would be interesting to apply variational techniques to this problem, which will be tackled in future work.We do expect, however, that the delay time will affect the speed of fronts even in the nonlinear case, which does not apply to the problem analyzed in the present paper but may be relevant in other applications.

V. CONCLUDING REMARKS
It had been previously shown that hyperbolic reactiondiffusion equations can be useful in the description of human population expansions ͓26͔.However, it was also shown that the corrections with respect to the classical, parabolic approach ͑based on Fisher's equation͒ can be relevant indeed ͓26͔.Thus the validity of a hyperbolic, second-order approach used was not clear.Here we have developed a model that allows one to include terms of up to an arbitrarily high order.We have shown how an infinite number of such terms can be taken into account ͓see Eq. ͑26͔͒.We have applied this to perform estimations of the velocity of spread of the Neolithic expansion in Europe.The results obtained are in good agreement with the empirical evidence from archaeology and anthropology that is available at present.This generalizes previous models of time-delayed reaction-diffusion ͓24-26͔, ͓36͔ and puts them on a more rigorous basis.
Before closing, we would like to mention that there are several additional fields of application of the results here reported.On one hand, chemically-reacting systems ͓21,22,53,54͔, superconductors ͓55͔, liquid crystals ͓56͔, and solidification ͓57͔ are topics in which extensive simulations are being performed in order to determine wave-front velocities: the methods presented here could be a useful analytical approach in situations such that the effect of a delay time could be important.On the other hand, here we have focused our attention on the Neolithic transition in Europe simply because the quantity of archeological observations is higher in this case than for other human expansions.However, there is evidence that similar processes took place in Africa ͓58͔, America ͓59͔, Polynesia ͓60͔, and China ͓61͔.The approach we have presented should, in our opinion, be applied to these expansions as soon as sufficient and reliable empirical data become available.Finally, we stress that the model presented is not restricted to a specific system and could thus be useful in the study of animal and plant expansions ͓62͔ and other biophysical topics in which reaction-diffusion is of utmost importance, such as the spread of epidemics ͓29,23͔, nerve conduction ͓29,63͔, cellular sensitivity ͓64͔, growth of bacterial colonies ͓65͔, and models of mitochondrial tissue ͓66͔.

extent these
When Eq. ͑5͒ is combined with the mass balance equation, one obtains the well known result ͑see, e.g., p. 236 in Ref. ͓29͔͒ ‫ץ‬p ‫ץ‬t ϭDٌ 2 pϩF, ͑6͒ which is Fisher's equation.Here FϭF(p) is the source function corresponding to reactive processes in the system.approximations are reliable.It seems reasonable to try to conserve the simplicity of the model in Ref. ͓26͔ as far as it is possible to do so.Thus we will simply keep an arbitrary number N of terms in the expansions above, i.e., we rewrite Eq. ͑9͒ as