Detecting regime shifts in artificial ecosystems

Ecosystems are subjected to a range of perturbations that have the potential to induce relatively sharp transitions in states. These can be referred to as regime shifts or critical transitions. They may be driven by perturbations that vary over a wide range of spatial and temporal scales, from responses to deforestation within a small field to responses to the gradual increase of carbon dioxide in the Earth’s atmosphere. Here we investigate potential early warning signals that may presage regime shifts in model ecosystems. We hypothesise and model a relationship between biodiversity and community structure that influences ecosystem structure. We argue that Artificial Life methodologies have potential to make substantial contributions to efforts searching to predict large changes in ecosystems and other elements in the Earth system, as there is a recognised limitation in empirical data and ability to conduct experiments in the real-world. Consequently simulation and exploration of the low-level mechanisms that give rise to regime shifts in artificial in-silico ecosystems represents a useful line of enquiry.


Introduction
The relationship between complexity and stability is a well established topic in ecology (McCann 2000).Much of this discourse stems from May's seminal paper that found an inverse relationship between diversity and stability in simple model ecosystems (May 1973).As the number of linear connections between species increased, the probability that the ecosystem would be stable decreased.
Here we explore to what extent the complexity of artificial systems, and by extension real-world ecosystems, can be used to provide early warning signals of impending large changes in their states -sometimes referred to as regime shifts.That is, we do not explore the (causal) relationships between complexity and stability, rather we seek to leverage the maximum amount of information from very often complex systems in order to give an indication of how stable the system is to perturbations and how close they are to collapse.
Dynamical systems theory has been applied to a range of real-world systems in order to determine to what extent they exhibit bistability and the potential to undertake catastrophic fold bifurcations, and thus rapid changes from one attractor to another (Scheffer et al. 2001).One example of this phenomenon is the eutrophication of a lake.When a lake is subject to nutrient loading, changes may occur from a state of clear water with rich submerged vegetation to a turbid one without submerged vegetation and loss of animal biodiversity.At early stages, water clarity seems little affected by increased nutrient concentrations until a critical threshold is reached.Once this threshold is passed, the lake abruptly shifts to a turbid stable state.
Given the range of perturbations facing ecosystems and an increasing awareness of their contribution to human well being in the terms of the services they provide, a great deal of research effort is being devoted to formulate robust early warning signals (Scheffer et al. 2012;Lenton et al. 2012).Leading indicators of impending regime shifts have been identified that operate regardless of the underlying processes driving the change (Scheffer et al. 2009;Carpenter and Brock 2006;Guttal and Jayaprakash 2008).Some studies however, have identified the potential for false positives or even an absence of signals (Hastings and Wysham 2010), and studies of the drivers of regime shift reveal other signals in community composition (Scheffer and Nes 2007).The sooner a warning can be detected, the more time is available to act on it.The challenge is to find signals that warn at the earliest opportunity but do not provide false alarms.
Our ability to conduct experiments on real-world ecosystems is often very limited.For example, studies of regime shifts in the Lake Erhai, Yunnan Province, in China involved gathering then analysing sediment cores along with the integration of large amounts of socio-economic data (Wang et al. 2012).At this scale, it is effectively impossible to manipulate ecosystem variables let alone ensure control systems are available to compare experimental results.One way to increase our understanding in this area is to perform simulation experiments on artificial ecosystems.Such simulations can be understood as conceptual models that produce data appropriate for statistical analysis and so allow the development of proof of concepts and initial hypothesis testing.dance of keystone species.A keystone species is a species that has large effects on its environment and consequently plays a key role in the maintenance of the structure and integrity of its ecosystem.The loss of a keystone species will be expected to produce changes in a large number of other species such that the entire ecosystem may transition to a new state: a regime shift.
If we seek to undertake interventions and so avoid future regime shifts, waiting to observe changes in keystone species may be too late as the transition is already in progress.We propose instead monitoring properties of ecosystems that includes measures of complexity and structure.The novelty of our method is in using as detectors the species most vulnerable to extinction, which we term the 'ecosystem canaries', whose response to the perturbation precedes the approach to a critical tipping point.Provided they are correctly identified as sensitive to the driver of environmental change, they are not prone to false alarms.Provided they are not keystone species, their signal precedes those of the leading indicators for imminent critical transition.

Nestedness and biodiversity analyses
We hypothesise that biodiversity has a relationship with ecosystem structure that depends on important properties of the dynamical system.In this section we define ways to measure both the biodiversity and the structure of a community.Then we correlate these measurements over time to investigate the relationship between these quantities under different stress regimes.

Hill's biodiversity analyses
Species diversity has two main components: abundance of each species and species richness (Magurran and Magurran 1988).Hill's biodiversity indices incorporate these two components in a single value and are easy to interpret ecologically.The index value increases when either the number of species or the evenness of the community increases.In this study, we use Hill's biodiversity index N 2 , which is the inverse of the commonly used Simpson's diversity index S i=1 p 2 i (Simpson 1949).Here S is the number of species and p i is the proportion of the ith species in the community (Hill 1973).
Simpson's diversity index indicates the probability of finding an already observed individual by randomly choosing it from the community (Simpson 1949).Hence, the larger the value of Simpson index, the less diverse the community will be.Contrastingly, the higher the value of Hill's biodiversity index N 2 , the more diverse the community will be.

Nestedness analyses
If we assume that a community is an assemblage of metacommunities, where each meta-community contains a subset of the species found in richer meta-communities, then we can ask: how much information about the structure of the community can we obtain when we observe a single species?
The answer is potentially a surprising amount, if we also understand how nested the community is.
The concept of nestedness was first introduced as an explanation of insular faunal structure where species abundances decrease with distance of islands from a continent and species in distant islands are a subset of species in proximate ones (Patterson and Atmar 1986).Latterly, the concept has been applied extensively to terrestrial communities (Cutler 1991;Fischer and Lindenmayer 2005).
Nestedness is a characteristic of the interconnectedness of species in a community (Ulrich et al. 2008).If a community is completely nested, the community structure and composition may be predicted entirely from the presence of the least abundant species; if it is completely un-nested, its composition is unpredictable, as in the case of communities with a high flux of species.
In general, nestedness studies focus on the analysis of spatial data (Rodríguez-Gironés and Santamaría 2006), or spatial analyses replicated through time (Timi and Poulin 2008;Heino et al. 2009).However, nested temporal assemblages can occur when most species respond similarly to inter-annual variation in conditions.In contrast, assemblages might be non-nested when different sets of species occur in different years (Elmendorf and Harrison 2009).
Here we develop a novel approach in which measures of nestedness are computed for communities that change only temporally and not spatially (see Fig. 1).

Mainland
shows decreasing species abundance S with distance from mainland.Island S(ln) contains a subset of the species found in previous islands so S(ln) ∈ S(ln−1) ∈ ... ∈ S. Analogously, a given community S(tn) at t = tn can be seen as a subset of the community S(t0) at t = 0, such that S(tn) ∈ S(tn−1) ∈ ... ∈ S(t0).

ECAL 2013
One way to measure the nestedness of a community is via its incidence temperature (Patterson and Atmar 1986).This temperature ranges from 0 (completely nested) to 100 (completely un-nested), and provides a measure of the species richness across non-chronological times with respect to their incidence.The use of the term "temperature" can be understood from physical analogy: as the temperature of a collection of molecules increases then this additional energy will typically lead the molecules to occupy a larger volume.Knowing the position of a single molecule in a low temperature systems provides more information about the location of all the other molecules than higher temperature systems.
Originally, the nestedness temperature was presented as a measure of the degree of uncertainty in species extinction order and was linked in perfect analogy with the entropy of a system (Atmar and Patterson 1993).Recently, it has been argued that the nestedness temperature does not increase with extinction disorder (Almeida-Neto and Ulrich 2011) and therefore it cannot be explained as the inverse of the entropy (Almeida-Neto et al. 2007).However, nestedness temperature does gives a measure of the community structure dependent on the community distribution.In this sense, nestedness temperature may be measured in nondimensional units of entropy where a completely structured community would have the lowest entropy while completely disordered communities would have maximum entropy values.

Hypothesis
According to Montoya et al. (2006) the interaction between species in an ecosystem may be so complex as to be impossible to understand.Here we argue that, although interactions between species in a community are very complex, the underlying mechanism of their dynamics are reflected in the structure and biodiversity of the system.If we think of an ecosystem as a complex network where each species may (or may not) be associated to other species, the structure of this network will certainly affect the function of the system (Strogatz 2001).Inter-species connections may be trophic in the sense of who eats whom, or competitive in the sense of who competes for whom over a food or resource.They may also emerge from ecosystems engineering effects (Jones et al. 1994) where one species alters the environment for another.As we are motivated to produce a general hypothesis in which if the structure of a community at a given time is affected by its biodiversity, we deliberately do not prescribe on the type of interactions among species.Instead, we propose that a correlation between nestedness temperature and Hills Biodiversity index may exist.
We assume that a community is composed of different types of species among which we can find keystone species, interacting species and "canary species".The keystone species are those with a large effect on the system, the interacting species will also be linked to other species in the community, while the presence or absence of canary species will have little impact.Canary species are also the most vulnerable to extinction.Consequently they will be rare in the community and will have little or no interaction with other species (see Fig. 2).For a potentially broad range of ecosystems that are characterised as having well-mixed or homogeneous environmental variables, we propose three key stages towards a regime shift (see Fig. 3).

Keystones Interacting Canaries
(a) Long before the regime shift: Temperature and biodiversity both fluctuate without obvious trend in biodiversity and temperature.However there is a weakly positive correlation (i.e., more biodiversity correlates with more disorder).This represents a healthy community that adds and loses canary species more or less at random.
(b) Prior to the regime shift up to its cusp: Temperature and biodiversity fluctuate without obvious trends in either, but now with a strongly negative correlation (i.e., more biodiversity correlates with less disorder).Upward fluctuations of biodiversity add competitively dominant (keystone) and fast-fugitive (weedy) species at the expense of canaries, thereby also raising order.In effect, the out-ofphase fluctuations of biodiversity and temperature ratchet out the canaries and replace them with more strongly connected species, reflecting a tightening of the community in response to stress.
(c) During and after the regime shift: Temperature falls dramatically with biodiversity in a strongly positive correlation, as keystone species are lost, leaving only the few most abundant and robust (weedy) species that have always been present.
These hypotheses and assumptions are consistent with empirical results obtained from Lake Erhai, Yunnan Province of China.A dataset for changes in a diatom community over time was obtained from a 63-cm sediment core in Lake Erhai, representing about 500 years of sedimentation prior to a critical transition in 2001, and 8 years posttransition.We developed a temporal analysis of the diatom community to obtain information about the community composition and structure before and after a critical transition that took place in 2001, according to results found by (Wang et al. 2012).
Before the critical transition, the community shows relatively high temperature (low levels of nestedness), though decreasing gradually towards the tipping point.During and after the critical transition, temperature decreases drastically.We correlated nestedness temperature with community biodiversity, finding that the sign of the correlation switched from positive to negative at about 50 years prior to the tipping point, and then back to positive immediately after the tipping point.We interpret these changes as a potential signal of ecosystem stress, which leads to the community tightening up as it loses first the canary species and finally the keystone species as it goes over the critical transition.
These empirical results provide a benchmark for the generation and analysis of the artificial ecosystems presented below.With limited capacity for further progress through analyses of empirical datasets, because of the financial and time costs involved in obtaining them, we perform numerical simulations of communities similar to the ones found empirically, and we test our assumptions on them.Such simulations can be understood as an initial evaluation of our hypothesis.A next step would be the development of more detailed agent based models that would explicitly capture the important processes that we identified in (a), (b) and (c) above.

Simulation experiments
We generated artificial ecosystem matrices, that simulate community distributions of species under the influence of different types of stressors.The objectives of these simulations are: (1) to test for variations in the nestedness temperature in response to stressors, and (2) to identify relationships between nestedness and biodiversity through time.We anticipate that the results of these numerical models may point towards a robust early warning signal for regime shifts in community composition.

Calculating nestedness
Nestedness temperature measures how much the incidence matrix departs from perfect nestedness.To calculate this metric, the list of species present in a series of times is summarised in an incidence matrix of presence-absence.The rows and columns of the incidence matrix are reordered so that nestedness is maximised, using the algorithm proposed by Rodríguez-Gironés and Santamaría (2006).The matrix is re-arranged to show species presences on the top left corner and the absences away from the top left corner.This creates a matrix where the columns rank species rarity (increasing from left to right) and the rows rank species richness (increasing bottom to top).Then an isocline of perfect nestedness is calculated to show the expected distribution of presences if the matrix were perfectly nested.
Absences to the top and left of the isocline are defined as unexpected, and so are presences below and to the right (see Fig. 4).The matrix nestedness temperature is calculated as the sum of squared deviations from the isocline of unexpected presences and absences divided by the maximum value possible for the matrix, multiplied by 100.Thus, the temperature is a non-dimensional index measuring how much the matrix departs from the perfectly nested state (Rodríguez-Gironés and Santamaría 2006;Almeida-Neto et al. 2007).

Generating power-law population distributions
Many populations follow power-law frequency distributions with respect to their abundance across time (Mitzenmacher 2004;Allen et al. 2001).Accordingly, we generated power-law frequency distributions of theoretical communities which we used as the basis of numerical experiments.We generated a matrix with 150 rows and 150 columns, where the M (i, j) represents the species in column j at time i.Each row is an independently generated power-law community distribution of 150 species.Once a core matrix was generated, we applied a number of forcing functions to it.By doing this we introduced what would be analogous to stressors of change in the environment (e.g.changes in nutrients, salinity, pH).The proposed functions had an immediate effect in the change of nestedness temperature across the nested blocks over the temporal scale.Equally, these changes in the community structure, affected Hill's biodiversity index.The main objective of these numerical simulations was to obtain correlations between these quantities to compare to the empirical data on diatoms from Lake Erhai.

Core matrix analysis
Using the method proposed in the previous section, we generated 100 matrices with power law distributions.Each core matrix was transformed into a presence-absence matrix from which we found an average nestedness temperature of T ≈ 87.57. Figure 5 shows the nestedness temperature for a matrix with these characteristics.The high temperature of the community matrix does not show any type of structure.However, as we will see in the next section, when the system was subjected to external perturbations, the community acquired structure and the nestedness temperature decreased.

Forcing functions
We have limited knowledge of the type of stressors acting on a real ecosystem.Consequently, we formulated a series of forcing functions that reduce the number of species in the community proportional to the value of parameter, z.The different forcing functions affect the abundance of species in richness Species rarity Figure 5: An example of a high incidence temperature system that shows little structure.A dataset with 150 species following a power law distribution for 150 times depicts an incidence temperature of 87.5.The continuous black line represents the isocline for a perfectly nested structure.The empty squares above the isocline denote surprising absences, while the red squares below the isocline denote surprising presences.different ways.For example one forcing function will preferentially remove species with high abundance, another will reduce species with low abundance.We are motivated to do this as we seek to understand how the correlation between nestedness temperature and biodiversity changes under different perturbations and stressors.
Max Eliminates the species with the highest number of individuals in each row: − z.These species are located in the peak of the distribution, and must exceed because the number of species with such high abundances is very small.
Min Eliminates the species with the lowest number of individuals in each row: These species are located in the tail of the distribution and comprise a large percentage of the community.
Middle Eliminates the most abundant and the least abundant species, retaining those in the middle of the distribution: Outflux Decreases the number of individuals of each cell in the matrix by a percentage z given by a fixed parameter In-Out This function simulates the influx and outflow of individuals f (z) = M ij − z + i.The population in each cell increases with time (row number i) and decreases with the size of z.
Here we have assumed that M i,K corresponds to the ith row for all the columns of the matrix M , M ij correspond ECAL -General Track

629
ECAL 2013 to the cell (i, j) of the matrix and z is a scalar parameter of fixed value.These functions are not an exhaustive selection; they give general indications of what responses stressors might induce in an originally well-mixed community.
In order to calculate how nestedness temperature changed over time in response to the different forcing functions, each matrix was divided in 135 consecutive nested subsets, where each subset had 15 rows (i.e.rows 1 to 15, then rows 2-16, 3-17, and so on).The nestedness temperature was then obtained for each subset.

Results
The nestedness temperature varied with respect to the forcing function and the size of the forcing parameter z (see Fig. 6).In the cases of functions Max, Min, Middle and Outflux, the temperature decreased as the intensity of the forcing function increased (as z increases), with differently decreasing curves depending on the forcing function.This is expected since as the number of species decreases loosing the canary species, the community becomes more nested.In the case of function In-Out, we observed that the nestedness temperature increased from 0 to 34.57.The reason for this is that the competing in-and out-fluxes balance each other as the parameter z increases.When z = 0 the population in each row increases with the row number, e.g. each species in row i will increase in i% adding more individuals to the community.As z increases, individuals leave the row community i, balancing out the fluxes into and out of the community.Consequently, when z ≥ 100, the population only decreases until it eventually vanishes.
As a preliminary analysis, we choose three matrices from the whole set displayed in Fig. 6 and investigate the possi-ble correlations between the community structure and biodiversity.We analysed relatively nested communities with nestedness temperature of about 25 units for three forcing functions in order to determine how a drop/increase in temperature and Hill's biodiversity changes the correlation between these two quantities.The results for these analyses are summarised in Table1.
We chose functions Min and Middle since they best account for a linear loss of canary species, while In-Out shows drastic changes in nestedness temperature.Our analysis shows that a linear increase or decrease in nestedness temperature maintains the sign of the correlation between temperature and biodiversity.This implies that changes in community structure will maintain the same type of dynamics on the biodiversity.Conversely, we found that a sharp decrease or increase in nestedness temperature changes the correlation sign.The best example for this change of sign in the correlation is given by function In-Out displayed in Fig. 7.A sharp change of temperature and Hill's diversity index results in a change of correlation sign from negative (while the biodiversity starts increasing the nestedness remains more or less constant due to the species flow into and out of the system) to positive (a sharp increase in biodiversity results in a higher temperature due to the large influx of individuals into the community), then a sharp drop in the temperature and Hill's diversity index results in a positive correlation (as biodiversity decreases sharply the temperature decreases as well as only few structured species remain).Notice that after this decrease in temperature and biodiversity, the correlation remains positive.This is due to the continued decrease of biodiversity and temperature after the shift.Notice also that, when the drop of the temperature is abrupt, the correlation is stronger.1.

Discussion
Our study of the biodiversity and nestedness of an artificial ecosystem aims to identify sensitive detectors of changes in environmental processes that drive regime shifts.Our hypothesis assumes that "canary species" (those most susceptible to extinction) function as detectors of drivers of change, preceding the approach to a critical tipping point.
We generated artificial populations with power law distributions and subjected them to a number of drivers of change.We correlated nestedness with biodiversity and established which conditions of external stress produced which correlations.We found that large changes in biodiversity and nestedness can lead to significant changes in the sign and magnitude of the correlation of the two measures.Such a change is qualitatively similar to that observed in a real-world lake ecosystem.We hypothesise that the change in correlations in the lake system are produced by the loss of ecosystem canary species in response to external perturbations.As these vulnerable species are lost the community 'tightens'.Increasing the intensity of perturbation would lead to further changes such that keystone species were affected.This would not only change the correlation between biodiversity and nestedness, but also lead to an imminent collapse in the systems as major ecosystem level properties and structure would rest.
Our analysis presented here was performed for a small number of systems.In order to further explore the utility of an early warning signal of an impending regime shift that is generated by the correlation of biodiversity and nestedness we propose the development of a series of agent based models that would build on this initial work.These models would allow the inclusion of traits such as growth rate and competitive ability.In doing so we would be able to explore the interactions between competitively dominant but slow growing keystone species, fast growing but competitively week species and canary species that are both poor competitors and slow growing.If our initial assumptions are correct, then what would emerge in response to progressive perturbations on the population would be a robust indicator of a regime shift.

Figure 2 :
Figure 2: A community of species represented by their interactions with each other.Red dots represent keystone species, which are the species with the most links in the community.Blue and black dots are interacting species with fewer links to other species.Green dots are canary species, which are relatively isolated from the community and have very few links with other species.The lines between dots correspond to the links between species.Modified fromKrebs (2012).

Figure 3 :
Figure 3: Increase of phosphorus influx into a lake drives the system from an oligotrophic state (clear water) towards a regime shift characterised by a eutrophic state (turbid water) passing through three stages (a) -(c).We propose that a negative correlation of nestedness temperature to biodiversity arises in stage (b) because rising stress on the community (phosphorus influx in this case) causes species accumulation rather than turnover with the loss of canaries and consolidation of strong competitors.During and after the regime shift (stage (c)), both the biodiversity and the nestedness temperature decrease and the system becomes a compact, nested community with very few species left.This produces a positive correlation.

Figure 4 :
Figure4: Incidence matrix of m species identities (ranked from most to least often present), in n samples (ranked from most to least species rich, bottom to top).The isocline curve corresponds to the theoretically perfect nested structure of the system.All the red squares are presences and white squares denote absences.Temperature is calculated from the sum of all squared deviations of unexpected presences/absences u ij = (dij/Dij) 2 .Modified from Oksanen 2012.

Figure 6 :
Figure6: Changes in nestedness temperature of an artificial community of 150 species, over 150 time steps in response to increasing intensity (z) of each of six stressors.The stressors that eliminate the rare species first (e.g., Min-Outflux) tend to reduce nestedness temperature most rapidly, reflecting an increased ordering and predictability in species composition.Function In-Out simulates a continuous flux of species through the community, which raises nestedness temperature from zero, reflecting the reduced predictability of composition induced by the turnover of species.

Figure 7 :
Figure 7: Nestedness temperature and running average of Hill'sdiversity index N2 for the numerically generated matrix In-Out with a stress intensity of z = 30.Both, the temperature and the Hill's index increases and decreases very drastically between row 8 and row 30.The arrows between vertical dotted lines show the regions analysed to obtain the correlation between these two nonsmoothed quantities.These correlations are summarised in Table1.

Table 1 :
Trends in nestedness temperature (T) and Hill's biodiversity (H) and their correlation coefficient (r and significance P).Trends are + or ++ for weakly or strongly positive, − or −− for weakly or strongly negative, n.t. for no trend.Correlations between nestedness temperature and Hill's biodiversity index show that linear changes maintain the sign of the correlation, while drastic changes result in a correlation shift.