Nonequilibrium Dynamics in Conservation Biology: Scales, Attractors and Critical Points

1. Introduction

The future of our biosphere is likely to be decided within this century. Potential paths to sustainable growth are to some extent in conflict with ongoing human population growth, overexploitation of resources and the resulting impacts on climate. Different studies suggest that decisive actions might be required within the next few decades (Solé and Levin, 2022). Because of its relevance, biodiversity emerges here as the cornerstone of any relevant path to Earth stewardship (Chapin et al., 2022). What are the right and effective ways of dealing with conservation? Alternatives range over a very vast landscape of possibilities, from the expansion of protected areas to ecosystem bioengineering. No matter what kind of path is followed, we need to deal with the complex nature of ecosystems and the biosphere (Levin ,1998; Brown et al., 2002; Green et al., 2005; Solé and Bascompte, 2012; Oro and Martinez-Abrain, 2023; Riva et al., 2023). In this context, an obvious problem is the presence of multiple scales that we naturally appreciate when looking at ecosystems organization and how they appear to be somewhat disconnected. Information concerning genomic-level description of species or communities might be very useful to map genetic diversity, but will tell little about population fluctuations and their responses to environmental challenges. Species-level properties, on the other hand, will not be enough to explain the spatiotemporal properties of communities (Levin, 1992; Levin et al., 1997). One obvious example is pattern formation in drylands, where large-scale patterns (sometimes highly regular) arise from short-scale interactions among individual plants.

Such scales are, all of them, relevant to our understanding of the biosphere and the threads associated to climate change (Peñuelas et al., 2013). Genomics can play a leading role in conservation efforts (Breed et al., 2019) while data-driven monitoring of global ecosystems by means of remote sensing can help define a resilience index (Lenton et al., 2022). However, against the traditional views of physics and biology, each scale is not reducible (in mechanistic terms) to the properties of the lower levels. This is at the core of the concept of emergence, which percolates complexity in living systems (Solé and Bascompte, 2012; Levin, 1999). At the heart of this complexity lies a common thread: the presence of nonlinear thresholds that can trigger sudden and profound shifts among alternative states within ecosystems. These tipping points represent critical junctures where small changes in environmental conditions lead to disproportionately significant and irreversible consequences (Scheffer et al., 2001; Lenton et al., 2008; Scheffer, 2009; Dakos et al., 2019). Recognising the existence and implications of these tipping points demands a paradigm shift in our perception of risk and calls for innovative and adaptable management strategies (Suding and Hobbs, 2009).

Different factors influence the preservation of healthy communities and how we approach conservation efforts. From the local level, challenges arise from habitat loss, fragmentation, and the disruption of ecological interactions. Concurrently, at larger scales, the intricate web of ecological complexities and the influence of extreme events further compound the conservation landscape. Understanding the linkages among these diverse scales is crucial to crafting effective conservation practices.

In this paper I present several key ideas and case studies that illustrate the presence of emergent phenomena in ecosystems with relevant implications for biological conservation. This is not a systematic overview nor a collection of examples, some of which are shown in Figure 1a–f. Instead, our aim is to clarify the nature of some universal properties that appear to be common to a broad range of situations. As part of the analysis, I want to show why (and how) very simple models can be helpful to answer crucial questions about the nature of nonlinear behaviour and how they pervade habitat degradation, extinction or defaunation dynamics. Three main topics will be explored in terms of the simple mathematical models that can account for their dynamics. All have important consequences to our potential to preserve biodiversity and define conservation strategies. Specifically, I have chosen these: (a) habitat loss and fragmentation, (b) facilitation and (c) criticality and forest fire dynamics. These three case studies, and the mathematical models used to make sense of them, illustrate the need of theory as a robust approach to the different scales involved and their origins.

A common feature of complex systems is the presence of negative and positive feedbacks (DeAangelis and Waterhouse, 1987). Both kind of feedbacks (or the combination of them) introduce the nonlinearities that determine the stability of the system and the emergence of complex dynamics. Typically, highly stable states allow fluctuations to be kept at bay thus allowing a robust homeostatic state despite the noise. This was one of the key contributions of cybernetics: how feedback control creates stability. The former have a rather different role. They involve interactions that enhance growth and as such are amplifiers of fluctuations. Their relevance in science is illustrated by their presence in different fields, from neuroscience to Earth sciences.

Negative feedbacks are easily identified in the classical models of population growth. This is the case of the so called logistic growth model, where a single-species population expands by reproduction while competition for available, finite resources, limits further growth (Case,1999). It is such an implicit limitation what describes the negative feedback: as soon as population size starts to be comparable with the maximum population size achievable, further growth starts to be inhibited and a slowdown in growth rates will be inevitable. The presence of negative feedback is behind part of the ability of ecosystems to respond homeostatically to external changes. But the interaction of growth and death can also generate thresholds that separate two "phases" associates to community persistence and extinction, respectively.

Positive feedbacks, on the other hand, are a natural consequence of mutualistic interactions. As will be shown below, these positive exchanges are responsible for the presence of multiple alternative states (May, 1977). Oftentimes, a stable population can suddenly jump into a degraded state under small variations of parameters (Scheffer et al., 2001). Finally, the spatial dimension naturally introduces further levels of complexity and new emergent properties (Levin, 1992; Rietkerk and Van de Koppel, 2008). Positive and negative exchanges can, under dispersal, generate stable self-organised spatial structures far beyond the short scale of local interactions (Margalef, 1968; 1979; Nicolis and Prigogine, 1989; Bascompte and Solé, 1995; Solé and Bascompte, 1998; Murray, 2003) but also highly dynamical states, such as forest fires, which act as global herbivores (Bond and Keeley, 2005). All these examples share a common feature that has significant implications to conservation and engineering practices. On the one hand, universal phenomena pervade their nature and deeply constrain their emergence. On the other, this universality is what makes possible to predict future scenarios of change as well as the potential impacts of interventions.

2. Ecological Attractors

As Robert May aptly noted, simple mathematical models can exhibit highly intricate dynamics (May, 1976). As will be shown in the following examples, the nonlinear dynamics of complex ecosystems might be properly captured by seemingly simplistic models. Prior to delving into our illustrative examples, let us provide a concise overview of essential definitions and properties that will be utilised here and will help to follow the logic of simple models and their unexpected predictions. Readers already acquainted with these concepts can readily proceed to the subsequent section.

The most elementary class of mathematical models encompasses one-dimensional population models, where the fraction of individuals denoted as x follows a dynamic equation:

$$\frac{dx}{dt}=f_{\mu }\left(x\right)$$

(1)

In this equation, the rate of change of our population, denoted by the derivative $dx/dt$, is determined by a population-dependent function $f_{\mu }\left(x\right)$. The subscript $\mu$ serves as a reminder that all processes typically rely on one or more parameters. Occasionally, continuous variations in these parameters can give rise to discontinuous changes in x.

Two fundamental concepts, straightforward yet pivotal to our analysis, will be employed below. The first is the notion of an equilibrium (or fixed) point. We designate a fixed point $x^{*}$ as a population value for which

$${\left(\frac{dx}{dt}\right)}_{x^{*}}=0$$

(2)

The underlying idea is that, when the population reaches this value, it remains unchanged. However, the nature of this equilibrium can vary widely. This distinction is often depicted through a mechanical analogy, as illustrated in Figure 1g–j. Equilibrium points can be likened to marbles positioned in an abstract landscape that somewhat represents potential available states of the given ecosystem. Under a given set of conditions (where parameters are fixed) we will observe the state of the community (defined by some variables) as located at the bottom of one valley ($S_1$ in Figure 1g). As parameters change (such as increased aridity, acidity or reduced precipitation) the landscape is deformed and the original state becomes unstable and a new alternative state is reached in a catastrophic way ($I\to j$).

In a general scenario, with several dimensions involved, the underlying landscape will be much more complex, as illustrated in (Figure 1k). Here several stable attractors are shown (1,2,3). A given initial condition within the basin of a given valley (such as 4) will roll down towards the stable state 1. Moreover, we also display another state (5) that could be stable: if we place a marble here, it might be that it remains there. But any slight perturbation will break the balance either towards 1 or 2. This constitutes an unstable fixed point. A distinct yet noteworthy (and vital, as explained later) case is that of marginal fixed points. These arise when perturbing $x^{*}$ with a slight displacement simply transitions the state to another (analogous to a flat landscape). Determining whether a fixed point is stable, unstable, or marginal is a key consideration.

There are two ways in which we can answer this question. One is to determine the stability of $x^{*}$ is to find the sign of the derivative of the growth function $f_{\mu }\left(x\right)$ against x and evaluated at the given $x^{*}$. In other words, calculate the parameter:

$$\lambda \left(x^{*}\right)={\left(\frac{d{f}_{\mu }\left(x\right)}{dx}\right)}_{x^{*}}$$

(3)

which is a real number. It can be shown (see appendix) that a given, small deviation y from equilibrium (i.e., pushing the system slightly away from $x^{*}$ to a new state $x^{*}+y$) will evolve in time following

$$y\left(t\right)\sim e^{\lambda t}$$

(4)

if $\lambda <0$, the fixed point will be stable (since $x\to 0$) while $\lambda >0$ indicates an unstable (since y will grow). Marginal points are characterised by a zero value.

One final component of this short review of simple dynamical systems concerns the mechanical analogy described above. Although the qualitative picture is a compelling one, a rigorous mathematical definition of this landscape would be desirable. Although there is no general result for dynamical systems with dimension two or larger, it is possible to write down a well-defined expression for the so called potential function $\mathcal{V}\left(x\right)$. Specifically, the system is said to derive from a potential function $\mathcal{V}\left(S\right)$ (Strogatz, 1994; Solé, 2011) if we can write the dynamical system as:

$$\frac{dx}{dt}=-\frac{d\mathcal{V}}{dx}$$

(5)

From the previous equation, it can be shown (see appendix) that the potential reads:

$$\mathcal{V}\left(x\right)=-{\int }^x{f}_{\mu }\left(x\right)dx.$$

(6)

One of the main messages gathered from the study of emergent phenomena is that complex dynamics can be captured by very simple models. More importantly, there are several fundamental properties that behave in universal ways. Many of the key results developed within theoretical models of species extinction are connected with (or equivalent to) well-established models of the physics of strongly interacting systems (Chopard and Droz, 1998). This is in fact one of the most relevant domains of statistical physics, where the concept of universality is central to our understanding of complexity (Sethna, 2021). In a nutshell, nonlinear systems with strongly interacting parts will typically exhibit low-dimensional dynamics as a few degrees of freedom dominate the macroscopic patterns (Haken, 1975, 1980; Nicolis and Prigogine, 1989; Gutenkunst et al., 2007). This is the basis of our success from using simple mathematical models displaying emergent behaviour.

3. Habitat Loss and Extinction Thresholds

Habitat destruction has been long recognised as one of the most impactful sources of biodiversity loss (Laurance, 2010; Haddad et al., 2015). It provides a first example of the relevance of simple models allowing to understand the presence of threshold phenomena. As we will see, despite the apparently simple nature of the process, some of the implications of this class of models is sometimes unexpected and contrary to linear thinking.

A very simple approach to habitat loss is provided by an extension of Levins metapopulation model (Levins, 1969). In the original formulation, a given finite habitat is considered as made of a collection of connected patches. The connectivity is defined by accessibility through dispersal: two local patches are connected if they can access each other by means of the movement of individuals. One given patch might contain one or many individuals, or be empty. If $1\ge x\ge 0$ is the fraction of occupied habitat, Levins model reads:

$$\frac{dx}{dt}=cx(1-x)-ex.$$

(7)

Here a set $\mu =\{c,e\}$ of two parameters are involved. One is the probability c that a given empty patch is colonised by an occupied one, while e is the probability that a site experiences an extinction event. We can easily describe these two conflicting rates in terms of probabilities of colonisation and decay on a given abstract habitat, as shown in Figure 2a, where occupied and empty patches are indicated as filled and empty squares, respectively. A crucial consequence of this model Is that there are thresholds to extinction. Since the fixed points obtained from $dx/dt=0$ are $x_0^{*}=0$ (extinction) and $x_1^{*}=1-e/c$, it is easy to see that the second (nontrivial) point $x_1^{*}>0$ only if $e<c$. How is this model extended to incorporate habitat loss?

Habitat loss can be easily introduced by means of a single parameter D that gives the fraction of patches that have been destroyed and cannot be colonised. In its simplest form, we just need to extract D from the available habitat, which now will be $1-D-x$. It is easy to see that this corresponds to a slightly different equation, namely:

$$\frac{dx}{dt}=cx(1-D-x)-ex$$

(8)

The new fixed points are now $x^{*}=0$ and

$$x^{*}=1-D-\frac{e}{c}$$

(9)

It can be shown that extinction will occur if $\lambda (x^{*}=0)<0$ which gives a threshold:

$$D>D_c=1-\frac{e}{c}$$

(10)

which is our first key, relevant result. It is summarised in Figure 2b, where we plot the so called phase space of our model on the $(D,C)$ plane. The shaded area corresponds to the parameter combinations allowing for a stable metapopulation, whereas the white area is the "phase" where extinction is expected.

When we plot $x^{*}$ against both c and D, we can see that, for a fixed colonisation, $x^{*}$ decreases from $x^{*}(D=0)=1-e/c$ in a linear fashion (white line) until it reaches extinction at the corresponding threshold value of destruction. And here comes a fundamental result. Once the critical threshold has been reached, there is still a potential available habitat given by $1-D$. The linear view would have predicted that a smaller, but nevertheless non-zero fraction of patches should be occupied. However, and this is a consequence of the original Levins model, the species-specific parameters c and e are also important and impose boundaries to metapopulation persistence. The take-home message is clear: even if some significant part of a landscape remains suitable for a given species, once we cross the minimal requirements imposed by the previous model, persistence is not an option. Mounting evidence has shown how habitat loss can trigger the extinction of endangered species. Examples include, among others, the Northern Spotted Owl (Lande, 1987) that has become a canonical example of how a few key parameters associated to dispersal and demography allow defining the critical conditions for survival.

Another important property of these thresholds is the presence of a marked increase of fluctuations as we approach $D_c$. Some further analysis of the model can help to show that this phenomenon must be expected. Consider our original Equation (8) for our metapopulation model and suppose that we start from a very small metapopulation size, i.e., $x\ll 1$. In that case, the model becomes much simpler:

$${\left(\frac{dx}{dt}\right)}_{x\ll 1}\sim c\left(D_c-D\right)x$$

(11)

which can be solved to give a dynamical Equation of x close to criticality:

$$x\left(t\right)=x_0{e}^{c\left(D_c-D\right)t}$$

(12)

How long it takes, from a given $x_0$, to reach extinction? The time $T(ϵ,D)$ needed to achieve this can be found from the last equation. By replacing $t=T(ϵ,D)$ and the final state $x\left(T\right(ϵ,D\left)\right)=ϵ$ we obtain

$$T(ϵ,D)=\frac{ln(ϵ/x_0)}{c\left(D_c-D\right)}\sim {(D_c-D)}^{-1}$$

(13)

This result has an enormous implications. It tells us that the relaxation time to equilibrium diverges close to the critical point. In practical terms, it means that, because the marble is rolling now on a rather flat valley, population size can fluctuate with a high variance. The exact number of individuals becomes less predictable while the time required to achieve extinction can be very long. This is in fact a well-known result from the statistical physics of phase transitions and is known as critical slowing down (Solé, 2011; Ganghi et al., 1998). It is a universal property of systems that display continuous phase transitions where conflicting tendencies between order and disorder are present.

What are the implications of these long transients? One is that we can actually use this feature to detect so called warning signals: as we get close to a critical points where some undesirable shift can occur, we should observe that population time series should display increasing variance (Scheffer et al., 2009; Clement and Ozgui, 2018). Another consequence (to be discussed in the next section) is that high variances might translate into confusing results concerning the outcome of habitat loss. Although this process should in principle correlate with biodiversity decays, some studies have pointed out that the opposite might sometimes be observed. Here long-term fluctuations might obscure the predicted decay and be a byproduct of the slowing down close to threshold points. In this context, generalisations of the previous model can be derived, including a S-species competitive community. This is the case for the following model (Solé et al., 2004):

$$\frac{d{x}_i}{dt}=c_i{x}_i\left(1-D-\sum _{j=1}^S{x}_j\right)-\varphi \left(c_i\right)x_i$$

(14)

with $I=1,...,S$ and $\varphi \left(c_i\right)$ introduces a trade-off between colonisation and extinction so that $c/e$ is the same for all species. Using $e_i=\varphi \left(c_i\right)=\alpha c_i$, this deterministic, quasi-neutral model behaves in similar ways as the single-species scenario, with all species getting extinct for $D_c=1-\alpha$. However, the spatially-extended counterpart displays a maximum of species diversity, as measured by Shannon’s entropy, right on that point. This is an important problem that needs to be investigated. Additionally, other models involving a ranked structure of competitors shows that extinction can occur in a delayed manner: even when given thresholds have already been crossed, it can take a long interval to observe the final demise of a given species (Tilman et al., 1994). As pointed out by some researchers, this poses a challenge to biological conservation policies (Kussaari et al., 2009).

A final but not less important phenomenon needs to be considered. Our model so far is a mean field model: we completely ignore the existence of local interactions. This is not a minor matter. The introduction of space in an explicit fashion can generate novel phenomena, from stable regular patterns to complex waves (Murray, 2003). These are emergent properties: they are system-level phenomena out of local, microscopic interactions, but cannot be reduced to the lower-level components (Artime and De Domenico, 2022).

It was shown in the early times of spatial ecology that introducing space in an explicit manner can deeply modify the expectations from well-mixed systems (Durrett and Levin, 1994; Bascompte and Solé, 1995). In Figure 3a–c we show what happens when space is taken into account. To build such a spatial model, as the one outlined in Figure 2a for the Levins model, we consider a set of probabilistic rules of patch colonisation and extinction (see Bascompte and Solé 1997 for a detailed explanation). In a nutshell, a given empty site $\mathit{i}\in \Omega$ can be colonised, with probability c, by any of the occupied sites in the neighbourhood ${\Gamma }_i$ (Figure 2a, zoomed area). Similarly, an occupied site can get extinct with a probability e.

The metapopulation occupancy resulting from these set of rules is shown in Figure 3a for different habitat destruction values using a 100 lattice (gray curve) and compared with the theoretical prediction from Levins model (black line). A marked difference is observed: the occupied habitat becomes much less than the one expected from the well-mixed one. What is going on? We will be back to this when we explore the problem of percolation, but let us advance that an emergent phenomenon occurs at a given destruction threshold, indicated in Figure 3b, where we plots the size of the largest connected set of patches for each value of D. Despite the smooth character of the increased habitat loss (depicted in Figure 3c for three D values), a phase transition takes place at $D\approx 0.4$. The system becomes fragmented afterwards, becoming broken in multiple patches. Because of this phenomenon, different isolated clusters of healthy habitat will sustain disconnected metapopulations, that will be more and more affected by fluctuations. This combination of smaller domains and persistence thresholds translates into the observed displacement of the extinction threshold.

Field observations and theoretical models have been instrumental in developing a conservation strategy to fight the effects of habitat fragmentation. In order to recover the connectivity required to increase the effective habitat size (as would be defined by the largest connected habitat) is to design biological corridors (Anderson and Jenkins, 2006; Pimm and Jenkins, 2019). Using models allow to address the problem of optimal reserve design, often using concepts of graph theory (Williams et al., 2005; Niebuhr et al., 2015; Rocha et al., 2021) particularly within the context of future climate change scenarios (McGuire et al., 2016).

4. Facilitation and Tipping Points

In the previous section we have considered models where species responses are somewhat isolated from the rest of the community. This is the case of Levins model and its extensions to habitat loss scenarios. In those models, the predicted critical point that separates a stable metapopulation from its demise applies to a case study where only one species is present. Single-species metapopulation models that address the conditions for extinction in a given landscape have received much attention (Hanski, 1999). However, this cannot be a satisfactory picture if we deal with species interactions, which can introduce new kinds of nonlinearities. To explore their impact, let us now consider a simple, but extremely important phenomenon known as facilitation, that is very relevant in dryland ecosystems (Rietkerk and van de Koppel 1997, Scanlon et al. 2007, Kéfi et al., 2007a-b, Solé 2007) where warming, steady declines in rainfall and increased grazing are likely to promote a sudden shift from a vegetated ("Green") to a desert state (Solé, 2017). The analysis and modelling of spatial patterning in semiarid habitats consistently supports the suggestion that rapid shifts might occur in a near future (Kéfi et al. 2007a-b, Solé 2011).

Facilitation refers to positive (cooperative) interactions between individuals that confer benefits upon at least one of the participating entities (Brooker et al. 2008). A particularly relevant instance encompasses instances where one organism enhances the local environment, rendering it more conducive for another, irrespective of whether they belong to the same species or not. This affirmative influence can manifest directly, such as through shading mechanisms that alleviate water or nutrient stress, or indirectly, like the removal of competitors or the deterrence of predators (Bruno et al. 2003).

In order to push forward the idea that the addition of qualitative properties to a given ecological model leads to new emergent dynamics, I will use as a case study a straightforward extension of the previous habitat fragmentation models including facilitation in its simplest form (Sardanyes et al., 2019). With this simple model, we will be able to address several important questions: What implications does habitat destruction bear on an ecological system that incorporates facilitation? How does this impact interact with the nonlinear attributes associated with facilitation?

Let us start with some complexity and see how can can reduce an initial multispecies model to a single-variable model. Consider the available space as a finite resource that we assume (for convenience) normalized to one. Let us also assume that part of this habitat has been lost due to habitat destruction. Those affected areas (either by the presence of barren soil or contaminated patches) are no longer amenable for colonization. Consider a population of S species where abundances are indicated as $x_1,x_2,...,x_S$. Let us assume that they display both competition and cooperation. The model now involves multiple dimensions (species diversity) and reads:

$$\frac{d{x}_i}{dt}=\sum _{j=1}^S{c}_{ij}{x}_i{x}_j\left(1-D-\sum _{w=1}^S{x}_k\right)-e_i{x}_i,$$

(15)

Here $c_{ij}$ stands for the mutual help (facilitation) that a given plant species j provides to the individuals belonging to species i. The cooperaive nature of facilitation is introduced by the cross-interaction terms $x_i{x}_j$ indicating that both populations are needed to allow growth. As with Levins model, the amount of available area is limited to a maximum $1-D$. This is a high-dimensional system with many variables, but relevant results can be obtained if we consider a homogeneous parameter set. If $c_{ij}=c$, $e_i=e$ for all species, then we have:

$$\frac{d{x}_i}{dt}=c{x}_i\sum _{j=1}^S{x}_j\left(1-D-\sum _{w=1}^S{x}_k\right)-e_i{x}_i,$$

(16)

Now assuming that, since all species are parametrically identical, and indicating by $x={\sum }_i{x}_i$ the total population, we expect $x_k\sim x/S$. We can now write:

$$\frac{d{x}_i}{dt}=\frac{1}{S}\frac{dx}{dt}=c\left(\frac{x}{S}\right)x(1-D-x)-e\frac{x}{S},$$

(17)

Finally, the mathematical model describing the metapopulation occupancy can be written as follows:

$$\frac{dx}{dt}=c{x}^2(1-D-x)-ex,$$

(18)

If $e=D=0$, this model reduces to a simple logistic growth model where the population would eventually occupy the whole area.

The model is almost the same as Levins model with habitat loss, with just a nonlinear term $x^2$. But this linear term, indicating that interactions among individuals are taken into account, changes everything. Using the steady state condition:

$${\left(\frac{dx}{dt}\right)}_{\{c,D,e\}}=c{x}^2(1-D-x)-ex=0$$

(19)

for a given parameter combination $\{c,D,e\}$, three possible fixed points are found, namely $x_0^{*}=0$ and the pair

$$x_{\pm }^{*}=\frac{1}{2}\left(1-D\pm \sqrt{{(1-D)}^2-4r}\right),$$

where $r\equiv e/c$. These possible states are summarized in the bifurcation diagram shown in Figure 3a. The upper branch $x_{+}^{*}$ is a stable state (while $x_{-}^{*}$ is unstable) and as we can see, increasing levels of destruction gradually reduce the size of the vegetation cover. However, (and here comes a big point) there is a critical threshold:

$$D_c=1-2\sqrt{r}.$$

associated with a catastrophic shift. For $D>D_c$ the two equilibria $x_{\pm }^{*}$ do not exist in the biologically-meaningful domain, and the system faces extinction in a discontinuous manner. Using a spatially extended version of this model, we can observe a replacement process between two community states, namely ${\mathcal{C}}_R\to {\mathcal{C}}_D.$ After a given transient, the slight increase in D has produced a dramatic outcome.

We can use again the potential function to visualize the system’s dynamics against changing levels of habitat loss. The potential function now reads:

$$\mathcal{V}\left(x\right)=-c{x}^2\left(-\frac{r}{2}+\frac{1-D}{3}x-\frac{1}{4}{x}^2\right)$$

(20)

and is displayed in Figure 3b for several values of habitat destruction D. For $D<D_c$ (solid lines) two wells are found, corresponding to the two alternative stable states ($x_0^{*}$ and $x_{+}^{*}$) resulting in bistability, that will be achieved depending on the initial conditions. As we can appreciate, once the $D_c$ state is crossed, a rapid roll down of the marble will take place towards the only available state, located at the bottom of the desert state valley. But in this case, by contrast with Levins model, where we could consider a recovery strategy based on restoring habitat quality and available patches, another important phenomenon is present: there is an irreversible change, and improvements in habitat availability will not allow the system for getting back to the vegetated state again. In order to do so, restoration can only occur if we add the adequate amount of individual plants in order to put urselves above the unstable point $x_{-}^{*}$. Otherwise, the dynamics drives us again back to the zero-vegetation state. The implications of this phenomenon are obvious: any restoration effort once we cross the threshold is likely to require great investments and active efforts, which can have staggering costs.

One could argue that a realistic approach to drylands needs to consider many other variables beyond vegetation cover. Such a multivariate formulation has been developed, based on the concept of multifunctionality (MF). As implicit in the term, MF encapsulates multiple ecosystem functions. Using remote sensing and field data (Berdugo et al., 2017) a picture of global drylands emerges that gives aupport to the alternative state scenario depicted by our simple model. However, there are several potential thresholds associated to drylands degradation. In this context, a global analysis comprising $5\times {10}^4$ data points spanning different scales (from individuals to ecosystems) revealed three well-defined shifts (Berdugo et al, 2020) resulting from increasing aridity levels. In Figure 4d we summarize these results and the corresponding interpretation in terms of potential functions. As aridity $\alpha$ is increased, the thresholds include:

1.

For $\alpha =0.54$, a decay in photosynthetic activity and nitrogen storage that result from the shrinkage of leaf area. The shift is observable in particular as a marked reduction in plant productivity.

2.

For $\alpha =0.7$, a deep weakening of plant-soil interactions occurs, with a soil disruption involving a shift from stress-tolerant communities (STC) towards infertile shrublands (US). These transition has been studied within the context of shrub encroachment in arid grasslands, and is known to involve a sudden transition (D’Odorico et al., 2012; 2013).

3.

As we approach $\alpha =0.8$ we get back to the green-desert transition with a systemic breakdown. The shift is marked by drastic reductions in plant cover as well as changes in leaf traits associated with stress avoidance.

Can our understanding of these tipping points and their underlying emergent patterns be exploited to prevent them from occurring? This is largely an open question, but theoretical models suggest that their nonlinear properties could be exploited to avoid undesirable shifts. Using tools from statistical physics, it has been shown that changes in the patterns of fluctuations or dispersal could change the nature of the transition from discontinuous to continuous (Vila Martin et al., 2015). Similarly, using synthetic biology it would be possible to modify the dynamics of the microbiome of soil crusts in drylands (Belnap and Lange, 2001; Pointing and Belnap, 2012), creating synthetic soil crusts that could push away catastrophic shifts (Vidiella et al., 2020).

Finally, we should mention that drylands models also predict the presence of ghosts associated to the very long transient times to collapse close to tipping points. By making small-scale, regular interventions, these ghost dynamics can be exploited (even beyond the critical point) to maintain the vegetated state (Vidiella et al., 2018). In all these cases, it is important to point out that interventions are aimed at enhancing bottom-up properties of communities, such as biodiversity, by contrast with standard, top-down engineering (Solé, 2015; Solé and Levin, 2022).

5. Percolation: What Makes A System

Ecology is one of the few sciences where the systems-level view has played a central role since the early beginnings of the discipline (Odum, 1983; Wulf, 2015). Several common threads, such as space, hierarchy, information or network structure have a long tradition. But one could ask a naive question: what defines a "system"? (Keller, 2005). It might seem a rather vague question, perhaps leading (as it occurs with many fundamental questions) to qualitative answers and little predictive power. But here too, there is a powerful concept that provides an answer with deep implications for biodiversity conservation. Before we present the theoretical framework, let us consider three different case studies spanning three very different time scales, as shown in Figure 5a–c.

The first example is a spatial set of patches obtained from the canopy tree plot from Barro Colorado Island (Welden et al., 1991; Hubbell et al., 1999). The map contains a $250\times 250$ meters domain where black squares (with a $5\times 5$ meters width) indicate locations of a gap in the canopy, indicating a height smaller than $h_c=15$ meters (Figure 5a). This pattern changes over time, given the highly dynamic nature of the gap formation process (Brokaw, 1985; Whitmore, 1990; Scnitzer and Carson, 2001; Wright, 2002), but its statistical properties remain the same. As we can appreciate, clusters of canopy points are highly heterogeneous. In fact, if we look at the frequency of clusters of size s, it is found that it follows a power law distribution (Solé and Manrubia, 1995), namely:

$$P\left(s\right)=\frac{1}{Z}{s}^{-\gamma }$$

(21)

where $\gamma \sim 1.8$ is the so called scaling exponent and

$$Z={\int }_{s_0}^{s_m}{s}^{-\gamma }ds$$

(22)

a normalization constant, defined from the minimum gap size $s_0$ to the maximum observed ($s_m$). Additionally, the spatial organization of this area has been shown to display a fractal structure: if we zoom over this plot, the same pattern looks the same on different scales1. How do these scale-free structures emerge? One potential explanation of the patterns behind the BCI data is a critical state. An example would be the critical habitat loss $D_c$ depicted in Levins model. As we already discussed abpve, at the critical point several unexpected phenomena occur and available evidence strongly suggests that the rainforest might live close to criticality (see Schlicht and Iwasa, 2004).

The second example, which works on a spatial scale of Kilometres, is also tied to the problem of habitat connectivity in a heterogeneous landscape. The plot in Figure 5b shows a number of suitable habitat areas where a given species can thrive. If individuals of a given species are able to move from one patch to another, we say that these two sites (open circles in our plot) are connected. Formally speaking, two sites $S_i$ and $S_j$ will be connected if their Euclidean distance $d_{ij}$ is lower than the characteristic dispersal distance $d_c$ of the species under consideration (here the Northern Spotted Owl). This threshold condition is all we need to build the landscape connectivity graph $\mathcal{G}=(\mathbf{S},\mathbf{A})$, defined as the set of nodes $\mathbf{S}$ and a $n\times n$ matrix $\mathbf{A}$ (the adjacency matrix) where $a_{ij}=1$ if the pair $\{S_i,S_j\}$ is connected through dispersal. The key problem addressed by this graph is to measure the total number of connected sites. The graph shown here is a special one, obtained when the critical distance is $d_c=1500$. Below this value, the graph is disconnected into many separated sites. For $d>d_c$ becomes hyperconnected. Right on $d=d_c$, a sudden jump in the number of connected sites occurs. This threshold provides a powerful tool to predict the conditions for species persistence (Gardner et al., 1987; Urban and Keith, 2001; Keith et al., 1997).

Our third case study involves a global survey of 130 million forest fragments from three continents (Taubert et al., 2018). This study revealed two important things. The first was the presence of a robust scaling law in the size distribution of forest fragments. Specifically, it was found that it scales as described by Equation (20), now with a scaling exponent $\gamma \approx 2$ (see Figure 5c for the Americas set).

The three examples described above share the a threshold condition that is associated to the problem of percolation (Stauffer and Aharony, 2018; Christensen and Moloney, 2005). Percolation is not only the simplest example of a phase transition. It perfectly illustrates the idea that universal patterns emerge close to critical points and that they have crucial implications to fundamental features of very diverse types of systems. In fact, percolation defines a threshold condition that allows us to sharply separate a connected, coherent system from a set of many disconnected parts. In other words, percolation is at the core of what defines a complex system. The idea itself is very simple, and is fully illustrated by a two-dimensional toy model.

Consider a lattice $\Omega$ initially made of $L\times L$ empty sites. Let us now introduce "trees" by choosing random sites. A given empty site will be occupied by a tree with probability p. We want to define a measure of connectivity that tells us something about how easy is to move from one tree to its neighbours. As before, we will consider the four nearest sites as the neighbourhood of a given site. Using this definition, we can count the size of connected clusters. For small p values, we cannot expect to observe large patches of connected trees, while for large p we should be able to move across the whole landscape jumping from tree to tree. What about intermediate values? To address this question, let us define $P(p,\infty )$, i.e., the fraction of sites that belong to the largest connected cluster. This quantity tells us how much of the trees in our lattice belong to a single system, and compare it with the size of the whole lattice. When we plot this fraction against the occupation probability p, a rather unexpected result emerges, as shown in Figure 6a. Instead of a continuous increase, two well-defined phases are observed. The two phases are separated through a well-defined critical point $p_c\sim 0.59$. For very large L, it can be shown that $P(p,\infty )$ behaves as

$$P(p,\infty )\sim {(p-p_c)}^{\beta }$$

(23)

for $p>p_c$ (and zero for $p<p_c$) with $\beta \approx 1.16$.

If we measure the number of clusters of different size s, a power law distribution is obtained (Figure 6b) and the structure of the largest connected cluster is very complex, defining a fractal object (Figure 6c). The distribution is very similar to the one reported in (Taubert et al., 2018). This likely indicates that the observed forest fragmentation pattern is close to criticality. As deforestation proceeds, we should expect the system to cross the critical point, thus moving from an essentially connected system to a disconnected one (Figure 6d,e). What is the importance of these findings? Since deforestation rates and future estimates can be determined (within some bounds), the model provides an estimated time window to reach the transition to $p_c$.

Using different versions of the percolation model, the authors forecast for the America’s (where deforestation and reforestation rates are estimated) these time window would be around 2050. For each continent, estimated of expected deforestation and planned or required reforestation allow to exploit the robust theoretical predictions while guiding mitigation plans. In this case, by contrast with the transition to extinction in levin’s or facilitation models) there is a reversibility intrinsic to the nature of the percolation threshold. If reforestation is performed close to $p_c$ and within a time scale where fragmentation might not cause major population effects, recovery of connected landscapes can be achieved.

6. Fires and Ecological Time Asymmetries

Our last example deals with the long-time recognised existence of a separation of time scales in ecological dynamics. It has been known since a long time ago that slow-dynamical regimes are sometimes followed by sudden changes occurring in a very short time scale (Levin, 2000). These have been traditionally connected to the presence of multiple states separated by breakpoints (May, 1977). Our case study here is fire, which is known to be a very important environmental factor in many biomes. It was so long before the ascent of our species, but has become heavily influenced by human activity through the Anthropocene (Bowman et al., 2011; Hantson et al., 2015).

Fire is one particularly relevant example of this separation between slow and fast time scales that departs from our previous examples. It is connected with a phenomenon that, despite its global impacts, is seldom found in textbooks of conservation ecology (Bond and Keeley, 2005). Along with soil properties and environmental factors, fire has played a major role in shaping ecological communities (Bond et al., 2005), by acting as a "global herbivore". Fires fully illustrate the importance and implications of the asymmetry that pervades them. As the Catalan ecologist Ramon Margalef once said: "A forest needs one hundred years to get built, but burns in one day". But this only reveals one aspect of fires. They are self-organized, highly dynamical systems, coevolving with diverse range of species that became adapted to fire regimes.

What makes forest fires different? The answer needs to be found in the way they occur and how their emergence interacts with vegetation. Because they need vegetation cover to spread, a process of forest growth is first required, which takes a long time, while the effects are -in comparison- almost instantaneous. Since fires can get started anytime at some location and spread until they run out of fuel, what kind of dynamical state results from this dynamics? One important clue to answer this question is provided by fire size distributions. If s indicates the size of a given fire event, it has been shown (Pueyo, 2007) that the frequency of fires of size s follows a power law distribution (Figure 7a), which we have already described within the context of forest fragments i.e.,

$$P\left(s\right)=\frac{1}{Z}{s}^{-\beta }$$

(24)

where the normalization constant reads

$$Z=\left(\frac{1}{\beta -1}\right)s_0^{-\beta +1}$$

(25)

Here $s_0$ is the minimum forest fire size (usually in units of area) and we assume that largest fire sizes $s_m$ are such that $s_m\gg s_0$. As it occurs with other systems exhibiting power laws, the distribution has no characteristic scale (there are no "typical" fires with some expected, average size). While fires will be typically small, very large ones with scales close to the system’s size, are expected to occur. Since we have seen in previous examples that power laws are expected close to critical points, one could ask: what is the parameter that we need to tune in this case? The answer is that there is no tuning: the system self-organizes itself to criticality (Muñoz, 2018). The concept of selforganized criticality (SOC) was first proposed by Danish physicist Per Bak (Bak et al., 1987; Bak and Chen, 1991; Bak, 2013; Jensen, 1998) and is a keystone in our understanding of criticality in ecology and evolution (Solé et al., 1999; Levin, 1999). The intuitive example is provided by a sandpile: if we slowly add grains of sand to a surface (one by one), a sandpile if slowly formed with an increasing slope. At the beginning, the grains fall down and their behaviours is capture by standard physics. However, once the pile reaches a given critical slope, a collective process starts to unfold: one grain can trigger an avalanche. Most of the time, these avalanches are small, involving just a few grains, but very large (and fast) avalanches will also occur. The self-similar pattern of avalanches is an emergent phenomenon resulting from the constant driving followed by reset events.

A simple cellular automaton model (Chopard and Droz, 1998) of fire propagation in a virtual forest illustrates the idea within our context. Consider the so called Forest Fire Model (Bak et al., 1990, Drossel and Schawbl, 1992a, 1992b; Clar et al., 1996) where a discrete description of fire spreading is applied to the states of each site on a two-dimensional lattice $\Omega$. Each site $\mathbf{r}\in \Omega$ can be in three possible states, namely $S\left(\mathbf{r}\right)\in \Sigma =\{E,T,F\}$, where E denotes an empty cell, T a tree cell and F a fire cell. The state of the system will be updated by means of three probabilistic events, namely: (1) spontaneous burning of a tree, i.e., a transition $T\to F$ at a rate f, leading to a burning site (fire cell); (2) growth of new trees from empty sites, i.e., with a probability p we have $E\to T.$ (3) The last rule allows fire propagation: if a given tree has a burning neighbour, it burns too. This means a (deterministic) transition $T\to F$.

Despite its simplicity, the FFM exhibits a very complex behaviour. In particular, for $p,f\ll 1$ and $f\ll p$ the dynamics of fire propagation exhibits a self-organized critical (SOC) state that includes extreme events. An example of the time series generated by the model is displayed in Figure 7b,c. While fires have a marked bursting dynamics, tree cover reveals a continuous growth pattern followed by a reset that correlates with fire peaks. One snapshot of the model is also shown in Figure 7d, associated with a peak in the fire spreading. Trees, fires and empty sites are indicated as green, yellow and black colours, respectively.

Despite of its simplicity, versions of the FFM and variations of it have been successfully adapted to model the statistical patterns of the actual fires (Pueyo, 2007). In this context, it has been shown that the observed exponents in woodland forests depart from those predicted by the FFM. In particular, the sudden appearance of trees in random locations is clearly unrealistic. In fact, it has been shown that the accumulation of fuel after previous fires plays a key role. This qualitatively relevant observation can be easily incorporated (Pueyo, 2007) by a simple rule that accounts for fuel increase. This is done by introducing a propagation rate that would suddenly decay after fire and then slowly increasing. The resulting model (which also takes into account environmental factors) successfully accounts for the observed statistical patterns. Importantly, in this case the observed exponent can change in different geographic locations due to the influence of weather (Pueyo, 2007; Hantson et al., 2014).

The work by these authors reveals that forest fire dynamics combines features of the SOC mechanism (the separation of time scales) along with heterogeneities associated to fuel accumulation, which plays a role of an actual tunable parameter. By looking all these factors, a global map of fire size distribution can be obtained, with different areas marked by a characteristic, estimated exponent $\beta$ (Hantson et al, 2014). Since this parameter is affected by climate change, it offers a unique possibility predict the future impacts of global warming on fire sizes.

7. Discussion

The preservation and restoration of biodiversity stand as formidable challenges in an increasingly interconnected world. As we face the looming threat of crossing planetary boundaries in the coming decades, it becomes imperative to address these challenges with a comprehensive and multiscale approach. Biodiversity conservation necessitates considering the intricate interplay of multiple scales of complexity, and embrace theoretical models as a powerful approximation to reality. Such reality often departs from the linear intuition, and makes difficult to create a narrative that fits our intuitions. Since the success of conservation policies is not only a matter of good science science but also public understanding, this is not a minor problem.

Delving into the emergent properties arising from ecosystem dynamics offers valuable insights into developing successful conservation strategies. By identifying these emergent phenomena, which differentiate the various scales involved, we can design targeted interventions that tackle conservation challenges more efficiently. Moreover, uncovering the universal properties that underpin ecosystem dynamics allows us to develop simple yet powerful models for each scale, aiding in the formulation of evidence-based policies and management decisions. In this paper I have summarised a few well-established models that share a common simplicity in their formulation. This simplicity, however, does not mean that they fall short in explaining observed patterns. Instead, simple models of nonlinear systems can fully capture their crucial properties.

A summary of the models presented here, the underlying case studies and their relationships is provided in Figure 8. From bottom ti top, the different kinds of predicted transitions described by the previous models and the mathematical models are displayed. This scheme also helps emphasising the idea that the introduction of extra rules can generate novel, emergent phenomena. When dealing with single-species studies (a,b), species are considered as isolated entities whose populations change as a consequence of a balance between dispersal (immigration) and extinction. In these models, intra-specific competition is the dominant interaction, although the estimated rates are of course somehow coarse-graining interspecific traits. If habitat loss is introduced (c,d), an extra parameter D modifies the previous transitions, but have a common continuous behaviour until extinction occurs at a given critical habitat destruction level.

If an extra layer of complexity associated with species interactions (beyond competition) is introduced, the nature of the transition dramatically changes. Here we have considered facilitation (e) where local interactions have a cooperative nature that translates into a quadratic growth term. This means that increasing changes in the key parameter (habitat destruction) trigger a catastrophic shift after the threshold level is reached. Finally, on a very large scale (f) the process of habitat loss requires a model that takes space into account in an explicit manner. In this case, percolation is the conceptual framework that provides robust predictors for fragmentation thresholds. Each of these models is able to guide conservation efforts at a given scale, but cannot be extrapolated as we move to the next level. Each level of complexity is irreducible. The selection of the models and associated case studies showcased here is based on their deliberate simplicity. This choice serves as a demonstration of the efficacy of well-defined models. Due to their straightforward nature, these models enable us to methodically analyze the significance of various variables and ascertain their relative importance, or lack thereof. These models can also be expanded upon to encompass additional aspects of population dynamics and their interaction with the environment, as demonstrated by Bauch et al. (2016), Weinberger et al. (2017), and Farahbakhsh et al. (2022).

A notable example pertains to the issue of habitat fragmentation arising as a consequence of agricultural expansion. In the work by Bengochea Paz et al. (2022), the authors present a model that intertwines agricultural expansion with human population growth. The latter is dependent on the former, leading to increased habitat loss and fragmentation. Consequently, this results in the depletion of ecosystem services provided by these fragmented habitat patches and a decrease in overall productivity. This, in turn, can drive land managers to further expand crop areas, exacerbating habitat loss. This feedback loop eventually pushes the system towards a state of percolation. The model’s predictions indicate a transition from sustainable states to a phase of collapse. The precise boundary between these two phases hinges on the delicate balance between land use under intensive agriculture and the rate at which demographic changes occur.

Embracing a complex systems view empowers us to approach biodiversity preservation with greater foresight, acknowledging the complex web of dependencies and feedback loops that sustain life on Earth. In such a world, where the fate of our planet hangs in the balance, it is vital to integrate sound scientific knowledge, forward-thinking management strategies, and global cooperation to safeguard and restore biodiversity. Such integration will require a collective effort across disciplines where our understanding of different scales will benefit from the universal patterns displayed by simple models of ecological complexity.