Modeling post-Pleistocene megafauna extinctions as complex social-ecological systems

Abstract

The role of human hunting behavior versus climate change in the mass extinction of megafauna during the Late Quaternary is much debated. To move beyond monocausal arguments, we treat human–megafauna–environment relationships as social–ecological systems from a complex adaptive systems perspective, to create an agent-based model that tests how human hunting may interact with environmental stress and animal life history to affect the probability of extinction. Using the extinction of Syncerus antiquus in South Africa at 12–10 ka as a loose inspirational case study, we parameterized a set of experiments to identify cross-feedbacks among environmental dynamics, prey life history, and human hunting pressure that affect extinction probability in a non-linear way. An important anthropogenic boundary condition emerges when hunting strategies interrupt prey animal breeding cycles. This effect is amplified in patchy, highly seasonal environments to increase the chances of extinction. This modeling approach to human behavior and biodiversity loss helps us understand how these types of cross-feedback effects and boundary conditions emerge as system components interact and change. We argue that this approach can help translate archaeological data and insight about past extinction for use in understanding and combating the current mass extinction crisis.

Introduction

We are currently amid a highly accelerated anthropogenic extinction event conservatively estimated to be one hundred times higher than expected background rates (Ceballos et al., 2015). It is estimated that 48% of animal species are currently in decline (Finn et al., 2023) and that 50% of the plant and animal species on Earth may be lost by 2100 (Braje and Erlandson, 2013, p. 21). Biologists have identified four recent periods of significant acceleration in the extinction rate of land animals (63.8–32.2 ka, 16–9.5 ka, 2300–600 yr, and 180–120 yr) (Andermann et al., 2020). There is a general temporal correlation between the dispersal of anatomically modern humans and the demise of faunal species around the world (Braje and Erlandson, 2013; Boivin et al., 2016), but the role of human-driven effects versus climate-driven effects in these extinctions continues to be debated (Barnosky et al., 2004; Koch and Barnosky, 2006; Braje and Erlandson, 2013).

In our opinion, archaeological data about past extinction is not being utilized to its fullest potential in combating the current mass extinction. Analyses of specific extinction events are typically reductive and attempt to identify and test the validity of specific and often monocausal explanations for specific past cases of population decline, extirpation, or extinction (Supplemental Text, Section 1). The coarse-grained and incomplete nature of archaeological, paleoenvironmental, and paleoclimatological data often combines with the highly variegated nature of extinction events across time and space to provide fodder for endless debate about causal mechanisms for extinctions in various parts of the world (Lyons et al., 2004; Turney et al., 2008; van der Kaars et al., 2017; Saltré et al., 2019). If we want these findings to be relevant to current-day issues (i.e., “action archaeology” sensu Sabloff, 2016), it is growing clearer that multifactor explanations are needed. Additionally, the combination and relative influence of the potential factors driving extinction need to be contextualized at a higher level than individual cases of past extinction events. But reductionist inductive approaches can at best only show correlations in the proxy record of the past, which are sometimes weak, sometimes spurious, and often inconclusive (Faith et al., 2020). The dynamics of how extinction events played out over time are completely missing from even the best-preserved proxy records, and even the most detailed analyses of fossil data sets are flawed by the sampling biases introduced by time, taphonomy, and by how we find, retrieve, and analyze material remains of the past.

We intend to show that a more fruitful way for archaeologists to create actionable knowledge relevant to addressing contemporary issues is by using paleontological and archaeological evidence to parameterize experiments within formal, generative, deductive models derived from theories of human and animal behavior as well as theories of climate and environmental change. To do so, we employed the agent-based modeling (ABM) formalism (Epstein, 1999; Bankes, 2002) informed by a complex adaptive systems (CAS) perspective about social–ecological systems (SES) (Preiser et al., 2018) to create a deductive model of interactions among human hunting preferences, environmental change, and animal demography that is used to investigate extinction dynamics in large-bodied, k-selected, mammals (i.e., megafauna >45 kg).

ABM is a computer simulation formalism (Supplemental Text, Section 1) centered around multiple interacting software agents that can sense stimuli in their environment and respond according to behavioral rules (Epstein, 1999; Bankes, 2002; Romanowska, 2015; Wilensky and Rand, 2015). The CAS framework identifies aggregates of interrelated variables that form systems, such as the individual cells of the nervous system, or the individuals and institutions of an economic system (Daems, 2021). The components of these systems change and reorganize themselves to adapt to their surroundings (Holland, 1992). There is no central all-controlling mechanism that guides the system's behavior (Mitchell, 2009). Rather, complex, self-organizing behavior emerges out of simple operational rules and complex feedback interactions (Kohler, 2012). SES are a special type of CAS where human social system components connect to and interact with ecosystem components (Walker et al., 2006; Glaser et al., 2012; Preiser et al., 2018).

In our model, we consider humans, prey species, and the environment to be key components of an interconnected, adaptive SES, in which feedback occurs between and among individuals and communities as the system changes over time. We demonstrate the usefulness of this approach by presenting a loose case study inspired by megafauna extinction dynamics that may have played out in the greater South African Cape Floristic Region across the Pleistocene–Holocene transition. All models are abstractions of the real world and are by definition limited in scope, realism, and detail (Box and Draper, 1987), but deductive models are flexible and are preferred for the initial exploration of the capacity for emergence within a general problem domain (Robinson et al., 2007; Rounsevell et al., 2012; Railsback and Grimm, 2019). Specificity and complexity are added to the model as required by our understanding of the behavior, or as the model problem domain evolves, which is the approach that we have taken here (Railsback and Grimm, 2019; Romanowska et al., 2021). We think that these methods may allow us to take lessons from past extinction dynamics and apply them to present SES by identifying extinction thresholds and the ways in which such systems move towards or away from them.

Methods

The megafauna hunting pressure model

We created our ABM simulation model, which we call the “Megafauna Hunting Pressure Model” (MHPM), using the NetLogo modeling platform (Version 6.3.0; Wilensky, 1999). The MHPM is described in detail using the Overview, Design Concepts, and Details (ODD) protocol (Grimm et al., 2006, 2010, 2020) in the Supplementary Text, and summarized here following the Summary ODD format outlined by Grimm et al. (2020).

The MHPM is a model of large-bodied ungulate population dynamics with human predation in a simplified, but dynamic grassland environment. The overall purpose of the model is to understand how environmental dynamics and human predation preferences interact with ungulate life history characteristics to affect ungulate population dynamics over time. The model considers patterns in environmental change, human hunting behavior, prey profitability, herd demography, herd movement, and animal life history as relevant to this main purpose. The most important design concept of the model is the way that we have composed and connected the following functional areas in regard to prey animal extinction dynamics: (1) human and prey animal energetics and foraging, (2) herd animal demography and life history, (3) human and prey mobility and movement, and (4) environmental dynamics. The MHPM contains individual prey animals, forager bands, and landscape patches as entities. Each entity type has several state variables and attributes, and the model has temporal processes that involve each type of entity. Table 1 lists all of the model's input state variables and maps them to the model entities and to the basic functional areas of the model.

Table 1. List of MHPM model variables with mappings to variable type, model entity, and general functional area(s) within the model

Prey animals can be male or female, and prey animal sex is related to payoff and costs to human foragers who may prey on them. Prey animal population change is a key process in the model. Animals can mate and produce offspring, and female prey animals carry offspring through a gestational and weaning period, after which offspring become independent animal entities. The foraging payoff of in utero or unweaned offspring is added to that of their mothers during this period. Prey animals must graze on vegetated landscape patches to maintain energy levels. Energy is expended when moving, and animals will move toward nearby vegetated patches first, in the same direction as most nearby animals second, or in a random direction third. Animals die from old age, starvation (low energy levels), intrinsic causes of death (illness, accident, animal predation), or by human predation.

Forager entities are considered to be forager bands (coresidential groups of related human hunter–gatherers) that collectively have an energy level that is expended when moving and replenished by hunting prey animals. Energy level affects forager hunting strategy and, optionally, low energy can lead to forager starvation and death. Foragers will move toward nearby prey animals first, or in a random direction if no animals are nearby. Foragers must hunt to obtain food energy from prey animals, but they also expend energy to pursue and process prey. Hunting prey animals is one of the most important key processes in the model. After a successful hunt, foragers may pause movement until all excess food energy is consumed, or they can leave the excess behind and continue to forage. Forager entities employ a foraging decision engine to choose whether to, and which animal to hunt. This decision logic is inspired by, but does not exactly replicate, classical optimal foraging theory (OFT) forging models. OFT has long been employed in archaeology as a baseline for understanding how human foraging behavior would play out under controlled circumstances with perfect knowledge of conditions and complete rationality (Gremillion, 2002; Codding and Bird, 2015). These restrictions have been fodder for many critiques (e.g., Keene, 1983), but we use the ABM formalism to mitigate the strict assumptions of classical OFT by embedding optimality decision logic in a changing, dynamic environment where foragers have limited knowledge about current conditions.

Landscape patch entities exist in an abstract square 211 × 211 gridded landscape universe (44,521 patches) that wraps at the edges (like a game of Pac-Man). Patches can be vegetated or unvegetated, and vegetation can be removed from vegetated patches by the grazing action of animal entities. Vegetation regrowth is another key process of the model. Vegetation can regrow on grazed patches after a specified number of time steps. Optionally, regrowth can be constrained to only a portion of the year using a seasonality variable. The proportion and spatial patterning of vegetated to non-vegetated patches can be altered, and initially non-vegetated patches can be set as persistent non-vegetated areas.

The temporal and spatial scales of the model are not hardcoded and are determined in relation to the way that parameters are set in relation to input data about prey animal life histories, environmental variability and dynamics, and human foraging strategies. Time units used in NetLogo are “ticks,” and it is most useful to scale these to a finite proportion of an annual cycle of vegetation growth. Although the landscape patches employed in the model are unitless, because animals move from one patch to an adjacent patch, it is most useful to consider their scale based on approximations of daily distance traveled by the particular prey species being modeled.

Simulating extinction dynamics with the MHPM

Background: the extirpation of Syncerus antiquus in the Cape Floristic Region of South Africa

We conducted a set of general modeling experiments with the MHPM loosely inspired by the case of the extinction of Syncerus antiquus in the Greater Cape Floristic Region in South Africa during the Late Quaternary. This case is a useful starting place for employing our modeling approach, because it provides an opportunity to simulate the confluence of interesting human and environmental factors related to megafauna extinction dynamics. Syncerus antiquus remains have been found in southern, eastern, and northern Africa, dating as far back as 1 Ma (Rossouw, 2001). The population of S. antiquus in South Africa became extirpated from the Cape Floristic Region by 12–10 ka (Faith, 2014) along with a handful of other large, grazing ungulates, including the Cape zebra (Equus capensis), a giant wildebeest (Megalotragus priscus), and a springbok (Antidorcas australis) (Klein, 1980; Rossouw, 2001; Faith, 2011a; Steele and Klein, 2013). Within the Cape Floristic Region, S. antiquus was present for ca. 1 Ma prior to its extirpation at the end of the Pleistocene (Faith, 2013; Steele and Klein, 2013). Although there is no definitive consensus on the drivers of S. antiquus extinction, hypotheses about this extinction have revolved around human exploitation of the species as a prey animal, or significant climatic and environmental shifts that occurred around the same time as the species was extirpated (Klein, 1980; Faith, 2014).

Anthropogenic explanations of this extinction event began with Klein's (1975, 1976) analysis of S. antiquus remains at the Klasies River Mouth site, which revealed a pattern of newborn (or late gestational) individuals of S. antiquus alongside mature adults. This stands in contrast to other large bovids at the site, whose fossil records show a more evenly distributed survivorship curve (Klein, 1975, p. 266). Klein hypothesized that this represented an intentional hunting pattern that targeted vulnerable individuals to buffer costs associated with hunting larger, dangerous males (Klein, 1975, 1976, 1978), and suggested that, “such a pattern could have endangered the survival of the species” (Klein, 1975, p. 267). The slower life histories of megafauna in general contribute to the plausibility of this hypothesis (e.g., Brook and Johnson, 2006).

Evidence of hafted microliths and spears appear early in the archaeological record of the region (Lombard, 2005; Wurz and Lombard, 2007; Wilkins et al., 2012; Schoville et al., 2016), and large bovids are frequently present in faunal assemblages at archaeological sites across time (Faith, 2011a; Reynard and Henshilwood, 2017), including direct evidence of human predation of mature S. antiquus from the Middle Stone Age layers at the Klasies River Mouth site (Milo, 1998). Therefore, hunters living during the Pleistocene–Holocene transition at the South African Cape were likely capable of taking large game such as S. antiquus (Reynard and Wurz, 2020). Additionally, human population density likely grew over the Pleistocene–Holocene transition (Clark and Kandel, 2013; Hallinan and Parkington, 2017), which would have amplified any negative effects of human predation on the S. antiquus population.

Other analyses, however, place greater weight on the significant climatic shifts that occurred across the Pleistocene–Holocene transition as part of a process extending back 1 Ma (Faith, 2011a, 2014). Widespread evidence suggests there was a decline of grassy forage in the region at the end of the Pleistocene that likely coincided with increased aridity and a decline in the productivity of the remaining grasslands (Klein, 1978; Schweitzer and Wilson, 1982; Cowling et al., 1999; Parkington et al., 2000; Faith, 2013; Klein and Cruz-Uribe, 2016; Chase et al., 2017, 2019; Faith et al., 2019). These climatic events likely had significant effects on the size and density of S. antiquus herds, and regression analysis of ungulate abundance data from late-middle Pleistocene through the Early Holocene in the Cape Floristic Region confirm this (Faith, 2011b). There was an overall reduction in biorichness due to changes in primary production, which may have combined with rising sea levels and localized declining precipitation to greatly diminish the paleo-Agulhas plain ecosystem. During the MSA, this expansive grassland ecosystem likely supported large numbers of grazing ungulates traveling north–south in response to seasonal rainfall patterns (Copeland et al., 2016). The loss of this plain was likely detrimental to ungulate biorichness in the Cape Floristic Region (Faith, 2011a; Reynard and Wurz, 2020). In other parts of the Cape, there was a significant reduction in seasonality at this time, characterized by a shift to year-round rainfall regimes (Sealy et al., 2016, 2020; Chase et al., 2017; Faith et al., 2019).

Comparisons to their still-extant close relative, the Cape buffalo (Syncerus caffer), and analysis of co-occurring species in faunal assemblages suggest that S. antiquus was a social grazer that sought grassy fodder in expansive, open environments (Reynard and Wurz, 2020). Reynard and Wurz (2020) suggested that the significant differences in ungulate diversity and richness among browsers, mixed-feeders, and grazers (like S. antiquus) observed in the archaeological record across time indicate ecological shifts among grassy, bushy, and mosaic environments. This suggests that S. antiquus may have been vulnerable to ecological changes that constrained the availability of open grasslands, causing them to either die out or migrate. Syncerus caffer prefers patches with higher green-leaf content, which meets their protein (nitrogen) needs (7–8% of diet) (Hildebrandt, 2014). This is particularly common when lower-quality food is available (Winnie et al., 2008). However, as large ruminants, S. caffer often engage in bulk grazing—and appear less picky—when foraging on higher quality patches that are “good enough,” suggesting their diets are more constrained by the protein content of their fodder (Prins, 1989; Ryan, 2006; Winnie et al., 2008). Additionally, it is important to note that S. caffer mortality tends to be influenced more by resource limitation than predation (Sinclair et al., 2003). Multiple lines of paleodietary evidence indicate that S. antiquus was a grazer with a similar diet to S. caffer (Faith, 2014). More details regarding estimations of the life history of S. antiquus are available in the Supplemental Text, Section 3, but this reconstruction combined with the potential types of environmental changes outlined above suggests that the extirpation of S. antiquus from this region could also plausibly have been associated with longer (but less productive) growing seasons and reduced availability of grassy fodder that decreased its grazing niche. Environmental and human hunting hypotheses are not necessarily at odds with one another, however, which presents an opportunity to employ a simulation approach to clarify how human hunting preferences and climate-driven environmental dynamics in this ancient SES could have worked, either in tandem or individually, to push or pull the South African population of S. antiquus towards an extinction threshold.

Formulating experiments in the MHPM

We parameterized a set of 24 modeling experiments (Table 3) inspired by the types of potential extinction dynamics that may have occurred in the greater Cape Floristic Region during the Late Quaternary using different combinations of values for environmental patchiness, growing season length, male prey handling cost, and forager population size. Three forager population configurations include scenarios with no hunters (init-foragers = 0), scenarios with a moderate hunter population (init-foragers = 10), and scenarios with a larger population of hunters (init-foragers = 20). Four different environmental configurations are parameterized by two possibilities for the intrinsic patchiness of the environment (grass-proportion = 0.33, grass-proportion = 0.66) and two possibilities for the length of the growing season (seasonality = 0.4 yr, seasonality = 0.8 yr). Two possibilities for male prey handling costs include aggressive males with higher handling costs than females (processing-cost-males = 6) or less-aggressive males with handling costs equal to those of females (processing-cost-males = 3). We set the other parameters of the model (Table 2) in relation to an approximate reconstruction of S. antiquus life history as well as archaeological and paleoecological data from the Late Quaternary of the region to help situate these hypothetical experimental scenarios in a reasonably realistic setting. We explain these parameterizations in the remainder of this section.

Table 2. Values of all variables used in the experiments presented here. Baseline variable values were not changed, but experimental variables were allowed to change through the sets of values listed here for a total of 20 functionally distinct experiments

The MHPM can be computationally expensive, so we needed to balance the spatio-temporal complexity of the model runs to achieve an acceptable level of detail in a reasonable amount of runtime. We set the temporal scale of our experiments to a daily granularity (ticks-per-year = 365 ticks). The spatial scale of the model is relative. If S. antiquus displayed average walking speeds and daily distance movements like those reported for S. caffer (~3.4 km/day, Supplemental Text, Section 3) (Ryan and Jordaan, 2005), then the spatial scale of a single patch in our model would be roughly equivalent to ~911 km², and the entire world would be ~500,000 km². These values are relative and speculative, however, because the actual daily movement speed of S. antiquus is not known. This size of the model world, however, provides ample space for large herds to develop and self-organize within distinct seasonal home ranges. We set the forager search radius to twelve patches, which, using the relative values from above, would translate to a “real world” value of not more than 40 km. It is important to reiterate that the “search radius” of the MHPM is not the same as the “foraging radius” of the central place foraging model (e.g., Kelly, 1995), but instead delineates the local environment for which a forager agent would have some general current knowledge about animal densities while making foraging decisions.

Syncerus antiquus was larger than S. caffer, with estimates of adult weight varying by methodological approach, including 1300–1800 kg (Klein, 1976), ~1,000 kg (Smith et al., 2003), or ~850 kg (Bibi and Kiessling, 2015). Syncerus antiquus was likely a highly sexually dimorphic species with males standing ~2 m tall with broad horns spanning 3 m or more in width (Klein, 1976). Recalling that food value and handling costs are highly generalized in the MHPM, we tuned these parameters to reflect differences in food value and processing costs inherent to a sexually dimorphic species (like S. antiquus) and held male prey handling values as an independent variable related to their potential for aggressiveness. Recalling that energy units are arbitrary in the MHPM (see the ODD protocol in the Supplemental Material for detailed equations and descriptions of variables and units) and that the maximum forager-energy-value is 100, in these experiments we considered that a large-bodied prey animal would provide substantial payoff. Thus, male prey contribute an intrinsic food value of 30 with a processing cost of either 3 or 6, for a total payoff of either 27 or 24 energy units. Harvested female prey contribute an intrinsic food value of 27 with a fixed processing cost of 3, for a total payoff of 24 energy units. Pregnant or recently postpartum female prey provide an additional food value of 3 when harvested, which means that females with calves will have a combined payoff of 27 energy units, which is equal to or greater than the payoff for males, depending on male processing costs. This means that the initial intrinsic prey rankings, as well as the outcome of the instantaneous prey choice decision model, will differ across experiments. Recalling that forager agents are (1) conceptually bands rather than individuals, (2) we intentionally kept the forager agents simple, (3) other food items are ignored in the model, and (4) forager energy value is measured on an abstract scale of 0–100, a single successful hunt provides about a quarter of the maximum forager energy value. The forager movement cost is 1, which means that forager bands need to harvest at least one animal per month in order to maintain peak energy values. Forager agents will remain stationary after a successful hunt until the entire food value is consumed. These parameterizations are comparable to data on big game hunting styles and frequencies of successful kills among contemporary African hunter–gatherers (Hawkes et al., 1991).

Prey life history variables were scaled to the temporal granularity of 365 ticks per year and parameterized to be roughly equivalent to the reconstruction of S. antiquus life history through comparison to S. caffer and information about other extinct large bovines (Supplemental Text, Section 3). We set the prey age of senescence to 7300 ticks (20 years) and the prey maximum lifespan to 10,220 ticks (28 years) (~30% longer than values reported for S. caffer, Supplemental Text, Section 3). We set the prey birth spacing interval to 730 ticks (two years), which approximates a seasonal reproductive schedule where females will, on average, produce a new calf every two to three years (~30% longer than the average birth spacing interval for S. caffer, Supplemental Text, Section 3; Supplemental Table 1, Supplemental Fig. 2). We set the initial population of prey as 400 individuals split evenly between males and females, providing approximately 125 ha per individual (Supplemental Text, Section 3). We set an external mortality rate of 0.0025 (~1 death per 400 animals), and an internal mortality of 0.025 (~1 death per 40 animals). The maximum animal energy is 100 in the MHPM. We set the animal starvation threshold to 25 energy units, the animal movement energy cost to 1, and the energy gain for a grazed patch of grass to 1.5 energy units. This means that a prey agent with a full energy value of 100 could become resource stressed after 75 days of continuous movement without grazing returns but could regain full energy in about half a year of continuous grazing returns. Altogether, these dietary and life history parameters interact to produce reasonable population dynamics, including senescence curves and population growth rates. We set the length of each experiment to 73,000 ticks, allowing for population dynamics to play out over 200 years, which is ~10 sequential lifetimes of prey.

We conducted 40 repeated runs of each scenario (see Table 3) to characterize the inter-run variability of each scenario (Romanowska, 2015). This resulted in a total of 960 model runs. We used the BehaviorSpace tool in NetLogo to set up and execute these experiments on a high-performance computer with 88 independent compute nodes and 270 gigabytes of RAM. Temporal dynamics in male, female, and postpartum-female populations were saved to delimited text files for each run, and we used the Pandas and Seaborn data organization, analysis, and visualization packages in Python to parse, analyze, and plot results across all repetitions for each experiment.

Table 3. Matrix of the 16 experiment scenarios implemented in this research. Each scenario is labeled by a unique numerical code that combines the slider values for the four foraging scenarios and the four types of environmental conditions, as denoted by the row and column descriptors

Results

We measured the probability that each experimental scenario could lead to localized extinction, which we defined as a prey animal population of 0 at any point before the full 73,000 ticks. Sensitivity analysis, reported in Supplemental Text Section 4, indicates that in the absence of human predation, prey populations are stable unless the environment is extremely patchy (grass-proportion < 0.2, Supplemental Figs. 3–5). Demographic instability in these conditions is exacerbated by a shortened growing season, particularly when the growing season is less than half a year (seasonality < 0.5, Supplemental Figs. 3–5). Adding human foraging pressure to the simulation changes these dynamics. Figure 1 shows the percentage of repetitions in each model scenario that ended in an extinction event, Figure 2 plots the timing of those events, and Table 4 summarizes these results. Across environmental scenarios, the main determinant of the timing of extinction events is male prey processing costs. When males are low cost (foraging scenarios 10-3 and 20-3, Figs. 1a–d and 2a–d), prey populations endure for significantly longer times regardless of forager population size. With small forager populations and low male handling costs (foraging scenario 10-3, Figs. 1a–d and 2a–d), no extinctions occurred in the less-patchy grassland environment scenarios (environmental scenarios 0.66-0.4 and 0.66-0.8, Figs. 1c, d and 2c, d), and few extinctions occurred in the patchy grassland environment scenario when the growing season was long (environmental scenario 0.33-0.8, Figs. 1b and 2b). When male prey handling costs are high (foraging scenarios 10-6 and 20-6, Figs. 1a–d and 2a–d), extinction always occurs quickly, regardless of environmental conditions. These poor environmental conditions, exemplified by high patchiness coupled with a short growing season (environmental scenario 0.33-0.4, Figs. 1a and 2a), elicited the greatest number of extinction events, which also occurred earlier in the simulation timeline than other scenarios. No extinctions occurred in the scenarios with no hunters (foraging scenarios 0-3 and 0-6, Fig. 1a–d and 2a–d).

Figure 1. Percentage of scenario repetitions ending in an extinction event for the different foraging scenarios tested in this paper. Subplots in the top row represent scenarios with sparse patchy grasslands; subplots in the bottom row represent scenarios with fuller, more continuous grasslands. Subplots in the left column represent scenarios with a short growing season; subplots in the right column represent scenarios with a long growing season. Foraging scenarios appear along the x-axes and are coded as a combination of forager population size (number to left of dash) and male buffalo processing cost (number to right of dash; see Table 3 for full experiment codes).

Figure 2. Timing of extinction events in each modeled scenario of human predation of prey across all repetitions of each simulation scenario. Subplots in the top row represent scenarios with sparse patchy grasslands; subplots in the bottom row represent scenarios with fuller, more continuous grasslands. Subplots in the left column represent scenarios with a short growing season; subplots in the right column represent scenarios with a long growing season. Foraging scenarios appear along the x-axes and are coded as a combination of forager population size (number to left of dash) and male prey processing cost (number to right of dash; see Table 3 for full experiment codes). Each dot represents the time (in number of ticks) where prey populations dropped to zero in a simulation run. The median value and 99% confidence interval (determined by bootstrap) of time to extinction in each scenario is overlain in black.

Table 4. Number and timing of extinction events occurring in each model scenario. See Table 3 for explanation of the Experiment ID codes

We tracked demographic variability in the prey populations within each scenario by recording the proportion of males to females at each tick (Fig. 3). This illustrates how male processing costs are the primary factor driving differences in hunting strategies. When males have low processing costs (foraging scenarios 10-3 and 20-3, Fig. 3b, c, upper row of subplots) hunters prefer males but exploit all three prey “types” (males, females, and postpartum females, Supplemental Figure 6), resulting in a female-skewed sex ratio (Fig. 3b, c), these scenarios only lead to extinction under poor environmental conditions (Figs. 1a–d and 2a–d, Table 4) because female-skewed herds still maintain high reproductive capacity. This pattern is reminiscent of a traditional “overkill” scenario, where the overall quantity of animals killed is the main anthropogenic driver of extinction. When males have high processing costs (foraging scenarios 10-6 and 20-6, Fig. 3b, c, lower row of subplots), hunters will always prefer postpartum females because the payoff is the highest. This results in a male-skewed sex ratio (Fig. 3b, c, lower row of subplots), which, in our model, always leads to extinction because of the concomitant reduction in reproductive capacity of the prey population that diminishes the ability of the herd to replace individuals lost to human predation. This differs from a purely “overkill” scenario because extinction is driven by feedback between prey demography and hunting strategies that culminate in rapid extinction when one sex is preferred over the other. Scenarios with no hunting (Fig. 3a) resulted in evenly balanced sex ratios across all environmental scenarios.

Figure 3. Scatterplots showing variability in the ratio of male prey to female prey in each simulated foraging scenario. Each dot represents a specific time-step (tick) in one of the repetitions of each scenario. The black diagonal dashed line represents parity between the number of males and females. Points falling above this line indicate time-steps with an excess of females; points falling below this line indicate time-steps with an excess of males. The left two columns represent scenarios with sparse patchy grasslands; the right two columns represent scenarios with fuller, more continuous grasslands. The first and third columns represent scenarios with a short growing season; the second and fourth columns represent scenarios with a long growing season. Each row represents a different simulated foraging scenario, categorized as a combination of forager population size (number to left of dash) and male prey processing cost (number to right of dash; see Table 3 full experiment codes). A 1% sample of all data points is displayed so that the smaller nuances of the data patterns are less obscured by plot marker overlap.

We estimate the environmental effect of prey animal herds with an “overgrazing coefficient,” which is the ratio of the number of grazed patches to the total possible number of grass-covered patches for each scenario (Fig. 4). Human hunting keeps prey populations below the emergent carrying capacity of the equivalent environmental conditions when there is no human predation (Supplemental Fig. 4). This results in lower overgrazing coefficients for all scenarios with any amount of human hunting pressure. Scenarios with small forager populations and low-cost males (foraging scenario 10-3) resulted in average overgrazing coefficients like those of the control runs, especially in the highly seasonal but well-vegetated environmental scenario (environmental scenario 0.66-0.4). In general, the highly seasonal environments (environmental scenarios 0.33-0.4 and 0.66-0.4) resulted in more highly overgrazed landscapes compared to the less-seasonal ones (environmental scenarios 0.33-0.8 and 0.66-0.8), indicating that seasonality reduces an important environmental buffer that limits the resilience of herds recovering from human predation.

Figure 4. The “grazing coefficient” (ratio of grazed patches to total initial grass patches) for each model scenario and control model. Red solid line is the mean value for all ticks of all repetitions. The dotted lines represent an error envelope at the 1σ level. The dashed line represents the maximum overgrazing coefficient recorded for each scenario. The x-axis is labeled by experiment codes, which are explained in Table 3.

Additional results are available in Supplemental Text, Section 5, where we report the number and sex of animals foraged per tick (including a separate tally of foraged postpartum female animals, Supplemental Fig. 5), the total number of prey animal births and deaths that occurred in each tick of each run (Supplemental Fig. 6), the total number of males, females, and postpartum females at each time-step of each run (Supplemental Figs. 7, 8), and the average animal energy levels at each time-step (Supplemental Fig. 8). The patterns in these additional results align with those reported above. When male processing cost is high, postpartum females become the preferred prey (Supplemental Fig. 6). Scenarios that include human hunting have an overall lower prey animal birth-to-death ratio when compared to scenarios without humans (Supplemental Fig. 6). In these situations, prey populations are also smaller (Supplemental Figs. 8, 9), and seasonal energy dynamics in the overall size of prey populations are dampened (Supplemental Fig. 8).

Discussion

Feedback and boundary conditions

The results of our modeling experiments show that when modeled as a complex set of feedbacks, extinction is not necessarily a forgone conclusion of any specific human behaviors or of any specific sets of environmental conditions. Rather, we can characterize combinations of socio-environmental conditions where the probability of extinction is more and less likely. Extinction itself occurs when feedback processes self-amplify within the system to set prey animal populations on precipitous pathways of decline without the possibility of recovery. In other words, under many system configurations, extinction can happen, but does not always happen. But extinction is much more likely to occur under some system configurations than others. In CAS theory, the size and configuration of areas of uncertain outcome are related to the adaptive capacity of the SES (Preiser et al., 2018). This adaptive region is often circumscribed by permeable, dynamically derived, boundary conditions, across which the probability of different system trajectories (e.g., extinction) are more and less definite, but which are not necessarily determinable a priori because they tend to be emergent properties of the system (Preiser et al., 2018). Here, we explore the sets of environmental, anthropogenic, and cross-level feedbacks that may control extinction dynamics and the emergence of boundary conditions under the set of social-environmental conditions that we modeled.

Environmental patchiness and seasonality interact to influence prey animal dispersal and probability of starvation. These dynamics affect the overall likelihood of finding mates and successfully reproducing across the different scenarios, thereby influencing the resilience of the population of prey animals to additional deaths from human hunting. When the environment is patchy and highly seasonal (scenario 0.33-0.4), it becomes harder to find mates because the dwindling population is spread thinner across the landscape. We consider this to be strong negative feedback that increases extinction risk over time, especially if deaths due to starvation begin to outweigh births. Comparing the results of the sensitivity analysis (Supplemental Fig. 3) with the experiment results (Figs. 1, 2), we see that this value (scenario 0.33-0.4) is lower without human predations (grass-proportion = 0.2 when seasonality < 0.6 yr). This is an emergent boundary condition that increases the probability of extinction due to environmental risk within only a few years.

On the other hand, when vegetation is patchy, but the growing season is longer (environmental scenario 0.33-0.8), the prey population is scattered and it is easier for individual animals to find grazable patches over most of the year. In this environmental scenario, fewer animals die of starvation, and they can still find mates with enough frequency to maintain a stable, if not maximally large population. We consider this to be “neutral” environmental feedback that neither increases nor decreases extinction risk over time.

In a less-patchy environment with a short growing season (environmental scenario 0.66-0.4), prey populations tend to coalesce into emergent herds that follow the seasonal regrowth of grasses across the landscape. Although grazable grass can become seasonally scarce in this scenario, it occurs in contiguous patches and most animals can consume enough fodder to prevent excessive deaths due to starvation. The density of animals in the emergent “herds” means that they are also able to easily locate potential mates and reproduce. This maintains a large, stable, but not maximal population. We consider this environmental scenario to have a slightly positive feedback effect.

Finally, when grasslands are extensive and the growing season is long (environmental scenario 0.66-0.8), fodder is always easy to find, and animals rarely enter into starvation conditions. Mates are easy to find because the overall population remains high, even if not always spatially clumped. Here, maximum population sizes are achieved, and we consider this to be a positive feedback effect. Again, comparing the experiment results (Figures 1, 2) to the sensitivity analysis results (Supplemental Fig. 3), we see that in the absence of human predation, prey populations stabilize when grass-proportion > 0.2 and seasonality ≠ 0, which we consider to be another environmental boundary condition, across which the risk of extinction from environmental pressure is minimal.

Foraging scenarios that led to preferential hunting patterns targeting postpartum females (foraging scenarios 10-6 and 20-6) enhance the chance and speed of extinction events across multiple replications, indicating a strong negative feedback process (Figs. 1a–d, 2a–d). This is an anthropogenic boundary condition of the SES that emerges quickly and leads to extinction under all but the most favorable environmental conditions. In large k-selected species such as S. antiquus and other megafauna, interruptions to breeding can have drastic short- and long-term effects. Their slow life histories (particularly their low fecundity and slow maturation) make them particularly vulnerable to human predation (Johnson, 2002; Brook and Johnson, 2006). When the time between subsequent births is long, the loss of even a single breeding female reduces the overall reproductive rate significantly, especially when the population is small. This becomes a strong negative feedback effect when breeding females are targeted systematically as part of a human adaptive strategy. Scenarios with large forager populations (foraging scenarios 20-3 and 20-6) also enhanced the probability of extinction, but did so because of over-harvesting, rather than (or in addition to) disrupting the breeding cycle (Figs. 1a–d, 2a–d). This is another anthropogenic boundary condition, which aligns closely to the classic “overkill” hypotheses of human-caused extinction (e.g., Martin, 1990). On the other hand, when a small forager population size combined with hunting strategies where males and females were equally targeted (foraging scenario 10-3), human hunting had a minimal effect on population dynamics unless the environment was both highly seasonal and patchy (Figs. 1a, 2a). This indicates a neutral to weak negative feedback loop and suggests the presence of another boundary condition across which human predation of large prey could be considered sustainable.

In a SES, socio-environmental feedback can alter the boundary conditions in ways that cannot be predicted if the social and natural components are only considered independently. When we examine how feedbacks from human hunting pressures combine with those from environmental variation across the modeled scenarios, we see distinct patterns in extinction rates. This suggests that a new set of cross-feedbacks emerged from these socio-natural interactions that are more than simply the sum of the individual environmental and anthropogenic effects. For example, our experiments show how the addition of human hunting in general may reduce the range of sustainable environmental configurations for prey herds (compare Figs. 1 and 2 with Supplemental Fig. 3). This potentially introduces new boundary conditions, or extends or foreshortens those that already exist. In other words, cross-system feedback produces non-linear combination effects that amplify or dampen existing extinction risks. Alongside other research (e.g., Braun et al., 2021), these results contribute to our understanding of deep-time socio-natural dynamics that were keystone to human evolution.

Insights into the dynamics of Late Quaternary megafauna extinctions

Our model can contribute to broader understandings of Late Quaternary megafauna extinction by providing a view of how extinction dynamics potentially play out in changing conditions. As we discussed in the introduction to this paper, the inductive models that archaeologists generally build often struggle to adequately fill in these dynamics, which are never directly preserved in the archaeological record. Formal modeling approaches, such as the ABM approach we have employed here, offer a robust alternative that stands to help move long-standing debates forward. For example, our experiments show that overkill hypotheses, which place blame squarely on human hunters, are possible, but are only likely to be correct when a specific set of boundary conditions are surpassed. Recent empirical research corroborates this insight, but also begs the need for additional formal modeling studies. For example, Bradshaw et al. (2021) and Stewart et al. (2021) both failed to confirm simple correlations among human population growth, hunting susceptibility, and extinction. Louys et al.'s (2021) examination of early insular megafauna extinction events conceded that humans could have played a role in island extinctions, but their deep-time analysis painted a more complex picture of human–animal interactions. Similar work by Boivin et al. (2016) and Jukar et al. (2021) showed that large animals may potentially coevolve with humans in areas with prolonged human–megafauna interactions, such as Africa and India, where longer time periods of interaction correlate with a greater number of surviving species of megafauna. In another region with a very long history of human occupation, Hocknull et al. (2020) suggested that the extinction of northeastern Australian megafauna 40 ka ago resulted from complex interactions of environmental variables (hydroclimate deterioration combined with prolonged environmental change). And yet, Brook and Johnson's (2006) examination of megafauna decline in Australia suggested that the slow life history of the largest known marsupial, Diprotodon optatum, meant that even a slow rate of exploitation of juveniles (0.2 kills/yr per capita) would be sufficient to catalyze extinction within a few centuries. Similarly, in southeastern Australia, Saltré et al. (2019) found that hunting pressure from small human groups in search of drinkable water may have interacted with environmental pressure accelerating localized extinctions of megafauna. Again, these inductive approaches provide tantalizing glimpses into past worlds, but do little to help us understand how, why, when, and how likely it is that extinctions will occur in different SES across time and space.

Our modeling results also inform the relative effects of climatic and environmental fluctuations in megafauna extinction dynamics. In the case study of S. antiquus, the correlation of extinction of species with the decline of grasslands in the greater Cape Floristic Region has led many researchers to assume that environmental change was a major causal factor (Klein, 1980; Faith, 2014). During the Pleistocene–Holocene transition in the Cape Floristic Region, it does appear that grasslands contracted while closed habitats expanded (Klein, 1974; Faith, 2011b). At the same time, however, the ratio of C4 to C3 grasses was increasing, and isotopic analysis of S. antiquus teeth recovered from Nelsons Bay shows that they were increasingly focusing on C4 plants across this transition (Sealy et al., 2020). Sealy et al. (2020) suggested the shift from C3- to C4-dominant grasslands in the area was due to a change from year-round rainfall at the last glacial maximum to more seasonal (likely summer) rainfall in the Early Holocene.

These different possible environmental conditions are captured in the four environmental scenarios used in our experiments, which show that it is possible that contracting habitat alone may lead to prey animal extinction. But, again, this is not a foregone conclusion. Instead, when deleterious anthropogenic and environmental conditions amplify due to cross-feedback, the probability of crossing an extinction boundary condition increases rapidly. In a complimentary example, Seersholm et al. (2020) found that DNA dating to the Younger Dryas period in Hall's Cave, Texas demonstrated that although the contracting grasslands correlated with decreased biodiversity, taxa that faced direct human predation did not survive the Pleistocene–Holocene transition (including several species of megafauna). The same is not true of taxa that were sheltered from human predation (such as rodents and many plants). These empirical results largely align with the behavior of our model.

Conclusion

Extinctions are multifaceted, complex events that are highly influenced by risk factors and feedback processes that emerge as a consequence of system dynamics over time. Any successful explanation for a particular extinction event must be able to satisfactorily answer three basic questions: (1) Why did this animal go extinct, and not others? (2) Why did extinction occur at this time, and not at other times? (3) Why did extinction happen in this location, and not others? The proxy nature of the archaeological and paleoecological records makes conclusive answers to these questions difficult to find. But, even if they could be answered conclusively, the answers may not be particularly useful if we want to translate archaeological knowledge and data about past extinctions towards understanding (and perhaps even preventing) current and future anthropogenic extinctions. Formal deductive models, such as the MHPM, can, however, elucidate the boundary conditions of systems and illuminate the feedback pathways that make extinction more or less likely occur in past SESs. When carefully designed and parameterized in reasonable and meaningful ways, modeling can help us learn about the inevitability (or not) of extinction in different scenarios, which is knowledge that could be directly put to use in current conservation situations (Johnson et al., 2005; Hjelle et al., 2012; Dearing et al., 2014; Barton et al., 2016).

But much like the earlier calls for interdisciplinary efforts to engage archaeologists and conservationists around the concept of historical ecology (Balée, 2006; Szabó, 2010; Hjelle et al., 2012), modelers and empiricists must work together to make this happen. We have shown how models and modeling experiments can be inspired by real archaeological cases. To complete the loop, hypotheses about potential past boundary conditions that emerge from modeling experiments should inspire new archaeological questions and research. In our case study of the early Holocene extinction of S. antiquus, four important boundary conditions in the greater Cape Floristic Region during the Pleistocene–Holocene transition warrant renewed investigation by field scientists: (1) the extent of grassland fragmentation, (2) the extremity of the changes in rainfall seasonality, (3) the sex ratios of species found in faunal assemblages, and (4) the density of human foragers in the region. By combining new empirical data with the insights from our modeling experiments, we may be able to characterize these boundary conditions with observed data so that, following the work of Nielsen et al. (2021), they can become actionable knowledge that can be used to combat the current widespread loss of biodiversity. In this way, we can move beyond the practice of telling archaeological “just so” stories about the past, and towards a future of archaeology as part of a multidisciplinary approach to solving contemporary and future problems.