A neutral model of species’ replacement

We suggest that a model of migration and neutral species drift can explain the replacement of Neanderthals by Moderns and is in line with the evidence to date. Our model assumes no selective differences between the two species; that is, the competitive interaction between individuals or groups from the different species is identical to the competitive interaction between individuals or groups within the same species. Thus the two species are equivalent to two non-interbreeding subgroups of a single species. Although the two species did interbreed to some extent4, 30, 48, for simplicity we do not incorporate introgression into our model. In it, the only trait of interest is the species’ identity of individuals or groups; this formally equates the model with a simple, well-studied scenario: two selectively neutral alleles segregating at a genetic locus32. Perhaps the most fundamental property of such a scenario (in the absence of mutation) is that random drift will ultimately lead to the fixation of one allele and the extinction of the other. Applied to species, the analogous process has been termed “species drift”76.

This portrayal of the Neanderthals–Moderns situation is already sufficient to explain why one of the species had to eventually disappear, and is in line with the archaeological evidence that points to a period of co-occurrence of the two species in Europe and the Levant. However, in order to understand how and why the two species’ history would necessarily result in the Neanderthals’ extinction, we must take into account geographic and demographic aspects of the two species’ populations at the time. To do so, we model Europe and the Levant (deme 1, for simplicity, referred to henceforth as “Europe”) and Africa (deme 2, “Africa”) as separate demes with migration between them. The two demes have constant but possibly different hominin carrying capacities. See Supplementary Note 3 for details.

For realism and simplicity, we consider the dynamics of bands of individuals. That is, the entities whose fate is tracked in our model are small groups of individuals: such a band may die out by chance and be replaced by a propagule from another band (similar to the “propagule pool” model described in ref. 83, see also refs. 34,35,36, a propagule should be regarded as a copy of its band of origin). The carrying capacity of bands that reside in deme x (x = 1,2) is the constant N x . The probability of outgoing migration of a propagule from deme x per time step is denoted M x , and is proportional to a parameter m x and to N x . The rate of migration is assumed to be small enough that at most a single propagule can migrate per time step; accordingly, if m x ⋅ N x > 1, we set M x = 1, in which case migration occurs with probability 1 at every time step.

The population dynamics are those of a birth-death process akin to a Moran process with migration: at every time step a band chosen at random (regardless of its species’ identity) dies out, and is randomly replaced by a propagule from one of the other bands in its deme (Africa or Europe) or by a migrant propagule that had arrived from the other deme during the most recent time step. These stochastic death-and-replacement dynamics are a result of local extinctions due to environmental fluctuations, stochasticity in reproductive success, or competition among bands. This competition, which is assumed to be identically harsh between any two bands regardless of their species’ identity, may include inter-group competition for resources (e.g., territory, food, dwelling sites) and even direct inter-group violence. Importantly, the result—which of the competing bands will die out—is independent of the bands’ species’ identity (see also Supplementary Note 3, section 4). Note that one implication of these dynamics is that the population within each deme is fully mixed, i.e., this model does not account for the geographic spatial structure within each deme. Below we present a simple model that considers spatial structure.

We use the term “establishment” to describe the case in which a propagule migrated and was chosen to replace a band that died out. The probability of a migrant propagule’s establishment in deme x after arriving from deme y is thus 1/N x . Only propagules of existing bands migrate; thus migration has no effect on the population in its deme of origin. These dynamics ensure that the “population size” in each deme is constant, equal to the “carrying capacity” of that deme, N x . Except where noted otherwise, both terms refer henceforth to the number of bands in a deme.

The following set of transition probabilities constitutes a full mathematical description of this model. Let i x denote the number of bands of Moderns in deme x, and \({P_{{i_x} \to {i_x} + 1}}\) denote the probability per time step, at time t, of an increase by 1 of i x , \({P_{{i_x} \to {i_x} - 1}}\) the probability of decrease by 1 of i x , and \({P_{{i_x} \to {i_x}}}\) the probability of no change to i x . As indicated by the lower subscript of i x , Eqs. (1)–(3) describe the dynamics in deme 1, representing Europe, and equations for the dynamics in deme 2, representing Africa, are obtained by replacing every 1 by 2 and 2 by 1 in Eqs. (1)–(3).

$${{P_{{i_1} \to {i_1} + 1}} = \frac{{{N_1} - {i_1}}}{{{N_1}}}\left[ {\left( {1 - {M_2}} \right) \cdot \frac{{{i_1}}}{{{N_1} - 1}} + {M_2} \cdot \left( {\left( {\frac{{{i_2}}}{{{N_2}}}} \right) \cdot \frac{{{i_1} + 1}}{{{N_1}}} + \left( {\frac{{{N_2} - {i_2}}}{{{N_2}}}} \right) \cdot \frac{{{i_1}}}{{{N_1}}}} \right)} \right]}$$ (1)

$${{{P_{{i_1} \to {i_1} - 1}} = \frac{{{i_1}}}{{{N_1}}}\left[ {\left( {1 - {M_2}} \right) \cdot \frac{{{N_1} - {i_1}}}{{{N_1} - 1}} + {M_2} \cdot \left( {\left( {\frac{{{i_2}}}{{{N_2}}}} \right) \cdot \frac{{{N_1} - i_1}}{{{N_1}}} + \left( {\frac{{{N_2} - {i_2}}}{{{N_2}}}} \right) \cdot \frac{{{N_1} - {i_1} + 1}}{{{N_1}}}} \right)} \right]}}$$ (2)

$${P_{{i_1} \to {i_1}}} = 1 - \left( {{P_{{i_1} \to {i_1} + 1}} + {P_{{i_1} \to {i_1} - 1}}} \right)$$ (3)

These equations are derived as follows: In Eq. 1, \(\left( {\frac{{{N_1} - {i_1}}}{{{N_1}}}} \right)\) is the probability that in the current time step, a Neanderthal band in deme 1 dies out (otherwise an increase in the Modern’s population in this deme during this time step is impossible). The terms \(\left( {1 - {M_2}} \right)\) and \({M_2}\), respectively, represent the probabilities that migration from deme 2 did not occur and that it did occur. The term \(\frac{{{i_1}}}{{{N_1} - 1}}\) is the probability that the propagule chosen to replace the one that died out in deme 1 (Europe) is Modern, given that no migration occurred; thus the number of candidate propagules that can act as a replacement is N 1 −1. \(\frac{{{i_2}}}{{{N_2}}}\) and \(\frac{{{i_1} + 1}}{{{N_1}}}\) represent, respectively, the probabilities that the migrant propagule to Europe is Modern, and that a Modern propagule is chosen to replace the band that died out. Another possibility that increases the Modern population in Europe, given that migration had occurred, is represented by \(\left( {\frac{{{N_2} - {i_2}}}{{{N_2}}}} \right) \cdot \frac{{{i_1}}}{{{N_1}}}\), i.e., the migrant to Europe is Neanderthal, and yet a Modern propagule is chosen to replace the band that died out. Eq. 2 is composed of analogous constituents, whose interpretation is analogous to the description above. This model is similar (but not identical) to the Moran process with mutation that is studied, for example, by Ewens (ref. 32, p. 106), if one of the migration probabilities is zero. The expected number of time steps required for species’ exclusion, in the process described by Eqs. (1)–(3), can be calculated using equations 2.144 and 2.160 in ref. 32, and is not provided here.

We limit the scope of our exploration to conditions in which one of the demes (deme 1, representing Europe) is initially populated only by bands of Neanderthals, and deme 2, representing Africa, is initially populated only by bands of Moderns. We first analyze the case in which migration occurs only from Africa to Europe (deme 2 to deme 1), the scenario that is widely believed to have taken place near to and during the interaction between the Neanderthal and Modern populations, based on the lack of evidence so far that would support Neanderthals’ existence in Africa6,7,8,9,10,11,12. For completeness, we then report simulation results for the case in which migration occurs in both directions, with both equal and unequal migration rates in the two directions.

Methods for comparing model results to empirical evidence

The replacement of Neanderthals by Moderns seems to have occurred surprisingly fast when compared to archaeological and evolutionary time scales of the two species’ existence. This may be why many scholars assume that the process was necessarily driven by selection, but whether the process should be considered as having been rapid depends on properties of the model assumed, such as the Neanderthal population size and the expected duration of replacement. Our model, which takes into account major aspects of the two species’ demography at the time but does not include selection, may be regarded as a baseline, or a null model, for such an evaluation. To assess whether this null model can be rejected in favor of a selection scenario, we study the distribution of replacement durations produced by the model and compare it to the replacement duration suggested by empirical evidence. Such an attempt faces a number of obstacles.

First, it is notoriously hard to correlate the time units in evolutionary models with the time span of real-life scenarios. This is also the case in our model, which entails the selectively neutral replacement of bands, whose empirical rate is hard to gauge and for which there are no clear estimates. Second, the species’ replacement duration in our model depends critically on the parameter N 1 , the number of bands in Europe. Estimates of hominid population sizes in Europe near the end of the middle Paleolithic vary over more than an order of magnitude and—according to reconstructions of paleo-climate during this era—are likely to have changed significantly while the replacement was taking place15. Finally, the way in which replacement duration is estimated may affect the result by more than an order of magnitude. Comparison between the simulations and the archaeologically estimated period of coexistence of the two species should take into account the time point at which species’ coexistence is likely to be evident in the archaeological record. That is, the appropriate duration to be compared should not be the full duration of each model simulation, but the period during which both species are likely to have a demonstrable presence. This would be based on archaeological findings that can be clearly associated with the identity of the species that produced them, and could, for example, be the period from the initial crossing of some frequency threshold by the Moderns until the Neanderthals constitute less than this threshold in the overall European population, or the period between the last crossings of these thresholds (see Supplementary Note 1 for a discussion of alternatives).

Accordingly, we have conducted statistical analyses for a range of combinations of parameters. We attempt to assess, under a range of possible choices and assumptions, whether the time required for species replacement according to our model is within a reasonable range compared to the empirical evidence. Because our model serves as a null model relative to models of species replacement that are not neutral (for example, that assume selection), we test, for each combination of parameters, whether our model can be rejected at the p = 0.05 significance level. For this, we ask whether the archaeologically supported duration of species coexistence falls, for example, in or below the range of the 5% shortest coexistence durations produced by our model, in which case our model should be rejected, or outside of this range (if one wishes to conduct this test with a lower rate of false rejections, i.e., using a more stringent p-value, the 5% threshold should be decreased accordingly). In our analysis, the considered duration of species coexistence is the time from the last simulation time step in which Moderns comprised 10% of the population in Europe until the last time that they reached 90% of it (see Supplementary Note 1 for a discussion of this choice and alternatives). The rates are described in units of the probability of band replacement per generation (25 years).

Each of our model’s time steps is equal to a stochastic replacement of one band by another. The rate of band replacement, the characteristic band size, and the hominid carrying capacity of Europe in the middle Paleolithic, are variables that are not necessary in order to derive the numerical results of our model, but an estimate of these is required, post-simulation, in order to compare the results to the period of species’ coexistence in the archaeological record and determine whether the model should be rejected. We conduct simulations for a broad range of carrying capacities in Europe, ranging from 10 to 500 bands. This covers the full range of estimates suggested in the literature, with the size of bands, or identity-groups, ranging from 50 to 100034,35,36,37,38,39,40 and with population sizes ranging from 5000 to 70,00037, 33, 15. Considering the rate of band replacements in units of the per-generation probability that a band is replaced, we determine whether our model should be rejected for the full range of possible values of this variable, with the per-generation replacement probability ranging from 0 to 1.

A spatial model of species’ replacement

Any model, particularly one such as ours that aims to describe the null expectation, is a highly simplified version of reality. The aspect in which this simplification deviates most significantly from the reality of the inter-species dynamics is that it does not account for the spatial structure of the dynamics beyond the split into two major demes. A full treatment of the spatial complexity of the dynamics requires complex modeling and a large number of strong assumptions, making the model less general and more sensitive to specific details, over which there is no agreement among researchers. Such treatment is well beyond the scope of our current study. However, to check whether spatial structure would lead to qualitatively different results from those of the non-spatial model, we have examined a simple model that captures the effect of a spatial component in the inter-species dynamics. We designed it such that it would be easily tractable and intuitive in its structure, and would require a minimal number of assumptions.

We model the spread of Moderns into Eurasia as a stepping-stone model84, where the inter-species dynamics play out on a shifting front of contact, i.e., by reducing the 2-D spatial structure of Eurasia (henceforth “Europe”) to a 1-D series of single-band territories, as illustrated in Fig. 7, where a species range increases or decreases via replacement of the foremost band of one species by the other. As previously, we consider a Moran model, where at each time step one band stochastically dies out (chosen regardless of its species identity, reflecting the null assumption that the species are selectively equivalent) and is replaced by a randomly chosen neighboring band. As illustrated in Fig. 7, in this model, each band has only two neighboring bands.

Fig. 7 An illustration of the 1-D stepping-stone model. Each deme is occupied by a single band, whose species identity is denoted by M (Modern) or N (Neanderthal). When a band stochastically dies out, it is replaced by a propagule from one of its immediate neighbors. Migration of Moderns from Africa occurs to the rightmost deme within Europe, and is assumed to occur immediately if the band in this deme dies out Full size image

We use this model to study only the scenario of unidirectional migration out of Africa. We focus on the duration of co-habitation of parts of Europe by bands of the two species; the period of interest always begins once both species are found within Europe, and so the migration dynamics need not be modeled explicitly. This is realized by assuming that the first (henceforth, “rightmost”) stepping stone (henceforth, “deme”) would be habitable only by Moderns; Neanderthals cannot establish in this deme, and if the band in it happens to die out, it is immediately replaced by a Modern band, representing a successful establishment of a migrant band from Africa. The initial conditions of each simulation run are that all other demes are initially occupied by Neanderthals.

As in the model of unidirectional migration described previously, each simulation is determined to end in full Neanderthal replacement, because the system has a single attracting state: the rightmost deme is always populated by Moderns, which may spread from there leftwards via drift but can never go extinct, while Neanderthals initially populate all other demes. If the drift dynamics gradually lead to all of the demes being populated by Moderns, the Neanderthal’s population cannot be replenished, and they go extinct. The question of interest is primarily how long this process of neutral species replacement via random drift is expected to be. Because the average time between band replacements is a parameter in the analysis, which should be derived from empirical evidence and that does not affect the model’s result when time is measured in units of mean number of band replacements per territory, this question is analogous to the question of how many times, on average, each deme changes hands between bands until full replacement. In addition, it is particularly interesting to explore how many of these band replacement events are within-species and how many are between species.

The spatial structure and the initial conditions described above lead to dynamics in which there is always a single front of interactions between the two species, to the right of which all demes are populated by Moderns, and to the left—by Neanderthals. The dynamics of the front of species interaction follow those of a random walk in 1-D with a single-absorbing state.

The number of demes in the model must reflect the carrying capacity of Europe, in terms of bands, transformed into a linear sequence of neighboring territories. This is the model’s primary point of weakness: it is unclear what should be the correct transformation in order to best simulate reality in this respect. One extreme possibility would be to assume that the length of the series of demes equals the total number of bands estimated to exist in Europe. This is analogous to assuming that the spread of Moderns into all of Europe had to occur strictly along a single trajectory, as illustrated in Fig. 8a. A more realistic assumption would be that the model captures one of multiple such trajectories that occur in parallel (Fig. 8b; also a significant simplification of reality, as discussed previously). The results would be highly sensitive to the number of these assumed independent trajectories. Since the correct choice of this number is unknown, we do not rigorously analyze the space of parameter combinations. Instead, we explore whether a range of reasonable choices for the number of demes leads to replacement durations that are of a similar order of magnitude to those found for the non-spatial model described previously, and compare the dynamics of inter-species replacement between the two models in terms of the average number of band replacements in which a band of one species replaces a band of the other.

Fig. 8 Alternative models of Moderns’ spatial trajectory in Europe. Schematic illustrations of modeling alternatives that differ in the number of trajectories along which Moderns may spread into Europe. One trajectory in a, 10 trajectories in b. Each trajectory is modeled as a single 1-D series of demes. The model realizes the dynamics along a single trajectory Full size image

In the reported simulations, we set the number of independent trajectories across Europe, somewhat arbitrarily, to ten (Fig. 8b), thus determining that the length of the stepping-stone chain for this model to be a tenth of Europe’s carrying capacity (measured in the number of bands it supports).

Code availability

All simulation code and scripts for further analysis are available at https://github.com/orenkolodny.

Data availability

All data supporting the findings are found in the paper and in its Supplementary Notes. Unprocessed simulation outputs are available from O.K. upon request.