Coupled Societies are More Robust Against Collapse: A Hypothetical Look at Easter Island. Ecological Economics

Inspired by the challenges of environmental change and the resource limitations experienced by modern society, recent decades have seen an increased interest in understanding the collapse of past societies. Modelling efforts so far have focused on single, isolated societies, while multi-patch dynamical models representing networks of coupled socio-environmental systems have received limited attention. We propose a model of societal evolution that describes the dynamics of a population that harvests renewable resources and manufactures products that have positive effects on population growth. Collapse is driven by a critical transition that occurs when the rate of natural resource extraction passes beyond a certain point, for which we present numerical and analytical results. Applying the model to Easter Island gives a good ﬁt to the archaeological record. Subsequently, we investigate what effects emerge from the movement of people, goods, and resources between two societies that share the characteristics of Easter Island. We analyse how diffusive coupling and wealth-driven coupling change the population levels and their distribution across the two societies compared to non-interacting societies. We ﬁnd that the region of parameter space in which societies can stably survive in the long-term is signiﬁcantly enlarged when coupling occurs in both social and environmental variables.


Introduction
Modern society is pushing beyond the safe operating boundaries of its global environment (Rockström et al., 2009). Resource depletion, species extinction, deforestation and climate change are symptoms associated with passing these boundaries. A starting point for managing modern societies response to these problems is understanding the dynamics governing human-environment interactions, and studying past societies offers the possibility of gaining insight into these dynamics. There are many examples throughout history of societies that have collapsed: Mesopotamia, the Minoan and Mycenaean Civilisations, the Western Roman Empire, the Lowland Classic Maya and the Chacoans to name a few (Tainter, 1988). Understanding and quantifying the conditions which lead to collapse can provide a guide for present society concerning how to cope with resource constraints and environmental challenges, and avoid a potential collapse.
Societal collapse can be defined as a "rapid, significant loss of an established level of sociopolitical complexity" (Tainter, 1988, p. 4).
Under what conditions do human societies face a collapse or, in particular, a population crash? This question has been a recurrent issue in the history of mathematical and computational modelling. Arguably, the earliest study addressing it can be considered to be the work of Malthus (1798) which raised wide concerns over population growth and land availability. Much later, a more sophisticated model was developed in Limits to Growth (Meadows et al., 1972) that showed a decline in population levels in the second half of the 21st century. Both works were misunderstood and faced harsh criticism (Neurath, 1994), mostly because the quantitative reasoning and modelling they employed were unfamiliar at the time to large portions of the educated population, especially in the case of Limits to Growth (Bardi, 2011). In recent times, interest in the mathematical modelling of societal dynamics has revived (Anderies, 2000;Turchin, 2008), with several models quantifying historical cases of societies that have experienced collapse (Axtell et al., 2002;Brander and Taylor, 1998;Heckbert, 2013).
The focus of most recent models has been on the dynamics of single societies, most prominently Easter Island (Brander and Taylor, 1998;DaltonandCoats,2000;GoodandReuveny,2006;PezzeyandAnderies, 2003;Reuveny and Decker, 2000). While, like Easter Island, many cases of collapse occurred in isolation, the dynamics of multiple societies is no less important to historical knowledge and understanding. Given the highly interconnected nature of the modern world, studying the behaviour of a network with many coupled, interacting socio-environmental systems is of much interest for the present environmental debate (Rockström et al., 2009). Furthermore, it is likely that it can shed light on past cases of collapse. The movement of populations, along with the products they carry, have had a significant effect on social dynamics throughout human history (Manning and Trimmer, 2013). Historically, commerce and the establishment of trade routes proved crucial to the growth and development of many complex societies, in particular the Roman empire (Van der Leeuw and de Vries, 2003). The building of road networks and the development of exchange systems contributed to its growth (Dowdle, 1987). Subsequently, the migration of Germanic tribes had a large impact on the later stages of the empire (Halsall, 2007). Seafaring has been equally important to migration and trade, e.g., the Roman Empire with Egypt, but also in the society of the Maya (McKillop, 1996). The systematic movement of people and goods over large distances and the building of the necessary infrastructure (roads, bridges etc.) to facilitate transit and transport are characteristic signs of high levels of societal complexity (Tainter, 1988).
In this paper we present a dynamical systems model of societal development, which we apply to Easter Island. Then, we take a first step in the direction of quantifying the dynamics of networks of societies by starting with the simplest possible case, namely two coupled societies interacting by various forms of movement of people, goods or resources. To add historical context to our study, we again parametrise the model for Easter Island. By investigating the situation of interacting societies, we are essentially addressing the hypothetical situation: What if Easter Island was not an isolated island? What if people migrated and traded goods with an adjacent island? Would such interactions have prevented, alleviated or accelerated the collapse?
The paper is organised as follows. In Section 2 we review the recent literature on models of the long-term dynamics of societies. We find that the modelling methodology regarding dynamical systems is separated into two schools of thought, one focusing on economic principles underlying societal dynamics, while the other uses a more heuristic ecological approach to model building.
In Section 3 we propose a model aligned to the ecological methodology for the dynamics of a single society. In Subsection 3.1 we describe the model and apply it to reproduce the archaeological data of Flenley and Bahn (2003) about the collapse of Easter Island. Furthermore, a comparison is made between the results of the proposed model and the model of Brander and Taylor (1998), that also focuses on Easter Island. In Subsection 3.2 we analyse the critical transitions in the model when important parameters are changed. Appendix A contains a more detailed mathematical treatment.
Having validated the model for a single society, we proceed in Section 4 to investigate scenarios of two societies coupled together through exchanges of population, resources or goods under different conditions. We initially investigate simple diffusion in Subsection 4.1, and then in Subsection 4.2 we study targeted, wealthdriven migration, in which migration occurs from poor to rich societies driven by a difference in the wealth per capita. The results pinpoint when the coupling is beneficial or not for the sustainability of the entire system.

Related Literature
The earliest arguments regarding long-term societal sustainability or collapse that were amenable to quantitative modelling were in the work of Malthus (1798). Concerns regarding humanenvironment interactions started to re-emerge in the 1960s, when sustainability science was born, arguably through the publication of the work of Carson (1962). A work that sparked intense debate and interest in mathematical modelling was Limits to Growth (LTG) (Meadows et al., 1972), which proposed an aggregated world model formulated in terms of a set of equations describing the dynamics and time evolution of several key aspects of society. Such a set of equations constitutes a dynamical system that attempts to capture important feedbacks present in socio-environmental systems. Recent research by Turner (2008) has compared the trajectories from Meadows et al. (1972) with real data over a span of 30 years. The comparison shows a reasonably good fit despite the initial modelling effort not aiming to be predictive but only precautionary.
In recent times, modelling the long-term development of societies, in particular the possibility of their collapse, has received renewed attention. The paper that re-sparked interest in this topic was by Brander and Taylor (1998). By appealing to neo-classical utility maximisation arguments, along with a set of functional forms widely used in ecology (e.g., logistic growth) Brander and Taylor (1998) arrive at a set of predator-prey type equations used to describe the evolution of the human population and renewable resources on Easter Island. Following this work a stream of papers appeared that adopted the methodology of Brander and Taylor (1998), expanded on it or applied it to other cases. We can broadly categorise these models as "economic type models", meaning that they represent people as utility maximising, rational agents.
There is a second class of models that aim to capture societal development that we label "ecologically inspired models." In this case the choice of dynamical system is made heuristically to capture the observed real-world dynamics of the society, while respecting modelling principles of population biology (Turchin, 2003a), but rationality of individuals is not enforced. An early model following this approach was developed by Anderies (1998) to capture the social dynamics of the Tsembaga of New Guinea. The wider diversity of assumptions that underlie the ecological style of modelling means that the endeavour tends to lack the conceptual unity of economic models. The theoretical appeal of the more unified economic framework is understandable, with modelling efforts that were initially ecological in style (Anderies, 1998;Janssen et al., 2003), soon joining the economic camp (Anderies, 2000;Janssen and Scheffer, 2004). Nevertheless, the reliability of the rationality hypothesis has been questioned many times (Nell and Errouaki, 2013), including by Janssen and Scheffer (2004).
The underlying difference between the two modelling schools is in the theoretical framework from which they start the modelling process. The economic type models use a narrower set of assumptions, typically including utility maximisation as a driver of human behaviour along with a decision on the global institutional policy. The ecologically flavoured models are less restrictive in their theoretical underpinnings, more attentive to characteristic features of the society and try to account for emergent social phenomena that can contrast with or even contradict economic rationality, e.g. sunk-cost effects (Janssen et al., 2003) or war rituals (Anderies, 1998).
We proceed to give a short overview of the different strands of literature that deal with societal modelling, starting with the economic type models. Dalton and Coats (2000) extend the model of Brander and Taylor (1998) for Easter Island to account for possible institutional reforms and how these could effect feast-famine cycles. Reuveny and Decker (2000) investigate other possible solutions to the "Malthusian Trap" for Easter Island in the form of technological progress and population management. Pezzey and Anderies (2003), along with considerations of institutional adaptation, also add a resource subsistence requirement to preferences. Dalton et al. (2005) analyse Easter Island by taking into account the possibility of economic growth. Good and Reuveny (2006) look at limits to resource management institutions and conclude that even if the people of Easter Island had a complete assignment of property rights, implemented optimal resource management with an infinite horizon, or had a social planner to perform such a task, a boombust cycle would still occur. D' Alessandro (2007) investigates what happens when two different renewable resources are present which cannot recover once below a critical threshold. A comprehensive survey of these types of models is to be found in Nagase and Uehara (2011) and Reuveny (2012).
We now turn to the ecological strand of models which have arguably received less attention. Anderies (1998) modelled the social structure of the Tsembaga, which employ a simple swidden (slash-and-burn) agricultural system that also features domesticated animals, predominantly pigs which are part of is a periodic war ritual called the Kaiko. Other research attempts to quantify conceptual theories of societal dynamics, such as the work of Turchin (2003b) who, among several ideas, also develops models reflecting Ibn Klaldun's theory of collective solidarity ("asabiya"). Janssen et al. (2003) introduce a model that aims to illustrate what are called sunkcost effects, which refer to people taking into consideration prior investments when deciding what course of action to take. A group "may not suggest abandoning an earlier course of action because this might break the existing unanimity" (Janssen et al., 2003, p. 722). Turchin (2009) takes into account the demographic-structural theory proposed by Goldstone (1991) in developing a model that reproduces the dynastic cycles seen in many societies. Motesharrei et al. (2014) propose the Human and Nature Dynamics (HANDY) model for societal dynamics that attempts to capture the relation of economic inequality to collapse. The model suggests inequality can be a significant factor leading to collapse.
Regarding Easter Island, Erickson and Gowdy (2000) point out discrepancies between Brander and Taylor (1998) and the archaeological record for Easter Island reconstructed by Bahn and Flenley (1992) and propose the addition of a capital stock to account for continued population growth despite resources being depleted. While a more abrupt collapse is observed than compared to Brander and Taylor (1998), the overall trajectories for the population and resources levels do not match the quoted historical record. Basener and Ross (2004) build a two-dimensional model that aims to reproduce a discrete set of data points of the population levels, particularly around the population peak. While the model achieves the desired fit, the results show a near instantaneous collapse of the population after 1700 which is unrealistic. Bologna and Flores (2008) develop a very similar model to Basener and Ross (2004) but with a Lotka-Volterra term for resource extraction. While no instant collapse occurs, the results of Bologna and Flores (2008) show a much earlier peak in the population compared to data Flenley and Bahn (2003), similar to Brander and Taylor (1998). Basener et al. (2008) extend the model from Basener and Ross (2004) to account for rat infestation on Easter Island, while Brandt and Merico (2015) also consider rat infestation along with an epidemic component. With the higher dimensional model, Brandt and Merico (2015) investigate multiple scenarios regarding the collapse that reflect different theories regarding it. Plausible results are obtained under each scenario and none can be ruled out, which is not a surprising fact given the large number of parameters of the model and their uncertainty. We aim for a simpler treatment and propose a three dimensional model, that uses a similar parametrisation to Brander and Taylor (1998), provides analytic insights and has results that match very well with the time series by Flenley and Bahn (2003).
The examples in the above categories develop dynamical systems that describe a single, isolated society. Follow-up studies by Quirin et al. (1977) and Mesarovic and Pestel (1974) extended the model in LTG by disaggregating the world into two and ten regions, respectively. The finer resolution allows for an analysis of societal development that moves past the homogeneity assumption that went into building the model behind LTG. For example, both Quirin et al. (1977) and Mesarovic and Pestel (1974) take into consideration the differences in internal development between the various regions, the imports and exports among them and variation in industrial activity.
While the work of Mesarovic and Pestel (1974) increased spatial resolution, it decreased the temporal horizon to 50 years and shows increasing trends in population and resource usage and no global collapse, much like the results of Meadows et al. (1972) up to the middle of the 21st century. Quirin et al.'s (1977) standard run also shows steadily increasing production of services, industrial output and population but for the next 300 years, completely disagreeing with Meadows et al. (1972). A significant role in this outcome is likely played by the annual rate of technological growth, which Quirin et al. (1977) assumed to be at least 5% and has a beneficial effect on all aspects of the model.
Instead of using a model of high complexity such as the one in LTG, here we opt for lower dimensional models and focus on systematically understanding their dynamics. Research on (relatively) simple models of coupled dynamical oscillators has a rich history in theoretical ecology, the closest related literature to our purposes being the Lotka-Volterra (LV) predator-prey models that explore diffusion between two patches. Early research on the topic was pursued by Levin (1974), who finds predator and prey can co-exist provided the migration rate is sufficiently low, while higher rates make the system act like a single patch and co-existence is no longer possible. Nisbet et al. (1992) and Jansen (1995) find that, for identical patches linked together by dispersal at a constant rate, the only possible equilibrium is symmetric with equal densities on each patch, which is a general result highlighted in a review by Briggs and Hoopes (2004). On two heterogeneous patches with LV predation, Murdoch et al. (1992) find that different prey birth rates have a stabilising effect, while Jansen (1995) concludes that temporal variability in predator death rates decreases the amplitude of fluctuations in the system when it is perturbed by noise.
While the ecological literature provides us with useful mathematical tools, there are nevertheless specific features that predator-prey models do not capture when applied to societal dynamics. Examples of this include the building and use of products derived from natural resources (such as tools, food reserves and infrastructure) or specific incentives regarding migration, like economic prosperity, both of which are key aspects that we are modelling. However, methods from population biology modelling have been used to asses the impact of human population growth and migration on the carrying capacity of the environment, and to understand the feedbacks associated with these processes. In particular, the Prehistoric U.S. Southwest has been researched from this perspective, with early work on the problem by Zubrow (1971), and a more recent study by Anderies and Hegmon (2011) who propose a heterogenous, threepatch model that captures population and resource dynamics, with net migration occurring towards a region if its resource stock is sufficiently large compared to neighbouring patches. Anderies and Hegmon (2011) find that the equilibrium states of regions with and without migration are very similar.
In the present work we present a thought experiment wherein two regions similar to Easter Island interact through diffusion or through targeted migration due to differences in wealth per capita. The regions not only exchange people, but also resources and manufactured goods can interchange between patches. The coupling of the patches significantly changes the equilibrium states obtained in the regions, a result which is in opposition to the findings of Anderies and Hegmon (2011).
Societal dynamics has also been explored with agent based models. Representative examples for the work done in this area are Axtell et al. (2002), who develop a model of the Anasazi in north-eastern Arizona (U.S.), Heckbert (2013), who focuses on the Maya civilisation and Turchin et al. (2013), who aim to understand how intense competition, in particular warfare, contributed to the emergence of large-scale complex societies. Agent based models, like the examples mentioned, incorporate detailed features such as spatially extended terrain and multiple forms of social interaction, and account for heterogeneity in agent behaviour and in the environment. The focus of this paper is, instead, to obtain a description of the time evolution of aggregate measures of societal macro-features, such as the population and resource consumption. Following previous research in the area (Basener and Ross, 2004;Bologna and Flores, 2008;Brander and Taylor, 1998;Brandt and Merico, 2015;Erickson and Gowdy, 2000), dynamical systems are an adequate tool for our purposes, which we use to formulate and specify our model.

Model Specification
Here, we opt for an ecological type of model for the dynamics of a single society that we apply to Easter Island. Using, where possible, the same parameter choices as Brander and Taylor (1998), we compare the solution of the model to the historical data provided by Flenley and Bahn (2003), who present an exhaustive review of the archaeological research of Easter Island. For comparison with an economic approach, we also show the model of Brander and Taylor (1998) and its solution. Subsequent papers, that extended the treatment of Brander and Taylor (1998), did not attempt to match the historical data, hence, our focus for comparison is on the original model.
By matching model output to real data we, at least partially, validate the proposed model. This means that by reproducing the observed historical trends, we can consider the relationships that the model describes to be plausible mechanisms at play within the socioenvironmental system. Hence, the model moves beyond a thought experiment, and makes a better connection to reality which lends validity to later extensions of the model for coupled societies. However, it is beyond the scope of this paper to explain the collapse of Easter Island or quantify all relevant factors that led to it, such as rats (Hunt, 2007) or European contact (Stevenson et al., 2015).
We propose the following model for the dynamics of a society: We regard the x variable as a stock of population, y as a stock of renewable natural resources and z as the cumulated goods produced from extracting natural resources, which we also refer to as wealth. In more concrete terms we could think of y as crops, wood, animal stock, etc., and z as products derived from these, like food stocks, shelter, tools, etc.; these are characteristic elements of past agricultural societies.
The dynamics governing the population change is consistent with the early phases described by demographic transition theory: in a developing society (which has relatively low, but growing levels of wealth per capita) the death rates drop quickly due to improvements in food supply and material conditions (Kirk, 1996;Landry, 1934). Data for modern societies indicates that infant death rates have fallen significantly over time (Chesnais, 2001), sometimes exponentially. We do not include a mechanism for the reduction of birth rates because the maximum net birth rate in the model is 0.2% per year, comparable to net growth rates in the High Middle Ages in Europe (Russell, 1972), which we judge to already be quite low.
The demographic mechanism is implemented as follows: if the wealth per capita z/x is high compared to a threshold q then the population grows exponentially at the (maximum) rate b, while if Table 1 Parameters for model (1) and the model of Brander and Taylor (1998)  z/x is close to 0 then the population decreases exponentially at the (minimum) rate b − d, where d is the maximum raw death rate. The presence of material wealth lowers the death rate, hence it has a positive impact on net growth rates, consistent with the first stages hypothesised by demographic transition theory.
As has been common in many models of societal dynamics, we consider natural resources to recover logistically at a rate r, with a maximum capacity of K. Another classical feature is that the resource extraction is governed by a predator-prey term axy, with people acting as predators and resources as prey. The extraction rate a is the fraction of the resource base that can be extracted by one person over a year. Thus, the equation for the rate of change of resources is analogous to that used by Brander and Taylor (1998).
Material goods, i.e. the wealth z, are produced at the same rate as the resources are extracted, which reflects conservation of matter. The consumption of goods is proportional to the population and the rate of consumption per capita is equal to the subsistence requirement s when wealth per capita is high. We expect that, if the average wealth per capita decreases, then the average consumption per capita also decreases. Hence, the consumption per capita tends to zero when z/x is low, which is similar to the dynamics proposed by Motesharrei et al. (2014). To account for wear-and-tear, manufactured products decay exponentially at a rate c, a feature previously considered by Erickson and Gowdy (2000). A summary of the parameters and their values is found in Table 1.
We show that by solving the equations of model (1) it is possible to reproduce to high accuracy the archaeological record for Easter Island as presented in Flenley and Bahn (2003). To achieve this, the choice of parameters in Eq. (1) is made to coincide with the values (and notation) used in model (2) proposed by Brander and Taylor (1998): where L is the population and S is the renewable resource stock. We set most of the parameters in Eq.
(1) to have the same values as the choices Brander and Taylor (1998) made for model (2). The maximum birth rate is b = 0.002, while the maximum death rate is d = 0.012, so that b − d = −0.01 as in (Brander and Taylor, 1998). The regeneration rate (per year) of resources is r = 0.004, the maximum resource level is K = 12, 000 and the extraction rate is a = 10 −6 (per person, per year). Some parameter choices for model (1) differ from model (2) of Brander and Taylor (1998). We choose the subsistence requirement per person per year to be s = 0.004. Despite being mentioned, no explicit value is quoted by Brander and Taylor (1998) for the subsistence requirement. The largest difference between Eqs. (1) and (2) is how the resource extraction is modelled and how this impacts birth rates.
For Brander and Taylor (1998), the parameter b is the fraction the population involved in resource extraction, so ab is the effective extraction rate per capita in model (2). The parameter a was assigned the value 10 −6 so that at maximum resource levels a "household could provide its subsistence requirements in about 20 percent of available labour time . . . A value of 0.4 for b is probably in the reasonable range" (Brander and Taylor, 1998, p. 128). In model (1) we eliminate the redundancy of the extra b parameter and take the effective extraction rate per capita to be a = 10 −6 , such that at maximum resource levels a person can meet his subsistence requirement in 33% of available time. While debatable, this value for the extraction rate allows the solution of Eq. (1) to match the archaeological record, see Fig. 1.
The parameter 0 = 4 scales the effect of resource consumption on birth rates in model (2). Resources do not directly impact birth rates in model (1) so, there is no 0 parameter in our case. We made the assumption in Eq. (1) that wealth per capita affects demographic rates and the parameter q = 0.1 sets the scale at which the impact of wealth cumulation on population growth is significant, corresponding to q/s = 25 years of material requirements (food reserves, tools etc.). There is no equivalent parameter for the decay rate c in model (2) so, we take c = 0 in model (1) if we are considering Easter Island. By choosing c = 0, we are assuming that wealth decreases mainly through human consumption. Fig. 1 (a) shows the historical evolution of population and resource levels on Easter Island according to Flenley and Bahn (2003). For initial conditions of 1100 people, 12,000 units of natural stocks and no wealth, the output of model (1) is shown in Fig. 1 (b). The initial conditions for model (2) are set to 40 people, while resources are at 12,000 units, with trajectories in Fig. 1 (c).
The population peak in Fig. 1 (c) occurs several centuries earlier than in the data of Fig. 1 (a), as noted by Erickson and Gowdy (2000). Also, the decline in population and resource levels in Fig. 1 (c) is much less abrupt than Flenley and Bahn (2003) show in Fig. 1 (a). Furthermore, the initial population value of 40 is much lower than estimates of the minimum viable population sizes, which are placed between 100s and 10,000s for humans (Smith, 2001). A higher initial population leads to an even earlier population peak in model (2).
Model (1) does not suffer from these shortcomings and its output in Fig. 1 (b) fits Fig. 1 (a) visibly better, both in population and aggregate resource trajectories, while using almost identical parameter estimates to Brander and Taylor (1998). The population peak, its timing, the collapse profile and the significant decline in resources, as indicated in Fig. 1 (a) are accurately reproduced in Fig. 1 (b). Thus, the results in Fig. 1 partially validate model (1). The resource level in Fig. 1 (b) should be seen as an aggregated mean of renewable stocks (trees, livestock etc.) which are not all included in Fig. 1 (a). This explains why the minimum resource level in Fig. 1 (b), which occurs close to the time when the population reaches its maximum level, does not precisely match the resource minima in Fig. 1 (a).
Why do we obtain a better fit using model (1)? Part of the reason lies in the delay induced by the additional equation describing the accumulated amount of manufactured products (wealth), which leads to a longer feedback loop within system (1) compared to other models.
The fit between the output of model (1) and the historical record shows us that the structure of the model and the feedbacks captured by it are viable choices to describe the dynamics of a simple agricultural society. Moreover, by fitting the model to a real worldcase, we obtain reasonable parameter estimates, which can inform further modelling efforts and make them more realistic. Given this, the model can serve as a starting point to explore new dynamics, like that of multiple interacting societies. Nevertheless, it is worth pointing out that the results do not fully explain or quantify the scope of the collapse of Easter Island, as many more factors were at play in  (Flenley and Bahn, 2003, p. 201). We see the population reaches a peak of approximately 10,000 people around the year 1600 CE. The resources were continually declining as the population was increasing. (b) The population and resource levels for Easter Island as determined by model (1) from the year 400 CE to 2000 CE assuming no outside influence. The predictions are valid up to the late 19th century when significant interaction between Easter Island and the outside world started taking place. We see a good fit between the output of Eq. (1) and the real data.
(c) The solution of the model by Brander and Taylor (1998) shown for comparison. The timing of the population peak as well as the overall shape of the model output is inconsistent with the data in Fig. 1 (a). The initial population value is below the lower bound for the minimum viable population size for humans. the collapse and model (1) does not incorporate all the complexities of reality.

Critical Transitions
In this section we focus on the critical transitions that model (1) exhibits. Fig. 2 shows the bifurcation diagram for system (1)  varying the extraction rate a. The critical transition points are highlighted through vertical lines that separate the different regimes. The different regimes correspond to different long-term outcomes for our model society.
The mathematical analysis of the equilibrium points and transitions undergone by the system (1) is provided in Appendix A. The model has a single interior equilibrium point that we denote by E = (x e , y e , z e ) where: The first transition occurs at a = (s(1 − b/d) + qclogd/b)/K, which we take as a reference value for the other transitions. When the extraction rate a is varied beyond a it leads to a stable steady state (at the point E) or to oscillations, see Fig. 2. This indicates the presence of a critical transition (a Hopf bifurcation) in the system at a certain threshold value of the extraction rate, which we denote as a c . That value separates regimes with long-term stable outcomes from boom-bust cycles. The precise determination of the transition point is of interest, as it dictates the eventual state of the system. The structure of system (1) allows for the exact analytic determination of the critical value a c at which the Hopf bifurcation occurs: where .
In case of a < a the equilibrium population is zero, which means that the extraction rate is too low for any human population to meet their subsistence requirements. For a < a < a c there exists an attractive fixed point E, as the single branch of the bifurcation diagram in Fig. 2 indicates. This regime corresponds to a sustainable society in the sense that it equilibrates to a long term steady state. For a > a c the system oscillates and we can interpret this as a sequence of boom-bust cycles the society goes through.
For very large extraction rates, exceeding a threshold we call a h , the natural resources are depleted quickly and they reach low enough levels for long enough that the population dies off and all wealth disappears. Without population pressure, natural resources recover and the system stabilises at a fixed point with no human population or wealth, but with maximum amount of natural resources.
In practice, if a > a c , the larger amplitude cycles would likely manifest as a collapse due to the low population values that are reached at the lower bound of the periodic orbit.
Relation (4) allows us to determine the behaviour of the critical threshold a c when varying other parameters. In particular, we notice that the relative extraction rate a c /a , at which the transition to the boom-bust regime occurs, has no dependence on the maximum resource level K. The dependence of a c /a on the value of r, the regeneration rate of natural resources, is shown in Fig. 3 (a), while the dependence on s, the subsistence requirement, is presented in Fig. 3 (b). As we see, the higher the regeneration rate r the higher the critical ratio a c /a at which the transition to the unsustainable regime takes place. A similar trend is seen for the subsistence rate s. Also, with increasing values of the regeneration rate r, the threshold a h /a is shifted upwards, whereas for the case of subsistence requirement s it is roughly constant.
With an increasing regeneration rate r we would expect an upward shift in the points of critical transition as the natural resources replenish faster. An increase in the subsistence requirement means goods are consumed at a higher rate which tends to lower the cumulated wealth, and hence the wealth per capita which leads to lower birth rates. We would thus expect the system to reach a steady state over a wider range of extraction rates, which is what Fig. 3 (b) confirms.
In the case of Easter Island we have that a 0 /a = 3.6 which is larger than the critical ratio a c /a = 2.8 so, the system is in a boombust regime that in real terms translates into a collapse. Thus, the critical value a c of the extraction rate serves to separate the sustainable regime (that leads to a steady state) from unsustainable regimes (where collapse can occur).  Table 1.
We can give model (1) and its behaviour a more concrete interpretation by identifying the feedbacks at play. If the extraction rate is above the critical value than the feedback loop dominance in the model changes between two competing mechanisms: resource extraction and wealth consumption. Initially, the resource extraction is the dominant feedback (an increasing population leads to a decrease in resources and an increase in wealth, which implies an increase in wealth per capita and, thus, a further increase in the population). Once the resources are depleted, the cumulated wealth reserves are no longer produced but only consumed, which leads to a decrease in wealth per capita, which implies a decrease in the population. A smaller population means a slower rate of wealth consumption. As can be seen, resource extraction is part of reinforcing feedback loop, while wealth consumption is in a negative feedback loop. The feedback loops balance out when a < a c , while if a > a c the feedbacks alternate in strength.

Diffusive Coupling
Having established and understood a model that reproduces the historical development of Easter Island, we now proceed with an analysis of two coupled socio-environmental systems. This is a first step at trying to quantify and understand the dynamics of networks of multiple, interacting societies, which is a situation that characterises the modern, industrial world but was also typical of many ancient societies. We aim at obtaining general insights on what new features the coupling of two socio-environmental systems gives rise to.
Human settlements and their surrounding environment, e.g. islands, form a natural partition of the overall socio-environmental landscape. Flows between different regions, such as migration or trade, provide additional sources (or sinks) for the respective stocks within regions. Accounting for these flows adds an important degree of realism to any model that tries to quantify inter-societal interactions.
There are many different types of possible interactions, and hence flows, that can occur between two systems. We start by considering possibly the simplest of such interactions, namely we model flows between systems through the mechanism of diffusion. If we assume people, resources and goods move or are moved randomly and unguided between different regions then a coarse description of such flows is given by a diffusion term in the equations. A further motivation for our choice is that diffusion is one of the most common types of coupling employed in the study of coupled dynamical systems (Pikovsky et al., 2001). Also, diffusion terms are widely used in spatial ecology to represent population dispersal from one region to another (Briggs and Hoopes, 2004). We do not model the geographical structure behind the flow of people, resource or goods. Here we are only considering how the average flows between two regions (or land masses, islands) impact the stocks levels within the different systems.
More specifically, we consider two instances of the system (1) which we couple by adding diffusion terms to all variables. The equations for the first socio-environmental system (referred to as S1) are given by: where the first terms are from model (1) and the additional terms represent the diffusive coupling, with s x , s y , s z acting as coupling constants. The equations for x 2 , y 2 , z 2 of the second system (called S2) are analogous to Eq. (6). For concreteness and mathematical convenience, the two societies share all the same parameter values except for the extraction rate and initial conditions. S1 has the same initial conditions (x 0 , y 0 , z 0 ) as Easter Island, while S2 has half those values, namely (x 0 /2, y 0 /2, z 0 /2). Hence, we can think of the second society as a sister island to Easter Island, which shares most of its features but starts off with fewer people and resources.
The diffusion in the population x can be interpreted as migration.
If x 1 > x 2 there will be a net influx of people from S1 into S2. A similar mechanism holds for the other diffusion terms. The diffusion of wealth z along with migration can be seen as people moving with their share of goods. A value of 0.1 for a diffusive coupling constant is quoted with respect to one year, just like the rest of the parameters in model (1). So, if s x = 0.1 then 10% of population difference can move between the societies (or islands) in the span of one year, e.g., if one island has 10,000 people and the other 5000, this means roughly 500 people moving per year between the two landmasses. We consider natural resources to be less mobile and when they do diffuse, only 1% of their difference is allowed to move between the regions in a year. . 4. Graphs for two societies coupled via simple diffusion. Horizontal and vertical (white) lines at the critical ratios divide panels into quadrants, labelled from I to IV. Top: The relative population change in population for system 1 (S1) when: (a) uncoupled from system 2 (S2), with (sx, sy, sz) = (0, 0, 0), or coupled via: (b) migration, with (sx, sy, sz) = (0.1, 0, 0), (c) migration and resource diffusion, with (sx, sy, sz) = (0.1, 0.01, 0) and (d) migration with wealth and resource diffusion, with (sx, sy, sz) = (0.1, 0.01, 0.1). The gray dotted curves indicate the bifurcation boundary determined using the first-order approximation of the interior equilibrium point for the coupled system, see Appendix B. Middle (e)-(h): The average population for S1 in the same regimes as above. Bottom (i)-(l): The ratio between the average population in S1, S2 and overall system in the coupled and uncoupled regimes along the black dotted line where a 2 /a = 1.8. The mesh size in the numerical simulations is Da/a = 0.1.
As we saw in the case of one society, the factor determining the long-term behaviour (steady state or collapse) is the critical value of the extraction rate a c . We are interested in investigating if coupling societies can expand the volume of the parameter space that leads to sustainable outcomes. Hence, we will focus on the equilibrium states of the system rather than the transient dynamics.
To help us identify the bifurcation boundary we first quantify the decline in the population by looking at the relative change in the population over the long term, i.e., the relative distance between the maximum and the minimum of the oscillations with the amplitude measured in the population variable.
Thus, we define: where lim inf t→∞ x(t) is the smallest value the population reaches in the long term and lim sup t→∞ x(t) is the largest. If a system reaches a steady state (sustainable outcome) then the relative change is 0%, with no decline in the population. On the other hand, if there are large amplitude oscillations (collapse) the relative change in the population is close to 100%. The use of the measure (7) allows us to represent the 3-dimensional bifurcation diagram for the population of S1 on a 2-dimensional diagram, see Fig. 4. Other relevant measures are the long-term average population l of a system: and the average population relative to the uncoupled regime, which we write as g: An index indicates what society we are referring to, e.g. d 1 for S1, while no index refers to both societies. We are interested in the long-term behaviour of the system for relative extraction rates above and below the critical ratio a c /a = 2.8 for a single society. To explore the parameter space we take the relative extraction rates of the two societies a 1 /a and a 2 /a , and vary them in the range 1 to 5. No population exists for an isolated society if a/a < 1, which sets a natural lower bound for the interval of exploration. The choice of the upper bound forms an almost symmetrical interval around the critical ratio. Fig. 4 (a)-(d) shows the relative change in the population of S1 when diffusion occurs in progressively more dimensions of the system. Similarly, Fig. 4 (e)-(h) shows the long-term average population of S1. The corresponding figures for S2 are given by reflecting the graphs in the two top rows of Fig. 4 with respect to the second diagonal.
Horizontal and vertical (white) lines divide the panels into 4 quadrants, such that both relative extraction rates are above the critical ratio in quadrant II and below it in quadrant III. Only one of the extraction rates is above the critical ratio in quadrants I (S2) and IV (S1). Fig. 4 (i)-(l) indicates the ratio between the average population size of S1, S2 and the overall system to its value in the uncoupled regime along the dotted black line, where a 2 /a = 1.8. Fig. 4 (a), (e) and (i) shows the relative population change, the average population and the relative population (along the black dotted line) in the uncoupled case and can be taken as a reference to compare the results in the coupled regime. We can see that Fig. 4 (a) reflects the critical transition seen in Fig. 2 (a) that separates the steady state and oscillatory regimes. The gray dotted line in Fig. 4 (a)-(d) represents the bifurcation boundary computed using an approximation of the stable equilibrium point of the system, see Appendix B. The determination of the boundary becomes more accurate due to better mixing of the system when more couplings are present, as can be seen in Fig. 4 (d) that includes diffusion in all the variables.
A steady state is a sustainable outcome, hence we refer to the region of parameter space where the relative population change is 0% as the sustainability region. By contrast, the region with a 100% change in the population is the collapse region. In all cases of Fig. 4 (a)-(d) quadrant III lies in the sustainability region and quadrant II in the collapse region. So, when both societies have a low extraction rate the overall system reaches a steady state, whereas if both societies intensively extract resources we have a case of mutual collapse. Fig. 4 (b)-(d) is symmetric along the second diagonal, so we only need to discuss quadrant IV. In contrast to Fig. 4 (a) all coupled regimes show the sustainability region extending to quadrant IV, where a 1 /a is above the critical threshold. In Fig. 4 (b) and (d) the sustainability area is similar, occurring below the main diagonal but in Fig. 4 (c) it extends to most of quadrant IV.
When a 1 = a 2 , which occurs along the second diagonal in Fig. 4 (b)-(d), the societies share the exact same parameters and they follow identical trajectories in the long run. Thus, they synchronise completely. The populations of the two societies continue to synchronise throughout most of the parameter space that we consider, which is consistent with the symmetrical stable equilibria seen in predator-prey patch models with diffusion (Briggs and Hoopes, 2004). For the sustainability region, synchronisation is illustrated by the population levels in Fig. 5 (c) and (d), (e) and (f), (g) and (h) that match closely. For the collapse region, synchronisation implies that the collapse occurs in both societies simultaneously.
The region with the highest average population in Fig. 4 (e)-(h) is shown in white and, in all cases, it occurs within the sustainability region. The average population then gradually decreases in the direction of the collapse region. In case of migration in Fig. 4 (f), the region with the highest average population does not change significantly from the uncoupled case of Fig. 4 (e). By adding couplings in more dimensions the area increases noticeably, as seen in Fig. 4 (g), (h). Fig. 4 (j)-(l) shows that, on the black dotted line where a 2 /a = 1.8, the long-term average population of each society in the coupled regime is equal to its population in the uncoupled regime, provided a 1 /a is smaller than the critical ratio a c /a = 2.8. For values above the critical ratio, the average population of S1 is higher in the coupled regime than when isolated, while for S2 it is lower. Along the black dotted line, the population of the overall system of coupled societies is greater or equal to the sum of the populations of the isolated societies, with the only exception occurring in Fig. 4 (j) when a 1 /a 5 and diffusion takes places only in the population variable. Repeated numerical simulations show that the results in Fig. 4 (a)-(k) are robust over changes in initial values and coupling constants. If the coupling constant for product diffusion is greater than that for population diffusion, namely s z > s x , the results are similar to those in Fig. 4 (d), (h).
We next look at a specific scenario: we take a 1 /a = 3.6 which corresponds to the case of Easter Island and a 2 /a = 1.8. The results are shown in Fig. 5 (a)-(h) for the different couplings. We see that in the uncoupled case in Fig. 5 (a), (b) S1 collapses and S2 reaches a steady state.
In the coupled regimes both systems stabilise so, Easter Island could potentially have been saved provided an almost identical sister island was present where the population extracted resources at a smaller rate. The most rapid approach to equilibrium occurs in Fig. 5 (e), (f) where migration and resource diffusion take place but no wealth is transported between the two societies. In the case of population diffusion in Fig. 5 (c), (d) and full diffusion in Fig. 5 (g), (h) there is a long-term, oscillatory transient until equilibrium is reached.
The increased sustainability regions and the new areas of higher population that occur in the coupled regime compared to the uncoupled one, allows us to conclude that diffusive coupling makes societies more robust against collapse. By this we mean a larger area of the parameter space can be explored without the risk of collapse or significant loss of population. We can broadly interpret these results as follows: the diffusion of people alleviates the demographic pressure from the more populous society, while the diffusion of goods lowers the wealth per capita and hence the birth rate in the society with higher overall wealth. This allows both societies to reach a steady state as long as one of the societies has a low enough extraction rate of natural resources.

Wealth-Driven Coupling
In the previous section we explored what happens when two societies are coupled in a simple, diffusive manner. We can explore alternative ways to couple societies and this is the goal here. What if people decided to move from poor to rich societies, i.e., from low to high wealth per capita?
We can write a dynamical system that captures this possibility as follows: Fig. 6. The graph of the step function h in terms of the difference in wealth per capita. As a reference, we indicate the zero point with a vertical line.
The function h is a smoothed out step function shifted to the right by one unit, as seen in Fig. 6. If z 1 /x 1 ≤ z 2 /x 2 then there will be no movement of people from S2 to S1. The greater z 1 /x 1 is relative to z 2 /x 2 the greater the influx to S1, which is shown in the monotone increase of the function h. So, if the wealth per capita in S1 is greater than in S2, then there will be an influx of people from S2 into S1 which corresponds to the intuition that people will move to the society with higher wealth per capita. Fig. 7 shows the same analysis as performed previously in Fig. 4 but now for the case of system (10).
For convenience, Fig. 7 (a) and (e) again shows the relative population change and average population of S1 in uncoupled regime, which is to be taken as a reference case. We see that in all cases of Fig. 7 (b)-(d) quadrants I to III are predominantly in the sustainability region, whereas quadrant II is mostly in the collapse region. Again, if a 1 = a 2 the two societies synchronise completely but, in contrast to the case of diffusive coupling, synchronisation no longer extends beyond the second diagonal, as is shown by the asymmetry in the average population levels in Fig. 7 The average long-term population of S1 is increased in the new sustainability regions of the population diffusion case, as Fig. 7 (f) shows. Compared to the other graphs, we notice some new features present in Fig. 7 (g) and (h) regarding the long-term average population of S1. The novelty of these scenarios consists in the very low (almost zero) population in S1 above the second diagonal, and the clustering of almost the entire population of both societies into S1 when 2 < a 1 /a < 3 and a 2 /a is close to 1. Both cases We can see that resource diffusion is important because in the region of parameter space with the highest population of Fig. 7 (g), (h) essentially most of the population moves into S1 and lowers the level of natural sources. In S2 on the other hand, the population is low while the resource level is high and this ensures a high influx of resources into S1, allowing it to maintain a higher population. Fig. 7 (i)-(l) shows the ratio of the average population in S1, S2 and the overall system to its population in the uncoupled regime for . 7. Graphs for two societies coupled via wealth-driven diffusion. Horizontal and vertical (white) lines at the critical ratios divide panels into quadrants, labelled from I to IV. Top: The relative population change in population for system 1 (S1) when: (a) uncoupled from system 2 (S2), with (sx, sy, sz) = (0, 0, 0), or coupled via: (b) migration, with (sx, sy, sz) = (0.1, 0, 0), (c) migration and resource diffusion, with (sx, sy, sz) = (0.1, 0.01, 0) and (d) migration with wealth and resource diffusion, with (sx, sy, sz) = (0.1, 0.01, 0.1). Middle (e)-(h): The average population for S1 in the same regimes as above. Bottom (i)-(l): The ratio between the average population of S1, S2 and the overall system in the coupled and uncoupled regimes along the black dotted line where a 2 /a = 1.8. The mesh size in the numerical simulations is Da/a = 0.1.
the subset of the parameter space that lies along the black dotted line in Fig. 7 (e)-(h). For diffusion occurring only between populations, the average population of S1 is equal to that of the uncoupled case if its relative extraction rate is below the critical ratio a c /a = 2.8, after which the population increases, as Fig. 7 (j) shows. For S2, the population is fairly similar to the uncoupled regime. In Fig. 7 (k), (l) at low extraction rates the population of S1 is 0, but then steadily increases at higher extraction rates. In S2 the situation is reversed, with the population high at low extraction rates (for S1) and then decreasing to 0 at higher rates. This contrast is to be expected from the asymmetry in Fig. 7 (g), (h). The overall average population of the system in the coupled regimes in Fig. 7 (j)-(l) matches the average population in the uncoupled case.
We again look at the specific scenario when S1 has the extraction rate corresponding to Easter Island, namely a 1 /a = 3.6 and S2 has half that value, a 2 /a = 1.8. The results are presented in Fig. 8 and we see that the coupling again stabilises both systems. In the case of migration in Fig. 8 (c), (d) both societies reach the same population level in the steady state. In the other coupled regimes of Fig. 8 (e)-(h) the population of S2 is 0, having moved completely into S1.
In all coupled cases we see an increased sustainability region in Fig. 7, noticeably larger than the corresponding simple diffusion cases of Fig. 4. Why is this the case? Whenever one of the societies has a higher wealth per capita than the other, people will migrate to it. The influx of people increases demographic pressure in the richer society but lowers its wealth per capita. The lower wealth per capita leads to a decreased population growth rate that helps stabilise the system. This different operating principle makes wealth-driven diffusion more effective at generating sustainable outcomes than simple diffusion, while maintaining similar overall levels of the population. Similarly, we can conclude that the coupling via wealth-driven diffusion makes the overall system of societies more robust against collapse.

Conclusion
In this paper we proposed and analysed a model to describe the simultaneous dynamics of the population, resources and manufactured goods of a society. We can match most of the parameters in the model with a well-known model of Easter Island by Brander and Taylor (1998) and by using similar values we obtain a faithful reproduction of the archaeological record provided by Flenley and Bahn (2003). The model captures crucial feedbacks present in most ancient agrarian societies, making it more general than the case of Easter Island.
We found the equilibrium points of the system along with the critical transitions (bifurcations) it undergoes when the rate of resource extraction is varied. At values below a particular extraction rate a c the system reaches a steady state (sustainable outcome), whereas for value above it the system displays oscillatory behaviour (unsustainable outcome). Large amplitude oscillations can be identified with collapse, and this is the regime that Easter Island is found to be in. It is important to mention that for ecological systems, fold bifurcations provide a general mechanism by which ecosystem transitions can be understood and modelled (Scheffer et al., 2001). Given the prevalence of persistent oscillations in data and models of societal dynamics (Anderies, 1998;Turchin, 2009;Turchin and Nefedov, 2009), we can hypothesise that Hopf bifurcations, as the one illustrated by the model we propose, could serve a similar role when modelling critical transitions, such as collapse.
We then investigated the case of two societies coupled through simple diffusion or wealth-driven diffusion. We set up experiments with two societies parametrised to match Easter Island, except for the extraction rates which we systematically varied. In contrast to the results of Anderies and Hegmon (2011), we generally found that population migration, along with resource and wealth diffusion, makes a significant change to the equilibrium states the subsystems reach compared to uncoupled, non-interacting societies.
In particular, coupling leads to new or extended regions of parameter space where sustainable outcomes occur. This means that societies that would have undergone large amplitude fluctuations and collapsed if they were isolated, can now reach a stable steady state, provided they are coupled with other societies with less resource intensive practices. The differences in resource extraction rates between the two societies help stabilise the overall system, which echoes findings that heterogeneity in two patch predator-prey models is stabilising (Briggs and Hoopes, 2004;Murdoch et al., 1992;Murdoch and Oaten, 1975). Furthermore, the regions of parameter space where a given society achieved the highest average population increased. The more couplings are added (in population, resource etc.) the stronger the effect.
The present world is in a similar situation where, because of competing power centres of comparable strength, any collapse that could occur is more likely to be simultaneous and global (Tainter, 1988). These results indicate that a coupled network of societies could prove more robust against collapse provided that at least some societies maintain lower extraction rates of natural resources to dampen the oscillations in the rest of the system. Precisely quantifying the regime in which overall system sustainability is achieved is of obvious importance and is the future focus of our research efforts. In all cases we see a reasonable match between analytic (dashed lines) and numerical results (dots). The results are quoted with respect to the equilibrium values for an isolated society with extraction rate 2a , which is taken as a reference.