Bounds on Transient Instability for Complex Ecosystems

Stability is a desirable property of complex ecosystems. If a community of interacting species is at a stable equilibrium point then it is able to withstand small perturbations to component species’ abundances without suffering adverse effects. In ecology, the Jacobian matrix evaluated at an equilibrium point is known as the community matrix, which describes the population dynamics of interacting species. A system’s asymptotic short- and long-term behaviour can be determined from eigenvalues derived from the community matrix. Here we use results from the theory of pseudospectra to describe intermediate, transient dynamics. We first recover the established result that the transition from stable to unstable dynamics includes a region of ‘transient instability’, where the effect of a small perturbation to species’ abundances—to the population vector—is amplified before ultimately decaying. Then we show that the shift from stability to transient instability can be affected by uncertainty in, or small changes to, entries in the community matrix, and determine lower and upper bounds to the maximum amplitude of perturbations to the population vector. Of five different types of community matrix, we find that amplification is least severe when predator-prey interactions dominate. This analysis is relevant to other systems whose dynamics can be expressed in terms of the Jacobian matrix.


Introduction
From the perspective of local stability analysis, if an ecosystem is close to a stable equilibrium point then the effect of a small perturbation, such as the loss of individuals from a population, will eventually decay and the system will return to its original equilibrium point [1,2]. But if the ecosystem is at an unstable equilibrium point then the perturbation will lead to the system settling at a new equilibrium point, possibly with fewer individuals or even species [3,4]. In theory, ecosystems with large numbers of species and interactions are more difficult to stabilise [5]. However, many ecosystems contain vast biodiversity [6,7]. Reconciling this finding with local stability analysis has motivated ecologists for over 40 years [8].
Recently, stability criteria were extended from randomly-assembled communities to include those with more realistic compositions of mutualistic, competitive and predator-prey interactions [9]. These criteria indicate that communities in which predator-prey interactions dominate are more likely to be stable. It was then shown, using empirical food webs, that the distribution and correlation of interaction strengths has a greater effect on stability than topology: how species interact with one another is more important than who they interact with [10,11].
Stability is a long-term concept: it indicates whether a system will, at some point in the future, return to the same state as before a perturbation [12]. Reactivity, on the other hand, indicates how a system will respond immediately after a perturbation has been applied [13][14][15][16][17]. A stable system can be non-reactive, meaning that a perturbation to species' abundances dies down immediately, or reactive, meaning that a perturbation is first amplified before eventually decaying (whether a particular perturbation is amplified in practice depends on which species are perturbed and by how much [13]). Reactivity criteria for large ecosystems indicate that communities on the verge of instability exhibit reactive dynamics [18], and identifying a system as reactive has been proposed as an early-warning signal for population collapse [19][20][21][22][23].
The starting point for deriving criteria for both stability and reactivity is the community matrix [24]. A spectral decomposition of the community matrix provides information on the asymptotic behaviour of the system for stability (t ! 1) and reactivity (t ! 0). But so far, little information has been extracted from the community matrix regarding transient dynamics: how the system evolves after a perturbation and before it either returns to equilibrium or becomes unstable [25][26][27].
Reactive dynamics are not possible if the community matrix M is normal, i.e., MM † = M † M, where M † is the adjoint of M [28]. But if M is a non-normal matrix, as is usually the case in analyses of realistic ecosystems, then transient dynamics may substantially differ from the asymptotic behaviour suggested by the eigenvalues of M. In addition, small changes to the entries of non-normal M can cause an otherwise stable matrix to become unstable [28]. In such cases, the dynamics implied by non-normal matrices are better described by pseudospectra, which detail the neighbourhood of eigenvalues in the complex plane for different average changes to the entries in M [29].
Here we formalise the transition from stability to instability in terms of pseudospectra. Using this approach, we consider the effect on dynamics of two kinds of perturbation: more commonly studied perturbations to the equilibrium abundance of species (to the population vector) and less commonly studied perturbations to the entries in M (which could be interpreted as uncertainty in, or small changes to, species' interaction strengths [30]). We describe critical values for community properties separating three regimes: stable and non-reactive dynamics, stable and reactive dynamics-'transient instability'-and unstable dynamics. We show that system dynamics at the boundary between non-reactive stability and transient instability can be affected by perturbations to entries of the community matrix. And, given a perturbation to the equilibrium abundance of species, we provide upper and lower bounds to the maximum amplification of such perturbations during transient instability. This allows us to sketch out the transient dynamics of complex ecosystems using only information from the community matrix. Finally, we compare the properties of community matrices representing ecological communities with five different types of interaction structure: random, mutualism, competition, mixture of mutualism and competition, and predator-prey.

Local stability analysis
Here we consider an ecological community of S species for which their population densities at time t are given by the vector Y(t), as in Tang & Allesina [18]. The dynamics of the population vector Y can be described by a system of coupled differential equations where f = [f 1 , f 2 Á Á Á, f S ] T is a vector of linear or nonlinear functions. An ecologically-relevant equilibrium point is a non-negative vector Y Ã such that The community matrix M is defined as which is the Jacobian matrix evaluated at an equilibrium point [24]. It is well known that an equilibrium point is (locally and asymptotically) stable if any infinitesimally small deviation, ΔY(0), eventually decays to zero, i.e., lim t!1 ΔY(t) = 0 [24]. In the vicinity of an equilibrium point, the time evolution of a perturbation can be described by Therefore, the spectrum of the community matrix M is clearly relevant for determining local stability. If Λ(M) is the set of eigenvalues of M, then an equilibrium point is stable if all eigenvalues have negative real part, i.e., Re(λ) < 0 8 λ 2 Λ(M) [5,9].

Generative models for community matrices
We parameterise community matrices using four quantities: S, C, μ and σ; where S, as above, is the number of species, C is the connectance (the fraction of realised interactions among species), μ is the strength of intraspecific interactions and σ is the standard deviation of the strength of interspecific interactions [9]. We assume that populations are self-regulating and so M ii = −μ, where μ > 0. Non-normal community matrices with different types of interactionrepresenting different types of ecological community-are generated by sampling off-diagonal entries (M ij , interspecific interactions) from different bivariate distributions. Having specified a particular distribution, stability criteria can be expressed in terms of S, C, μ and σ. Based on these criteria, it has been shown that predator-prey community matrices are the most stable, followed by random, competition, mixture and mutualism [9]. Generative models for these community matrices are described below.
Random. Each off-diagonal entry is sampled independently from a normal distribution N ð0; sÞ with probability C, and otherwise M ij = 0 with probability 1 − C.
Mutualism. Each off-diagonal pair (M ij , M ji ) is sampled from a half-normal distribution jN ð0; sÞj with probability C, and both entries are zero otherwise. These community matrices have a (+,+) sign structure for off-diagonal pairs.
Competition. Each off-diagonal pair (M ij , M ji ) is sampled from a half-normal distribution ÀjN ð0; sÞj with probability C, and both entries are zero otherwise. These community matrices have a (−,−) sign structure for off-diagonal pairs.
Mixture of mutualism and competition. Each off-diagonal pair (M ij , M ji ) is sampled from a half-normal distribution jN ð0; sÞj with probability C/2 or ÀjN ð0; sÞj with probability C/2, and both entries are zero otherwise. These community matrices have a (+, +) or (−, −) sign structure for off-diagonal pairs.
Predator-prey. The first entry in an off-diagonal pair is sampled from a half-normal distribution jN ð0; sÞj and the second entry from ÀjN ð0; sÞj with probability C/2, or with the halfnormal distributions reversed with probability C/2, and both entries are zero otherwise. These community matrices have a (+, −) or (−, +) sign structure for off-diagonal pairs.

Pseudospectra and transient instability
In general, the eigenvalues of M satisfy the following definition: or, equivalently, meaning that if z is an eigenvalue of M then by convention the norm of (zI − M) −1 is defined to be infinity (see Chapter I.1 in [29]). The '-pseudospectrum' has several comparable definitions which describe the eigenvalues of a matrix whose entries have been subject to noise of magnitude (in the sense of the matrix norm) [28]. We use the following definition: If a matrix is normal then its -pseudospectrum (henceforth just 'pseudospectrum') consists of closed balls of radius surrounding the original eigenvalues of M (see Theorem 2.2 in [29]). As mentioned earlier, normal matrices cannot exhibit reactive dynamics: perturbations of the population vector for a stable system decay immediately and with exponential profile as the system returns to its original equilibrium point. But with non-normal matrices, pseudospectra can be much larger and more intricate and reactive dynamics are possible: perturbations of the population vector for a stable system first increase in magnitude and reach a maximum amplitude before eventually decaying (Fig 1). This behaviour motivates a description of local stability analysis for community matrices in terms of pseudospectra. (Besides non-normal matrices and reactivity, it is worth noting that pseudospectra are still relevant for understanding the consequences of small changes to entries in normal matrices.) Local asymptotic stability is determined in the same way for normal and non-normal matrices. The 'spectral abscissa' of M is defined as where the supremum (sup) selects for the largest (real-part) of the rightmost eigenvalue in the set Λ(M). Stability is guaranteed for α(M) < 0. If M is normal, then ||e Mt || = e α(M)t and dynamics are completely described by α(M) see Eq (4). Otherwise, the dynamics implied by M can be more complicated: where the columns of matrix V are the eigenvectors of M, and κ(V) = ||V|| Á ||V −1 || is known as the conditioning of V [32][33][34][35]. The conditioning provides a bound from above-an upper bound-to the maximum amplitude of a perturbation of the population vector (it is worth noting that κ(V) does not provide any information about the time at which the perturbation reaches its maximum amplitude). In complement to stability is reactivity, which describes the behaviour of a system close to t = 0, at the application of a perturbation. The 'numerical abscissa' of M is defined as where H ¼ MþM t 2 [13][14][15][16][17]. The numerical abscissa is the maximum initial amplification rate following an infinitesimally small perturbation to the population vector. Dynamics are non-reactive if ω(M) < 0 and may be reactive if ω(M) ! 0. A stable system can be either reactive or non-reactive, but an unstable system is necessarily reactive.
With non-normal matrices, perturbations to the entries of M can affect whether a system is stable and non-reactive or stable and reactive. In other words, perturbations to the entries of M can affect how a system responds to perturbations to the population vector. The effect of such perturbations to M is not covered by Eq (10). However, we can study the pseudospectrum of a community matrix to better understand system dynamics between the limits of reactivity and stability. In what follows, we use the theory of pseudospectra to relate uncertainty in, or small changes to, the entries of M to bounds on the amplification of perturbations of the population vector.
The '-pseudospectral abscissa' of M is defined as a ðMÞ ¼ sup and therefore the function is useful for understanding transient dynamics. Eqs (12) and (13) are also valid for bounding normal matrices with positive spectral abscissa. As ! 0, α (M) converges to the spectral abscissa. If M has a positive spectral abscissa, then lim !0 α (M)/ ! 1, which confirms that the norm is unbounded and the equilibrium point is unstable. In the literature on pseudospectra, sup !0 f M ðÞ KðMÞ is known as the Kreiss constant [32,34]. Eqs (11) (12) and (13) are useful because they relate perturbations to the matrix norm -small changes to the elements of the community matrix as described by the noise parameter -to the effect of perturbations to the population vector (compare Eqs (8) and (11)). For a given community matrix, as the size of a matrix perturbation is increased from zero there may be some critical value Ã at which f M (*) = 1. In the pseudospectrum, this is illustrated by the Ã -contour crossing the imaginary axis (Fig 1). At this point, perturbations to the equilibrium population vector begin to be amplified. For a stable and non-reactive system, perturbations to the population vector are not amplified and the system always returns to its original equilibrium point. For an unstable and necessarily reactive system, perturbations are amplified and the system may move to a new equilibrium point. But for a stable and reactive system, perturbations are first amplified before the system eventually returns to its original equilibrium point-this is transient instability. Now that we can compute upper Eq (9) and lower bounds Eq (12) for amplifications, we are in a position to compare the transient dynamics of different types of ecological community as described by non-normal community matrices.

Results
We generated multiple sets of community matrices with C = 0.1, μ = 1 and various combinations of S and σ for the five generative models. We first consider lower and upper bounds to the maximum amplitude of perturbations to the population vector for random community matrices, before turning our attention to the other types of interaction. The data required to reproduce the plots in this article are available at [36].

Lower bound for random community matrices
We numerically evaluated the -pseudospectral abscissa using the recently proposed subspace method [37]. Consider an ensemble of community matrices generated with random interaction type and S = 100 and σ = 0.3, which is just below the threshold for instability (s c ¼ m ffiffiffi ffi SC p ¼ 1 ffiffiffi 10 p % 0:31). We found that the average value of f M () Eq (13) monotonically increases as a function of and eventually saturates. It is worth noting that although the average value of f M () monotonically increases, the average value was calculated over 100 matrices so this may not be the case for f M () for a single matrix. This is for instance the case for the Monte Carlo simulations we have performed.
The key result of this paper is that at Ã % 0.085 the curve crosses one, at which point perturbations are amplified and transient instability may be observable. The function f M () converges for all asymptotically stable community matrices considered here.
In general, we identify regions of stability, transient instability and instability by plotting sup !0 a ðMÞ Eq (12); in practice, we plot f M () for large values of ) as σ is varied (Fig 2). Similar regions can be identified as S is varied while σ is held constant (results not shown). In the stable region, there is no perturbation to the community matrix large enough (that can still be considered infinitesimally small) such that sup !0 a ðMÞ > 1, and so perturbations are never amplified. At some critical point, σ ti , there is a level of matrix noise = Ã above which perturbations to the population vector are amplified before decaying. As σ increases in the region of transient instability, Ã decreases until it reaches zero at σ c . At this point, system dynamics are guaranteed to be asymptotically unstable and any infinitesimally small perturbation to the population vector is amplified (without necessarily returning to the original equilibrium point). In the unstable region, f M () diverges and corresponding values for the lower bound should be treated with caution.
The critical point for transient instability with S = 100 is σ ti % 0.22. This is very close to the value given by reactivity criteria based on the numerical abscissa: 18]. Indeed, both approaches determine whether perturbations to the population vector are amplified based on eigenvalues related to M. As a point of difference, however, the pseudospectral approach allows for an additional treatment of uncertainty in, or small changes to, entries of the community matrix. For a given set of parameters, the numerical abscissa only indicates whether amplification is possible, whereas the pseudospectrum, through the -pseudospectral abscissa, also indicates whether amplification is possible given small changes to the strengths of interactions among species in the community.

Upper bound for random community matrices
We plot the frequency distribution of κ(V) Eq (9) for various combinations of S and σ to investigate the upper bound to the maximum amplitude of perturbations of the population vector. In general, distributions are strongly peaked and fat-tailed (Fig 3). This indicates that very large amplification is possible even for very small perturbations. The location of the peak changes very little as σ increases, but shifts rightwards as S increases (results not shown). The slope of the tail does not change much as either S or σ is varied. With S = 100 and σ = σ c = 0.31, the peak in the distribution of upper bound values is UB peak (σ c ) % 95 and the maximum value in the tail is UB tail (σ c ) * 1000. When a power law is fit to the tail, f(x) / x − α , the exponent is α % 2.9.

Community matrices with different types of interaction
The region of transient instability varies for different types of interaction, as do lower and upper bounds for amplification (Table 1). Transient instability becomes observable with smallest σ ti with mutualism, followed by mixture, competition, random and predator-prey. This order is the same as with the threshold for instability, σ c . However, the size of the region of transient instability, σ c − σ ti , has a different order: predator-prey is largest, followed by random, mutualism, competition and mixture. The pattern is similar if S is varied while σ is held constant (results not shown). As expected, these findings are consistent with earlier results based on the numerical abscissa and the correlation between off-diagonal entries in a community matrix [18]. Predator-prey community matrices are relatively stable and exhibit the largest range of parameter values for transient instability. The lower bound to the maximum amplitude of perturbations of the population vector also reaches its largest value among the five types of interaction for predator-prey community matrices. However, the peak in the distribution of upper bounds is at lower amplification and the slope of the tail is steeper (Table 1). This implies that perturbations are typically amplified less severely compared to the other types of interaction and the very largest possible amplitudes are not as large.
Mutualism (+, +) and competition (−, −) have different critical points for transient instability and instability, but similar bounds to the maximum amplitude of perturbations of the population vector. Interestingly, the peak in the distribution of upper bounds is at lower amplification for community matrices with a mixture of these two interaction types. The largest upper bound, UB tail (σ c ), however, is similar to mutualism and competition, so the exponent α is shallower.

Discussion
Here we described transient instability for non-normal community matrices using local stability analysis and pseudospectra. We showed how the shift from stable and non-reactive dynamics to transient instability changes if perturbations are applied to the community matrix. We also characterised how perturbations of the population vector are amplified during periods of transient instability for different types of interaction. We found an early, sharp and severe transition between stability and instability with mutualism, mixture and competition, but a later, longer and less severe transition with predator-prey community matrices.
In this study, we assumed a random topology of interactions between species. Although the correlation between interaction strengths-and therefore the predominant type of interaction in a community matrix-may be more important than topology for stability [10,11], it remains to be seen whether this is the case with transient instability. Nevertheless, it is likely that the particular trajectory of a perturbed system is sensitive to topology, and, of course, the direction of initial perturbation of the population vector. Understanding transient dynamics at this level of detail requires analysis of pseudoeigenvectors in addition to pseudoeigenvalues (see Chapter I.4 in [29]).
Local stability analysis is only one approach to understanding the capacity for ecosystems to withstand external shocks [38,39]. It will be informative to compare how the time evolution of the same shock to the same system is assessed under different approaches to measuring the 'stability', 'persistence' or 'resilience' of ecosystems [12].
Stability, in principle, promises a degree of certainty that biodiversity will not be lost [1,2]. Reactivity has been suggested as a possible early-warning signal for the onset of instability [19][20][21][22][23]. Transient instability not only fills the gap between these two concepts, but also highlights new consequences of rapid environmental change. The longer the period of transient instability and the larger the amplification of perturbations of the population vector, the more susceptible an ecosystem is to multiple perturbations. One perturbation may drive a stable system into a period of transient instability that eventually dissipates; but two or three perturbations in quick succession may force the system to a new, unknown equilibrium point that may correspond to a loss of species and biodiversity. Pseudospectra can be used to investigate which ecosystems are at risk of instability, and what could be done to mitigate that risk.