Elsevier

Physics Reports

Volume 1088, 17 October 2024, Pages 1-41
Physics Reports

Stability of ecological systems: A theoretical review

https://doi.org/10.1016/j.physrep.2024.08.001Get rights and content

Abstract

The stability of ecological systems is a fundamental concept in ecology, which offers profound insights into species coexistence, biodiversity, and community persistence. In this article, we provide a systematic and comprehensive review on the theoretical frameworks for analyzing the stability of ecological systems. Notably, we survey various stability notions, including linear stability, sign stability, diagonal stability, D-stability, total stability, sector stability, and structural stability. For each of these stability notions, we examine necessary or sufficient conditions for achieving such stability and demonstrate the intricate interplay of these conditions on the network structures of ecological systems. We further discuss the stability of ecological systems with higher-order interactions.

Keywords

Stability
Ecological systems
Lyapunov theory
Network structure
Random matrix theory
Generalized Lotka–Volterra model
Consumer-resource models
Higher-order interactions

1. Introduction

Ecological systems are usually subject to continual disturbances or perturbations. Here, disturbances often refer to large, infrequent events that can cause significant and lasting changes to the ecological system. They can reshape the ecological landscape and lead to a new equilibrium. Perturbations typically refer to small temporary changes that test the system’s response and resilience. They generally do not result in long-term changes to the overall structure of the ecological system. The responses of ecological systems to disturbances and perturbations are often characterized quantitatively as stability, a fundamental concept rooted in dynamical systems and control theory [1], [2], [3], [4], [5], [6], [7], [8], [9]. Loosely speaking, an ecological system is stable if its state represented by a vector of its species abundances, do not change too much under small perturbations. Therefore, understanding the stability of ecological systems is particularly important due to its implications for species coexistence, biodiversity, and community persistence [10], [11], [12], [13], [14], [15], [16], [17]. More importantly, these implications can contribute to the development of effective interventions and policies that maintain the overall health of ecological systems [18], [19], [20], [21].
Numerous mathematicians, physicists, and ecologists have explored the stability of ecological systems. In his pioneer work, May linearized an unspecified nonlinear dynamical system around an equilibrium to perform the local linear stability analysis of ecological systems [22], [23], [24]. The local linear stability is completely characterized by the eigenvalues of the Jacobian matrix (of the unspecified dynamical model), which is often called the community matrix in the ecological contexts. If all the eigenvalues of the community matrix have negative real parts, the ecological system is locally stable around the equilibrium. By incorporating random matrix theory, May successfully derived a simple condition to characterize the stability of ecological systems. While there have been intense debates about the relationship between stability and complexity over the past five decades [25], [26], [27], [28], [29], [30], [31], [32], May’s result offers a powerful foundation for analyzing the stability of ecological systems with diverse structural characteristics, such as interaction types [33], degree heterogeneity [34], self-regulation [35], modularity [36], and more.
May’s framework does not need to know the dynamical model of the ecological system. Instead, it focuses on the Jacobian matrix of the unspecified dynamical model. Hence, the framework is model-implicit. In the past, actually many ecological models, such as the generalized Lotka–Volterra (GLV) model (1910) [37], [38], Holling type II model (1959) [39], [40], and MacArthur’s consumer-resource model (1970) [41], [42], have been proposed to study the dynamics of ecological systems. Remarkably, the GLV model is the most commonly used model in the stability analysis of ecological systems due to its simplicity [43], [44], [45], [46], [47], [48], [49], [50]. Most of these findings rely on Lyapunov theory, which allows us to determine the stability of a system without explicitly integrating the differential equation. Lyapunov theory encompasses two fundamental methods. One is based on linearization (e.g., May’s approach), which is often referred to as the Lyapunov indirect method, while the other is called the Lyapunov direct method. The essence of Lyapunov direct method is to construct the so-called Lyapunov function, an “energy-like” scalar function whose time variation can be viewed as the “energy dissipation”. It is not restricted to small perturbations, and in principle can be applied to all dynamical systems. However, we lack a general theory to find a suitable Lyapunov function for an arbitrary dynamical system. Instead, we have to rely on our experience and intuition. Notably, certain diagonal-type Lyapunov functions are useful for the stability analysis of ecological systems described by the GLV model [51].
The feasibility of an ecological system, which concerns whether a specific set of conditions or parameters can support the system over time, can also impact its stability [28], [52], [53], [54], [55], [56], [57]. Intriguingly, the relationship between feasibility and stability is noteworthy when dealing with ecological systems described by the GLV model or MacArthur’s consumer-resource model. It has been proved that under certain conditions, the feasibility of such systems can imply their stability [28], [57], [58]. Moreover, feasibility is closely related to a concept known as structural stability, which reflects the ability of an ecological system to qualitatively maintain its dynamics and overall behavior under small perturbations of the system model itself [59]. A great amount of metrics have been proposed to quantify the structural stability of ecological systems through feasibility analysis [55], [59], [60].
  1. Download: Download high-res image (224KB)
  2. Download: Download full-size image

Fig. 1. Pairwise interactions versus higher-order interactions in ecological systems. (a) Pairwise interspecies interactions. (b) Higher-order interspecies interactions. (c) Consumer-resource interactions naturally implies higher-order interspecies interactions.

Recently, there has been a growing interest in higher-order interactions within ecological systems, which refers to the intricate relationships that extend beyond pairwise interspecies interactions [61], [62], [63], [64], [65], see Fig. 1a, b. Mathematically, these higher-order interactions emerged in an ecological system can be naturally represented as a hypergraph, where its hyperedges can connect an arbitrary number of nodes [66]. The resulting dynamics can be expressed in the form of polynomials [67]. The analysis of hypergraph dynamics often involves the application of tensor theory [67], [68], [69], [70], which deals with multidimensional arrays generalized from vectors and matrices [71], [72], [73]. Lyapunov theory can also be leveraged to analyze the stability of polynomial dynamical systems [74]. Furthermore, the consumer-resource model can be implicitly considered as an instance of higher-order interactions, with consumer species as nodes and resources as hyperedges, see Fig. 1c. Increasing theoretical evidence has revealed that higher-order interactions play a significant role in the stability of ecological systems [61], [63], [64].
Several valuable reviews have explored specific aspects of ecological stability analysis. For example, Cui et al. investigated the connections between ecology and statistical physics, focusing on techniques like the cavity method and random matrix theory [124]. In addition, Akjouj et al. examined the intersection of theoretical ecology and large random matrix theory [125]. Furthermore, Landi et al. addressed the relationship between complexity and stability in different types of models and real ecological systems [126]. However, these reviews primarily delve into specific areas of stability, making it challenging to gain a comprehensive understanding of the diverse stability criteria used in ecological research and the progress made within each area. This article presents a systematic and inclusive analysis of the theoretical frameworks employed to assess various stability types in ecological systems, along with their relationships to underlying network structures. In Section 2, we briefly review the fundamentals of Lyapunov theory, which acts as a foundation for the stability analysis of ecological systems. From Section 3 to 9, we comprehensively survey various stability notions, including linear stability, sign stability, diagonal stability, D-stability, total stability, sector stability, and structural stability, see Table 1. For each of these stability notions, we examine necessary or sufficient conditions that are required to establish such stability and illustrate the complicated relationships between these conditions and the network structures of ecological systems. We further delve into the stability of ecological systems with higher-order interactions in Section 10. Finally, we discuss the future prospects of the stability notions in Section 11 and conclude in Section 12.

Table 1. Summary of stability notions that have been used to study ecological systems.

NotionMain modelCharacteristicsKey references
Linear
Stability
Linear model (3),
GLV model (9)
Characterizing the eigenvalue spectrum of the community matrix M through random matrix theory.[22], [23], [33], [34], [35], [36], [57], [75], [76], [77], [78], [79], [80], [81], [82], [83], [84], [85], [86], [87], [88]
Sign
Stability
Linear model (3)Assessing stability based solely on the signs of the elements in the community matrix M without considering their magnitudes.[89], [90], [91], [92], [93], [94], [95], [96]
Diagonal StabilityGLV model (9)The diagonal stability of the interaction matrix A, i.e., there exists a diagonal matrix P such that AP+PA0, implies the Lyapunov stability of the GLV model.[45], [46], [48], [49], [51], [55], [79], [97], [98], [99], [100], [101], [102], [103], [104], [105], [106], [107]
D-StabilityGLV model (9)The D-stability of the interaction matrix A, i.e., the matrix XA is stable for any positive diagonal matrix X, implies that any feasible equilibrium of the GLV model is locally stable.[59], [79], [108], [109], [110], [111]
Total
Stability
GLV model (9)The total stability of the interaction matrix A, i.e., any principal sub-matrix of A is D-stable, implies that a new feasible equilibrium of the GLV model remains locally stable after the removal or extinction of certain species.[95], [112], [113], [114]
Sector
Stability
Generic population dynamics model (25)A semi-feasible equilibrium point is called sector stable if every solution that starts within a non-negative neighborhood remains in the same or an even larger non-negative neighborhood and eventually converges to that equilibrium point.[58], [115]
Structural StabilityGLV model (9)Measuring the capacity to qualitatively maintains the dynamics under small perturbations through feasibility analysis of the GLV model.[55], [59], [60], [116], [117], [118], [119], [120], [121], [122], [123]

2. Preliminaries

In this section, we review the fundamentals of Lyapunov theory, which plays a significant role in the various stability analyses of ecological systems. We adapt most of the notations and results from the comprehensive work of [127], [128], [129], [130], [131].

2.1. Definitions of Lyapunov stability

Consider a general S-dimensional nonlinear system (1)ẋ(t)=f(x(t),t)withx(t0)=x0,x(t)RS.In most cases, we are interested in the stability of the system near an equilibrium point xRS that satisfies f(x,t)=0. Without loss of generality, we can shift the origin of the system and assume that the equilibrium point of interest occurs at x=0. Let M be the Jacobian matrix evaluated at x. Denote n, n0, and n+ be the numbers of eigenvalues (counting multiplicities) of M with negative, zero, and positive real part, respectively. An equilibrium point is called hyperbolic if n0=0, i.e., there are no eigenvalues on the imaginary axis. A hyperbolic equilibrium point is called a hyperbolic saddle if nn+0. Since a generic matrix has n0=0, equilibrium points in a generic dynamical system are typically hyperbolic.
An equilibrium point x=0 is stable (in the sense of Lyapunov) at t=t0 if ϵ>0, δ(ϵ,t0)>0 such that x(t0)δ(ϵ,t0)x(t)ϵ,tt0>0.Otherwise, x is unstable, see Fig. 2. An equilibrium point x=0 is asymptotically stable at t=t0 if it is stable and locally attractive, i.e., δ(t0)>0 such that x(t0)δ(t0)limtx(t)=0.An equilibrium point which is stable but not asymptotically stable is called marginally stable, see Fig. 2. Both stability and asymptotic stability are defined at a time instant t0. In practice, it is often desirable for a system to have a certain uniformity in its behavior. Uniform (asymptotic) stability requires that the equilibrium point is (asymptotically) stable for all t0>0. Note that notions of uniformity are only relevant for time-varying or non-autonomous systems. For autonomous or time-invariant systems, (asymptotic) stability naturally implies uniform (asymptotic) stability. The definition of asymptotic stability does not quantify the rate of convergence. The notion of exponential stability guarantees a minimal rate of convergence. An equilibrium point x=0 is exponentially stable if δ,α,λ>0 such that x(t0)δx(t)αx(t0)eλ(tt0),tt0.Exponential stability is a very strong form of stability because it implies uniform asymptotic stability.
The above definitions of stability, asymptotic stability, and exponential stability are local and only describe the behavior of a system near an equilibrium point. We call an equilibrium point globally (asymptotically or exponentially) stable, if it is (asymptotically or exponentially) stable for all initial conditions x0RS. Though global stability is very desirable, it is very difficult to achieve in many systems. Note that linear time-invariant systems are either asymptotically stable, or marginally stable, or unstable. Moreover, linear asymptotic stability is always global and exponential, and linear instability always implies exponential blow-up. The above refined notions of stability are explicitly needed only for nonlinear systems.
  1. Download: Download high-res image (212KB)
  2. Download: Download full-size image

Fig. 2. Geometrical implication of stability. For nonlinear systems ẋ(t)=f(x(t),t), due to the complex and rich behavior of nonlinear dynamics, various types of stability, e.g., stability, asymptotic stability, and global asymptotic stability, can be discussed. Intuitively and roughly speaking, if all solutions of the system that start out near an equilibrium point x stay near x forever, then x is stable (in the sense of Lyapunov). More strongly, if x is stable and all solutions that start out near x converge to x as t, then x is asymptotically stable. x is marginally stable if it is stable but not asymptotically stable.

This figure was redrawn from [131].

2.2. Lyapunov’s indirect method

Lyapunov’s indirect method is based on linearization, and is concerned with the local stability of a nonlinear system and is based on the intuition that a nonlinear system should behave similarly to its linearized approximation for small perturbations The linearization can be used to give a conservative bound on the domain of attraction of the equilibrium point for the original nonlinear system.
For a general nonlinear system (1) with continuously differentiable f(x(t),t) around the equilibrium point x=0, the system dynamics can be rewritten as ẋ(t)=M(t)x(t)+fh.o.t(x(t),t),where M(t)=f(x,t)x|x=x is the Jacobian matrix of f(x,t) with respect to x, evaluated at the equilibrium point x, and fh.o.t(x,t) represents the higher-order terms in x. We typically require fh.o.t(x,t) approaches zero uniformly, which is obviously true for an autonomous system. Then ẋ(t)=M(t)x(t) is called the uniform linearization of the original nonlinear system (1) at the equilibrium point x=0. If this uniform linearization exists and the Jacobian matrix M(t) is bounded, then the uniform asymptotic stability of the equilibrium point x=0 for the linearization implies its uniform local asymptotic stability for the original nonlinear system.
For an autonomous system ẋ(t)=f(x(t)), the linearization is simply an linear time-invariant system ẋ(t)=Mx(t). Denote the eigenvalues of M as λ(M). We have the following results: (i) if the linearization is strictly stable (i.e., all the eigenvalues of M lie in the closed left half of the complex plane), then the equilibrium point x=0 is asymptotically stable for the original nonlinear system (note that a real square matrix M that satisfies Re[λ(M)]<0 is often called Hurwitz stable); (ii) if the linearization is unstable (i.e., at least one eigenvalue of M has positive real part), then the equilibrium point is unstable for the original nonlinear system; (iii) if the linearization is marginally stable (i.e, all eigenvalues of M are in the left-half complex plane, but at least one of them is on the imaginary axis), then one cannot conclude anything from the linearization—the equilibrium point may be stable, asymptotically stable, or unstable for the original nonlinear system.

2.3. Lyapunov’s direct method

Lyapunov’s direct method (also called the second method of Lyapunov) can be considered as a mathematical extension of a fundamental physical observation. If the total energy of a mechanical or electrical system is continuously dissipated, then the system must eventually settle down to an equilibrium point, no matter whether it is linear or nonlinear. Thus, we may conclude the stability of a system by studying the energy change rate of the system, without explicitly solving the original differential Eq. (1).
Let Bϵ be a ball of size ϵ around the origin x=0, i.e., Bϵ={xRS:x<ϵ}. A scalar continuous function V(x):RSR is locally positive definite (LPD) if V(0)=0 and x0V(x)>0 in the ball Bϵ. If V(0)=0 and the above property holds for the whole state space, then V(x) is globally positive definite (GPD). Similarly, a scalar continuous function V(x) is locally (or globally) positive semi-definite if V(0)=0 and x0V(x)0 in the ball Bϵ (or for the whole state space). A scalar continuous function V(x,t):RS×R>0R is locally positive definite if V(0,t)=0 and there exists a time-invariant LPD function V0(x) that is dominated by V(x,t), i.e., V(x,t)V0(x),tt0. If V(0,t)=0 and the above property holds for the whole state space, then V(x,t) is globally positive definite. A scalar continuous function V(x,t) is called decrescent if V(0,t)=0 and there exists a time-invariant LPD function V1(x) that dominates V(x,t), i.e., V(x,t)V1(x),tt0. A scalar continuous function V(x,t) is called radially unbounded if limxV(x,t) uniformly on t.
Let V(x,t) be a non-negative function with derivative V̇(x,t) along the trajectories of the system, i.e., V̇(x,t)=Vt+Vxf(x,t).Lyapunov’s theorems of equilibrium point stability are summarized in Table 2. Note that they are all sufficiency theorems. If for a particular choice of Lyapunov function candidate V, the condition on V̇ is not met, we cannot draw any conclusions on the system’s stability. Many Lyapunov functions may exist for the same system. Specific choices of Lyapunov functions may yield more precise results than others.
For an LTI system ẋ(t)=Mx(t), we often use a quadratic Lyapunov function candidate V=xPx, where P is a symmetric positive definite matrix (i.e., xPx>0,x0), denoted as P0. It is easy to derive that V̇=ẋPx+xPẋ=xQx,where (2)Q=Q=(MP+PM).The above equation is often called the Lyapunov equation of the linear time-invariant system. If Q is positive definite, then the system is globally asymptotically stable. However, if Q is not positive definite, then no stability conclusion can be drawn. Fortunately, Lyapunov proved a necessary and sufficient condition for a linear time-invariant system to be strictly stable. For any symmetric positive definite matrix Q, the unique matrix solution P of the Lyapunov Eq. (2) is symmetric positive definite. Therefore, we can start by choosing a simple positive definite matrix Q (e.g., the identity matrix I), then solve for P from the Lyapunov equation, and finally verify whether P is positive definite.

Table 2. Summary of Lyapunov stability theorems on equilibrium point x=0 [129].

V(x,t)V̇(x,t)Stability (in the sense of Lyapunov)
LPDLPSDstable
LPD, decrescentLPSDuniformly stable
LPD, decrescentLPDuniformly asymptotically stable
GPD, decrescent,
radially unbounded
GPDglobally uniformly asymptotically
stable

3. Linear stability

As a simple application of Lyapunov’s indirect method, linear stability analysis has been extensively applied to ecological systems, helping us better understand the intricate relationships between stability and biodiversity [22], [23], [33], [126], [132], [133], [134], [135], [136]. More significantly, under some special conditions, local stability determined by linear stability analysis implies global stability for certain ecological systems [137].
In his pioneer work, May considered an ecological system of S species coexisting at a feasible equilibrium point x of a dynamical system ẋ(t)=f(x(t)) that describes the time-dependent abundance vector x(t) of the S species [22], [23]. Suppose that one species is subjected to a small but sudden population increase or decrease. As a result, other species populations may show immediate changes away from the equilibrium. The manipulated species itself, if with initially increased (or reduced) abundance, may begin to decline (or increase) toward the equilibrium because of self-regulation or interspecific interactions. These immediate, direct changes can be referred to as first-order effects and described by the linearized equation of the original dynamical system, i.e., (3)ż(t)=Mz(t),where z(t)=x(t)xRS denotes the deviation from the equilibrium, and the Jacobian matrix M=f(x)x|x=xRS×S is often referred to as the community matrix. The diagonal elements of the community matrix Mii=fi(x)xi|x=x represent the self-regulation of species i, while the off-diagonal elements Mij=fi(x)xj|x=x capture the impact that species j has on species i around the equilibrium point x.
Numerous studies have explored the stability of the linearized Eq. (3), most of which employ random matrix theory to characterize the eigenvalue spectrum of the community matrix [22], [33], [75], [76], [77], [79]. These studies can be broadly classified into two categories. The first category, referred to as model-implicit approaches, focuses on analyzing the stability of the linearized equation without requiring specific knowledge about the underlying dynamics f(x). The second category, referred to as model-explicit approaches, incorporates additional information about f(x), e.g., the GLV model, to provide insights into the stability analysis.

3.1. Model-implicit approaches

Since the empirical parameterization of the exact functional form of f(x) is difficult for an ecological system, model-implicit approaches directly utilize the linearized representation to perform stability analysis.

3.1.1. May’s classical result

May considered that Mij are randomly drawn from a distribution with mean μ=0 and variance σ2 with probability C and are 0 otherwise. Hence, σ represents the characteristic interspecies interaction strength and C is the ratio between actual and potential interactions in the ecological system (often referred to as the connectance). For simplicity, the diagonal elements are chosen to be the same with d=1, representing the intrinsic damping timescale of each species, so that if disturbed from equilibrium, it would return with such a damping time by itself. May found that for random interactions drawn from a Gaussian distribution N(0,σ2), a randomly assembled system is stable (in the sense that all the eigenvalues of the community matrix have negative real parts) if the so-called “complexity” measure (4)σCS<d=1.This implies that more complexity (i.e., larger CS) tend to destabilize community dynamics [22], [23].
May’s classical result was inspired by Ginibre’s work on the circular law for random matrices with Gaussian distribution of entries [138] and Wigner’s work on the semi-circle law for random symmetric matrices [139]. In the 1980s, Girko developed a method to establish the circular law for more general distributions [140]. Hence, in random matrix theory, the circular law is also referred to as Girko’s Circular Law [141]. Consider a random S×S real matrix M with entries independent and taken randomly from a normal distribution N(0,σ2). Then as S, the eigenvalues of M/Sσ2 are uniformly distributed in the unit disk centered at (0,0) in the complex plane [140]. Sommers et al. considered possible correlations between the off-diagonal elements Mij and Mji and proved that the eigenvalues of M are uniformly distributed in an ellipse with the real and imaginary directions 1+ρ and 1ρ (where ρ=E[MijMji]/Var(Mij) for ij), respectively, known as Sommers’ Elliptical Law [142].

Table 3. Stability criteria for ecological systems with different interaction types. In each case, the stability criterion is derived for a large community matrix M. Diagonal elements Mii, representing self-regulation, are set to be d. For random interactions, off-diagonal elements Mij are randomly drawn from a normal distribution N(0,σ2) with probability C and are 0 otherwise. Predator–prey interactions come in pairs with opposite signs, i.e., Mij>0 and Mji<0. With probability C, we sample Mij from |N(0,σ2)| and Mji from |N(0,σ2)|. With probability (1C), both Mij and Mji are set to be 0. Community matrices with mutualistic and/or competitive interactions can be constructed similarly [33].

Interaction typeStability criterion
RandomσSC<d
Random with correlation (ρ)σSC(1+ρ)<d
Random with degree heterogeneity (ξ)σ(S1)Cξ<d
Predator–PreyσSC<dππ2
Mixture of mutualismσSC<dππ+2
and competition
Mutualismσ(S1)C2π<d
CompetitionσSCπ+24Cπ(π2C)+C2π<d
May’s result continues to be influential almost five decades later, not because it asserts that ecological systems must be inherently unstable, but rather because it highlights the importance of specific structural characteristics that enable real ecological systems to maintain stability despite their inherent complexity [143]. In other words, nature must adopt some devious and delicate strategies to cope with this stability-complexity paradox.

3.1.2. Impacts of interaction types and correlation

One of the specificity is the existence of well-defined interspecific relationships observed in nature, e.g., predator–prey, competition, mutualism, and a mixture of competition and mutualism. In 2012, by leveraging Sommers’ Elliptical Law [142], Allesina et al. refined May’s result and provided stability criteria for all these interspecific interaction types [33], as shown in Table 3 and Fig. 3a, b, and c. They found remarkable differences between predator–prey interactions, which are stabilizing, but mutualistic and competitive interactions, which are destabilizing. Additionally, the correlation between interaction strengths constrains the relative proportion of interaction types. Tang et al. incorporated correlation into the stability analysis by deriving a new stability criterion for large ecological systems: (5)max{σS(1+ρ),Sμ}μ<d,where the connectance of the network is implicitly involved in σ, ρ represents the overall correlation between pairs of interactions, and μ denotes the mean of the off-diagonal elements in the community matrix M [86]. In particular, if μ0, which is typically the case in food webs (with predator–prey interactions) and competitive networks, the stability criterion (5) reduces to σS(1+ρ)μ<d. Therefore, the authors concluded that the effect of the correlation of interaction strengths substantially influences the stability of large food webs compared to other network structural properties. Other studies also demonstrated that the presence of correlation between interactions of species can significantly influence the locations of the eigenvalues of the community matrix and the resulting stability of ecological systems [33], [144], [145].
  1. Download: Download high-res image (626KB)
  2. Download: Download full-size image

Fig. 3. Eigenvalue distributions of 10 community matrices (color) with d=1 on the diagonal and off-diagonal elements following the random (a), predator–prey (b), and mixture of competition and mutualism interactions (c), respectively. S=250, C=0.25, σ=1. The black ellipses are analytical results. Impact of degree heterogeneity on the stability of ecological systems with random interactions (d), predator–prey (e), and mixture of competition and mutualism interactions (f). The dots represent the results from numerical simulations on 3-modal networks with different connectances. Higher Re[λ1] indicates lower stability. Each error bar represents the standard deviation of 100 independent runs. S=1200, σ=1, and d=0. Panel d was drawn in the log–log scale.

Panels a, b, and c were redrawn from [33]. Panels d, e, and f were redrawn from [34].

3.1.3. Impacts of degree heterogeneity

Many of the “devious strategies” adopted by nature can now be tested with the revised formula as a reference point. For example, one can further study the impact of degree heterogeneity on the stability of ecological systems. Degree heterogeneity measures the variability of the number of interactions associated with each species. Yan et al. found that for ecological systems with random interactions or a mixture of competition and mutualism interactions, increasing the degree of heterogeneity always destabilizes ecological systems [34], see Fig. 3d, f. For ecological systems with predator–prey interactions constructed from either simple network models or a realistic food web model (cascade model), high heterogeneity is always destabilizing, yet moderate heterogeneity is stabilizing, see Fig. 3e. Furthermore, they obtained a stability criterion for large ecological systems with random interactions by considering the factor of degree heterogeneity, i.e., (6)σ(S1)Cξ<d,where ξ denotes the degree heterogeneity [34]. Similarly, Allesina et al. found that broad degree distributions tend to stabilize food webs by approximating the real part of the leading (“rightmost”) eigenvalue of the community matrix [75]. Recently, Baron derived a closed-form expression for the eigenvalue spectrum of a general directed and weighted network [146]. The findings of this study suggest that network heterogeneity appears to be a destabilizing influence in most circumstances, except when the interactions are very asymmetric (e.g., very asymmetric predator–prey interactions). These results are consistent with what were reported by Yan et al. [34].

3.1.4. Impacts of self-regulation

Self-regulation (reflected as negative diagonal elements of the community matrix) is also a key factor for stability in real or random ecological systems. For random ecological systems, Barabás et al. utilized the quaternionic resolvent method [147] to determine the minimum fraction of self-regulating species required for local stability [35]. They identified two key factors: the interaction correlation ρ, and the scaled strength of self-regulation d/(σS). Remarkably, their results show that local stability cannot be achieved unless either a significant portion of species strongly regulate themselves or ρ approaches 1. The theoretical findings for random ecological systems remain applicable to empirical food webs, indicating that stability is attainable only when the majority of species demonstrate substantially strong self-regulation. In addition, Tang et al. investigated the effect of non-constant diagonal on the eigenvalue distribution of the community matrix M [86]. They found that a moderate variance of the diagonal elements minimally affects the distribution of the eigenvalues. Nevertheless, when the variance significantly surpasses that of the interspecific interactions, the impact on stability depends on the specific patterns of the diagonal elements. When the self-regulation is stronger for species with fewer interactions, the impact of a considerable variance on stability remains negligible.

3.1.5. Impacts of modularity

The stability of ecological systems can also be influenced by network modularity. A modular network can be divided into different modules or subsystems, and the interactions within each subsystem are much more frequent than those between subsystems [148]. The modularity of a network can be defined as Q=LwE[Lw]Lw+Lb,where Lw is the observed number of interactions within the subsystems, and Lb is the number of inter-subsystem interactions. Grilli et al. studied the effect of modularity Q on the stability of ecological systems consisting of two subsystems or modules [36]. They found that modularity exhibits a moderate stabilizing effect when the subsystems have similar sizes and the overall mean interaction strength is negative. In particular, the stabilizing effect becomes stronger when there are negative correlations between interaction strengths. Conversely, anti-modularity is highly destabilizing, except for the case where the overall mean interaction strength is close to zero. This near-zero mean value of interaction strengths corresponds to the parameter regime yielding stability in many foundational diversity-stability studies [22], [33], [47]. The authors further investigated the effect of modularity in food webs through numerical simulations and found that the results remain qualitatively unchanged.

3.1.6. Impacts of dispersal

Spatial flows (e.g., exchanges of individuals, energy, and material) among local ecosystems are ubiquitous in nature [149]. Gravel et al. incorporated dispersal, i.e., spatial movement of species among local ecological systems, into the stability analysis of meta-ecological systems [80]. The community matrix M of a meta-ecological system can be expressed as the sum of three matrices, i.e., M=D+Q+K,where D is a diagonal matrix that accounts for intraspecific density dependence, Q is a matrix representing dispersal among patches with diffusion coefficient q, and K is a block diagonal matrix that contains the local community matrices (random matrices with S species, connectance C, and interspecific interaction strength σ). By assuming the number of local systems and S are large, the authors obtained the following stability criteria: (7)σC(S1)/Se<dfor sufficiently large q(q),σC(S1)<d+qfor small q,where Se=S/[1+(S1)ρ]. Therefore, the effect of dispersal can promote stability in meta-ecological systems for both large and small q. Dispersal shifts most of the eigenvalues of the community matrix towards more negative values and reduces the range of the remaining eigenvalues in proportion to the number of effective patches for large q, while it acts as a negative intraspecific feedback, which is known to be stabilizing, for small q. Additionally, Baron et al. discovered that the introduction of dispersal can lead to Turing-type instability, where the equilibrium point becomes unstable with respect to spatial disturbances of wavelengths in the abundances of species for an ecological system with a trophic structure [76]. While the inclusion of trophic structures often enhances stability in large ecological systems (e.g., food webs) [150], [151], it can lead to a peak in the real part of the maximum eigenvalue of the community matrix in the presence of dispersal, rendering the equilibrium point unstable.

3.1.7. Impacts of time delay

Empirical evidence has indicated that species interactions often exhibit time lags, rather than occurring instantaneously [152], [153]. Thus, time delays can hold significant implications for stability and coexistence. Notably, Pigani et al. investigated delay effects on the stability of large ecological systems [83], which can be captured by the following linearized equation: (8)ż(t)=Mz(t)+Mdelayz(tτ),where Mdelay is the community matrix with delay. For simplicity, the community matrix M is set to dI where d0, which guarantees stability for a sufficiently small community matrix delay. In other words, the current intraspecific interactions are always stabilizing. Mdelay is supposed to be a random matrix with a constant diagonal entry ddelay0 and off-diagonal elements normally distributed with zero mean, standard derivation σdelay, and connectance Cdelay. The authors proved that (8) is stable if all the roots of the characteristic equation, defined as H(z)=zλλdelayezτ where λ and λdelay are the eigenvalues of M and Mdelay, respectively, have negative real parts. Importantly, these roots can be solved implicitly as a function of τ, along with the parameters d, ddelay, σdelay, and Cdelay. They found that an increasing delay tends to destabilize the system, and if a system is already unstable for τ=0, the delay cannot stabilize the system. Furthermore, the authors determined the critical delay τ as the minimum value of τ above which the system becomes unstable. Finally, distributed delay was considered, i.e., the second term in (8) is modified to Mdelay0exp{ττ̃}/τ̃z(tτ)dτ, which can be derived from logistic or resource competition models [154]. They found that the system becomes more and more unstable as τ̃ increases. However, if τ̃ is large enough, the system eventually goes back to the stable regime.

3.1.8. Limitations of model-implicit approaches

May’s approach and the follow-up studies offer a valuable theoretical framework for understanding the stability of ecological networks with various structures. Yet, these model-implicit approaches heavily rely on randomly sampling interactions without an underlying dynamical model, which results in a lack of biologically realistic representation of species interactions. Consequently, the analytical results produced by these approaches are independent from the underlying model from which they are hypothetically derived. While informative, these results may not accurately capture the nuanced and complex interactions that characterize the stability of real ecological systems. One way to enhance this framework is to integrate underlying ecological models into the linear stability analysis, as detailed in the next subsection.

3.2. Model-explicit approaches

The most commonly used model to describe the dynamics of ecological systems is the generalized Lotka–Volterra (GLV) model [37], [38] defined as (9)ẋi(t)=xi(t)[ri+j=1SAijxj(t)]for i=1,2,,S, where ri is the intrinsic growth rate of species i and ARS×S is the interaction matrix whose off-diagonal elements Aij represent the effect that species j has upon species i. By assuming the existence of a feasible equilibrium point x (i.e., xi>0 for all species), the community matrix can be computed as M=XA, where XRS×S the diagonal equilibrium matrix such that Xii=xi. Hence, the community matrix can be regarded as the scaled interaction matrix. Similarly to model-implicit approaches, we can investigate the stability properties of the GLV model based the linearized Eq. (3).
  1. Download: Download high-res image (449KB)
  2. Download: Download full-size image

Fig. 4. Eigenvalue distributions of the three matrices M, Q, and J. The off-diagonal elements of A are sampled from a normal bivariate distribution with identical marginal μA=5/S, σA=5/S, and correlation ρA=0.5. The diagonal elements of A are fixed to μd=1, while the diagonal elements of X are sampled from a uniform distribution on [0,1]. The two matrices Q and J are defined as Q=X[(μdμA)I+μA1] and J=X(μdI+B), respectively, such that the bulk of the eigenvalue distributions of M are the same as these of J, and the outlier eigenvalue of M is the same as that of Q. Here, B is a random matrix with zero diagonal elements whose off-diagonal elements have mean zero and variance σA2. S=500.

This figure was redrawn from [79].

3.2.1. Impacts of equilibria

Equilibrium points play a significant role in the local stability analysis of the GLV model, as its community matrix is directly related to the equilibrium points of the system. Gibbs et al. explored the eigenvalue spectrum of the community matrix with effect of population abundances [79]. Assume that the diagonal elements of the equilibrium matrix X are drawn from an arbitrary distribution with positive support, the diagonal elements of the interaction matrix A are drawn from an arbitrary distribution with support in the negative axis, and each off-diagonal pair of A is drawn from independently from a bivariate distribution. Denote μX, μd, and μA as the mean of the diagonal of X, the diagonal of A, and the off-diagonal of A, respectively. The authors proved that if μA0, the eigenvalue spectrum of the community matrix M consists of a bulk of eigenvalues and an outlier. The mean of the bulk eigenvalues is μX(μdμA) determined by the eigenvalues of the matrix J=X(μdI+B), and the outlier is computed as λoutlier=μX[μd+(S1)μA],determined by the largest eigenvalue of the matrix Q=X[(μdμA)I+μA1] (1 is the all-one matrix), see Fig. 4. Based on their finding, the authors concluded that for mutualistic ecological systems (μA>0), if the interaction matrix is stable, then the community matrix will also be stable. The authors also showed that this property holds in the case where μA=0 and ρA=0 (correlation of the off-diagonal elements of A) using the cavity method. It is important to note that Gibbs et al. made the assumption that the distribution from which the equilibrium is drawn is independent of the interaction matrix, which does not typically hold in reality. Nevertheless, Liu et al. discovered that if the equilibrium point follows a specific distribution related to the elements of the interaction matrix A, Gibbs et al.’s assumption remains valid [155]. Around the same time, Stone demonstrated that for large ecological systems described by the GLV model, the stability of the interaction matrix implies the stability of the community matrix, under the condition that all species have positive equilibria [57]. Note that for small ecological systems, this result is not necessarily true [156].

3.2.2. Impacts of extinction

Extinction is generic in the GLV model with random interactions [77], [157], [158]. Recent studies have highlighted that an equilibrium point of the GLV model is locally stable if and only if all the eigenvalues of the reduced interaction matrix (i.e., the interaction matrix between the species in the surviving sub-community) have negative real parts [57], [77], [159], [160]. Baron et al. focused on the eigenvalue spectrum of the reduced interaction matrix in the GLV model, where the spectrum of the reduced interaction matrix also consists of a bulk set of eigenvalues and an outlier [77]. The authors demonstrated that the universality principle [161] (i.e., the eigenvalue spectrum of a random matrix often only depends on the first and second moments of its elements) holds for the bulk region following Ellipse Law but not for the outlier eigenvalue. The outlier eigenvalue can be solved from the generating functional approach. Prediction of extinction boundary and new equilibrium after a primary extinction were also discussed in [157], [158], respectively. Furthermore, the stability of ecological systems with extinction is closely related to the concept of sector stability [115], see Section 8.
  1. Download: Download high-res image (141KB)
  2. Download: Download full-size image

Fig. 5. Eigenvalue distributions with time delay, where the blue teardrop-shaped area represents the stability region with time delay. (a) The eigenvalues (represented by red dots) are within the teardrop-shaped area, indicating that the corresponding GLV model with time delay is locally stable. (b) The eigenvalues are outside the teardrop-shaped area, indicating that the corresponding GLV model with time delay is unstable. τ=1.

The figure was redrawn from [88].

3.2.3. Impacts of time delay

Similar to Pigani et al.’s work as discussed in Section 3.1.7, Yang et al. investigated the stability of the time-delayed GLV model [88], which is defined as (10)ẋi(t)=xi(t)[ri+j=1SAijxj(tτ)].The linearized equation of the time-delayed GLV model is given by (11)ż(t)=Mdelayz(tτ),where Mdelay=XA is the community matrix with delay. For simplicity, X is set to I such that each species has unit abundance. However, the validation of this assumption remains an open problem. They proved that (11) is stable if all the roots of the characteristic equation, defined as H(z)=zλdelayezτ where λdelay are the eigenvalues of Mdelay, have negative real parts. This equation implies that all the eigenvalues of Mdelay are required to be located in a teardrop-shaped region to ensure stability, see Fig. 5. The boundary of the teardrop-shaped region determined by time delay can be computed as τ=tan1(x/y)/x2+y2 where x and y are the real and imaginary part of eigenvalues, respectively. Since the distribution of the eigenvalues of Mdelay (i.e., the interaction matrix A is this case) with random, mutualistic, or competitive interactions is well-established, it becomes straightforward to determine stability. Based on the theoretical findings, the authors found that time delay plays a critical role in ecological systems, where large delay is often destabilizing, but short delay can considerably enhance community stability.
Similar conclusions were also obtained by Saeedian et al. when investigating the impact of time delay on the emergent stability patterns with the time-delayed GLV model (10) [162]. They further determined the existence of a Hopf bifurcation (i.e., the two complex conjugate eigenvalues of the community matrix, with non-zero imaginary part, simultaneously cross the imaginary axis into the right half-plane) at τc, which is computed as τc=minλiλ(Mdelay)1λiartan|Re[λi]Im[λi]|for the time-delayed GLV model with a large number of species. By increasing τ above τc, the amplitude of the trajectories also increases. However, when τ significantly exceeds τc, the amplitude of the oscillations becomes so large that it leads to numerical instabilities, causing the population trajectories to numerically diverge.

3.2.4. Impacts of stochastic noises

Real ecological systems are inherently stochastic with constant external perturbations and internal fluctuations [163]. Krumbeck et al. developed a theoretical framework for the analysis of temporal stability of ecological systems with stochastic noises [82]. They provided analytical predictions for the power spectral density of the stochastic GLV model defined as (12)ẋi(t)=xi(t)[ri+j=1SAijxj(t)]+1Wηi(t),where W is the area of the species living domain, and ηi(t) are Gaussian noises. The power spectral density is a statistical measure that can capture various aspects of temporal stability (e.g., the height of the spectrum gives information about the magnitude of stochastic fluctuations; the locations of nonzero peaks correspond to quasi-cyclic signals; a peak at zero indicates baseline wander). Specifically, the authors utilized the linearized equation of the stochastic GLV model (12), i.e., (13)ż(t)=Mz(t)+ζ(t),where ζ(t) is a vector of Gaussian white noises with correlation matrix F. The power spectral density of fluctuations Ψ(ω) in the frequency domain can be computed as Ψ(ω)=(MiωI)1F(M+iωI)1,where i here is the imaginary number. By analyzing the mean-field power spectral density, i.e., E[Φii], they showed that different network structures (random, mutualistic, and competitive) have unique signatures in the spectrum of fluctuations. These fluctuations are characterized by a few key parameters: the mean, variance, and correlation of the entries in the community matrix, as well as the noise correlator. They further investigated the effect of trophic structures and identified a gap in the power spectral density, which indicates that high-level trophic structures contribute to enhanced long-term temporal stability.

3.2.5. Impacts of evolved system size

Previous work on the linear stability of the GLV model is concerned with ecological systems of predetermined and fixed sizes. Galla explored the stability of ecological systems with evolved system size [78], using methods from statistical mechanics and the theory of disordered systems [164], [165]. Specifically, he exploited generating functionals to derive the effective dynamics for the GLV model (with ri=1) computed as (14)ẋ(t)=x(t)[1x(t)+μAM(t)+ρAσA20tG(t,t)x(t)dt+η(t)],where M(t) is the average species concentration, G(t,t) is the response function, and η(t) is the Gaussian noise. As described previously, μA, σA and ρA are the mean, standard derivation, and correlation of the off-diagonal elements of the interaction matrix A, respectively. The effective process characterizes the dynamics of a single representative species, denoted by x(t), and captures the statistical behaviors of the ecological system. The linearized effective dynamics can be computed as (15)ż(t)=x[z(t)+ρAσA20tG(t,t)z(t)dt+v(t)+ζ(t)],where z(t) denotes fluctuations about the equilibrium x, v(t) is the deviation of the noise in the effective process, and ζ(t) is the Gaussian white noise of unit amplitude. By performing Fourier transform with a focus on the long-time behavior of perturbations (ω=0), the author established that the GLV model has stable equilibrium points in the limit of large population for (16)σA<2(1+ρA)2.This result implies that predator–prey relationships enhance stability, while variability in species interactions promotes instability, which aligns with prior findings as reported in [34], [86].
Poley et al. applied Galla’s approach to analyze the stability of the GLV model with hierarchical interactions (i.e., species are ordered in a natural cascade or hierarchy such that any species can prey on only those below it and can be preyed on by those above it in the hierarchy, also known as the cascade model) [84]. The model is defined as (17)ẋia(t)=xia(t)[ria+b=1Bj=1SbAijabxjb(t)]for i=1,2,,Sa, where B is the total number of local systems, and a,b=1,2,,B denote the indices of the local systems with the associated system size Sa and Sb. Their findings indicate that a strong hierarchical structure promotes stability but also leads to increased species extinctions and lower abundances in the surviving community. Similarly, Sidhom and Galla adapted Galla’s approach to study the GLV model with nonlinear feedback defined as (18)ẋi(t)=xi(t)[ri+g(j=1SAijxj(t))]for i=1,2,,S, where the function g(u)=2aua+2|u| denotes the nonlinear feedback (a is the saturation parameter) [85]. This form of feedback was initially introduced to describe the predator’s growth rate when interacting with prey. It is plausible that the benefits from extra prey will eventually saturate as prey numbers become large. The authors concluded that the stability and diversity of ecological systems improves with the introduction of nonlinear feedback.
  1. Download: Download high-res image (300KB)
  2. Download: Download full-size image

Fig. 6. Stability of real and permuted food webs in relation to complexity. These permuted food webs were constructed from real food webs by removing some non-random features, namely the row structure, topology, pairwise correlation, and interaction strength. The more the “rightmost” eigenvalues, i.e., Re[λ1], close to zero, the more stable of the food webs (because the diagonal elements of the community matrix are set to zero).

This figure was redrawn from [81].

3.2.6. Applications to real ecological systems

Last but not least, while the stability analysis of random ecological systems offers powerful insights for the study of ecological dynamics and the understanding of how species interactions shape the stability and diversity of ecological communities, empirically applying the theory to real ecological systems poses a formidable challenge [32], [81], [166], [167], [168], [169], [170]. Notably, Jacquet et al. performed the local stability analysis of 116 real food webs sampled worldwide from marine, freshwater, and terrestrial habitats using the GLV model [81]. The community matrix M was constructed by multiplying the interaction matrix A with species biomass for each food web. The interaction coefficients of A can be translated from the parameters of the corresponding Ecopath model, a widely-used ecological modeling software that creates mass-balanced models of ecological systems, providing a way to quantify energy flows and the trophic interactions between different species [171]. The authors then measured the stability of food webs using the real part of the dominant eigenvalue of the community matrix M. Using randomization tests, they found that negative correlation between interaction strengths with high frequency of weak interactions is a strong stabilizing property in real food webs, see Fig. 6. Their findings reveal that empirical food webs exhibit some non-random characteristics that lead to the absence of a complexity–stability relationship.

4. Sign stability

The notion of sign stability is of particular interest for ecology, economics, chemistry, and engineering [93], [95], [172], [173], [174], [175], [176], [177]. In an ecological system, the community matrix M (or the interaction matrix A) associated with certain models (e.g., the GLV model) might be only known qualitatively in the sense that the signs of the elements Mij can be determined with reasonable confidence, but the actual magnitudes may be very difficult to determine, see Fig. 7. Such matrices are often referred to as sign matrices. A sign matrix M is called sign stable (or sign semi-stable) if each of its eigenvalues has negative real part (or non-positive real part, respectively) for all numerical matrices of the same sign pattern [92], [94], [95], [96], [174].

4.1. Characterizations of sign stable matrices

A sign matrix MRS×S is sign semi-stable if and only if the following conditions are all satisfied: (i) Mii0 for i=1,2,,S; (ii) MijMji0 for i,j=1,2,,S and ij; (iii) there exists no simple cycle of length 3 in the digraph generated by M [94]. It is also known that if Mii<0 for i=1,2,,S, then (ii) and (iii) are necessary and sufficient for M to be sign stable. Quirk and Ruppert further claimed that (i)–(iii) as well as (iv) Mii<0 for some i and (v) there exists a non-zero term in the expansion of det(M) are both necessary and sufficient for M to be sign stable [95]. However, this was proved to be incorrect later by Jeffries et al. [92]. Later on, Yamada demonstrated that such an exceptional case is very rare and proved that the conditions (i)–(v) proposed by Quirk and Ruppert is actually necessary and sufficient for a system to be generically sign stable, i.e., sign stable for almost all parameter values except for some pathological cases with measure zero [96].
Quirk and Ruppert’s conditions (i)–(v) can be translated into ecological terms to provide a comprehensive characterization of sign stable patterns (although it may not be ecologically meaningful) [93]. According to (i) and (iv), a sign stable ecological system should not include self-promoting species and must have at least one species with self-regulation. Condition (ii) necessitates the absence of both competitive and mutualistic relations. Condition (iii) states that there are no closed directed loops (i.e., simple cycles) with more than two links in the community structure. For instance, in the rock–paper–scissors dynamics with three species i, j and k, each species has a direct competitive interaction with another species, creating a cycle: ijki. This forms a simple cycle of length 3 in the directed graph of interactions. Finally, condition (v) requires that the digraph generated by M must include a specific number of non-overlapping simple cycles that encompass all the nodes. This stems from the fact that any permutation can be represented as the composition of a set number of cyclic permutations, each acting within non-overlapping subsets of the set {1,2,,S}.
  1. Download: Download high-res image (160KB)
  2. Download: Download full-size image

Fig. 7. Examples of signed ecological networks and their corresponding sign community matrices. (a) Self-regulated species 1 forms a commensal relationship with each species within the prey–predator pair of 2 and 3. (b) Linear trophic chain in which each successive species is preyed upon the subsequent one and species 3 is self-regulated.

This figure was redrawn from [93].
Jeffries et al. introduced two additional conditions, distinct from (iv) and (v), called the coloring and matching conditions [92]. The two conditions, combined with (i)–(iii), are necessary and sufficient for sign stability. Let RM={i|Mii0}. An RM-coloring of an undirected graph is a partition of its nodes into two sets, black and white, such that each node in RM is black (one of which can be empty), no black node has exactly one white neighbor, and each white node has at least one white neighbor. The coloring condition then states that in every RM-coloring of the undirected graph generated by M, all nodes are black. Additionally, a (VRM)-complete matching in an undirected graph is a set M of disjoint edges such that an exact cover of the node set V can be obtained using the pairs in M and certain singletons from RM. The matching condition then asserts that the undirected graph generated by M admits a (VRM)-complete matching. However, the ecological systems characterized by Quirk and Ruppert or Jeffries et al. are not commonly observed in nature.

4.2. Applications to ecological systems

The notion of sign stability has been applied to various ecological systems through different approaches. Dambacher et al. introduced two qualitative metrics — weighted feedback and weighted determinants based on the Hurwitz criterion [90], which can be recast into two conditions: (i) the coefficients of the characteristic equation of M must have the same sign; (ii) the corresponding Hurwitz determinants must all be positive. The two metrics offer a practical mean to identify the relative degree to which stable parameter space can be constrained based on system structure and complexity. Remarkably, Haraldsson et al. utilized the weighted feedback and weighted determinants to investigate the sign stability of social-ecological systems, which is an important tool to understand human-nature relations [91].
Moreover, Allesina and Pascual studied the sign stability of random and real food webs through a fairly empirical approach, which does not strictly adhere to the exact definition of sign stability [89]. Given a community matrix M (either randomly generated or empirical), first determine the stability based on its eigenvalues. If it is stable (in terms of Re[λ(M)]0), generate 100 matrices that have the same sign pattern with random magnitude. Lastly, measure the percentage of the random generated matrices that are stable. The percentage can tell whether the stability is due to a particular combination of coefficients or the sign pattern of the network itself. Based on this approach, the authors demonstrated that predator–prey interactions can promote stability, highly robust to perturbations of interaction strength, in real ecological systems.

5. Diagonal stability

The notion of diagonal stability was first introduced by Volterra in the 1930s [178]. It has been particularly useful for the stability analysis of ecological systems and other networked systems [51], [179], [180], [181], [182], [183], [184], [185]. A matrix A is called diagonally stable if there exists a diagonal matrix P0 (i.e., P is positive definite) that renders AP+PA=Q0,where the notion 0 means a negative definite matrix. We use the notation D(P) for the class of diagonally stable matrices. The positive diagonal matrix P is often called Volterra multiplier in literature [104]. In many cases, the necessary and sufficient conditions for the Lyapunov stability of nonlinear systems are also the necessary and sufficient conditions for the diagonal stability of a certain matrix associated to the nonlinear system. This matrix naturally captures the underlying network structure of the nonlinear dynamical system.
Diagonal stability has been successfully applied to various types of dynamical systems [186], [187], [188], [189], [190], [191]. It has nice “structural consequences”. For example, the principal sub-matrices of a diagonally stable matrix are also diagonally stable, implying that all the corresponding “principal sub-systems” of a given diagonally stable system are diagonally stable [51]. In this section, we begin by introducing the theoretical framework for determining the Lyapunov stability of the GLV model via diagonal stability analysis. Subsequently, we delve into the characterizations of diagonal stability of the interaction matrix that is associated with special interconnection or network structures, which offer an effective mean of determining the stability of ecological systems.

5.1. Lyapunov stability of the GLV model

Before talking about the stability of the GLV model, we first introduce the notion of Persidskii-type systems. Persidskii-type systems are typical examples that admit diagonal-type Lyapunov functions [192], [193], [194], [195]. We need to introduce a few concepts to define Persidskii-type systems. A function f:RSRS:xf(x) is called diagonal if, the ith component of f, i.e., fi, is a function of xi alone. A function f:RR is said to be in sector [a,b] if xR, f(x) lies between ax and bx, i.e., (f(x)bx)(f(x)ax)0.For example, sector [1,1] means |f(x)||x|. Sector [0,] means f(x) and x always have same sign. The class of infinite sector nonlinear functions S is defined to be functions in sector [0,] that satisfy 0xf(τ)dτ as |x|. Typical examples of infinite sector nonlinear functions are f(x)=x, f(x)=x3, and f(x)=tanh(x) [51]. A dynamic system ẋ(t)=f(x(t)) with x(t)RS and f:RSRS is said to be of Persidskii-type if it has the following form: (19)ẋi(t)=j=1SAijfj(xj(t))for i=1,2,,S, where fjS for all j=1,2,,S. In other words, f is diagonal and fj is in the class of infinite sector nonlinear functions. Note that (19) can be used to describe a wide range of complex networked systems, where Aij capture the weighted wiring diagram.
For Persidskii-type systems, we can introduce a diagonal-type Lyapunov function of the following form: V(x)=12i=1Spi0xifi(τ)dτ,where pi denotes the ith diagonal of P. The equilibrium point x=0 of the Persidskii-type system is globally asymptotically stable (in the sense of Lyapunov) if AD(P). This can be seen by computing V̇(x) along the trajectory of (19), yielding V̇(x)=f(x)(AP+PA)f(x).Since f is diagonal and fjS, if AD(P), therefore V̇(x) is negative definite. Moreover, the functions fiS ensure the radial unboundedness of V(x). Hence, according to Lyapunov’s theorems of stability, x=0 is globally asymptotically stable (in the sense of Lyapunov).
By assuming the existence of a non-trivial equilibrium point x (i.e., xi>0 for all species) and defining zi(t)=logxi(t)xiandgj(zj(t))=xj(ezj(t)1),we can bring the GLV model (9) into the form of Persidskii-type dynamics: żi(t)=j=1SAijgj(zj(t)),which admits the following diagonal-type Lyapunov function: V(z)=i=1Spi0zigi(τ)dτorV(x)=i=1Spixixixilogxixi.Using the above diagonal-type Lyapunov function, Goh first proved that for the GLV model (9), if the interaction matrix A is diagonally stable, then the non-trivial equilibrium point x in the positive orthant is globally asymptotically stable (in the sense of Lyapunov) [48]. Therefore, diagonal stability allows the GLV model to lie in the unique fixed-point phase, as depicted by Bunin [47]. Here, we provide an example of a 3-by-3 diagonally stable matrix and the corresponding vector field plot of the GLV model, see Fig. 8. Clearly, the non-trivial equilibrium point in the positive orthant is globally asymptotically stable.
The diagonal stability analysis has been applied to analyze the stability of various GLV-derived models. Following Goh’s findings, Wörz-Busekros offered a sufficient condition for the global stability of the GLV model with continuous-time delay via diagonal stability analysis [107]. Later on, Beretta and Takeuchi considered the global asymptotic stability of diffusion models with multiple species heterogeneous patches, in which each patch is governed by the GLV model with continuous-time delay [46]. Concurrently, Beretta and Solimano generalized the Wörz-Busekros’s outcome by considering a non-negative linear vector function of the species [45]. In addition, Kon exploited the diagonal stability theory to determine the stability of the GLV model with an age structure [49]. The notion of diagonal stability was further applied to more general ecological models such as Kolmogorov systems [102], where the GLV model is encompassed as a specific instance. Significantly, the stability properties of many quasi-polynomial dynamical systems, often used to represent many biochemical processes, can be studied through an equivalent GLV model that has a much simpler form [101], [103]. This is based on the fact that the former can be transformed into the latter with some appropriate changes of variables.
  1. Download: Download high-res image (456KB)
  2. Download: Download full-size image

Fig. 8. Example of a diagonally stable matrix. (a) Diagonally stable interaction matrix A. It is easy to check that the eigenvalues of (A+A)/2 are all negative. (b) Corresponding ecological network associated with A. (c) Vector field plot of the corresponding GLV model with the non-trivial equilibrium point indicated by the red star.

5.2. Characterizations of diagonally stable matrices

A general characterization of the diagonally stable interaction matrix in the GLV model remains elusive for more than three species [113], [196], though there exist efficient optimization-based algorithms to numerically check if a given matrix is diagonally stable, e.g., polynomial-time interior point algorithms [197]. In cases where the dimension is three or less, the diagonal stability of A can be determined by examining the signs of its principal minors [196]. For high-dimensional matrices under very special structural assumptions, one can derive necessary and sufficient conditions for diagonal stability. For example, if A is Metzler (i.e., Aij0 for ij), then A is diagonally stable if and only if all principal minors of A are positive [198]. More special examples are discussed in [104]. Moreover, there are approaches which can reduce the problem of determining whether ARS×S is diagonally stable into two simultaneous problems of (S1)×(S1) matrices, but the method becomes intractable for large S.
Grilli et al. imposed a condition of negative definiteness on A (i.e., all the eigenvalues of A+A are negative) to ensure the diagonal stability of A for studying the feasibility and coexistence of large random ecological systems [55]. It is known that a negative definite matrix is also diagonally stable, and the condition is much easier to verify and characterized for random matrices. Later on, Gibbs et al. estimated the “rightmost” eigenvalue of (A+A)/2 under the condition μA=0 (note that a negative μA only produces an equal shift in the rightmost eigenvalue) and proved that if (20)μd+2SσA2(1+ρA)<0,then A is diagonally stable, where μd is the mean of the diagonal elements of A, and σA2 and ρA are the variance and correlation of the off-diagonal elements of A [79].
Recently, necessary and sufficient diagonal stability conditions for matrices associated with special interconnection or network structures were studied [97], [98], [99], [100], [105], [106]. If an ecological system described by the GLV model exhibits these network structures, it will be effective and efficient to determine its Lyapunov stability through the diagonal stability of the interaction matrix A. We review these network structures as follows.
  1. Download: Download high-res image (132KB)
  2. Download: Download full-size image

Fig. 9. Negative feedback cyclic structure with the corresponding interaction matrix A, where αi>0, βi>0 for i=1,2,,S.

5.2.1. Negative feedback cyclic structure

In the negative feedback cyclic structure, the intermediate species play a facilitating role for the subsequent species, while the final species exerts inhibitory effects on the initial species, see Fig. 9. It has been shown that the interaction matrix A is Hurwitz, i.e., Re[λ(A)]<0, if it satisfies the so-called secant criterion [199], [200], [201], i.e., (21)β1β2βSα1α2αS<secπSS.When αi are equal, the secant criterion (21) is also necessary for A to be Hurwitz. Surprisingly, this secant criterion derived for linear stability is also a necessary and sufficient condition for diagonal stability of the corresponding class of matrices [99].
Note that a simple necessary condition for A to be diagonal stable is that all the diagonal elements Aii be negative. Since scaling the rows of A by positive constants does not change its diagonal stability, one can safely assume Aii=1. Moreover, a reducible matrix A can always be transformed into an upper-triangle form with a suitable permutation, and permutation does not change diagonal stability. Considering these two points, the secant criterion can be further generalized as follows [97]. Any matrix A that can be transformed via a suitable permutation P to the form of PAP=100β̃Sβ̃1100β̃21000β̃S11is diagonally stable if and only if |γ|Φ(sgn(γ),S)<1 where γ=β̃1β̃2β̃S0 is the cycle gain and Φ(sgn(γ),S)=cosS(π/S) for γ<0; 1 for γ>0.

5.2.2. Cactus structure

Arcak further generalized the secant criterion to multiple cycles when the weighted digraph of the interaction matrix A, denoted as G(A), possesses a “cactus” structure, i.e., any pair of distinct simple cycles have at most one common node [97]. Apparently, the negative feedback cyclic structure corresponds to a single cycle and is just a special case of the cactus structure. The digraph G(A) is defined to represent the off-diagonal entries of A that has S nodes and there is a directed edge (ji) with weight Aij if and only if Aij0. Self-loops corresponding to the diagonal entries Aii are excluded from G(A).
Arcak made two assumptions about A without loss of generality: (i) Aii=1 for i=1,,S; (ii) G(A) is strongly connected, or equivalently, A is irreducible. Arcak then defined an undirected graph for describing which cycles of the digraph G(A) intersect. Subsequently, he constructed a spanning tree in the undirected graph and simultaneously assigned directions to its edges to create an arborescence, i.e., a digraph in which every node can be reached from the root by one and only one path, see Fig. 10a, b, c. Using the hierarchy of the cycles of G(A) established by this arborescence presentation, he sequentially generated inequalities for the gains of each simple cycle.
Consider there are J simple cycles in G(A). The length of cycle j is denoted as nj for j=1,2,,J. Define Ij={ij(1),ij(2),,ij(nj)} and Ji={1jJ|iIj} as the set of nodes traversed by cycle j and the set of cycles that node i belongs to, respectively. Denote γj as the gain for cycle j. Then the stability condition can be summarized as follows: the interaction matrix A satisfying the above two assumptions (i–ii) is diagonally stable if and only if there exist constants Θj(i)>0 such that (22)iIjΘj(i)>|γj|Φ(sgn(γj),nj) for j=1,2,,JjJiΘj(i)=1 for i=1,2,,S.Arcak then outlined a systematic procedure for constructing Lyapunov functions based on the above stability condition. Notably, he illustrated the procedure with the GLV model for ecological systems [97].
  1. Download: Download high-res image (251KB)
  2. Download: Download full-size image

Fig. 10. Cactus and connected circle structures. (a) Digraph G(A) corresponding to the matrix A with simple cycles {1,4,1},{1,2,3,1},{1,7,8,1}, and {4,5,6,4} labeled as 1, 2, 3 and 4, respectively. Note that any pair of the cycles have at most one common node. (b) Undirected graph that describes which cycles in G(A) intersect. (c) Arborescence constructed on the undirected graph in (b) according to the broadcasting algorithm with cycle 1 selected as the root (shown in cyan). (d) Digraph with the connected circle structure where any pair of the cycles have at most one common edge or one common node.

Panels a, b, c were redrawn from [97].
Later on, Wang and Nešić extended the small gain condition to a more general circle structure called connected circles, where each pair of distinct simple cycles have at most one common edge or a common node [106], see Fig. 10d. In addition to the two assumptions (i) and (ii) made by Arcak, the authors further assumed that (iii) A has non-negative off-diagonal elements. Define Lj={ej(1),ej(2),,ej(nj1)} as the set of edges traversed by cycle j. Then the stability condition can be modified as follows: the interaction matrix A satisfying the above three assumptions (i–iii) is diagonally stable if and only if there exist constants Θj(i)>0 such that (23)iIjΘj(i)>γj for j=1,2,,JjJiΘj(i)=1 for i=1,2,,S,where γj=(p,q)LjApq is the gain for cycle j.

5.2.3. Rank-one structure

Recently, Simpson-Porco and Monshizadeh explored necessary and sufficient conditions for the diagonal stability of A with a rank-one network structure defined as A=D+S, where D=diag(d1,d2,,dS) is a positive diagonal matrix, and S=xy is a rank-one matrix for some x, yRS [105]. Suppose that y is a non-negative vector. The authors proved that the interaction matrix A with rank-one structure is diagonally stable if and only if (24)i=1S1di[xiyi]+<1,where xi and yi are the ith elements of x and y, respectively, and []+=max(,0) is the maximum operator with respect to zero. Significantly, they provided a theoretical stability analysis of automatic generation control in an interconnected nonlinear power system based on the above condition. While the rank-one structure may not naturally occur in real ecological networks, it might be useful in synthetic or specific ecological scenarios.

6. D-stability

The notion of D-stability was originally introduced by Arrow and McManus [108] and Enthoven and Arrow [110] in the late 1950s. D-stability can be considered as a weaker version of diagonal stability, see Fig. 11. A matrix A is called D-stable if for any positive diagonal matrix X, the matrix XA is stable (in terms of Re[λ(XA)]0). Clearly, the definition of D-stability is particularly relevant to the GLV model, where its community matrix is represented as M=XA. In this section, we first investigate the characterizations of D-stable matrices and then discuss existing results regarding D-stability in the context of the GLV model.

6.1. Characterizations of D-stable matrices

Explicitly characterizing D-stable matrices is only known for dimension less than or equal to four [113]. For example, a matrix AR4×4 is D-stable if and only if all the principal minors of A are positive and f(x,y,z)>0 for all positive x, y, and z, where f(x,y,z)=E1(DA)E2(DA)E3(DA)E3(DA)2E1(DA)2E4(DA),and D=diag(1,x,y,z). Here, Ei() denotes the sums of the order i principal minors of a given matrix. However, the problem becomes highly challenging as the dimension grows larger [79], [104], [113], [182]. Nonetheless, numerous sufficient conditions have been proposed [111], [202], [203]. Some of the better-known ones are:
  • If A is diagonally stable, then it is D-stable [108]. In fact, it is essentially the condition that Arrow and McManus offered.
  • If A is Metzler (i.e., Aij0 for ij) and all the principal minors of A are positive, then it is D-stable [204].
  • If there exists a positive diagonal matrix D such that AD=B satisfies Bii>ji|Bij| for i=1,2,,S, then A is D-stable [205]. The matrix A is referred to as quasi-dominant diagonal.
  • If A is triangular with Aii<0 for i=1,2,,n, then it is D-stable [111]. This is the most straightforward condition for D-stability.
  • If A is sign stable, then it is D-stable [111].
  • The element-wise product of P and A is stable for each positive definite symmetric matrix P [206].
Regrettably, none of these conditions is necessary for D-stability. More sufficient conditions can be found in [111].
  1. Download: Download high-res image (264KB)
  2. Download: Download full-size image

Fig. 11. Schematic containment plot of matrix stability notions. Detailed proofs of the relations among these stability notions can be found in [113].

This figure was redrawn from [162].

6.2. Results with the GLV model

Here, we present a series of findings and speculation regarding D-stability in the context of analyzing the GLV model. Chu offered a solvable Lie algebraic condition for the equivalence of four stability notions — linear stability, D-stability, total stability (defined in Section 7), and diagonal stability for the GLV model [109]. Given a Lie algebra L, let us define the inductive sequence as L(0)=L,L(i+1)={ABBA|A,BL(i)}.The Lie algebra L is referred to as solvable if there exists a positive integer k>0 such that L(k)={0}. The author proved that given the GLV model (9), if the two matrices, defined as (A+A)/2 and (AA)/2, generate a solvable Lie algebra, then the linear stability, D-stability, total stability, and diagonal stability of the interaction matrix A are equivalent. Moreover, although diagonal stability is not equivalent to Lyapunov stability in general, based on intensive numerical simulations, Rohr et al. conjectured that for mutualistic ecological systems captured by the GLV model, if the interaction matrix A is stable (in the sense of Lyapunov), then A is D-stable [59]. This conjecture has an important consequence for modeling mutualistic ecological systems with the GLV model, implying that if A is stable (in the sense of Lyapunov), any feasible equilibrium point is globally stable, corresponding to the unique fixed-point phase, as illustrated by Bunin [47]. Similarly, Gibbs et al. observed that large random matrices are D-stable almost surely (i.e., the set of positive diagonal matrices leading to instability has measure zero), so feasible unstable equilibrium points are very unlikely for the GLV model with random interactions [79].

7. Total stability

The notion of D-stability is also closely linked to another stability notion termed total stability. A matrix A is called totally stable if any principal sub-matrix of A is D-stable [95], [207]. Any principal sub-matrix of a totally stable matrix is also totally stable. In addition, total stability is a stronger notion compared to D-stability, but is weaker than diagonal stability, as illustrated in Fig. 11. Consequently, any totally stable matrix is inherently D-stable, and the class of totally stable matrices is closed under transposition and multiplication by a positive diagonal matrix. This implies that total stability holds simultaneously for both the interaction matrix A and the community matrix M of the GLV model [112]. A direct way to characterize a totally stable matrix is to check the D-stability of all its principal sub-matrices, which is computationally expensive. However, explicit characterizations are only known for S3 [112], [114], similar to diagonal stability and D-stability. A necessary condition of total stability is that all the principal minors of A of odd orders are negative, while those of even orders are positive [112]. Here, we provide an example of a 3-by-3 totally stable matrix that is not diagonally stable as follows (borrowed from [112]): A=10.320.9511.920.490.4951.It can be shown that all principal sub-matrices of A are D-stable. For example, the real parts of the eigenvalues of the following product XA=x100x210.30.951are equal to x1/2x2/2. On the other hand, (A+A)/2 contains a positive eigenvalue, so A is not diagonally stable.
Ecologically, the principal sub-matrix of the interaction matrix refers to the interactions between the subset of species that remains after removal or extinction of certain species. Therefore, total stability is related to the so-called “species-deletion stability” concept introduced in [208], implying that the property is preserved after the removal or extinction of any group of species from the initial composition [113]. In the context of the GLV model, the total stability of the interaction matrix suggests that the remaining species will reach a new locally stable equilibrium after the removal or extinction of certain species.

8. Sector stability

The notion of sector stability, proposed by Goh, focuses on the stability properties of semi-feasible equilibrium points (i.e., xi=0 for some i) for studying interesting processes in ecology, e.g., succession and extinction [58], [115]. A semi-feasible equilibrium point is called (locally) sector stable if every trajectory of the system that starts within a non-negative neighborhood remains in the same or an even larger non-negative neighborhood and eventually converges to that equilibrium point. The definition is analogous to that of (local) asymptotic stability. However, sector stability restricts the trajectories within a non-negative part of an open neighborhood of the equilibrium point.
Let us consider a generic population dynamics model defined as (25)ẋi=xifi(x1,x2,,xS)for i=1,2,,S, where fi have continuous partial derivatives at every finite point in the state space. Let I0={i|xi=0} and I1={i|xi>0}. Goh proved the semi-feasible equilibrium point x of the generic population model (25) is locally sector stable if all the eigenvalues of the community matrix defined as Mij=xifixj|x=x have negative real parts and fi(x)<0 for all iI0 [115]. Global sector stability results were also established by Goh through Lyapunov theory. Define U={x|xi>0 for iI1 and xi0 for iI0}. A semi-feasible equilibrium point is called globally sector stable if it is sector stable relative to the set U. The semi-feasible equilibrium point x of the generic population model (25) is globally sector stable if there exists positive constants d1,d2,,dS such that at every point in U, the function V(x)=i=1Sdi(xixi)fi(x)0,and it does not vanish identically along a solution of (25) except for x=x [115]. The result establishes valuable conditions for global sector stability in the GLV model (9): (i) there exists a positive diagonal matrix D=diag(d1,d2,,dS) such that DA+AD is negative semi-definite (i.e., A is diagonally semi-stable); (ii) the expressions ri+j=1SAijxj0 for all iI0; (iii) the function 12(xx)(DA+AD)(xx)+iI0dixi(ri+j=1SAijxj) does not vanish identically along any solution of (9) except for x=x.
However, the above global sector stability condition is difficult to verify in practice (especially when S>2). Thus, Goh came up with a conservative but simpler result. Suppose that there exists a constant matrix G such that (26)fi(x)xi|x=xGii<0for i=1,2,,S,|fi(x)xj|x=x|Gijfor ijin the set U. If all the leading principal minors of G are positive and fi(x)0 for all iI0, the semi-feasible equilibrium point x is globally sector stable [115]. As an illustrative example (borrowed from [115]), consider the following GLV model: ẋ1=x1(11.74x10.2x20.1x3)ẋ2=x2(1.20.8x1x20.2x3)ẋ3=x3(32x1x22x3).The model has four semi-feasible equilibrium points, which are (0,1,1), (0,6/5,0), (0,0,3/2), and (117/40,0,0). Let Gii=Aii for i=1,2,3 and Gij=|Aij| for ij such that (26) is satisfied. All the leading principal minors of G are positive. At the equilibrium point (117/40,0,0), f2(x)<0 and f3(x)<0. Therefore, this semi-feasible equilibrium point is globally sector stable with respect to U={x|x1>0 and x2,x30}, see Fig. 12.
  1. Download: Download high-res image (178KB)
  2. Download: Download full-size image

Fig. 12. Trajectories of the GLV model with three initial conditions closely positioned to the other three semi-feasible equilibrium points. The globally sector stable semi-feasible equilibrium point is marked by the red star.

The feasibility requirement for equilibrium points in complex ecological models such as the GLV model drastically restricts their parameter space, especially when assuming random interactions. Consequently, most equilibrium points are semi-feasible, which highlights the importance and necessity of sector stability. Furthermore, Goh’s results indicate that a semi-feasible equilibrium point will likely be globally sector stable if the strength of self-regulating interactions surpasses that of interspecific interactions [115].
  1. Download: Download high-res image (226KB)
  2. Download: Download full-size image

Fig. 13. Structural stability. (a) Andronov’s definition of structural stability. (b-d) Phase portraits of structurally unstable planar systems.

This figure was redrawn from [209].

9. Structural stability

The notion of structural stability of dynamical systems was first introduced by Andronov and Pontryagin under the name coarse (or rough) systems [116]. Different from the previous notions of stability which consider perturbations of initial conditions for a fixed dynamical system, structural stability concerns whether the qualitative behavior of the system trajectories will be affected by small perturbations of the system model itself [210], [211], [212], [213], [214]. In this section, we first discuss the mathematical definition of structural stability in dynamical systems. Then, we present various metrics that can be used to quantify the structural stability of ecological systems.

9.1. Mathematical definition

To formally define structural stability, we introduce the concept of topologically equivalence of dynamical systems. Two dynamical systems are called topologically equivalent if there is a homeomorphism h:RSRS mapping their phase portraits, preserving the direction of time. Consider two smooth continuous-time dynamical systems (i) ẋ(t)=f(x(t)) and (ii) ẋ(t)=g(x(t)). Both systems (i) and (ii) are defined in a closed region DRS, see Fig. 13a. System (i) is called structurally stable in a region D0D if for any system (ii) that is sufficiently C1-close to system (i) there are regions U, VD, and D0U, D0V such that system (i) is topologically equivalent in U to system (ii) in V, see Fig. 13a. Here, the systems (i) and (ii) are C1-close if their “distance”, defined as d=supxDf(x)g(x)+df(x)dxdg(x)dx,is small enough.
Andronov and Pontryagin offered sufficient and necessary conditions for a two-dimensional continuous-time dynamical system to be structurally stable [116]. A smooth dynamical system ẋ(t)=f(x(t)) with x(t)R2, is structurally stable in a region D0R2 if and only if (i) it has a finite number of equilibrium points and limit cycles in D0, and all of them are hyperbolic; (ii) there are no saddle separatrices returning to the same saddle, see Fig. 13b and c, or connecting two different saddles in D0, see Fig. 13d. This is often called the Andronov–Pontryagin criterion, which gives the complete description of structurally stable systems on the plane. It has been proven that a typical or generic two-dimensional system always satisfies the Andronov–Pontryagin criterion and hence is structurally stable [117]. In other words, structural stability is a generic property for planar systems. Yet, this is not true for high-dimensional systems. Later on, Morse and Smale established the sufficient conditions for an S-dimensional dynamical systems to be structurally stable [121], [122]. Such systems, often called Morse-Smale systems, have only a finite number of equilibrium points and limit cycles, all of which are hyperbolic and satisfy a transversality condition on their stable and unstable invariant manifolds.

9.2. Structural stability metrics

The notion of structural stability has been explored in diverse ecological systems. For instance, Recknagel investigated the structural stability of aquatic ecological systems by using the catastrophe theory and applied his approach as an aid in decision-making for water quality management [119]. In addition, the concept of structural stability have been heavily used in soil ecological systems [215], [216], [217], [218], [219]. Below, we survey various measures that can be used to quantify the structural stability of ecological systems.
Rohr et al. introduced a mathematical framework based on the concept of structural stability to elucidate the influence of network architecture on community persistence with the GLV model [59]. They proposed that an ecological system becomes more structurally stable as the area of the parameter space of the model expands, resulting in both a dynamically stable and feasible equilibrium. In particular, the authors investigated the range of conditions necessary for the stable coexistence of all species in mutualistic systems and showed that numerous observed mutualistic network architectures tend to maximize the volume of parameter space under which species coexist. This implies that having both a nested network architecture and a small mutualistic trade-off is one of the most preferable structures for community persistence.
Following Rohr’s framework, Grilli et al. studied the range of conditions necessary for feasible coexistence of large ecological systems using random matrix theory [55]. They quantified the structural stability of an ecological system using the GLV model as the volume of the feasibility region when varying intrinsic growth rates, which can be approximated by (27)Δ(1+1πCAμA(2μdSCAμA)μdSCA2μA2)S,where CA is the connectance of the interaction matrix A, μA is the mean of the off-diagonal elements of A, and μd is the mean of the diagonal elements of A. The authors further analytically predicted the range of coexistence conditions in more than 100 empirical ecological systems. Recently, using the above approximation of structural stability, Portillo et al. explored the correlation between structural stability and various network measures, such as centrality and modularity, and illustrated that optimal modularity has a negative impact on biological diversity (structural stability) in empirical ecological systems [220].
During the same time, Saavedra et al. proposed novel metrics analogous to stabilizing niche differences and fitness differences with the GLV model, which can measure the range of conditions compatible with multi-species coexistence, i.e., structural stability [120]. The structural analogs of the niche and fitness differences are defined as (28)Ω(A)=|det(A)|π/2SR0SexAAxdxandϑ=arccos(rrcrrc),respectively, where A is the interaction matrix, r is the vector of growth rates, and rc is the centroid of the feasibility domain defined as rc=1S(a1a1+a2a2++aSaS).Here, ai denotes the ith column of A. Therefore, feasible solutions can be fulfilled as long as r is inside the cone defining the domain of feasibility D(A)={rR>0S such that A1r>0}. In other words, the structural analog of the fitness difference θ is small enough relative to the structural analog of the niche difference Ω(A). The authors further applied their structural approach to a field system of annual plant competitors occurring on serpentine soils.
Saavedra and his coauthors have utilized similar approaches to study the structural stability of various ecological systems [60], [123], [221], [222], [223], [224], [225], [226], [227]. For example, they quantified structural stability using the quantity (1cos2(ϑ))/(cos2(ϑ)), which measures how big the deviations are from the structural vector (geometric centroid) compatible with a positive stable equilibrium [60]. They further proved that the smaller the level of global competition, the broader the conditions for having feasible solutions. In addition, Song and Saavedra proposed a measure of structural stability using the quantity Ω(A)1/S with a transformation 2AA=Σ1, which offers an approximation to the level of external perturbations tolerated by an ecological system [123]. The authors further found that their measure is the only consistent predictor of changes in species richness in empirical ecological systems among different ecological and environmental variables.
Lately, Pettersson et al. developed a metric called instability to quantify the proximity to collapse and level of structural stability with the GLV model, which provides deep insights into the dynamics and limits of stability and collapse of ecological systems [118]. The instability metric is defined as (29)ϰ(s)=σAσf(Spred(s))σc(Spred(s))σf(Spred(s)),where σA, σf, and σc are the standard derivations of interspecific interaction strengths, the first extinction boundary, and the collapse boundary, respectively. Here, Spred(s) is the predicted initial biodiversity in terms of the actual biodiversity s. The metric ϰ[0,1], and it can be computed from observable quantities of an ecological system. Notably, a higher ϰ value indicates lower structural stability of a system, which increases the likelihood of collapse due to perturbations or external pressures.

10. Stability analysis for systems with higher-order interactions

Various complex networks, such as ecological networks, social networks, biological networks, and chemical reaction networks, exhibit higher-order interactions, in which the interactions go beyond pairwise relationships and involve multiple nodes or entities at the same time [228], [229], [230], [231], [232], [233], [234], [235], [236]. In particular, in ecological networks, species interactions frequently occur in higher-order combinations, where the relationship between two species is influenced by one or more additional species [15], [61], [237], see Fig. 14. From the modeling perspective, higher-order interactions can be explicitly modeled by adding higher-order terms to the classical GLV model. Alternatively, they can be implicitly captured through consumer-resource models, i.e., consumer species associated with the same resource form a higher-order interaction. While the significance of higher-order interactions has been recognized, their comprehensive impact on the stability of ecological systems has not been fully understood.
Pioneering research has shown that one type of third-order interaction, in which one species mitigates the negative interactions between two others, can stabilize well-mixed ecological systems with particular network configurations [238]. In addition, increasing the order of interactions among species can weaken the destabilizing effect and amplify the variance of species abundances at the equilibrium [239], [240]. In this section, we first introduce the GLV model with higher-order interactions and present various results related to the stability of ecological systems with higher-order interactions. We subsequently outline two potential approaches for analytically determining the stability of the GLV model with higher-order interactions via polynomial systems theory. Finally, we discuss various consumer-resource models for implicit higher-order interactions and their stability properties. Note that the stability results for ecological systems with higher-order interactions involve linear stability analysis, Lyapunov analysis, or pure simulations.
  1. Download: Download high-res image (187KB)
  2. Download: Download full-size image

Fig. 14. Examples of pairwise, third-order, and fourth-order interactions in ecological systems. For pairwise interactions, species 1 could produce an antibiotic, inhibiting the growth of species 2. For third-order interactions, species 3 might degrade the antibiotics produced by species 1, thus alleviating the inhibitory effect of species 1 on species 2. For fourth-order interactions, the activity of the antibiotic-degrading enzyme by species 3 may in turn be inhibited by compounds produced by species 4.

This figure was redrawn from [61].

10.1. GLV model with higher-order interactions

The dynamics of ecological systems with higher-order interactions are often described by the GLV model with higher-order interactions [241], [242] defined as (30)ẋi(t)=xi(t)[ri+j=1SAijxj(t)+j=1Sk=1SBijkxj(t)xk(t)+j=1Sk=1Sl=1SCijklxj(t)xk(t)xl(t)+]for i=1,2,,S. Here, BRS×S×S is a third-order tensor whose off-diagonal elements represent the effect that species j and k has upon species i. CRS×S×S×S is a fourth-order tensor whose off-diagonal elements represent the effect that species j, k, l has upon species i. In fact, the GLV model with higher-order interactions (30) belongs to the family of polynomial dynamical systems [71]. AlAdwani and Saavedra [241] determined the number of non-trivial equilibrium points of the GLV model with higher-order interactions (30) based on Bernshtein’s theorem [243] from algebraic geometry. However, the stability of those non-trivial equilibrium points remains unknown. So far, most of the stability analysis of higher-order interactions depends on linearization or numerical simulations of the GLV model with higher-order interactions or other similar models.
Bairey et al. found that the classical relationship between diversity and stability is reversed when considering high-order interactions by simulating the dynamics of ecological systems using a replicator model [61], which is defined based on the GLV model with higher-order interactions. Denote the bracket part in the model (30) as fi(x). The replicator model is defined as (31)ẋi(t)=xi(t)[fi(x)j=1Sxjfj(x)],where the extra term ensures x to be solved in the form of relative abundances. The authors further derived a stability criterion by examining the effective pairwise interactions through linearizing the higher-order interactions around an equilibrium point. Suppose that the elements of the interaction matrix A, third-order interaction tensor B, and fourth-order interaction tensor C are randomly drawn from Gaussian distributions with means zero and variances ϖ, ϱ, and ς, respectively. For simplicity, assume that third-order elements Bijk are non-zero for i>j>k (similarly with other higher-order terms). An ecological system modeled by the replicator model (31) is stable with high probability if (32)ϖ+ϱS+ςS2+<1S.If the assumption on non-zero elements is removed, the higher-order interaction strengths in (32) are multiplied by constant factors, but the overall scaling remains unaffected. Therefore, interactions of any given order destabilize systems, and the impacts of different orders are additive, scaling differently with diversity. The combination of these interactions establishes both a lower and an upper bound on the number of species, which can be computed as Smin,max=1ϱ±(1ϱ)24ϖς2ϖ.When the discriminant vanishes, i.e., 1ϱ=2ϖς, the two bounds narrow a defined optimal number of species S=ς/ϖ.
Later on, Grilli et al. explored the role of higher-order interactions in competitive ecological systems [64]. Particularly, the authors simulated different competitive ecological systems with only third-order interactions of the form: (33)xi̇(t)=xi(t)[j=1Sk=1SBijkxj(t)xk(t)]for i=1,2,,S, which can be viewed as a special case of the higher-order GLV model. The interaction tensor B is defined from pairwise interactions and can be computed as Bijk=2HijHikHjiHjkHkiHkj where the first term represents the probability of species i beats both species j and k, and the remaining two terms represent the probabilities that either species j or k dominate. Here, the matrix H encodes the dominance relationships among species, i.e., Hij represents the probability that species i outcompetes species j. Compared to the pairwise dynamics, the equilibrium point of (33) is unchanged, and more importantly, it is globally stable. Therefore, they concluded that incorporating higher-order interactions into competitive ecological systems enhances stability, making species coexistence resilient to perturbations of both population abundance and parameter values [64].
Singh and Baruah derived general rules for species coexistence modulated by higher-order interactions by simulating the higher-order GLV model [242]. They demonstrated that negative higher-order interactions can promote coexistence if these interactions strengthen intraspecific competition more than interspecific competition, while positive higher-order interactions can stabilize coexistence across a wide range of fitness differences, disregarding differences in strength of interspecific and intraspecific competition. Recently, Gibbs et al. also exploited the higher-order GLV model and found that when interspecific higher-order interactions became excessively harmful compared to self-regulation, the stability of coexistence in diverse ecological systems will be compromised [63]. Furthermore, they found that coexistence can also be lost when these interactions are weak and mutualistic higher-order effects became dominant [63].

10.2. Two potential approaches

As mentioned, the GLV model with higher-order interactions (30) belongs to the family of polynomial dynamical systems. Thus, polynomial systems theory can be leveraged to analytically study the stability properties of the model. We summarize two recent approaches that can analytically determine the stability of homogeneous polynomial dynamical systems [71], [74]. For simplicity, let us assume that a higher-order ecological system with S species only contains pth-order interactions. The resulting GLV model with higher-order interactions thus becomes homogeneous (of degree p), similar to the form (33). In fact, every polynomial dynamical system can be homogenized by introducing additional variables. Therefore, the two approaches can be in principle applied to analyze the stability of the GLV model with higher-order interactions.

10.2.1. Lyapunov approach

Ali and Khadir proved that the existence of a rational Lyapunov function, defined as the ratio of two polynomials, is both necessary and sufficient for asymptotic stability of a homogeneous polynomial dynamical system at the origin [74]. The rational Lyapunov function takes the form of V(x)=q(x)(i=1Sxi2)r,where r is a non-negative integer and q(x) is a homogeneous positive definite polynomial of degree 2r+2. Since the Lyapunov inequalities on both the rational function and its derivative have sum of squares certificates, V(x) can always be found by semi-definite programming [244].
The authors further proved that an S-dimensional homogeneous polynomial dynamical system of degree p is asymptotically stable (in the sense of Lyapunov) at the origin if and only if there exist a non-negative integer r, a positive even integer s, with 2r<s, and symmetric matrices PI and QI, such that n(x),Qn(x)=2x2J(m(x))Pm(x),f(x)+2rm(x)Pm(x)x,f(x),where m(x) and n(x) denotes the vector of monomials in x of degree s2 and s+p+12, respectively, and J(m(x)) denotes the Jacobian matrix of m(x) [74]. However, solving this hierarchy of semi-definite programs requires trying all possible combinations of s and r, and the degree s cannot be bounded solely based on the values of S and p.

10.2.2. Tensor decomposition approach

Chen exploited tensor algebra to determine the stability of homogeneous polynomial dynamical systems [71]. First, the author proved that every S-dimensional homogeneous polynomial dynamical system of degree p can be equivalently represented by a (p+1)th-order S-dimensional tensor, denoted by ARS×S×p+1×S. Suppose that A is orthogonally decomposable, i.e., (34)A=i=1Sλivivip+1vi,where λi are often referred to as the Z-eigenvalues of A with the corresponding Z-eigenvectors vi [245]. The stability properties of the homogeneous polynomial dynamical system can be readily obtained through the Z-eigenvalues, similar to linear stability. Let the initial condition x0=i=1Saivi. The equilibrium point x=0 is (i) stable if and only if λiaip10 for all i; (ii) asymptotically stable if and only if λiaip1<0 for all i; (iii) unstable if and only if λiaip1>0 for some i. When p is odd, the initial condition can be ignored, and the stability conditions only depends on the Z-eigenvalues (exactly same as the case of linear stability).
The tensor decomposition approach provides straightforward criteria for determining the stability of homogeneous polynomial dynamical systems compared to the Lyapunov approach. Notably, Chen successfully employed this approach to analytically determine the stability of an ecological system described by the GLV model with higher-order interactions (considering third-order interactions only) with supply rates [71]. The model is defined as (35)ẋi(t)=xi(t)[j=1Sk=1SBijkxj(t)xk(t)]+si,where si is the supply rate for species i. The stability of the model can be obtained for orthogonal decomposable B. However, not all tensors can be decomposed in the form of (34), which limits the applicability of this approach [245]. Nevertheless, we believe that both approaches (the Lyapunov and tensor decomposition approaches) have a significant potential for analytically understanding the stability of higher-order dynamics in ecological systems.

10.3. Implicit higher-order interactions

Implicit higher-order interactions have been considered in consumer-resource models [246]. MacArthur’s consumer-resource model is one of the most commonly used models for studying consumer-resource interactions [41], [42]. It is defined as (36)ẋi(t)=xi(t)[θil=1NCliyl(t)νi]ẏk(t)=yk(t)[ιk(1yk(t)Kk)j=1SCkjxj(t)]for i=1,2,,S and k=1,2,,N, where xi(t) is the abundance of species i, yk(t) is the abundance of resource k, θi0 is the efficiency rate of species i for converting the consumed resources into biomass, Cli is the rate at which species i consumes resource l, νi is the mortality rate for species i, and ιk and Kk are the maximal growth rate and carrying capacity of resource k, respectively.
In MacArthur’s original formulation where the consumers and resources have different timescales, it can be shown that a feasible equilibrium point implies global stability by using a minimization principle [53], [247], [248]. Later on, the global asymptotic stability of MacArthur’s consumer-resource model with consumers and resources evolved in the same timescale has been proved by Case and Casten [58]. Notably, they rewrote MacArthur’s consumer-resource model into the GLV form (9) with x̃=x1xSy1yN,r̃=ν1νSι1ιN,Ã=00θ1C11θ1CN100θSC1SθSCNSC11C1Sι1/K10CN1CNS0ιN/KN,where x̃, r̃, and à are the composite vector (or matrix) of species abundances, intrinsic growth rates, and interactions, respectively. Based on the formulated GLV model, they successfully constructed a Lyapunov function to demonstrate that any feasible equilibrium point is globally asymptotically stable.
Aparicio et al. explored the feasibility domain of MacArthur’s consumer-resource model [53]. Based on Case and Casten’s result, the feasibility domain can be automatically treated as the stability domain in MacArthur’s consumer-resource model. Specifically, the authors utilized the structural stability metrics as defined in (28) to compute the feasibility domain. Their findings unveil that the feasibility of MacArthur’s consumer-resource model diminishes with the pool size of the consumers increases, but the expected fraction of feasible consumers rises in this scenario [53]. Conversely, if resources exhibit linear growth, an increase in the resource pool reduces the feasibility of the model and the expected fraction of feasible consumers. However, if the resources increase logistically, the trend is reversed.
In the following, we present several results on the stability of consumer-resource models with various structure.

10.3.1. Impacts of cascade structure

Cascade structure plays a critical role in the stability of ecological systems [249]. Schreiber investigated the global stability of consumer-resource cascades (serially arranged containers with a dynamic consumer population that receives a flow of resources from the preceding container) [250]. The cascade model is defined as (assuming the numbers of consumers and resources are the same) (37)ẋi(t)=xi(t)[θ(xi(t))h(yi(t)/xi(t)b)ν(xi(t))]ẏi(t)=χ(yi1(t)yi(t))xi(t)h(yi(t)/xi(t)b)for i=1,2,,S, where θ is the assimilation efficiency of the consumer dependent on xi, ν is the mortality rate of the consumer dependent on xi, χ is the flow rate, and h(yi/xib) is the functional response for b[0,1] [251]. To ensure both biological realism and mathematical tractability, it is assumed that (i) h:RR is C1 with h(0)=0 and h(x)>0 for all x0; (ii) The limit h=limxh(x) exists and is finite; (iii) limxxh(x)=0; (iv) 0<θ1 and ν>0 for all xi0; (v) θ(xi)0 and ν(xi)0 for all xi; (vi) θ(0)hν(0)>0. Under these assumptions (i)–(vi), the author proved that there exists an equilibrium point of the model which is globally stable. Hence, the ratio-dependent functional response can lead to persistence of the consumer population in all containers [250].

10.3.2. Impacts of mutualistic interactions

Microbial species capable of both resource consumption and production have the potential for mutualistic interactions [252], [253]. Butler and O’Dwyer explored the stability of a consumer-resource model with mutualistic interactions by explicitly considering the resources that microbes both consume and produce [254]. Their model is defined as (assuming the numbers of consumers and resources are the same) (38)ẋi(t)=xi(t)[θj=1SCjiyi(t)j=1SRjiνi]ẏi(t)=ιiyi(t)j=1SCijxj(t)+j=1SRijxj(t)for i=1,2,,S, where Rij is the rate of the production of resource i by consumer j, and other notations defined similarly as those in MacArthur’s consumer-resource model. For simplicity, they assumed that C=cI meaning each consumer specializes on a single resource. In addition, the equilibrium points for the consumers and resources are simplified as x1 and y1, where 1 is the all-one vector, achieved by tuning the influx and mortality rates. This setting corresponds to the most general case of purely mutualistic interspecific interactions. The authors proved that feasible equilibrium points of the model are locally stable if (39)(jiRij)2<cxθ(cyRiicx4θ)for all i=1,2,,S, see Fig. 15. Consequently, if microbes consume resources without producing them, any feasible equilibrium point will be locally stable. Yet, in the presence of cross-feeding, stability is no longer assured, where it can be obtained only when mutualistic interactions are sufficiently weak or when all pairs of taxa reciprocate each other’s assistance.
  1. Download: Download high-res image (315KB)
  2. Download: Download full-size image

Fig. 15. Eigenvalue distributions of the Jacobian matrix of the consumer-resource model with mutualistic interactions. (a) The randomly sampled production matrix R does not satisfy the condition (39), so there are positive eigenvalues. (b) The randomly sampled production matrix R satisfies the condition (39), so all eigenvalues are negative.

This figure was redrawn from [254].

10.3.3. Impacts of warming temperature

Temperature can affect the stability of ecological systems by causing changes in metabolic demands, ingestion rates, growth rates, and others [255], [256], [257], [258], [259]. Synodinos et al. used the classical Rosenzweig–MacArthur model [260] with a Holling type II functional response to study how temperature will alter stability [261]. The model is defined as (assuming one consumer and one resource) (40)ẋ(t)=x(t)[θoy(t)1+owy(t)ν]ẏ(t)=y(t)[ι(1y(t)K)ox(t)1+owy(t)],where o is the consumer attack rate, w is the handling time, ν is the mortality rate of the consumer, ι is the intrinsic growth rate of the resource, and K is the carrying capacity of the resource. They analyzed the impacts of temperature on stability using the tendency for intrinsic oscillations in population abundances, which can be defined as φ=owK, i.e., the product of the three temperature-dependent parameters. The larger φ, the lower the stability of the ecological system. The authors found four different warming-stability patterns from various empirical ecological systems: stability increases, decreases, is hump-shaped or U-shaped with temperature.

10.3.4. Impacts of dynamic switching

The switching of the distribution of the rates of consumption between two resources is common in nature [262], [263], [264]. Gawronski et al. introduced a dynamic switching mechanism in the classical Rosenzweig–MacArthur consumer-resource model with one consumer and two resources [265]. The model is formulated as (41)ẋ(t)=x(t)[oθ1y1(t)+θ2y2(t)1+w(θ1y1(t)+θ2y2(t))ν]ẏ1(t)=y1(t)[(1L11y1(t)L12y2(t))oθ1x(t)1+w(θ1y1(t)+θ2y2(t))]ẏ2(t)=y2(t)[(1L21y1(t)L22y2(t))oθ2x(t)1+w(θ1y1(t)+θ2y2(t))],where Lij represents the limitation of growth of resource i imposed by resource j. The dynamic switching is driven by maximizing the consumer’s population x, i.e., θ̇1=(1/Ϝ)dx/dθ1, where Ϝ is the characteristic time. Through numerical simulations, the authors found that oscillations of the consumer and the mutually synchronized two resources, which emerge at θi=0.5, become unstable for large or small θi. Additionally, the model converges to a stable equilibrium point when either θ1>0.5 and y1<y2 or the opposite. This suggests that the consumer is unable to switch its preferred resource once it has been chosen.

10.3.5. Impacts of periodic environments

Periodic fluctuations and signals are widespread and strongly influence the structure and function of ecological systems [266], [267]. Bieg et al. investigated the role of periodic environmental forcing on consumer-resource interactions [268]. They used an extension of the classic Rosenzweig–MacArthur consumer-resource model [260] defined as (assuming one consumer and one resource) (42)ẋ(t)=x(t)[θoy(t)ψ+y(t)ν]ẏ(t)=y(t)[ι(1y(t)Kmean+Kforce(t))ox(t)ψ+y(t)],where ψ is the half-saturation coefficient, Kmean is the average resource carrying capacity, and Kforce(t) is a sinusoidal function capturing periodic forces defined as Kforce (t)=Asin(2Pπt).Here, A is the amplitude of the forcing, and P is the forcing speed. Through local stability analysis, the authors found that the forcing speed significantly affects the stability of the model, with fast environmental forcing having a stabilizing effect. Their results further suggest that periodic fluctuations from climate change may cause sudden shifts in ecological dynamics and stability.

10.3.6. Impacts of consumption threshold

The idea of consumption threshold is borrowed from the basic reproduction in epidemiological models [269]. Duffy and Collins introduced this notion to determine the stability of the classical Rosenzweig–MacArthur consumer-resource model [270], which is defined similarly as (assuming one consumer and one resource) (43)ẋ(t)=x(t)[θoy(t)ψ+y(t)ν]ẏ(t)=y(t)[ι(1y(t)K)ox(t)ψ+y(t)].The consumption threshold is defined as C0=θoKν(ψ+K). Ecologically, it represents the parameter combination that guarantees the minimum resource consumption required for the consumer to survive. The authors proved that the trivial equilibrium point (i.e., x=y=0) is always unstable, and the semi-feasible equilibrium point (i.e., x=0 and y=K) is locally and globally asymptotically stable if C01. Furthermore, They also demonstrated that the model undergoes a Hopf bifurcation at C0=Kψ+K(1+(ν+θo)ψνK),and the non-trivial equilibrium point (i.e., x>0 and y>0) exists for C0>1 and is locally stable if C0<C0. In summary, when the consumption threshold C01, only the resource can survive, while when C0>1, both the consumer and resource can coexist.

10.3.7. Impacts of delayed age structure

The introduction of delay in the consumer-resource models plays a significant role in their stability. Akimenko studied the stability of a delayed age-structural consumer-resource model with N resource patches and consumers that can move freely between the patches [271]. The resource dynamics of the model is defined as (44)ẏi(t)=yi(t)[ιi(1yi(t)Ki)ϕiXˆ(t)1+ψiXˆ(t)]for i=1,2,,N, where ϕi>0 is the search rate, ψi0 is the saturation coefficient, and Xˆ(t)=0adκ(a)Υ(a,t)da is the weighted quantity of consumers at time t. Here, κ(a) represents the age-specific consumer’s preference, and Υ(a,t) is the age-specific density of consumer population at age a and time t. On the other hand, the consumer dynamics Υ(a,t) is described by the delayed McKendrick–Von Foerster’s age-structured model defined as (45)xt+xa=ν(a,r(tτ))Υ(a,t),where ν(a,r) is the age- and calorie intake rate-dependent consumer mortality rate [272], [273], [274]. Combined with initial and boundary conditions, the author proved that the trivial equilibrium point (i.e., yi=0 and x(a)=0) of the model is unstable for all τ>0 and the semi-feasible equilibrium point (i.e., yi=Ki and x(a)=0) is locally asymptotically stable if the consumer’s basic reproduction number R(r)<1 and is unstable if R(r)1 for all τ>0. The consumer’s basic reproduction number dependent on the calorie intake rate w is defined as R(r)=amatamaxΞ(a,r)exp(0aν(a,r)da)da,where Ξ(a,r) is the age- and calorie intake rate-dependent consumer fertility rate, and amat and amax are the age of maturation and maximum age of reproduction, respectively. Furthermore, the non-trivial equilibrium point (i.e., yi>0 for some i and x(a)>0) exists if and only if R(r)=1. Define I0={i|yi=0}. The non-trivial equilibrium point is unstable for all τ>0 if I0 and at least one of the following conditions is satisfied: (i) there exists iI0 such that ψiϕi/ιi; or (ii) ψi<ϕi/ιi for all iI0 and XˆmaxiI0(ϕi/ιiψi)1. It is locally asymptotically stable for all τ>0 if one of the following conditions is satisfied: (iii) I0, ψi<ϕi/ιi for all iI0, and Xˆ>maxiI0(ϕi/ιiψi)1; (iv) I0=. Akimenko’s findings show that the digestion period (i.e., the delay parameter) τ does not destabilize the consumer population at the trivial, semi-feasible, or non-trivial equilibrium points [271].
Recently, El-Doma improved the stability criteria proposed by Akimenko in analyzing the delayed age-structural consumer-resource model (44), (45) [275]. Notably, the author proved that the non-trivial equilibrium point (with I0) of the model is locally asymptotically stable if X(ϕi/ιiψi)>1 and is unstable if X(ϕi/ιiψi)<1 for iI0. Clearly, the conditions are simpler than Akimenko’s.

10.3.8. A case study of larch moth interactions

Plant’s quality plays an important role in the population dynamics of consumers [276], [277]. Din and Khan explored the dynamics between budmoths and the quality of larch trees located in the Swiss Alps [278]. They proposed an improved discrete-time consumer-resource model to capture the interactions between the plant quality index x and the moth population density y defined as (46)x(t+1)=(1)+x(t)y(t)ψ1+y(t)y(t+1)=y(t)exp(ι(x(t)ψ2+x(t)y(t)K)),where 0<<1 is the plant vulnerability, is the moth’s maximal uptake rate of the plant, ψ1 and ψ2 are the half-saturation coefficients, ι is the growth rate of moths, and K is the carrying capacity of moths. The authors proved that the model has a non-trivial equilibrium point, and it is locally asymptotically stable if and only if (47)ψ1ψ2ιy(ψ2+x)2(ψ1+y)2+(+1)(2Kιy)K>0,+ψ1ψ2ιy(ψ2+x)2(ψ1+y)2ιyK<1.The model also undergoes a period-doubling bifurcation at the non-trivial equilibrium point. Significantly, Din and Khan validated their theoretical findings with experimental data.

11. Discussion

The six stability notions we reviewed, including linear stability, sign stability, diagonal stability, D-stability, total stability, and structural stability, are not mutually exclusive but rather complement each other, which offers a multifaceted approach to explore the stability of ecological systems. Next, we discuss the pros and cons of each stability notion and outline its potential future prospects in the context of higher-order interactions.
Linear stability analysis is the most widely used approach to study the stability of ecological systems due to its simplicity. It has been heavily used in ecological systems with various network structures, as presented in Section 3. However, it cannot precisely capture the global dynamics and complex behaviors of ecological systems that are highly nonlinear in reality. Nevertheless, it will be interesting to perform the linear stability analysis on more complex models, such as the GLV model with higher-order interactions. For example, how to express the community matrix M in terms of the interaction matrix A and the higher-order interaction tensors, and how to establish stability conditions based on higher-order interaction strengths and correlations with various higher-order network structures?
Sign stability analysis focuses on the qualitative nature of interactions in ecological systems, disregarding of their magnitudes. It is particularly advantageous when quantitative data on interaction strengths are missing. However, the conditions described in Section 4 are too stringent to realize in real ecological systems. Hence, it will be intriguing to explore simpler and more realizable conditions that are associated to the underlying network structure for the sign stability of an ecological system. Generalizing the notion to nonlinear models, e.g., the GLV model or the GLV model with higher-order interactions, also presents exciting avenues for further work.
Diagonal stability, D-stability, and total stability analyses contribute to enhancing the understanding of the stability of the GLV model through characterizing the interaction matrix A. Moreover, they provide an insightful exploration of how the network structure encoded in the interaction matrix impact the overall stability of an ecological system. Nevertheless, characterizing diagonally stable, D-stable, or totally stable matrices is increasingly difficult for high-dimensional systems. Therefore, it will be worthwhile to further explore various network structures that can be used to establish necessary or sufficient conditions for achieving diagonally stability, D-stability, or total stability with the GLV model. Extending the current results to the GLV model with higher-order interactions also holds promising avenues for future research.
Sector stability analysis is particularly useful when dealing with semi-feasible equilibrium points, which are prevalent in modeling ecological systems. Local sector stability can be readily assessed, but a definitive algorithmic approach for verifying global sector stability remains elusive. For example, how can we efficiently compute the crucial constant matrix G defined in Section 8? Moreover, it will be worthwhile to explore simplified conditions on the sector stability of other nonlinear ecological models such as the GLV model with higher-order interactions and MacArthur’s consumer-resource model.
Structural stability analysis assesses how the overall behavior of an ecological system is maintained under small perturbations, which is different from all the previous stability notions. Structural stability can offer unique insights into the feasibility and robustness of ecological systems. An interesting avenue for future research would involve comparing various structural stability metrics and understanding their implications for the underlying network structure of an ecological system. Nonetheless, all the current structural stability metrics are developed based on the GLV model. Thus, extending these metrics to more complicated models, such as the GLV model with higher-order interactions, would also be a valuable exploration.

12. Conclusion

In this article, we presented a comprehensive and systematic review of the literature concerning the stability of ecological systems. Notably, we inclusively surveyed various stability notions that arise in ecological systems, encompassing linear stability, sign stability, diagonal stability, D-stability, total stability, sector stability, and structural stability. We further explored the stability of ecological systems with higher-order interactions.
Additionally, applying those various stability analysis to empirical ecological systems requires accurate estimation of the species interaction matrix or community matrix. For example, the trophic flows between consumers and resources in a food web with the Ecopath model can be translated into the interaction coefficients of the GLV model [81]. However, inferring the interaction network for large ecological systems, e.g., microbial communities, is challenging. Alternative frameworks, such as sensitivity testing of species abundance in response to the presence of additional species [32], may prove valuable for assessing the connectance and interaction strength of the interaction matrix or community matrix. To facilitate research in this field, we have compiled a summary of the most commonly used databases for real ecological systems in Table 4.
As mentioned earlier, empirically applying this theory to real ecological systems poses a significant challenge. First, most of the stability criteria are derived from the random matrix theory, which cannot accurately capture real ecological systems. For instance, the food web topology studies have demonstrated the non-randomness of ecological systems, rendering to incorporate more realistic network structures [126]. Second, the dynamics governing real ecological systems may not adhere to the GLV model, or even unknown in many cases. Thus, the stability analysis of real ecological systems may not align with the predictions of the theory.

Table 4. Summary of the most commonly used databases for real ecological systems. Links to the databases are provided.

DatabaseDescriptionReferences
Web of LifeA graphical user interface for visualizing and downloading data on various ecological networks of species interactions.[279]
Network RepositoryA network repository containing hundreds of real-world networks including ecological networks.[280]
GloBIAn extensible and open-source infrastructure for importing, searching, and exporting species-interaction data.[281]
NEONMonitoring ecological systems, including freshwater and terrestrial systems across the United States.[282], [283]
Global WebAn online collection of food webs.N/A
DoPIComprising ecological networks of British pollinator–plant interactions from the published scientific literature or submitted datasets.[284], [285], [286]
Finally, these stability analyses could be applied to other systems beyond ecological systems, such as power grid and gene regulatory networks. For example, the stability of power-grid networks in the homogeneous state, where the voltage frequencies of generators are all equal to a constant reference frequency, requires the largest eigenvalues of the Jacobian matrix to be negative [287]. Similarly, the stability of gene regulatory networks in regimes where the lifetime of mRNAs is much shorter than that of proteins depends solely on the eigenvalues of the Jacobian matrix. Furthermore, in the case of random gene regulatory networks, the underlying system becomes unstable when either the system size exceeds a critical threshold or the regulation strength becomes too strong, as determined using random matrix theory, which aligns with May’s theory [288].

Nomenclature

NotationDefinition
x(t)Species abundance (state vector)
xi(t)Species abundance for species i
f(x,t)/f(x)Unknown nonlinear dynamics
STotal number of species (states)
MCommunity matrix (Jacobian matrix)
x0Initial abundance (initial condition)
xEquilibrium point
Frobenius norm
λ()Eigenvalues of a matrix
Re[]Real part of an eigenvalue
Im[]Imaginary part of an eigenvalue
BϵBall of size ϵ
V(x,t)/V(x)Lyapunov functions
0/0Positive/negative definiteness
IIdentity matrix
z(t)Deviation from equilibrium point
N(0,σ2)Normal distribution with mean zero and standard derivation σ
E[]Expectation of a random variable
Var()Variance of a random variable
μMean of off-diagonal elements of M
σStandard derivation of off-diagonal elements of M
CConnectance of M
dDiagonal elements of M
ρCorrelation of off-diagonal elements of M
ξDegree heterogeneity of a network
QModularity of a network
qDiffusion coefficient of dispersal
MdelayCommunity matrix with time delay
τTime delay
ddelayDiagonal elements of Mdelay
σdelayStandard derivation of off-diagonal elements of Mdelay
CdelayConnectance of Mdelay
λdelayEigenvalues of Mdelay
rIntrinsic growth rate vector
riIntrinsic growth rate of species i
AInteraction matrix
XDiagonal equilibrium matrix
μXMean of diagonal elements of X
μdMean of diagonal elements of A
μAMean of off-diagonal elements of A
ρACorrelation of off-diagonal elements of A
CAConnectance of A
1all-one vector or matrix
λoutlierOutlier eigenvalue
τcHopf bifurcation of time delay
WArea of the species living domain
ηi(t)Gaussian noise for species i
ζ(t)Gaussian white noise
FCorrelation matrix of ζ(t)
Φ(ω)Power spectral density of fluctuations with frequency ω
x(t)Species abundance in effective process
M(t)Average species concentration in effective process
G(t,t)Response function in effective process
η(t)Gaussian noise in effective process
z(t)Derivation from equilibrium point in effective process
v(t)deviation of noise in effective process
ζ(t)Gaussian white noise of unit amplitude in effective process
xia(t)Abundance of species i in local system a
BTotal number of local systems
SaTotal number of species in local system a
Det()Determinant of a matrix
DClass of diagonally stable matrices
SClass of infinite sector nonlinear functions
G(A)Digraph of A
JTotal number of simple cycles
njLength of cycle j
IjSet of nodes traversed by cycle j
JiSet of cycles that node i belongs to
γjGain of cycle j
LjSet of edges traversed by cycle j
Ei()Sums of order i principal minors
LLie algebra
ΔVolume of feasibility region
ΩStructural niche difference
ϑStructural fitness difference
rcCenter of feasibility domain
ϰInstability metric
σfStandard derivation of first extinction boundary
σcStandard derivation of collapse boundary
SpredPredicted initial biodiversity
BThird-order interaction tensor
CFourth-order interaction tensor
ϖVariance of pairwise interactions
ϱVariance of third-order interactions
ςVariance of fourth-order interactions
sisupply rate for species i
NTotal number of resources
yk(t)Abundance of resource k
θi/θ(xi)/θEfficiency rate of species i
CjiConsumption rate at which species i consumers resource j
νi/ν(xi)/νMortality rate for species i
ιk/ιMaximal growth rate for resource k
Kk/KCarrying capacity of resource k
χFlow rate
h()Function responses
RijProduction rate of resource i by consumer j
oConsumer attack rate
wHandling time
Lijlimitation of growth of resource i imposed by resource j
ϜCharacteristic time
ψi/ψHalf-saturation coefficient for species i
KmeanAverage resource carrying capacity
Kforce(t)Sinusoidal function capturing periodic forces
C0Consumption threshold
ϕiSearch rate for species i
κ(a)Age-specific consumer’s preference
Xˆ(t)Weighted quantity of consumers
XˆEquilibrium point of weighted quantity of consumers
Υ(a,t)Age-specific density of consumer population
r(t)Calorie intake rate
R(r)Consumer’s basic reproduction number
Ξ(a,r)Age- and calorie intake rate-dependent consumer fertility rate
Plant vulnerability
Moth’s maximal uptake rate of plant

CRediT authorship contribution statement

Can Chen: Writing – review & editing, Writing – original draft, Visualization, Methodology, Investigation, Conceptualization. Xu-Wen Wang: Writing – review & editing, Writing – original draft, Visualization, Methodology, Investigation, Conceptualization. Yang-Yu Liu: Writing – review & editing, Writing – original draft, Supervision, Methodology, Investigation, Funding acquisition, Conceptualization.

Declaration of Generative AI and AI-assisted technologies in the writing process

During the preparation of this work, the authors used ChatGPT in order to improve readability and language. After using this tool, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work was supported by the National Institutes of Health, United States R01AI141529. Xu-Wen Wang acknowledges the funding support from National Institutes of Health (K25HL166208).

References

Cited by (0)

View Abstract