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., Reference Ceballos, Ehrlich, Barnosky, García, Pringle and Palmer2015). It is estimated that 48% of animal species are currently in decline (Finn et al., Reference Finn, Grattarola and Pincheira-Donoso2023) and that 50% of the plant and animal species on Earth may be lost by 2100 (Braje and Erlandson, Reference Braje and Erlandson2013, 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., Reference Andermann, Faurby, Turvey, Antonelli and Silvestro2020). 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, Reference Braje and Erlandson2013; Boivin et al., Reference Boivin, Zeder, Fuller, Crowther, Larson, Erlandson, Denham and Petraglia2016), but the role of human-driven effects versus climate-driven effects in these extinctions continues to be debated (Barnosky et al., Reference Barnosky, Koch, Feranec, Wing and Shabel2004; Koch and Barnosky, Reference Koch and Barnosky2006; Braje and Erlandson, Reference Braje and Erlandson2013).
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., Reference Lyons, Smith, Wagner, White and Brown2004; Turney et al., Reference Turney, Flannery, Roberts, Reid, Fifield, Higham and Jacobs2008; van der Kaars et al., Reference van der Kaars, Miller, Turney, Cook, Nürnberg, Schönfeld, Kershaw and Lehman2017; Saltré et al., Reference Saltré, Chadoeuf, Peters, McDowell, Friedrich, Timmermann, Ulm and Bradshaw2019). If we want these findings to be relevant to current-day issues (i.e., “action archaeology” sensu Sabloff, Reference Sabloff2016), 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., Reference Faith, Rowan, Du and Barr2020). 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, Reference Epstein1999; Bankes, Reference Bankes2002) informed by a complex adaptive systems (CAS) perspective about social–ecological systems (SES) (Preiser et al., Reference Preiser, Biggs, De Vos and Folke2018) 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, Reference Epstein1999; Bankes, Reference Bankes2002; Romanowska, Reference Romanowska2015; Wilensky and Rand, Reference Wilensky and Rand2015). 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, Reference Daems2021). The components of these systems change and reorganize themselves to adapt to their surroundings (Holland, Reference Holland1992). There is no central all-controlling mechanism that guides the system's behavior (Mitchell, Reference Mitchell2009). Rather, complex, self-organizing behavior emerges out of simple operational rules and complex feedback interactions (Kohler, Reference Kohler2012). SES are a special type of CAS where human social system components connect to and interact with ecosystem components (Walker et al., Reference Walker, Gunderson, Kinzig, Folke, Carpenter and Schultz2006; Glaser et al., Reference Glaser, Ratter, Krause, Welp, Glaser, Ratter, Krause and Welp2012; Preiser et al., Reference Preiser, Biggs, De Vos and Folke2018).
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, Reference Box and Draper1987), 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., Reference Robinson, Brown, Parker, Schreinemachers, Janssen, Huigen, Wittmer, Gotts, Promburom and Irwin2007; Rounsevell et al., Reference Rounsevell, Robinson and Murray-Rust2012; Railsback and Grimm, Reference Railsback and Grimm2019). 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, Reference Railsback and Grimm2019; Romanowska et al., Reference Romanowska, Wren and Crabtree2021). 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, Reference Wilensky1999). The MHPM is described in detail using the Overview, Design Concepts, and Details (ODD) protocol (Grimm et al., Reference Grimm, Berger, Bastiansen, Eliassen, Ginot, Giske and Goss-Custard2006, Reference Grimm, Berger, DeAngelis, Polhill, Giske and Railsback2010, Reference Grimm, Railsback, Vincenot, Berger, Gallagher, DeAngelis and Edmonds2020) in the Supplementary Text, and summarized here following the Summary ODD format outlined by Grimm et al. (Reference Grimm, Railsback, Vincenot, Berger, Gallagher, DeAngelis and Edmonds2020).
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, Reference Gremillion2002; Codding and Bird, Reference Codding and Bird2015). These restrictions have been fodder for many critiques (e.g., Keene, Reference Keene, Moore and Keene1983), 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, Reference Rossouw2001). The population of S. antiquus in South Africa became extirpated from the Cape Floristic Region by 12–10 ka (Faith, Reference Faith2014) 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, Reference Klein1980; Rossouw, Reference Rossouw2001; Faith, Reference Faith2011a; Steele and Klein, Reference Steele and Klein2013). Within the Cape Floristic Region, S. antiquus was present for ca. 1 Ma prior to its extirpation at the end of the Pleistocene (Faith, Reference Faith2013; Steele and Klein, Reference Steele and Klein2013). 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, Reference Klein1980; Faith, Reference Faith2014).
Anthropogenic explanations of this extinction event began with Klein's (Reference Klein1975, Reference Klein1976) 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, Reference Klein1975, 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, Reference Klein1975, Reference Klein1976, Reference Klein1978), and suggested that, “such a pattern could have endangered the survival of the species” (Klein, Reference Klein1975, p. 267). The slower life histories of megafauna in general contribute to the plausibility of this hypothesis (e.g., Brook and Johnson, Reference Brook and Johnson2006).
Evidence of hafted microliths and spears appear early in the archaeological record of the region (Lombard, Reference Lombard2005; Wurz and Lombard, Reference Wurz and Lombard2007; Wilkins et al., Reference Wilkins, Schoville, Brown and Chazan2012; Schoville et al., Reference Schoville, Brown, Harris and Wilkins2016), and large bovids are frequently present in faunal assemblages at archaeological sites across time (Faith, Reference Faith2011a; Reynard and Henshilwood, Reference Reynard and Henshilwood2017), including direct evidence of human predation of mature S. antiquus from the Middle Stone Age layers at the Klasies River Mouth site (Milo, Reference Milo1998). 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, Reference Reynard and Wurz2020). Additionally, human population density likely grew over the Pleistocene–Holocene transition (Clark and Kandel, Reference Clark and Kandel2013; Hallinan and Parkington, Reference Hallinan and Parkington2017), 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, Reference Faith2011a, Reference Faith2014). 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, Reference Klein1978; Schweitzer and Wilson, Reference Schweitzer and Wilson1982; Cowling et al., Reference Cowling, Cartwright, Parkington and Allsopp1999; Parkington et al., Reference Parkington, Cartwright, Cowling, Baxter and Meadows2000; Faith, Reference Faith2013; Klein and Cruz-Uribe, Reference Klein and Cruz-Uribe2016; Chase et al., Reference Chase, Chevalier, Boom and Carr2017, Reference Chase, Boom, Carr, Chevalier, Quick, Verboom and Reimer2019; Faith et al., Reference Faith, Chase and Avery2019). 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, Reference Faith2011b). 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., Reference Copeland, Cawthra, Fisher, Lee-Thorp, Cowling, le Roux, Hodgkins and Marean2016). The loss of this plain was likely detrimental to ungulate biorichness in the Cape Floristic Region (Faith, Reference Faith2011a; Reynard and Wurz, Reference Reynard and Wurz2020). 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., Reference Sealy, Lee-Thorp, Loftus, Faith and Marean2016, Reference Sealy, Naidoo, Hare, Brunton and Faith2020; Chase et al., Reference Chase, Chevalier, Boom and Carr2017; Faith et al., Reference Faith, Chase and Avery2019).
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, Reference Reynard and Wurz2020). Reynard and Wurz (Reference Reynard and Wurz2020) 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, Reference Hildebrandt2014). This is particularly common when lower-quality food is available (Winnie et al., Reference Winnie, Cross and Getz2008). 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, Reference Prins1989; Ryan, Reference Ryan2006; Winnie et al., Reference Winnie, Cross and Getz2008). Additionally, it is important to note that S. caffer mortality tends to be influenced more by resource limitation than predation (Sinclair et al., Reference Sinclair, Mduma and Brashares2003). Multiple lines of paleodietary evidence indicate that S. antiquus was a grazer with a similar diet to S. caffer (Faith, Reference Faith2014). 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, Reference Ryan and Jordaan2005), then the spatial scale of a single patch in our model would be roughly equivalent to ~911 km2, and the entire world would be ~500,000 km2. 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, Reference Kelly1995), 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, Reference Klein1976), ~1,000 kg (Smith et al., Reference Smith, Lyons, Ernest, Jones, Kaufman, Dayan, Marquet, Brown and Haskell2003), or ~850 kg (Bibi and Kiessling, Reference Bibi and Kiessling2015). 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, Reference Klein1976). 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., Reference Hawkes, O'Connell and Blurton Jones1991).
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, Reference Romanowska2015). 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., Reference Preiser, Biggs, De Vos and Folke2018). 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., Reference Preiser, Biggs, De Vos and Folke2018). 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, Reference Johnson2002; Brook and Johnson, Reference Brook and Johnson2006). 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, Reference Martin1990). 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., Reference Braun, Faith, Douglass, Davies, Power, Aldeias and Conard2021), 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. (Reference Bradshaw, Johnson, Llewelyn, Weisbecker, Strona and Saltré2021) and Stewart et al. (Reference Stewart, Carleton and Groucutt2021) both failed to confirm simple correlations among human population growth, hunting susceptibility, and extinction. Louys et al.'s (Reference Louys, Braje, Chang, Cosgrove, Fitzpatrick, Fujita and Hawkins2021) 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. (Reference Boivin, Zeder, Fuller, Crowther, Larson, Erlandson, Denham and Petraglia2016) and Jukar et al. (Reference Jukar, Lyons, Wagner and Uhen2021) 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. (Reference Hocknull, Lewis, Arnold, Pietsch, Joannes-Boyau, Price and Moss2020) 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 (Reference Brook and Johnson2006) 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. (Reference Saltré, Chadoeuf, Peters, McDowell, Friedrich, Timmermann, Ulm and Bradshaw2019) 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, Reference Klein1980; Faith, Reference Faith2014). During the Pleistocene–Holocene transition in the Cape Floristic Region, it does appear that grasslands contracted while closed habitats expanded (Klein, Reference Klein1974; Faith, Reference Faith2011b). 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., Reference Sealy, Naidoo, Hare, Brunton and Faith2020). Sealy et al. (Reference Sealy, Naidoo, Hare, Brunton and Faith2020) 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. (Reference Seersholm, Werndly, Grealy, Johnson, Keenan Early, Lundelius and Winsborough2020) 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., Reference Johnson, Kohler and Cowan2005; Hjelle et al., Reference Hjelle, Kaland, Kvamme, Lødøen and Natlandsmyr2012; Dearing et al., Reference Dearing, Wang, Zhang, Dyke, Haberl, Hossain and Langdon2014; Barton et al., Reference Barton, Ullah, Bergin, Sarjoughian, Mayer, Bernabeu-Auban, Heimsath, Acevedo, Riel-Salvatore and Arrowsmith2016).
But much like the earlier calls for interdisciplinary efforts to engage archaeologists and conservationists around the concept of historical ecology (Balée, Reference Balée2006; Szabó, Reference Szabó2010; Hjelle et al., Reference Hjelle, Kaland, Kvamme, Lødøen and Natlandsmyr2012), 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. (Reference Nielsen, Marteau, Bauer, Bradbury, Broad, Burgess and Burgman2021), 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.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/qua.2024.6
Acknowledgments
We would like to thank C.M. Barton for making his Diet Breadth NetLogo model freely available for reuse and modification. M.C.K. and I.U. conceived and designed the research; I.U. coded the model and conducted the experiments; I.U. conducted the statistical analysis of model output; M.C.K. and I.U. interpreted the results; and M.C.K. and I.U. authored the paper. This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors. Computational resources in the Computational Archaeology Laboratory at San Diego State University were used to develop the MHPM and to run and analyze the simulations.
Software and Data Availability
The Megafauna Hunting Pressure Model was created in NetLogo 6.3.0 (Wilensky, Reference Wilensky1999) and is made freely available under the GNU GPL version 3. The code for this project is hosted at https://github.com/isaacullah/MegafaunaHuntingPressure. The static release of version 1.1 of the model was used for the experiments presented here and can be referenced with the following: https://doi.org/10.5281/zenodo.7686273. The simulation output data and the Python scripts used for statistical analysis for the experiments presented in this paper can be accessed at the following permalink, https://osf.io/sbau2/, and referenced with https://doi.org/10.17605/osf.io/sbau2.
Open access
