Scott Heckbert (2013)
MayaSim: An Agent-Based Model of the Ancient Maya Social-Ecological System
Journal of Artificial Societies and Social Simulation 16 (4) 11
Received: 20-Aug-2012 Accepted: 02-Jun-2013 Published: 31-Oct-2013
Bookmark: Share on email Share on print Share on gmail Share on facebook Share on linkedin Share on reddit Share on researchgate Share on stumbleupon Share on xing Share on twitter Share on baidu Share on orkut
This paper presents results from the MayaSim model, an integrated agent-based, cellular automata, and network model representing the ancient Maya social-ecological system. The model represents the relationship between population growth, agricultural production, soil degradation, climate variability, primary productivity, hydrology, ecosystem services, forest succession, and the stability of trade networks. Agents representing settlements develop and expand within a spatial landscape that changes under climate variation and responds to anthropogenic impacts. The model is able to reproduce spatial patterns and timelines somewhat analogous to that of the ancient Maya, although this proof-of-concept model requires refinement and further archaeological data for calibration. This paper aims to identify candidate features of a resilient versus vulnerable social-ecological system, and employs computer simulation to explore this topic, using the ancient Maya as an example. Complex systems modelling identifies how interconnected variables behave, considering fast-moving variables such as land cover change and trade connections, meso-speed variables such as demographics and climate variability, as well as slow-moving variables such as soil degradation.
Social-Ecological System, Archaeology, Cellular Automata, Network Model, Trade Network, Agent-Based Model
Few topics gain as much cross-disciplinary interest as the rise and fall of ancient civilisations. The story of development and demise in complex societies contains narratives of the human endeavour threatened by devastating droughts, greedy rulers, foreign imperialism, and overuse of natural resources, among others. Societies are, however, a set of interacting elements which as a whole express characteristic features, interpreted as emergent properties of underlying processes at multiple scales. Designing a holistic approach to understanding social-ecological systems requires methods which simultaneously observe patterns in many dimensions, a kind of observation for which van der Leeuw (2012) argues that traditional Western science is not very well equipped. An analogy is the example of solving a Rubik's Cube, in that one cannot get the cube 'in order' by dealing first with one side, then the next, and so forth. The only way to arrive at order is by looking at the patterns on all sides simultaneously, and not favouring any particular one at any time (van der Leeuw 2012). This paper presents a method to identify candidate features of a resilient versus vulnerable social-ecological system, and employs complex systems science, using computer simulation to explore this topic using the ancient Maya as an example.
A number of research questions are presented for exploration:
- What dynamics lead to the development of the densely populated and interconnected human geography of the ancient Maya?
- Is it possible to use computational social science to 'grow' the three Maya temporal periods of the Preclassic (1000 BC - AD 250), Classic (AD 250 - 900), and Postclassic (AD 900 - 1500)?
- How does the simulated social-ecological system develop and respond to changing conditions, and what modelled indicators warn of vulnerability?
In order to explore these research questions, a simulation model was designed and calibrated for the landscape of Central America. Model runs produce temporal and spatial patterns that can be understood through examining the underlying assumptions of the different integrated components of the model. MayaSim is a combined agent-based, cellular automata, and network model that represents the ancient Maya social-ecological system. Agents, cells, and networks are programmed to represent elements of the historical Maya civilisation, including demographics, trade, agriculture, soil degradation, provision of ecosystem services, climate variability, hydrology, primary productivity, and forest succession. Simulating these in combination allows patterns to emerge at the landscape level, effectively growing the social-ecological system from the bottom up. This approach constructs an artificial social-ecological laboratory where different theories can be tested and hypotheses proposed for how the system will perform under different configurations.
The model is able to reproduce spatial patterns and timelines somewhat analogous to that of the ancient Maya's history. This proof of concept model requires refinement and further archaeological data for calibration to improve results, although it is noted that there is little empirical evidence by which to validate such models, and such evidence is generally site-specific and discontinuous through time.
The purpose of the model is to better understand the complex dynamics of social-ecological systems and to test quantitative indicators of resilience as predictors of system sustainability. An integrated agent-based, cellular automata, and network model was constructed using the software Netlogo (Wilensky 1999). The full model, code and documentation is available in Heckbert (2012) via the www.openabm.org website, and further description of the model in the context of Maya archaeological literature is presented in Heckbert et al. (in press).
The MayaSim model represents settlements as agents located in a gridded landscape. The model is constructed using the software Netlogo (Wilensky 1999). The software interface, shown in Figure 1, presents the spatial view of the model with graphs tracking model data and a user interface for interacting with the model. The view can be changed to observe different spatial data layers within the model. The model operates at a spatial extent of 516,484 km2 with a 20 km2 resolution. Temporal extent is 650 times steps, each representing roughly 2 years.
Upon model initialisation, base GIS layers are imported using the Netlogo GIS extension. Static cell variables are set therein, dynamic variables are reset to default values set on the user interface, and settlement agents are randomly initialised in the spatial landscape. Imported spatial data include elevation and slope (Farr et al. 2007), soil productivity (FAO 2007), and temperature and precipitation (Hijmans et al. 2005). Data for soil productivity was slightly smoothed using the Netlogo function diffuse to make differences along polygon boundaries less stark, but all other data is unchanged from the cited sources. Data is resampled at a 20km2 resolution using the Netlogo GIS extension (see Wilensky 1999). All GIS data is high resolution and deemed to be of high quality except possibly the soil productivity dataset. The latter dataset can be replaced with finer resolution regional or country-specific data, but because the region extends across several modern countries with inconsistent data sets, the best globally consistent dataset available was used. Future work will aim to replace the existing dataset with the Harmonized World Soil Database (FAO et al. 2009) for finer spatial resolution and consistency across country borders.
Figure 1. MayaSim model interface with interactive controls, spatial view, and graphs tracking model data. Agents operate on a cellular landscape and are connected by links within a network.
The simulation begins with calculations of biophysical variables for water flow and net primary productivity, and these are further used to calculate forest succession, agricultural production, and ecosystem services. Settlement agents interact with the spatial landscape to generate agricultural yield through cropping, derive benefit from local ecosystem services, and generate trade benefits within their local trade network. The combined benefits of agriculture, ecosystem services, and trade drives demographic growth including migration. Simulating the integrated system reveals how the social-ecological system functions through time. Additional agents include a 'migrant' agent who settle new locations, a 'raindrop' agent which routes hydrological surface flow, and a 'network search agent' who traverses links between connected settlement nodes to calculate network statistics. The following section describes biophysical and anthropogenic processes that are tracked in the model.
Biophysical functions: climate, hydrology, soil, primary productivity, forests and ecosystem services
Spatial data for precipitation and temperature (Hijmans et al. 2005) representing current conditions (1950 to present day) is adjusted within the model using assumptions that mimic the paleoclimatic variation presented in Prufer et al. (2011). Variation in precipitation is more pronounced towards the northwest of the case study area. This diagonal north-south / east-west variation in rainfall is created over the current-day GIS data for precipitation, using the function:
where Rj,T is precipitation [mm] for cell j at initial time step T, and CLn is a localized rainfall effect due to the presence of cleared land on neighbouring cells n = 1 … 8 with weighting parameter δ determining the strength of this effect. DFj is the distance [km] of each cell from the top northwest corner of the map and is the furthest distanced cell from this point. RCt cycles from + 20% to -10% linearly over a 56 time step cycle, and t = 1 … 650. This function serves to reduce and increase rainfall cyclically, with a more pronounced effect further towards the northwest.
A surface-flow hydrological model was devised to recreate past conditions based on elevation and rainfall. These data are used to calculate surface flow and location of potential seasonal standing water, consistent with Reaney (2008). Each cell generates a mobile 'raindrop' agent containing precipitation volume Rj,T. The raindrops follow the elevation data, repeatedly moving to the adjacent cell with the lowest elevation (and considering the summed volume [mm] of raindrops already at that location). If raindrops cannot move (i.e., a location is flooded), the raindrops 'pool' and form river and lake patterns, as shown in Figure 2 depicting simulated water flow under climate variation assumptions. This flow pattern can be validated against hydrological processing functions in GIS software. The function serves to move water based on elevation, and can generate the spatial distribution and surface water flow as precipitation varies across the climate cycle. See equations in model code readily explored in Heckbert (2012). The result of these equations is a finely-resolved spatial pattern that is able to recreate the potential river, lake, and wetland systems under climate variation assumptions. An important caveat of the model assumptions, however, is the lack of heterogeneity in infiltration rates across the region. In reality, the surface water in many areas drains rapidly to the aquifer owing to the karst limestone bedrock, however spatial data to inform this process was not able to be sourced, and as such a constant rate of infiltration is assumed across the landscape, and is a topic for further research effort.
Figure 2. Simulated surface water flow [m3] (blue on white background, darker blue being greater flow) based on climate variation assumptions of a) -15%, b) 0%, and c) +15% change in annual precipitation in equation 1.
The base GIS layers for rainfall and temperature are used to calculate net primary productivity, which in turn is used to calculate forest succession. GIS soil data is used to calculate agricultural productivity, and all these combine to calculate provision of ecosystem services. While we recognize that these relationships are complex, simplifying assumptions allow initial representations as methodologies are further refined.
Forest succession operates as a cellular automata model, where the state of a cell is dependent on internal conditions and is influenced by the condition of neighbouring cells. In this proof-of-concept stage, cells take on one of three general forest states that represent climax forest, secondary regrowth, and cleared/cropped land, referred to as state 1, 2 and 3 respectively. The forest state is decremented for 3.5% of randomly selected cells, to represent natural disturbance. The disturbance rate is linearly amplified by population density of nearby settlements to represent local wood harvesting, to a maximum of 15%. Cells advance in their forest state based on the time since last disturbance and the relative net primary productivity of the cell. Once the time since last disturbance is above a threshold (40 * NPPmax j, t / NPPj, t years for secondary regrowth and 100 * NPPmax j, t / NPPj, t years for climax forest, to account for spatial variation in net primary productivity) the forest converts to the new state. For conversion to climax forest, a cellular automata function is applied that requires a number of neighbouring cells to also contain climax forest. This rule represents the need to have local vegetation for seed dispersal.
Net primary productivity NPPj, t [gC m2 -1 yr-1] is a function of precipitation and temperature, calculated based on the Miami model (Lieth 1972) as:
where Rj,t is precipitation [mm] and Tj is temperature [degrees C].
For each cell, agricultural productivity AGj,t is calculated as:
where SPj is soil productivity (FAO 2007) [index 1-100], Sj is slope [%], WFj is water flow calculated as the sum volume of water agents traversing any given cell j, as depicted in Figure 2, and SDj,t is soil degradation [% loss of productivity]. Soil degradation occurs at a constant rate of 1.5% per time step for each cropped cell. This constant rate of soil degradation is obviously a simplifying assumption, as soil management and maintenance of soil productivity is in itself a complex social-ecological system. The β parameters are weightings for calibration, and presented in Heckbert (2012).
Ecosystem services are modelled by quantifying the availability provisioning services (as defined in MA 2005; TEEB 2010) relating to water, food, and raw materials. This is a subset of ecosystem services and does not include a full set of indicators which would incorporate supporting services (for example erosion prevention), habitat services (such as maintenance of genetic diversity) or cultural services (such as inspiration for culture, art, and design). The current ecosystem services equation incorporates a subset of four important services provision based on arable soils, precipitation, access to available freshwater, and timber resources. Ecosystem servicesESj are calculated as:
where AGj,t is taken from equation 3, Rj,tis taken from equation 1, WFj,tis the simulated water flow volume, and Fj,tis the forest state [1-3], ESDj,tis a catch-all proxy variable for all other ecosystem services degradation [%] as a function of population density.
Anthropogenic functions: agriculture, trade, and demographics
Each settlement agent i maintains at least one cell j for generating agricultural yield. Settlements perform an agriculture benefit-cost assessment considering the costs of production, travel cost given the distance of the cell from the settlement site, and with larger settlements achieving economies of scale, modelled as;
where BCAj,t is the total benefit provided from agriculture yield, κj, α, and φ, are crop yield and slope parameters, AGj,t is again taken from equation 3, γ is the establishment cost of agriculture (annual variable costs), Oj is the agriculture travel cost as a function for distance from the city and a per km cost parameter, and Pj,t is population of the settlement.
The benefit-cost of agriculture function generates yields that are spatially distributed based on individual conditions of the cells and the location of settlements. Costs of production, including distance from settlements, results in adding cropped cells, generating yield and increasing population, which in turn add more cropped cells, but causes soil degradation. The system adjusts over time in response to the spatially-explicit agricultural benefit-cost.
A series of functions represent trade within a spatially connected network of agents. It is assumed that through the process of specialization, settlements that are connected to one another within a network will generate benefits from trade. It is assumed a larger network produces greater trade benefits, and also the more central a settlement is within the network, the greater the trade benefits for that individual settlement. To model these benefits, settlements are connected via a network of links that represent trade routes. As a simplifying assumption of how they connect together, it is assumed when a settlement reaches (or drops below) a certain size, they will add routes (or allow routes to degrade) to nearby settlements within a 40 km radius. At each time step, the size of the local network is calculated as well as each settlement's centrality within that local network, further discussed below.
Combining the functions for agriculture, ecosystem services, and trade benefit, total real income per capita RIi is calculated as:
where Nj,t is the network size [# nodes], Cj,t is the centrality [degree] and TCi is the travel cost, and ϑ parameters are prices for agriculture, ecosystem services, and trade, respectively. Benefits from agriculture are calculated only for cells under cropping production AGJi,t = 1 … n whereas ecosystem services are calculated encompassing the entire 'area of influence' of each settlement IJi,t = 1 … m which is based on the population size of the settlement, increasing linearly to a maximum of 40 km in diameter for the most populous settlements (those with populations greater than 15000 people), as taken from Heckbert et al. (in press) and interpreted from Chase and Chase (1998). Travel cost measures the relative 'friction' of different land cover types, and is represented as:
where Sj is slope and WFj,t is simulated water flow volume, both described in previous equations, resulting in areas of higher slope being relatively more costly to travel through, mitigated by the presence of flowing water for canoe transport.
The terms NJi,t network size, and Ci,t centrality in equation 6 are calculated using a network search algorithm custom-built for this application. The network search algorithm is the most detailed procedure in the model, and uses a recursive search function using 'network walking agents' who travel every possible combination of routes along a network path. These report to the network nodes (settlements) the total size of the local cluster, and the position of the node within that cluster (degree). The network search algorithm is the lengthiest procedure in terms of amount of code written, is the most computationally expensive procedure in the model (more than doubles the total run time of a simulation of 650 time steps). The outcome of the procedure is depicted in Figure 3 for time step 250, with the network links coloured black and the range of values on the visual scale halved compared to following figures of the same indicator, for visual purposes.
Figure 3. Simulated 'trade strength' [last term of equation 6] (red on white background, darker red being greater trade strength) based on centrality, size of local network cluster, and travel cost of each settlement (temple icon), connected via trade routes (black links). A local cluster is identified for visual purposes with a blue oval over the same cluster, with a) and b) being zoomed-in views of the wider system depicted in c).
We see in Figure 3 an area highlighted by a blue oval, zoomed-in at three different scales. In Figure 3a), we see that this local cluster is not connected to adjacent clusters, and Ni,t network size [# of nodes in network] is small. In Figure 3b) we see that the trade strength (red colouring) is greater for larger-sized local clusters given the larger value of Ni,t. Lastly, in Figure 3c) we observe that central nodes in larger clusters (those located in the darkest red regions), have greater trade strength compared to nodes at the periphery of the same cluster.
The value of each contribution to RIi,t is determined by the weighting ϑ, which is static for ϑAG and ϑES; however, the value of trade ϑTRt is dynamic. Specifically, ϑTRt increases each time rainfall decreases, according to the climate variation assumptions. The assumption here is that settlements specialize production within an overall trade network to increase the value of trade goods relative to other commodities. The effect is to linearly increase the value of trade each time the climate cycle is in decline. This assumption is obviously rudimentary, and further formulations are explored.
After determining RIi,t, settlement demographics account for births, deaths, and migration. The birth rate is assumed to remain constant at 15%, while death rate and out-migration decrease linearly with increased RIi,t per capita, with a maximum out-migration rate of 15% and a maximum death rate of 25% per annum. Settlements with a population below a minimum number required to maintain subsistence agriculture are deleted. Settlements that register out-migration above a minimum threshold of the number of people required to maintain subsistence agriculture create a 'migrant agent'. The migrant agent uses a utility function (Heckbert et al. 2010) to select locations to create a new settlement. The migration utility function is calculated as:
where λ parameters are weightings for travel cost and ecosystem services, and ESj,t is taken from equation 4, and Dj is the distance from the origin settlement to the potential new settlement site.
Model outcomes are dependent on parameterization of equations 1-8. This stage of model development is offered for proof of concept, for transparency, and to receive feedback on model assumptions. Model code and parameterization is available for download via the http://www.openabm.org website in Heckbert (2012).
This section describes results and sensitivity analysis from simulation runs. This proof-of-concept stage model reveals familiar patterns analogous to the Maya historical timeline and offers a quantitative basis for these assessments. The model was initialised according to default values (Heckbert 2012). Simulation runs were conducted with a temporal extent of 650 time steps with a spatial resolution of 20km2. Figures 5-7 describe model results, depicting mean values for 20 runs (with confidence intervals [a=0.05] where relevant), and reporting on social-ecological indicators.
Figure 4 presents spatial outcomes for five indicators, at 200 time step intervals. Population density, forest condition, settlement trade strength, soil degradation, the condition of ecosystem services and forests each contain a narrative describing the development and reorganisation of the simulated social-ecological system. By time step 200, settlements have expanded into all regions, first occupying areas with greater ecosystem services and progressively growing with agriculture development. Population densities are higher in areas where settlements have clustered and formed local trade connections. By time step 400, as the value of trade increases, the population dramatically increases, extending local trade connections to 'global' connectivity. A notable fringe exists between the connected network and the periphery . The centre of the global trade network is approximately in the broad region where ancient Maya captials of Tikal and Caracol existed. The condition of the forest is markedly changed, with only small patches of climax forest remaining in agriculturally unsuitable areas, forming ecological refugia within the near-completely settled landscape. Dramatically, by time step 600 the trade network has disintegrated, the centre of the most densely populated areas is nearly entirely abandoned, leaving only a small number of locally connected settlements of any notable size in what was once the fringe of the globally connected network. Abandoned cropland and significantly decreased fuelwood harvesting allows broad-level secondary regrowth, and climax forest eventualy expands out from its refugia to an extent similar to pre population expantion levels. However, soils and ecosystem services remain severely impacted and large scale resettlement of degraded areas is not possible.
Figure 4. Population density, forest condition, settlement trade strength, soil degradation and the condition of ecosystem services are shown at 200 time step intervals. Darker colouring shows increased a) population density (blue), c) trade strength (red), d) soil degradation (red), and e) ecosystem services (green). Forest condition b) depicts three states of cleared / cropped cells (yellow), secondary regrowth (light green) and climax forest (dark green).
The model reports quantitiative indicators through time and can be used to drill down in order to explore the dynamcis of the development and reoganisation observed. Figure 5 presents the total population of all simulated settlements and contributions to real income by ecosystem services, agriculture and trade, respectively. In the first _ of the simulation, ecosystems services provide the majority of value. Ecosystem services values are superseded by agriculture by time step 150, and both are superceded by trade around time step 350. Neither recover as trade values decrease in the latter half of the simulation run, and population adjusts accordingly.
Figure 5. Total population [# people, primary axis] of all simulated settlements over time, and contributions to real income by ecosystem services, agriculture and trade [# proxy value units, secondary axis]. Ecosystem services values are superceded by agriculture by time step 150, and both are superceded by trade around time step 350. Neither recover as trade values decrease in the latter half of the simulation run.
The rapid change in the value of trade can be explained by examining Figure 6, which depitcs the total number of settlement nodes, the number of nodes within the largest cluster, and total natural capital, which is represented as the total sum of ecosystem services taken from equation 4. The network grows from local clusters to a near-globally connected system. Periodic perturbations give the clusters structure which inevitably form the 'skeleton' of the global structure. However, when natural capital has reached its lowest level, perturbations result in cascading failure in the network. Although natural capital recovers, the number of settlements does not due to the lack of trade network structure.
Figure 6. Total number of settlement nodes and number of nodes within the largest cluster [primary axis], and total natural capital [total sum of ecosystem services values, secondary axis]. The network grows from local clusters to a near-globally connected system through growth in link connections and periodic perturbations which give the clusters structure. However, when natural capital has reached its lowest level, perturbations result in cascading failure in the network.
Figure 7 depicts soil degradation and forest condition by three states of cleared / cropped cells (state 1), secondary regrowth (state 2) and climax forest (state 3). The initial period, corresponding to roughly the first third of the simulation run, is descibed by accelerated decline of climax forest and inhibited regrowth is a result of cropping and timber harvesting for construction and fuelwood as a result of increasing population levels. The following period, corresponding to roughly the second third of the simulation run, shows continued repressed regrowth, and increased areas of cleared/cropped land as population pressure results in marginal lands being put under agricultural production. As a result, soil degradation increases at its highest rate during this period. The last third of the simulation run shows a rapid decline in cleared/cropped land with corresponding large scale secondary regrowth, and eventual succession into climax forest which recovers to near pre-development levels.
Figure 7. Forest condition [% of total area, primary axis] by three states of cleared / cropped cells (state 1), secondary regrowth (state 2) and climax forest (state 3), and soil degradation [total proxy soil productivity units lost, secondary axis]. Accelerated decline of climax forest and inhibited regrowth is a result of cropping and timber harvesting. Increased cleared/cropped area results in higher soil degradation rates, with soils beginning to recover after climax forest again becomes dominant.
The timeline for each indicator presented in Figures 5-7 reveals the underlying dynamics of the modelled social-ecological system. Table 1 presents an integrated narrative of model indicators at 100 time step intervals, akin to the analogy of examining a Rubik's Cube by looking at the patterns on all sides simultaneously (van der Leeuw 2012). This complex systems perspective highlights that a reductionist method of identifying factors contributing to resilience and vulnerability cannot be used to identify a single cause (or a linear progression of single causes) of the unravelling of the social-ecological system.
Table 1: Narrative description of modelled indicators for population, forests, soils, and trade in 100 time step intervals
PopulationNumber of settlements grows based on migration to unsettled areas with high ecosystem servicesPopulation increases with agriculture developmentPopulation increases with trade valuePopulation reaches height, significant perturbations with climate variabilityPrecipitous decline in populationPopulation does not recover to prior levels and settlements of any significant size are located at the former fringe of the most densely populated areas
ForestsGradual decline of climax-dominated forest and gradual increase of secondary regrowth and cleared / cropped areas.Roughly even distribution of three forest states.Climax forest reduced to 10% of landscape area, timber harvesting and cropping suppresses secondary regrowth.Cleared / cropped land accounts for 70% of land surface, climax forest reduced to refugia.Broad-scale secondary regrowth.Climax forest begins to emerge from refugia to nearly regain pre-development levels.
Soils and AgricultureAgricultural expansion into prime soils.Increase of cropped areas in broadly settled landscape.Expansion of agriculture into marginal and more distant locations, soil degradation rate increases.Broad-scale cropping albeit with diminishing returns relative to land area under production. Soil degradation reaches height. Abandonment of crops and loss of economies of scale in agriculture production, legacy soil degradation inhibits re-sowing crops.Soils remain unproductive and areas suitable for agriculture are predominantly in the former fringe of densely cropped areas.
Trade and network structureLittle connectivity of trade routes, with small network cluster sizes of