Time-Delayed Theory of the Neolithic Transition in Europe

The classical wave-of-advance model of the neolithic transition (i.e., the shift from hunter-gatherer to agricultural economies) is based on Fisher’s reaction-diffusion equation. Here we present an extension of Einstein’s approach to Fickian diffusion, incorporating reaction terms. On this basis we show that second-order terms in the reaction-diffusion equation, which have been neglected up to now, are not in fact negligible but can lead to important corrections. The resulting time-delayed model agrees quite well with observations. [S0031-9007(98)08286-6]

In reactive systems, Fickian diffusion leads to parabolic reaction-diffusion (PRD) equations.Parabolic reactiondiffusion equations have been applied to the spread of advantageous genes [8], dispersions of biological populations [9], epidemic models [10], etc.But, if information propagates at a finite speed, linear flux-force laws-and thus PRD equations-do not hold [1].This unphysical feature can be avoided by making use of hyperbolic reaction-diffusion (HRD) equations [11,12].Thus HRD equations, instead of the usual PRD equations, should be regarded as the first choice from a conceptual perspective.Hyperbolic reaction-diffusion equations have been very recently applied to the spread of epidemics [13], forest fire models [14], and chemical systems [15].
An interesting application of PRD equations arose after archaeological data led to the conclusion that European farming originated in the Near East, from where it spread across Europe.The rate of this spread was measured [16], and a mathematical model was proposed according to which early farming expanded in the form of a PRD wave of advance [17].Such a model provides a consistent explanation for the origin of Indo-European languages [18], and also finds remarkable support from the observed gene frequencies [19].However, this PRD model predicts a velocity for the spread of agriculture that is higher than that inferred from archaeological evidence, provided that one accepts those values for the parameters in the model that have been measured in independent observations [17].Here we will analyze this problem by means of a HRD model.
Let p͑x, y, t͒ stand for the population density (measured in number of families per square kilometer), where x and y are Cartesian coordinates and t is the time.We assume that a well-defined time scale t between two successive migrations exists.We begin, as usual [20], noting that, between the values of time t and t 1 t, both migrations and population growth will cause a change in the number of families in an area differential ds dx dy, i.e., ͓p͑x, y, t 1 t͒ 2 p͑x, y, t͔͒ds ͓p͑x, y, t 1 t͒ 2 p͑x, y, t͔͒ m ds 1 ͓p͑x, y, t 1 t͒ 2 p͑x, y, t͔͒ g ds , where the subindices m and g stand for migrations and population growth, respectively.We denote the coordinate variations of a given family during t by Dx and Dy.The effect of migrations on the evolution of p͑x, y, t͒ will be derived here by means of a simple extension of Einstein's model of Fickian diffusion [21].The migration term in Eq. (1) where F͑x, y, t͒ is the change, due to births and deaths, in the local population density per unit time.
From Eqs. ( 1), (2), and ( 5) we find p͑x, y, t 1 t͒ 2 p͑x, y, t͒ `p͑x 1 Dx, y 1 Dy, t͒f͑Dx, Dy͒dDx dDy 2 p͑x, y, t͒ We now assume that t, Dx, and Dy are small enough in comparison with the measured time intervals and distances (t ø t, Dx ø x, and Dy ø y) so that both sides of Eq. ( 6) may be approximated by their second-order Taylor expansions, Making use of Eqs. ( 3) and ( 4), this reduces to ≠p ≠t 1 t 2 where we have defined with D p Dx 2 1 Dy 2 .In the special case of a vanishingly small delay time, t !0, Eq. ( 7) reduces to the usual PRD equation [8].In the limit F ! 0 and t !0, Fickian diffusion is recovered [21].Our derivation of the HRD equation (7), in contrast to previous ones [11,12], makes it possible to note that in the application we are dealing with the interpretation of t is completely different from that in a reacting-diffusing gas mixture (where collisions have a negligible duration in comparison with the travel time between two subsequent collisions).The time interval between two successive migrations is the sum of the time of travel (which may of the order of some days or weeks) and the time of "residence," i.e., the time interval between the arrival of a family and the subsequent migration (presumably of about one generation [22]).Therefore, in our case, t is approximately the time of residence.
We will assume that the population growth can be described by where a is the initial growth rate and p max is the saturation density.Equation ( 9) is the logistic growth function, and it compares favorably to a wealth of experimental results [23].It has been recently shown that Eq. ( 7) leads to wave fronts with asymptotic speed of propagation [12] y This result would not change if we considered an exponential growth function.We note from Eq. ( 10) that the delay time t slows down the wave front, as was to be expected intuitively.In the limit t !0, Eq. ( 10) becomes which is the basis of the classical wave-of-advance theory of the neolithic transition in Europe [17].
In a different context, wave fronts for the autocatalysis A 1 B !A 1 A have been studied for a reaction term of the form F kp͑p 0 2 p͒, where p is the density of A, p 0 is the total density, and k is the rate constant.In this case, Eq. ( 10) can be written as y r!0 y 1 1 t 2 kp 0 , i.e., an additional term proportional to the rate constant is predicted.Such a term has been proposed before [15], although no explicit form for it has been derived.
We return to the problem under consideration.Equation (8) may be written as in fact, a well-known result for two-dimensional diffusion.Here t is the mean square displacement during the time interval t that separates two successive migrations.Previous approaches did not take the factor 1 4 in Eq. ( 12) into account but relied on the approximation [17] The front velocity for the expansion of agriculture can be predicted from Eq. ( 10) provided that independent estimations for the values of t, a, and D are available.As in Ref. [17], we consider twenty-five years for the average generation, t 25 yr.Let us, for the moment, assume an initial growth rate of 3%, i.e., a 0.03 yr 21 [19] and t 1700 km 2 ͞generation [17].
In Fig. 1 we reproduce the archaeological data for the spread of the neolithic transition in Europe [17].Use of the values above into the approximation (13) (which was applied in Ref. [17]) and Eq. ( 11) yields a front velocity of y t!0 2.86 km͞yr.In contrast, Eqs. ( 11) and ( 12) yield a velocity of 1.43 km͞yr.The dasheddotted and dashed lines in Fig. 1 are the best fits with these slopes, respectively.The front speed implied by Eq. ( 13) is much higher than that inferred from the data.The prediction from Eq. ( 12), on the other hand, shows the usefulness of the wave-of-advance theory (provided that the factor 1 4 in Eq. ( 12) is included, as has been done here for the first time).Still, a lower front speed is clearly FIG. 1.Comparison between empirical evidence and theoretical models for the spread of agriculture in Europe.The points are the data already analyzed in Ref. [17], distances being measured as great circle routes from Jericho (the presumed center of diffusion).Dates are conventional radiocarbon ages in years before present (BP).The solid line is the regression by Ammerman and Cavalli-Sforza (correlation coefficient r 0.89) [17] implied.This may also be seen from the regression (solid line) in Fig. 1 [16,17]: Its slope yields a velocity of 1.00 km͞yr.When the same values for the parameters as above are introduced into Eqs.( 10) and ( 12), one obtains a prediction of y 1.04 km͞yr.Therefore, the prediction from Eq. ( 7) agrees very well with observations (see the values of x 2 in Fig. 1).
Changing the value of the average generation leads to much the same result (e.g., we might consider t 28 yr [17] in Eqs. ( 10) and ( 12), which yields y 0.95 km͞yr).In Figs. 2 and 3, the curves labeled with number 1 give the possible values of a and D 2 t compatible (according to each one of the three models) with the observed speed of 1 km͞yr.In fact, the principal axis was preferred [16,17] to perform a fit to the data (solid line in Fig. 1), but the two regressions (distances versus dates and vice versa) imply a velocity range between about 0.8 and 1.2 km͞yr [16].In Figs. 2 and 3 we also include curves corresponding to these velocities.The hatched regions in these figures correspond to likely ranges of the parameters, and have been obtained as follows.Birdsell [24] was able to collect detailed evolution data of two human populations which settled in empty space.A fit of these data (either to an exponential or to logistic curve) yields a 0.032 6 0.003 yr 21 , with 80% confidence level (C.L.).Values for t have been derived [17,25] from observations of Ethiopian shifting agriculturalists and Australian aborigines.A statistical analysis of these values yields   t , according to the available empirical evidence.Each curve is labeled with a number, which gives the considered front velocity (in km͞yr).The three lower curves correspond to the approximation (13), which was used in Ref. [17].The three upper curves correspond to the same model, but taking into account Eq. ( 12) for twodimensional diffusion.it is seen that the use of this correction (three upper curves), instead of the classical approach (lower curves), makes the nondelayed model marginally consistent with the experimental range of a and t values (hatched rectangle).This shows that the first-order, or PRD, approach is compatible with the demographic data, but only assuming extreme values for the parameters.In the second-order, or HRD, model, Eq. ( 10) holds, from which Fig. 3 follows.From this figure, we conclude that the agreement between the available empirical data and our time-delayed model is very satisfactory, in spite of the simplicity of the latter.This shows that it is not necessary to assume extreme values of the measured parameters, provided that one is willing to accept the hypothesis that a time interval of about one generation separates successive migrations.
Finally, it is worthwhile to take into account that simulations of the neolithic transition in Europe for the values a 0.02 yr 21 and D 2 t 973.4 km 2 ͞generation yield a front velocity of y 1.09 km͞yr [25].Making use of Eq. ( 11), the two-dimensional result (12) and the above values for a and D 2 t , we obtain y t!0 0.88 km͞yr.In fact, in Ref. [25] an irregular lattice was considered.Since Eq. ( 12) can be easily generalized into D 1 2n t , where n is the number of dimensions and 2n is the number of neighbors in a simulation lattice, we note that a decrease in the number of neighbors corresponds to an enhancement of diffusion, as was to be expected intuitively.If we take into account that, for the irregular lattice used in Ref. [25], the mean number of neighbors was 3.4, making use of Eq. ( 11) we estimate y t!0 0.96 km͞yr, which is closer to the value 1.09 km͞yr observed in simulations.Since the simulations in Ref. [25] included corrections due to the acculturation of hunter-gatherers by farmers (which could lead to a more rapid agricultural expansion), we conclude that there is reasonable agreement between our results and the simulations.Of course, the more general result (10) cannot be compared to those in Ref. [25], simply because that study is based on Eq. ( 7) in the limit t !0.
We would like to stress that, unless one is willing to accept an instantaneous propagation of signals, HRD equations are much more natural than the usual PRD equations.
Physical models could also contribute to the modeling of genetic gradients, to the possible contributions of higherorder terms, etc.

4 D 2 t 4 D 2 t 2 t
FIG.1.Comparison between empirical evidence and theoretical models for the spread of agriculture in Europe.The points are the data already analyzed in Ref.[17], distances being measured as great circle routes from Jericho (the presumed center of diffusion).Dates are conventional radiocarbon ages in years before present (BP).The solid line is the regression by Ammerman and Cavalli-Sforza (correlation coefficient r 0.89)[17].The other three lines are least-squares fits with slopes calculated from the classical wave-of-advance model with D ഠ D 2 t

D 2 t
FIG. 2. Predictions of the classical wave-of-advance model [Eq.(11)].The hatched regions correspond to the range of allowed values for a and

FIG. 3 . 2 t
FIG. 3. The same as Fig. 2, but for the time-delayed waveof-advance model we have presented [Eq.(10)] and an average generation of t 25 yr.It is seen that, for almost all of the likely values of the parameters a and D 2 t (hatched rectangle), the predictions of the model are consistent with the observed front velocity ͑1.0 6 0.2 km͞yr͒.