Anisotropic dispersion, space competition and the slowdown of the Neolithic transition

The front speed of the Neolithic (farmer) spread in Europe decreased as it reached Northern latitudes, where the Mesolithic (hunter-gatherer) population density was higher. Here, we describe a reaction–diffusion model with (i) an anisotropic dispersion kernel depending on the Mesolithic population density gradient and (ii) a modified population growth equation. Both effects are related to the space available for the Neolithic population. The model is able to explain the slowdown of the Neolithic front as observed from archaeological data.


Introduction
The spread of the Neolithic, one of the most important socioeconomic changes in human history, has been widely studied using physical models in recent years (for a review, see [1]). The Neolithic expansion has been tackled from different approaches such as age-structured population models [2], population spread along rivers [3] and settlement formation [4].
Here, we will focus on the fact that the spread of the Neolithic in Europe was not homogeneous from the macroscopic point of view. Archaeological observations show that, as the front propagated from the Near East across Europe, its speed slowed down as higher latitudes were reached [5].
This decrease of the front speed can be intuitively seen from figure 1, which shows the arrival time of the Neolithic across Europe. The arrow on the map represents the average direction along which the expansion from the Near East to the Baltic Sea took place (within the rectangle). In figure 1, it can be seen that the distance advanced during 500 years is lower on reaching northern latitudes (a quantitative analysis will be presented in section 5).
Although it would seem that the more intuitive reason for the decrease in speed is the time needed by crops to adapt to temperate climates, evidence exists that this effect was, in fact, minimal [6]. Indeed, when establishing their settlements in colder regions, Neolithic populations just cultivated the more adaptable crops and dropped the rest.
From archaeological studies, one of the most accepted reasons for the presence of a gradient in the front speed when spreading to the North of Europe is the presence of Mesolithic hunter-gatherer populations [7], which had higher densities at Northern latitudes. Thus, motivated by the observational data, in this paper we extend a homogeneous model [8] to study how the presence of indigenous Mesolithic populations affects the speed of the Neolithic invasion front.
We describe a reaction-diffusion model for Neolithic population density with a directiondependent dispersion kernel determined by the space dependence of the Mesolithic population density. We also introduce in this model the effect of the presence of Mesolithic populations on the Neolithic population growth process. We compare the results from the model with archaeological data [9].

Anisotropic dispersion kernel
In the case we assumed that the spread of the Neolithic front took place in a homogeneous space, it would be reasonable to consider that the probability φ to jump would be the same in all directions; thus, mathematically we would have [8,10] φ(x, y; θ, ) that is, the jump probability could be expressed as a function ψ that depends only on the jump distance, , and is independent of the jump direction θ or the position in space (x, y). We have assumed that ∞ 0 ψ( )d = 1. However, Neolithic individuals do not move in a homogeneous space, since the density of Mesolithic individuals they encounter depends on the position and direction they move. Then, for a given position (x, y), the Neolithic individuals will preferably move in the direction along which they encounter a lower Mesolithic population density, i.e. along the direction where more free space is available.   [9]. The arrow corresponds to the y-direction in our model. Thus, we can assume that, in this situation, the jump distance probability distribution (1) will be modulated by the available space, s, at the final jump point (x + x , y + y ), in each direction θ = tan −1 ( y / x ) and for every jump distance = 2 x + 2 x . Thus, the dispersion kernel is now of the form where α is a normalization constant. We now need a mathematical expression for the available space s(x + x , y + y ). If M max is the carrying capacity for Mesolithic populations and M(x, y) the actual density of Mesolithic individuals at the position (x, y), then the fraction of occupied space at this point can be expressed as Thus, the fraction of space available for Neolithic settlements is and the space-dependent jump distance probability (2) can be written as 4 For simplicity, we assume that the variation in Mesolithic population density takes place mainly in one direction, y, in figure 1, whereas it remains approximately constant along the x-direction, i.e.
Now, Taylor-expanding the term within square brackets in equation (6), we obtain that the spacedependent jump distance probability is approximately Normalizing equation (7), we obtain that the normalization constant α is and the jump distance probability becomes We can see from equation (9) that if the Mesolithic (indigenous) population density M increases along the direction y, then the probability of Neolithic invaders to jump forward (θ = π/2) is minimum and the probability to jump backwards (θ = 3π/2) is maximum.

Population growth
In population dynamics, a commonly used expression to describe the first-order variation in population density due to population growth (reproduction minus deaths) is the logistic growth equation [8,11,12], where a is the initial growth rate, N max the carrying capacity and N the density of the Neolithic population. The subindex g stands for population growth, i.e. for variations in population density N due to births and deaths (but not to dispersal). The logistic equation (10) describes an exponential growth for low values of population density, whereas it is self-limiting for large densities, saturating at N max . Note that the limiting term (within brackets) in equation (10) is similar to expression (4) for the available space that we have used in the previous section. Therefore, one can say that population growth, according to equation (10), is limited by the fraction of available space [12]. Now, equation (10) corresponds to a single population reproducing without external competition. But when we have a second population using the same space and resources, the presence of this additional population must also contribute to limiting the growth process. Thus, we can modify equation (10) so that the growth function of the Neolithic population, N , also includes the effect of the fraction of space occupied by Mesolithic populations, M [12]. A population density M occupies a fraction M/M max of the space available, in addition to that occupied by N . Therefore, N within parentheses in equation (10) should be replaced by (N + (M/M max )N max ). Then, Growth functions similar to (11) have been applied to competing micro-organisms [12].

Evolution equation
We can describe the variation in the Neolithic population density N , during a generation time T , as the sum of the variations due to the dispersion process and due to population growth (see [8] for details), where, as in equation (10), the subindex g stands for population growth (as opposed to dispersal, which corresponds to the first two terms on the right-hand side). Now, if we Taylor-expand equation (12) up to first order in time and to second order in space, we find that We could also have Taylor-expanded equation (12) up to second order in time [8], finding in this case slightly lower speeds; however, the conclusions we find here would not change.
The direction-dependent diffusion coefficients D x and D y , for our kernel (9), are where T is the generation time, and the mean value of a variable, for example 2 x , is defined as The advection terms U x , U y and U x y , using the kernel (9), are 6 As could be expected from the fact that the jump probability distribution (9) depends only on y = sin θ but not on x , we have obtained advection only in the y-direction, equation (18).

Front speed
As usual, we apply that for t → ∞ the front can be considered locally planar; thus for y → ∞ we can consider the variation in the x-direction negligible [13], and the evolution equation (13) for y → ∞ becomes where we will use equation (11) for the growth function F(N ).
As usual [13], we look for constant-shaped solutions to equation (20) with the form N = N 0 exp[ − λ(y − ct)] for N 0. As the Neolithic population density N at the leading edge of the front is low, F(N ) in equation (11) can be linearized, and we obtain from equation (20) In order for λ to be real, the term within the square-root must be non-negative, so the front speed c is Equation (22) can be also obtained, without the need for equation (21), by noting that equation (20) is simply Fisher's equation with (i) a modified growth term (11), which after linearization leads to a modified initial growth rateã = a(1 − m(y)), and (ii) an advection velocity v = 2D ∂m/∂ y 1−m(y) . Thus, the speed of front solutions to equation (20) must be Fisher's, namely 2 √ a D, minus the advection velocity, in agreement with equation (22). Moreover, the front speed obtained from the linear analysis described above is a lower bound to the front speed c. However, it is easy to apply variational analysis [14] and derive an upper bound with the same result, so equation (22) is the exact front speed for equation (20).
From equation (22), we see that if the Mesolithic population density increases with y, then the front speed decreases because of two effects: (i) the higher the gradient of the reduced Mesolithic density m, the higher the correction on the front speed; (ii) the speed also changes if there is less available space for the Neolithic population, i.e. for lower values of s = (1 − m(y)) (if s = 1, this second effect disappears).
To see the actual behavior of the front speed, equation (22), we need an expression for the variation of the Mesolithic population density M with y. However, the precise function M(y) is unknown because published data on Mesolithic settlements are scarce and restricted to very specific local areas, and also because the estimation of population densities from archaeological data relies on assumptions that are difficult to test and cause important methodological problems [15]. However, as explained in the introduction, we do know that the Mesolithic density did increase at northern latitudes [7]. Thus, we apply equation (22) To estimate the anthropological parameters a and D appearing in equation (22), we apply that the initial growth rate for preindustrial populations has a mean value of a = 0.028 year −1 [10], the mean-squared jump distance is 2 = 1531 km 2 [8] and the mean generation time is T = 32 years [16] 2 As we expected from equation (22), we see in figure 3 that each of the four test functions leads to a decrease in the front speed along the y-direction. To better compare the results with archaeological data, in figure 3 we have plotted c/c max , where the maximum speed from equation (22) is given by Fisher's value [11], In fact, this should be corrected due to a time-delay effect [8]. This would further complicate our equations, so we will not include this effect because, rather than comparing to the absolute . Symbols: observed front speeds calculated from archaeological data [9].
value of the maximum front speed (which we already analyzed in [8]), here we are interested in focusing our attention on the slowdown effect. This is simpler by considering the relative speed from equation (22), In figure 3 we compare the results obtained from equation (25) (curves) with Neolithic front speed data (symbols). The latter was obtained by computing the areas within isochrones separated 250 years inside the rectangle in figure 1 (such isochrones are shown in figure 1 every 500 years for clarity) 3 . Comparing the results from equation (25) to those from archaeological data in figure 3, we see that, even though none of the four test functions reproduce exactly the behavior of the archaeological data (which is not surprising for such a complex phenomenon), they do give a good approximation to the general behavior (especially m 4 ). Thus, a simple physical model can explain qualitatively the decrease in the front speed during the Neolithic expansion range in Europe. Therefore, physical models are useful to explain not only the average Neolithic front speed [8], but also its gradual slowdown in space.
The reaction-diffusion model presented in this paper could be applied to many examples of invasion fronts in which the indigenous population and the invasive one compete for space in a single biological niche, both in natural habitats [17,18] and in microbiological assays [12,19].