Stochastic model for population migration and the growth of human settlements during the Neolithic transition

We present a stochastic two-population model that describes the migration and growth of semisedentary foragers and sedentary farmers along a river valley during the Neolithic transition. The main idea is that random migration and transition from a sedentary to a foraging way of life, and backwards, is strongly coupled with the local crop production and associated degradation of land. We derive a nonlinear integral equation for the population density coupled with the equations for the density of soil nutrients and crop production. Our model provides a description of the formation of human settlements along the river valley. The numerical results show that the individual farmers have a tendency for aggregation and clustering. We show that the large-scale pattern is a transient phenomenon which eventually disappears due to land degradation.


I. INTRODUCTION
The wave of colonization by migrating farmers and the establishment of farming communities in Europe between 5000 and 3500 B.C. is currently a topic of great interest in prehistoric archaeology, linguistics, and anthropology ͓1-3͔.Ammerman and Cavalli-Sforza developed a model for the expansion of farming as a demic diffusion which spread into Europe in the form of a wave of advance ͓4͔.Using radiocarbon dates, they found that farmers spread at an average rate of about 1 km a year ͓5,6͔.Interest in simulation and spatial modeling of the spread of agriculture has been growing rapidly in the last decade, especially in the physics community ͓7-9͔.One of the main reasons for this is that the geographical spread of population can be effectively described by the classical Fisher-Kolmogorov-Petrovskiy-Piskunov ͑KPP͒ equation and its various generalizations ͓10,11͔.These models have recently attracted considerable interest in physics and biology, because of the huge number of potential applications.Cohen ͓12͔ developed a model that takes account of a varying environment.He derived a modified Fisher-KPP equation in which the growth rate depends on food production, land fertility, etc. ͑see also ͓13͔͒.Fort and Mendez applied a time-delayed theory for the Neolithic transition ͓7͔ which involves a hyperbolic correction to the Fisher-KPP equation.The transition from hunter-gathering to farming did not happen in a uniform way, and that is why Davison et al. ͓14͔ took into account both advection and spatial variation in diffusivity and carrying capacity, in the framework of the Fisher-KPP equation.Aoki et al. ͓15͔ studied the spread of farmers into an area occupied by huntergatherers in terms of the system of reaction-diffusion equations.Ackland et al. ͓13͔ analyzed a wave-of-advance model involving the interaction of three populations: Neolithic farmers, hunter-gatherers, and converts.
The main objective of those works was to reproduce the observed rate at which agricultural expansion takes place in Europe.Despite the interest in the establishment of farming communities in Europe, there remains no models explaining the spatial structure formed behind the wave of advance.The main challenge of our work presented below is to set up a stochastic model which provides an explanation to not only the propagation of farming, but also the formation of human settlements.We develop a model that demonstrates the tendency of the distribution of population to form clusters.Furthermore, we show that this large-scale pattern in population evolution is a transient phenomenon, which disappears due to degradation of land in the form of an extinction wave.Our specific motivation is the successive migration of settlements along parallel river valleys in the Tripolye-Cucuteni system, which has been thoroughly investigated and documentede.g., ͓16͔.This system can be considered as being one dimensional and so is particularly amenable to numerical investigation.

A. Two-population model for semisedentary foragers and sedentary farmers
In our model the population consists of semisedentary foragers and sedentary farmers who share the same territories.The semisedentary foragers are the population of individuals randomly moving from place to place along a river valley and searching for food and other resources.An implicit consequence of this behavior is the foundation of new settlements ͑large localized values of population density͒ and an interchange between farming and foraging populations.On the contrary, the sedentary farmers are individuals who do not migrate.They live in small villages scattered near cultivated land in the valley of a major river.Their main activities are the cultivation of soil and crop production.In this paper we are interested in the total population density n͑x , t͒ at location x along the river at time t.We define this density as n = n 1 + n 2 , where n 1 is the density of semisedentary foragers and n 2 is the density of sedentary farmers.
We assume that there is not a strict distinction between foraging and sedentary lifestyles.There are always random transitions from a sedentary to a foraging way of life and vice versa; these transitions depend strongly on the local food supply.This is one of the main features of our random walk model.The key feature of the movement of semiseden-tary foragers is that they do not jump from place to place completely randomly.Unlike Brownian particles in physics, in general the migration of people cannot be explained by a standard diffusion law, in which the flux is proportional to the gradient of number density of individuals.There are many extensions of the standard diffusion techniques-e.g., ͓17͔.For example, Cohen introduced a forced flux corresponding to the population tendency to migrate toward more fertile land ͓12͔.This idea is very similar to chemotaxis, where bacteria and other organisms migrate up a concentration gradient of chemoattractant.To describe the random migration of semisedentary foragers and random transitions from one lifestyle to another, we adopt a biased random walk whose statistical characteristics depend on the local food supply.In our model, the probability of a random migration event resulting in a jump z in the time interval from t to t + ⌬t is ⌬t ͓17͔.The probability of a transition from a foraging lifestyle to farming is ␣ 1 ⌬t.The probability of the conversion of farmers to semisedentary foragers is ␣ 2 ⌬t.Thus we introduce a new variable, the local crop production per individual per year q͑x , t͒, such that the frequency of jumps and transition rates ␣ 1 and ␣ 2 depend on the crop production: It is natural to assume that the frequency ␣ 2 ͑q͒ is a decreasing function of q.That is, when sedentary farmers are not able to produce enough food to sustain their population, some of them start to migrate from their neighborhood at the rate (q͑x , t͒).

B. Balance equations for population densities
We now set up the balance equation for the density of semisedentary foragers, n 1 ͑x , t͒.According to our random walk model, the density of foragers n 1 at location x at time t + ⌬t can be written as follows: This equation is the conservation law for foragers.The first term on the right-hand side represents those foragers who stay at location x and do not move during time ⌬t and do not become sedentary farmers.The second term gives the number of foragers who arrive at x during time ⌬t from different places x − z, where the jump distance z is distributed by dispersal kernel ͑z͒.The last term gives the number of sedentary farmers who become semisedentary foragers during time ⌬t.For simplicity in this preliminary investigation we have assumed that the reproduction term for foragers is negligible.Inclusion of such a term in ͑1͒ would be technically straightforward.
The sedentary farmers do not migrate, and their density n 2 ͑x , t͒ obeys the balance equation involving logistic growth and lifestyle transitions: where the last term gives the number of foragers who become the farmers during the time ⌬t.Here r is the growth rate of the sedentary population.Typical values of the parameter r are 0.01-0.03yr −1 ͓14͔.The carrying capacity K is, in principle, an increasing function of the local crop production q.In numerical simulations, K is assumed to be a constant.Note that Cohen ͓12,13͔ derived a modified Fisher-KPP equation in which the growth rate r depends on food production and land fertility.
In the limit ⌬t → 0, from ͑1͒ and ͑2͒ we obtain two differential equations for n 1 ͑x , t͒ and n 2 ͑x , t͒: In this paper we are interested in the evolution of population density n͑x , t͒ on characteristic time scales around of 100 yr.Therefore we can adopt a local equilibrium for the two populations, describing the balance between the nomadic and sedentary ways of life in proportions p and 1 − p.We write where the proportion of foragers p is a function of the crop production q-i.e., p = p͑q͒.͑6͒ We can derive an explicit expression for the function p͑q͒ in terms of the transition rates ␣ 1 ͑q͒ and ␣ 2 ͑q͒.Equations ͑3͒ and ͑4͒ involve four characteristic rate parameters: the growth rate r, the jump rate , and the two transition rates ␣ 1 and ␣ 2 .If we assume that both transition rates are large compared to the growth and jump rates, then the "fast" transition process can be averaged.The fast local dynamics of population densities n 1 and n 2 is governed by the reduced equations For the initial densities n 1 ͑0͒ and n 2 ͑0͒, the solution is When ͑␣ 1 + ␣ 2 ͒t is large, the second terms in ͑8͒ and ͑9͒ tend quickly to zero.In other words, for a large intermediate time If we take into account "slow" migration and population growth, the density of population splits locally into the numbers ␣ 1 +␣ 2 n of sedentary farmers and Therefore the dependence of p on the transition rates ␣ 1 ͑q͒ and ␣ 2 ͑q͒ is One of the advantages of this formula is that instead of estimating two functions ␣ 1 ͑q͒ and ␣ 2 ͑q͒, we need to estimate only the single p͑q͒.
The evolution equation for biased migration, Eq. ͑3͒, and population growth, Eq. ͑4͒, can be rewritten as one balance equation for the total population density.By adding Eqs.͑3͒ and ͑4͒ using ͑5͒, we obtain a single equation for n͑x , t͒:

͑13͒
Note that in ͑13͒ the number of semisedentary foragers moving depends on the local crop production and thus the number density of population, n.Thus we are describing an essentially nonlinear movement of population.Now we need to find an approximate form for the function p͑q͒.It is natural to assume that for a large crop production q, the population mostly consists of sedentary farmers-that is, p close to zero-while for a low crop production the population mostly is semisedentary foragers-that is, p is close to 1.A simple approximation for p might be where H͑x͒ is the Heaviside step function: H͑x͒ =0 if x Ͻ 0; H͑x͒ =1 if x ജ 0. Thus farmers do not migrate if the crop supply is greater that the critical value q min .A better approximation which we use in our numerical simulations is a linearly decreasing function in the interval between the minimum value q min and the maximum value q max : p"q͑x,t͒… = Ά p max , q ഛ q min , − aq + b, q min Ͻ q Ͻ q max , p min , q ജ q max , · where a = −p max +p min q max −q min and b = p max +p min 2 + a 2 ͑q max + q min ͒ ͓see Fig. 1͑a͔͒.Based on ͓18͔, we take q min , q max = 300, 736 kg per person per year, respectively.For the results discussed below we set p max = 0.95 and p min = 0.05.It should be noted that the exact values for q min and q max are not crucial to our results.
Regarding the dispersal kernel ͑z͒, we assume that it has an exponential form for large jumps: where l is the characteristic length scale.Since farmers do not travel very small distances, the function ͑z͒ should also satisfy the condition ͑0͒ = 0.That is why in numerical simulations we introduce the concept of the most probable value for jumps z m that corresponds to the point where ͑z͒ has a maximum value.The function ͑z͒ is shown in Fig. 1͑b͒, where z m is just below 10 km.Clearly the choice of a functional form for is rather arbitrary.We have just attempted to choose pragmatically a form with the desired characteristics such as ͐͑z͒dz = 1 and ͑z͒ → 0 as z → ϱ.Note that there is some evidence for long-ranged kernels such as given by a power law with infinite second moments ͓19͔.The most probable value z m can be estimated from archaeological evidence from the Tripolye-Cucuteni culture that settlement separations were typically 10-20 km ͓16,18͔.We believe that typical jumps would be of comparable magnitude.For computational reasons we introduce a maximum cutoff for ͑z͒.Since the typical length of the river is around 500 km, then it seems reasonable to assume that the maximum cutoff is around 50 km.Note that we allowed "forward" jumping only, but verified that our results were essentially unaltered when a small probability of "backward" jumping was also allowed.

C. Equation for crop production
We assume that the human population derives most of their food from the cultivation of land and thus suggest the following formula for the local food production q͑x , t͒: where F͑x , t͒ is the density of soil nutrients, n 0 is the constant, ␣ is the production rate coefficient, and ␤ is the parameter that determines how the yield depends on the nutrients.This equation describes how the rate of food supply q FIG. 1.The functions ͑a͒ p͑q͒ and ͑b͒ ͑z͒.
increases due to the increase in the population density n and how the degradation of land ͑the decrease of soil nutrients F͒ leads to a decrease of food production through the factor 1 − e −␤F͑x,t͒ .Note that the factor n͑x,t͒ n 0 +n͑x,t͒ describes a tendency toward group solidarity that increases the efficiency of food production.Numerical experimentation for a typical solution showed that in the absence of this term, the clustering phenomenon reported in this paper did not occur.The main difference between our approach and Cohen's model ͓12͔ is that we use the effect of group solidarity while Cohen used the food production function which is inversely proportional to the population density.There is no clustering phenomenon for this function.We know from existing data that the yield-i.e., the production of food per unit area-can be up to 0.0736 kg/ m 2 yr ͓18,20͔.We estimate that the area cultivated by people is approximately 10 4 m 2 per person.Thus, we have an estimate for ␣ of 736 kg per person per yr.In general, it is difficult to model how the crop production depends on the nutrients, since this relation is strongly dependent on the environmental conditions.One of the existing models is given by the Mitscherlich-Baule yield response ͓21͔, which takes into account productivity losses due to soil erosion.Following this model, in ͓22͔ the factor 1 − e −␤F͑x,t͒ has been proposed.The value ␤ Ӎ 890 m 2 kg −1 is found for the corn cultivation.However, other studies ͓23͔ seem to suggest lower values ͑␤ Ӎ 65 m 2 kg −1 ͒.

D. Land degradation equation
We adopt the following model for nutrient depletion and the corresponding land degradation.We assume that equation for the density of soil nutrients F͑x , t͒ has the form ‫ץ‬F͑x,t͒ ‫ץ‬t = 1 − 2 F͑x,t͒ − ␥q͑x,t͒n͑x,t͒, where 1 is the rate at which soil nutrients regenerate naturally, 2 is the rate of nutrient depletion due to environmental reasons ͑erosion, flooding, etc.͒, and ␥ is the rate at which nutrients are depleted due to the harvests.
Here we assume that natural nutrient depletion is much slower than depletion due to human activity.In order to estimate ␥, we can use the data from ͓20͔.These data correspond to prehistoric agriculture in Hawaii; however, similar values are cited for Sub-Saharan agroecosystems and there are no data available ͑as far as we know͒ for other regions.So using that study we estimate the rate of nutrient depletion, ␥ = 5.43ϫ 10 −3 kg P / kg, where kg P denotes kilograms of phosphorus in the soil.However, other studies mention that the total concentration of nutrients in crops can be up to 3% ͓23͔.It would suggest that the value ␥ = 0.03 kg P / kg can be taken as a maximum threshold-cf.͓14͔.Of course, other nutrients are also important, but phosphorous can be regarded as a proxy for them.Soil regeneration is a very complex process and depends on many environmental parameters.Nevertheless, in some agroeconomics models regeneration is modeled as a constant ͓24͔ as we have assumed in our model.

III. NUMERICAL RESULTS
Numerical simulations of the nonlinear integral equation ͑13͒ together with Eqs.͑15͒ and ͑16͒ reveal a very interest- FIG. 2. n͑x , t͒ for the model with jump rate = 0.1 and depletion rate ␥ = 0.01 discussed in the text.Panels ͑a͒-͑h͒ show successively the solution at times 0 -700 yr at intervals of 100 yr.
ing dynamical behavior: the emergence of a large-scale pattern in the population density.This phenomenon can be interpreted as the formation of human settlements along a river valley.Thus our model provides an explanation for the formation of settlements as a dynamical phenomenon.The individual farmers have a tendency for aggregation and clustering as a result of nonlinear random migration.Moreover, our model describes not just a population clustering, but subsequent decay of these clusters ͑settlements͒ due to land degradation.
Our model has a rather large number of parameters, but fortunately most of them can be estimated from existing sources, as indicated above.For simplicity, the value of is assumed to be a constant and determined by the requirement that the velocity z m be of order of the observed migration speed of about 1 km yr −1 .For numerical work we use as unit of time 1 yr and measure densities in units of m −2 .We take the values r = 0.03, ␣ = 736, ␤ = 200, K =10 −4 , 1 =10 −4 , 2 =0, q min = 300, q max = 736, and l = 10.Note that in this initial study we take K to be a constant, rather than a function of q.The value of K can be roughly estimated from Roman and medieval evidence for densities of a little more than 1 hectare per person ͓14͔.Our choice of r is also guided by ͓14͔. 2 is set to zero for simplicity and to minimize the number of parameters in the absence of strong experimental estimates.We explore ranges 0.01ഛ ␥ ഛ 0.03 yr −1 and 0.05 ഛഛ0.3 yr −1 .Our initial conditions were rather arbitrary: n = 0.1 K in 0 ഛ x ഛ 15 km, n = 0 elsewhere, although from experiments with other initial conditions with n͑x ,0͒ nonzero in x Շ 100 it appears that results are quite insensitive to details of the initial state.Again, for simplicity, we set the initial value F͑x ,0͒ to be uniform and equal to F 0 .Our "stan- FIG. 3. n͑x , t͒ for the model with jump rate = 0.2 depletion rate ␥ = 0.01 discussed in the text.Panels ͑a͒-͑l͒ show successively the solution at times 0 -1100 yr at intervals of 100 yr.
dard" choice was, again rather arbitrarily, F 0 = 0.01, and we also explored a larger value.If the range of x were to be extended to negative values, the evolution of n͑x , t͒ would be approximately symmetric about x =0.
We show in Fig. 1 the functions p͑q͒ and ͑x͒.In Figs. 2 and 3 we present results n͑x , t͒ for typical cases = 0.1 and 0.2, respectively, with ␥ = 0.01, at intervals of 100 yr.We see clear evidence of clustering of population ͑settlements͒ and of their subsequent decay.Panel ͑a͒ shows the "arrival" of farmers at location 0 Ͻ x Ͻ 15 along the river.After 300 yr, one can see clear evidence of clustering of population ͑settle-ments͒.Panel ͑f͒ in Fig. 3 shows their subsequent decay and reappearance of clusters behind the extinction wave.The entire population decays over about 1000 yr, depending to some extent on the exact choice of parameters.
A general feature of our modeling is that most population clusters grow and decay in situ, without significant movement.This is illustrated in Fig. 4, where we show at increased time resolution the evolution of population density for the model of Fig. 3 over a limited time interval, and in Fig. 5 we show the evolution of population density n for two selected population clusters ͓i.e., n͑x , t͒ at fixed values of x͔ in the chosen model.Although the model is one dimensional, we can make an order-of-magnitude estimate of the total population of a cluster by taking the linear extent of a cluster as the diameter of a circular settlement.This gives typical figures of order 10 4 individuals.
We also explored the dependence of the speed of advance of the population and the separation of the clusters as a function of the parameters.The speed of advance is rather ill defined, but for simplicity we monitored the movement of the main peak of the distribution.͑Arguably we could have, e.g., studied the position of the cluster farthest from the origin.͒We found the speed of advance to depend approximately linearly on both and ␥.The separation of clusters is independent of both and ␥.If Շ0.05, the clustering phenomenon does not occur.While we have concentrated our efforts on studying variations among a manageable subset of parameters while retaining fixed plausible values for the others, we also made a small study of the effects of varying F 0 and 1 .Increasing the value of F 0 , up to 0.1 from our canonical value of 0.01, prolongs the time scale of the population evolution by a factor of up to about 2. Values of F 0 significantly smaller than 0.01 do not give clustering for typical choices of the other parameters.For larger values of 1 , տ10 −4 , the regeneration is so strong that, although clustering occurs initially much as for 1 =10 −4 , the final state is a spatially uniform population.For smaller values of 1 -e.g., 1 =10 −5 -we still see strong clustering, but the overall phenomenon persists for a significantly shorter time, as the episodes of cluster regeneration seen when 1 =10 −4 do not now occur.Physically, we can deduce that for soils with low fertility clustering does not occur.Soils are exhausted more rapidly for lower regeneration rates, and for sufficiently rapid regeneration soil exhaustion does not occur.While these findings are intuitively unremarkable, it is essential that our model should display these properties.

IV. DISCUSSION AND CONCLUSIONS
Our stochastic model gives an explanation for the formation of human settlements as a dynamical phenomenon.The individuals have a tendency for aggregation and clustering as a result of nonlinear random migration.The key feature is that random migration and the transition from a sedentary to a foraging way of life and backward is strongly coupled with local crop production.Moreover, our model describes the subsequent decay of the clusters due to land degradation and corresponding food crisis in the form of an extinction wave.We took our inspiration and guidance from the successive south to north migration of the Tripolye-Cucuteni cultures along parallel river valleys ͓16͔.For reasonable choices of parameters and without fine-tuning, our solutions exhibit a distinct clustering of population, at intervals of 10-30 km.Plausibly these clusters are estimated to have around 10 4 individuals.These estimates are all consistent with what is known about the Tripolye-Cucuteni cultures ͓16,18͔.Overall, we feel that our quite naive modeling captures some of the essence of the clustering seen in population migration and points the way for more sophisticated modeling.We can think of a number of significant developments in the future, including taking K = K͑q͒ and allowing for a latency effect, in which the population in a cluster ͑settlement͒ has a reduced probability of moving, representing an attachment to the investment of building a house and clearing land.It is easy to allow for additional resources available from, e.g., forest or river, which are maybe harder to exhaust than the fertility of the land and so can act as a reservoir of resource.It would be interesting to apply our model to understand societal collapses to which environmental problems such as habitat destruction, soil degradation, and overpopulation contribute ͓25͔.
FIG. 4. Details of the evolution of the model illustrated in Fig.3between times 200 and 300 yr.Panels ͑a͒-͑e͒ show n͑x , t͒ at intervals of 25 yr between t = 200 yr and t = 300 yr.