# Phase-separation physics underlies new theory for the resilience of patchy ecosystems

Edited by Pablo Marquet, Pontificia Universidad Catolica de Chile, Santiago, Chile; received February 19, 2022; accepted November 22, 2022

Commentary

January 25, 2023

## Significance

Human-induced environmental changes push ecosystems worldwide toward their limits. Therefore, there is a growing need for indicators to assess the resilience of ecosystems against external changes and disturbances. We highlight a novel class of spatial patterns in ecosystems for which resilience indicators are lacking and introduce a new indicator framework for these ecosystems, akin to the physics of phase separation. Our work suggests that aerial imagery can be used to monitor patchy ecosystems and highlights a link between physics and ecosystem resilience.

## Abstract

Spatial self-organization of ecosystems into large-scale (from micron to meters) patterns is an important phenomenon in ecology, enabling organisms to cope with harsh environmental conditions and buffering ecosystem degradation. Scale-dependent feedbacks provide the predominant conceptual framework for self-organized spatial patterns, explaining regular patterns observed in, e.g., arid ecosystems or mussel beds. Here, we highlight an alternative mechanism for self-organized patterns, based on the aggregation of a biotic or abiotic species, such as herbivores, sediment, or nutrients. Using a generalized mathematical model, we demonstrate that ecosystems with aggregation-driven patterns have fundamentally different dynamics and resilience properties than ecosystems with patterns that formed through scale-dependent feedbacks. Building on the physics theory for phase-separation dynamics, we show that patchy ecosystems with aggregation patterns are more vulnerable than systems with patterns formed through scale-dependent feedbacks, especially at small spatial scales. This is because local disturbances can trigger large-scale redistribution of resources, amplifying local degradation. Finally, we show that insights from physics, by providing mechanistic understanding of the initiation of aggregation patterns and their tendency to coarsen, provide a new indicator framework to signal proximity to ecological tipping points and subsequent ecosystem degradation for this class of patchy ecosystems.

### Sign up for PNAS alerts.

Get alerts for new articles, or get an alert when an article is cited.

Theoretical models and observations (1–8) suggest that the large-scale regular spatial patterns found in many ecosystems change in a predictable way when approaching a potential tipping point, which is considered an important indicator for assessing ecosystem resilience (9–12). This indicator framework is based on reaction-diffusion models that generate spatially periodic (regular) patterns through a so-called Turing instability (13–15). These self-organized patterns form due to a combination of short-range facilitation and long-range inhibition, generally referred to as scale-dependent feedbacks (16–18). Such feedbacks result, for instance, from the concentration of a limiting resource (e.g., water) that is flowing through the system (12, 13, 19, 20), or the production of a growth-inhibiting substance (14, 21, 22). However, many ecosystems exhibit patterns that are far less regular and lack clear scale-dependent feedbacks (Fig. 1). These include non-regular patterns generated by self-propelling disturbances (23), or patterns of coexisting patches resulting from positive feedback processes (11). Non-regular patterns have for instance been found of in grazing systems (24–28), biogeomorphological systems (29–33), such as seagrass beds and dune fields, and nutrient-poor systems (34–39), such as peatlands, meaning that the leading indicator framework does not apply to an important group of patchy ecosystems.

Fig. 1.

Here we identify a new generic mechanism behind self-organization of ecosystems, involving the aggregation of an abiotic or biotic species into non-regular patches. We will refer to this new mechanism as density-dependent aggregation, and we propose it as a possible explanation for the ubiquitous irregular patterns observed in ecosystems that apparently lack scale-dependent feedbacks. This mechanism involves aggregation of a mobile species or substance that is approximately conserved within the ecosystem, meaning that its fluxes in or out of the ecosystem are negligible compared to local redistribution processes. We demonstrate that these ecosystems can be expected to exhibit fundamentally different dynamics compared to ecosystems governed by scale-dependent feedbacks and introduce a novel indicator framework for ecosystem resilience that is based on phase-separation theory (40–44). First, we highlight a number of ecosystems where density-dependent aggregation may be responsible for the formation of self-organized spatial patterns. Second, we present a generalized, universal model of these ecosystems, discuss the feedbacks that trigger spatial pattern formation in this class of ecosystems, and demonstrate that they are expected to exhibit principal dynamics that are akin to phase separation, a concept stemming from physics and increasingly highlighted for biological systems (45–52). Finally, building on insights from physics, we introduce a new indicator framework to signal if ecosystems are approaching tipping points at which dramatic shifts in ecosystem states occur.

## Regular Patterns and Scale-Dependent Feedbacks

Arid vegetation patterns are a well-known example of spatial self-organization by means of scale-dependent feedbacks. In arid ecosystems, the growth of plants is locally enhanced through the concentration of water by plants (called small-scale positive feedback). This concentration process leads to a deprivation of water on a larger spatial scale (large-scale negative feedback; Fig. 2

*A*), leading to dense patches of bush interspersed with bare soil (13, 53–55). Similarly, regular spatial patterns in intertidal mussel beds have been explained by scale-dependent feedbacks resulting from blue mussels that decrease their own mortality locally, while competing for algae on larger spatial scales (17, 56).Fig. 2.

There are many other examples of spatial, scale-dependent feedbacks generating regular patterns in ecology (16, 18, 57, 58). A common characteristic of these patterns is that positive feedback is dominant at short spatial scales, while negative feedback prevails at greater distances. A crucial condition for scale-dependent feedbacks to occur, and hence for these patterns to have a characteristic scale, is that resources are locally replenished and depleted (e.g., acting as a dissipative system). For example, water enters arid ecosystems through rainfall and leaves through evaporation (5, 18), and carbon enters and leaves mussel beds through algal inflow and respiration by mussels (17, 59), respectively. There are, however, ecosystems where spatial self-organization is caused by feedbacks between a sessile engineering species and a mobile element (for instance a nutrient, or grazer, see examples in Table 1) that are not locally replenished and depleted, but instead are largely conserved within the system over relevant timescales (weeks, months, years). In these ecosystems, patchiness is driven by a conserved ecosystem component, either a resource or a consumer, whose mobility is determined by a sessile species. This class of ecosystems includes irregularly patterned grazing systems, biogeomorphic ecosystems, and nutrient-poor ecosystems, as we have summarized in Table 1. We highlight these classes in the following sections.

Table 1.

Ecosystem | Variables | Fluxes | Feedbacks |
---|---|---|---|

Grazing systems | $V$: Foraging herbivores $U$: Grazing herbivores | A: Herbivores switch to foragingB: Herbivores switch to grazing | 1: Grazing herbivores reduce plant biomass. Less attractive high biomass causes herbivores to stop grazing and start foraging. 2: Grazing herbivores reduce plant biomass. More attractive low biomass causes herbivores to stop foraging and start grazing. |

Biogeomorphic systems | $V$: Aeolian/suspended sand $U$: Deposited sand | A: ErosionB: Deposition | 1: Elevated areas enhance productivity. Plants reduce erosion. 2: Elevated areas enhance productivity. Plants enhance deposition. |

Oligotrophic ecosystems | $V$: Dissolved resource $U$: Sessile species | A: Release of resource through mortalityB: Uptake, trapping, and demobilization of resource | 1: $U\stackrel{-}{\to}A$ Reduced release of resource with increasing biomass of sessile species. 2: $U\stackrel{+}{\to}B$ Enhanced uptake, trapping, and demobilization of resource with increasing biomass of sessile species. |

Variables, fluxes, and feedbacks refer to Fig. 2

*B*. Variables $U$ and $V$ represent the less mobile and highly mobile state of the conserved species. $A$ and $B$ are the fluxes between $U$ and $V$. Feedbacks 1 and/or 2 depicted in Fig. 2*B*are established through the interaction between $U$ and a sessile species, as depicted in the right-hand column. In nutrient-poor ecosystems, $U$represents the ecosystem engineer, and the feedbacks result from facilitation (decreased mortality). Examples of the listed feedbacks are discussed in the main text.## Ecosystems with Aggregation-driven Patchiness

### Grazing Systems: Two-Phase Mosaics in Terrestrial Rangelands and Intertidal Seagrass Meadows.

Many natural and domestic grazing systems are characterized by strong patchiness in the form of alternating short-grazed low biomass areas and high biomass patches that are hardly grazed. In these so-called two-phase mosaics, herbivores typically prefer grazing in the low biomass areas, thereby maintaining this state of low biomass in these areas [sometimes referred to as grazing lawns. (24, 60)]. There are various reasons for this preference. In shrublands, plants become woody as they mature or develop thorns or spines to guard against grazing. In grasslands, grasses become less nutritious with their increase in biomass (Fig. 1

*A*). In intertidal seagrass meadows, sediment accumulates in dense seagrass patches, which stabilize the sediment and prevent erosion, hampering waterfowl in applying their preferred grazing strategy (dabbling) during low tide, as the dense patches are then elevated above the water, leaving only the water-logged hollows for grazing (see refs. 24, 25, and 61 and Fig. 1*B*). In these grazing systems, herbivores are a mobile agent that distributes depending on the spatial distribution of the plants, the sessile species. Plant density is in turn controlled by the grazing pressure of the herbivores. Although the local plant biomass affects the movement behavior of herbivores, it does not affect their reproduction and mortality rates, at least not at the same timescales. This means that, unlike arid ecosystems or mussel beds, the production and depletion of a resource, and the resulting scale-dependent feedbacks, do not play a role. Instead, the aggregation of herbivores in the short-grazed patches is what drives pattern formation.### Biogeomorphic Systems: Coastal Dunes and Nebkha Landscapes.

In desert and beach systems, characterized by mobile sand, vegetation can be found to create hummocked structures. Examples include relatively small nebkha dunes, observed in many desert ecosystems (Fig. 1

*C*), and extensive beach dunes covered by marram grass (Fig. 1*D*). These landscapes are formed through feedbacks between vegetation and sand. Vegetation decreases the mobility of sand by enhancing deposition rates and/or decreasing erosion rates, thereby demobilizing aeolian sand and preventing (re)mobilization of deposited sand, respectively. Sand, in turn, enhances the growth of plants through various mechanisms. It provides a fresh pathogen-free substrate for many dune grasses (31, 62, 63), enables vegetation to escape erosive conditions and elevated salinity levels, and provides nutrients and acts as water storage in desert systems (33, 64). Although vegetation traps sand in these ecosystems, it does not consume sand or remove it in another way, as is the case with water in arid ecosystems and algae in mussel bed ecosystems. Instead, plants in dune landscapes control the mobility of sand. Hence, the aggregation of sand in plant patches is the mechanism creating patterned dune landscapes, rather than scale-dependent feedbacks.### Nutrient-Poor Ecosystems: Coral Reefs, Peat Lands, and Subtidal Seagrass Meadows.

Many nutrient-poor ecosystems exhibit patchiness, as the limited amount of resources within these systems does not allow for a uniform cover to be attained (Fig. 1

*E*and*F*). In these systems, resources are typically either suspended in (soil) water or have been taken up by organisms. Uptake by sessile organisms demobilizes the resource, while resources are again released to the water during or after an organism’s life. In coral reefs for example, filter-feeding corals and sponges passively or actively filter particulate organic matter out of the water column, an important source for inorganic nutrients. The nutrients are again released to the water due to shedding of cells by sponges or the secretion of mucus by corals (the latter is converted back to particulate form by sponges, also known as the “sponge loop”, refs. 65–69). This feedback loop, in combination with the mobility of suspended particles, is likely to amplify patchiness in many coral reefs. Hence, facilitative interactions play a key role in pattern formation in coral reef ecosystems, as it reduces the mortality-induced loss of resources. Similarly, patterning of subtidal seagrass meadows has been associated with facilitative effects, including stabilization of fine-grained sediment by the dissipation of wave energy (70). Here, seagrasses grow better when there is more sediment, while entrapment of sediment increases with seagrass density creating a local positive feedback loop. Besides demobilization of resources through uptake, resources can be demobilized in other ways. For example, phosphorous-limited Siberian peatlands exhibit patterning (Fig. 1*E*) which is linked to the vascular plants that draw down the water table through evapotranspiration during the growth season, preventing phosphorous to leave the ridges dominated by dwarf shrubs, beech, and pine trees (34, 71, 72). Hence, phase separation could provide a viable alternative explanation for pattern formation in peat lands.## Results

### A Generic Conceptual Framework.

The ecosystems described above have in common that their patchiness is caused by an ecosystem component, either a resource or a consumer, whose mobility is determined by a sessile species. In these ecosystems, the dynamics of the mobile species can be described as being either in a highly mobile state $V$ (e.g., aeolian sand, foraging herbivores, and dissolved resources) or in a less mobile state $U$ such as deposited sand, grazing herbivores, and absorbed resources, with the sessile species controlling the fluxes between these two states (Fig. 2

*B*). The sessile species can be modeled explicitly as*U*in case of a plant consuming a mobile nutrient, or implicitly in case of an arid shrub immobilizing sand, transferring from loose aeolian sand*V*to immobilized, fixed sand*U*under its roots.The reversibility of the resource and consumer states is a crucial element of the model framework. This can be captured by a system of reaction-diffusion equations:

$$\frac{\partial U}{\partial T}=F\left(U,,,V\right)+{D}_{U}\left(\frac{{\partial}^{2}U}{\partial {X}^{2}},+,\frac{{\partial}^{2}U}{\partial {Y}^{2}}\right),$$

[1a]

$$\frac{\partial V}{\partial T}=G\left(U,,,V\right)+{D}_{V}\left(\frac{{\partial}^{2}V}{\partial {X}^{2}},+,\frac{{\partial}^{2}V}{\partial {Y}^{2}}\right),$$

[1b]

where $U(X,T)$ and $V(X,T)$ represent the mobile species in its less mobile and its highly mobile form, respectively, $T$ represents time, and $X$ and $Y$ represent space. These equations describe the conversion of mobile species $U$ into sedentary species $V$ and vice versa through reaction terms $F(U,V)$ and $G(U,V)$, as well as the lateral movement through diffusion. Since the mobile species (herbivores, sand, or resources) is not locally produced, nor destroyed and because it is assumed to reside either in the form of $U$ or $V$, we write

$$G\left(U,V\right)=-CF\left(U,V\right),$$

[2]

where $C$ is a positive constant converting $U$ into $V$, in real-world terms for instance the incorporation of nutrients into biomass. Substituting and scaling gives (see

*SI Appendix*, Texts S1, S2, and S5 for details)$$\frac{\partial u}{\partial t}=f\left(u,v\right)+\delta \left(\frac{{\partial}^{2}u}{\partial {x}^{2}},+,\frac{{\partial}^{2}u}{\partial {y}^{2}}\right),$$

[3a]

$$\frac{\partial v}{\partial t}=-f(u,v)+\left(\frac{{\partial}^{2}v}{\partial {x}^{2}},+,\frac{{\partial}^{2}v}{\partial {y}^{2}}\right),$$

[3b]

with $\delta $ representing the ratio between the rate of diffusion of $V$ and $U$, and hence, δ < 1 is required for spatial patterns to develop. In this paper, we adopt $f(u,v)$ equals $uv-\frac{\alpha u}{u+\beta}$ as a dimensionless generic model (see

*Materials and Methods*). As Eq.**3**might prove abstract and unspecific, we provide an example model of a nutrient-poor ecosystem in Table 2. It is worth noting that Eq.**3**does not follow Turing’s mechanism for pattern formation because its dynamics effectively follow a spinodal phase separation (73, 74), as we discuss later.Table 2.

**3**

Mathematical equation | Parameter |
---|---|

$\begin{array}{c}\frac{\partial P}{\partial t}\phantom{\rule{-1pt}{0ex}}=\phantom{\rule{-1pt}{0ex}}{c}_{1}NP-\frac{{c}_{2}}{{c}_{3}+P}P+{D}_{P}\left(\frac{{\partial}^{2}P}{{\partial x}^{2}},+,\frac{{\partial}^{2}P}{{\partial y}^{2}}\right)\hfill \end{array}$ $\frac{\partial N}{\partial t}={c}_{4}\left(-,{c}_{1},N,P,+,\frac{{c}_{2}}{{c}_{3}+P},P\right)+{D}_{N}\left(\frac{{\partial}^{2}N}{{\partial x}^{2}},+,\frac{{\partial}^{2}N}{{\partial y}^{2}}\right)$ | c_{1} : Plant uptake coefficientc_{2}, c_{3} : Plant mortality constantsc_{4} : Nutrient to biomass ratioD_{i} : Diffusion constant of P or N |

The model represents a plant $P$ that consumes a very mobile nutrient $N$ that is conserved at large spatial scales in the ecosystem. Plant growth here is determined by its consumption of nutrient and described as ${c}_{1}NP$. Plant mortality decreases with plant biomass due to a self-facilitation process and is described by ${c}_{2}P/({c}_{3}+P)$. Nutrients are removed from the nutrient pool by plant uptake and is replenished by plant mortality. The ratio of nutrient to biomass in plant tissue is given by ${c}_{4}$.

The densities of $u$ and $v$ can vary in space and in time; however, the total amount of the mobile species within the modeled domain remains constant in the absence of external sources. The conserved, spatially averaged density of the mobile species ($u+v$) will be referred to as the global density $\langle \tau \rangle $ :

$$\langle \tau \rangle =\langle u+v\rangle =v+\delta u,$$

[4]

where $\langle \ast \rangle $ takes the spatial average (

*SI Appendix*).Reaction term $f(u,v)$ includes both fluxes from $v$ to $u$ (demobilization, for instance deposition) and fluxes from $u$ to $v$ (mobilization, for instance erosion), as depicted with arrows

*A*and*B,*respectively, in Fig. 2*B*. In*SI Appendix*, Figs. S1–S3, we derive that for pattern formation to occur, $u$ needs to reduce mobilization of $u$ into $v$ (dashed thin arrow 1 of Fig. 2*B*), and/or enhance demobilization of*v*into*u*(dashed thin arrow 2 of Fig. 2*B*), amplifying its own local density. Such feedbacks can emerge directly, for instance when a mobile nutrient is incorporated in plant tissue, or through interactions with a sessile engineering species with slower dynamics, for instance when aeolian sand is fixed by dune or desert vegetation, as discussed in the previous sections and given schematically in Table 1.When in a patchy state, the net diffusion of $u$ and $v$ occurs in opposite directions, as illustrated by the horizontal arrows of Fig. 2

*B*. For the equilibrium states of Eq.**3**, it is possible to derive the conditions of spatial mass redistribution as defined by the diffusion ratio (diffusion of $u$ to diffusion of $v$, as shown in,*SI Appendix*, Fig. S7), as a function of $u$. Fig. 2*D*shows that the net diffusion can be negative for sufficiently high local densities of $u$, leading to aggregation and subsequently an increase in $u$. Hence, in Eq.**3**, density-dependent aggregation drives pattern formation and dynamics by determining the spatial fluxes of $u$ and $v$ between neighboring low-density and high-density domains (74).### Phase Separation and Ostwald Ripening.

A number of fundamental differences exist between pattern formation arising from scale-dependent feedback (e.g., Eqs.

**8**and**9**in*Methods*) and density-dependent aggregation (Eq.**3**). When scale-dependent feedbacks describe biological or physical processes with clearly defined spatial scales, they generate a stable periodic pattern that peaks at regular intervals (Fig. 3*A*). Patterns caused by density-dependent aggregation are defined by two states, or “phases,” coexisting in space and separated by a transition or interface. This interface between the lower and higher plateau is mobile, and consequently patches in systems with density-dependent aggregation mechanism remain dynamic. As a result, density-dependent aggregation is characterized by patterns defined by a range of spatial scales (but not scale-free), following a log-normal distribution. Hence, density-dependent aggregation is characterized by spatial patterns consisting of coexisting phases (resembling alternative stable states), as shown in the model simulation of Fig. 3*C*and*D*. In physics, the emergence of such alternating phases out of uniformity is referred to as phase separation (75, 76). Phase separation finds its origin in the unmoving of two immiscible fluids like oil and water as put forward by Cahn and Hilliard (75, 76), but has also been identified as a mechanism behind cell polarization (51, 77–79) and has been suggested to be the principle mechanism behind the aggregation of animals (24, 42, 80–82) in recent years.Fig. 3.

Movie S1.

Movie S2.

Non-regular and scale-free patterns have been put forward by a wide range of modeling approaches (see refs. 6, 7, 27, 83, and 84 for examples). Yet, the density-dependent aggregation mechanism has further properties separating it from scale-dependent feedback, especially when scale-dependent feedbacks are limited to realistic spatial scales. The lack of a characteristic patch size leads to a pattern that is not only irregular [it may self-organize into a disordered hyperuniform state (85, 86)], but that also becomes coarser as time progresses, even though the fractional cover remains constant. This resembles a phenomenon known in physics as Ostwald ripening (43, 87, 88), which describes how small particles dissolve and contribute to the growth of larger particles. Similarly, in the density-dependent aggregation model, the mean patch size (radius) becomes larger over time, while the patch-size distribution remains the same. In fact, patch-size distributions at different moments in time collapse to the same distribution when the patch radius is scaled by dividing by the mean patch size, which is in line with Ostwald ripening theory

*.*Phase-separation theory furthermore puts forward that the average wavelength of the patterns is expected to increase following a coarsening power law, as is predicted by the Cahn–Hilliard model (89), and observed in ecosystems (24), and sorted patterned ground (90). To the contrary, the wavelength of the patterns predicted by scale-dependent feedback models is typically centered around a constant value, although a multitude of wavelengths can develop under different conditions (11, 91, 92).### Resilience to Local Disturbances.

A crucial characteristic of the patchiness created by density-dependent aggregation is that a local disturbance beyond a critical threshold (open circle in Fig. 2

*D*, and dashed line Fig. 3*B*), can tip a patch from a state of expansion with persisting aggregation to a state of decline where the patch is being dispersed. As a consequence, this patch will disappear even when the disturbance has passed. Moreover, besides the bistability within patches, dynamics between patches can further amplify local disturbances. Fig. 4 shows that, in the absence of disturbances, larger patches tend to grow (green region of Fig. 4*A*) at the expense of smaller patches that tend to shrink (red region). If a growing patch is reduced in size due to some local disturbance, pushing it below a critical size threshold (the boundary of the red and green areas in Fig. 4*A*), it will start to shrink, eventually leading to the complete loss of the patch (Fig. 4*C*). Hence, the resilience of aggregation patterns to local disturbances is low due to a combination of bistability of patches and between-patch competition. Yet, the model highlights that due to a water bed-like effect, patches will expand on other locations as they absorb the mobilized substances (resources or species). Therefore, resilience to local disturbances is much higher on system-wide scale (cf. Fig. 4*C*and*D*).Fig. 4.

Similar to what is observed in systems with scale-dependent feedback (5, 11, 15, 93), gradual changes in environmental conditions can trigger rapid and hard to reverse shifts in ecosystems governed by density-dependent aggregation. Such gradual changes can for example be caused by the inflow or outflow of resources from elsewhere or the migration or hunting of herbivores, expressed in our model as changes in the global density of the aggregating species or substance ($\langle \tau \rangle \equiv \langle u+v\rangle $). Both global density

*τ*and diffusive ratio δ are essential degrees of freedom because they are the only control variables that dynamically affect the emergent properties of ecosystems (Fig. 5*A*). We present the bifurcation diagram of Eq.**3**showing how the phase-separation patterns and their stability properties are affected by the global density $\tau $ by means of a combination of theoretical analysis (*SI Appendix*, Texts S4 and S6) and numerical simulation, as shown in Movies S3.Fig. 5.

Movie S3.

Fig. 5

*A*shows that, depending on the global density $\langle \tau \rangle $, the model resides in either a patchy state or in a uniform state. The region in parameter space in which only a patchy state is stable is referred to as the spinodal region in physics. Here patches form out of small inhomogeneities in a uniform ecosystem state, a process known as spinodal decomposition (76), similar to diffusive instability in the context of reaction-diffusion models (15). Gradual changes, corresponding to pattern formation at the border of the spinodal region, lead to large-scale patterns, whereas more rapid changes, that allow patterns to form in the center of the spinodal region, lead to smaller-scale patterns. The spinodal region shares its border with a region referred to as the binodal region where, instead of through infinitesimal inhomogeneities of a uniform state, patches form through nucleation. Here, pattern formation requires nucleation sites, that is, large preexisting local inhomogeneities. It can be shown for Eq.**3**that a binodal region always exists on both sides of the spinodal region (*SI Appendix*, Texts S4 and S6), meaning a transition to a patchy state can be triggered by the introduction of nucleation sites (for instance, local restoration measures) under conditions where patches cannot establish by themselves. Due to the existence of both a patchy and a uniform ecosystem state (bistability) in the binodal region, both pattern formation and the loss of patterns are characterized by abrupt shifts and hysteresis associated with ecosystem degradation. This behavior is similar to what is observed in systems with scale-dependent feedback.Box 1. A Physics-Derived Indicator Framework

## Discussion

Self-organized patterns have been proposed as a key indicator for imminent catastrophic shifts in ecosystems, as well as an adaptation of ecosystems evading dramatic changes, acting as a hallmark of ecological resilience (1, 8, 11, 27, 83, 97, 98). Thus far, this concept has mainly been applied to systems with regular patterning, where pattern formation is driven by scale-dependent feedback processes that follow Turing’s activator–inhibitor principle (16, 53, 93). The specificity of this mechanism limits the wide application of spatial patterns as an ecological indicator system. Here, we identify a different process for self-organized patterning that is akin to the physical process of phase separation (73, 76, 89). This type of self-organization applies well to ecosystems where pattern formation is determined by more or less fixed amounts of biotic or abiotic agents such as nutrients, sediment, or consumers which become aggregated in patches, determining patch cover (53). Exchanges of these abiotic or biotic resources among the patches can result in apparent, asymmetric competition between neighboring patches (74). This results in a patch coarsening process, where larger patches absorb more and more resources, while smaller patches wither away (Figs. 3 and 4). Such coarsening dynamics are a defining characteristic for phase-separating dynamics in ecosystems (24, 42, 80, 90).

Density-dependent aggregation and the resulting dynamics have crucial implications for the resilience as well as the management of patterned ecosystems. While systems with regular patterns driven by scale-dependent feedback are highly resilient in returning to a patterned state while perturbed (10–12), perturbations have very different effects in ecosystems with phase-separation dynamics (cf. Fig. 4

*C*and*D*). Any localized perturbation is likely to amplify the competitive interactions between patches and may lead to local collapse of patches in the landscape. Yet, the resources contained in this patch will redistribute over the landscape, increasing patch growth elsewhere, creating a water bed effect. For example, sediments lost from degrading seagrass patches may end up in other seagrass patches, facilitating growth in these locations. Hence, phase separation will strongly decrease local resilience, hampering any restoration efforts, but maintaining ecosystem resilience at larger spatial and temporal scales. If human restoration efforts focus on on-the-spot restoration, for instance in very small reserves, phase-separation dynamics might limit recovery of localized patches of a mobile species. If restoration efforts involve large, landscape scales, phase-separation dynamics are beneficial, as local collapse will stimulate recovery elsewhere in the landscape. Hence, populations or ecosystems governed by phase-separation dynamics need to be managed in a fundamentally different way, focusing on landscape-wide, ecosystem-based management, rather than on localized efforts to compensate for human-induced disturbances.Human-induced global changes push ecosystems worldwide toward and over their limits (99, 100). Therefore, there is a growing need for indicators to assess the resilience of ecosystems against external changes and disturbances. So far, existing studies have pointed out that self-organized patterns can serve as excellent hallmark of ecological resilience (1, 4, 6, 27, 101). However, these patterns are based on the well-defined Turing-like patterns. Hence, the patterns are stationary in time. In contrast, here we have introduced a physics-based indicator framework to determine the resilience of another important subset of self-organized patchy ecosystems, including grazing systems, biogeomorphic systems, and nutrient-poor ecosystems. Our new generic framework also displays spatial self-organized patterns, but the patterns display coarsening behavior where patch size increases over time, which results in a more patchy pattern with less regularity, comparable to so-called coexistence patterns (11, 91, 92) or disturbance-driven patterns (23). However, individual aggregation patches display either shrinkage or growth depending on their size relative to the critical patch size (Fig. 4), as well as a relatively constant patch cover, a dynamic not previously been found in self-organized ecosystems. We emphasize again that the basic mechanism for patterning presented here does not depend on the specific feedbacks chosen for $f(u,v)$ in Eq.

**3**(i.e., different from the well-studied scale-dependent feedback mechanism). They have the common property that they lead to aggregation of a conserved, mobile substance or species into a patchy, less mobile state.Our study here provides a generic theoretical framework for ecological systems that approximate mass conservation (42, 73). This generic theory provides an explanation for observations in some grazing systems that are not well understood. For example, Berg et al. (26) showed that for a constant herbivore density, the fractional cover of less digestible tall salt-marsh vegetation remains the same (47%), but that the pattern becomes coarser over time. Kéfi et al. (27) showed that patch-size distributions in Mediterranean rangelands start to deviate from a right-tailed distribution [power-law or log-normal distributions; (28, 102)] as grazing pressure increases, as predicted by our model. Similarly, Xu et al. (103) observed a decrease in the skewness of patch sizes in a global dataset of irregularly patterned arid rangelands as cover decreases, in line with our theoretical findings. Our theory generalizes these observations from grazing systems, meaning that similar trends are expected to be found in nutrient-poor systems and biogeomorphic systems with a limiting nutrient or substrate, respectively, that is contained in the landscape.

Scale-dependent feedback and density-dependent aggregation processes are not mutually exclusive in causing pattern formation in ecosystems, as both types of processes can interact to sculpt the spatial distribution of species (33, 104). For instance, in desert plants (bromeliad

*Tillandsia landbeckii*) surviving on advective fluxes of humidity (fog) form banded patterns through interception of fog-water, creating a scale-dependent feedback. In turn, vegetation aggregates aeolian sand to create nebkha-like landscapes, which follow the density-dependent aggregation principle (33, 104). Similarly, in mussel beds, large-scale patterns form through density-dependent feedback due to the interplay of facilitation and competition between mussels (10, 17, 105), while aggregative movement within mussels generates pattern formation at small spatial scales (<1 m) nested within the large-scale patterns (42, 56). Model analysis highlighted that the nested patterns have a strong positive effect on the resilience of the mussel beds. A third example is the formation of bacterial colony patterns in a particular strain of*Bacillus subtilis*(106, 107), where several mechanisms determining movement interact at the same time. Finally, physical and geological processes can impose patterns and other forms of heterogeneity upon ecosystems, for instance in the form of stone aggregation driven by frost upheaval in arctic ecosystems (90). Large-scale landscape heterogeneity, for instance in the slope along hillsides, can generate a multitude of pattern wavelengths to occur close to each other. Moreover, recently insights into multistability of patterns (3) in ecosystems with scale-dependent feedback highlight that a multitude of stable, alternate pattern wavelengths is possible even without landscape variability. Hence, the interplay of biological and physical/geological processes would create landscapes that exhibit a variety of patterns with different lengths and shapes. In many real-world ecosystems, multiple pattern-forming processes interact to determine ecosystem shape and functioning.Non-regular, patchy patterns have been observed in a wide range of systems. Many of them are caused by processes different from aggregation of substances or species and have been described by models that follow different principles (28, 108) than the one proposed in this paper. Small-scale, self-propelling disturbances can generate large-scale, scale-free patterns in forest prone to fire (23, 84), or mussel beds on rocky shores (109), a pattern-forming process that is referred to as self-organized criticality in the literature (109). Scale-free patterns also have been predicted to occur in relation to intense grazing (27), or fast water redistribution in arid ecosystems (7). So-called coexistence patterns may develop without mass conservation in systems where strong local positive feedbacks create two-phase patchiness resulting from bistability (11, 91). The current study provides a number of handles to distinguish the above patterns from patterns caused by density-dependent aggregation. First, patch-size distributions caused by density-dependent aggregation are not scale-free but follow a log-normal distribution. Second, the patterns caused by density-dependent aggregation are predicted to coarsen in time and hence do not form a stationary patch distribution. This feature separates them from scale-dependent and coexistence patterns. Yet, coarsening might be difficult to observe in real-world systems, as the coarsening process slows down as time progresses, and moreover, coarsening can cease due to decreased mobility of the organisms, as was observed for instance in mussels (42). Moreover, coarsening can be caused by various physical processes, as described in this paper and in other studies, and hence additional observations or experiments are needed to further solidify potential underlying processes. Identifying pattern-forming processes based on the dynamic properties of spatial patterning will provide an interesting but challenging avenue of research for future studies of pattern formation in ecology and beyond.

In this paper, we have extended the use of phase separation as a pattern-generating process from single-species examples [e.g., aggregation of mussels or ungulate herbivores (42, 80)] to universal two-species ecosystems, potentially expanding the applicability of phase separation within ecology. Our study demonstrates that within-patch bistability (e.g., patch collapse following a disturbance), between-patch competition, and hysteresis cause intrinsic vulnerability of patchy ecosystems governed by density-dependent aggregation. We have derived the feedbacks required for density-dependent aggregation and highlighted ecosystems with aggregation patterns, which could help identify these vulnerable ecosystems (Fig. 4

*C*and*D*). Statistics on size and spacing of patches could provide additional ways to distinguish between aggregation patterns and other types of patterns (11, 28, 108), even without knowledge of the exact pattern-forming processes. One of the most significant features of phase-separation patterns is the continuing coarsening, which distinguished this class of patterns from other forms of spatial topography, such as heterogeneities of the soil. Density-dependent aggregation generates a diversity of patterns because of the coarsening behavior. Those patterns share properties similar to at least a subset of coexistence patterns (11), but differ substantially in their dynamics. Whether and how these patterns are different and how they are mechanistically linked remains an outstanding question for future research. For instance, theory on Ostwald ripening (88, 110) and hyperuniformity (85, 86, 111) could guide such disordered-like statistics in future studies on the diversity of self-organized patterns.To summarize, our study contributes to a more complete indicator framework based on self-organized patchiness in ecology, providing a thorough understanding of irregularly patterned ecosystems that lack both scale-dependent feedbacks associated with spatially periodic Turing patterns (5, 16, 93) and self-propelling localized disturbances associated with scale-free patterns that typify self-organized criticality (112). Our extension of the theoretical framework of spatial indicators for tipping points in ecology can be used to provide tools for researchers and managers monitoring ecosystem degradation (61, 97).

## Materials and Methods

### Mathematical Analysis of Eq. **3**.

The steady states of Eq.

**3**comprise spatially uniform states and patterned states. The uniform steady states are given by combinations of $u$ and $v$ which satisfy $f(\stackrel{-}{u},\stackrel{-}{v})=0$ and $\overline{u}+\overline{v}=\langle \tau \rangle $, where $\stackrel{-}{u}$ and $\stackrel{-}{v}$ denote the steady states. Patterns form through diffusive instabilities of the uniform state at $\langle \tau \rangle ={\tau}_{D}$, which corresponds to$$-1<{\left.\frac{df\left(u,,,v\right)}{\mathit{du}}\right|}_{\left(\stackrel{-}{u},,,\stackrel{-}{v}\right)}<-\delta .$$

We provide a detailed linear stability analysis of Eq.

**3**on the patterned formation in*SI Appendix*, Texts S1–S3.### Numerical Analysis.

Most of the results presented in this paper are based on mathematical analysis of the generic model 3 (Eq.

**3**). For the analysis of the patch dynamics and patch statistics, we numerically solved a model with the specified reaction term $f(u,v)$:$$f\left(u,,,v\right)=uv-\frac{\alpha u}{u+\beta}.$$

The model was compared with a scale-dependent feedback model, similar to the model (5, 17)

$$\frac{\partial u}{\partial t}=uv-\frac{\alpha u}{u+\beta}+\delta \left(\frac{{\partial}^{2}u}{\partial {x}^{2}},+,\frac{{\partial}^{2}u}{\partial {y}^{2}}\right),$$

[8]

$$\frac{\partial v}{\partial t}=\gamma -v-uv+\left(\frac{{\partial}^{2}v}{\partial {x}^{2}},+,\frac{{\partial}^{2}v}{\partial {y}^{2}}\right).$$

[9]

### Model Implementation.

The models were implemented by discretizing $t$, $x,$ and $y$. The dimensionless time $t$ ran from 0 to 2,500, with a timestep of $\mathrm{\Delta}t=0.005$. Dimensionless space was implemented on a grid of $1,024\times 1,024$, representing a domain of $100\times 100$.

The models were solved in MATLAB R2019B and PyOpenCL using finite difference methods. All codes to produce the figures are available in the repository GitHub https://github.com/liuqx315/PS-Resilience-in-Ecosystems (113). The rates of change in $u$ and $v$ were calculated from forward differences, and the diffusion rates of $u$ and $v$ were calculated using the second-order central differences. The grids were loaded on a graphics processing unit (Nvidia Quadro K6000) to accelerate the computations.

Periodic boundaries were used to minimize boundary effects. The uniform steady state was used as initial condition for $u$ and $v$ for all model runs. To trigger pattern formation, a small amount of additive, spatially uncorrelated, normally distributed noise with a standard deviation of $\sigma ={10}^{-4}$ was added to both $u$ and $v$ at $t=0$.

The following parameter values were used for the model runs: $\alpha =4$, $\beta =1$, and $\delta =1/16$ for the density-dependent aggregation model (Eq.

**3**) with function $f(u,v)$, but $\beta =1$, $\gamma \in [2.1,3.5]$, and $\delta =0.01$ for the scale-dependent feedback model (Eqs.**8**and**9**). The global density of the mobile species was set to $\langle \tau \rangle =3$, except for the model runs of Fig. 5*C*and*D*.### Patch-Size Distributions.

The sizes of the patches in the model were determined using $u$ at $t=250$. A gridcell was said to be in a patch if $u{\stackrel{~}{u}}_{-}$, yielding a binary matrix. The areas of patches were determined by counting the number of connected grid cells and multiplying by the area of a single grid cell $\mathrm{\Delta}x\times \mathrm{\Delta}y$, where connected grid cells share a Von Neumann neighborhood, while acknowledging the periodic boundaries.

From the patch areas ${a}_{i,}$ a patch radius ${r}_{i}$ was estimated to enable comparison with theoretical predictions: ${r}_{i}=\sqrt{{a}_{i}/2\pi}$. Notice that, although this conversion only makes sense for circular patches, it does not affect the fit of the log-normal distribution, which best described the patch-size distributions for the high global densities at which non-circular patches were found. This is because if $a$ is log-normally distributed (as suggested by theory and empirical studies), then $\sqrt{a}$ is too, since $\mathrm{ln}\sqrt{a}=1/2\phantom{\rule{1pt}{0ex}}\mathrm{ln}\phantom{\rule{1pt}{0ex}}a$. In Figs. 3

*F*(*Inset*) and Fig. 5*C*, the patch sizes were scaled by dividing by the mean of the estimated patch radius.Two probability density functions were fitted to the obtained patch-size distributions. The first is the log-normal distribution:

$$f\left(r,;,\mu ,,,\sigma \right)=\frac{1}{r\sigma \sqrt{2\pi}}\mathrm{exp}\left(-,\frac{{(\mathrm{ln}r-\mu )}^{2}}{2{\sigma}^{2}}\right).$$

The model was parametrized using the maximum likelihood estimates:

$$\widehat{\mu}=\frac{1}{n}\sum _{i=0}^{n}\mathrm{ln}\phantom{\rule{2pt}{0ex}}{r}_{i},$$

[10]

$${\widehat{\sigma}}^{2}=\frac{1}{n}\sum _{i=0}^{n}{(\mathrm{ln}\phantom{\rule{2pt}{0ex}}{r}_{i}-\widehat{\mu})}^{2}.$$

[11]

The second model is a distribution that follows from Lifshitz–Slyozov–Wagner (LSW) theory for Ostwald ripening (87–89, 110):

$$g\left(r,;,\mu \right)=\frac{4{r}^{2}}{9{\mu}^{3}}{\left(\frac{3\mu}{3\mu +r}\right)}^{7/3}{\left(\frac{3\mu}{3\mu -2r}\right)}^{11/3}{e}^{\frac{2r}{2r-3\mu}},$$

[12]

for $0r3/2\mu $ and $g\left(r,;,\mu \right)=0$ for $r\ge 3/2\mu $, with $\mu $ being the mean patch radius. Following Tiryakioglu et al. (110), this model was parametrized using the maximum likelihood estimate of $\mu $, rather than the sample mean. The maximum likelihood estimates of $\mu $ are given by the solution of:

$$\begin{array}{c}\frac{\mathit{dL}}{d\mu}=-3n+\frac{7}{3}\sum _{i=0}^{n}\frac{{r}_{i}}{3\mu +{r}_{i}}-\frac{11}{3}\sum _{i=0}^{n}\frac{2{r}_{i}}{3\mu -2{r}_{i}}\hfill \\ \phantom{\rule{3.4em}{0ex}}+3\mu \sum _{i=0}^{n}\frac{2{r}_{i}}{{(2{r}_{i}-3\mu )}^{2}}=0.\hfill \end{array}$$

[13]

As the global density $\langle \tau \rangle $ decreases, the modeled patch-size distributions become left skewed (negative skewness, see Fig. 5

*C*). The log-normal distribution always has a positive skewness:$$\mathrm{\Gamma}=({e}^{{\widehat{\sigma}}^{2}}+2)\sqrt{{e}^{{\widehat{\sigma}}^{2}}-1},$$

and is not able to reproduce the change in skewness. The LSW distribution on the other hand has a constant negative skewness, $\mathrm{\Gamma}=-0.92002$, which approximates the simulated skewness at low global densities.

Both models were further compared by calculating weighted AIC values (wAIC), following Wagenmaker and Farrell (114). The wAIC can be interpreted as the probability that a distribution function is the best model, given the data and the set of candidate models.

## Data, Materials, and Software Availability

All codes to produce the figures are available in the repository GitHub https://github.com/liuqx315/PS-Resilience-in-Ecosystems. All study data are included in the article and/or

*SI Appendix*.## Acknowledgments

We thank Holly Moors, Bas Arens, and Maarten Eppinga that authorized the license of the pictures in Fig. 1

*A*,*C*, and*E*. We also thank Nigel Goldenfeld, Maarten Eppinga, and Mara Smeele for fruitful discussions and contributions. The research was supported by the EU Horizon 2020 project MERCES (689518), the project “Coping with deltas in transition” within the Program of Strategic Scientific Alliances between China and The Netherlands (PSA), financed by the Royal Dutch Academy for Arts and Sciences (Ref. PSA-SA-E-02), the Chinese Ministry of Science and Technology (2016YFE0133700), and the National Natural Science Foundation of China (32061143014, 41676084).### Author contributions

K.S., Q.-X.L., and J.v.d.K. designed research; K.S., Q.-X.L., and V.R. performed research; K.S., Q.-X.L., V.R., and A.D. analyzed the model; and K.S., Q.-X.L., V.R., T.v.d.H., M.R., A.D., C.B., and J.v.d.K. wrote the paper.

### Competing interest

The authors declare no competing interest.

## Supporting Information

Appendix 01 (PDF)

- Download
- 6.31 MB

Movie S1.

Dynamic behavior of spatial patterns results from scale-dependent feedback model Eqs.(8-9), where stationary patterns with saturated wavelength emerge for a long temporal evolution.

- Download
- 5.52 MB

Movie S2.

Dynamic behavior of spatial patterns results from the density-dependent aggregation model (Eq. 3), where the wavelength of patterns are unstable and coarsening with the temporal evolution.

- Download
- 5.94 MB

Movie S3.

The spatial dynamic behavior and spinodal bifurcation of phase-separation patterns within the (⟨τ⟩, δ)-parameter space.

- Download
- 6.43 MB

## References

1

M. Rietkerk, S. C. Dekker, P. C. de Ruiter, J. van de Koppel, Self-organized patchiness and catastrophic shifts in ecosystems.

*Science***305**, 1926–1929 (2004).2

V. Deblauwe, P. Couteron, O. Lejeune, J. Bogaert, N. Barbier, Environmental modulation of self-organized periodic vegetation patterns in Sudan.

*Ecography***34**, 990–1001 (2011).3

R. Bastiaansen, A. Doelman, M. B. Eppinga, M. Rietkerk, The effect of climate change on the resilience of ecosystems with adaptive spatial pattern formation.

*Ecol. Lett.***23**, 414–429 (2020).4

H. de Paoli et al., Behavioral self-organization underlies the resilience of a coastal ecosystem.

*Proc. Natl Acad. Sci. U.S.A.***114**, 8035–8040 (2017).5

C. A. Klausmeier, Regular and irregular patterns in semiarid vegetation.

*Science***284**, 1826–1828 (1999).6

E. Meron,

*Nonlinear Physics of Ecosystems*(CRC Press, Boca Raton, 2015), https://doi.org/10.1201/b18360.7

J. von Hardenberg, A. Y. Kletter, H. Yizhaq, J. Nathan, E. Meron, Periodic versus scale-free patterns in dryland vegetation.

*Proc. Roy. Soc. Biol. Sci.***277**, 1771–1776 (2010).8

S. Sankaran, S. Majumder, A. Viswanathan, V. Guttal, Clustering and correlations: Inferring resilience from spatial patterns in ecosystems.

*Methods Ecol. Evol.***10**, 2079–2089 (2019).9

M. Scheffer et al., Early-warning signals for critical transitions.

*Nature***461**, 53–59 (2009).10

Q. X. Liu et al., Pattern formation at multiple spatial scales drives the resilience of mussel bed ecosystems.

*Nat. Commun.***5**, 5234 (2014).11

M. Rietkerk et al., Evasion of tipping in complex systems through spatial pattern formation.

*Science***374**, eabj0359 (2021).12

L. X. Zhao et al., Fairy circles reveal the resilience of self-organized salt marshes.

*Sci. Adv.***7**, eabe1100 (2021).13

M. Rietkerk et al., Self-organization of vegetation in arid ecosystems.

*Am. Nat.***160**, 524–530 (2002).14

H. Meinhardt,

*The Algorithmic Beauty of Sea Shells*(Springer-Verlag, Berlin, Heidelberg, New York, 2009), p. 269, https://doi.org/10.1007/978-3-540-92142-4.15

J. D. Murray, “Mathematical biology II: Spatial models and biomedical applications”.

*Interdisciplinary Applied Mathematics*(Springer, New York, NY, 2003),**vol. 18**, p. 1.16

M. Rietkerk, J. Van de Koppel, Regular pattern formation in real ecosystems.

*Trends Ecol. Evol.***23**, 169–175 (2008).17

J. van de Koppel, M. Rietkerk, N. Dankers, P. M. J. Herman, Scale-dependent feedback and regular spatial patterns in young mussel beds.

*Am. Nat.***165**, E66–E77 (2005).18

R. Lefever, O. Lejeune, On the origin of tiger bush.

*Bull. Math. Biol.***59**, 263–294 (1997).19

S. Getzin et al., Discovery of fairy circles in Australia supports self-organization theory.

*Proc. Natl. Acad. Sci. U.S.A.***113**, 3551–3556 (2016).20

J. von Hardenberg, E. Meron, M. Shachak, Y. Zarmi, Diversity of vegetation patterns and desertification.

*Phys. Rev. Lett.***87**, 198101 (2001).21

J. van de Koppel, C. M. Crain, Scale-dependent inhibition drives regular tussock spacing in a freshwater marsh.

*Am. Nat.***168**, E136–E147 (2006).22

T. van der Heide et al., Spatial self-organized patterning in seagrasses along a depth gradient of an intertidal ecosystem.

*Ecology***91**, 362–369 (2010).23

M. Pascual, F. Guichard, Criticality and disturbance in spatial ecological systems.

*Trends Ecol. Evol.***20**, 88–95 (2005).24

Z. Ge, Q.-X. Liu, Foraging behaviours lead to spatiotemporal self-similar dynamics in grazing ecosystems.

*Ecol. Lett.***25**, 378–390 (2022).25

T. van der Heide et al., Ecosystem engineering by seagrasses interacts with grazing to shape an intertidal landscape.

*PLoS One***7**, e42060 (2012).26

G. Berg, P. Esselink, M. Groeneweg, K. Kiehl, Micropatterns in Festuca rubra-dominated salt-marsh vegetation induced by sheep grazing.

*Plant Ecol.***132**, 1–14 (1997).27

S. Kéfi et al., Spatial vegetation patterns and imminent desertification in Mediterranean arid ecosystems.

*Nature***449**, 213–217 (2007).28

M. E. Ritchie, H. Olff, Spatial scaling laws yield a synthetic theory of biodiversity.

*Nature***400**, 557–560 (1999).29

T. K. Michaels, M. B. Eppinga, J. D. Bever, A nucleation framework for transition between alternate states: Short-circuiting barriers to ecosystem recovery.

*Ecology***101**, e03099 (2020).30

X. Dong, A. B. Murray, J. B. Heffernan, Ecohydrologic feedbacks controlling sizes of cypress wetlands in a patterned karst landscape.

*Earth Surf. Process. Landforms***44**, 1178–1191 (2019).31

V. C. Reijers et al., A Lévy expansion strategy optimizes early dune building by beach grasses.

*Nat. Commun.***10**, 2656 (2019).32

E. J. Weerman et al., Changes in diatom patch-size distribution and degradation in a spatially self-organized intertidal mudflat ecosystem.

*Ecology***93**, 608–618 (2012).33

C. Latorre et al., Establishment and formation of fog-dependent Tillandsia landbeckii dunes in the Atacama Desert: Evidence from radiocarbon and stable isotopes.

*J. Geophys. Res. Biogeosci.***116**, G03033 (2011).34

M. B. Eppinga, P. C. de Ruiter, M. J. Wassen, M. Rietkerk, Nutrients and hydrology indicate the driving mechanisms of peatland surface patterning.

*Am. Nat.***173**, 803–818 (2009).35

X. Dong, N. B. Grimm, J. B. Heffernan, R. Muneepeerakul, Interactions between physical template and self-organization shape plant dynamics in a stream ecosystem.

*Ecosystems***23**, 891–905 (2020).36

L. X. Zhao, C. Xu, Z. M. Ge, J. van de Koppel, Q. X. Liu, The shaping role of self-organization: Linking vegetation patterning, plant traits and ecosystem functioning.

*Proc. R. Soc. Biol. Sci.***286**, 20182859 (2019).37

J. A. Castillo Vardaro et al., Resource availability and heterogeneity shape the self-organisation of regular spatial patterning.

*Ecol. Lett.***24**, 1880–1891 (2021).38

K. E. Anderson, F. M. Hilker, R. M. Nisbet, Directional biases and resource-dependence in dispersal generate spatial patterning in a consumer–producer model.

*Ecol. Lett.***15**, 209–217 (2012).39

D. R. Foster, G. A. King, P. H. Glaser, H. E. Wright, Origin of string patterns in boreal peatlands.

*Nature***306**, 256–258 (1983).40

D. Horváth, V. Petrov, S. K. Scott, K. Showalter, Instabilities in propagating reaction-diffusion fronts.

*J. Chem. Phys.***98**, 6332–6343 (1993).41

A. J. Bray, C. L. Emmott, Lifshitz-Slyozov scaling for late-stage coarsening with an order-parameter-dependent mobility.

*Phys. Rev. B***52**, R685–R688 (1995).42

Q. X. Liu et al., Phase separation explains a new class of self-organized spatial patterns in ecological systems.

*Proc. Natl. Acad. Sci. U.S.A.***110**, 11905–11910 (2013).43

D. Zwicker, A. A. Hyman, F. Jülicher, Suppression of Ostwald ripening in active emulsions.

*Phys. Rev. E Stat. Nonlin. Soft Matter Phys.***92**, 012317 (2015).44

F. Liu, N. Goldenfeld, Dynamics of phase separation in block copolymer melts.

*Phys. Rev. A***39**, 4805–4810 (1989).45

A. A. Hyman, C. A. Weber, F. Jülicher, Liquid-liquid phase separation in biology.

*Annu. Rev. Cell Dev. Biol.***30**, 39–58 (2014).46

J. A. Riback et al., Composition-dependent thermodynamics of intracellular phase separation.

*Nature***581**, 209–214 (2020).47

J. Stenhammar, A. Tiribocchi, R. J. Allen, D. Marenduzzo, M. E. Cates, Continuum theory of phase separation kinetics for active Brownian particles.

*Phys. Rev. Lett.***111**, 145702 (2013).48

K. Takae, T. Kawasaki, Emergent elastic fields induced by topological phase transitions: Impact of molecular chirality and steric anisotropy.

*Proc. Natl. Acad. Sci. U.S.A.***119**, e2118492119 (2022).49

A. Klosin et al., Phase separation provides a mechanism to reduce noise in cells.

*Science***367**, 464–468 (2020).50

A. Radja, E. M. Horsley, M. O. Lavrentovich, A. M. Sweeney, Pollen cell wall patterns form from modulated phases.

*Cell***176**, 856–868.e810 (2019).51

S. Alberti, A. Gladfelter, T. Mittag, Considerations and challenges in studying liquid-liquid phase separation and biomolecular condensates.

*Cell***176**, 419–434 (2019).52

J. Wang et al., A molecular grammar governing the driving forces for phase separation of prion-like RNA binding proteins.

*Cell***174**, 688–699.e616 (2018).53

R. Lefever, N. Barbier, P. Couteron, O. Lejeune, Deeply gapped vegetation patterns: On crown/root allometry, criticality and desertification.

*J. Theor. Biol.***261**, 194–209 (2009).54

E. Meron, Pattern-formation approach to modelling spatially extended ecosystems.

*Ecol. Modell.***234**, 70–82 (2012).55

J. A. Sherratt, Using wavelength and slope to infer the historical origin of semiarid vegetation bands.

*Proc. Natl. Acad. Sci U.S.A.***112**, 4202–4207 (2015).56

J. van de Koppel et al., Experimental evidence for spatial self-organization and its emergent effects in mussel bed ecosystems.

*Science***322**, 739–742 (2008).57

C. E. Tarnita et al., A theoretical foundation for multi-scale regular vegetation patterns.

*Nature***541**, 398–401 (2017).58

R. M. Pringle, E. T. Tarnita, Spatial self-organization of ecosystems: Integrating multiple mechanisms of regular-pattern formation.

*Annu. Rev. Entomol.***62**, 359–377 (2017).59

H. Jansen, L. van den Bogaart,

*Blue Carbon by Marine Bivalves*(Wageningen University, The Netherland, 2020), https://doi.org/10.18174/537188.60

S. J. McNaughton, Grazing lawns: Animals in herds, plant form, and coevolution.

*Am. Nat.***124**, 863–886 (1984).61

I. Stavi, H. Yizhaq, Y. Osem, E. Argaman, Positive impacts of livestock and wild ungulate routes on functioning of dryland ecosystems.

*Ecol. Evol.***11**, 13684–13691 (2021).62

W. H. van der Putten, P. W. T. Maas, W. J. M. Van Gulik, H. Brinkman, Characterization of soil organisms involved in the degeneration of Ammophila arenaria.

*Soil. Biol. Biochem.***22**, 845–852 (1990).63

O. Durán Vinent, L. J. Moore, Barrier island bistability induced by biophysical interactions.

*Nat. Clim. Change***5**, 158–162 (2014).64

J. Quets, “A spatiotemporal analysis of desert ecosystems with phytogenic mounds (nebkhas),” PhD thesis, Universiteit Antwerpen (2015).

65

J. M. d. Goeij et al., Surviving in a marine desert: The sponge loop retains resources within coral reefs.

*Science***342**, 108–110 (2013).66

Z. Dubinsky, I. Berman-Frank, Uncoupling primary production from population growth in photosynthesizing organisms in aquatic ecosystems.

*Aquat. Sci.***63**, 4–17 (2001).67

C. E. Nelson, A. L. Alldredge, E. A. McCliment, L. A. Amaral-Zettler, C. A. Carlson, Depleted dissolved organic carbon and distinct bacterial communities in the water column of a rapid-flushing coral reef ecosystem.

*ISME J.***5**, 1374–1387 (2011).68

A. F. Haas et al., Effects of coral reef benthic primary producers on dissolved organic carbon and microbial activity.

*PLoS One***6**, e27973 (2011).69

C. Wild et al., Coral mucus functions as an energy carrier and particle trap in the reef ecosystem.

*Nature***428**, 66–70 (2004).70

D. Ruiz-Reynes et al., Fairy circle landscapes under the sea.

*Sci. Adv.***3**, e1603262 (2017).71

É. G. Devoie et al., Mechanisms of discontinuous permafrost thaw in peatlands.

*J. Geophys. Res. Earth Surf.***126**, e2021JF006204 (2021).72

M. Rietkerk, S. C. Dekker, M. J. Wassen, A. W. M. Verkroost, M. F. P. Bierkens, A putative mechanism for bog patterning.

*Am. Nat.***163**, 699–708 (2004).73

F. Brauns, J. Halatek, E. Frey, Phase-space geometry of mass-conserving reaction-diffusion dynamics.

*Phys. Rev. X.***10**, 041036 (2020).74

F. Brauns, H. Weyer, J. Halatek, J. Yoon, E. Frey, Wavelength selection by interrupted coarsening in reaction-diffusion systems.

*Phys. Rev. Lett.***126**, 104101 (2021).75

J. W. Cahn, J. E. Hilliard, Free energy of a nonuniform system. I. Interfacial free energy.

*J. Chem. Phys.***28**, 258–267 (1958).76

J. W. Cahn, Phase separation by spinodal decomposition in isotropic systems.

*J. Chem. Phys.***42**, 93–99 (1965).77

J. Halatek, E. Frey, Rethinking pattern formation in reaction–diffusion systems.

*Nat. Phys.***14**, 507–514 (2018).78

Y. Mori, A. Jilkine, L. Edelstein-Keshet, Asymptotic and bifurcation analysis of wave-pinning in a reaction-diffusion model for cell polarization.

*SIAM J. Appl. Math.***71**, 1401–1427 (2011).79

N. Verschueren, A. Champneys, A model for cell polarization without mass conservation.

*SIAM J. Appl. Dyn. Syst.***16**, 1797–1830 (2017).80

Q. X. Liu et al., Phase separation driven by density-dependent movement: A novel mechanism for ecological patterns.

*Phys. Life Rev.***19**, 107–121 (2016).81

L. McNally et al., Killing by type VI secretion drives genetic phase separation and correlates with increased cooperation.

*Nat. Commun.***8**, 14371 (2017).82

Y. Chen, J. E. Ferrell Jr.,

*C. elegans*colony formation as a condensation phenomenon.*Nat. Commun.***12**, 4947 (2021).83

T. M. Scanlon, K. K. Caylor, S. A. Levin, I. Rodriguez-Iturbe, Positive feedbacks promote power-law clustering of Kalahari vegetation.

*Nature***449**, 209–212 (2007).84

S. Kéfi et al., Robust scaling in ecosystems and the meltdown of patch size distributions before extinction.

*Ecol. Lett.***14**, 29–35 (2011).85

M. Huang, W. Hu, S. Yang, Q. X. Liu, H. P. Zhang, Circular swimming motility and disordered hyperuniform state in an algae system.

*Proc. Natl. Acad. Sci. U.S.A.***118**, e2100493118 (2021).86

D. Hexner, P. M. Chaikin, D. Levine, Enhanced hyperuniformity from random reorganization.

*Proc. Natl. Acad. Sci. U.S.A.***114**, 4294–4299 (2017).87

C. Wagner, Theorie der alterung von niederschlägen durch umlösen (Ostwald-Reifung).

*Zeitschrift für Elektrochemie, Berichte der Bunsengesellschaft für physikalische Chemie***65**, 581–591 (1961).88

J. H. Yao, K. R. Elder, H. Guo, M. Grant, Theory and simulation of Ostwald ripening.

*Phys. Rev. B***47**, 14110–14125 (1993).89

I. M. Lifshitz, V. V. Slyozov, The kinetics of precipitation from supersaturated solid solutions.

*J. Phys. Chem. Solids***19**, 35–50 (1961).90

A. Li et al., Ice needles weave patterns of stones in freezing landscapes.

*Proc. Natl. Acad. Sci. U.S.A.***118**, e2110670118 (2021).91

R.-H. Wang, Q.-X. Liu, G.-Q. Sun, Z. Jin, J. van de Koppel, Nonlinear dynamic and pattern bifurcations in a model for spatial patterns in young mussel beds.

*J. R. Soc. Interface***6**, 705–718 (2009).92

K. Siteur et al., Beyond Turing: The response of patterned ecosystems to environmental change.

*Ecol. Complex***20**, 81–96 (2014).93

A. M. Turing, The chemical basis of morphogenesis.

*Phil. Tran. Roy. Soc. Lond. Series. B***237**, 37–72 (1952).94

W. F. Smith, J. Hashemi,

*Foundations of Materials Science and Engineering*(McGraw-Hill Series in Materials Science, vol.**509**, McGraw-Hill Higher Education, Boston, ed. 4, 2006), p. 21.95

D. V. Alexandrov, I. V. Alexandrova, From nucleation and coarsening to coalescence in metastable liquids.

*Trans. R. Soc. A Math. Phys. Eng. Sci.***378**, 20190247 (2020).96

I. V. Alexandrova, D. V. Alexandrov, E. V. Makoveeva, Ostwald ripening in the presence of simultaneous occurrence of various mass transfer mechanisms: An extension of the Lifshitz-Slyozov theory.

*Philos. Trans. R. Soc. A Math. Phys. Eng. Sci.***379**, 20200308 (2021).97

A. Génin et al., Monitoring ecosystem degradation using spatial data and the R package spatial warnings.

*Methods Ecol. Evol.***9**, 2067–2075 (2018).98

S. Majumder, K. Tamma, S. Ramaswamy, V. Guttal, Inferring critical thresholds of ecosystem transitions from spatial data.

*Ecology***100**, e02722 (2019).99

J. Rockström et al., A safe operating space for humanity.

*Nature***461**, 472–475 (2009).100

W. Steffen et al., Planetary boundaries: Guiding human development on a changing planet.

*Science***347**, 1259855 (2015).101

B. R. Silliman et al., Facilitation shifts paradigms and can amplify coastal restoration efforts.

*Proc. Natl. Acad. Sci. U.S.A.***112**, 14295–14300 (2015).102

F. Meloni, C. R. F. Granzotti, S. Bautista, A. S. Martinez, Scale dependence and patch size distribution: Clarifying patch patterns in Mediterranean drylands.

*Ecosphere***8**, e01690 (2017).103

C. Xu et al., Can we infer plant facilitation from remote sensing? A test across global drylands.

*Ecol. Appl.***25**, 1456–1462 (2015).104

A. I. Borthagaray, M. A. Fuentes, P. A. Marquet, Vegetation pattern formation in a fog-dependent ecosystem.

*J. Theor. Biol.***265**, 18–26 (2010).105

J. A. Sherratt, Q. X. Liu, J. van de Koppel, A comparison of the “reduced losses” and “increased production” models for mussel bed dynamics.

*Bull. Math. Biol.***83**, 99 (2021).106

T. Vicsek, Ecological patterns emerging as a result of the density distribution of organisms.

*Phys. Life Rev.***19**, 139–141 (2016).107

M. E. Cates, D. Marenduzzo, I. Pagonabarraga, J. Tailleur, Arrested phase separation in reproducing bacteria creates a generic route to pattern formation.

*Proc. Natl. Acad. Sci. U.S.A.***107**, 11715–11720 (2010).108

M. E. Ritchie,

*Scale, Heterogeneity, and the Structure and Diversity of Ecological Communities*(Princeton University Press, 2009).109

F. Guichard, P. M. Halpin, G. W. Allison, J. Lubchenco, B. A. Menge, Mussel disturbance dynamics: Signatures of oceanographic forcing from local interactions.

*Am. Nat.***161**, 889–904 (2003).110

M. Tiryakioğlu, G. Ökten, D. Hudak, R. T. Shuey, J. P. Suni, On evaluating fit of the Lifshitz–Slyozov–Wagner (LSW) distribution to particle size data.

*Mater. Sci. Eng. A***527**, 1636–1639 (2010).111

Z. Ma, S. Torquato, Random scalar fields and hyperuniformity.

*J. Appl. Phys.***121**, 244904 (2017).112

R. V. Solé, J. Bascompte,

*Self-Organization in Complex Ecosystems*(Monographs in Population Biology 42, Princeton University Press, Princeton, 2006), p. 373.113

Siteur et al., “Phase-separation physics underlies new theory for the resilience of patchy ecosystems.” GitHub. https://github.com/liuqx315/PS-Resilience-in-Ecosystems. Deposited 2 June 2022

114

E.-J. Wagenmakers, S. Farrell, AIC model selection using Akaike weights.

*Psychon. Bull. Rev.***11**, 192–196 (2004).## Information & Authors

### Information

#### Published in

#### Classifications

#### Copyright

Copyright © 2023 the Author(s). Published by PNAS. This open access article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

#### Data, Materials, and Software Availability

All codes to produce the figures are available in the repository GitHub https://github.com/liuqx315/PS-Resilience-in-Ecosystems. All study data are included in the article and/or

*SI Appendix*.#### Submission history

**Received**: February 19, 2022

**Accepted**: November 22, 2022

**Published online**: January 3, 2023

**Published in issue**: January 10, 2023

#### Keywords

#### Acknowledgments

We thank Holly Moors, Bas Arens, and Maarten Eppinga that authorized the license of the pictures in Fig. 1

*A*,*C*, and*E*. We also thank Nigel Goldenfeld, Maarten Eppinga, and Mara Smeele for fruitful discussions and contributions. The research was supported by the EU Horizon 2020 project MERCES (689518), the project “Coping with deltas in transition” within the Program of Strategic Scientific Alliances between China and The Netherlands (PSA), financed by the Royal Dutch Academy for Arts and Sciences (Ref. PSA-SA-E-02), the Chinese Ministry of Science and Technology (2016YFE0133700), and the National Natural Science Foundation of China (32061143014, 41676084).##### Author Contributions

K.S., Q.-X.L., and J.v.d.K. designed research; K.S., Q.-X.L., and V.R. performed research; K.S., Q.-X.L., V.R., and A.D. analyzed the model; and K.S., Q.-X.L., V.R., T.v.d.H., M.R., A.D., C.B., and J.v.d.K. wrote the paper.

##### Competing Interest

The authors declare no competing interest.

#### Notes

This article is a PNAS Direct Submission.

### Authors

## Metrics & Citations

### Metrics

#### Citation statements

#### Altmetrics

### Citations

If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

#### Cited by

Loading...

## View Options

### View options

#### PDF format

Download this article as a PDF file

DOWNLOAD PDF### Get Access

#### Login options

Check if you have access through your login credentials or your institution to get full access on this article.

Personal login Institutional Login#### Recommend to a librarian

Recommend PNAS to a Librarian#### Purchase options

Purchase this article to get full access to it.