September 22, 2017, Friday
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at Boulder CU 
Search Links

MBW:A Multispecies Overkill Simulation of the End-Pleistocene Megafaunal Mass Extinction

From MathBio

Jump to: navigation, search

Part 1


In John Alroy’s 2001 paper: A Multispecies Overkill Simulation of the End-Pleistocene Megafaunal Mass Extinction [1], he uses mechanistic models to argue that human hunting can accurately simulate the End-Pleistocene extinction of large mammals.


It should be noted that archaeologists and paleontologists have great disdain for the overkill hypothesis on account of an extreme paucity of evidence that humans hunted North American megafauna. There are known kill sites where humans certainly hunted Mastodon and Mammoth, but these are just two of many species to go extinct during the period. Due to high rates of kill sites found in Europe, the likelihood of large scale human hunting of extinct North American megafauna is vanishingly slim. [2]

This emboldened me to try and review the paper in a critical light. The paper presents no underlying mathematical analysis and provides insufficient implementation details to attempt to recreate the experiment. Given that the author is trying to prove an unlikely hypothesis, the experiment ought to be clearly reproducible.

This made reviewing and explaining the model highly difficult as it's not clear what is actually being modeled.

This page will attempt to explain some of the underlying mathematical behavior that might cause a species to go extinct in the North American overkill model, and to point out many of the model's weaknesses.


"Megafauna" is a term used to denote a large animal. Before the Pleistocene to Holocene transition, North America was populated by 41 species of herbivorous large mammals. The body masses of mammals used in this study range from the massive 5826.6 kg Mammuthus columbi to the minute (and extinct), 21 kg Capromeryx minor.

Columbian mammoth.JPG

(Columbian Mammoth skeleton, image from Wikipedia)

30 of these 41 species go extinct very rapidly (on an evolutionary timescale) around 12,000 years before present.

Near the time of these extinctions is the arrival of the first humans across the Bering Straight. The earliest dated evidence of humans in North America is dated to 13,400 years ago. This creates a temporal overlap which suggests the plausibility of humans having a hand in their extinction.

History of Model

The father of the North American overkill model, Paul S. Martin, proposed a situation where a wave of human hunters sweep across North America after crossing the Bering Land Bridge and drove most large mammal species to extinction. Alroy’s goal is to explore this event and show that population dynamics of interacting species with realistic parameters can effectively explain the extinctions, without having to rely on climate effects.

The debate here is whether climate change OR human hunting caused these megafauna to go extinct. However it is highly unnecessary to create a model where human hunting alone causes extinction, as this model attempts to do. Rapid warming and cooling periods put a large amount of stress on animal populations. Many megafaunal species around the world come under pressure during the period and some go extinct in places humans haven't even reached. [3] (p. 5). Thus a model which assumes healthy initial megafauna populations (at their time averaged carrying capacity) is highly unrealistic. However in order to model the extinctions, "without invoking climate change," this is what Alroy does.

Simple Model

An example of a simple version of this model is to combine North American mammals into a single generic prey species and ignore any spatial effects. This isn’t very useful in this circumstance as the feedback from a dwindling prey population on the human population make it unlikely that the prey population can be driven to extinction.

Under the homogenized prey assumption, and the assumption that humans are primarily feeding on this prey, we essentially get a Lotka-Volterra predator-prey model.

LVfeedback s.jpg

In addition to this feedback making extinction unlikely, Alroy notes that several species of large mammals don't go extinct, so we shouldn't expect the homogenized model to result in extinction.

Full Model

Alroy's model is a synthesis of:

A logistic growth equation to model prey population growth, a hunting model with type II functional response, and a human population growth modeled as a function of nutritional intake.

In order to make this model more realistic Alroy simultaneously models the population of 41 large herbivore species, the human population, and human hunting. He independently models the populations in 754 geographic cells, each spanning 1 degree of longitude and latitude.

Initial Population Setup

The cells are populated by 41 different megafaunal species, of which 30 are now extinct. Each species i has an estimated mass m_{i} (kg). The range of each species is estimated based on fossil finds. A minimum convex polygon is drawn around fossil finds and cells touching the polygon are included in that population’s range.

Initial population is an estimate extrapolated from the number fossil collections which contain that species. This initial population number is assumed to be the population’s carrying capacity k

An initial population of 100 humans are placed at 49 North 114 West to represent the humans crossing the Bering Strait.

The model starts with t=0 representing 14000 years before present to allow the humans to spread out and grow enough to be present in large enough numbers for the Clovis finds dated to 13600 B.P.

Prey Population Model

To model prey population i in a cell we use Ricker Logistic Growth (scramble competition limit of Hassel equation):


r is the intrinsic rate of growth, determined by


m_{i} and k are defined above as mass and carrying capacity of the i th species.

An example of unperturbed logistic growth for initial population 100 and carrying capacity 1000 is shown below.


The relationship between mass and intrinsic growth rate r is shown below.

IntrinsicGrowth s.jpg

This spells trouble for large mammals as their intrinsic growth rates are dwarfed by those of small mammals. The extinct Colombian Mammoth, with a mass of 5826.6 kg gets an estimated reproductive ratio of 0.0141. The still extant 30 kg collared peccary meanwhile gets a ratio of 0.0988, just over seven times that of the mammoth.

I'll describe hunting more later, but this model assumes that humans hunt indiscriminately. It also assumes humans hunt each prey species in proportion to their population. Thus we can slightly non-dimensionalize this problem by only considering the percent each prey species takes up of its carrying capacity. Let n=N/ k be this non-dimensionalized proportion.

If we simplify the equations by assuming that the human population is independent of the prey populations (which is somewhat reasonable for low levels of hunting and humans getting the majority of their nutrition from gathering and hunting small game), then we can assume an indiscriminately hunting human population kills some small H proportion of each species per year.

We can alter the logistic growth equation to then include a "killed by hunting" proportion:


SimpleHunt1.JPG naturally suggests the differential equation:


The 1-n term is at most 1 near low populations. This makes the maximum reproductive ratio r-H, which clearly must be greater than 1 for the species to survive. An example with H=0.02 is given below.


Because we model r as a function of mass and H is a result of indiscriminate hunting, we can easily model all animals larger than a certain mass going extinct based on these assumptions.

Alroy acknowledges that "predicting survival strictly on the basis of body mass by declaring all species heavier than 180 kg to be extinct would identify 23 of 30 extinct and 7 of 11 surviving species, for a prediction success rate of 73%"

While his model: "correctly predicts the fates of 32 out of 41 species (78%). The exceptions are six “surviving” species that actually are extinct and three “extinct” species that actually survive."

However he doesn't acknowledge the the genesis of this effect being a simple consequence of his assumptions.

Prey Competition

We can further increase the accuracy of the model by including competition for resources between species. Competition between species is modeled by having the Ricker logistic function for each species in a cell use its own density plus the population of all other species times a competition coefficient c. Thus when c=0 the species will not compete. When c=1 all species in a cell behave as an “undifferentiated ecological unit.”

In this context we reinterpret the carrying capacity of each species i as the energy (or grazing) needs of that species. Population sizes are weighted by Energy.JPG to get them to reflect energy intake, based on Kleiber's law.

Alroy doesn't further specify the implementation of this effect, but from chapter 2.5 of Britton's Essential Mathematical Biology and some guesswork and c=1 (the preferred parameter value) we get the universal energy capacity for the cell


So in each cell, the new k value is the weighted energy needs of all S species in the cell.

The competition term in the Ricker logistic growth equation changes to make the equation for the growth of the i'th species:


I set up a toy example problem where a 5000 kg mammal (mammoth-ish mass) and a 500 kg mammal (elkish mass) share a biome and exist in the same niche (c=1). Giving them each 100,000 "units" of energy corresponds to 946 elk and 168 mammoth.

Because of the symmetry involved (The species have the same c value in each of their equations), the populations will stay in their initial balance in the absence of an external perturbation:


However introducing just a tiny human kill rate of H=0.002 slowly drives the larger species to extinction.

PerturbedCompetition s.jpg

Both species die off from initial hunting, but the more fecund (smaller) species can expand and deal with losses, expanding into the energy left vacant by the larger mammal before the mammoth can do so.

Higher levels of hunting can create a distinct competitive advantage for mammals that are less drastically different in size. With a 400 kg mammal, a 600 kg mammal, and 2% hunting, the 600 kg mammal is driven to extinction while the smaller one survives.

PerturbedCompetition2 s.jpg

Note that above I showed even a 1000 kg animal can survive 2% hunting without competition with other herbivores. c=1 competition, along with slight human pressure is highly prone to getting extinction as a model result. The indiscriminate human hunting assumption essentially demands that all but the smallest species in an ecosystem go extinct. However if these assumptions build a mechanistic model of hunter gatherer influence on populations, we would expect to be able to extrapolate these results to other hunter gatherer populations outside of North America. However large mammals weren't all wiped out all over the world. The issue here is that Alroy sells these assumptions as being the most conservative because they make no assumptions about human hunting preferences. However the assumption of indiscriminate killing is in and of itself a massive assumption, and one guaranteed to make the model predict extinctions for large mammals. This is true regardless of the other complicating factors which might seem to lend the model more credibility.

It is also interesting to note that after determining the ranges of each species, any cell containing only one species was omitted from the model: "Grid cells were excluded if their midpoints fell in an ocean or if they contained only one species, indicating an extreme lack of data." While it is reasonable to not want to include cells with unreasonable initial conditions in the study, this is somewhat equivalent to throwing out data points which don't agree with the hypothesis. The mechanism of extinction modeled here largely relies on this competition effect. The way Alroy presents his data, I'm not entirely sure that the model doesn't often simply leave only the smallest mammal in each cell and all others going extinct (the results for individual cells or individual species is never given). If this were the case, the omission of cells is a highly artificial way of getting the model to predict the desired extinctions.

Alroy does run the simulation for lower c values after fixing the preferred hunting ability parameter. For c=0 only 11 species in total go extinct while for c=1, 27 do. This demonstrates how vital competition is to the model. Alroy understates this by saying it plays a "minor role" and then leaves the parameter value at c=1. He justifies this via: "Higher values may be more realistic, because large-sized terrestrial herbivores are known to compete for food resources even with rodents. Therefore, the other simulations assume full competition." He should note that this choice is coincidentally by far the best for producing the extinctions from hunting that his model needs. It is also highly unreasonable to put all species at the same levels of competition. Almost all species in the model have well known relatives in other parts of the world or just smaller versions of themselves. Empirical data, comparisons of their diets, or comparisons of fossilized dung samples could be used to produce much more accurate competition parameters.

With species all occupying the same ecological niche, the principle of competitive exclusion demands this to be resolved by extinction, which makes c=1 an absurd assumption in a model trying to determine if extinction is a plausible outcome.

Human Population Model

Alroy models human population growth with a nutrition based exponential growth model. In each cell, calories are brought in from hunting and gathering. The calories provided determine the growth of the population in the cell.

Baseline nutrition, B=2200 kcals/person/day, is the amount of food needed to keep the population stable (r=1)

Let actual nutrition provided be A. A is capped at 1.4 times B to represent the maximum (~3150 kcals) a person can consume in a day. At most, the population can grow at a rate of 3% per year (r=1.03), which is a reasonable estimate given observed hunter gatherer population growth. The population function is designed to match this, so humans are given the variable growth ratio HumanGrowth.JPG

It isn't specified if this results in a final discrete population equation of the form of either:


or something like:


but the growth is slow enough that the linear and the exponential form are negligibly different.

Human Gathering Model

Alroy's comments on secondary food resources are limited to:

"Productivity of plant and small game foods was assumed to be 200 g/m^2 for all geographic regions because areas with higher net primary productivity than this probably yield relatively less edible plant resources (9). Of this amount, 2% was assumed to be edible (9). Regrowth was modeled with logistic dynamics, with the intrinsic rate of increase set to 22% the initial value (9). Caloric values were set to 3 kcal/g (13), which is similar to estimates used in other simulations (10)."

The assumptions lead us to the calculation that for each square meter there are 200*(0.02)*3=12 kcals of edible food. In a square kilometer this scales up to 12 million kcals. Alroy gives population densities in 100 km^2, so this scales to 1.2 billion kcals per 100 km^2.

To put this in perspective, the modeled population consumes at baseline 2200 kcals/day. So in, "baseline consumption person years", this is:

1.2 billion * 1/(365*2200)=1500 (person years).

Let G(t) be the person years of food available in year t.

"Logistic dynamics" suggest a regrowth formula:


Where r=1.22 and k=1500.

The differential form:


Will tell us the maximum number of supportable humans if we find its maximum.

Differentiating and then finding the root numerically suggests a maximal growth amount when "G"=729 person years.

As expected from the shape of a logistic growth function, this is near half of the carrying capacity.

Plugging 729 in to find dG gives 87 person years / year.

i.e. every 100 square kilometers can support 87 persons on secondary food resources.

This is far above the realized densities in Alroy's model, but he doesn't give any explicit equations for how secondary resources are gathered, so I don't know how to do more analysis on it. However his referenced material may make things clearer.

By looking at the table of results, we see that humans with high hunting ability have the highest max and final population densities, yet they have quickly driven all prey species to extinction and have a "meat in diet" proportion of 0.000. Thus the hunting ability must also be determining foraging ability. This is highly misleading as Alroy states:

The simulation results are unambiguous (Fig. 1 and Table 1). Human population growth and hunting almost invariably leads to a major mass extinction. In fact, it is hard to find a combination of parameter values that permits all species to survive. These few scenarios (trials 1 through 3 and 18 through 21) require very low final human population densities of, 0.13 people/100 km^2, an order of magnitude below the observed range for modern hunter-gatherers (4).

The simulations with low "hunting ability" are much more accurately described as very low foraging ability. However because the two ability are totally intertwined in this model, realistic human population levels can only be reached by increasing this one parameter. Due to this, Alroy sells trials in which humans hunt most of the megafauna to extinction as the most realistic.

To do this honestly, foraging ability and hunting intensity undoubtedly need to be controlled by separate parameters. There isn't really a logical reason to combine the two, and many reasons to keep them separate. Much of the criticism of the overkill model is its envisioning early humans as overly risky and highly focused on hunting big game. With higher foraging ability but lower hunting ability the model might be able to predict reasonable human population levels, while also predicting much less hunting impact, which would do a better job of explaining the lack of empirical evidence for hunting most megafauna.

Human Hunting Model

To model hunting we first calculate baseline hunting success h' , (grams of food) by summing across all prey species in cell:


m_{i} = body mass of i'th species,

N_{i} = population size of the i'th species,

a = hunting ability constant (low:24, preferred:38, high:54, (in millionths))

Nh = human population size.

To find realized success “h” (grams), we use a type 2 function response with the equation: h = (3.5h'/hk)/(1 + 2.5h'/hk) , which converges to 1.4 for high h' values.

hk= hunting rate at equilibrium (amount of food needed to sustain population).

This suggests h represents a nutritional proportion.

The baseline success rate is somewhat bizarre as the hunting ability parameter a is only involved in computing the interference effect when the human population is too large. "Human population size" naturally seems to mean the population of the cell. We can plug in some reasonable numbers to see how the equation behaves.

The preferred hunting ability constant is 38 x 10^-6. A reasonable population density from the model is 5 people per 100 square km. The size of a cell of 1 degree longitude and 1 degree latitude depends on your proximity to the equator. At Boulder Colorado the size is 85 km by 110 km, for 9350 square kilometers in a cell. This produces a reasonable population estimate of 5 times 93.5, i.e. 467.5.

However e^467.5 is a nonsense number for us to involve in the calculation. Furthermore, the hunting ability parameter a is clearly negatively correlated with hunting success, which is nonsensical given what the parameter represents.

Ecological Response

The modulation of hunting via a type II ecological response is actually a curious choice and deserves some examination. When I first started this project, the initial idea that excited me was that it might be possible to abstractify human ability to drive animals to extinction based on human versatility.

Vertebrate predators are generally specialized and relatively focused on just one or a couple prey species. Due to this, it might be much harder for these predators to drive their prey species to extinction. We would see a feedback effect as prey species dwindle and predators die off from lack of food resources.

On the other hand, omnivorous and fast learning humans hunter gatherers might be much more able to switch between prey species or foraging resources as their preferred food sources dwindle.

One way of exploring these relationships are the Holling functional response equations.

Let H be the hunting capture rate and N be prey density.

A type I functional response describes a linear relationship between between prey density and hunting rate:


alpha is just a positive constant of proportionality. This simple relationship is seen in the LV equations, but in practice it is most applicable to filter feeders like sponges or baleen whales as these animals have no real prey handling time [4]

The type II functional response on the other hand describes a saturating relationship. At high levels of prey density more of the predator's time is taken up by prey handling and feeding while less is spent searching for prey.


Taking the limit as prey density approaches infinity reveals that 1/t represents the hunting rate at saturation.

Dividing the numerator by the denominator changes the form of the equation to:


Which, for |alpha*t*N|<1, can be simplified by representing the right part of the equation as a geometric Taylor series:


Which, for low prey density simplifies to:


This is nice as we see the expected return to a linear relationship as prey density approaches zero. This is expected as the predator will then be spending all of its time searching for prey and virtually no time handling kills.

This functional response is useful for characterizing a human saturation with food, however it is lacking when it comes to human relationships with individual food sources. This is on account of the fact that vertebrates are typically characterized by type III functional responses.

The type III response:


1/t still gives us saturation, but now values of the power n>1 now dramatically change behavior for low prey concentrations.

Applying the same logic as for the type II response, the behavior for low prey densities is:


I.e. for low prey densities the capture rate grows as the function alpha*N^n. n>1 Implies that for very low prey density the capture rate is minuscule, however it grows at a faster than linear rate.

The three functional responses, with n=2 can be compared graphically:


Zooming in on low prey densities:

Lowfunctional s.jpg

This is somewhat confusing, but it makes sense if thought of as a learning requirement for complex predators. At very low levels of prey density they switch almost entirely to an alternative source of food. As the density increases they get better at hunting it. However eventually saturation occurs, just as in type II, and all time is spent processing kills.

This relationship is highly interesting for extinctions as a type III functional response creates a sort of natural defense mechanism from overkill. If there are other abundant species of prey the predator might hold a rare species at low densities in a "predator pit" [5].

To model the functional response of a generalist predator the appropriate functional response form is:


S is the number of prey species modeled

i and j are indexes representing prey species and the parameters are all prey specific.

The alpha values represent a preference value for each species.

The n values are critical for determining the behavior of the system.

Estimating all of these parameter values for modern ecological systems is difficult. Accurately doing so for extinct ecological systems may be impossible. However this does not justify using a universal n value of one. While hunting preferences may always be guesswork, the behavior of complex generalist predators demands that n>1.

By using a universal type II response Alroy makes this assumption. Despite this assumption being massively extinction predicting, Alroy doesn't justify it. He also doesn't mention it as a weakness of his model, stating that,

"Many complicating factors are ignored by the model. Several involve aspects of human biology and history: selective human hunting of individual prey species; human-induced habitat change; ... By and large, adding these factors to a model should increase extinction rates"

While it is certainly true that assuming humans preferred to hunt all extinct species could produce a model that predicts extinction of extinct species, it is absolutely unreasonable to state that adding most complicating factors is in the hypothesis's favor.

Spatial effects

Human movement is modeled as a diffusion from cells of high concentration to cells of low concentration.

Animal movement is represented by (for most trials) having the animal population redistributed across its range each year. Some trials are run without this effect and each cell involves no movement of prey species.

An important effect of the 754 cells is that it allows the model to have a variety of overlapping ranges of prey species. Due to competition, other mammals have the opportunity to take up some of the carrying capacity freed by the death of some larger species. The presence of other mammals can also have an effect by providing the humans with enough food to sustain a higher population. This larger human population will then have a higher nutritional need and hunt more.

However it creates the substantial drawback of making the experiment virtually impossible to recreate. With 754 cells there is potentially now 754 initial condition setups for the experiment. Each cell can have its own range of diversity of species and carrying capacity.

Doing an experiment in a single cell and presenting population graphs for this cell makes what is going on tremendously clearer.


Here is an excerpt from Alroy's table of results:

Table s.jpg

The hunting ability constant a is very closely related to human population density, and the extinction of species. This is the independent variable here.

The densities are given in people per 100 km.

Meat in diet represents the proportion of food the diet that comes from the hunting modeled.

Energy use represents the amount of their carrying capacity the prey use at the end of the simulation.

Actually extinct refers to the number of the 30 species that actually are extinct that go extinct in the trial.

Actually surviving refers to the number of the 11 species that are not extinct that go extinct in the trial.

True Fates is the number of species that the trial "gets right" and correctly predicts extinction or survival.

The preferred trial for the experiment is trial 8, on account of it having the highest number of True Fates, and reasonable numbers for human population density and the time at which species go extinct. This trial uses the competition parameter c=1, dispersal of prey species about their range each year, and an initial human population of 100 starting from Beringia.

Alroy also gives the a logistic overall North American population map for this trial:

Timelines s.png

The thick black line represents human population. The thin black and thin grey lines are the represent surviving and extinct species respectively.

The fact that the trial is preferred primarily because it backs up the hypothesis that humans caused the extinctions is highly troubling. It is understandable to vary model parameters until they agree most closely with empirical data, however which species go extinct is only useful empirical data here if hunting is certainly the mechanism which causes them to to extinct.

Trying to make the simulation predict reasonable hunter gatherer population densities is a much better approach. This trial apparently does line up relatively well with expected population densities, but as I mentioned in the gathering model section, hunting ability and gathering ability are conflated in Alroy's model. As you can see from the bizarre relationship between hunting ability, meat in diet, and the population densities, the better hunters quickly wipe out all the animals yet don't have a population crash and their diet approaches zero percent meat. There are no trials in which humans are efficient gatherers but super efficient indiscriminate hunters.

In addition to this, to describe the trial in which humans are predicted to drive all of these mammals to extinction borders on an absurd rejection of empirical evidence. In, A requiem for North American Overkill [6], Grayson and Meltzer write:

How many of those genera can be shown to have been human prey during Clovis times? The answer is two–mammoth and mastodon- and there are only 14 sites that securely document this relationship. As has long been known, this is not a sampling fluke. There are more late Pleistocene occurrences of horse than there are of mammoth or mastodon, and nearly as many for camel as for mastodon, yet there are no demonstrable kill sites for horse or camel or for any of the remaining genera. This is not for want of looking. Given the high archaeological visibility of the remains of extinct Pleistocene mammals, and their great interest to archaeologists and Quaternary paleontologists alike, if such sites were out there, they would surely be found. Indeed, there is a strong bias in the Clovis archaeological record toward just such sites.

The rarity of megafaunal kill sites is such an evident feature of the late Pleistocene archaeological and paleontological records of North America that Martin has had to address it. After all, other parts of the world—Late Pleistocene Europe, for example—are littered with sites that document human predation on large mammals.

From this context, varying the "hunting" parameter until it predicts all of these extinctions from hunting and labeling this the most likely value for the parameter is highly illegitimate. Choosing the hunting ability parameter to best fit extinction events only makes sense if the extinction events were caused by hunting. Because the hypothesis of the paper is effectively, "The megafaunal mass extinction was caused by hunting," choosing this parameter to agree with the hypothesis is highly circular logic.

Even when varying the parameters for the most true fates, the best the model ever does is predict 32/41 correct results. This is a much better than random result, but simply killing off every species in the model gives an accuracy rating of 30/41. Alroy also never gives results for individual species, other than to say that it has a difficult time prediction some many species to go extinct without always driving Moose, Elk, and Bison to extinction as well. This isn't surprising, as these animals all have very large body masses and the simulation is inherently designed to predict extinction of large mammals.


Grayson and Meltzer also specifically address Alroy's model, (referencing specifically the paper I'm reviewing) but their concerns are over the empirical lack of evidence for hunting. My concern is more that Alroy claims to have tried to make a highly conservative model that cant help but predict hunting causing a mass extinction.

Alroy presents this as a clear and inevitable conclusion:

"The improved, ecologically realistic model outlined here challenges the common-sense notion that no amount of overkill could have resulted in a true megafaunal mass extinction. The simulations demonstrate not merely that overkill scenarios are plausible, but that an anthropogenic extinction was unavoidable given the facts of ecology and the fossil record—even assuming that human predation was limited and nonselective."

If Alroy presented his model as a way of arguing that humans could cause extinctions under just the right circumstances it wouldn't deserve so much scrutiny and scorn. Instead he has subtly rigged up his model not to be accurate, but to be an asset in this absurd argument over what killed the megafauna.

While the model does a slightly better job of describing extinctions than simply choosing to kill off all mammals heavier than some mass (30 true fates vs 32), Alroy never actually demonstrates where he gets these extra two correct results, or even which species they are. Due to how enormously favorable the model is to the survival of very small mammals, my prediction is that some smaller species which actually go extinct share a range with the minuscule Capromeryx minor or Collared Peccary, which effectively inherit this models world. Alternatively, some larger species that actually survive may be blessed with a range free of dangerously small mammals.

While some models aren't as realistic as possible as a trade-off for mathematical elegance and ability to probe them analytically, this model is neither. It's both inaccurate and highly opaque. Given its disdain for both reproducibility and empirical evidence, it doesn't belong in the realm of science.


Species Info

Body masses (kg) used in study:


30 Pecari tajacu

60.6 Rangifer tarandus

68 Antilocapra americana

91 Oreamnos americanus

91 Ovis canadensis

106.85 Odocoileus virginianus

118.1 Odocoileus hemionus

286 Ovibos moschatus

422 Bison bison

457 Alces alces

500 Cervus elaphus (now Elk)


21 Capromeryx minor

45 Oreamnos harringtoni

52.5 Platygonus compressus

52.5 Stockoceros conklingi

53.7 S. onusrosagris

60.6 Tetrameryx shuleri

73.5 Mylohyus fossilis

222.6 Navahoceros fricki

238.1 Hemiauchenia macrocephala

244.7 Palaeolama mirifica

306 Equus conversidens

312 Holmesina septentrionalis

324 Tapirusveroensis

368 Equus francisi

439 Equus complicatus

485.6 Cervalces scotti

498.8 Euceratherium collinum

522.8 Bison priscus

533.4 Equus niobrarensis

555 Equus scotti

574 Equus occidentalis

614 Nothrotheriops shastensis

665.8 Glyptotherium floridanum

753 Bootherium bombifrons

995 Camelops hesternus

1320 Megalonyx jeffersonii

1990 Paramylodon harlani

3174.2 Mammuthus primigenius

3297.7 Mammut americanum

5826.6 Mammuthus columbi

There is a clear gap in the species used in the simulation. There are 14 species under 120 kg. Half of these go extinct and half survive. The next smallest mammal is 222.6 kg. Of these mammals larger than 220 kg, 4 survive and 23 go extinct.

Part 2


In this part my aim is to create a create a multispecies overkill model with a careful examination of underlying assumptions and their influence on the simulation. In some aspects this will be a substantial simplification of the Alroy model. In order to effectively communicate the effects of assumptions it is much easier to work in a single, well mixed compartment.

It also may prove useful to model human population as independent from megafaunal populations. For a megafauna population nearing extinction, the species must provide a negligible amount of nutrition to the human diet, and thus have a negligible effect on the human population. One possible model for human population dynamics are limit cycles \cite{Belovsky}, which we will examine in order to shed light on the impact of human population boom/bust on mammalian populations.

In some respects, my model will need to follow in the the footsteps of the Alroy model simply due to the constraints of knowledge and expertise. While it is easy to find examples of fault and scientific inaccuracy in a model, it is impossible to create a model combining all of human scientific knowledge and expertise. Thus some assumptions I found unacceptable in my previous paper have to be used so the model/paper can be produced at all.

Models Drawn From

Alroy Model

Many of the ideas for this model are built upon those in the Alroy model. It is very difficult to obtain virtually any desired parameter for an ecology that no longer exists. While it is possible to make highly specific estimates for each species, these are still just estimates and would require a model to veer wildly off course. This leads me to adopt many of the strategies used in the Alroy paper in order to be able to focus on the relationships I intend to study.

For intrinsic fecundity I will use the same formula. Which, for working in kilograms is:


Another idea borrowed from Alroy is the use of a universal competition parameter. Again,while possible to do better, i.e. specifying the parameter for each animal pair, this growsunwieldy incredibly quickly. The number of relationships is the factorial of the number ofspecies involved.

Multispecies Ecological Response Models

One parameter I find far too simplified in the Alroy model is the ecological response power and the grouping of all hunted animals into the same ecological response. The generic functional response equation is given by:


Alroy's model assumes this holds universally across hunting and gathering food sources and with n=1.

This is dangerous because food gathering is a highly different process from hunting, and because humans are generalist predators.

For generalist predators prey switching can occur whereby at low densities of preferred prey the predator can switch to an alternative source of food. To model this we use the multispecies functional response:


Which tracks consumption of each prey item.

In "The Functional Response of a Generalist Predator" Smout et al. study the relationship between a predator, Hen harriers, and three primary sources of food during breeding season: meadow pipits, voles and red grouse chicks. Using data on both predation rates and prey populations, she calculates n values of 2.51, 1.14, and 1.18 for grouse, vole, and pipit hunting respectively.

Higher n values correspond to more prey switching when that prey item becomes scarce.

In, "The roles and impacts of human hunter-gatherers in North Pacific marine food webs," Dunne finds evidence that while humans can certainly have a devastating effect on ecosystems, a prey switching generalist human predator promotes ecological integrity through the prey switching effect.

Hunter-Gatherer Population Model

Belovsky's, "An optimal foraging-based model of hunter-gatherer population dynamics" mechanistically models human population in relation to the primary productivity of a landscape. Belovsky's model is widely cited and provides effective predictions of hunter gatherer population densities. It is also an ubiquitous source for the population modeling in Alroy's paper.

Despite the paper's success, it presents the somewhat bizarre conclusion of predicting violent limit cycles in population for large primary productivity. E.g. (figure 7 in Belovsky's paper).


However these cycles provide intriguing circumstances under which animal populations may become substantially strained as humans attempt to ward o�ff starvation during a population crash. His paper is also useful in providing estimation tools for primary productivity based on the biomass of large herbivores.



Megafaunal Population

Population estimates for large mammals in the late Pleistocene are very difficult and different methods obtain vastly different results. Determining accurate populations and the probabilities associated with their accuracy is a mathematical biology project in and of itself. For populations within a specific time frame (Late Pleistocene/Early Holocene transition), the possible sample size of data is minuscule.

Using estimates based on numbers of fossil finds and estimated "live:dead" ratios, Hanes gets rough sizes of the four most prevalent species of North American extinct species of megafauna:

Mammuthus, N~174,100

Mammut, N~72,500

Equus, N~60,000

Camelops, N~59,000

Essentially we assume that for each find there where 1000 individuals alive at any one time.

Using these numbers would perhaps make sense if we assumed all of North America was a well mixed container, but this is highly unreasonable. The origin of these populations being relatively small is that it is inaccurate to think of North America as a continent wide mammoth steppe teeming with life. Instead, what is borne out from the fossil record is a collection of relatively small "sources" surrounded by large "sinks" (Hanes). The megafaunal populations within the source can sustain themselves and grow while in the sinks they mostly decline. Thus there is a significant amount of clustering of the populations around the sources.

This is explained in part by the rapid climate changes. While almost all of the continental US was no longer under a glacier, plant life lags behind climate changes significantly due to the difficulty of ecological succession.

Regardless, this suggests not modeling the populations continent wide, but instead within a source.

This is an important change from the Alroy model because it will likely greatly increase the ratio of megafauna to humans. By spreading the megafauna populations across a convex polygon of their entire range, it pulls population from the sources into the vast sinks. Most of human nutrition probably came from gathering, but it is likely that regions unsuitable to other mammals were also unsuitable to humans.

For this project I used FAUNMAP to view the locations of fossil sites and chose to examine a cluster of fossil sites making up most of southern Missouri.

Cluster s.png

I draw a convex polygon around cluster of late Pleistocene/early Holocene megafauna finds. Area: 127000 km^2

After searching through the sites for megafaunal herbivores we get the minimum number of individuals required to create the fossils:

Deer: 23

Bison: 6

Elk: 5

Mastodon: 4

Goat: 2

Muskox: 2

Mammoth: 1

Horse: 1

Giant sloth: 1

Moose: 1

While this is a too measly a sample to put any real faith into, it provides a vague indication of the megafauna in the region; mostly deer, elk, bison, and mastodon. While there is only one Mammoth find, this is within the range of Mammuthus columbi, and at 5800 kg per animal they can still take up a huge chunk of total biomass.

It is worth noting that including carnivores in this population estimation technique yields an impossibly high proportion of carnivores due to many of the sites being caves. This calls into question whether these are relative numbers of herbivores, or relative numbers of herbivores hunted by cave dwelling predators.

Studying just Mammoth, mastodon, bison, elk, and deer provides a wide range of population sizes and body masses, while preventing excess clutter.

The estimated x1000 relationship between finds and living specimen gives us an initial population state vector:

Initial populations: N_0=1000x[23, 5, 6, 4, 1]

for parameter vectors:

Species: [Deer, Elk, Bison, Mastodon, Mammoth]

Mass in kg: M=[110, 250, 470, 3300, 5800]

Unfortunately however, these numbers result in a prey biomass for the region that is unreasonably low, 0.2 g/m$^2$. While we didn't include many species present in the area, we kept the primary contributors to biomass, and want these species to represent a reasonable ecology. For the generally expected harvestable primary productivity average of 200 g/m$^2$ used by Alroy and Belovsky, we would expect the biomass to be approximately 1 g/m^2. As this is a region of high megafaunal density, we would expect to at least meet this average benchmark. Thus we scale up the estimated live:dead ratio for prey species by 5 to 5000.

Initial populations: N_0=5000 x [23, 5, 6, 4, 1]

As these are grazers, simplifying the inter-species competition scheme to full competition is not nearly as unreasonable as it was with mammoth fully competing with collared peccary. We approximate the ecological carrying capacity with the populations and the M^0.75 relationship to energy via Kleiber's law.


For the currently occupied carrying capacity, take


Applying the continuous logistic growth equation (as Alroy does, except continuous instead of discrete) with a system wide competition term, this gives a growth equation:


As a simple check of the model thus far, we model the effect of a sudden tripling of the available carrying capacity:


Biomass scaled logarithmic growth. While mastodon and mammoth dominate initial population, their lower intrinsic growth rates result in substantially smaller gains than those of deer.

Human Population

Within the convex polygon there are also a minimum 3 human finds across 2 sites dating to the the time period in question. This is quite substantial as there are only 18 total Homo sapiens sites within the time period on FAUNMAP.

Into the subsequent Holocene epoch, the Missouri region becomes dominant in terms of human population.

Holo s.png

In keeping with the ratios used before, three Homo finds suggests a population of 3x5000=15000 humans in the region. This produces a population density of 12 humans /100 km^2. This is surprisingly well matched to the population seen in the Belovsky model for PROD=200 g/m^2.

In order to achieve similar effects to a human population near its carrying capacity, we hold the human population constant at its expected rate.

I originally wanted to be able to simulate the human population as well, however this simply adds many more assumptions that have to be made in order to produce reasonable effects on the rest of the system.

Ecological Response

In the general multispecies functional response equation we use biomass instead of true population numbers as human consumption is limited by the amount of biomass they can consume.

We use the nutritional estimates of 3 kcal/gram biomass and limiting human consumption of 3150 kcal/day, which translates to 3150x365x(1/3)(1/1000)=380 kg/year.

To model prey switching (Dunne) I used an equation of the form:


The value returned, F_i gives the biomass hunted of the ith species as a proportion of saturation.

The significant drawback of this method is it draws attention and importance to somewhat arbitrary distinctions between prey items. If plant gathering is taken together as a lumped biomass, it will dwarf all other sources of food. Using the scalar to diminish its influence requires making including it in the first place pointless.

Observed foragers (Belovsky) consume about 69% vegetables, 31% meat. So we can reduce the problem to looking at prey alternatives for the 31% of a maximized diet (3150 kcals). Making the new biomass consumption 118 kg/year.

Parameter Variation

To study the effect of the ecological response we want to create similar hunting levels while varying the parameter in question. An approximate solution for this is to also take the small attack coefficient to the (1+q) power.

We then solve the system for q=[0,1/2,1]

In addition to this, For some trials I introduce massive 100 year period limit cycles in human population.


The relationship between the variables ended up being significantly less profound than was hoped. The first trial:


This displays the similar unlikely behavior of the Alroy model, increasing the strength of the type 3 ecological response simply pulls all of the final populations in the model closer together. However increasing the response inherently caused the equations to put much more pressure on all of the populations. In order to see the prey switch behavior the pressure put on all of the populations had to be continually tuned down, which substantially reduces any claim that I came to a natural result mechanistically.

Second trial:


Third Trial:


Also largely inconsequential were the large human population limit cycles. Due to the variation in pressure from the human population, I had speculated that this could allow for advances by the rapidly reproducing populations against the slower ones, but both century long massive cycles with huge in population and improbably fast and violent cycles caused no qualitative change other than a slight ripples in the prey populations:


q=1, Type 3 Response with limit cycle


An enormous portion of what I learned in trying to put this model together is the absurdity involved in the process of speculative modeling in a historical context. In trying to develop accurate population numbers within a region of high megafauna density, I ended up getting ridiculous numbers and had to scale them massively in order to return to normal values.

However ultimately I largely lost sight of whether the overarching mathematical relationships would be either correct or insightful. What I did here is essentially what Alroy did in his paper except opposite in that I made some sweeping assumptions which inherently led to it being impossible to produce extinction outcomes.

A possible important alternative method of modeling megafauna population is a structured population model. This is due to the massive preference predators have for hunting baby mammoth/mastodon, as they were presumably much less dangerous to hunt. Humans also have been shown to have a large preference for hunting baby elephants. At the Lehner kill site for example there are 9 mammoth fossils, and all of them are from immature animals (Haury). This suggests examining these species in a structured manner, as hunting may effect their population in a way strongly biased towards killing youngsters.

However humans also hunted adult mammoth, which another interpretation of the structured population may show to be more influential than juvenile deaths. If humans could hunt adult mammoth while other predators could not, that may have suddenly put a large amount of pressure on their population.


Matlab Code


%main %clc %clf

global M r K A q;

TIME=1000;  %years tspan=[0,TIME];


%Humans h=15000;

%Human Type 3 Strength q=1;

%Human Appetite A=118*h;

%Species S=['Deer', 'Elk', 'Bison', 'Mastodon', 'Mammoth']; %S=['Elk', 'Bison', 'Mastodon', 'Mammoth']; %S=['Bison', 'Mastodon', 'Mammoth'];

%Masses M=[110;250;470;3300;5800]; %M=[3300;5800];

%Growth Rates r=0.3467*M.^(-0.37);

%Initial Populations N0=5000*[23;5;6;4;1]; %initial state vector %N0=5000*[4;1];

%Carrying Capacity K=sum(N0.*M.^(0.75));

%SOLVER [T,N] = ode45(@f,tspan,N0);  %Refer to f.m for equations

%PLOTTING clf; hold on;

for n=1:5 plot(T,M(n)*N(:,n)) end

xlabel('Time (years)'); ylabel('Total Biomass (kg)') legend('Deer', 'Elk', 'Bison', 'Mastodon', 'Mammoth','Rel capacity filled') title('Type 3 Response') shg


function [ dN ] = f( t, Nt ) global M r K A

%Occupied carrying capacity Kt=sum(Nt.*M.^(0.75));

%Human hunting and gathering H=A*FR(Nt)./M;

%DIFFERENTIAL EQUATIONS dN=zeros(5,1); for n=1:5 dN(n)=r(n)*Nt(n)*(1-Kt/K)-H(n); end

%dN=[r.*N.*(1-Kt/K)-H]; end


function [ Rate ] = FR(Nt) global M A q %q Type 3 strength n=1+q;%exponent Bo=A/2;%Half saturation B=Nt.*M;%Biomasses





J. Alroy, A multispecies overkill simulation of the end-Pleistocene megafaunal mass extinction, Science 292 (2001) 1893–1896.


Details (some) at: [8]

A PDF of the paper can be found at [9]

G. E. Belovsky, "An Optimal Foraging Based Model of Hunter-Gatherer Population Dynamics" J. Anthropol. Archaeol. 7, 329 (1988).

Britton, N. F. Essential Mathematical Biology. London: Springer, 2003.

Dunne, J. A. et al. The roles and impacts of human hunter-gatherers in North Pacific marine food webs. Sci. Rep. 6, 21179; doi: 10.1038/srep21179 (2016).

FAUNMAP Working Group. 1994. FAUNMAP: a database documenting late Quaternary distributions of mammal species in the United States. Illinois State Museum Scientific Papers 25(1-2):1-690.

Grayson, D. K. Meltzer, D. J. "A requiem for North American overkill." Journal of Archaeological Science 30 (2003) 585–593 [10]

Hanes, Gary. "Late Pleistocene megafaunal extinctions." Royal Tyrrell Museum Speaker Series 2012. Lecture.

Haury, Emil W., E. B. Sayles, and William W. Wasley. ``The Lehner Mammoth Site, Southeastern Arizona”. American Antiquity 25.1 (1959): 2–30. Web

Smout, Sophie et al. “The Functional Response of a Generalist Predator.” Ed. Frederick R. Adler. PLoS ONE 5.5 (2010): e10761. PMC. Web. 30 Mar. 2016. [11]