Biased dispersal can explain fast human range expansions

Some human fronts spread faster than expected by models based on dispersal and reproduction. The only explanation proposed so far assumes that some autochthonous individuals are incorporated by the expanding populations, leading to faster front speeds. Here we show that simple models without this effect are also consistent with the observed speeds of two fronts (a Khoi-khoi expansion of herders and a Bantu expansion of farmers), provided that the dispersal of individuals is biased (i.e., more probable) in directions closer to the front propagation direction. The physical models presented may also be applied to other kinds of social phenomena, including innovation diffusion, rumor propagation, linguistic fronts, epidemic spread, diffusion in economic space and the evolution of cooperation in spatial systems. They can be also adapted to non-human systems with biased dispersal, including biological invasions, cancer tumors and virus treatment of tumors.


Biased dispersal in human population fronts
and Kolmogorov et al. 5 simultaneously derived the following reaction-diffusion equation 2 where t is the time, N the population density (of Neolithic individuals in our case), D its diffusion coefficient and F N ( ) its growth function, i.e. the net reproduction (births minus deaths per unit time). Usually a logistic function is used for F N ( ) because it is well-known to agree with many human population data (see 16 and references therein), i.e., where a is the initial growth rate and K the carrying capacity (maximum population density). Fisher's Eq. (1) can be obtained by performing Taylor expansions in the so-called non-cohabitation equation [Eq.
(3) below], which is derived by applying simply that the variation in population density at position x y ( , ) is due to individuals arriving minus those leaving plus the net reproduction, i.e 16 .
x y x y x y where T is the generation time (mean time interval between the migration of a parent and one of her/his children) 16,31 . The dispersal kernel φ ∆ ∆ ( , ) x y can be defined so that φ ∆ ∆ ∆ ∆ d d ( , ) x y x y is the probability that an individual born at − ∆ − ∆ x y ( , ) x y and time t will have a child born between (x, y) and + ∆ + ∆ x d y d ( , ) x y at time + t T (i.e., that a jump in the random walk will have x-coordinate between ∆ x and d x x ∆ + ∆ and y-coordinate between y ∆ and ∆ + ∆ d y y ) 16 . Also in Eq. (3), R N x y t [ ( , , )] T is the effect of logistic net reproduction, i.e 17,32 . Note that Eqs. (3) and (4) without dispersal yield the same result as integrating the logistic differential equation[Equations (1)- (2), also without diffusion], between times t and t T + . In the limit T 0 → , Eq. (4) tends to N x y t ( , , ), and so does the first term on the right-hand side of Eq. (3) because dispersal distances tend to zero, i.e., ( , ) x y φ ∆ ∆ tends to a Dirac delta centered at ( , ) (0, 0) x y ∆ ∆ = . By performing Taylor expansions up to second order in space and first order in time in the non-cohabitation Eq. (3) and assuming an isotropic kernel, it is easy to obtain the Fisher Eq. (1) 16 , which was applied by Ammerman and Cavalli-Sforza to model the spread of the Neolithic in Europe 6,7 . However, in contrast to many physical and chemical systems, it has been shown that for human populations the speed of front solutions to Fisher's Eq. (1) is not a valid approximation to the more precise Eq. (3) 17 . Moreover, newborn humans cannot survive separated from their parents during some years and for this reason, the non-cohabitation Eq. (3) should be replaced by the so-called cohabitation equation, namely 17,33,34  but the first term on the right of Eq. (3) implies that parents that have migrated away from x y ( , ) will have left their children there (see Fig. 1 in 17 ). This is appropriate for some biological species but not for humans 33 . In contrast, according to Eq. (5) parents and children will be at the same final position. It is important to note that in the cohabitation Eq. (5) we have assumed that parents reproduce after migrating, but the speed of front solutions is the same if they migrate after reproducing and even if they reproduce during migration (for a detailed discussion of this point, see ref. and 33 and Fig. 17 in ref. 3 ). It is well-established that for human populations it is necessary to use Eq. (5) rather than (3) or (1) 3,17,33 , so the rest of this paper is based on the cohabitation Eq. (5).
As shown in previous work 17,33 , the speed of front solutions to Eq. (5) can be derived by considering the leading edge of the front, i.e., assuming that < N x y t K ( , , ) in Eq. (4), choosing the positive x-axis along the local front propagation direction and using the ansatz N x y t N e ( , , ) (5), where 0 λ > and c is the front speed. Then 17,33 cT aT 0 cos where arccos x θ = ∆ ∆ is the angle between the individual displacement ( , ) x y ∆ ∆ and the front propagation direction ( > x 0). For later use, we have assumed that the dispersal kernel can be written as x y ( ) is the probability of migration as a function of distance ∆ and θ Φ( ) is the probability of migration as a function of direction θ. This is justified if the migration distance ∆ is statistically independent of the migration direction θ.
In recent years, the framework summarized above has been generalized by adding cultural transmission, i.e., an indigenous population (e.g., hunter-gatherers) the individuals of which can join the invading population (e.g., herders or farmers) 12 . In contrast, here we will introduce an alternative model that does not assume cultural transmission but, instead, allows for the possibility that individuals disperse with different probabilities in different directions. This seems reasonable because humans in a population invading a landscape will likely tend to migrate to places where the population density is still low, i.e., the probability of dispersal or migration will increase for directions closer to that of the front propagation (lower values of θ) and this should lead to faster fronts than isotropic (i.e., non-biased) dispersal. It is important to distinguish the effect that we will analyze in this paper (biased dispersal) from a different effect that leads to similar equations and is called correlated dispersal [35][36][37][38][39] . In the latter case individuals also perform a random walk, but the probability of the direction of motion does not depend on its angle relative to a fixed direction (as in biased dispersal) but on the angle relative to the direction of motion during the previous step of the random walk (i.e., during the previous generation in human expansions). We do not deal with this case (correlated dispersal) because we see no reason to think that the direction between, e.g., the birthplace of an individual (A) and one of his children can be related to the direction between the birthplace of the first individual (A) and that of his mother or father. For a detailed comparison between biased and correlated random walks, see ref. 37 .
As in previous work 2,7,12,18,19 , we can use ethnographic data from preindustrial populations for the probability Ψ ∆ ( ) of migration as a function of distance ∆. Such data are usually reported as histograms, i.e., a set of distances and the corresponding probabilities, where r i is the mean distance of the i-th bin in the histogram, p i is its probability A simple way to derive the right-hand side of Eq. (8) is to consider Eq. (7) and the normalization conditions Unfortunately, in contrast to Ψ ∆ ( ), for preindustrial populations there are no data at present that allow to measure the probability θ Φ( ) of migration as a function of the angle θ between the direction of individual dispersal ∆ ∆ ( ) , x y and that of the front propagation (although such data will hopefully become available in the following years, as explained in the Conclusions section below). Therefore we will introduce and compare three increasingly complex models of biased (or non-isotropic) dispersal, i.e., three different functions ( ) θ Φ . This will be enough to fulfill our aim, namely to determine whether or not biased dispersal is a viable mechanism to explain the fastness of the Neolithic fronts of herders 18 and famers 19 mentioned above. Intuitively we expect that ( ) ( ) θ θ Φ − = Φ , i.e. that dispersal depends only on the magnitude of the angle θ between the dispersal jump and the front propagation direction (θ = 0), so the three models that we consider have this property.
Obviously geographical features (such as rivers and mountains) can introduce further effects, but the front speed is measured using data separated by very large distances (up to several thousand kilometers) and the high values of the correlation coefficients for the corresponding linear fits of times and distances (up to r 0 85 = . for the Khoi-khoi expansion of herders 18 and = . r 0 87 for the southern Bantu expansion of farmers 19 ) make it reasonable to develop models in homogeneous space. It is also worth to note that at present there are not enough data from independent observations to parametrize non-homogeneous models. In any case, here we are interested in the simplest possible physical models to answer the question that has motivated this study, namely whether non-isotropic dispersal can explain why the speed of some human range expansions is too fast to agree with demic models based on isotropic dispersal 18,19 .

Model 1.
A very simple model is to assume that individuals have probability p to migrate forward (i.e., outwards from the invaded range), namely with 0 , and with probability p 1 − to migrate backward (i.e., inwards to the invaded range), namely with 0 . This simple assumption has been previously applied in microbiology, specifically to the spread of some brain tumors in which cell dispersal is biased outwards from the tumor core 29 . However, in that case the non-cohabitation Eq. (3) was used because the progeny of tumor cells can disperse away from their parent cells immediately after reproduction (i.e., there is no cohabitation effect). In contrast, here we have to use the cohabitation Eq. (5) because we are dealing with human populations 3,17,33,34 . As explained in the Introduction, the problem that we want to solve is that some front speeds are faster than those predicted by isotropic dispersal 18,19 . For this reason, obviously we have to consider the case p 0 5 > . . Model 1 is shown in Fig. 1 for = . p 0 5 (isotropic dispersal) and = p 1 (which corresponds to the maximum possible dispersal anisotropy for this model). Using the normalization condition (10), we find that in model 1 the angular dispersal distribution is Using this distribution and Eq. (8) in Eq. (6), we find that the front speed in model 1 is is the modified Bessel function of the first kind and order zero, and is the modified Struve function of the fist kind and order zero. We have used Eqs. 3.339 and 8.551-2 in ref. 40 . We have also assumed that the minimum possible speed is that selected by the front, a well-established result not only for differential 2 but also for integro-difference 12 equations. In the isotropic case ( = p 1/2) Eq. (12) reduces to Eq. (14) in 33 or Eq. (9) in 17 , as it should. In Fig. 2 we apply Eq. (12) to the spread of Khoi-khoi herders in Southwestern Africa (2,300-1,100 yr BP), for which the observed speed is 1.4-3.3 km/yr 18 (horizontal hatched rectangle in Fig. 2a). This is a fast speed, if compared to that of the the spread of the most well-known case, namely the spread of the Neolithic in Europe (0.9-1.3 12 ). The parameter values, already used in ref. 18 , are . ≤ ≤ . a 0 023 0 033 yr −1 for the initial growth rate (estimated archaeologically and ethnographically for several small preindustrial populations that settled into empty space), www.nature.com/scientificreports www.nature.com/scientificreports/ bound (1.4 km/yr) and sufficiently biased dispersal, namely ≥ .
p 0 9. This range is shown as a vertical hatched rectangle in Fig. 2a,b. Figure 2b shows the effect of the dispersal bias, namely the speed c from Eq. (12) minus the corresponding speed for non-biased dispersal (p 1/2 = ), divided by the former and multiplied by 100. We see that according to model 1, consistency between the observed and predicted speeds (hatched vertical rectangle, i.e., ≥ . p 0 9) implies that the bias effect is between 18% and 23% (Fig. 2b). In Fig. 3 we apply the same model 1, but to a different range expansion which was also rather fast, namely that of Bantu farmers in southeastern Africa (3,400-1,400 yr BP), for which the observed speed southwards from the Great Lakes area is 1.5-2.3 km/yr 19 (hatched horizontal rectangle in Fig. 3a), using as in ref. 19 12,17,18 , see Supplementary Text S1). We see in Fig. 3a that model 1 is again viable, but only for p 0 9 ≥ . (vertical hatched rectangle) and sufficiently slow speeds. Interestingly, this is the same conclusion as for the Khoi-khoi expansion (Fig. 2a), in spite of the fact that the dispersal kernels are different. For the Bantu expansion the bias effect is 19-25% (Fig. 3b), very similar to the range 18-23% for the Khoi-khoi (Fig. 2b).

Model 2.
Model 1 is consistent with the Khoi-khoi and Bantu data, but only assuming that the speed was very close to its minimum possible value suggested by the archaeological data (Figs. 2a and 3a, dashed and dashed-dotted curves). It is thus reasonable to ask if a different non-isotropic model can agree with a wider range of the observed speed. Moreover, although model 1 is very useful conceptually, it has the obvious pitfall of assuming the same probability (p) for all forward directions ( π θ π − ≤ ≤ /2 /2) and, similarly, the same probability ( p 1 − ) for all backward directions (see Eq. (11)). There is thus a discontinuity at θ π = /2 (as seen in Fig. 1, curve p 0 55 and = p 1). For model 3, the bias parameter p has no upper limit and the case p → ∞ corresponds to a Dirac delta centered at θ = 0, i.e., all individuals moving in the front propagation direction. (2020) 10:9036 | https://doi.org/10.1038/s41598-020-66045-2 www.nature.com/scientificreports www.nature.com/scientificreports/ for model 1 and p 1 = ), which is a questionable feature intuitively. We introduce model 2 to avoid this limitation as follows, This model 2 has been used previously in human dispersal 34 but only in the diffusive approximation (i.e., a diffusive coefficient or second-order moment instead of the complete dispersal kernel), and that approximation is nowadays known to be invalid for human populations 17 . As in model 1, the case = p 1/2 corresponds to isotropic or non-biased dispersal (uniform distribution θ π Φ = ( ) 1/(2 )) and the case = p 1 to the maximum possible bias allowed by this model (because we must have ( ) 0 θ Φ ≥ for all values of θ). Model 2 (Eq. (13)) is shown for p 1/2 = and = p 1 in Fig. 1. In this figure, comparison of model 2 to model 1 (both of them for p 1 = ) shows that  18 . For each model, the range between the two curves is its predicted front speed. It is seen that models 1 and 2 are consistent with the observed front speed only if the latter is very close to its lower bound (1.4 km/yr). In contrast, model 3 is consistent with the observed front speed for a substantially wider range of the latter (up to 1.89 km/yr). This happens if p 0 53 ≥ . (this consistency range is not shaded for clarity), which implies that the bias does not need to be very strong (see Fig. 1 and the text below Eq. (17)). (b) Effect of the dispersal non-isotropy on the same Khoi-khoi front speed. For each value of the dispersal bias p, this effect is defined as the difference between the corresponding front speed minus the speed for non-biased dispersal ( = . www.nature.com/scientificreports www.nature.com/scientificreports/ model 2 is much more acceptable intuitively, because the discontinuity is no longer present. We next determine how using model 2 instead of model 1 changes the conclusions concerning the compatibility of predicted and observed speeds for the Khoi-khoi and Bantu Neolithic expansions. Using the distribution (13) is the modified Bessel function of the first kind and order one, and we have used Eqs. 3.339, 3.915-2, 8.406-3 and 8.447-2 in ref. 40 . Figure 2a,b include, for the same parameter values as for model 1 above, the results of model 2 (full and empty circles). It is seen that the speeds predicted by model 2 are very similar to those by model 1. Thus, interestingly, Figure 3. (a) Shows the speed of fronts for the Bantu expansion that transformed southeastern Africa from a land of hunter-gatherers into one of farmers from about 3,400 until 1,400 yr BP. The observed range for the front speed, estimated from archaeological data, is shown as a horizontal hatched rectangle (1.5-2.3 km/yr) 19 . For each model, the range between the two curves is its predicted front speed. In spite of the fact that the dispersal kernel is different from that used in Fig. 2 (because here we are dealing with farmers rather than herders), the conclusions are very similar. Indeed, models 1 and 2 agree with the observed front speed only if the latter is very close to its lower bound (1.5 km/yr in this case). In contrast, model 3 is consistent with the observed front speed for a substantially wider range of the latter (up to 1.96 km/yr), which happens if ≥ .
p 0 53, so the bias does not need to be very strong (see Fig. 1 and the text below Eq. (17)). (b) The effect of the dispersal bias p on the front speed of the same southern Bantu expansion is 19-25% for model 1 ( ≥ . p 0 9), 18-23% for model 2 (same consistency range of p) and 19-42% for model 3 ( ≥ . p 0 53).

Model 3.
Both models 1 and 2 yield predicted speeds similar to those observed only if the latter are close to their minimum possible values (lower side of the hatched horizontal rectangle in Fig. 2a,b). As mentioned at the beginning of the previous subsection, we are interested to know if there is some other biased-dispersal model that agrees with the observations for a wider range of the front speed. Intuitively, it is easy to note a drawback of both models 1 and 2, namely that they cannot describe populations for which the bias is arbitrarily large. This can be seen from Fig. 1, where for the maximum possible bias ( = p 1) we see that a large portion of individuals move in directions that are not close to that along which the front propagates (θ = 0), and for model 2 even towards the space that has been already invaded by other individuals ( /2 θ π < − or /2 θ π > ). This is unlikely to describe all possible non-isotropic dispersal behaviors, because it is perfectly reasonable that individuals of an invading human population may migrate only in directions close to that of the invasion front (θ ≈ 0). For this reason, we introduce model 3 with a Gaussian dispersal distribution, where we have have used θ 2 rather than θ in the exponent to ensure that ( ) θ Φ − = ( ) θ Φ (see the text below Eq. (10)). The normalization condition (10) , where x erf( ) is the error function and we have used Eqs. 3.321,2 in ref. 40 . In order to compare to models 1-2, it will be useful to introduce the following bias parameter = + . p q 10 20 (16) As in models 1-2, the case p 1/2 = (or = q 0) corresponds to isotropic dispersal. Other forms of p with this property are possible, but Eq. (16) has the advantage that we will be able to compare the speeds from models 1-3 in a single figure. Before doing so, note from Fig. 1 that the Gaussian distribution (15) (model 3) can describe situations in which many individuals move close to the front propagation direction (θ ≈ 0), whereas models 1-2 cannot. It is also seen that higher values of p correspond to more strongly biased dispersal (i.e., a larger fraction of individuals migrating in directions close to 0 θ = ). This fraction can be arbitrarily large in model 3, and this corresponds to the fact that, in contrast to models 1-2, there is not any upper limit to the value of p in model 3.
Using the distribution (15) and Eq. (8) in Eq. (6), we find that the front speed in model 3 is In contrast to models 1-2, it does not seem possible to solve this integral analytically. Therefore, we follow a different approach. For a given set of parameter values (a, T , {p j } and {r j }), we plot the quotient in the right-hand side of Eq. (17) and find its minimum numerically. In this way, we have computed the speeds for model 3 in Figs. 2a and 3a. We see that model 3 is consistent with the observed speed for a substantially wider range of the latter (up to about 1.9 km/yr) than models 1-2, and that this happens for > .
p 0 53, both for the Khoi-khoi expansion of herders (Fig. 2a) and the Bantu expansion of farmers (Fig. 3a), in spite of the fact that both populations have different dispersal kernels. Remarkably, p 0 53 = . is not an extremely isotropic case but, as shown in Fig. 1, it is intermediate between = .
p 0 5 (non-biased case) and p 0 55 = . (which is mildly biased, in the sense that some individuals move even backwards, i.e. with /2 θ π < − or θ π > /2, in contrast to, e.g., model 3 with p 1 = in Fig. 1). Thus we conclude that biased dispersal can explain the speed of both Neolithic fronts. According to model 3, the bias effect is 17-41% for the Khoi-khoi (using the consistency range above > .
p 0 53 in Fig. 2b) and 19-42% for the Bantu ( > . p 0 53 in Fig. 3b). Finally we note that in the most extremely biased case (p→∞), the angular distribution (15) becomes a Dirac delta centered at θ = 0 and Eq. (17) reduces to the new, very simple result conclusions In this paper we have shown that when a population expands its range, the front speed can be substantially faster if the dispersal of individuals is biased towards the front propagation direction, which seems a very reasonable assumption because individuals will presumably tend to move towards places with more free space (lower population densities) than to places with less free space (higher population densities). Moreover, for two human expansions (one of herders and one of farmers) we have shown that this effect can explain the range of the front speed estimated from archaeological data.
The present paper has crucial implications on a very important issue in Archaeology, Anthropology and Genetics, namely the relative importance of demic and cultural diffusion. For the two case studies considered in www.nature.com/scientificreports www.nature.com/scientificreports/ this work (Figs. 2-3), it has been previously shown that cultural transmission (i.e., the conversion of hunter-gatherers into farmers) 12 is a possible explanation of the fact that the observed front speeds are faster than those predicted by an isotropic, purely demic model ( = . p 0 5 in Figs. 2-3) 18,19 . However, here we have shown that biased dispersal may also explain the speeds of both range expansions (Figs. 2-3, especially model 3). Thus understanding Neolithic front speeds is by no means trivial, in the sense that they have these two possible explanations (and perhaps many others, see Supplementary Texts S1 and S2). Here we have shown that at least two different explanations may agree with the observed data. We would like to stress, however, that both explanations may partly apply, with different relative importance depending on the front considered. The problem of estimating the relative importance of biased dispersal and cultural transmission may be solved in the future as follows. It is highly likely that in some years it will become possible to measure dispersal kernels of prehistoric individuals directly (not using ethnographic analogues as in this paper and all previous work), because ancient genetics has been already used to detect a few prehistoric parent-child pairs buried at different places 41,42 . When many such pairs are detected, it should become possible to determine the anisotropy in the dispersal kernel (e.g., parameter p in our 3 models) in addition to the dispersal distances and probabilities (r i and p i in Eq. (8)). Then, the difference between the observed front speed and that predicted by a purely demic model with biased dispersal (e.g., models 1-3 in this paper) would be an estimation of the importance of cultural transmission (i.e., the conversion of hunter-gatherers into herders or farmers). This will hopefully make it possible to disentangle the roles of both complex mechanisms (cultural transmission on one hand, and biased dispersal on the other) in human range expansions.
The interest of our results is not limited to Neolithic inland spread, because there are also some data on Neolithic coastal 13,15 , Paleolithic 43 and modern 8 front speeds for which the dispersal bias effect could be of importance. Besides human range expansions, our approach can be applied to other social applications of statistical physics including innovation diffusion 22 , rumor propagation 23 , linguistic fronts 24 , diffusion in economic space 26 and the evolution of cooperation in spatial systems 27 , because it seems reasonable to expect that such human systems may exhibit higher dispersal probabilities in directions closer to that of the front propagation (for the reason given at the beginning of this section). Our 3 models are easy to adapt to some other fields for which the cohabitation Eq. (5) must be replaced by the non-cohabitation one (3), including ecology (biological invasions 28 and epidemic spread with biased dispersal due, e.g., to wind or water currents 25 ), microbiology (growth of bacterial colonies and cancer tumors displaying biased dispersal 29 ), medicine (virus treatment of tumors 30 ), etc.

Data availability
No datasets were generated or analyzed during the current study.