Consider an ecological community composed of S species, whose population densities at any time t are denoted as a vector Y(t) of length S. Species dynamics are typically modeled by a system of non-linear autonomous ordinary differential equations:

where f = [f₁, f₂, …, f_S]^T is a set of functions describing the dynamics of the S species. A feasible equilibrium point of system (1) is a non-negative vector Y^* = [Y^*₁, Y^*₂, …, Y^*_S]^T such that:

The evolution of a small perturbation y₀ = ΔY(0) applied to the equilibrium Y^* can be approximated by the linear equation:

where M is the so-called “community matrix” (Levins, 1968; May, 1972, 2001): the Jacobian matrix evaluated at the equilibrium, i.e.,

An equilibrium is called “(locally asymptotic) stable” if any infinitesimal perturbation decays to zero eventually: for any ΔY(0) sufficiently small, we have lim_{t → ∞}ΔY(t) = 0. Moreover, if all eigenvalues of M, denoted by λ^M_i, i = 1, 2, …, S, have negative real parts, the equilibrium point is stable (May, 1972, 2001). Here we always order the eigenvalues of M such that ℜ(λ^M₁) ≤ ℜ(λ^M₂) ≤ ··· ≤ ℜ(λ^M_S). Thus if ℜ(λ^M_S) < 0, the equilibrium is stable.

The stability of an equilibrium defines the long-term response to small perturbations. However, the real parts of M's eigenvalues do not provide information on the instantaneous response to perturbations. Among stable equilibria, the instantaneous dynamics can differ dramatically, even if all perturbations will eventually die out.

Neubert and Caswell (1997) firstly introduced to ecology the concept of reactivity, which quantifies the maximum amplification rate of perturbations. For reactive equilibria, some perturbations will immediately grow in magnitude before eventually decaying. Conversely, if all small perturbations decay immediately, such equilibria are non-reactive. As such, reactivity is defined as the maximal initial amplification rate, which can be calculated from:

where ||·|| is the norm operator. Neubert and Caswell also showed that the maximum initial amplification rate is simply λ^H_S, the largest eigenvalue of H, where H is the symmetric part of the community matrix M, defined as H=M+MT2, where M^T is the transpose of M. Note that H is symmetric and therefore all its eigenvalues are real. An equilibrium is reactive if λ^H_S > 0 (i.e., positive instantaneous amplification rate) and non-reactive if λ^H_S < 0 (i.e., negative instantaneous amplification rate) (Neubert and Caswell, 1997). For any matrix M, we have that λ^H_S ≥ ℜ(λ^M_S) (Snyder, 2010): the largest real part for the eigenvalues of M is bounded from above by the largest eigenvalue of the symmetric matrix H=M+MT2, and the largest imaginary part is bounded from above by the largest eigenvalue of the symmetric matrix H′=M−MT2 (Wolkowicz and Styan, 1980). As such, unstable equilibria (ℜ(λ^M_S) > 0) are certainly reactive (λ^H_S > 0), whereas stable equilibria (ℜ(λ^M_S) < 0) can be either reactive or non-reactive.

As when studying stability (Tang et al., in press), the diagonal entries of M are assumed to have a negative expectation 𝔼 (M_ii) < 0; and we denote 𝔼 (M_ii) = −d for some positive number d. The diagonal entries are sampled independently from a normal distribution with mean −d and finite variance σ²_d. The off-diagonal entries of M are sampled as independent pairs according to a bivariate distribution of (X₁, X₂) in the following way: (1) with probability C/2, (M_ij, M_ji)_{i > j} is sampled from the distribution of (X₁, X₂); (2) with probability C/2, (M_ij, M_ji)_{i l > j} is sampled from the distribution of (X₂, X₁); (3) with probability (1 − C), (M_ij, M_ji)_{i > j} = (0, 0). The symmetric part of M is computed from the definition H=M+MT2. Note that M constructed in this way can describe different types of communities when we vary the underlying bivariate distribution. For example, when (X₁, X₂) is chosen to follow a bivariate mean-zero Gaussian distribution, with two independent components, M can represent the collection of randomly-assembled communities studied by May (1972). If the distribution is defined such that the two components X₁ and X₂ always have opposite signs, then M can represent a food web, containing only consumer-resource interactions, as in Allesina and Tang (2012).

For a random matrix M whose entries are independently sampled from a statistical distribution, five quantities are essential for determining the largest real part of its eigenvalues: (1) S, the dimension of M, which equals the number of species in the network, (2) E, the expectation of the off-diagonal entries, describing inter-specific interaction strengths, (3) −d, the expectation of the diagonal entries, describing intra-specific interaction strengths, (4) V, the variance of the off-diagonal entries, and (5) ρ, the pairwise interaction correlation:

In Allesina and Tang (2012) and Tang et al. (in press), we showed that there is always one eigenvalue of M (denoted by λ^M_R) whose real part is close to the expected row sum of M:

When S is sufficiently large, the other (S − 1) eigenvalues of M are approximately uniformly distributed on an ellipse centered at (−E −d, 0), whose horizontal and vertical axes are about 2SV(1 + p) and 2SV(1 − p), respectively. Thus the rightmost one of these (S − 1) eigenvalues (denoted by ℜ(λ^M_EL)) can be estimated using the “center” plus the semi-length of the horizontal axis:

Depending on the sign and the magnitude of ℜ(λ^M_R), stability may be determined by either ℜ(λ^M_R) or ℜ(λ^M_EL), i.e., the largest real part for all eigenvalues of M is just the larger one of these two eigenvalues:

When E > 0, ℜ(λ^M_R) grows linearly in S whereas ℜ(λ^M_EL) grows only sub-linearly, so if S is sufficiently large, for example, S≫V(1+ρ)2E2, we have:

In this case, λ^M_R lies on the right of the ellipse formed by the other (S − 1) eigenvalues, and the largest real part of the eigenvalues of M is thus ℜ(λ^M_R), i.e., for sufficiently large S,

When E ≤ 0, we always have

which means λ^M_R lies on the left of the ellipse. Hence the rightmost eigenvalue of M is λ^M_EL. In this case, for sufficiently large S, we have

A random, symmetric matrix, such as H, is a special case of general random matrices. As such, we may use the same approach to estimate the largest eigenvalue of H. We thus compute the corresponding statistical quantities for H, namely,

In particular, we would like to express these five quantities in terms of the statistics of M. First of all, the expectations of the diagonal and off-diagonal entries do not change by taking the symmetric part of M, because

However, the variance, the expectation of the products of pairwise interactions, and the pairwise correlation do differ in M and H. These quantities for H can be computed according to the relationship between M and H:

Similar to the case of stability, two eigenvalues of H are essential in determining the reactivity of the underlying equilibrium, one corresponding to the expected row sum (denoted by λ^H_R), and the other one corresponding to the largest one among the other (S − 1) eigenvalues forming an “ellipse.” In the case of symmetric matrices, since ρ^H = 1, the vertical axis of the ellipse is zero:

This means the ellipse actually degenerates to a line segment on the real axis. The density of the eigenvalues of a random symmetric matrix is described by Wigner's semicircle law (Wigner, 1958), so we denote the largest eigenvalue among the S − 1 eigenvalues (excluding λ^H_R) by λ^H_SC, where “_SC” stands for “semicircle.” Similarly to the computation of ℜ(λ^M_R) and ℜ(λ^M_EL), λ^H_R and λ^H_SC can be computed as follows:

These equations are analogous to those determining stability (Equations 6, 7): the only difference is that instead of using E, V, and ρ for matrix M, we take the corresponding values for matrix H.

Again, the largest eigenvalue of H is the larger one of λ^H_R and λ^H_SC:

If E > 0, for sufficiently large S, namely S≫2V(1+ρ)E2, we have λ^H_R > λ^H_SC

whereas if E ≤ 0, the largest eigenvalue of H is λ^H_SC and thus

We use numerical simulations to verify the derivation of the reactivity criteria. We choose S = 250, C = 0.2, and standard normal distribution (0, 1) for the marginal distributions of X₁ and X₂. We consider three different cases of the joint distribution (X₁, X₂). The first case is called the Random case, where X₁ and X₂ follow independent (0, 1) distributions and (M_ij, M_ji)_{i > j} can be (+, −), (−, +), (+, +), or (−, −). In this case, we have E = 0, V = C, and ρ = 0. The second case is called the Predator-Prey case, where we constrain the two entries in each non-zero pair (M_ij, M_ji)_{i > j} to have different signs, either (+, −) or (−, +), and this also gives E = 0, V = C, but a negative pairwise correlation ρ=E2V=−C𝔼|X1|𝔼|X2|C=−2π<0. The last case is the Mixture of mutualism and competition case (or “Mixture” case for short), where each non-zero pair (M_ij, M_ji)_{i < j} contains two entries of the same sign, either (+, +) or (−, −). In the Mixture case, we still have E = 0, V = C, but a positive pairwise correlation ρ=E2V=C𝔼|X1|𝔼|X2|C=2π>0. The joint distributions of (X₁, X₂) for the three cases are illustrated in Figure 1. To simulate the phase transition for stability and reactivity, we increase d from 0 to 15 (i.e., 𝔼 (M_ii) is between 0 and −15), for three levels of the variances of the diagonal entries σ²_d = 0, 0.5² and 1². For each combination of d and σ_d, we construct 1000 community matrices for each of the three cases, and for each matrix, we check stability and reactivity and estimate the probabilities of stability and reactivity, respectively.

Similar to the derivation for the stability criterion, in which the goal is to estimate the largest real part of the eigenvalues of M, the derivation of the reactivity criterion relies on estimating the largest eigenvalue of H when S is sufficiently large. The reactivity criteria are summarized and compared to their stability counterparts in Table 1. We see that the five quantities affect the largest eigenvalue of H in the same direction as they affect the largest real part of the eigenvalues of M: if changing one quantity while the others are held constant increases (or reduces) ℜ(λ^M_S), such change also increases (or reduces) λ^H_S.

If E > 0, for sufficiently large S, ℜ(λ^M_S) and λ^H_S are both of the order (S − 1)E − d, implying that the critical point defining the stability-instability boundary (i.e., when ℜ(λ^M_S) = 0) is asymptotically the same as that for the non-reactivity-reactivity boundary (i.e., when λ^H_S = 0). If E ≤ 0 and S is sufficiently large, ℜ(λ^M_S) and λ^H_S are approximately SV(1 + ρ) − E − d and 2SV(1+ρ) − E − d, respectively. Note that, in this case, since −1 ≤ ρ ≤ 1, we have

which yields ℜ(λ^M_S) ≤ λ^H_S. It means we necessarily have ℜ(λ^M_S) < 0 (“stable”) whenever λ^H_S < 0 (“non-reactive”). Our estimate for ℜ(λ^M_S) and λ^M_S are consistent with the fact that a non-reactive equilibrium is always stable (Snyder, 2010).

The assumption of large S is important in determining the rightmost eigenvalue when E > 0. Similar to the study of stability, the reason here is that although λ^H_R grows linearly in S and is eventually larger than λ^H_SC, which grows only sub-linearly, if E is small and S is not sufficiently large, we may still have λ^H_R < λ^H_SC and thus the largest eigenvalue of H is then λ^H_SC instead of λ^H_R. In practice, for a system of a given size S and mean interaction strength E, we need to consider the estimates for both eigenvalues, λ^H_R and λ^H_SC. Similar to the stability criteria derived in Tang et al. (in press), we thus combine the criteria for both E > 0 and E ≤ 0 cases and state the criteria in terms of the maximum of λ^H_R and λ^H_SC (Table 1, row 3).

Note that the expressions for ℜ(λ^M_S) and λ^H_S obtained above are the estimates for their average values. In other words, for fixed S, C, d, and the underlying bivariate distribution (X₁, X₂) (which determines E, E₂, V, and ρ), if one constructs random community matrices M repeatedly, Equations (6–9) provide estimates for the average positions of ℜ(λ^M_S) and λ^H_S among all constructions. For one single realization of M and H, the actual values of ℜ(λ^M_S) and λ^H_S will deviate from their mean estimates due to the randomness. Furthermore, the variations of ℜ(λ^M_S) and λ^H_S around their means decrease with the size of the network: when S is larger, the estimates given in Equations (6–9) are more accurate. Finally, if S ≫ 1 and E is small compared to SV(1 + ρ), we may simplify the expressions in Table 1 accounting for these facts (1) S − 1 ≈ S, (2) SV(1 + ρ) − E ≈ SV(1 + ρ), and (3) 2SV(1+ρ) − E ≈ 2SV(1+ρ).

From Table 1, we see that the inequalities define the boundary conditions for stability and reactivity. Intuitively, since non-reactive equilibria are necessarily stable but not vice versa, we may think non-reactivity as a stronger form of stability in the sense that all infinitesimal perturbations not only decrease to zero eventually, but also do so immediately.

For fixed S, V, E, and ρ, the negative mean intra-specific interaction d may be viewed as an internal source to stabilize the equilibrium. If d = 0, implying that the eigenvalues of M are centered at zero, we always have unstable and hence reactive equilibrium. When we increase d, which is equivalent to moving the center of the eigenvalues to the left, it is possible to achieve stability and non-reactivity. We denote by d_crit,1 and d_{crit, 2} the minimum (i.e., critical) average strengths of intra-specific interactions needed to achieve stability and non-reactivity, respectively. When E > 0 and S is sufficiently large so that ℜ(λ^M_R) and λ^H_R determine stability and reactivity, respectively, the minimum d required to stabilize or make the equilibrium non-reactive is d_crit,1 = d_crit,2 = (S − 1)E. However, in real ecological networks, the case E < 0 is more common (Tang et al., in press), especially when consumer-resource interactions are preponderant. This is because the positive effects of resources on consumers are usually smaller in magnitude than the negative effects of consumers on resources due to imperfect conversion efficiency. Moreover, adding a higher proportion of competitive interactions compared to mutualistic interaction can make the mean interaction strength even more negative. For this case, the minimum average strength of intra-specific interaction required to stabilize the equilibrium is d_{crit, 1} = SV(1 + ρ) − E, whereas to make the equilibrium non-reactive we need average strength of inter-specific interaction being at least d_{crit, 2} = 2SV(1+ρ) − E, greater than d_crit,1.

Our criteria are quite accurate even for communities of moderate sizes (Figure 2), as confirmed by numerical simulations. In our simulations, when we increase the mean intra-specific interaction strength d, we observe two sequential phase transitions as d passes through the two critical points. First, near d = d_crit,1, the probability of stability rapidly increases from 0 to 1, marking the first phase transition from instability to stability (Figure 2, upward triangles). After that, as d continues to increase, the probability of reactivity rapidly drops from 1 to 0 near d = d_crit,2, marking the second phase transition from reactivity to non-reactivity (Figure 2, downward triangles). As S becomes larger, the behavior for the eigenvalues of M and H will be closer and closer to that estimated according to random matrix theory, and, consequently, the ℜ(λ^M_S) and λ^H_S will be closer and closer to their mean positions (Tang et al., in press). These phase transitions will be sharper for larger S, and the behaviors of the transitions will be more like sudden jumps occurring at d_crit,1 and d_crit,2, respectively. The critical points predicted using our criteria, i.e., d_crit,1 and d_crit,2, are marked by the dashed vertical lines in each case in Figure 2. Conversely, if we start with a non-reactively stable network and gradually push it toward the unstable region, for example, by reducing the mean strength of the intra-specific interactions, the network will first become reactively stable and then unstable. This makes intuitively sense, since if all small perturbations would decay immediately (i.e., non-reactive, stable), then any small perturbation is equivalent to a smaller perturbation in a slightly different direction, which will also decay immediately. Therefore, if a network has a non-reactive, stable equilibrium, it has to first cross the “reactively stable” phase before reaching the unstable phase, assuming its state is changing continuously. Thus, if reactivity were to be detectable in empirical data, it would be an excellent candidate for an early-warning signal for the approach of a bifurcation.

From the inequalities in Table 1, we see that negative pairwise correlation is beneficial for both stability and non-reactivity. Given the potential of reactivity for signaling the approach of a bifurcation, it is important to know how “early” reactivity can warn us when the system is approaching the unstable state. This is actually measured by the distance between the phase transitions for reactivity and stability, i.e., the absolute difference δ = λ^H_S − ℜ(λ^M_S). When E > 0 and S is sufficiently large, ℜ(λ_S) ≈ λ^H_S and reactivity is asymptotically equivalent to instability. In other cases, either E ≤ 0 or E > 0 but S is not large enough so that we have ℜ(λ^M_S) ≈ SV(1 + ρ) − E − d and λ^H_S ≈ 2SV(1+ρ) − E − d, the difference is

depending on S, V, and ρ only. Here, we can see that the absolute distance δ changes monotonically with S and V, indicating that a higher biodiversity results in a larger distance between the two phase transitions. However, δ does not change monotonically with ρ. Computing the first derivative of δ with respect to ρ and equating it to zero, we get:

The second derivative at ρ^* is negative, which implies δ reaches its maximum when ρ=−12. Thus, for fixed S and V, the distance from the non-reactivity-reactivity phase transition to the stability-instability is large when the interactions between pairs of species maintain a moderate negative correlation (Figure 3).

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