Coupled land use and ecological models reveal emergence and feedbacks in socio‐ecological systems

Acknowledgements: This work was supported by an EPSRC Doctoral Training Centre grant (EP/G03690X/1). Supplementary material (Appendix ECOG‐04039 at ). Appendix 1.


Introduction
Interactions between human and natural systems have profound global consequences, making the identification of sustainable practices a key priority and challenge. The complexity of these interactions means that they are difficult to fully understand or predict, producing unexpected consequences from well-intentioned interventions (Alberti et al. 2011, Malawska et al. 2014. Nonlinear responses, feedback mechanisms, time lags and shocks are all characteristic of socioecological systems (Liu et al. 2007, Filatova and Polhill 2012, Filatova et al. 2016. Understanding these dynamics requires the integration of research and management approaches that have typically adopted distinct methods to focus on either social or ecological issues in isolation (Voinov andShugart 2013, Liu et al. 2015).
Integration of approaches emerging in both ecological and socio-economic research offers substantial potential for improved representation, exploration and prediction of socio-ecological dynamics. In particular, process-based modelling approaches that have become widespread in recent years share key characteristics in focusing on fundamental system behaviours via 'bottomup' model architectures, in principle allowing for the coupling of models of different sub-systems.
For instance, individual-or agent-based modelling (hereafter ABM) has become a common method of studying underlying interactions in complex social and ecological systems on the basis of individual-level characteristics and behaviours within varied environmental, social and economic contexts (Matthews et al. 2007, Grimm and Railsback 2013, Synes et al. 2016. In ecology, ABMs are increasingly used to study animal movements (Tang andBennett 2010, Watkins et al. 2015), help identify threats to populations (Wiegand et al. 2004), test the efficacy of conservation strategies (Synes et al. 2015, Aben et al. 2016) and predict the impacts of land-use scenarios (Gimona et al. 2015). In social and land-use science, ABMs have been used to study the responses of land managers, communities and societies to environmental change (Batty et al. 2012, Brown et al. 2017, to explore the evolution of social behaviours (Macy andWiller 2002, Smith andConrey 2007) and to contribute to political or practical strategies for sustainable development (Pahl-Wostl 2002, Berger et al. 2006).
However, these parallel fields currently do little to consider feedback dynamics from each other.
Few ecological studies have considered temporally dynamic landscapes, and those that have done so have used pre-defined changes in land management and the environment (Bithell et al. 2008, Imron et al. 2011, Gimona et al. 2015, Synes et al. 2015, Marshall et al. 2018, thus not including feedbacks between ecological and environmental processes (Keith et al. 2008, Franklin et al. 2014. Similarly, models of social or land-use systems and their environmental impacts tend to incorporate a simplistic or static representation of the environment (Veldkamp and Verburg 2004, Luus et al. 2013, Brown et al. 2017, in which humans act either as a driver or a user of the environment, but rarely both (Matthews and Selman 2006, Filatova et al. 2013, Schulze et al. 2017). The few existing models of feedback dynamics suggest a largely unrealised potential for uncovering important emergent effects through integrated modelling (e.g. Mathevet et al. 2003, Monticino et al. 2007, Rebaudo et al. 2011, Carrasco et al. 2012, Le Page et al. 2013, Polhill et al. 2013).
One reason that integrated modelling has not yet achieved its potential for providing insights into socio-ecological dynamics is that model development is both conceptually and practically challenging. Existing models often provide advanced, validated and well-understood representations of particular sub-systems. However, existing models also often include hidden inconsistencies in assumptions, complexities that are unnecessary for an integrated application, or exogenous parameters that represent factors endogenous to the coupled system (Grimm et al. 2005, Marohn et al. 2012, Luus et al. 2013. Another key consideration is the potential for differences in the temporal and spatial scales at which the interacting systems (and hence models) operate (Janssen andOstrom 2006, Malawska et al. 2014). These complexities mean that model coupling must be carried out with considerable caution if meaningful results are to be produced.
Here, two existing ABMs are coupled to explore the feasibility and utility of integrated modelling of farmer decision making with wildlife responses in a single socio-ecological system. We use a hypothetical application to compare system dynamics with and without two-way links between social and ecological systems, and with and without an explicit representation of species population dynamics as well as habitat.

Case study
A key issue affecting the sustainability of socio-ecological systems is the interaction between agriculture and pollination. Approximately 35% of global food crop production relies on pollination (Klein et al. 2007), as do almost 90% of wild flowering plants (IPBES 2016). Global declines in pollinator populations are therefore of great concern, especially in the context of growing food demands and consequent agricultural expansion and intensification that contribute to those population declines (Potts et al. 2010, Whitehorn et al. 2013, van der Sluijs et al. 2015. The ability of models to simulate the feedbacks between agricultural production and pollinator ecology is therefore of great importance. We developed a hypothetical system for a virtual pollinator species (illustrative of real pollinator species, and part of the virtual ecologist approach: Zurell et al. 2010) in an agricultural landscape under increasing pressure from demand for food, in order to explore the sensitivity of the system to processes and factors controlling socio-ecological interactions.

Model coupling
RangeShifter (Bocedi et al. 2014), an ABM of animal demography and dispersal, was integrated with CRAFTY (Murray-Rust et al. 2014), an ABM of land-use dynamics, through loose coupling, i.e. the models interacted through file-sharing.
An important feature of CRAFTY for this application is that each pixel of land has a number of "capitals" representing the land's potential for each included ecosystem service, rather than landuse types being proxies for ecosystem services (Murray-Rust et al. 2014). Natural capital is the stock of natural assets from which humans can derive ecosystem services (Costanza et al. 1997), and the representation of natural assets has previously been suggested as a key feature to study feedbacks in socio-ecological systems (Luus et al. 2013), leading to greater realism by allowing land parcels with the same land-use type to differ in terms of productive potential (as defined by natural capital) and hence ecosystem service levels. Here, we adopt the use of the term 'capitals' to denote specific components of both natural and anthropogenic productive potential: namely, the properties of land parcels relating to crop productivity, timber productivity, and grassland productivity (Murray-Rust et al. 2014, Harrison et al. 2015. Social agents compete for land on the basis of their ability to utilise capitals and provide ecosystem services, for which societal demands are defined. Therefore, ecosystem service levels are the basis for agents' decisions about land management, thus making CRAFTY suitable for integration with models of natural systems, allowing feedbacks between changes in land use and ecology. For this case study, the simulated pollinator population has a direct influence over crop productivity (reduced in locations without pollinators), while other capitals remain static. A change in crop productivity then changes the productive potential of the land, and hence the competitiveness of land-use agents, potentially leading to land-use change.
RangeShifter operates on an annual time-step, whilst one time-step of CRAFTY incorporates a full set of agent decisions about land-uses that do not have a fixed timescale. RangeShifter was modified to pass updated capitals to CRAFTY, on the basis of which CRAFTY simulated one set of agent decisions (an initialisation step and one subsequent time-step) to represent a single year.
RangeShifter then loaded the new landscape provided by CRAFTY, and advanced to the subsequent year in its own simulation ( Figure 1). The timescale implications of this model coupling were that animal population dispersal and reproduction, and land-use change all occur annually.

Landscape
An artificial landscape was created to allow for the simple representation of a scenario where demand for food (meat and crops) is increasing but the production potential of the land has a finite limit. The landscape was defined in terms of a simplified set of three capitals representing the productive potential of each cell in the landscape for ecosystem services: crop productivity, livestock (grassland) productivity, and timber productivity. Capitals were created randomly for each cell from a uniform distribution (0 < x ≤ 1), specifically to avoid spatial autocorrelation that would influence land uses, as one potential outcome of interest was the emergent spatial configuration of land uses in the coupled model. The "crop" is defined as generic, since its reliance on pollination is one of the independent variables being tested; crop selection is not considered, with all crop farmers producing the same crop. In this landscape of 100 x 100 cells (equating to a per-cell land unit size of 25 ha at 500 m resolution), each cell could either be managed by a unique agent implementing one type of land use or left unmanaged. Agents belonged to the following types: high intensity crop farmers, low intensity crop farmers, high intensity livestock (pasture) farmers, low intensity livestock (pasture) farmers, foresters (managed forest). Land-use conversion costs and delays were not included; when a new agent takes over a cell, ecosystem services are immediately produced to the maximum potential of that cell. However, every agent type was assigned a 'giving-in' and a 'giving-up' threshold (see Supplementary Table S1). The 'giving-in' threshold determines the maximum competitive deficit an agent is willing to accept before allowing a new agent (and therefore land use) type to take over. The 'giving-up' threshold determines the minimum acceptable return, below which the cell is abandoned. Similarly to previous applications of CRAFTY (Brown et al. 2016), high intensity farmers were defined with higher productivity than low intensity farmers, but also with greater dependency on the capitals available in their cell (see Supplementary Table S2). The capacity of CRAFTY to assign unique productive abilities, asset sensitivities and threshold values to each agent were not used here, and social networks between agents were also not implemented. The agent population therefore differed only in sensitivities to asset levels, abilities to produce ecosystem services, and tolerance for competition and returns.

Species
The hypothetical insect pollinator populations were subject to carrying capacities which depended on land use type (Supplementary Table S3). Forests supported the highest carrying capacity for the pollinators, followed by unmanaged land, then cells managed by low intensity farmers. Cells managed by high intensity livestock and crop farmers were modelled as unsuitable, and the pollinators could not reproduce in these land-use types. A female-only population model was used, in line with research on pollinators that has found male production to be non-limiting in general and particularly when populations are declining (Beekman andvan Stratum 1998, Whitehorn et al. 2012). Here, an 'individual' represented a single colony of pollinators rather than individual insects.
The species' population dynamics were modelled at the cell scale, i.e. the individuals (colonies) present in each cell represent a distinct sub-population, and density-dependent emigration operated between cells. A number of population parameters were varied to study their impact on model results: maximum fecundity (Rmax), carrying capacity reduction factor (dK), and the inclusion or exclusion of long-distance dispersal (LDD) in the pollinator movement model. These parameters were selected for their substantial impacts on population dynamics, and the recognised need to explore their effects in population models through sensitivity analyses (Naujokaitis-Lewis et al. 2009, Lurgi et al. 2015. For the full RangeShifter parameter specification, see Supplementary Table S3.

Model type: coupled versus uncoupled
To test the effect of incorporating population dynamics and species ecology into modelled landscape change, two versions of each simulation were run. In one, the coupled models included bi-directional feedbacks between land-use and pollinators, meaning that land-use affected habitat suitability for pollinator populations and that pollinator distributions affected crop productivity, changing the competitiveness of land-use. In the other, the uncoupled models included only unidirectional feedback, meaning that land-use affected habitat suitability for pollinator populations, but pollinator distributions did not affect crop productivity. Instead, the uncoupled model took suitable habitat as a proxy for species presence; that is, all suitable habitat was assumed to contain pollinators.
Crops in each cell containing a pollinator population (in the coupled model), or containing suitable pollinator habitat (in the uncoupled model), and in the eight neighbouring cells were pollinated (i.e. assuming a maximum pollinator foraging distance of 500 m in line with empirical findings for different species of bumblebee (Osborne et al. 1999, Darvill et al. 2004). In both models, pollinated cells retained the full crop productivity capital value; in the absence of pollination, the crop yield of that cell was multiplied by a factor varied in separate model runs from 0.1 (reduced to 10%) to 0.9 (reduced to 90%) in increments of 0.2. We herein refer to this reduced crop yield as "crop yield in absence of pollination". These values were chosen to be representative of the variable dependence that different crops have on pollination. Pollination can be essential for some crops (without pollination, production is reduced to 10% or less), but for other crops pollination has little effect (without pollination, production is reduced to 90% or more) (Klein et al. 2007).
Including model coupling, population parameters and the variations in crop yield in absence of pollination, a total of five model parameters were varied across different simulations in a full factorial design (Table 1).

Simulations
The initial demand for crop produce was set at 2.5 x initial demand for livestock produce, approximately matching the proportions of world demand for crops and livestock (Valin et al. 2014). A 'spin-up' CRAFTY simulation was run for 20 time-steps with constant demand, allowing the agents to achieve a stable spatial distribution at initial demand levels. The resulting land-use map was used to initialise all of the main simulations. The initial demand level was below the productive capacity of the landscape, allowing for a mixture of high and low intensity farming, and for 757 of the 10,000 cells to be unmanaged.
Demand for services was defined exogenously to CRAFTY, from an assumed non-spatial human population. The same demand curve was used for every simulation, beginning with 10 years of constant demand to allow the pollinator populations and land-use agents to stabilise (required due to the differences in crop yield in absence of pollination, but omitted from the results). Every simulation was initialised with pollinator populations at carrying capacity and occupying every suitable cell in the landscape. This meant that both coupled and uncoupled models started on the basis that all suitable habitat contained pollinators -only the uncoupled model continued on the assumption that this was always the case. The constant demand period was followed by 50 years of linear annual increases resulting in an overall 74% increase in demand for both livestock and crop produce, matching the mean increase projected by Valin et al. (2014) for 2050. Demand for forests in the case study encompassed both demand for timber and the protection of forests for conservation, and decreased to zero by the end of the 50-year period. This represented a scenario in which forest protection is gradually reduced due to the increasing demand for food. It therefore assumed that no consideration was made (in terms of demand from the non-spatial population) for the natural capital of forests, and the ecosystem services they provide -including their potential role in supporting pollinator populations.
The spatial auto-correlation (Moran's I) of land-uses was monitored annually throughout the simulations to identify any emergent spatial effects from the interaction of (spatially explicit) pollination and (spatially explicit but independent) land use changes. Ten replicates of CRAFTY (with the same demand curves) were also run, independent of RangeShifter, to identify whether spatial effects could emerge from the modelled land use processes alone. The effects of ecological processes on spatial auto-correlation were tested by comparing species parameterisations and the uncoupled simulations. Both expected values (assuming the underlying spatial processes were random) and observed values were calculated and then a z-score and p-value were computed to test for a statistically significant difference.

Results
Both the coupled and uncoupled models demonstrated the dynamic of agricultural intensification driving declines in pollinator habitat, and the feedback cycle of declining pollinator habitat driving decreases in crop yields, leading to greater intensification and greater loss of available pollinator habitat (Figure 2a, b). In the context of the model, use of the term "agricultural intensification" refers to increases in the number of high intensity farmers -where crop productivity may be higher, but at the expense of natural habitat. As expected, the increasing demand for food in combination with a decreasing protection of forests also led to both the intensification and expansion of farming. However, this gradual intensification did not result in a steady increase in crop supply; instead, severe collapses in supply occurred as pollination decreased and reduced the mean value of the crop productivity capital (Figure 2c). This reduction in capital levels occurred for high intensity crop farmers only when crop yield in the absence of pollination was low (Figure 2d Figure S3). For all other land-uses, there was no difference in Moran's I values from those expected by chance. The variation in spatialautocorrelation for high intensity livestock farming in coupled simulations came from the ecological parameterisations, which had no effect in the uncoupled model (Figure 6a). Lower maximum fecundity and lower carrying capacity resulted in lower Moran's I (i.e. less spatialautocorrelation of high intensity livestock farmers), but the inclusion of long distance dispersal had no effect. In the coupled models when crop yield in the absence of pollination was low (0.1), most spatial aggregation occurred rapidly between simulation years 20 and 25 (when pollinator populations were decreasing most rapidly), regardless of the ecological parameterisations ( Figure   6b). This was driven by the lack of suitable habitat for pollinators in cells managed by high-intensity livestock farmers, reducing populations in neighbouring cells and therefore making those cells less suitable for crop production. After this rapid change, the level of spatial-autocorrelation diverged towards the values in Figure 6a, dependent on parameterisation.

Discussion
The coupled socio-ecological model presented here reveals important, emergent system dynamics that would be obscured by established modelling approaches that focus on either social or ecological processes. These emergent dynamics include the development of socio-ecological feedback loops, spatial clusters of different land-use types, and exposure of intrinsic uncertainty in potential outcomes reflecting the interplay of social and ecological characteristics.
The hypothetical case study used here, in which agricultural intensification reduces pollinator population sizes, pollination services and hence crop yields, is designed to reflect observed links between agricultural intensity and pollinator populations, or biodiversity more generally (Kremen et al. 2002, Potts et al. 2010, Lanz et al. 2018. At the same time, the simplified modelling approach allows for exploration of the effects of uncertainties about land-use impacts and species ecology.
At a general level, we find the key feedback between agricultural intensification and declining pollinator populations is robust to the modelled range of these uncertainties but not to model uncoupling, which breaks the link between agricultural decision-making, yields and species dynamics.
Our experimental design includes an increasing demand for food that requires expansion and intensification of agriculture. In the coupled model, pollinator populations decline as the land becomes dominated by high intensity farmers that cannot support populations, and the resulting decrease in crop supplies drives further, rapid intensification and expansion of farming. These threshold changes occur earlier and become more severe as crop yield in the absence of pollination decreases, demonstrating the importance of understanding both the impacts of land use on pollination and other ecosystem services, and the dependency of land use on those services. Furthermore, the divergence of coupled and uncoupled model results following variations in parameter values controlling species fecundity, carrying capacities and dispersal range shows that species ecology is a key aspect of these dynamics, and highlights the dangers of assuming that suitable habitat is a reasonable proxy for population presence (Fordham et al. 2012, Conlisk et al. 2013. We found that the modelled system had the greatest shortfall in food production when parameter values for pollinator carrying capacities and maximum fecundity were at their lowest, and the smallest shortfall when these parameter values were at their highest. Pollinator population dynamics are known to be sensitive to population sizes, growth rates, densitydependence, fragmentation and genetic diversity (Purvis et al. 2000, Henle et al. 2004, Melbourne et al. 2004, Whitehorn et al. 2011). In the case study used here, failure to consider the ecology and population dynamics of pollinators leads to substantial overestimation of the modelled system's ability to meet future food demand levels.
The emergence of spatial land-use patterns, particularly where crop yields were more reliant on pollination, occurred in both coupled and uncoupled simulations. Once again, the coupled model exhibited greater sensitivity to ecological parameters, with lower maximum fecundity and carrying capacities resulting in less clustering of high intensity livestock farmers. Independent simulations of CRAFTY (without RangeShifter) showed no clustering effect; high intensity livestock farmers remained marginally over-dispersed. This confirms it was the interaction between land-use and ecological processes that caused spatial land-use patterns to emerge in our case study simulations.
The presence of land uses that neither support nor require pollinators (cells managed by high intensity crop and livestock farmers as modelled in this study) decreases the potential sources of foraging pollinators for neighbouring land uses; this increases the likelihood that the productivity of adjacent land uses reliant on pollination will be reduced, and that they will therefore become less competitive and more likely to be displaced. Such a mechanism has the potential for a positive feedback loop similar to that demonstrated by the model of segregation of Schelling (1971). It is also consistent with patterns of agricultural intensification, from local to global scales, that both drive and are driven by biodiversity losses -although limits clearly exist to this process in reality (Lanz et al. 2018). The rate of cluster formation in our study may be affected by the initialisation of the models with spatially uncorrelated land uses. While this provides a neutral basis on which to compare the coupled and uncoupled approaches, it does not account for the effects that realworld spatial autocorrelation in land uses may have in either increasing or reducing the rate of clustering of particular land use types. Similarly, we do not simulate ecological interactions beyond those implicit in the ecological parameters used above.
This study demonstrates a system modelling approach that can capture important socio-ecological dynamics between pollinator ecology and agricultural intensification, and, in principle, between land use and ecosystem services more generally. It is essential that such dynamics are better understood and modelled if meaningful steps towards sustainable agricultural systems are to be made (Haberl et al. 2016, Rissman andGillon 2017), particularly in areas of high reliance on pollination or other ecosystem services (Lautenbach et al. 2012). However, much ground remains to be covered. In the case of pollinators, more precise and species-specific responses to land use (including pesticide usage) and environmental conditions would allow for practical assessments of population vulnerabilities and conservation plans (van der Sluijs et al. 2015, Aguirre-Gutiérrez et al. 2017. Dynamics similar to those we simulate have led to various farming adaptations in reality, including establishment of wildflower strips, the use of farmed pollinators or even artificial pollination (Potts et al 2016). We deliberately excluded such adaptations here in order to focus on underlying feedback effects, but the magnitude of impacts that emerge from the coupled models clearly forms a basis for simulation of adaptive processes in the future. More generally, feedbacks between ecosystem services and elements of biodiversity warrant increased attention, particularly where they relate to soil biodiversity (Vukicevich et al. 2016, de Valença et al. 2017, Delgado-Baquerizo et al. 2017. The timescales of such feedbacks can also have crucial implications for socio-ecological dynamics (Lafuite and Loreau 2017), as can relatively neglected aspects of genetic and functional diversity (De Palma et al. 2017). Such effects are also strongly modulated by the knowledge and attitudes of land managers, which we did not consider here, but which can be varied in the model framework we used (Brown et al. 2016, Pe'er et al. 2017).
In any case, future research is increasingly likely to involve the coupling of models to study interacting systems (Synes et al. 2016), and realistic couplings require that feedback mechanisms are implemented between the study systems. A model of animal movement and population dynamics will often therefore require a model of the changing landscape or environment in which the species lives. The environmental modelling requirement can vary greatly depending on the species, and may include models of land-use, climate, vegetation, hydrology, or even finer scale environments (Franklin et al. 2014). Such model pairs can be integrated simplistically by creating a sequence of landscapes to be loaded in a time-series by spatial models of animal movement.
However, a one-way integration such as this makes the assumption that the animal has no influence on the landscape upon which it exists, an assumption that is rarely if ever accurate. As this study has shown, the impact of pollinators on crop yield can radically change the spatial and temporal properties of agricultural intensification. Similarly, grazing animals interact with vegetation dynamics, the presence of endangered or protected species may lead to habitat designation, and the presence of invasive species can disrupt local biodiversity and vegetation.
Such feedbacks are a fundamental characteristic of socio-ecological systems, and currently constitute a missing link in the assessment of processes and impacts of global change.

Author contributions
NWS, PEO and JMJT initially conceived the study, and then designed it together with CB and KW. NWS led the modelling, supported by CB, GB and SCFP. NWS analysed the simulation output and wrote the first draft of the paper together with CB. All authors contributed to further drafts. Included -10% chance that the long distance dispersal kernel will be used Figure 1: Workflow for the coupled RangeShifter (Bocedi et al. 2014) and CRAFTY model (Murray-Rust et al. 2014) illustrating loose coupling whereby models exchange output at each time step, driving the dynamics of the coupled subsystem in the following timestep.       Carrying capacity reduction factor (dK) VARIED: 0%; 50%

Movement model Dispersal kernel
Kernel type VARIED: negative exponential (kernel I only); double negative exponential (kernels I and II) Mean distance I (metres) 500 Mean distance II (metres) 1500 Probability of kernel II 10%

Dispersal -settlement
If the arrival cell is unsuitable Die Figure S1: