Three recent papers consider the role of climate change in Neanderthal replacement. Timmerman (2020), on the basis of recent simulation experiments, suggests that Homo sapiens must have been “considerably more effective in exploiting scarce glacial food resources as compared to Neanderthals.” Timmerman found that in simulation at least, Homo sapiens must be able to sustain a higher population than Neanderthal in the same environment (i.e. have a higher carrying capacity) in order to produce the kind of rapid replacement we find in the archaeological record.

Timmerman doesn’t specify the source of that advantage but mentions the usual suspects: Technological innovation (Aurignacian vs Mousterian), greater mobility, larger, less fragmented groups (social advantage). Environmental change does play a role in the explanation, but it has be abrupt, and allows the faster growing, more mobile, and more adaptable sapiens population to replace Neanderthals in the wake of sudden regional declines in carrying capacity.

Meanwhile, Raia et al. (2020), using paleoclimate models, find that our various cousin branches became extinct mostly as a result of climate change, but that in the case of Neanderthals, it was “probably exacerbated by competition with H. sapiens”. Haws et al. (2020), focusing on Iberia, also conclude that environmental change was the major driver for sapiens expansion. It is a familiar situation in our field to have environmentally driven explanations competing with culturally driven ones.

Timmerman looks at species competition in a very constrained simulation of the real-world geography, climate change and archaeological record of the sapiens expansion. Likewise, Raia et al. use a highly specified real world climate model.

Building on those experiments, I take a more general look at the overall dynamics between population growth rate, mobility, effectiveness at mobilizing carrying capacity, social connectedness, and rate and magnitude of environmental change, in the context of population competition and replacement. What are the parameters under which replacement happens? What conditions make it more or less likely?

The simulation

I use Netlogo to build a simple environmental change and population expansion and diffusion model. Patches have a carrying capacity and can host two kinds of populations (APop and BPop). Users can set the roughness of the initial environment by assigning a probability that a patch will vary from the maximum carrying capacity for the environment, and the degree by which patches vary. They can set the probability that a patch will undergo carrying capacity change in a time step, and the maximum amplitude of that change.

Users can set the growth rate, carrying capacity multiplier, diffusion rate, network radius, and minimum number of individuals in a viable network for each of APop and BPop, as well as the initial probability that a patch will be populated.

Each time step, the carrying capacity of each patch has a given probability of increasing or decreasing up to a maximum proportion. Each time step, APop and BPop grow at a specified rate. If the combined population of APop and BPop on a patch reach carrying capacity, both populations search neighbouring patches within a given network radius for available carrying capacity and move accordingly. I refer to this as overflow.

Independent of spread through overflow, APop and BPop also have diffusion rates. Even if a patch is not at maximum carrying capacity, a proportion of APop and BPop will search for available carrying capacity within their network radius and move accordingly.

At the end of a time step, any population over the carrying capacity in a given patch is removed. APop and BPop can use carrying capacity differently through a carrying capacity multiplier. For example, a carrying capacity multiplier of 2 for APop means that APop is multiplied by 2 for purposes of carrying capacity calculations. In other words, APop can sustain half the population that BPop could on the same patch.

Additionally, the population of any patch whose APop or BPop cannot form a minimum viable network (as specified by the user) within the network radius is removed. If the user specified that the minimum viable network for APop is 50, and that the network radius is 3, any patch whose APop cannot connect with 50 APop (including itself) within a 3 patch radius will be removed.

Users can set the maximum duration of a run in time steps, or specify a set of population stability conditions that will end it. Users can define a range of variation of BPop proportion over a number of time steps that will end the run. A run always ends when BPop proportion reaches 0 or exceeds 99%.

Sound options

When a run begins, there is a two second trumpet sound. A synthesizer voice note indicates the proportion of BPop in the overall population. The pitch increases as BPop proportion increases. An underlying (slightly lower volume) flute note indicates the proportion of patches in which BPop are a majority of the population, showing the spatial extent of BPop areas. The pitch increases as the spatial extent of BPop increases. A one second tubular bell note sounds when a period of population stability ends. The longer the period of stability relative to the user-defined maximum, the higher the pitch. The soundscape has a bit of a delay as the simulation progresses, and unfortunately does not work very well when doing long series of runs under behaviourspace.

Display options

The default display shows patches with no population as black. A red component shows APop and a blue component shows BPop. Patches in which both APop and BPop are present will show a proportional shade of purple.

An alternative display shows the carrying capacity a patch as a green component. In this mode, APop will show as amber to gold, and BPop will show as light turquoise to deep blue.


Starting with 1% of patches hosting an APop and a single patch hosting a BPop (of 10k patches in a 100 x 100 grid), I explored conditions under which BPop could replace APop. In all experiments APop and BPop had the same growth rate of 1.02 and the maximum carrying capacity was 200 per patch. I defined population stability as less than 5% variation in BPop proportion over 1000 time steps. All initial environments had the same structure, with 10% probability of a patch deviating from the maximum carrying capacity by up to 25%.

In the first set of experiments, I varied carrying capacity change frequency and carrying capacity change amplitude between 10% and 100%. APop diffusion rate was set to 0. In subsets of experiments, BPop diffusion rate was set at 0.1, 0.3, or 0.5. BPop network radius was set at 1 or 4. There was no minimum viable network size. All runs had a maximum duration of 50 000 time steps, but ended if the BPop proportion of the overall population reached 0 or exceeded 99%.

In the second set of experiments, I explored the influence of network parameters under good replacement conditions (see below). I set carrying capacity change frequency (CCCFr) to 20% and carrying capacity change amplitude (CCCAm) to 70%, APop diffusion rate to 0 and BPop diffusion rate to 0.01. I set carrying capacity multipliers to 1 for both APop and BPop.

I set network radius to 1 for APop. It varied between 1 and 11 for BPop. Minimum viable network size varied between 0 and 150 for both APop and BPop.

For all experiments, I did 5 runs of each combination of parameters for a total of 3000 runs in the more exploratory first experiment, 210 runs in the second experiment, and 250 runs in the third experiment.


Figure 1: Average BPop proportion at end of run, carrying capacity change frequency vs amplitude

Combining all runs (3000) in experiment one (Figure 1), the highest BPop proportion values (0.62) are at 30% CCCFr and 100% CCCAm, and 40% CCCFr and 90% CCCAm, with an arc of high values (0.55-0.62) running from low frequency-high amplitude, to low amplitude-high frequency.

Figure 2: Maximum BPop proportion at end of run, carrying capacity change frequency vs amplitude

The highest values for individual runs (Figure 2) are produced by low frequency-high amplitude change. A run at 70% CCCFr and 20% CCCAm produced a BPop proportion of 0.98. Replacement level values (<95%) are found in both low CCCFr-high CCCAm runs and in high CCCFr-high CCCAm runs, but are more abundant when change frequency is low and change amplitude is high.

Figure 3: Average BPop proportion at end of run, BPop diffusion rate = 0.01, BPop network radius = 4, carrying capacity change frequency vs amplitude

The average values are pulled down by runs in which BPop diffusion was greater than BPop growth rate. Higher BPop network radius increased BPop proportion. Focusing on runs with the best replacement conditions (Figure 3), in which BPop diffusion rate was 0.01 and BPop network radius was 4 (500 runs), shows average BPop proportion values at replacement level between 10% and 70% CCCFr and between 40% and 80% CCCAm, with the highest values (0.97) at 20-30% change frequency and 60-80% change amplitude.

Figure 4: Average BPop proportion at end of run, CCCAm = 70, CCCFr = 20, BPop diffusion rate = 0.01, BPop network radius vs minimum viable network size

The second experiment, run under good replacement conditions (20% change frequency, 70% change amplitude, BPop diffusion rate 0.01) shows that any BPop network radius larger than 1 gives BPop a substantial advantage (Figure 4) and leads to replacement under a wide range of network radius and minimum viable network size values.


In general, population replacement happens more often when the environment changes abruptly but infrequently. This agrees with Timmerman’s basic result. However, replacement can also happen when the environment changes more subtly and more frequently. One of the many problems with the archaeological record is that it was produced by a single experiment, so contingency may have played a significant role. But we can say that low frequency-high amplitude change is more conducive to replacement.

Diffusion in addition to overflow helps with replacement as long as the diffusion rate is lower than the growth rate of the population, at least using the combinations of parameters tested here. A more mobile population that doesn’t wait for population pressure to force it into new areas has an advantage over one that stays in place as much as possible.

As Timmerman suggests, the mechanism at work here seems to be that when a patch becomes depopulated as a result of a large sudden decrease in carrying capacity and later sees a sudden large increase, the more mobile and exploratory population is in a better position to move back in right away. This mechanism also operates when change is frequent and small rather infrequent and large, but it seems to work less well.

A spatially wider social network only adds to the advantage of mobility under a wide range of environmental change regimes, especially when minimum viable networks are relatively large. It allows the population to mobilize regional-level carrying capacity more efficiently and in a more timely manner.

The carrying capacity multiplier is not necessary to produce replacement in the conditions tested here. I am currently exploring its impact in a third set of experiments. The other thing I haven’t explored at all here is the impact of the initial structure of the environment. More on that soon. As always with simulation, there may be various bugs or flaws, and you may disagree with some of the design choices. I encourage all feedback and discussion, and I look forward to hearing from you.



Haws et al. 2020. The early Aurignacian dispersal of modern humans into westernmost Eurasia, PNAS 117.

Raia et al. 2020. Past Extinctions of Homo Species Coincided with Increased Vulnerability to Climatic Change, One Earth 3.

Timmerman A 2020. Quantifying the potential causes of Neanderthal extinction: Abrupt climate change versus competition and interbreeding, Quaternary Science Reviews 238.

2 thoughts on “Thinking about Neanderthals: A simple simulation tool to explore population competition and replacement

  1. “One of the many problems with the archaeological record is that it was produced by a single experiment, so contingency may have played a significant role. ” Indeed…


Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s