Agent-based modelling, molluscan population dynamics, and archaeomalacology

Quaternary International xxx (2015) 1e14

Contents lists available at ScienceDirect

Quaternary International journal homepage: www.elsevier.com/locate/quaint

Agent-based modelling, molluscan population dynamics, and archaeomalacology Alex E. Morrison a, b, *, Melinda S. Allen a a b

Anthropology, School of Social Sciences, University of Auckland, Private Bag 92019, Auckland, New Zealand International Archaeological Research Institute, 2081 Young Street, Honolulu, HI 96826, USA

a r t i c l e i n f o

a b s t r a c t

Article history: Available online xxx

Agent-based modelling (ABM) is an emerging archaeological tool that offers insights into processes which are archaeologically invisible or difficult to detect. Here we illustrate the potential of ABM for archaeomalacology, posing two research questions and comparing ABM results with Pacific archaeological sequences. The first analysis considers how molluscan energetic return rates (ERR) and age of reproductive maturity (ARM), singularly or in combination, influence prey population resilience. The second analysis assesses how prey spatial structure affects foraging efficiency and prey susceptibility to resource depression. Consistent with expectations from evolutionary ecology and life history theory, the ABM results demonstrate that both ERR and ARM influence prey population resilience (or vulnerability). However, the analysis also demonstrates that ARM is the more important variable and taxa with high ERR (i.e., large-bodied) are disproportionately affected by human harvesting. Not only are efficient foragers more likely to target high ERR taxa, but these prey often have delayed ARM and un-foraged individuals are more likely to be smaller and immature, with disadvantages for population stability and recovery. In short, early-maturing taxa are highly resilient, while late-maturing organisms are more vulnerable; these outcomes also are observed archaeologically. The ABM analyses also demonstrate the effects of prey spatial structure on molluscan susceptibility to resource depression. High prey aggregation initially allows for high foraging efficiency, but prey abundance and encounter rates often rapidly decline. In contrast, when prey are dispersed, search time is greater, leading to lower encounter rates and reduced foraging efficiency, but greater prey population stability. Our ABM and archaeological examples further illustrate that while general principles can be derived, the resilience and spatial structure of specific prey populations, as well as foraging outcomes, are context dependent and continuously evolving. Finally, we note that model departures from theoretical expectations serve to stimulate further research, including use of additional parameters, consideration of novel contextual evidence, and/or investigation of social, technological or environmental hypotheses. © 2015 Elsevier Ltd and INQUA. All rights reserved.

Keywords: Agent-based modelling Marine mollusc resilience Life history traits Prey spatial structure Foraging dynamics Pacific Islands

1. Introduction Archaeologists analyze plant and animal remains discarded by past human societies with the aims of determining the process or set of processes responsible for creating the observed spatial and temporal patterns. However, in historical sciences (like archaeology, paleontology, and paleoclimatology), the phenomena that comprise the empirical record are the result of multiple processes occurring over different temporal scales, many of which are not

* Corresponding author. Anthropology, School of Social Sciences, University of Auckland, Private Bag 92019, Auckland, New Zealand. E-mail address: [email protected] (A.E. Morrison).

directly observable. Problematically, historical explanations, including those of archaeology, often suffer from equifinality, where a given pattern can be reached from more than one path or set of initial parameters or conditions (Premo, 2010, p. 31; Barton, 2013, p. 155; O'Sullivan and Perry, 2013, p. 52). Of interest here is the character of human molluscan foraging activities in prehistory and in particular the dynamic interactions with molluscan prey populations over time. One way that zooarchaeological studies assess these dynamics is by evaluating patterns of foraging efficiency (energy gained from prey in relation to energy spent in capture). Because search and handling time are archaeologically intractable, archaeological fauna are typically ranked using some proxy of profitability, most often body size or caloric return (e.g., Broughton, 1999; Broughton et al., 2011). The

http://dx.doi.org/10.1016/j.quaint.2015.09.004 1040-6182/© 2015 Elsevier Ltd and INQUA. All rights reserved.

Please cite this article in press as: Morrison, A.E., Allen, M.S., Agent-based modelling, molluscan population dynamics, and archaeomalacology, Quaternary International (2015), http://dx.doi.org/10.1016/j.quaint.2015.09.004

2

A.E. Morrison, M.S. Allen / Quaternary International xxx (2015) 1e14

abundances of high-ranked prey items, relative to less profitable taxa, are used to assess patterns of foraging efficiency over time and explain changes in prey abundances. Human predation often is concluded to be the process responsible for harvesting impacts or resource depression. However, many researchers acknowledge, and sometimes seek to exclude, other factors such as environmental change, technological innovations, and/or non-human predatorprey interactions (e.g., Smith, 1991; Madsen and Schmidt, 1998; Broughton, 1999; Butler, 2001; Wolverton, 2005; Whitaker, 2008a, 2008b, 2010; Faulkner, 2009, 2011; Jerardino, 2014). Many researchers also recognize that resource depression is not a singular process. Charnov et al. (1976) identified three distinct ways that resource depression may arise. Exploitation resource depression occurs when predators adversely impact prey populations directly through hunting, resulting in lowered capture rates. Prey numbers also may be reduced when animals alter their behavior to avoid detection by predators, a phenomenon Charnov et al. (1976) termed behavioral resource depression. This might include flocking behaviors or shifts in the timing of prey activities. Finally, microhabitat resource depression refers to locational shifts in prey to contexts where they are harder to detect or capture, such as offshore movements in fish or elevation shifts in ungulates, processes which may occur without any change in prey population sizes. Distinguishing between these processes can be important to understanding human foraging strategies, prey dynamics, and long-term outcomes, but this is a considerable challenge when working with archaeological assemblages alone. In this paper we use an agent-based model to address two primary research questions regarding long-term dynamics between molluscan prey and human populations: 1) How do differences in life history characteristics affect prey resilience and vulnerability to exploitation resource depression? and 2) What role does prey spatial structure play in vulnerability or resilience to exploitation resource depression? We examine these two questions through the development, application, and analysis of an agentbased simulation model. 1.1. Mollusks and resource depression Some mollusks, especially sessile taxa, may be especially susceptible to population fluctuations as a result of a variety of processes. These include not only human foraging pressures (e.g., Swadling, 1976; Anderson, 1981; Catterall and Poiner, 1987; Mannino and Thomas, 2001, 2002; Braje et al., 2007; Morrison and Hunt, 2007; Faulkner, 2009; Jerardino, 2010; Allen, 2012), but also habitat alterations such as climate-induced coral bleaching
foue €t et al., 2013) (Glynn, 1985, 1993; Przeslawski et al., 2008; Andre or increased water turbidity (Jerardino, 1997; Catterall et al., 2001; Shaw and Richardson, 2001; Morrison and Cochrane, 2008), and variation in prey spatial aggregation (Braje and Erlandson, 2009). Additionally, taxonomic variation in life-history characteristics, such as size at reproductive maturity, reproductive output, and prey body size, may make some taxa more susceptible to exploitation resource depression than others (e.g., Catterall and Poiner, 1987; Poiner and Catterall, 1988; Lasiak, 1991; Whitaker, 2008b; Codding et al., 2014; Whitaker and Byrd, 2014). Roff (2002), for example, notes that the age at which a prey taxon reaches reproductive maturity is negatively correlated with population growth and therefore maturity may be an indicator of the potential resilience of a taxon to exploitation. More specifically, prey that reach reproductive maturity at relatively young ages have more capacity for rapid recovery from intensive harvesting or other traumatic events (Reynolds et al., 2001; Denney et al., 2002; Reynolds, 2003; Hutchings and Reynolds, 2004). Another important life history character that may influence prey resilience is size-at-sexual-

maturity. Lasiak (1991) suggests prey taxa that reach reproductive maturity at sizes below those favored by human collectors are less vulnerable to exploitation resource depression, being less likely to be collected and thus have more opportunities to reproduce and contribute to population maintenance. Prey reproductive output, recruitment success, and intrinsic rate of natural increase also play a significant role in the long-term viability of a species. Pimm (1991, p. 32) notes, “How fast a population recovers from a severe decline in density could depend on its reproductive rate: the more young produced, the faster the population can recover its former level…it seems likely that high reproductive rates will translate into higher population resilience”. Archaeological and ethnographic studies also demonstrate that species with high reproductive rates are relatively impervious to over-harvesting (e.g., Lasiak, 1991; Whitaker and Byrd, 2014). In marine fauna, reproductive output, however, is mediated by cohort recruitment success and juvenile mortality, processes influenced by a variety of environmental factors, including sea surface temperature, substrate characteristics, predation, and competition (Gosselin and Qian, 1997; Hunt and Scheibling, 1997; Black et al., 2011). Problematically these variables are likely to differ from one habitat to the next and their impact can be difficult to evaluate archaeologically. Human ecologists also recognize the role that resource aggregation plays in foraging return rates (e.g., Harpending and Davis, 1977; Kelly, 1983; Cashdan, 1992). While life history traits such as mobility and reproductive strategies can be fundamental to the development of population spatial structures, spatial patterns also can arise from external conditions. For example, habitat change such as coral bleaching can create patchy, segregated prey populations which may be difficult for predators to find. Moreover, mollusk species which aggregate in dense and sessile colonies may be more susceptible to exploitation resource depression relative to those that are solitary, mobile, or dispersed (Catterall and Poiner, 1987; Braje and Erlandson, 2009). Some archaeological studies have explicitly examined the potential role that resource spatial structure plays in foraging return rates and prey susceptibility to over-exploitation (Broughton, 1999; Cannon and Meltzer, 2008; Braje and Erlandson, 2009; Wolverton et al., 2012). Broughton (2002) notes that variation in the spatial distribution of individual prey of different ages and sex leads to intra-species variability in sensitivity to over-exploitation. Specifically, species that form dense, seasonal breeding colonies, such as some waterfowl and pinnipeds, are prone to both exploitation and behavioral resource depression (Broughton, 2002, p. 65; Whitaker, 2010, p. 2569). Here we explore the potential of ABM modelling to evaluate the impact of both life history traits and spatial patterning on forager return rates and prey resilience (or vulnerability) over time. 2. The utility of agent-based models Agent-based modelling is a computational technique oriented towards understanding how population scale patterns arise from the behavior of autonomous agents (individuals or delimited groups) operating at lower scales (Lake, 2015, p. 4). These discrete units or agents vary according to behavioral, life history, or phenotypic traits such as feeding strategy, reproductive rate, age, or size. Agents are generally placed in a model environment which can be either hypothetical or realistic. Agents interact with the environment and with each other through a set of behavioral rules and scheduled actions that are specified at local scales (Railsback and Grimm, 2012, p. 10). Multiple types of agents with different goals and actions can be integrated into a single model. The archaeological record of past mollusk use is a time-averaged phenomenon influenced by a variety of processes operating at

Please cite this article in press as: Morrison, A.E., Allen, M.S., Agent-based modelling, molluscan population dynamics, and archaeomalacology, Quaternary International (2015), http://dx.doi.org/10.1016/j.quaint.2015.09.004

A.E. Morrison, M.S. Allen / Quaternary International xxx (2015) 1e14

sometimes quite different spatial and temporal scales. On onehand, archaeofaunal assemblages reflect prey selected by human foragers, based on both the taxa that are locally available and a set of principles regarding which organisms to pursue and which to ignore. The latter are often economical in focus but mediated by cultural preferences and restrictions. Temporal patterning in archaeofaunal records also may stem from other processes that affect faunal abundances, including climate variability, habitat change, or technological innovations. To further complicate matters, taphonomic processes (e.g., fragmentation, diagenesis, erosion, post-depositional intrusions), sampling, and recovery procedures can hinder our ability to deduce the set of processes that produced empirical patterning in the archaeological record (Schiffer, 1987; Lyman, 1994; Jerardino and Navarro, 2008; Faulkner, 2010). Many archaeological applications use foraging theory models derived from evolutionary ecology (Stephens and Krebs, 1986) to document declining resource return rates, variation in foraging efficiency, and general changes in the use of molluscan taxa through time. Determining the factors that affect foraging efficiency, and building an improved understanding of these processes, can help explain archaeological patterning and generate new archaeological hypotheses. In this respect, an important point is that changing patterns of foraging efficiency are largely determined by prey encounter rates, which can be influenced by varied processes, many which cannot be deduced from the empirical archaeological record alone. To this end, agent-based simulation models are becoming an essential tool for investigating complex ecological interactions between individual organisms and their environments across space and through time. A hallmark of these simulations is that they allow population scale characteristics to be understood as the outcome of interactions between individual organisms, the scale at which natural selection generally operates. These agent-based simulations avoid the problem of teleology which often affects ecosystem modelling, where macro-scale phenomena are treated as the result of macro-scale processes (Richerson, 1977). Here we distinguish between the individual scale, where we model individual prey and predator interactions,

3

and higher level, macro-scale features of our prey populations, which derive from aggregated individual results. While we have incorporated a number of realistic parameter values into the ABM, our approach is largely exploratory (following Premo, 2010; Davies et al., 2015). Exploratory models are built to examine the probability that a given pattern or outcome can be produced under different configurations of parameter settings. As a consequence, the model developed here is purposefully simple and therefore limited in scope to the few relevant variables we are interested in exploring. While the drawback of this approach might be a lack of precision in absolute metrics and predictive power (see O'Sullivan and Perry, 2013), an exploratory approach allows us to determine what different processes among alternatives are more or less likely to produce a given pattern. Furthermore, as a result of this generality, exploratory models are less context-dependent in application and can therefore be used to examine a variety of datasets across multiple settings (e.g., Davies et al., 2015). 3. General model description Our agent-based model was developed using NetLogo 5.1 (Wilensky, 1999) and the results were analyzed using R statistical software (R Core Development Team, 2014). The system we model is a predator-prey environment in which human foragers collect a variety of marine invertebrates; the latter are modelled after shellfish which are common in oceanic island archaeological deposits (details in Tables 1e3). We simulate a hypothetical marine environment consisting of a 101
101 set of regularly spaced square grids, represented as a torus to avoid edge effects (see Fotheringham et al., 2000). The marine environment is akin to a lagoon ecosystem with randomly dispersed coral reef substrate patches, interspersed with sandy bottom substrates. Iterations of the model correspond to one day and foraging trips only occur every three days. The decision to allow foraging once every three days is based on ethnographic accounts from Pacific Island settings which suggest that shellfish were typically one of several food resources and shellfish foraging was not necessarily a daily activity (e.g., Bird and Bliege Bird, 1997; Thomas, 2007; Pinca et al., 2009).

Table 1 State variables and parameter settings used in the agent-based model. State variable Predator Current location Success on last forage Prey captured Number of moves Prey Location No. prey on patch Age Maturity Substrate Reproductive age Reproductive output Energetic return rate (kcal/hr)

Description

Initial parameter values and ranges

Current x, y coordinates corresponding to the location of the agent Boolean variable indicating TRUE or FALSE Number of prey captured by each predator on a single foray Total number of forays an agent has participated in

e FALSE 0 0

X, Y coordinates corresponding to the patch a prey resides within Number of prey residing within a patch Age of the prey items If the prey age is greater than Mk then the value is TRUE Boolean variable indicating either hard substrate or soft substrate Age parameter (Mk) at which prey participate in the reproductive sub-model The number of successful offspring recruited to a hard substrate The caloric return rate of the taxa

n of hard substrate patches 0e2000 0 - 25 years TRUE or FALSE e 1, 2, 3, 4, 5 years 1, 2, 3, 4, 5 250, 500, 750, 1000

Table 2 Simulation symbols, description, and initial model values. Symbol

Description

Initial parameters & ranges

f1 f2 Ak Mk

Movement Strategy 1: area restricted search Movement Strategy 2: correlated random walk Prey age Prey age of reproductive maturity (ARM)

e e 0 to 25 years 1e5 years (continued on next page)

Please cite this article in press as: Morrison, A.E., Allen, M.S., Agent-based modelling, molluscan population dynamics, and archaeomalacology, Quaternary International (2015), http://dx.doi.org/10.1016/j.quaint.2015.09.004

4

A.E. Morrison, M.S. Allen / Quaternary International xxx (2015) 1e14 Table 2 (continued ) Symbol

Description

Initial parameters & ranges

P k Eh d Mr Pw R Rk Moran's I

Prey population starting value Carrying capacity Number of patches with hard substrate Density limit Reproductive output Prey mortality rate Predator population size Prey energetic return rate (kcal/hr) (ERR) Degree of spatial aggregation

15,000e40,000 d* Eh 0 - 10,201 2000/10m2 1e5 0.15 50 250e1000 0.004e0.321

Table 3 Modelled taxon values. Modelled species

ARM

ERR (kcal/hr)

Individual weight (gm)

1 2 3 4

5 4 2 1

1000 750 500 250

1250 16 25 8

The prey population consists of four hard substrate species that vary in terms of energetic return rate (ERR) and age of reproductive maturity (ARM; Table 3). These prey taxa are randomly placed within hard substrate patches. Prey population mortality and reproduction are scheduled once a year and prey age with each iteration of the model. In the model scenarios presented below, individual prey are reproductively viable once they reach a minimum maturity age, after which they can participate in a reproduction sub-model as described below. ARM is a flexible parameter and an important component of one of the research questions we examine in our first case study. To keep foraging pressure uniform we have simplified the model structure by holding human population growth and reproduction constant. We realize that a more realistic modeling exercise would allow human population sizes to rise and fall in accord with the availability of subsistence resources in the corresponding environment. Future models could be adjusted to incorporate Lotka-Volterra predator-prey dynamics (Lotka, 1925; Volterra, 1926). However, marine invertebrates, while important to human diets on Pacific Islands, were only one dietary component in an array of resources that influenced population growth. Accurate modelling of human population growth would add considerable complexity and our current models were designed to be exploratory (see above). 3.1. Mollusk foraging and collection Upon commencement of the model, human foragers participate in two foraging related sub-routines. The “movement” subroutine allows the forager to move according to two possible choices. If the forager was successful at finding prey on their previous foray, they undertake an area restricted search (ARS) (f1) in which the search direction is randomly chosen from a normal distribution, with a possible search radius of 90 in either direction that the agent is facing. If a forager is unsuccessful, they follow a correlated random walk (f2) where the turn angle is randomly chosen from a normal distribution with a range of 2 centered on the direction the agent is currently facing (see Bailleul

et al., 2013 for a similar approach). The foragers move one cell length per time step and check for the presence of molluscan prey items within their corresponding patch. Following foraging theory, the four taxa are ranked according to their energetic return rates (Table 3) and when a predator searches a given patch, the highest ranked items (those most profitable in terms of energetic returns) are taken upon encounter before moving on to collect lower ranked items. Based on ethnographic studies (e.g., Bird and Bliege Bird, 1997; Thomas, 2007), the model allows foragers to collect up to 10 kg of shellfish on a foraging trip. The model uses individual weights of 1250 gm for Species 1, 16 gm for Species 2, 25 gm for Species 3, and 8 gm for Species 4. While we recognize that these mean values may not be representative of every given foraging regime or habitat, they serve the primary purpose of modelling differences in collection numbers based on size variation in the individual taxa collected. We do not model field processing or differential transport of shell back to a central location, although we recognize the potential impact that these handling processes may have on shell patterning in the archaeological record (e.g., Bird and Bliege Bird, 1997).

3.2. Mollusk prey age, reproduction, and mortality In the model, each individual mollusk has a corresponding life history attribute of age (Ak). Mollusks age one day per iteration, in accordance with the general time step of the model. When individuals reach a reproductively mature age (Mk) they are able to reproduce. In our first case study, the reproductively mature age (Mk) parameter is flexible and varies across model runs, from a minimum age of one year to a maximum of five years (see Table 1). Additionally, prey items must be located on a patch with another reproductively mature organism. Reproductive output (RO) (Mr) also is a flexible parameter and varies across model runs. A reproductively viable mollusk produces from one to five offspring, depending on the reproductive output setting, which we vary across model runs to explore how differential reproduction may affect resilience (susceptibility). Successful offspring are considered to have survived the larval stage and become established on a hard substrate. Importantly we do not attempt to model larvae production, dispersal, or early juvenile mortality, due to the difficulties of determining realistic parameter settings for these variables. Notably, invertebrate reproductive output is challenging to measure in natural settings due to the small size of larvae and the rapid rate of early post-settlement mortality, both which make monitoring of these processes problematic (Hunt and Scheibling, 1997, p. 270). Most shellfish species have complex reproductive cycles that include the transport of larvae, through currents, to suitable recruitment sites. The complexity and often large spatial range over which these processes occur also make it

Please cite this article in press as: Morrison, A.E., Allen, M.S., Agent-based modelling, molluscan population dynamics, and archaeomalacology, Quaternary International (2015), http://dx.doi.org/10.1016/j.quaint.2015.09.004

A.E. Morrison, M.S. Allen / Quaternary International xxx (2015) 1e14

difficult to determine the relationship between the number of successfully recruited organisms and the reproductively viable population (Worrapimphong et al., 2010, p. 1337). Rather than attempt to model individual taxonomic differences in reproductive output, given the foregoing uncertainties, we use an exploratory approach and vary reproductive output from one to five individuals per reproductively viable organism across our model runs (see Worrapimphong et al., 2010 for a similar solution). In addition to reproduction, the prey population (P) is subject to a specified adult mortality rate (Pw). The mortality sub-model is run by calculating the prey population at t þ 1 or (P 0 ), which occurs once every 360 iterations.

P 0 ¼ P  ðP*Pw Þ Several researchers have noted that mollusk mortality and survivorship can be affected by a variety of factors, including age, habitat, and taxonomic characteristics (Black et al., 2011; Smith, 2011; Mekawy and Madkour, 2012; Mies et al., 2012). Based on these studies, we fix prey mortality at 15% in all of the models (following Black et al., 2011). We also include a density limit (d) of 2000 prey items per 10 m2 based on formal marine biology surveys and other investigations (e.g.,
foue €t Pradatsundarasar et al., 1989; Schulte et al., 2009; Andre et al., 2013). If the prey population on a patch grows larger than the density limit, the population within the patch is automatically reduced to d-1 by a random mortality procedure. The density-limit sub-model sets the prey population carrying capacity at (k) ¼ d * Eh, where Eh is the number of hard substrate patches available for prey. In our first case study, prey resource taxa are ranked according to post-encounter return rates (Rk), expressed as kcal/hr.

4. Mollusk life history and prey resilience: model results Using a combination of theoretical models drawn from evolutionary ecology, ABM, and archaeological observations, in our first case study we examine how variability in molluscan energetic return rates (ERR), reproductive output (RO), and age of reproductive maturity (ARM) contribute to prey resilience (or vulnerability), in the context of sustained human foraging. Foraging theory provides a set of initial expectations about how human foragers might rank and gather molluscan prey upon encounter, in the main on the basis of energetic return rates or kcal/hr (ERR). Life history theory further highlights the potential impact of ARM and RO, often linked to body size, on prey population stability (Swadling, 1976; Catterall and Poiner, 1987; Poiner and Catterall, 1988; Codding et al., 2014). Particularly relevant to large-bodied shellfish is a tendency for comparatively late reproductive maturity, low reproductive output, and/or long life spans (Peters, 1983; White et al., 2007). In general, these k-strategy life history characteristics are often correlated and therefore exacerbate the effects of human exploitation on high-ranked prey. Further, as noted above, ARM and RO are often negatively correlated with population growth (Roff, 2002, 2003), with entailments for population recovery following adverse impacts. ABM is used to investigate the effects of ERR, ARM, and RO at variable settings across multiple model runs for a set of taxa modelled after shellfish that are common in Pacific Island subsistence economies, now and in the

5

past (see also Morrison and Allen, 2016). We compare the ABM outcomes with Pacific archaeological records, identifying both similarities and departures from theoretical and model expectations. 4.1. Model details In this first case study, human foragers are represented as mobile agents and prey are characterized as immobile properties of patches. Tables 1e3 provide the state variables, their parameters, symbols, and parameter settings, and the ARM values and prey return rates are loosely based on values from economically important Pacific shellfish taxa (Table 3). In the models discussed below, the initial predator population size (R) is 50 and the prey population starting size is set at 10,000 for each of the four species, resulting in a total prey population (P) of 40,000. The model was allowed to run for 10,000 iterations for each combination of parameter settings investigated. Each of these parameter combinations was run 30 times and the reported results are the mean of these 30 runs. Modelled Species 1 reaches reproductive maturity at five years of age (ARM) and has the highest kcal/hr or ERR. As such it is ranked first in terms of profitability and in the ABM Species 1 is always taken upon encounter. Species 2 is ranked second in ERR and reproduces at four years of age. Species 2 individuals are taken only when Species 1 is not present within a given patch, or has been exhausted before reaching the 10 kg transport maximum. Species 3 is ranked third in ERR and has an ARM of two years. Species 3 is only taken in the absence of both Species 1 and 2. Finally, Species 4 is the lowest ranked taxon based on ERR or kcal/hr, has an ARM of one, and is only taken when the other taxa are absent. We also model reproductive output (RO) but, as discussed above, we do not attempt to incorporate realistic reproductive output values for the modelled taxa. 4.2. Results of life history and prey resilience analysis The ABM demonstrates the expected relationship between prey reproductive output (RO), age at reproductive maturity (ARM), and long-term resilience to foraging pressures (Figs. 1e2, Table 4). The modelled taxon with the youngest ARM (Species 4; ARM of one year) produces resilient populations (stable or increasing) when RO is set to parameter values of 2, 3, 4, and 5. Only when RO is set to one year of age does the prey population size diminish due to the combined effects of foraging pressure and low reproductive output. In contrast, the taxon with the highest ERR (Species 1) is always susceptible to resource depression, regardless of RO settings. The susceptibility of the modelled taxa is the result of both high ERRs and low ARMs (see also Section 4.3). Intermediate taxa (Species 2 and 3; ARM of four and two years respectively) are sometimes susceptible and sometimes resilient, depending on the specific values of RO. Fig. 2 is a bivariate sensitivity analysis contour plot based on ARM and RO. The Z value is the terminal prey population size at model completion. The results demonstrate that resilient prey populations are found in taxa that have young ARM and high RO. High RO settings increase prey resilience in all taxa, except that with the greatest ARM (Species 1, Table 3).

Please cite this article in press as: Morrison, A.E., Allen, M.S., Agent-based modelling, molluscan population dynamics, and archaeomalacology, Quaternary International (2015), http://dx.doi.org/10.1016/j.quaint.2015.09.004

6

A.E. Morrison, M.S. Allen / Quaternary International xxx (2015) 1e14

Table 4 The 20 different model parameter combinations and the terminal prey population sizes at the end of each model run. All models began with individual taxon population sizes of 10,000 and a total population of 40,000. Model run

Modelled taxon

Relative reproductive output

Reproductive age

Ending population size

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

4 3 2 1 4 3 2 1 4 3 2 1 4 3 2 1 4 3 2 1

1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 5

1 2 4 5 1 2 4 5 1 2 4 5 1 2 4 5 1 2 4 5

18 10 0 1 161,094 2687 39 16 275,567 40,281 50 58 3,290,198 230,559 1882 187 4,028,447 1,637,606 7102 636

4.3. Relationship between age of reproductive maturity (ARM) and energetic return rate (ERR) Notably, in the modelled taxa presented above there is a positive correlation between ERR and ARM; taxa with high ERR also are late-maturing while taxa with low ERR also are earlymaturing. This makes it difficult to assess how the two variables interact to produce resilient versus susceptible prey populations. To examine this issue, we explored outcomes when the relationships are reversed, that is when ARM and ERR are negatively correlated. In this inverse iteration, we modelled taxa so that those with the greatest ERR had the youngest age of reproductive maturity, and those with the lowest ERR the slowest maturation (i.e., high ARM values). While this is somewhat unrealistic, examining the inverse relationship between ERR and AMR allows the influence of ARM to be explored independent of ERR. Fig. 3 compares the two iterations of the model. The upper panels depict the model run with the positive correlation between ARM and ERR and the lower panel shows the results of the somewhat contrived inverse model. Comparison of these model iterations demonstrates that when ERR is set to the highest value, and the taxon also has low ARM, the prey population is resilient. This demonstrates that while ERR is an important variable in determining if prey taxa will be susceptible to exploitation resources depression, low values of ARM can mediate the effects of heavy predation and foster prey resilience. This would be especially true if low ARM values occurred in a high ERR species, since foragers only move to less profitable prey when higher ranked (i.e., high ERR) taxa are not encountered. Although not common, if a taxon with high ERR and early ARM was present, it would potentially be relatively resilient and effectively shield other more susceptible late-maturing taxa from exploitation resource depression. With respect to prey population collapse (three of the four upper panels in Fig. 3), this occurs when recruitment is no longer able to keep pace with forager harvesting. In our ABM this threshold also is partially structured by the model's specification that reproduction can only take place when there is another reproductively mature individual in close proximity. In real life there are numerous conditions that might contribute to a taxon's decline, especially if already in a weakened state from human

foraging (e.g., environmental change, intensive predation by other species, etc.). As such, the model outcomes presented here are probably “best case” scenarios. Decline is likely to be faster in latematuring taxa (Fig. 3, upper row, far left panel) both because they require more time to reproduce and because the individuals remaining after a foraging bout will most likely be smaller, reproductively immature individuals. The results of this ABM support the findings of Morrison and Allen (2016) that it is the combination of prey rank (based on kcal/hr) and ARM that produces prey susceptibility or resilience. Of note, prey rank is situational and varies according to the taxonomic make-up of a given locality, along with previous patterns of exploitation and other historical factors. In the next section we explore how the ABM results compare with long-term archaeological records from across Polynesia. 4.4. Comparisons with Pacific archaeological records Perhaps not unexpectedly, the ABM results conform to predictions from the prey choice model and life history theory regarding the effects of predation on prey of variable energetic return rates and reproductive ages. Recognizing that a large number of cultural, ecological, and/or environmental factors can affect prey populations in the real world, our ABM effectively serves as a null hypothesis about the potential influences of these key variables. These results can in turn be compared with empirical archaeological records of variability in prey abundances over centennial time frames. To this end, we examine several Polynesian archaeological assemblages from stratified sites where there are taxa similar to our modelled examples (Table 5). Three species that were common food sources in the past are considered: Tridacna, Turbo, and Nerita (Table 6). The expectation is that Tridacna, as a taxon with the highest ERR and ARM (akin to our modeled Species 1), will be the most susceptible to human predation, the small Nerita will be resilient (similar to modeled Species 4), and Turbo will be intermediate (broadly similar to Species 2). Taxa are discussed primarily in terms of rank order abundances, but occasionally in terms of relative abundances (percentages), as for example when rank orders are stable but changing patterns of use are nonetheless suggested.

Please cite this article in press as: Morrison, A.E., Allen, M.S., Agent-based modelling, molluscan population dynamics, and archaeomalacology, Quaternary International (2015), http://dx.doi.org/10.1016/j.quaint.2015.09.004

A.E. Morrison, M.S. Allen / Quaternary International xxx (2015) 1e14

7

Table 5 Details of the archaeological sites that provided the assemblages discussed in the text. Sequence duration refers to the main strata which provided the archaeomalacology assemblages. Sites are ordered by sequence length and approximately from west to east. Island & archipelago

Site name

Island size (km2)

Sequence duration

Molluscan analysis reference

Ofu, American Samoa Tikopia, Solomon Is. Tutuila, American Samoa Aitutaki, Cook Is. Kaua'i, Hawai'i Great Barrier, New Zealand

To'aga (AS-13-1) Sinapupu (Site TK1) Fatu-ma-futi (AS-25-055) Moturakau Shelter (MR1) Nualolo Kai (Site K3) Harataonga Beach (multiple areas)

7 5 142 18 1430 285

c. 2200 yrs 1800 yrs 1500 yrs 850 yrs 500 yrs 450 yrs

Nagaoka, 1993 Kirch and Yen, 1982 Morrison and Addison, 2008 Allen, 1992 Morrison and Hunt, 2007 Allen, 2012

Table 6 Energetic return rate (ERR) and age of reproductive maturity (ARM) for common archaeological taxa (from Codding et al., 2014 and sources therein). Taxon

Energetic return (kcal/hr)

Age of reproductive maturity

Nerita Turbo Tridacna

42 to 1106 520 to 606 2622 to 13,064

1e2 years 3e4 years 4e5 years

On the small islands of Tikopia and Aitutaki, Tridacna maxima clams generally are archaeologically well represented, despite ethnographic studies showing that the heavy shells are often discarded at sites of collection (e.g., Bird and Bliege Bird, 1997; Thomas, 2007). At the Sinapupu site (Tikopia Island) Tridacna is the 2nd or 1st ranked taxon at several points in the sequence but eventually slips to below 5th and in the last few centuries of occupation further declines in abundance (Kirch and Yen, 1982, Table 45, Figure 117). At Moturakau Rockshelter (Aitutaki Island), Tridacna ranks first in the oldest occupation layers but decreases in abundance over time to rank 2 or 3, (Allen, 1992, Appendix E, Tables E.4 and E.5). In both cases, Turbo becomes increasingly important while Tridacna decreases, consistent with our ABM model predictions that comparatively smaller-bodied and earliermaturing taxa will be more resilient. Also consistent with the ABM outcomes, Turbo, a taxon of intermediate body size and intermediate reproductive age, sustains foraging pressures for many centuries, with no statistically significant size changes observed in the Moturakau case (Allen, 1992, Figure 8.11). At To‘aga (Ofu Island), Turbo setosus is the highest ranked taxon by weight. While its rank is stable over time, the proportional contribution of Turbo increases across the approximately 2200-year sequence, as represented in the “main trench” (Kirch and Hunt, 1993; Nagaoka, 1993). Tridacna is the second highest ranked taxon on Ofu across nearly all layers but its relative abundance more than halves over time. At a second Samoan site, Fatu-ma-futi (Tutuila Island), there is more stability, especially if the quite small sample relating to initial use of the newly formed beach (Layer IV) is excluded (see Morrison and Addison, 2008). Tridacna maxima is a top ranked taxon throughout the archaeological sequence (varying between ranks 1 and 3), but declines in frequency over time, from 20 to 22% (Layers II and III) to 15% (Layer I) (Morrison and Addison, 2008, Table 1). Turbo (ranking between 1st and 2nd), fluctuates between 16 and 19% of the assemblage in any given layer. The small Nerita, in contrast, more than doubles over time, increasing from 5% (Layer III) to 27% in the most recent layer (Layer I). In the Hawaiian Archipelago, at Nualolo Kai (Kaua‘i Island) where Tridacna does not occur, the endemic Turbo sandwicensis was one of the largest species available to prehistoric foragers, reaching lengths of 90 mm (in Morrison and Hunt, 2007, 333). Here Turbo is the top-ranking taxon across all three analytical zones but it declines in both absolute and relative abundance over time, from 64% in the early Zone C to 19% in the more recent Zone A. The much smaller Nerita picea does not figure prominently in the overall

shellfish assemblage but as a component of the shoreline habitat group it increases from 37 to 45% over time. Moving south to subtropical New Zealand, at Harataonga (Great Barrier Island) the largest native turban (Cookia sulcata) was common in the earliest occupation layer varying from 1st to 2nd rank (39e21%), but declines to around rank 4 to 5 (5%) in later occupation layers, where it is eclipsed by the smaller Turbo smaragdus (Allen, 2012). These shifts are accompanied by moderately significant decreases in Turbo opercula sizes across the two main time periods. At the same time use of the very small Nerita atramentosa increases dramatically in late prehistory, rising to rank 2. The foregoing examples highlight the relative vulnerability of the large Tridacna clams to human predation and the process of prey switching to lower ERR and more resilient genera as vulnerable taxa decline. Overall, our results suggest that ERR and ARM are major drivers of foraging-induced resource depression. Further, because ARM tends to be correlated with body size, large taxa are not only disproportionately targeted by foragers but also more disadvantaged with respect to population recovery relative to young ARM taxa. 5. Prey spatial structure and foraging efficiency: model results In a second case study we use ABM to consider the specific processes that link differences in prey spatial structure with archaeological measures of foraging efficiency, prey population size, and exploitation resource depression. Aggregation is defined using Moran's I test for spatial autocorrelation (Moran, 1950) and in this analysis is a measure of the spatial proximity of the coral reef patches which are available for mollusk settlement. As outlined above, prey aggregation may stem from life history characters, such as colonial habits or larval dispersal ranges, while prey dispersion could be the result of not only exploitation by human predators but also habitat fragmentation, disease, or other factors. We focus here on the outcomes rather than specific causes, examining how spatial structure (aggregated or dispersed) and sustained harvesting interact to affect prey population sizes and prey persistence, along with human foraging efficiency over time. 5.1. Model details The location of patches available for mollusks is controlled by an algorithm that runs at initialization of the model. During the model setup procedure, a set of patches is used to generate the degree of spatial aggregation of prey resources by varying the radius (r) within which further patches are allowed to host prey (following Bailleul et al., 2013). An error parameter (e) also is used to impart a degree of stochasticity into the model. The overall result is that as r increases, so does the degree of aggregation. Spatial autocorrelation occurs when values in one patch trait influence corresponding trait values in nearby patches (Fotheringham et al., 2000, p. 12). Moran's I is used here to quantify

Please cite this article in press as: Morrison, A.E., Allen, M.S., Agent-based modelling, molluscan population dynamics, and archaeomalacology, Quaternary International (2015), http://dx.doi.org/10.1016/j.quaint.2015.09.004

8

A.E. Morrison, M.S. Allen / Quaternary International xxx (2015) 1e14

Fig. 1. Results of the ABM model runs (means of 30 runs for each parameter setting). Note that the Y-axis values for each individual model pane refer to prey population size and the X-axis values refer to the time step. All model individual taxon populations began at 10,000 and the ending population sizes are presented in Table 4.

the degree of spatial autocorrelation in prey resource abundance per patch, across all of the patches in the environment. This allows for more accurate quantification of how habitat structure and prey aggregation affect prey population size and forager return rates. We quantify foraging efficiency as the number of items captured per foraging trip at the population scale. The start environment for the different aggregation parameter values and model runs is shown in Fig. 4. The model used in this second case study follows the same basic structure as that used in case study 1, with some important differences. Rather than include four species, for simplicity forgers focus on only one taxon, our hypothetical Species 1. We also do not vary reproductive output (RO) and keep the parameter setting at a value of 2. In this model a total of 15,000 prey are randomly

distributed across the available hard substrate patches upon initialization. We set the error parameter so that 10% of the prey population will settle outside of the defined aggregation areas. 5.2. Results of prey spatial structure analysis The spatial analysis informs on two aspects of predator-prey interaction at five levels of aggregation as defined by Moran's I. First, looking at the relationship from the perspective of prey, the models show how prey population size is affected over time at different levels of spatial aggregation. When spatial aggregation is high (Moran's I ¼ 0.321 and 0.279), foraged prey populations decline over time and eventually become extinct (Fig. 5, red and blue lines (in the web version)). In contrast, foraged populations

Please cite this article in press as: Morrison, A.E., Allen, M.S., Agent-based modelling, molluscan population dynamics, and archaeomalacology, Quaternary International (2015), http://dx.doi.org/10.1016/j.quaint.2015.09.004

A.E. Morrison, M.S. Allen / Quaternary International xxx (2015) 1e14

9

are cyclical episodes of coral reef bleaching and die-off, particularly in response to above average sea surface temperatures, or periods
foue €t et al., 2013; of increased terrestrial erosion and run-off (Andre Bahr et al., 2015), processes which adversely affect the availability and quality of suitable substrates. More robust explanations of long-term prehistoric fluctuations in Tridacna, and in particular conclusive demonstration that ancient foragers were responsible, will ultimately require evaluation of the contemporaneous paleoclimate conditions, and local reef histories. 6. Discussion 6.1. Foraging and life history analysis

Fig. 2. Bivariate sensitivity analysis contour plot comparing age at reproductive maturity with reproductive output. The Z value is the population size at completion of the model. The values are the means of 30 runs for each parameter setting.

that are not aggregated (Moran's I ¼ 0.004 and 0.097) are stable or increase (Fig. 5, yellow and green lines (in the web version)). The analysis also informs on how prey spatial structure affects foraging efficiency (Fig. 6). Aggregated prey (Moran's I ¼ 0.321, and 0.279) initially allow for very high foraging efficiency (above 0.8), but resource depression occurs quite rapidly (red and blue lines (in the web version)), trending below outcomes for prey that are more dispersed. Alternatively, when aggregation is low (Moran's I ¼ 0.004 and 0.097), foraging efficiency is reduced or dampened, but foraging returns are relatively stable over time. The main result from this model is the demonstration that spatially aggregated prey initially allow for high foraging efficiency, essentially because concentrated resources reduce search and travel time. However, the benefits may be short-lived, with exploitation resource depression being a likely outcome. These results also can be compared with the foregoing archaeological sequences. Along with providing high ERR and having a late ARM, Tridacna maxima populations are further disadvantaged with respect to human foragers by their proclivity for highly clustered settlement patterns. Turbo, in contrast, are both mobile and less colonial, characteristics which provide a buffer against human foragers. The results of the spatial structure analysis also raise questions about the recurring regional patterns of Tridacna decline. While there are several reasons to expect that the observed archaeological patterning discussed above reflects exploitative resource depression, the ABM spatial analysis highlights how other processes which fragment or destroy marine habitats or substrates also can contribute to declines in mollusk populations over time. This is especially true of sessile taxa like Tridacna, as well as others like mussels, abalone, oysters and chitons which vary in ARM. Among the other potentially important processes of disturbance

Our agent-based model (ABM) successfully reproduced outcomes predicted by life history theory with respect to interactions between foragers, prey reproductive output (RO), and prey age of reproductive maturity (ARM). In general, high ERR prey are more vulnerable to exploitation resource depression because they are more frequently taken upon encounter. Similarly, molluscan taxa that reach reproductive maturity comparatively late in life, a common feature in large-bodied prey, are more likely to be sensitive to human foragers because they are more likely to be harvested (having higher ERR), and therefore have reduced opportunities to reproduce and contribute to population maintenance. High reproductive output (RO) increases the resilience of prey populations under predation pressure; nonetheless, our data show that taxa with both high ERR and high ARM are still susceptible to exploitation resource depression regardless of the reproductive output values modeled here. One unexpected finding of the current ABM is that taxa with high ERR that also are early-maturing (low ARM) would inadvertently contribute to stability in the other high ERR prey taxa with comparatively higher ARM values by deflecting predation pressures. Finally, we found that archaeomalacology sequences from across the Pacific region, derived from a variety of marine settings, variable time sequences, and presumably accumulated under different harvesting regimes, have resonance with our ABM results. With these findings, we are now in a position to assess the role of additional parameters, beyond ERR, ARM and RO through our ABM platform. For example, Catterall and Poiner (1987) predict that prey visibility and ease of collection are key variables affecting prey resilience and their potential for rebound following resource depression. Visibility can be affected by prey habits, such as burrowing in soft sediments or exposure during reproductive periods. Additionally, other life history traits that typify large-bodied taxa across many vertebrate groups should be evaluated for these and other economically important molluscan taxa. These include variation in growth rates, larval production and survival, and natural population densities. Predation effects related to other predators also need to be considered, especially in contexts where trophic relations have been altered by human foraging (e.g., Rick and Erlandson, 2009). The ABM also helps to identify missing contextual information and new research questions. For example, are the specific life history traits of the endemic Hawaiian Turbo a factor in its rapid decline at the Hawaiian site of Nualolo Kai? Why was Tridacna resilient on Tikopia Island for centuries but vulnerable on Aitutaki Island where there is a shorter cultural sequence but also a large, coral studded lagoonda high quality habitat for clams? One possibility is variation in Tridacna management strategies and mollusk harvesting techniques (e.g., Whitaker, 2008b). Differences in shell collection and discard also could play a role in structuring the archaeological records of these two localities. Alternatively, the size of Tikopia, coupled with cultural controls (see Kirch and Yen, 1982),

Please cite this article in press as: Morrison, A.E., Allen, M.S., Agent-based modelling, molluscan population dynamics, and archaeomalacology, Quaternary International (2015), http://dx.doi.org/10.1016/j.quaint.2015.09.004

10

A.E. Morrison, M.S. Allen / Quaternary International xxx (2015) 1e14

Fig. 3. Comparison of ABM runs where the relationship between ARM and ERR is positive (upper panels) and negative (lower panels). Reproductive Output is set to 2. All individual taxon prey populations began at 10,000. The values are the means of 30 runs for each parameter setting.

may have limited human population growth on this diminutive island. 6.2. Prey aggregation and resilience versus susceptibility Our spatial ABM model demonstrates that prey aggregation initially increases foraging efficiency but then rapidly leads to exploitation resource depression, declining foraging efficiency, and ultimately prey extinction in populations with highly aggregated spatial structures. These results corroborate a variety of empirical and theoretical models from behavioral ecology. For example, in a study of predator-prey dynamics involving prairie skinks that foraged on red-legged grasshoppers, Pitt and Ritchie (2002, p.161) found that as prey aggregation increased, so did predator consumption rates. They conclude that “predators could consume more aggregated prey than dispersed prey and increased consumption rate may result in predator versus resource limited prey populations dynamics...” (p. 162). ABM simulations on other types of foragers also are consistent with our ABM outcomes. Barraquand

et al. (2009, p 511), for example, found that when predators could not remember the locations where they had previously taken prey, dispersion of prey significantly limited the rate at which they were captured, while aggregation increased prey detection and consumption. ABM in combination with field analyses of beluga whale migratory patterns also is instructive (Bailleul et al., 2013). Prey aggregation affected both food intake and the timing of whale migrations, with the whales being more likely to prolong their stay during summer months when prey were moderately aggregated (Bailleul et al., 2013, pp. 2541-2). These studies, combined with our results, illustrate how prey aggregation can not only reduce predator search time but also affect prey consumption patterns. Environmental processes also can create dispersed distributions with implications for human foraging efficiency. In some cases the number of prey encountered may decrease relative to search time, even though the overall number of prey in the population remains stable. Our ABM prey outcomes could arise from any process that fragments large, continuous areas of marine habitat, whether climate change, disease, or human harvesting activities. In contrast,

Please cite this article in press as: Morrison, A.E., Allen, M.S., Agent-based modelling, molluscan population dynamics, and archaeomalacology, Quaternary International (2015), http://dx.doi.org/10.1016/j.quaint.2015.09.004

A.E. Morrison, M.S. Allen / Quaternary International xxx (2015) 1e14

11

Fig. 4. Start world for the agent-based model runs used in case study 2. Spatial aggregation increases from A) to E). Moran's I values ¼ A) 0.004, B) 0.097, C) 0.199, D) 0.279, E) 0.321.

prey with naturally aggregated spatial structures may initially provide for high encounter rates and high foraging efficiency. These situations can arise from a species' natural patterns, specific reproductive behaviors (Broughton, 2002), or gregariousness as a defensive mechanism (Wolverton et al., 2012). Increased encounter rates and foraging efficiency are especially anticipated when prey are stationary and their locations predictable, as is the case for many species of marine mollusks. When prey are aggregated and immobile, forager strategies that facilitate recall and communication of resource locations would be important adaptations, enhancing resource exploitation and adversely affecting prey population resilience (Morrison and Allen, 2014). Finally, it is worth noting that naturally aggregated species can vary considerably with

respect to ARM, from late-maturing genera like Tridacna to earlymaturing taxa like mussels or barnacles (Suchanek, 1981); this variation in ARM has the potential to differentially offset the generally adverse effects of aggregation. 6.3. Conclusions Agent-based modelling is a rapidly developing computational technique which is increasingly being used to gain insights into processes that might be archaeologically invisible or difficult to detect. Among the advantages are an explicitly defined model environment, well controlled state variables, and the possibility of repeated testing of both environmental and agent processes under

Please cite this article in press as: Morrison, A.E., Allen, M.S., Agent-based modelling, molluscan population dynamics, and archaeomalacology, Quaternary International (2015), http://dx.doi.org/10.1016/j.quaint.2015.09.004

12

A.E. Morrison, M.S. Allen / Quaternary International xxx (2015) 1e14

ABM is now firmly established within a range of ecological, economic, and social sciences. We anticipate that ABM will increasingly become a mainstream tool for the archaeological community as well, stimulating investigation into a wide range of both cultural and ecological processes (e.g., Axtell and Epstein, 1996; Grimm and Railsback, 2005; Kohler and Varien, 2012; Wren et al., 2014; Davies et al., 2015; Wurzer et al., 2015), and with results that will be of value to a diversity of disciplines. We suggest these simulations models are especially useful for zooarchaeological studies where the evolutionary outcomes of interacting agents, who have the abilities to make both patterned and random decisions, are a focus. Fig. 5. Effects of prey spatial structure on prey population size: values are the means of 30 runs for each parameter setting.

Acknowledgments This research, and Morrison's position as a Post-doctoral Research Fellow at University of Auckland, was supported by Marsden Fund, The Royal Society of New Zealand, Grant 11-UOA027 to M. S. Allen, A. M. Lorrey, and M. N. Evans. Ben Davies provided valuable advice on programming using lists in NetLogo and Beau DiNapoli commented on an early version of the ABM which helped to improve the final product. We thank two anonymous reviewers for comments that greatly improved the manuscript. Briar Sefton assisted with final production of the graphics. Zac McIvor's assistance with the bibliography is much appreciated. Finally, we thank the special issue editors Antonieta Jerardino, Patrick Faulkner, and Carola Flores for their invitation to be part of this special issue. References

Fig. 6. Effects of prey spatial structure on foraging efficiency: values are the means of 30 runs for each parameter setting.

variable settings. These features of ABM allow complex ecological processes, like predator-prey dynamics, to be examined simply and iteratively, one variable at a time, or more realistically through complex multivariate simulations, albeit with less control over the interaction of specific parameters. ABMs also offer the advantage of providing long-term, decadal, centennial, or even millennial scale outcomes which can be explored for emergent phenomena through multiple model runs. The two case studies examined here consider specifically how molluscan prey resilience and human foraging efficiency are affected by prey energetic returns, life history traits, and spatial patterning. The model results in many respects replicate predictions from evolutionary ecology and life history theory, and empirical archaeological sequences derived from varied marine environments, geographic and cultural contexts, and centennial to millennial time periods. We identify prey with high energetic returns, late reproductive ages, and spatially-clustered settlement structures, like Tridacna, as potentially highly vulnerable to human foragers. In contrast small-bodied, early-maturing prey, like Nerita, are often resilient to intensive foraging. Prey of intermediate ERR and ARM, such as some species of Turbo, represent forager tradeoffs between energetic return rates and sustainable food sources. Usefully, empirical archaeological departures from ABM outcomes are potentially as interesting as results which are congruent, often suggesting new avenues of research. These can include investigation of parameters not considered in the original simulation, as well as exploration of additional parameter values. In some cases model results and archaeological comparisons may point to the need for additional contextual evidence with which to better assess specific archaeological outcomes or, alternatively, testing of other social, technological or environmental hypotheses.

Allen, M.S., 1992. Dynamic Landscapes and Human Subsistence: Archaeological Investigations on Aitutaki Island, Southern Cook Islands. University of Washington. University Microfilms International, Ann Arbor (Ph.D. thesis). Allen, M.S., 2012. Molluscan foraging efficiency and mobility amongst agricultural foragers: a case study from northern New Zealand. Journal of Archaeological Science 39, 295e307. Anderson, A.J., 1981. A model of prehistoric collecting on the rocky shore. Journal of Archaeological Science 8, 109e120.
foue €t, S., Van Wynsberge, S., Gaertner-Mazouni, N., Menkes, C., Gilbert, A., Andre Remoissenet, G., 2013. Climate variability and massive mortalities challenge giant clam conservation and management efforts in French Polynesia atolls. Biological Conservation 160, 190e199. Axtell, R., Epstein, J.M., 1996. Growing Artificial Societies: Social Science from the Bottom Up. The MIT Press and Brookings Institution, Cambridge. Bahr, K.D., Jokiel, P.L., Rodgers, K.S., 2015. The 2014 coral bleaching and freshwater flood events in K aneʻohe Bay, Hawaiʻi. PeerJ 3, e1136. Barraquand, F., Inchausti, P., Bretagnolle, V., 2009. Cognitive abilities of a central place forager interact with prey spatial aggregation in their effect on intake rate. Animal Behaviour 78, 505e514. Barton, C.M., 2013. Stories of the past or science of the future? Archaeology and computational social science. In: Bevan, A., Lake, M. (Eds.), Computational Approaches to Archaeological Spaces. Left Coast Press, Walnut Creek, CA, pp. 151e178. Bailleul, F., Grimm, V., Chion, C., Hammill, M., 2013. Modeling implications of food resource aggregation on animal migration phenology. Ecology and Evolution 3, 2535e2546. Bird, D.W., Bliege Bird, R.L., 1997. Contemporary shellfish gathering strategies among the Meriam of the Torres Strait Islands, Australia: testing predictions of a central place foraging model. Journal of Archaeological Science 24, 39e63. Black, R., Johnson, M.S., Prince, J., Brearley, A., Bond, T., 2011. Evidence of large, local variations in recruitment and mortality in the small giant clam, Tridacna maxima, at Ningaloo Marine Park, Western Australia. Marine and Freshwater Research 62, 1318e1326. Braje, T.J., Erlandson, J., 2009. Mollusks and mass harvesting in the Middle Holocene. California Archaeology 1, 269e289. Braje, T.J., Kennett, D.J., Erlandson, J.M., Culleton, B.J., 2007. Human impacts on nearshore shellfish taxa: a 7,000 year record from Santa Rosa Island, California. American Antiquity 72, 735e756. Broughton, J.M., 1999. Resource Depression and Intensification during the Late Holocene, San Francisco Bay: Evidence from the Emeryville Shellmound Vertebrate Fauna. University of California Press, Berkeley. Broughton, J.M., 2002. Prey spatial structure and behavior affect archaeological tests of optimal foraging models: examples from the Emeryville Shellmound vertebrate fauna. World Archaeology 34, 60e83.

Please cite this article in press as: Morrison, A.E., Allen, M.S., Agent-based modelling, molluscan population dynamics, and archaeomalacology, Quaternary International (2015), http://dx.doi.org/10.1016/j.quaint.2015.09.004

A.E. Morrison, M.S. Allen / Quaternary International xxx (2015) 1e14 Broughton, J.M., Cannon, M.D., Bayham, F.E., Byers, D.A., 2011. Prey body size and ranking in zooarchaeology: theory, empirical evidence, and applications from the northern Great Basin. American Antiquity 76, 403e428. Butler, V.L., 2001. Changing fish use on Mangaia, southern Cook Islands: resource depression and the prey choice model. International Journal of Osteoarchaeology 11, 88e100. Cannon, M.D., Meltzer, D.J., 2008. Explaining variability in early Paleoindian foraging. Quaternary International 191, 5e17. Cashdan, E.A., 1992. Spatial organization and habitat use. In: Smith, E.A., Winterhalder, B. (Eds.), Evolutionary Ecology and Human Behavior: Foundations of Human Behavior. Aldine de Guyter, New York, pp. 237e266. Catterall, C.P., Poiner, I.R., 1987. The potential impact of human gathering on shellfish populations, with reference to some NE Australian intertidal flats. Oikos 114e122. Catterall, C.P., Poiner, I.R., O'Brien, C.J., 2001. Long-term population dynamics of a coral reef gastropod and responses to disturbance. Austral Ecology 26, 604e617. Charnov, E.L., Orians, G.H., Hyatt, K., 1976. Ecological implications of resource depression. American Naturalist 19, 247e259. Codding, B.F., O'Connell, J.F., Bird, D.W., 2014. Shellfishing and the colonization of Sahul: a multivariate model evaluating the dynamic effects of prey utility, transport considerations and life-history on foraging patterns and midden composition. The Journal of Island and Coastal Archaeology 9, 238e252. Davies, B., Holdaway, S., Fanning, P.J., 2015. Modeling the Palimpsest: an Exploratory Agent-based Model of Surface Archaeological Deposit Formation in a Fluvial Arid Australian Landscape. The Holocene (in press). Denney, N.H., Jennings, S., Reynolds, J.D., 2002. Life-history correlates of maximum population growth rates in marine fishes. Proceedings of the Royal Society of London. Series B: Biological Sciences 269, 2229e2237. Faulkner, P., 2009. Focused, intense and long-term: evidence for granular ark (Anadara granosa) exploitation from late Holocene shell mounds of Blue Mud Bay, northern Australia. Journal of Archaeological Science 36, 821e834. Faulkner, P., 2010. Morphometric and taphonomic analysis of granular ark (Anadara granosa) dominated shell deposits of Blue Mud Bay, northern Australia. Journal of Archaeological Science 37, 1942e1952. Faulkner, P., 2011. Late Holocene mollusc exploitation and changing near-shore environments: a case study from the coastal margin of Blue Mud Bay, northern Australia. Environmental Archaeology 16, 137e150. Fotheringham, A.S., Brunsdon, C., Charlton, M., 2000. Quantitative Geography: Perspectives on Spatial Data Analysis. Sage, Thousand Oaks, CA. ~ o-associated disturbance to coral reefs and post disturGlynn, P.W., 1985. El Nin bance mortality by Acanthaster planci. Marine Ecology Progress Series 26, 295e300. Glynn, P.W., 1993. Coral reef bleaching: ecological perspectives. Coral Reefs 12, 1e17. Gosselin, L.A., Qian, P.Y., 1997. Juvenile mortality in benthic marine invertebrates. Marine Ecology Progress Series 146 (1), 265e282. Grimm, V., Railsback, S.F., 2005. Individual-based Modeling and Ecology. Princeton University Press, Princeton. Harpending, H., Davis, H., 1977. Some implications for hunter-gatherer ecology derived from the spatial structure of resources. World Archaeology 8, 275e286. Hunt, H.L., Scheibling, R.E., 1997. Role of early post-settlement mortality in recruitment of benthic marine invertebrates. Marine Ecology Progress Series 155, 269e301. Hutchings, J.A., Reynolds, J.D., 2004. Marine fish population collapses: consequences for recovery and extinction risk. BioScience 54, 297e309. Jerardino, A., 1997. Changes in shellfish species composition and mean shell size from a late-Holocene record of the west coast of southern Africa. Journal of Archaeological Science 24, 1031e1044. Jerardino, A., 2010. Large shell middens in Lamberts Bay, South Africa: a case of hunteregatherer resource intensification. Journal of Archaeological Science 37 (9), 2291e2302. Jerardino, A., 2014. Variability in late Holocene shellfish assemblages: the significance of large shore barnacles (Austromegabalanus cylindricus) in South African west coast sites. Journal of Archaeological Science 52, 56e63. Jerardino, A., Navarro, R., 2008. Shell morphometry of seven limpet species from coastal shell middens in southern Africa. Journal of Archaeological Science 35, 1023e1029. Kelly, R.L., 1983. Hunter-gatherer mobility strategies. Journal of Anthropological Research 39, 277e306. Kirch, P.V., Yen, D.E., 1982. Tikopia. The Prehistory and Ecology of a Polynesian Outlier. Bishop Museum Press, Honolulu. Kirch, P.V., Hunt, T.L., 1993. The To'aga Site: Three Millennia of Polynesian Occupation in the Manu'a Islands, American Samoa. University of California, Berkley. Kohler, T.A., Varien, M.D., 2012. Emergence and Collapse of Early Villages: Models of Central Mesa Verde Archaeology, vol. 6. University of California Press. Lake, M.W., 2015. Explaining the past with ABM: on modelling philosophy. In: Wurzer, G., Kowerik, K., Reschreiter, H. (Eds.), Agent-based Modeling and Simulation in Archaeology. Springer International Publishing, pp. 3e35. Lasiak, T., 1991. The susceptibility and/or resilience of rocky littoral molluscs to stock depletion by the indigenous coastal people of Transkei, southern Africa. Biological Conservation 56, 245e264. Lotka, A.J., 1925. Elements of Physical Biology. Williams and Wilkins, Baltimore. Lyman, R.L., 1994. Vertebrate Taphonomy. Cambridge University Press, Cambridge. Madsen, D.B., Schmidt, D.N., 1998. Mass collecting and the diet breadth model: a Great Basin example. Journal of Archaeological Science 25, 445e455.

13

Mannino, M.A., Thomas, K.D., 2001. Intensive Mesolithic exploitation of coastal resources? Evidence from a shell deposit on the Isle of Portland (Southern England) for the impact of human foraging on populations of intertidal rocky shore molluscs. Journal of Archaeological Science 28, 1101e1114. Mannino, M.A., Thomas, K.D., 2002. Depletion of a resource? The impact of prehistoric human foraging on intertidal mollusc communities and its significance for human settlement, mobility and dispersal. World Archaeology 33, 452e474. Mekawy, M.S., Madkour, H.A., 2012. Studies on the Indo-Pacific Tridacnidae (Tridacna maxima) from the Northern Red Sea, Egypt. International Journal of Geosciences 3, 1089e1095. Mies, M., Braga, F., Scozzafave, M.S., Lavanholi de Lemos, E., Sumida, P.Y.G., 2012. Early development, survival and growth rates of the giant clam Tridacna crocea (Bivalvia: Tridacnidae). Brazilian Journal of Oceanography 60, 127e133. Moran, P.A., 1950. Notes on continuous stochastic phenomena. Biometrika 17e23. Morrison, A.E., Addison, D.J., 2008. Assessing the role of climate change and human predation on marine resources at the Fatu-ma-Futi site, Tutuila Island, American Samoa: an agent based model. Archaeology in Oceania 43, 22e34. Morrison, A.E., Allen, M.S., 2014. Investigating the Effects of Climate Variability on Long-term Predator/prey Dynamics in Marine Ecosystems. An Agent-based Model. In: Paper Presented at the 12th International Conference of Archaeozoology, September 22 e 27, 2014, San Rafael, Argentina. Morrison, A.E., Allen, M.S., 2016. Modelling prey population resilience in ancient
te
marine ecosystems: an agent-based model application. Bulletin de la Socie
historique française (in press). pre Morrison, A.E., Cochrane, E.E., 2008. Investigating shellfish deposition and landscape history at the Natia Beach site, Fiji. Journal of Archaeological Science 35, 2387e2399. Morrison, A.E., Hunt, T.L., 2007. Human impacts on the nearshore environment: an archaeological case study from Kaua'i, Hawaiian Islands. Pacific Science 61, 325e345. Nagaoka, L., 1993. Faunal assemblages from the To'aga site. In: Kirch, P.V., Hunt, T.L. (Eds.), The To‘aga Site: Three Millennia of Polynesian Occupation in the Manu‘a Islands, American Samoa. Contributions of the University of California Archaeological Research Facility, Berkeley, pp. 189e216. O'Sullivan, D., Perry, G.L.W., 2013. Spatial Simulation: Exploring Pattern and Process. Wiley-Blackwell, Chichester. Peters, R.H., 1983. The Ecological Implications of Body Size. Cambridge University Press, Cambridge. Pimm, S.L., 1991. The Balance of Nature?: Ecological Issues in the Conservation of Species and Communities. University of Chicago Press. Pinca, S., Awira, R., Kronen, M., Chapman, L., Lasi, F., Pakoa, K., Boblin, P., Friedman, K., Magron, F., Tardy, E., 2009. Cook Islands Country Report: Profiles and Results from Survey Work at Aitutaki, Palmerston, Mangaia and Rarotonga (February and October 2007). Secretariat of the Pacific Community, Coastal Fisheries Program, Noumea, New Caledonia. Pitt, W.C., Ritchie, M.E., 2002. Influence of prey distribution on the functional responses of lizards. Oikos 96, 157e163. Poiner, I.R., Catterall, C.P., 1988. The effects of traditional gathering on populations of the marine gastropod Strombus luhuanus Linne 1758, in southern Papua New Guinea. Oecologia 76, 191e199. Pradatsundarasar, A.O., Saichuae, P., Teerarkup, K., Gajaseni, N., 1989. Changing of Chemical Conditions in Mae Klong River and Ecosystem of Mae Klong River Mouth Around Don Hoi Lord, Samut Songkhram Province (Phase 2). Department of Biology. Chulalongkorn University, Bangkok. Premo, L.S., 2010. On the role of agent-based modeling in post-positivist archaeology. In: Costopoulos, A., Lake, M. (Eds.), Simulating Change: Archaeology into the 21st Century. Utah University Press, Salt Lake City, pp. 28e37. Przeslawski, R., Ahyong, S., Byrne, M., Woerheide, G., Hutchings, P.A.T., 2008. Beyond corals and fish: the effects of climate change on non-coral benthic invertebrates of tropical reefs. Global Change Biology 14, 2773e2795. Railsback, S.F., Grimm, V., 2012. Agent-based and Individual-based Modeling: a Practical Introduction. Princeton University Press, Princeton. R Core Development Team, 2014. R: a Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. http:// www.R-project.org. Reynolds, J.D., 2003. Life histories and extinction risk. In: Blackburn, T.M., Gaston, K.J. (Eds.), Macroecology: Concepts and Consequences. Blackwell, Oxford, pp. 195e217. Reynolds, J.D., Jennings, S., Dulvy, N.K., 2001. Life histories of fishes and population responses to exploitation. In: Reynolds, J.D., Mace, G.M., Redford, K.H., Robinson, J.G. (Eds.), Conservation of Exploited Species. Cambridge University Press, Cambridge, pp. 147e168. Richerson, P.J., 1977. Ecology and human ecology: a comparison of theories in the biological and social sciences. American Ethnologist 4, 1e26. Rick, T.C., Erlandson, J.M., 2009. Coastal exploitation. Science 325, 952e953. Roff, D.A., 2002. Life History Evolution. Sinauer Associates, Sunderland. Roff, D.A., 2003. Foraging Theory. Princeton University Press, Princeton. Schiffer, M.B., 1987. Formation Processes of the Archaeological Record. University of Utah Press, Salt Lake City. Shaw, E.A., Richardson, J.S., 2001. Direct and indirect effects of sediment pulse duration on stream invertebrate assemblages and rainbow trout (Oncorhynchus mykiss) growth and survival. Canadian Journal of Fisheries and Aquatic Sciences 58, 2213e2221. Schulte, D.M., Burke, R.P., Lipcius, R.N., 2009. Unprecedented restoration of a native oyster metapopulation. Science 325, 1124e1128.

Please cite this article in press as: Morrison, A.E., Allen, M.S., Agent-based modelling, molluscan population dynamics, and archaeomalacology, Quaternary International (2015), http://dx.doi.org/10.1016/j.quaint.2015.09.004

14

A.E. Morrison, M.S. Allen / Quaternary International xxx (2015) 1e14

Smith, E.A., 1991. Inujjuamiunt Foraging Strategies: Evolutionary Ecology of an Arctic Hunting Economy. Aldine de Gruyte, New York. Smith, S.D.A., 2011. Growth and population dynamics of the giant clam Tridacna €ding) at its southern limit of distribution in coastal, subtropical maxima (Ro eastern Australia. Molluscan Research 31, 37e41. Stephens, D.W., Krebs, J.R., 1986. Foraging Theory. Princeton University Press, Princeton. Suchanek, T.H., 1981. The role of disturbance in the evolution of life history strategies in the intertidal mussels Mytilus edulis and Mytilus californianus. Oecologia 50, 143e152. Swadling, P., 1976. Changes induced by human exploitation in prehistoric shellfish populations. Mankind 10, 179e194. Thomas, F.R., 2007. The behavioral ecology of shellfish gathering in western Kiribati, Micronesia 1: prey choice. Human Ecology 35, 179e194. Volterra, V., 1926. Fluctuations in the abundance of a species considered mathematically. Nature 118, 558e560. White, E.P., Ernest, S.K.M., Kerkhoff, A.J., 2007. Relationships between body size and abundance in ecology. Trends in Ecology and Evolution 22, 323e330. Whitaker, A.R., 2008a. The Role of Human Predation in the Structuring of Prehistoric Prey Populations in Northwestern California. University of California, Davis (Ph.D. thesis). Whitaker, A.R., 2008b. Incipient aquaculture in prehistoric California?: Long-term productivity and sustainability vs. immediate returns for the harvest of marine invertebrates. Journal of Archaeological Science 35, 1114e1123.

Whitaker, A.R., 2010. Prehistoric behavioral depression of cormorant (Phalacrocorax spp.) on the northern California coast. Journal of Archaeological Science 37, 2562e2571. Whitaker, A.R., Byrd, B.F., 2014. Social circumscription, territoriality, and the late Holocene intensification of small-bodied shellfish along the California coast. The Journal of Island and Coastal Archaeology 9, 150e168. Wilensky, U., 1999. NetLogo. Center for Connected Learning and Computer-Based Modeling, Northwestern University, Evanston, IL. http://ccl.northwestern.edu/ netlogo/. Wolverton, S., 2005. The effects of the hypsithermal on prehistoric foraging efficiency in Missouri. American Antiquity 91e106. Wolverton, S., Nagaoka, L., Dong, P., Kennedy, J.H., 2012. On behavioral depression in White-tailed Deer. Journal of Archaeological Method and Theory 19, 462e489. Worrapimphong, K., Gajaseni, N., Le Page, C., Bousquet, F., 2010. A companion modeling approach applied to fishery management. Environmental Modelling & Software 25 (11), 1334e1344. Wren, C.D., Xue, J.Z., Costopoulos, A., Burke, A., 2014. The role of spatial foresight in models of hominin dispersal. Journal of Human Evolution 69, 70e78. Wurzer, G., Kowerik, K., Reschreiter, H., 2015. Agent-based Modeling and Simulation in Archaeology. Springer International Publishing.

Please cite this article in press as: Morrison, A.E., Allen, M.S., Agent-based modelling, molluscan population dynamics, and archaeomalacology, Quaternary International (2015), http://dx.doi.org/10.1016/j.quaint.2015.09.004