Multiple Resource Use Strategies and Resilience of a Socio-Ecosystem in a Natural Protected Area in the Yucatan Peninsula, Mexico

As the world faces unprecedented ecological and social changes, there is a need to better understand the complex dynamics of social-ecological systems (SES) and the mechanisms that underlie their resilience. In Mexico, Natural Protected Areas (NPAs) constitute complex SES as they are generally established on territories that different peoples have historically inhabited and managed. They generally manage their resources following a multiple use strategy (MUS), which involves local traditional agricultural practices and has been proposed as a resilience-enhancing mechanism. In this paper we study the MUS as practiced by the Yucatec Maya communities that inhabit the NPA Otoch Ma'ax Yetel Kooh and its buffer zone in the Yucatan Peninsula, Mexico. Due to the restrictions imposed by the decree of the reserve and the growth of tourism in the region, some of these communities have started to abandon the MUS and specialize on tourism-related activities. To study the consequences of these changes and to better understand the mechanisms by which the MUS may enhance the resilience of this SES, we built an evidence-based dynamical computational model that allows us to explore different virtual scenarios. The model, through the incorporation of agent-based and boolean network modeling, explores the interaction between the forest, the monkey population and some productive activities done by the households (milpa agriculture, ecotourism, agriculture, charcoal production). We calibrated the model, explored its sensibility, compared it with empirical data and simulated different management scenarios. Our results support the hypothesis that the MUS enhances the resilience of this SES in terms of income and food availability, as it increases the system's response diversity and functional redundancy, thus reducing income variability and increasing the resistance to natural and anthropogenic disturbances. We also identify other possible mechanisms related with the MUS that provide a more nuanced understanding of the resilience of this SES. Our study, in addition to highlighting the importance of local management practices for resilience, also puts forward a novel integration of diverse mathematical formalisms and illustrates how computational modeling and a systems perspective are effective means of integrating and synthesizing information from different sources.

1. Introduction

Mexico, as one of the biologically megadiverse territories in the world, has adopted Natural Protected Areas (NPAs) as its main governmental conservation strategy. Nationwide, 10.87% of the terrestrial and inland water ecosystems are protected under this program (SEMARNAT, 2020). However, far from being “pristine areas,” NPAs in Mexico are generally established in the territories of indigenous groups with long-lived productive and cultural practices (Toledo, 2013), and as such, they constitute highly diverse and complex social-ecological systems (SESs) (Cumming and Allen, 2017). SESs are complex systems in which humans are considered as part of nature (Berkes et al., 2003). In contrast to traditional ecological and social research, studies with a SES perspective explicitly consider ecological and human components and their interactions (Liu et al., 2007). However, most of these studies lack a true integration of both social and ecological variables and their feedbacks (Herrero-Jáuregui et al., 2018). In this sense, to gain insights into the mechanisms underlying SES dynamics, it is necessary to use approaches and tools that account for the intrinsic complexity of these systems (Preiser et al., 2018).

As the world undergoes unprecedented global changes, SES research has been particularly interested in understanding the resilience of SESs. This is generally understood as the capacity of a system to maintain its structure, functionality, identity, and feedbacks in the face of disturbance and ongoing change (Walker et al., 2004). The concept of resilience has promoted new approaches to ecosystem and resource management (Holling and Meffe, 1996) and has become a “boundary concept” fostering communication and collaboration across disciplines (Brand and Jax, 2007; Baggio et al., 2015). However, the operationalization and put in practice of resilience in real-world systems has been difficult because of its conceptual vagueness (Brand and Jax, 2007), the difficulty in quantifying it (Standish et al., 2014) and the limited understanding on the mechanisms that underlie the dynamics of SESs (Schlüter and Pahl-Wostl, 2007; Egli et al., 2018). A way to overcome this latter difficulty is by systematically testing different mechanisms that may affect the response capacity of SESs in the face of disturbances using, for instance, dynamical models (Schlüter and Pahl-Wostl, 2007). As proposed by Egli et al. (2018), this kind of explorations can be guided by the theoretical characterization of resilience-enhancing principles, such as those proposed by Biggs et al. (2012): (1) maintain response diversity and functional redundancy, (2) manage slow variables and feedbacks, (3) manage connectivity, (4) foster an understanding of social-ecological systems as complex adaptive systems, (5) encourage learning and experimentation, (6) broaden participation, and (7) promote polycentric governance systems. So, to advance the practice of resilience, there is a need for a better mechanistic understanding of how these principles may apply in different contexts and how they may interact with each other (Biggs et al., 2012).

The long quest for resilience-enhancing mechanisms on SES has led to the revalorization of some traditional practices, largely overlooked and considered unproductive and environmentally damaging (Berkes et al., 2000). An example of these is the multiple-use strategy (MUS) on which indigenous communities inhabiting tropical forests in Mexico often base their resource management and family household (Toledo et al., 2003). A MUS involves a set of different productive activities (e.g., agriculture, agroforestry, gathering) on a diversity of land units (e.g., milpa, successional forest, mature forest). The core element of the MUS is the milpa, a polyculture which typically combines maize, beans and squash, along with other domesticated, semi-domesticated and tolerated species (Hernandez Xolocotzi et al., 1995; Nigh and Diemont, 2013). The milpa system is highly diversified and adapted to an ample variety of environments, playing a key role in the maintenance of biological diversity and food sovereignty (Altieri et al., 2012). In general terms, a MUS amplifies the subsistence options available for the households to ensure a continuous flow of goods and services, thus minimizing the vulnerability associated with different disturbances (Barrera-Bassols and Toledo, 2005). However, further research is needed to assess the coupled environmental and social outcomes of traditional resource management strategies (Chazdon et al., 2009).

In this paper we study the MUS practiced by a Maya community inhabiting Otoch Ma'ax Yetel Kooh (OMYK, “house of the spider monkey and the jaguar” in Yucatec Maya), a NPA located in the northeastern part of the Yucatan Peninsula, Mexico (Figure 1). This NPA was established in 2002 in response to a local initiative, after years of community-based conservation. It consists of a body of lakes surrounded by tropical forest in multiple successional stages that harbors a large population of spider monkeys and other ecologically important species. OMYK and its buffer zone are inhabited by Yucatec Maya communities, which, until recently, used to manage their resources following a MUS based on traditional swidden milpa agriculture (García-Frapolli et al., 2008). However, due to the growing tourism industry in the region and the restrictions imposed by the decree of the NPA, the local communities have experienced important changes in their management strategies. First, before the decree of the reserve, the community inside the NPA started an ecotourism project based on the observation of spider monkeys in their natural habitat. The MUS, by adding this activity to the traditional ones, promoted the resilience of the SES by diversifying the subsistence means of households (García-Frapolli et al., 2012). However, after the creation of the NPA a management plan was developed which prohibited fire use inside the reserve (CONANP, 2006), thus restricting the production of milpa agriculture in the traditional way, which depends on fire for the preparation of the plots. This situation, in addition to the fast growth of the Riviera Maya tourist corridor which led to an increase on the number of visitors to the reserve, has driven some households to abandon the MUS and specialize in the provision of ecotourism services (García-Frapolli et al., 2012; Bonilla-Moheno et al., 2020). These important changes lead us to ask: how do different productive strategies affect the resilience of this SES in the face of some typical disturbances in the region (e.g., hurricanes, fires, tourism-related fluctuations)? We address this question using a dynamical computational model, a novel approach to understand this SES. This approach allows us to formally and explicitly integrate some of the results of previous studies carried out in the area for a better understanding of the interactions and feedbacks of the social and ecological element of this SES, as well as for exploring different virtual scenarios to systematically investigate different resilience mechanisms (Schlüter and Pahl-Wostl, 2007).

FIGURE 1

Figure 1. Location of the Natural Protected Area Otoch Ma'ax Yetel Kooh (OMYK) showing the main locations mentioned in the text.

Dynamical models are a tool that complements empirical research by facilitating the representation of some key aspects of the complex dynamics of SES (Egli et al., 2018). These models allow researchers to make reproducible experiments in fully controlled virtual SESs, on relatively short periods of time and without directly affecting the people that inhabit them (Barreteau et al., 2001). Agent Based Models (ABMs) are a type of dynamical model that has been widely used in the study of SESs (An, 2012). In ABMs, a system is modeled as a set of individual interacting elements, or agents. Each agent owns a set of variables that describe its state and a set of rules that determine its behavior (Railsback and Grimm, 2012). Some examples of the use of ABM to investigate SESs include the study of ecological degradation scenarios on indigenous communities on Amazonian Guyana (Iwamura et al., 2014, 2016); the study of the impact of demographic human changes on the deforestation of the panda habitat on the Wolong reserve, China (An et al., 2001); and the study of the effect of different management strategies on “La Sepultura” reserve, Mexico (Braasch et al., 2018). Boolean Network Models (BNMs) are another type of dynamical models that provide a qualitative representation of a system. BNM are discrete models in which a system is represented as a directed graph where nodes represent the variables of the system and edges represent regulatory relationships between them (for a detailed description see Saadatpour and Albert, 2013). BNMs have been mainly used in the study of molecular and cellular scale systems and it has not been until recently that they have been used to study ecological and agroecological systems (e.g., Robeva and Murrugarra, 2016; Gaucherel et al., 2017; López-Martínez, 2017).

The objective of this study is to understand the mechanisms by which the MUS may enhance the resilience of the SES associated with OMYK. For this purpose we integrate currently available field evidence into a dynamical computational model that includes both ABM and BNM and use it to explore the effect of different productive and management strategies on the resilience of the SES. We spatially delimit the SES to the area encompassed by the reserve and we represent households, landscape units and important features of the ecosystem (biotic and abiotic). The choice of input variables and rules for their dynamics was based on the ecological and social studies that have been carried out in the area. By running simulations under different scenarios, we analyse the specific resilience (Carpenter et al., 2012) of three output variables (resilience of what?): mature forest area, monkey population size and average monetary value; to three different disturbances (resilience to what?): hurricanes, forest fires and tourism fluctuations. For achieving this we run simulations under different scenarios, where households adopt different resource use strategies or where the SES faces different disturbances. The average monetary value of a household's productive activities, integrates the outcomes of different productive activities and is thus used as an indicator of the resilience of the households' economy. Meanwhile the size of the monkey population and the area of the mature forest, are considered as indicators of ecological resilience. Since spider monkeys are large mammals that require relatively big and continuous extensions of forest with mature native fruiting trees (Ramos-Fernández and Ayala-Orozco, 2003), the viability of their population indirectly reflects the integrity of the local plant communities and the biogeochemical processes that allow for such plant communities to thrive (van Roosmalen and van Roosmalen, 1985; Urquiza-Haas et al., 2009). Similarly, the extent of mature forest provides complementary information regarding the viability of biological populations inhabiting the protected area (Fahrig, 2003). The model allowed us to test the hypothesis that a MUS promotes resilience of the system by exploring the mechanics of how some of the resilience-enhance principles proposed by Biggs et al. (2012) may apply in this particular SES.

2. Methods

2.1. Study Site

OMYK is located in the northeastern region of the Yucatan Peninsula, Mexico, and has an extension of 5,367 ha. Mean annual temperature in the region is 24.3°C and the mean annual precipitation is 1,127.4 mm (SMN, 2019). The region has a tropical wet and dry climate (Aw2) with a dry season from December to April and a wet season from May to November (SMN, 2017).

Dominant vegetation is medium semi-evergreen forest. This forms a heterogeneous landscape composed of a mosaic of vegetation in different successional stages as a result of many years of traditional milpa agriculture and multiple natural disturbances, such as hurricanes and forest fires (Bonilla-Moheno, 2008).

The site has a high faunistic diversity and is inhabited by multiple endangered species, such as spider monkeys (Ateles geoffroyi), jaguars (Panthera onca), pumas (Puma concolor), howler monkeys (Alouatta pigra), and others (CONANP, 2006). The spider monkey population inhabiting the reserve has been widely studied (Ramos-Fernández et al., 2018) and has been the main tourist attraction of the site.

The management plan of the reserve formally recognizes three user communities: Punta Laguna, Campamento Hidalgo, and Nuevo Yodzonot (CONANP, 2006). On this project we focus on Punta Laguna, which given its size and closeness to the reserve has experienced important changes after the decree of the NPA. This community is located in the southeast of the reserve over the roadway Cobá-Nuevo Xcan. Punta Laguna has 136 inhabitants distributed in 28 households. All the inhabitants are indigenous Yucatec Maya (Rivera-Núñez, 2014).

Before the management plan of the reserve entered into force in 2006, households used to manage their resources following a MUS. They used to manage five different land units (milpa, homegardens, secondary forest, old-growth forest, and aquatic systems) and implement 13 different productive activities (milpa agriculture, beekeeping, charcoal production, gather firewood, gather medicinal plants, gather wood for home construction, hunting, fishing, ecotourism, home gardening, sheep herding, scientific research assistance, and temporary work outside of the community; García-Frapolli et al., 2007). Nowadays, there is a tendency toward an abandonment of the MUS and of the traditional activities, and a growing trend toward a specialization on tourism-related activities (García-Frapolli et al., 2012; Bonilla-Moheno et al., 2020).

There has been continuous research on the site since 1996. The research lines are diverse and include: behavioral ecology of the spider monkeys (Ramos-Fernández et al., 2018), forest succession and restoration (Bonilla-Moheno, 2008, 2010; Bonilla-Moheno and Holl, 2010), local systems of nature management (García-Frapolli, 2006; García-Frapolli et al., 2008; Bonilla-Moheno et al., 2020), land cover change analysis (García-Frapolli et al., 2007; Bonilla-Moheno et al., 2020), conservation (Bonilla-Moheno and García-Frapolli, 2012; Rivera-Núñez, 2014; García-Frapolli, 2015) and local institution analysis (García-Frapolli et al., 2013; Rivera-Núñez, 2014). The results of these previous studies were the main sources of information for building the model.

2.2. Description of the Model

The model was built integrating an ABM with a BNM. As the BNM can be understood as part of the ABM (Figure 2; Supplementary Table 1; Supplementary Figure 1) we describe the whole model following the ODD (Overview, Design concepts, and Details) protocol for ABMs (Grimm et al., 2010).

FIGURE 2

Figure 2. Representation of the model. The model can be mainly understood as an agent based model composed of four different types of agents. The world where the agents interact is based on a land use and vegetation map. The model simulates eight different processes, or submodels, six of which occur every 2 months (each time step) and the last two only occur annually (each six time steps). The three submodels indicated with an asterisk (*) were modeled using boolean networks (Supplementary Figure 1). The state of the modeled system is monitored through three output variables.

2.2.1. Purpose

The purpose of the model is to explore how different productive and management strategies affect the resilience of the household's economy, the spider monkey population and the forest of this SES to some frequent disturbances in the region: hurricanes, forest fires and temporality of tourism. The model has been designed for scientists and managers, mainly those interested in conservation and forest management in the region, and it is aimed at exploring scenarios on the relationship between management strategies and the system's resilience, integrating and synthesizing some of the information that has been generated after more than 20 years of multidisciplinary research in the area.

2.2.2. Entities, State Variables, and Scales

The model consists of four agents: landscape patch, household, monkey, and fuel biomass.

Landscape patches represent a square of 3 ha that can be of one of five different types: forest, agriculture (milpa), burned, water or household. Patches of type forest, milpa and burned have a successional age that increases year by year, simulating the regeneration of the forest. When the milpa and burned patches reach a successional age of 2 years they change their type to forest. Patches of type forest and milpa can accumulate fuel biomass, which makes them susceptible to become burned. When a patch is burned or is converted to milpa, its successional age becomes zero and the fuel biomass accumulated on it is consumed. Biannually, households that practice agriculture can convert patches of type forest to milpa by following certain rules (see subsection Household activities). After a milpa has been used for 2 consecutive years, it is abandoned and allowed to regenerate, and the household opens a new plot (García-Frapolli et al., 2007). Water-type patches represent water bodies and the household-type patches represent the areas where households build their homes. Water and household-type patches remain fixed throughout the simulation.

Households have a fixed position in the virtual ABM world. They can work on four different productive activities: milpa agriculture, apiculture, charcoal production and provide ecotourism services. The time dedicated to each productive activity is distributed along the year simulating the annual household's dynamics and temporality of productive activities (Figure 3). The set of activities that households perform remain fixed throughout the simulation. Households obtain a monetary value from each of the productive activities they do. Disturbances, such as hurricanes, forest fires and fluctuations in tourism, can affect the productive activities done by the households, thus reducing their monetary value.

FIGURE 3

Figure 3. Attractor recovered by the submodels built with boolean network modeling. Rows represent each node of the network, columns are the attractor states and each color represents a different state of the node (0, inactive; 1, active). Each of the attractor states of this cyclic attractor of period six can be interpreted as a snapshot of the approximate state of the SES in a given bimester (see subsection Submodels Built With Boolean Network Modeling).

A monkey agent represents a spider monkey. These can live and move in the forest-type patches with a successional age ≥30 years. Monkeys can be born or die following a local logistic growth equation (see section Monkey population dynamics).

Fuel biomass represents an aggregate of decomposing biomass that can ignite a fire. Fuel biomass accumulates on patches after storms and hurricanes. Each fuel biomass unit has a duration time before it disappears. While the fuel biomass is present on a patch, it can cause it to burn.

One time step of the model represents 2 months (see section Climate). The model is run on a 37 × 94 grid where each grid cell represents 3 ha of landscape.

2.2.3. Process Overview and Scheduling

Within a bimester, a sequence of processes takes place in the following order. First, the weather conditions of the bimester are determined (i.e., rainy season or dry season). Second, depending on the weather conditions, different meteorological events may occur (e.g., tropical storms, hurricanes) and depending on the type of event a proportion of patches accumulates fuel biomass. Third, the flow level of tourists visiting the reserve is determined. Fourth, households carry out their productive activities and calculate how much monetary value they obtain. Fifth, forest fires can occur on patches with fuel biomass accumulated. Finally, monkeys can move from patch. Annually (i.e., every six iterations), two additional processes follow: patches regenerate and spider monkey population dynamics take place.

As simulations are performed along discrete time steps, all the state variables of the agents are updated synchronously.

2.2.4. Design Concepts

2.2.4.1. Emergence

The interaction and feedback between the social and ecological systems emerge as a result of land use changes and the productive activities done by the households.

2.2.4.2. Adaptation/Objectives/Learning/Prediction/Sensing

Agents do not present any adaptive behavior, objectives, learning, prediction nor sensing. The set of activities a household carries out makes up a household strategy, but households do not change their strategy along the simulation.

2.2.4.3. Interactions

Households interact indirectly with monkeys through the fragmentation and habitat loss caused by milpa agriculture. Monkeys also indirectly affect the monetary values of the households by influencing the tourist flow. Households also interact directly with the patches by cultivating them and by removing fuel biomass. Spider monkeys are also affected by high levels of tourist flow.

2.2.4.4. Stochasticity

The occurrence of the different meteorological events is stochastically determined based on occurrence probabilities (see sections below). The selection of the patches that accumulate fuel biomass is random. Tourist flow has a stochastic component. The selection of new parcels for cultivation, the site for placing the beehives and the fuel biomass used by the households is partly random. Forest fires are modeled using a stochastic percolation model (see section Forest fires below). And the movement of the monkeys is partly random.

2.2.4.5. Collectives

Contiguous patches of the same type and within the same range of successional age aggregate into patch neighborhoods. Monkeys living in the same patch neighborhood are aggregated into a local population with its own population dynamics.

2.2.4.6. Observation

The average annual monetary value of the productive activities done by the households (Supplementary Table 3) is observed to assess the economic state of the households. The mature forest area (successional age >50 years; the vegetation type that harbors the greater species richness and most rare species; Bonilla-Moheno, 2008), and the total number of spider monkeys in the area are observed to assess the ecological state of the reserve.

2.2.5. Initialization and Input

According to literature reports (Rivera-Núñez, 2014), the number of households is set to 28 and the number of monkeys is set to the maximum possible in each habitable patch neighborhood. To create the virtual world, we used a 2003 land use and vegetation map (García-Frapolli et al., 2007). This map was rasterized and overlaid on the landscape grid. The successional age of the patches is set to the minimum age of the range to which they belong on the base map (2–7, 8–15, 16–29, 30–50, >50 years). The model only considers the area that is occupied by the reserve, so all the grid cells outside of the polygon of the reserve are ignored. The model parameters were set to the values shown in Supplementary Table 2.

2.2.6. Input

The model does not use input data to represent time-varying processes.

2.2.7. Submodels

2.2.7.1. Climate

Weather conditions are modeled using BNM. This submodel is based on the submodel of the same name of López-Martínez (2017) and is used to give a temporality to the model. The submodel is composed of three nodes: temperature, pressure and precipitation. The interpretation of the state of the nodes, regulation functions and their explanations are shown in Supplementary Table 1. The dynamics of this submodel generates a periodic attractor of length 6 where the node precipitation remains in one same state for three consecutive states of the attractor and in a qualitatively different state during the next three states (Figure 3). This dynamics is interpreted as an annual climatic dynamic with 6 months of rainy season and 6 months of dry season. This climatic dynamic is characteristic of the region, so each state of the attractor is mapped to one bimester of the year. This mapping allowed us to locate different events throughout the year.

2.2.7.2. Storms

Tropical storms and hurricanes are a frequent disturbance in the region. The strong winds of hurricanes generate defoliation, snapping and uprooting of trees (Bonilla-Moheno, 2010). So, after these disturbances there is a high accumulation of fire biomass that can generate forest fires in the subsequent dry seasons. To simulate this, a proportion of affected patches (P(s)) is calculated as a function of the speed of the wind (s). We suppose that the relation between these two variables follows a sigmoid function:

$P(s)=100\frac{1}{1+e^{p_1(-s+p_2)}}$

where p₁ and p₂ are parameters that determine the shape and position of the inflection point, respectively. The speed of wind (s) takes a value depending on the type of event that occurs on a bimester: no storm (10 mph), tropical storm (20 mph), and hurricanes from category one to five (74, 96, 111, 130, and 157 mph). In the rainy season, storms and hurricanes can take place with certain occurrence probability estimated with data from NOAA (2019) (Supplementary Table 2). When an event occurs a unit of fuel biomass is created on an amount of patches of type milpa or forest randomly selected equivalent to the proportion of affected patches (P(s)). When created, fuel biomass units are assigned a time of duration from normal distribution (mean duration of fuel biomass, Supplementary Table 2).

2.2.7.3. Tourism

The tourist flow to the reserve is modeled using BNM. This submodel is composed of two nodes, tourism and tourismH, that together represent three levels of tourists flow: no tourists, low season flow and high season flow. The interpretation of the state of the nodes, regulation functions and their explanation are shown in Supplementary Table 1. With these rules we obtain a dynamic with two bimester of high season and four of low season that is similar to the tourist dynamic in the region.

As hurricanes generally reduce the flow of tourists (García-Frapolli et al., 2012), we assume that if a hurricane occurs on a bimester, then no tourists arrive (i.e., the nodes tourism and tourismH turn off). Likewise, as the main tourist attraction at the site are spider monkeys (García-Frapolli et al., 2013), we assumed that if there are not at least 20 monkeys at a ratio of 3 km from the households location then no tourists arrive.

Finally, a stochastic factor was added to simulate the effect of economic and social disturbances that also reduce the number of tourists. This was simulated through a parameter that with certain probability avoids a high flow of tourists (i.e., the node tourismH is kept turned off; probability of high flow of tourists, Supplementary Table 2).

2.2.7.4. Household activities

Household activities are also modeled using BNM. This submodel is composed of seven nodes: openMilpa, plantMilpa, youngMilpa, adultMilpa, harvestMilpa, harvestApiculture, and charcoalProduction. The interpretation of the state of the nodes, regulation functions and their explanation are shown in Supplementary Table 1. With these rules we obtained a dynamic that allows us to locate the realization of the households' productive activities throughout the year (see section Submodels Build With Boolean Network Modeling). In contrast to the other submodels built with BNM, whose nodes are global (i.e., one set of nodes for all the system), the nodes that constitute this submodel are unique for each household considered in the simulation (i.e., each household has its own set of nodes) (Supplementary Figure 1).

To spatially simulate the production of milpa agriculture, when a household opens a milpa it chooses a plot equivalent to 3 ha within a radio of 3 km from its location. Milpas are only opened in patches with successional age ≥5 and <50 years.

To simulate the effect of the hurricanes on the agriculture it is supposed that if a hurricane occurs and the milpa plot of a household is affected by it (i.e., accumulates fossil fuel), then the plants and harvests are lost (i.e., youngMilpa, adultMilpa, and harvestMilpa nodes turn off). Likewise, to simulate the effect of hurricanes on apiculture, it is supposed that a household places their beehives in a patch in a radius of 3 km from its location, and if a hurricane occurs and this patch is affected, then the beehives are damaged and there will be no harvests during the next 2 years. Finally, it is also supposed that households that produce charcoal consume one unit of fuel biomass in a radius of 3 km from its location.

Each bimester, households calculate their total bimonthly monetary value from the values shown in Supplementary Table 3.

2.2.7.5. Forest fires

This submodel uses a modified version of the percolation model of the model library of NetLogo 6.0.4 (Wilensky, 2006). In the model forest fires can only occur with certain occurrence probability during the dry season estimated with data from CONABIO (2019) (Supplementary Table 2). If a forest fire occurs, then one of the patches with most fuel biomass accumulated ignites, and with certain probability (burning probability of patches, Supplementary Table 2) the four neighbor patches containing fuel biomass can also ignite. Then, each burned patch repeats the process.

2.2.7.6. Monkey movement

Spider monkeys can move with a certain probability that depends on the successional age of the patch they are on (probability of permanence, Supplementary Table 2). If a monkey moves, then it randomly shifts to a neighboring patch with successional age ≥30 years. If there are any monkeys on patches of type milpa or burned, then they move to one of their nearest patches with successional age ≥30 years.

2.2.7.7. Forest regeneration

For simplification purposes, we supposed that the regeneration of the forest is deterministic and that the only factor affecting the successional state of the patches is time. So, each year the patches of type milpa, burn and forest increase their successional age.

2.2.7.8. Monkey population dynamics

This submodel is based on the submodel named “Animal meta-population dynamics” of Iwamura et al. (2014). First, patch neighborhoods (set of neighbor patches within the same range of successional age) of late successional forest (30–50 years) and mature forest (>50 years) are formed. Then, for each of these neighborhoods the local population size (N_t) is estimated using the logistic growth equation:

$N_t=N_{t-1}(1+R(1-\frac{N_t-1}{k}))$

Where N_t−1 is the local population size in the time t − 1, R is the discrete intrinsic growth rate and k is the carrying capacity of the neighborhood. R was obtained from the literature, and k for the nth neighborhood, denoted k_n, is calculated in the next way:

$k_n=D_{max}{C}_n$

Where D_max is the maximum density of a patch and C_n is the number of patches that form the nth neighborhood. D_max values were obtained from the literature (Supplementary Table 2).

Finally, as it has been proposed that a high flow of tourists can have an effect on the monkey population, (for instance, the noise generated by tourists and the opening of new trails may distress the monkeys affecting their breeding areas and feeding habits; García-Frapolli et al., 2007), we supposed that the discrete intrinsic growth rate (R) diminishes by a sixth for each bimester that there is a high flow of tourists.

2.3. Implementation

The model was implemented on NetLogo 6.0.4. The code can be found on: https://github.com/laparcela/OMYKmodel.

2.4. Calibration

We made a manual calibration of the following four parameters related to the extent of the forest fires for which no information was founded on the literature:

1. burning probability of patches (values explored: 0.4, 0.5, 0.56, 0.6, 0.7),

2. sigmoid function parameter 1 (p₁) (values explored: 0.05, 0.045, 0.040, 0.035, 0.03, 0.025, 0.02),

3. sigmoid function parameter 2 (p₂) (values explored: 74, 96, 111), and

4. mean duration of fuel biomass (values explored: 6, 9, 12).

We explored a total of 315 different treatments that represent all the possible combinations of the selected values for each parameter (Supplementary Figure 2). We run each treatment 100 times over a time span of 12 years (from 2003 to 2015) in which households do not produce milpa agriculture (as it was prohibited in 2006). We selected the combination of values for the parameters by applying the following criteria:

1. That the total burned area registered by the model was similar to the total area burned observed by Rangel-Rivera (2017) from 2003 to 2015 (i.e., that the empirical data was inside the interquartile range of the model results).

2. That the dispersion of the model results was from the ten smallest.

3. That the burning probability of patches was the least possible. This was done with the objective of avoiding the fires that affect the complete reserve that appear in longer simulations.

4. That the percentage of affected patches when no storm occurs was the lowest possible. This was done to avoid the Forest Fire and Storms submodels to work independently from each other.

After applying this criteria hierarchically (i.e., first applying criteria 1, next criteria 2, and so on) we chose the values for the parameters that are shown in Supplementary Table 2.

2.5. Sensitivity Analysis

We carried out a sensitivity analysis on six parameters as a way of examining the consistency (verification) and the general behavior of the model. The explored variables and the explored values for each were:

1. patch size: 1 and 3 ha;

2. burning probability of patches: from 0 to 1, by increments of 0.1;

3. sigmoid function parameter 1 (p₁): from 0.01 to 0.05, by increments of 0.005;

4. sigmoid function parameter 2 (p₂): 74, 96, 111, and 130;

5. mean duration of fuel biomass: 6, 9, 12, and 15 bimesters;

6. standard deviation of duration of fuel biomass: 1, 3, 6, 9, and 12 bimesters;

7. probability of high flow of tourists: from 0 to 1, by increments of 0.1;

8. minimum number of monkeys for tourism: from 10 to 35, by increments of 5.

We followed the “one at a time” method, in which one parameter is varied at a time while keeping the other fixed on the values shown in Supplementary Table 2. We ran 30 simulations for each treatment during a virtual time span of 200 years (so that the output variables stabilize) in which households produce the four productive activities. We assessed the sensitivity visually looking for noticable changes in the model results.

2.6. Validation Test

Model results were compared with empirical data of the different vegetation and land use classes in 2015 (Rangel-Rivera, 2017) and estimated monkey population size of 2015 (Spaan, 2017). We ran 100 simulations over a virtual time span of 12 years (from 2003 to 2015) in which households do not produce milpa agriculture. Model parameters were fixed on the values shown in Supplementary Table 2. The values of the four parameters for which no data was found in the literature and no calibration was made, were chosen arbitrarily. We look for values that would make simulations faster and that would not generate noticeable changes in the output variables according to the results of the sensitivity analysis.

2.7. Scenarios

To assess the effect of different management strategies on the output variables we explore the 16 different possible combinations of the four productive activities that can be implemented by the households (we only explored the cases where all the households implemented the same set of activities). Each treatment was run 30 times during a time span of 50 years. In this paper we followed Toledo et al. (2003) in their usage of the word strategy as the “intellectual construction or internalized plan among indigenous producers.” To facilitate the analysis of all different simulated scenarios in terms of the household's management strategies, the 16 combinations of productive activities were classified and averaged following the typology of strategies of García-Frapolli et al. (2007, 2008) (Supplementary Table 4):

1. Traditional: combinations where households produce agriculture but no ecotourism.

2. Mixed: combinations where households produce both milpa and ecotourism.

3. Service Oriented: combinations where households dedicated to ecotourism but no agriculture.

4. Other: combinations where households do not produce milpa nor ecotourism (households can only produce charcoal and/or apiculture).

To explore the effect of different management strategies on the resilience of the SES to some frequent disturbances we ran the next set of scenarios:

1. We increased the frequency of hurricanes and tropical storms by arbitrarily multiplying the probabilities of occurrence of tropical storms and of each type of hurricane by 3.

2. We increased the frequency of forest fires by arbitrarily multiplying their probability by 3.

3. We reduced the flow of tourism by arbitrarily reducing the probability of high flow to 0.3.

Each scenario was run 30 times during a time span of 50 years.

3. Results

3.1. Submodels Built With Boolean Network Modeling

The sequence of system states (periodic attractor) recovered by the three submodels built with BNM (just considering a single household) is shown in Figure 3. Each of the states of the attractor can be interpreted as a still image of the approximate conditions of the SES in a given bimester. For instance, the second state of the attractor may be interpreted as the March-April bimester, a bimester of the dry season (node precipitation turned off) in which a high flow of tourists can arrive (nodes tourism and tourismH turned on); in this bimester a household generally prepares the plot for the milpa (node openMilpa turned on); and also in these months a household can obtain harvests of apiculture and produce charcoal (nodes harvestApiculture and charcoalProduction turned on). Similarly, the fifth state of the attractor may be interpreted as a rainy season bimester (September-October), with a low flow of tourists and in which the milpa is almost ready to be harvested and in which households can produce charcoal.

3.2. Sensitivity Analysis

We observed the following behavior patterns for the three output variables:

1. Mature forest area. In most simulations this variable presents three events of remarkable growth as a result of the initial conditions of the model and of the wide length of the age categories used to classify the forest. These events take place at 21, 35, and 50 years and correspond to the moment when the patches that started at 30, 16, and 2 years, respectively, reach the mature forest category (>50 years). After these three events the variable stabilizes.

2. Number of monkeys. In most cases this variable decreases slightly during the first 21 years. Next, it gradually increases until it stabilizes.

3. Average monetary value. In most simulations this variable stabilizes since the beginning and fluctuates around certain equilibrium values. In some cases the equilibrium values gradually decrease until they reach a new, lower equilibrium.

Most output variables were sensible to the modification of the explored parameters. Mature forest area and the total number of monkeys were sensible to the burning probability of patches (Figure 4), the mean duration of the fuel biomass (Figure 5), the two parameters of the sigmoid function (Supplementary Figures 3, 4), and patch size (Supplementary Figure 5). Average monetary value was sensible to the burning probability of patches (Figure 4), patch size (Supplementary Figure 5), the probability of high flow of tourists (Supplementary Figure 6), and the minimum number of monkeys for tourism (Supplementary Figure 7).

FIGURE 4

Figure 4. Sensitivity analysis of the burning probability of patches. Each panel represents the change in an output variable for a 200 years period and 30 simulation runs for each value of the burning probability of patches: (A) Mature forest area in ha (forest >50 years of successional age). (B) Total number of monkeys. (C) Average monetary value of activities done by the households in daily minimum wages. Lines represent the mean and upper and lower shadows represent 5th and 95th quartiles, respectively, for each value of the burning probability of patches (different colors).

FIGURE 5

Figure 5. Sensitivity analysis of the mean duration of fuel biomass (bimesters). Each panel represents the change in an output variable for a 200 years period and 30 simulation runs for each value of the mean duration of fuel biomass. (A) Mature forest area in ha (forest with successional age >50 years). (B) Total number of monkeys. (C) Average monetary value of activities done by the households in daily minimum wages. Lines represent the mean and upper and lower shadows represent 5th and 95th quartiles, respectively, for each value of the mean duration of fuel biomass (different colors).

3.3. Validation Test

The model reproduced relatively well the area of forest of 2–15, 16–29, and 30–50 years (Figure 6). Nevertheless, the model overestimated the area of mature forest and the total number of monkeys.

FIGURE 6

Figure 6. Validation test results. Results of simulations and empirical values for five different variables. (A) Area of forest with successional age between 2 and 15 years. (B) Area of forest with successional age between 16 and 29 years. (C) Area of forest with successional age between 30 and 50 years. (D) Area of forest with successional age >50 years. (E) Total number of monkeys in the reserve. Boxplots are the model outputs after 12 years (from 2003 to 2015) for 100 simulation runs. Diamonds represent the empirical observed values for different forest types in 2015 and the total estimated number of monkeys in 2015.

3.4. Scenarios

The main findings of the simulations are summarized on Table 1. In the scenario without any extra disturbance (Figure 7 first column), the mature forest area was greater for Service oriented and Other strategies. However, there was a growth trend of this variable in all four strategies. Similarly, there was a growth trend of the monkey population in all the strategies, with a larger population increase in Other strategy, followed by the Service oriented and Traditional strategies, and finally by the Mixed strategy. Average monetary value was markedly higher for the strategies including ecotourism, though the variability was higher as the households depended more on this activity (Supplementary Table 5).

TABLE 1

Table 1. Summary of the main findings of the simulated scenarios.

FIGURE 7

Figure 7. Model results for different strategies under different scenarios. Each column of graphs represents a different scenario and each row a different output variable. The four scenarios shown are: a scenario with normal level of all disturbances, a scenario with an increase in the storms and hurricane occurrence probability, a scenario with an increase in the probability of occurrence of forest fires and a scenario with a decrease in the probability of the flow of tourists. First row: Mature forest area in ha (forest with successional age >50 years). Second row: total number of monkeys. Third row: average monetary value of activities done by the households in daily minimum wages. Lines represent the mean and upper and lower shadows represent 5th and 95th quartiles, respectively, for different strategies (in colors), each ran over a 50 years period and for 30 simulation realizations.

When the storm occurrence probability was increased there was just a slight reduction of the trajectories of the three output variables (Figure 7 second column). In the scenario of increased forest fires (Figure 7 third column) there was a considerable reduction of mature forest area and monkey population for all strategies, but there were just slight changes in the average monetary value. In contrast, when tourist flow was reduced (Figure 7 fourth column), there were no noticeable changes in the trajectories for the mature forest area and the number of monkeys, yet there was a remarkable reduction of the average monetary value for the strategies that included this activity when compared with the scenario with no extra disturbances as a reference state. The reduction of the aggregated monetary value under all disturbances was higher for the Service oriented strategy than for any other strategy (Figure 8, Supplementary Table 5).

FIGURE 8

Figure 8. Aggregate average monetary value for different strategies under different disturbance scenarios. The four scenarios shown are: a scenario with normal level of all disturbances, a scenario with an increase in the storms and hurricane occurrence probability, a scenario with an increase in the probability of occurrence of forest fires and a scenario with a decrease in the probability of the flow of tourists. Boxplots represent the yearly average monetary value of different strategies (in colors) for 30 simulations, each ran over a 50 years period.

4. Discussion

The model we have used here enabled the formal integration of ecological and social data and information collected over the past 20 years, and allowed us to test how different management strategies affect the socio-ecosystemic resilience in the long-term. To our knowledge, this is the first computational model for the study of SESs that explicitly integrates BNMs and ABMs. One of the main constraints of the latter is that they require large amounts of empirical data to be parameterized (Iwamura et al., 2014; Schulze et al., 2017). In contrast, BNMs require few or no parameters, making them a suitable tool for modeling scarcely characterized systems (Saadatpour and Albert, 2013). In this sense, the qualitative approach offered by the BNMs can facilitate and simplify the construction of models for the study of SES. In what follows, we first show some of the mechanisms by which the MUS practiced by the local communities may enhance the resilience of this SES. Specifically, we show how some of the resilience-enhancing principles proposed by Biggs et al. (2012) may apply in this particular SES. Next, we discuss the validation and sensitivity analysis results and highlight some limitations of the model.

When comparing the average monetary value obtained by Traditional and Other strategies with those obtained by the Service oriented and Mixed strategy, the significant economic benefits of the incorporation of ecotourism becomes evident. This is consistent with studies in the area that have shown that regular and relatively high incomes from this activity significantly reduce the short-term economic vulnerability of the households, as they reduce the need of the household members to search for temporary work outside their communities, improving family well-being (García-Frapolli et al., 2007; Rivera-Núñez, 2014). However, our results suggest that the recent tendency of specialization on the provision of ecotourism services and abandonment of the milpa agriculture, the central element of the MUS, reduces the resilience of this SES in several ways.

First, in our model the average monetary value of the Mixed strategy showed a greater resistance (i.e., the change of a state variable after a disturbance; Egli et al., 2018) and a lesser variability than the other three strategies in the face of all disturbances. This can be explained by the greater functional redundancy and response diversity that the Mixed strategy promotes (Biggs et al., 2012). In this context, functional redundancy is represented by the number of productive activities the household implements. However, only implementing multiple productive activities without assuring their disparity, i.e., that the elements are different from one another (Biggs et al., 2012), does not ensure resilience. For instance, the integration of charcoal production and provision of ecotourism services results in a low response diversity as these two activities strongly depend on the regional tourism industry. In contrast, a Mixed strategy enables two alternative pathways of resource fluxes: market-oriented and subsistence-oriented activities. This integration results in a higher response diversity as the failures and losses on market oriented-activities due, for instance, to the temporality of tourism or economic and health crises, are buffered by the goods produced by subsistence activities. Similarly, disturbances that mostly impact subsistence activities, such as droughts and hurricanes, may be buffered by the income of market oriented activities. However, as the model clearly shows, there is a noticeable economic trade-off between the productive specialization, represented by the Service oriented strategy, and the diversification, represented by the Mixed strategy. Specialization allows higher earnings (as households can spend more time in ecotourism related activities) with a less labor intensive activity, although compromising the household's economy in face of disturbances. On the contrary, diversification is a more labor intensive strategy that generates lower economic gains (as market forces unequally favor the ecotourism-related over traditional ones) but reduces the income vulnerability caused by the frequent disturbances in the region. In general, these results agree with studies that suggest that economic diversification helps rural households reduce poverty and vulnerability (Ellis, 2008; Thulstrup, 2015; Martin and Lorenzen, 2016) and studies that document the greater vulnerability of tourism specialization due to its strong dependency on a single highly temporal and variable activity (Tao and Wall, 2009; Su et al., 2016).

Yet another way in which a greater functional redundancy and response diversity enhance resilience is in the greater adaptation capacity they promote. The MUS is characterized by its flexibility, based on a second resilience-enhancing mechanism: active learning and experimentation (Biggs et al., 2012). The MUS practiced by the local households has allowed them to adapt to ongoing ecological, social and economic changes, as they have constantly integrated and abandoned different productive activities (García-Frapolli et al., 2012). Nevertheless, none of the previous adaptations had implied a decrease in the diversity of household activities as has the integration of ecotourism (García-Frapolli et al., 2008; Bonilla-Moheno et al., 2020).

Furthermore, specialization also threatens food sovereignty and biocultural diversity. As specialization on ecotourism increases, there is a greater dependency of communities on external agents and a loss of control over their subsistence. This means, for example, that communities are forced to buy products that they previously produced and decided how to produce. The abandonment of milpa agriculture is also worrying, given that it is a central component of the Yucatec Maya identity and, as such, it encompasses a great diversity of knowledge, practices and beliefs, which may be lost (Barrera-Bassols and Toledo, 2005; Lyver et al., 2019). Similarly, tourism dynamics have led to a process of acculturation, as shown by the way that the traditional ceremonies and rituals have been commodified into commercial products for the tourists (García-Frapolli et al., 2018). All these events, together, imply a challenge to a third resilience mechanism of the MUS: fostering an understanding of SESs as complex adaptive systems (Biggs et al., 2012). Which in this context is represented by the holistic cosmovision of the Yucatec Maya, that includes some features that parallel complex thinking, such as the recognition of surprises, of multiple scales and of strong interconnections between elements (Barrera-Bassols and Toledo, 2005).

For the ecological variables considered, our results show that the recovery tendency of the tropical forest and the growth tendency of the spider monkey population held in the face of all the different management strategies. This suggests that the forest and spider monkeys may be resilient to the disturbances generated by the kind of productive activities that have taken place in the protected area. This is particularly interesting in the case of the strategies that include milpa production. The small effect of milpa agriculture on the prevalence and increment of mature forest is due to the low rate of milpa-driven land use change and to the community agreement on not to use the mature forest surrounding the lake for agricultural activities (García-Frapolli et al., 2007). This last initiative, based on traditional local ecological knowledge, allowed the conservation of one of the main patches of mature forest used by the spider monkeys. As the model results suggest, this zonification, that did not eliminate completely milpa agriculture from the reserve, may have been a sufficient measure to conserve the spider monkey population in the long term and ensure the arrival of tourists. In addition, spider monkeys may be more flexible to the degradation of their habitat than previously thought (Spaan, 2017). In the reserve, spider monkeys are not restricted to mature forest as they also use and travel on late successional forest fragments (30–50 years old; Ramos-Fernández and Ayala-Orozco, 2003; Ramos-Fernandez et al., 2013). This allows them to exploit a wide range of environments in a fragmented habitat. More generally, in spite of the general perception that agricultural activities are a threat to primate biodiversity, there is evidence of primate presence on a diversity of agroecosystems (Estrada et al., 2006).

When increasing the probability of hurricane occurrence, both ecological variables maintain their growing trends. We didn't explicitly take into consideration any mechanism that may explain this behavior, but these results agree with studies that suggest that the forest in the region has evolved to become resilient to this frequent disturbance (Sánchez Sánchez and Islebe, 1999; Chazdon, 2003; Bonilla-Moheno, 2010). In contrast, the disturbance that most affected the mature forest and the monkey population was the increase in forest fires. The importance of this disturbance on the SES was also observed in the sensitivity analysis when modifying the burning probability of patches and the mean duration of biomass fuel. When varying the burning probability of patches we obtained the behavior characteristic of the percolation models, with a threshold that when exceeded prevents the growth of the mature forest and the monkey population. It is worth noting that in order to reproduce the empirically observed total burned area, in the calibration, it was found that this parameter has to be set near the critical point. Indeed, it has been proposed that some ecological systems organize near the critical points, particularly regarding the occurrence of forest fires (Ricotta et al., 1999). Similarly, when varying the mean duration of fuel biomass, we found a clear reduction of the mature forest area and of the spider monkey population. This was due to the greater extent of the forest fires given the greater proportion of affected patches at a given point of time, which exemplifies a fourth resilience mechanism: the management of connectivity (Biggs et al., 2012). This suggests that some traditional activities that involve the reduction of the volume of fuel material, such as the gathering of firewood and wood for construction and the controlled burning for traditional agriculture, may be important mechanisms to prevent large-scale forest fires as they may reduce the connectivity between sites with high probability of fires. However, as some of these activities are mainly seen as detrimental to the environment and some, as the controlled burning for agriculture, evidently involve a high fire risk, they have been prohibited in the reserve (CONANP, 2006). Nevertheless, regardless of these prohibitions, the lack of proper regulation over the charcoal production on neighborhood communities to the reserve, caused the large forest fires on the reserve between 2006 and 2011 after the 2005 hurricanes Emily and Wilma (Rangel-Rivera, 2017; Bonilla-Moheno et al., 2020).

Throughout most of the twentieth century, environmental policy on fire in most countries, including Mexico, was based on a fire suppression approach (Martínez-Torres et al., 2016). This approach is responsible for the widely negative perception of fire use in traditional agriculture, as it blamed the practices of rural and indigenous people as the main cause of forest fires. Currently, Mexico's environmental policy on fire is on a transition toward an integrated management approach which recognizes the role that fires could play in ecosystems (Gutiérrez Navarro et al., 2017). However, as environmental policy on fire responds mainly to ecological concerns, it continues to ignore the practices and knowledge on fire use of rural and indigenous people (Monzón-Alvarado et al., 2014; Martínez-Torres et al., 2016; Gutiérrez Navarro et al., 2017). In the case of the Yucatec Maya, “wind tenders” (yum ik'ob in Yucatec Maya) are the people who carry out the burns in the milpa production. Wind tenders have developed a variety of techniques to control fire, preventing it from expanding to the surrounding fields by placing firebreaks, while also assuring the complete conversion of biomass into ashes (Nigh and Diemont, 2013). Given the need of effective large-scale fire prevention strategies and the profound traditional ecological knowledge of the Maya people, it seems essential to revalue local practices and knowledge regarding fire management. We do not suggest that traditional or local practices are ideal, but that they should not be ignored and that the local communities' perspectives, needs and knowledge should not be excluded from fire management plans (García-Frapolli, 2015). In general, our findings suggest that forest fire prevention strategies in NPAs need to go beyond the recognition and prohibition of productive activities with possible detrimental effects on the environment, and also recognize the role that these activities may play in ecosystem management and resilience-enhancing.

Summing up, our results add to those of previous studies (García-Frapolli et al., 2007; Bonilla-Moheno et al., 2020) by suggesting that, contrary to the conventional belief, there is no trade-off between traditional milpa agriculture and biodiversity conservation. On the contrary, there may be a synergy between these two activities as traditional agriculture is strongly related to multiple practices and traditional knowledge that when implemented through a MUS enhance the resilience of this SES. However, the management plan of OMYK, instead of protecting and promoting the milpa and the MUS, has become the main driver of their abandonment, promoting the simplification of this complex SES. Instead of adopting reductionist approaches, management plans of NPAs should take into consideration the (sometimes hidden) trade-offs and synergies that characterize the complex dynamics of these SESs. As tourism keeps growing as a result of its promotion in the region, and as Mexican conservation policies keep facing important challenges on the inclusion of local population's needs and perspectives (García-Frapolli et al., 2009; Durand et al., 2014), the processes occurring in this particular SES may find parallels with other similar SESs in the region. Thus, without ignoring the risks of overgeneralizing, this SES may constitute a model system, providing lessons to inform and give insights and hypotheses into what may happen on similar, but less characterized SESs. For instance, the so-called Maya Train, one of the main infrastructure projects of the federal government planning to economically connect multiple southeastern states of Mexico (SEGOB, 2019), is expected to drive important changes on local communities as they integrate to the tourism and industrial sectors. The particular case of OMYK helps to highlight some of the potential threats of this project, such as the potential loss of some resilience-enhancing mechanisms, an increased dependency on outside agents and the loss of food sovereignty. Consequently, in the face of these projects, there is a greater need to protect and promote the MUS over the specialized strategies.

The results of the sensitivity analysis suggest that the model has an important degree of uncertainty associated with the modification of its parameters. Nevertheless, the behaviors observed in the output variables can be explained by the assumptions of the conceptual model (i.e., the mechanisms proposed in subsection 2.2.7), which increases the confidence in the model implementation. In the case of the validation test, the model was able to reproduce quite well the empirically observed area of successional forest categories, but failed to reproduce the observed mature forest area and estimated number of monkeys. This may be partly due to the unusually large forest fires that took place between 2006 and 2011, that damaged large areas at the north of the reserve, affecting particularly one of the main patches of mature forest (Rangel-Rivera, 2017; Bonilla-Moheno et al., 2020). Even when having used field data to calibrate the model, these events were not considered explicitly in the model simulations. Additionally, it is important to notice that the calibration and the validation test were done with data of the same study, which can explain the good fit of some of the model results. So, to better test the predictive capacity of the model it will be necessary to compare it with future independent studies. Yet, even if the model may not be correctly validated for all variables and even if it may not generate precise quantitative predictions, it reproduces fairly well the mid-term qualitative system behavior described in the literature. As such, the model proved to be a useful tool to explore the joint effect of the considered variables and gain new insights of this SES.

The model presented in this work is a simplified representation of a complex SES and has multiple limitations. Some of these are necessary assumptions of the modeling exercise, while others arise from the lack of data and the difficulty of representing political, cultural and historical processes and contingencies in mathematical models. Here we briefly discuss four limitations of our model. First, the model only considers the area occupied by the reserve, ignoring all the ecological and social processes that take place outside its limits. For instance, although the habitat of the monkeys may be recovering within and around the reserve (Bonilla-Moheno et al., 2020), other processes (e.g., intensification of agricultural activities or expansion of urban areas) may be occurring at a larger scale that may jeopardize the long-term permanence of monkey populations (e.g., loss of connectivity between habitats in the region). Second, the model does not consider the human demographic changes. Thus, we could not explore how population growth may promote the expansion of agriculture and urban areas and the overall consequences of such expansion. However, population growth in OMYK does not necessarily imply an expansion of the agricultural land given the existence of alternative productive activities (García-Frapolli et al., 2007). Third, the income distribution is rarely equitable and accessible to all, as assumed in the model. In OMYK, while the community of Punta Laguna has highly benefited from tourism, other communities have been excluded or do not have consolidated initiatives that allow them to enjoy the same benefits (Rivera-Núñez, 2014). Even among the members of Punta Laguna, there has been an exclusion of the households that have remained more attached to their traditional activities (García-Frapolli et al., 2018). Fourth, the model does not consider some of the main current drivers of change in this SES. Since the reserve's decree, members of the community of Punta Laguna have faced strong conflicts over land tenure with some members of the ejido assembly (a collective form of land tenure). If the conflict continues, it could enable the entry of external agents, quite likely from the tourism development and real estate industries.

In conclusion, in this study we show how different resilience mechanisms may apply in the SES associated with OMYK. In particular, we have shown how the MUS promotes a greater response diversity and functional redundancy and has key positive influences on resilience, such as the management of connectivity, the promotion of learning and experimentation, and the fostering of an understanding of SESs as complex adaptive systems. In addition, the model presented here shows how computational modeling and the SES perspective are effective means of integrating and synthesizing information from different sources. Our model can be used as a guide or support tool for future integrative studies in OMYK and other similar SESs. Likewise, the model can be used as a guide toward the development of support tools in participatory processes. This type of tools can help in the communication of knowledge, promote collective reflection and learning, and thus facilitate collective and adaptive decision-making for building sustainable management practices and effective governance (D'Aquino et al., 2003; Braasch et al., 2018). Finally, our work also illustrates some of the main limitations of a systems perspective and computational modeling on the study of SES, such as the difficulty in considering political, cultural and historical features.

Data Availability Statement

The model and figures code is available at https://github.com/laparcela/OMYKmodel. Simulations results datasets are available on request to the corresponding authors.

Author Contributions

All authors participated in an initial workshop where this study was first conceptualized. LG-J, MB, and GR-F conducted the research. EG-F, MB-M, and CR-R provided data and maps for the model. LG-J wrote the original draft. All authors reviewed and edited the paper.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer YM-U declared a shared affiliation, though no other collaboration, with several of the authors LG-J, MB, and GR-F to the handling editor.

Acknowledgments

We thank the community of Punta Laguna for their support of academic research over the years. We also thank Celene Espadas, Fernanda Figueroa, Fernanda Ríos-Beltrán, Tlacaelel Rivera-Núñez, and Ernesto Vega, members of La Parcela laboratory, Lorena Soto-Pinto and three reviewers for their commentaries on previous versions of this work. Ana María García helped proofread the manuscript. The Centro de Ciencias de la Complejidad (C3) at the National Autonomous University of Mexico hosted a workshop where some of the ideas that gave rise to this work were discussed. MB acknowledges financial support from UNAM-DGAPA-PAPIIT (IN207819). GR-F acknowledges support from CONACYT (grant 157353), UNAM-DGAPA-PAPIIT (IA200720), Instituto Politécnico Nacional, and EG-F acknowledges financial support form UNAM-PAPIIT (IN300520).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fsufs.2020.522657/full#supplementary-material

References

Altieri, M. A., Funes-Monzote, F. R., and Petersen, P. (2012). Agroecologically efficient agricultural systems for smallholder farmers: contributions to food sovereignty. Agron. Sustain. Dev. 32, 1–13. doi: 10.1007/s13593-011-0065-6

CrossRef Full Text | Google Scholar

An, L. (2012). Modeling human decisions in coupled human and natural systems: review of agent-based models. Ecol. Model. 229, 25–36. doi: 10.1016/j.ecolmodel.2011.07.010

CrossRef Full Text | Google Scholar

An, L., Liu, J., Ouyang, Z., Linderman, M., Zhou, S., and Zhang, H. (2001). Simulating demographic and socioeconomic processes on household level and implications for giant panda habitats. Ecol. Model. 140, 31–49. doi: 10.1016/S0304-3800(01)00267-8

CrossRef Full Text | Google Scholar

Baggio, J. A., Brown, K., and Hellebrandt, D. (2015). Boundary object or bridging concept? A citation network analysis of resilience. Ecol. Soc. 20:2. doi: 10.5751/ES-07484-200202

CrossRef Full Text | Google Scholar

Barrera-Bassols, N., and Toledo, V. M. (2005). Ethnoecology of the Yucatec Maya: symbolism, knowledge and management of natural resources. J. Latin Am. Geogr. 4, 9–41. doi: 10.1353/lag.2005.0021

CrossRef Full Text | Google Scholar

Barreteau, O., Bousquet, F., and Attonaty, J.-M. (2001). Role-playing games for opening the black box of multi-agent systems: method and lessons of its application to Senegal Rivere Valley irrigated systems. J. Artif. Soc. Soc. Simul. 4, 1–5. Available online at: http://jasss.soc.surrey.ac.uk/4/2/5.html

Google Scholar

Berkes, F., Colding, J., and Folke, C. (2000). Rediscovery of traditional ecological knowledge as adaptive management. Ecol. Appl. 10, 1251–1262. doi: 10.1890/1051-0761(2000)010[1251:ROTEKA]2.0.CO;2

CrossRef Full Text | Google Scholar

Berkes, F., Colding, J., and Folke, C. (2003). Navigating Social-Ecological Systems: Building Resilience for Complexity and Change. New York, NY: Cambridge University Press.

Google Scholar

Biggs, R., Schlüter, M., Biggs, D., Bohensky, E. L., BurnSilver, S., Cundill, G., et al. (2012). Toward principles for enhancing the resilience of ecosystem services. Annu. Rev. Environ. Resour. 37, 421–448. doi: 10.1146/annurev-environ-051211-123836

CrossRef Full Text | Google Scholar

Bonilla-Moheno, M. (2008). Forest recovery and management options in the Yucatan Peninsula, Mexico (Doctorado), University of California, Santa Cruz, Santa Cruz, CA, United States.

Bonilla-Moheno, M. (2010). Damage and recovery of forest structure and composition after two subsequent hurricanes in the Yucatan Peninsula. Carib. J. Sci. 46, 240–248. doi: 10.18475/cjos.v46i2.a12

CrossRef Full Text | Google Scholar

Bonilla-Moheno, M., and García-Frapolli, E. (2012). Conservation in context: a comparison of conservation perspectives in a mexican protected area. Sustainability 4, 2317–2333. doi: 10.3390/su4092317

CrossRef Full Text | Google Scholar

Bonilla-Moheno, M., and Holl, K. D. (2010). Direct seeding to restore tropical mature-forest species in areas of slash-and-burn agriculture. Restor. Ecol. 18, 438–445. doi: 10.1111/j.1526-100X.2009.00580.x

CrossRef Full Text | Google Scholar

Bonilla-Moheno, M., Rangel-Rivera, C. E., García-Frapolli, E., Ríos Beltrán, F., Espadas-Manrique, C., Aureli, F., et al. (2020). Protected Areas Are Effective Conservation Strategies But Promote the Simplification of the Social-ecological System: A Case Study on Land-Use, Vegetation Cover, and Resource Management Strategies.

Braasch, M., García-Barrios, L., Cortina-Villar, S., Huber-Sannwald, E., and Ramírez-Marcial, N. (2018). TRUE GRASP: actors visualize and explore hidden limitations of an apparent win-win land management strategy in a MAB reserve. Environ. Model. Softw. 105, 153–170. doi: 10.1016/j.envsoft.2018.03.022

CrossRef Full Text | Google Scholar

Brand, F. S., and Jax, K. (2007). Focusing the meaning(s) of resilience: resilience as a descriptive concept and a boundary object. Ecol. Soc. 12:23. doi: 10.5751/ES-02029-120123

CrossRef Full Text | Google Scholar

Carpenter, S., Arrow, K., Barrett, S., Biggs, R., Brock, W., Crépin, A.-S., et al. (2012). General resilience to cope with extreme events. Sustainability 4, 3248–3259. doi: 10.3390/su4123248

CrossRef Full Text | Google Scholar

Chazdon, R. L. (2003). Tropical forest recovery: legacies of human impact and natural disturbances. Perspect. Plant Ecol. Evol. Syst. 6, 51–71. doi: 10.1078/1433-8319-00042

CrossRef Full Text | Google Scholar

Chazdon, R. L., Harvey, C. A., Komar, O., Griffith, D. M., Ferguson, B. G., Martínez-Ramos, M., et al. (2009). Beyond reserves: a research agenda for conserving biodiversity in human-modified tropical landscapes. Biotropica 41, 142–153. doi: 10.1111/j.1744-7429.2008.00471.x

CrossRef Full Text | Google Scholar

CONABIO (2019). Puntos de Calor Detectados con Imágenes de Satélite. Mexico: CONABIO. Available online at: http://terra2.conabio.gob.mx/prueba_dom1/consultasAvanzadas.php

CONANP (2006). Programa de Conservación y Manejo Área de Protección de Flora y Fauna Otoch Ma'ax Yetel Kooh: México. Ciudad de México: Comisión Nacional de Áreas Naturales Protegidas (CONANP).

Cumming, G. S., and Allen, C. R. (2017). Protected areas as social-ecological systems: perspectives from resilience and social-ecological systems theory. Ecol. Appl. 27, 1709–1717. doi: 10.1002/eap.1584

PubMed Abstract | CrossRef Full Text | Google Scholar

D'Aquino, P., Le Page, C., Bousquet, F., and Bah, A. (2003). Using self-designed role-playing games and a multi-agent system to empower a local decision-making process for land use management: the selfcormas experiment in Senegal. J. Artif. Soc. Soc. Simul. 6. Available online at: http://jasss.soc.surrey.ac.uk/6/3/5.html

Google Scholar

Durand, L., Figueroa, F., and Trench, T. (2014). Inclusion and exclusion in participation strategies in the Montes Azules Biosphere Reserve, Chiapas, Mexico. Conserv. Soc. 12, 175–189. doi: 10.4103/0972-4923.138420

CrossRef Full Text | Google Scholar

Egli, L., Weise, H., Radchuk, V., Seppelt, R., and Grimm, V. (2018). Exploring resilience with agent-based models: state of the art, knowledge gaps and recommendations for coping with multidimensionality. Ecol. Complex. 40:100718. doi: 10.1016/j.ecocom.2018.06.008

CrossRef Full Text | Google Scholar

Ellis, F. (2008). The determinants of rural livelihood diversification in developing countries. J. Agric. Econ. 51, 289–302. doi: 10.1111/j.1477-9552.2000.tb01229.x

CrossRef Full Text | Google Scholar

Estrada, A., Saenz, J., Harvey, C., Naranjo, E., Muñoz, D., and Rosales-Meda, M. (2006). “Primates in agroecosystems: conservation value of some agricultural practices in Mesoamerican landscapes,” in New Perspectives in the Study of Mesoamerican Primates: Distribution, Ecology, Behavior, and Conservation, Developments in Primatology: Progress and Prospects, eds A. Estrada, P. A. Garber, M. S. M. Pavelka, and L. Luecke (Boston, MA: Springer US), 437–470. doi: 10.1007/0-387-25872-8_22

CrossRef Full Text | Google Scholar

Fahrig, L. (2003). Effects of habitat fragmentation on biodiversity. Annu. Rev. Ecol. Evol. Syst. 34, 487–515. doi: 10.1146/annurev.ecolsys.34.011802.132419

PubMed Abstract | CrossRef Full Text | Google Scholar

García-Frapolli, E. (2006). Conservation from below: socioecological systems in natural protected areas in the Yucatan Peninsula, Mexico (Doctorado), Universidad Autónoma de Barcelona, Barcelona, Spain.

Google Scholar

García-Frapolli, E. (2015). “Exclusión en Áreas Naturales Protegidas: Una Aproximación desde los Programas de Manejo,” in La naturaleza en contexto: Hacia una ecología política mexicana, eds L. Durand, F. Figueroa, and M. Guzmán (México: Centro de Investigaciones Interdisciplinarias en Ciencias y Humanidades), 221–236.

Google Scholar

García-Frapolli, E., Ayala-Orozco, B., Bonilla-Moheno, M., Espadas-Manrique, C., and Ramos-Fernández, G. (2007). Biodiversity conservation, traditional agriculture and ecotourism: Land cover/land use change projections for a natural protected area in the northeastern Yucatan Peninsula, Mexico. Landsc. Urban Plan. 83, 137–153. doi: 10.1016/j.landurbplan.2007.03.007

CrossRef Full Text | Google Scholar

García-Frapolli, E., Ayala-Orozco, B., Bonilla-Moheno, M., Espadas-Manrique, C., and Ramos-Fernández, G. (2012). “Conservación Comunitaria en Punta Laguna: Ecoturismo y Desarrollo Sustentable,” in Casos exitosos en la construcción de sociedades sustentables, Quehacer científico y tecnológico. Serie Hacia la sustentabilidad, Primera Edición Edition, eds E. Silva Rivera, M. d. C. Vergara Tenorio, and E. Rodríguez-Luna (Xalapa: Universidad Veracruzana), 413–430.

García-Frapolli, E., Ayala-Orozco, B., Oliva, M., and Smith, R. J. (2018). Different approaches towards the understanding of socio-environmental conflicts in protected areas. Sustainability 10:2240. doi: 10.3390/su10072240

CrossRef Full Text | Google Scholar

García-Frapolli, E., Bonilla-Moheno, M., and Ramos-Fernández, G. (2013). “Community conservation in Punta Laguna: a case of adaptive ecotourism management,” in Community Action for Conservation, eds L. Porter-Bolland, I. Ruiz-Mallén, C. Camacho-Benavides, and S. R. McCandless (New York, NY: Springer), 101–113. doi: 10.1007/978-1-4614-7956-7_7

CrossRef Full Text | Google Scholar

García-Frapolli, E., Ramos-Fernández, G., Galicia, E., and Serrano, A. (2009). The complex reality of biodiversity conservation through natural protected area policy: three cases from the Yucatan Peninsula, Mexico. Land Use Policy 26, 715–722. doi: 10.1016/j.landusepol.2008.09.008

CrossRef Full Text | Google Scholar

García-Frapolli, E., Toledo, V. M., and Martinez-Alier, J. (2008). Adaptations of a Yucatec Maya multiple-use ecological management strategy to ecotourism. Ecol. Soc. 13:31. doi: 10.5751/ES-02627-130231

CrossRef Full Text | Google Scholar

Gaucherel, C., Théro, H., Puiseux, A., and Bonhomme, V. (2017). Understand ecosystem regime shifts by modelling ecosystem development using Boolean networks. Ecol. Complex. 31, 104–114. doi: 10.1016/j.ecocom.2017.06.001

CrossRef Full Text | Google Scholar

Grimm, V., Berger, U., DeAngelis, D. L., Polhill, J. G., Giske, J., and Railsback, S. F. (2010). The ODD protocol: a review and first update. Ecol. Model. 221, 2760–2768. doi: 10.1016/j.ecolmodel.2010.08.019

CrossRef Full Text | Google Scholar

Gutiérrez Navarro, A., García Barrios, L. E., Parra Vázquez, M., and Rosset, P. (2017). De la supresión al manejo del fuego en la Reserva de la Biosfera La Sepultura, Chiapas: perspectivas campesinas. Región Soc. 29, 31–70. doi: 10.22198/rys.2017.70.a329

CrossRef Full Text | Google Scholar

Hernandez Xolocotzi, E., Baltazar, E., and Tacher, S. (1995). La milpa en Ycatán: un sistema de producción agrícola tradicional. Montecillo: Colegio De Postgraduados.

Google Scholar

Herrero-Jáuregui, C., Arnaiz-Schmitz, C., Reyes, M. F., Telesnicki, M., Agramonte, I., Easdale, M. H., et al. (2018). What do we talk about when we talk about social-ecological systems? A literature review. Sustainability 10:2950. doi: 10.3390/su10082950

CrossRef Full Text | Google Scholar

Holling, C. S., and Meffe, G. K. (1996). Holling and Meffe.pdf. Conserv. Biol. 10, 328–337. doi: 10.1046/j.1523-1739.1996.10020328.x

CrossRef Full Text | Google Scholar

Iwamura, T., Lambin, E., Silvius, K., Luzar, J., and Fragoso, J. (2014). Agent-based modeling of hunting and subsistence agriculture on indigenous lands: understanding interactions between social and ecological systems. Environ. Model. Softw. 58, 109–127. doi: 10.1016/j.envsoft.2014.03.008

CrossRef Full Text | Google Scholar

Iwamura, T., Lambin, E., Silvius, K., Luzar, J., and Fragoso, J. (2016). Socio-environmental sustainability of indigenous lands: simulating coupled human-natural systems in the Amazon. Front. Ecol. Environ. 14, 77–83. doi: 10.1002/fee.1203

CrossRef Full Text | Google Scholar

Liu, J., Dietz, T., Carpenter, S. R., Alberti, M., Folke, C., Moran, E., et al. (2007). Complexity of coupled human and natural systems. Science 317, 1513–1516. doi: 10.1126/science.1144004

PubMed Abstract | CrossRef Full Text | Google Scholar

López-Martínez, R. (2017). Modelos Dinámicos de Redes Ecológicas para un Sistema Productivo Complejo: la Milpa Mexicana (Maestría), Universidad Nacional Autónoma de México, Ciudad de México, Mexico.

Lyver, P. O., Timoti, P., Davis, T., and Tylianakis, J. M. (2019). Biocultural hysteresis inhibits adaptation to environmental change. Trends Ecol. Evol. 34, 771–780. doi: 10.1016/j.tree.2019.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, S. M., and Lorenzen, K. (2016). Livelihood diversification in rural Laos. World Dev. 83, 231–243. doi: 10.1016/j.worlddev.2016.01.018

CrossRef Full Text | Google Scholar

Martínez-Torres, H. L., Castillo, A., Ramírez, M. I., and Pérez-Salicrup, D. R. (2016). The importance of the traditional fire knowledge system in a subtropical montane socio-ecosystem in a protected natural area. Int. J. Wildl. Fire 25, 911–921. doi: 10.1071/WF15181

CrossRef Full Text | Google Scholar

Monzón-Alvarado, C., Waylen, P., and Keys, E. (2014). Fire management and climate variability: challenges in designing environmental regulations. Land Use Policy 39, 12–21. doi: 10.1016/j.landusepol.2014.03.003

CrossRef Full Text | Google Scholar

Nigh, R., and Diemont, S. A. (2013). The Maya milpa: fire and the legacy of living soil. Front. Ecol. Environ. 11, e45–e54. doi: 10.1890/120344

CrossRef Full Text | Google Scholar

NOAA (2019). IBTrACS–International Best Track Archive for Climate Stewardship. Silver Spring, MD: NOAA. Available online at: https://www.ncdc.noaa.gov/ibtracs/index.php?name=ibtracs-data

Google Scholar

Preiser, R., Biggs, R., De Vos, A., and Folke, C. (2018). Social-ecological systems as complex adaptive systems: organizing principles for advancing research methods and approaches. Ecol. Soc. 23:46. doi: 10.5751/ES-10558-230446

CrossRef Full Text | Google Scholar

Railsback, S. F., and Grimm, V. (2012). Agent-Based and Individual-Based Modeling: A Practical Introduction. Princeton, NJ: Princeton University Press.

Google Scholar

Ramos-Fernandez, G., Aguilar, S. E. S., Schaffner, C. M., Vick, L. G., and Aureli, F. (2013). Site fidelity in space use by spider monkeys (Ateles geoffroyi) in the Yucatan Peninsula, Mexico. PLoS ONE 8:e62813. doi: 10.1371/journal.pone.0062813

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramos-Fernández, G., Aureli, F., Schaffner, C. M., and Vick, L. G. (2018). “Ecología, comportamiento y conservación de los monos araña (Ateles geoffroyi): 20 años de estudio en Punta Laguna, México,” in La Primatología en Latinoamérica 2, volume Tomo II, eds B. Urbani, M. Kowalewski, R. G. Teixeira da Cunha, S. de la Torre, and L. Cortés-Ortiz (Caracas: Instituto Venezolano de Investigaciones Científicas), 13.

Google Scholar

Ramos-Fernández, G., and Ayala-Orozco, B. (2003). “Population size and habitat use of spider monkeys at Punta Laguna, Mexico,” in Primates in Fragments, ed L. K. Marsh (Boston, MA: Springer US), 191–209. doi: 10.1007/978-1-4757-3770-7_13

CrossRef Full Text | Google Scholar

Rangel-Rivera, C. E. (2017). Uso del hábitat del mono araña (Ateles geoffroi) y cobertura vegetal en un Área Natural Protegida de la península de Yucatán, México (Maestría), Instituto Politécnico Nacional, Centro Interdisciplinario de Investigación para el Desarrollo Integral Regional, Oaxaca, Mexico.

Google Scholar

Ricotta, C., Avena, G., and Marchetti, M. (1999). The flaming sandpile: self-organized criticality and wildfires. Ecol. Model. 119, 73–77. doi: 10.1016/S0304-3800(99)00057-5

CrossRef Full Text | Google Scholar

Rivera-Núñez, T. A. (2014). Conservacionismo biológico o agencia humana en el manejo ambiental? El caso del Área de Protección de Flora y Fauna Otoch Ma'ax Yetel Kooh (Maestría), Centro de Investigación de Estudios Avanzados del Instituto Politécnico Nacional, Mérida, Mexico.

Robeva, R., and Murrugarra, D. (2016). The spruce budworm and forest: a qualitative comparison of ODE and Boolean models. Lett. Biomath. 3, 75–92. doi: 10.30707/LiB3.1Robeva

CrossRef Full Text | Google Scholar

Saadatpour, A., and Albert, R. (2013). Boolean modeling of biological regulatory networks: a methodology tutorial. Methods 62, 3–12. doi: 10.1016/j.ymeth.2012.10.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Schlüter, M., and Pahl-Wostl, C. (2007). Mechanisms of resilience in common-pool resource management systems: an agent-based model of water use in a river basin. Ecol. Soc. 12:4. doi: 10.5751/ES-02069-120204

CrossRef Full Text | Google Scholar

Schulze, J., Müller, B., Groeneveld, J., and Grimm, V. (2017). Agent-based modelling of social-ecological systems: achievements, challenges, and a way forward. J. Artif. Soc. Soc. Simul. 20:8. doi: 10.18564/jasss.3423

CrossRef Full Text | Google Scholar

SEGOB (2019). Plan Nacional de Desarrollo 2019–2024. Mexico: SEGOB. Available online at: https://www.dof.gob.mx/nota_detalle.php?codigo=5565599&fecha=12/07/2019

Google Scholar

SEMARNAT (2020). Áreas Naturales Protegidas decretadas (hectáreas). Mexico: SEMARNAT. Available online at: http://dgeiawf.semarnat.gob.mx:8080/ibi_apps/WFServlet?IBIF_ex=D3_BIODIV04_17&IBIC_user=dgeia_mce&IBIC_pass=dgeia_mce&NOMBREANIO=*

Google Scholar

SMN (2017). Información Climatológica–estación Cobá. Mexico: SMN. Available online at: https://smn.conagua.gob.mx/tools/RESOURCES/Estadistica/23012.pdf

SMN (2019). Información Climatológica por estado. Mexico: SMN. Available online at: https://smn.conagua.gob.mx/es/

Sánchez Sánchez, O., and Islebe, G. A. (1999). Hurricane Gilbert and structural changes in a tropical forest in south-eastern Mexico. Glob. Ecol. Biogeogr. 8, 29–38. doi: 10.1046/j.1365-2699.1999.00317.x

CrossRef Full Text | Google Scholar

Spaan, D. (2017). The effects of habitat disturbance on the populations of Geoffroy's spider monkeys in the Yucatan Peninsula (Doctorado), Universidad Veracruzana, Mexico.

Google Scholar

Standish, R. J., Hobbs, R. J., Mayfield, M. M., Bestelmeyer, B. T., Suding, K. N., Battaglia, L. L., et al. (2014). Resilience in ecology: abstraction, distraction, or where the action is? Biol. Conserv. 177, 43–51. doi: 10.1016/j.biocon.2014.06.008

CrossRef Full Text | Google Scholar

Su, Y., Hammond, J., Villamor, G. B., Grumbine, R. E., Xu, J., Hyde, K., et al. (2016). Tourism leads to wealth but increased vulnerability: a double-edged sword in Lijiang, South-West China. Water Int. 41, 682–697. doi: 10.1080/02508060.2016.1179523

CrossRef Full Text | Google Scholar

Tao, T. C., and Wall, G. (2009). Tourism as a sustainable livelihood strategy. Tourism Manage. 30, 90–98. doi: 10.1016/j.tourman.2008.03.009

CrossRef Full Text | Google Scholar

Thulstrup, A. W. (2015). Livelihood resilience and adaptive capacity: tracing changes in household access to capital in central Vietnam. World Dev. 74, 352–362. doi: 10.1016/j.worlddev.2015.05.019

CrossRef Full Text | Google Scholar

Toledo, V. M. (2013). “Indigenous peoples and biodiversity,” in Encyclopedia of Biodiversity, ed S. A. Levin, (Waltham, MA: Elsevier), 269–278. doi: 10.1016/B978-0-12-384719-5.00299-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Toledo, V. M., Ortiz-Espejel, B. F., Cortés, L., Moguel, P., and Ordoñez, M. D. J. (2003). The multiple use of tropical forests by indigenous peoples in Mexico: a case of adaptive management. Conserv. Ecol. 7:9. doi: 10.5751/ES-00524-070309

CrossRef Full Text | Google Scholar

Urquiza-Haas, T., Peres, C. A., and Dolman, P. M. (2009). Regional scale effects of human density and forest disturbance on large-bodied vertebrates throughout the Yucatán Peninsula, Mexico. Biol. Conserv. 142, 134–148. doi: 10.1016/j.biocon.2008.10.007

CrossRef Full Text | Google Scholar

van Roosmalen, M. G. M., and van Roosmalen, M. G. M. (1985). Habitat preferences, diet, feeding strategy and social organization of the black spider monkey [Ateles paniscus paniscus Linnaeus 1758] in Surinam. Acta Amazon. 15, 7–238. doi: 10.1590/1809-43921985155238

CrossRef Full Text | Google Scholar

Walker, B., Holling, C. S., Carpenter, S. R., and Kinzig, A. P. (2004). Resilience, adaptability and transformability in social-ecological systems. Ecol. Soc. 9:5. doi: 10.5751/ES-00650-090205

CrossRef Full Text | Google Scholar

Wilensky, U. (2006). NetLogo Fire Simple Extension 1 Model. Available online at: http://ccl.northwestern.edu/netlogo/models/FireSimpleExtension1