Nonlinear analogue of the May−Wigner instability transition

Yan V. Fyodorov; Boris A. Khoruzhenko

doi:10.1073/pnas.1601136113

New Research In

Edited by Srinivasa S. R. Varadhan, Courant Institute of Mathematical Sciences, New York, NY, and approved May 9, 2016 (received for review January 22, 2016)

Significance

Complex systems equipped with stability feedback mechanisms become unstable to small displacements from equilibria as the complexity (as measured by the interaction strength and number of degrees of freedom) increases. This paper takes a global view on this transition from stability to instability. Our examination of a generic dynamical system whereby N degrees of freedom are coupled randomly shows that the phase portrait of complex multicomponent systems undergoes a sharp transition as the complexity increases. The transition manifests itself in the exponential explosion of the number of equilibria, and the transition threshold is quantitatively similar to the local instability threshold. Our model provides a mathematical framework for studying generic features of the global dynamics of large complex systems.

Abstract

We study a system of N≫1 degrees of freedom coupled via a smooth homogeneous Gaussian vector field with both gradient and divergence-free components. In the absence of coupling, the system is exponentially relaxing to an equilibrium with rate μ. We show that, while increasing the ratio of the coupling strength to the relaxation rate, the system experiences an abrupt transition from a topologically trivial phase portrait with a single equilibrium into a topologically nontrivial regime characterized by an exponential number of equilibria, the vast majority of which are expected to be unstable. It is suggested that this picture provides a global view on the nature of the May−Wigner instability transition originally discovered by local linear stability analysis.

Will diversity make a food chain more or less stable? The prevailing view in the midtwentieth century was that diverse ecosystems have greater resilience to recover from events displacing the system from equilibrium and hence are more stable. This “ecological intuition” was challenged by Robert May in 1972 (1). At that time, computer simulations suggested that large complex systems assembled at random might become unstable as the system complexity increases (2). May’s 1972 paper complemented that work with an analytic investigation of the neighborhood stability of a model ecosystem whereby N species at equilibrium are subject to random interactions.

The time evolution of large complex systems, of which model ecosystems is one example, is often described within the general mathematical framework of coupled first-order nonlinear ordinary differential equations (ODEs). In the context of generic systems, the Hartman−Grobner theorem then asserts that the neighborhood stability of a typical equilibrium can be studied by replacing the nonlinear interaction functions near the equilibrium with their linear approximations. It is along these lines that May suggested looking at the linear modeldyjdt=−μyj+∑k=1NJjk yk, j=1,…,N,[1]to study the stability of large complex systems. Here J=(Jjk) is the coupling matrix and μ>0. In the absence of interactions, i.e., when all Jjk=0, system [1] is self-regulating: If disturbed from the equilibrium y1=y2=…=yN=0, it returns back with some characteristic relaxation time set by μ. In an ecological context, yj(t) is interpreted as the variation about the equilibrium value, yj=0, in the population density of species j at time t. The element Jjk of the coupling matrix J, which is known as the community matrix in ecology, measures the per capita effect of species k on species j at the presumed equilibrium. Generically, the community matrix is asymmetric, Jjk≠Jkj.

For complex multispecies systems, information about the interaction between species is rarely available at the level of detail sufficient for the exact computation of the community matrix and a subsequent stability analysis. Instead, May considered an ensemble of community matrices J assembled at random, whereby the matrix elements Jjk are sampled from a probability distribution with zero mean and a prescribed variance α2. This is similar to the approach taken by Wigner in his description of statistics of energy levels of heavy nuclei via eigenvalues of random matrices, which proved to be very fruitful (3). A detailed review of May’s model in the light of recent advances in random matrix theory can be found in ref. 4.

The linear system [1] is stable if, and only if, all of the eigenvalues of J have real parts less than μ. Invoking Wigner’s arguments for studying eigenvalues of large random matrices, May claimed that, for large N, the largest real part of the eigenvalues of J is typically αN. Obviously, the model’s stability is then controlled by the ratio m=μ/(αN). For large N, system [1] will almost certainly be stable if m>1 and unstable if m<1, with a sharp transition between the two types of behavior with changing either μ, α, or N. In particular, for fixed μ,α, system [1] will almost certainly become unstable for sufficiently large N.

Despite the simplistic character of May’s model (5), his pioneering work gave rise to a long-standing “stability versus diversity” debate, which is not fully settled yet (4, 6 ⇓⇓–9), and played a fundamental role in theoretical ecology by prompting ecologists to think about special features of real multispecies ecosystems that help such systems to remain stable. Variations of May’s model are still being discussed nowadays in the context of neighborhood stability (see ref. 4 and references therein).

One obvious limitation of the neighborhood stability analysis is that it gives no insight into the model behavior beyond the instability threshold. Hence May’s model has only limited bearing on the dynamics of populations operating out of equilibrium. An instability does not necessarily imply lack of persistence: Populations could coexist thanks to limit cycles or chaotic attractors, which typically originate from unstable equilibrium points. Important questions posing serious challenges then relate to classification of equilibria by stability, studying their basins of attraction, and other features of global dynamics (4). Over the last few years, as computing power grew, nonlinear models have increasingly been used to investigate population dynamics on the global scale by means of numerical integration of the corresponding system of ODEs (10 ⇓⇓–13). Although such investigations captured a rich variety of types of behavior such as fold bifurcations when points of equilibrium merge/annihilate (14), or various types of chaotic dynamics (15, 16), they provide little analytic insight and are limited to small to medium-sized systems.

In this paper, we attempt to investigate the generic properties of the global dynamics of large complex multispecies systems by retaining in our model only the bare essentials—nonlinearity and stochasticity. Much in the spirit of May’s original approach, the model we propose is simple enough to allow for an analytic treatment yet, at the same time, is rich enough to exhibit a nontrivial behavior. In particular, our model captures an instability transition of the May−Wigner type, but now on the global scale. It also sheds additional light on the nature of this transition by relating it to an exponential explosion in the number of equilibria. Interestingly, despite the nonlinear setting of the problem, the random matrix ideas again play a central role in our analysis.

Similar to the May’s linear model, our toy model is likely to have rather limited practical significance for quantitative description of real ecosystems, but it might provide insight into the generic qualitative features of such systems and beyond, e.g., machine learning (17) or financial ecosystems (18, 19). The idea of destabilization by interaction is of relevance far beyond the mathematical ecology context, as applications of systems of many coupled nonlinear ODEs are vast [e.g., complex gene regulatory networks (20), neural networks (21, 22), or random catalytic reaction networks (23)].

Model

Consider a system of N coupled nonlinear autonomous ODEs of the formdxidt=−μxi+fi(x1,…,xN), i=1,…,N,[2]where μ>0 and the components fi(x) of the vector field f=(f1,…,fN) are zero mean random functions of the state vector x=(x1,…,xN). To put this model in the context of the discussion in the Introduction, if xe is an equilibrium of [2], i.e., if −μxe+f(xe)=0, then, in the immediate neighborhood of xe, system [2] reduces to May’s model [1] with y=x−xe and Jjk=(∂fj/∂xk)(xe).

The nonlinear system [2] may have multiple equilibria whose number and locations depend on the realization of the random field f(x). To visualize the global picture, it is helpful to consider first a special case of a gradient descent flow, characterized by the existence of a potential function V(x) such that f=−∇V. In this case, system [2] can be rewritten as dx/dt=−∇L, with L(x)=μ|x|2/2+V(x) being the associated Lyapunov function describing the effective landscape. In the domain of L, the state vector x(t) moves in the direction of the steepest descent, i.e., perpendicular to the level surfaces L(x)=h toward ever-smaller values of h. This provides a useful geometric intuition. The term μ|x|2/2 represents the globally confining parabolic potential, i.e., a deep well on the surface of L(x), which does not allow x to escape to infinity. At the same time, the random potential V(x) may generate many local minima of L(x) (shallow wells), which will play the role of attractors for our dynamical system. Moreover, if the confining term is strong enough, then the full landscape will only be a small perturbation of the parabolic well, typically with a single stable equilibrium located close to x=0. In the opposite case of a relatively weak confining term, the disorder-dominated landscape will be characterized by a complicated random topology with many points of equilibria, both stable and unstable. Note that, in physics, complicated energy landscapes are a generic feature of glassy systems with intriguingly slow long-time relaxation and nonequilibrium dynamics (see, e.g., ref. 24).

The above picture of a gradient descent flow is, however, only a very special case, because the generic systems of ODEs [2] are not gradient. The latter point can easily be understood in the context of model ecosystems. For, by linearizing a gradient flow in a vicinity of any equilibrium, one always obtains a symmetric community matrix, whereas the community matrices of model ecosystems are, in general, asymmetric. Note also a discussion of an interplay between nongradient dynamics in random environment and glassy behavior in ref. 25.

To allow for a suitable level of generality, we therefore suggest choosing the N-dimensional vector field f(x) as a sum of gradient and nongradient (solenoidal) contributions,fi(x)=−∂V(x)∂xi+1N∑j=1N∂Aij(x)∂xj, i=1,…,N,[3]where we require the matrix A(x) to be antisymmetric: Aij=−Aji. The meaning of this decomposition is that vector fields can be generically divided into a conservative irrotational component, sometimes called “longitudinal,” whose gradient connects the attractors or repellers, and a solenoidal curl field, also called “transversal.” As discussed in, e.g., ref. 26, such a representation is closely related to the so-called Hodge decomposition of differential forms and generalizes the well-known Helmholtz decomposition of the three-dimensional vector fields into curl-free and divergence-free parts to higher dimensions. Correspondingly, we will call V(x) the scalar potential and the matrix A(x) the vector potential. The normalizing factor 1/N in front of the sum on the right-hand side in Eq. 3 ensures that the transversal and longitudinal parts of f(x) are of the same order of magnitude for large N.

Finally, to make the model as simple as possible and amenable to a rigorous and detailed mathematical analysis, we choose the scalar potential V(x) and the components Aij(x), i<j, of the vector potential to be statistically independent, zero mean Gaussian random fields, with smooth realizations and the additional assumptions of homogeneity (translational invariance) and isotropy reflected in the covariance structure,〈V(x)V(y)〉=υ2ΓV(|x−y|2);[4]〈Aij(x)Anm(y)〉=a2ΓA(|x−y|2)(δinδjm−δimδjn).[5]Here the angular brackets 〈…〉 stand for the ensemble average over all realizations of V(x) and A(x), and δin is the Kronecker delta: δin=1 if i=n and zero otherwise.

For simplicity, we also assume that the functions ΓV(r) and ΓA(r) do not depend on N. This implies (27)Γσ(r)=∫0∞exp(−sr)γσ(s)ds, σ=A,V,where the radial spectral densities γσ(s)≥0 have finite total mass: ∫0∞γσ(s)ds<∞. We normalize these densities by requiring that Γσ′′(0)=∫0∞s2γσ(s)ds=1. The ratioτ=v2/(v2+a2), 0≤τ≤1,is a dimensionless measure of the relative strengths of the longitudinal and transversal components of f(x): If τ=0, then f(x) is divergence-free, and, if τ=1, it is curl-free.

Results

Determining and classifying all points of equilibria of a dynamical system with many degrees of freedom is a well-known formidable analytical and computational problem. In this paper, we shall focus our investigation on the simplest, yet informative, characteristic of system [2] by counting its total number of equilibria, that is, the total number Ntot of solutions of the simultaneous equations−μxi+fi(x1,…,xN)=0,i=1,…,N.[6]Certainly, finding Ntot is a good starting point of any phase portrait analysis.

Had we restricted ourselves to the gradient descent flows, Ntot would simply count the number of stationary points (minima, maxima, or saddle points) on the surface of the Lyapunov function L(x). The problem of counting and classifying stationary points of high-dimensional random energy landscapes of various types has attracted considerable interest in recent years (28 ⇓⇓⇓⇓–33). In particular, works (28, 29) study such energy landscapes generated by a potential equivalent to the above Lyapunov function. One of the main conclusions of that study is that, for large N, the topology of the Lyapunov function changes drastically with decrease of the strength of the confining term relative to that of the interaction term in L(x). The change manifests itself in the emergence of a multitude of equilibria, exponential in number. Such a transition is intimately connected to the spin-glass-like restructuring of the Boltzmann−Gibbs measure induced by the Lyapunov function when the latter is treated as an effective energy landscape.

We shall prove below that, for large N, the general autonomous system defined by Eqs. 2 and 3 exhibits a similar drastic change in the total number of equilibria when the control parameterm=μαN, where α=2υ2+a2,drops below the threshold value mc=1. As in the case of gradient systems, the proof involves the Kac−Rice formula as a starting point. However, performing the subsequent steps requires quite different mathematical techniques due to the asymmetry of the Jacobian matrix for nongradient systems.

The Kac−Rice formula (see, e.g., ref. 34) counts solutions of simultaneous algebraic equations. Under our assumptions (homogeneity, isotropy, and Gaussianity of V and A), this formula yields the ensemble average of Ntot in terms of that of the modulus of the spectral determinant of the Jacobian matrix (Jij)i,j=1N, Jij=∂fi/∂xj (see Materials and Methods),〈Ntot〉=1μN〈|det(−μδij+Jij)|〉,[7]thus bringing the original nonlinear problem into the realms of the random matrix theory.

The probability (ensemble) distribution of the matrix J can easily be determined in closed form. Indeed, the matrix entries of J are zero mean Gaussian variables and their covariance structure, at spatial point x, can be obtained from Eqs. 4 and 5 by differentiation,〈JijJnm〉=α2[(1+ϵN)δinδjm+(τ−ϵN)(δjnδim+δijδmn)],where ϵN=(1−τ)/N. Thus, to leading order in the limit N→∞,Jij=α(Xij+τδijξ),[8]where Xij, i,j=1,…,N are zero mean Gaussians with〈XijXnm〉=δinδjm+τδjnδim,[9]and ξ is a standard Gaussian, ξ∼N(0,1), which is statistically independent of X=(Xij). Note that, for the divergence free fields f(x) (i.e., if τ=0), the entries of J are statistically independent in the limit N→∞, exactly as in May’s model. On the other side, if f(x) has a longitudinal component (τ>0), then this implies positive correlation between the pairs of matrix entries of J symmetric about the main diagonal: 〈XijXji〉=τ if i≠j. Such distributions of the community matrix have also been used in the neighborhood stability analysis of model ecosystems (8). Finally, in the limiting case of curl-free fields (τ=1), the matrix J is real symmetric.

Representation [8] comes in handy, as it allows one to express [7] as a random matrix integral,〈Ntot〉=N−N2mN∫−∞∞〈|det(x δij−Xij)|〉XNe−Nt22dt2π/N,[10]where x=N(m+tτ) and the angle brackets 〈…〉XN stand for averaging over the real elliptic ensemble of random N×N matrices X defined in Eq. 9; see also Eq. 20. This one-parameter family of random matrices interpolates between the Gaussian Orthogonal Ensemble of real symmetric matrices (GOE, τ=1) and real Ginibre ensemble of fully asymmetric matrices (rGinE, τ=0) (see ref. 35 for discussions). Both rGinE and its one-parameter extension defined in Eq. 9 have enjoyed considerable interest in the literature in recent years (36 ⇓⇓⇓–40).

The matrix X is asymmetric (unless τ=1) and can have real as well as complex eigenvalues. The latter come in complex conjugate pairs. Their density, in the limit N→∞, is constant inside the ellipse with the main half-axis N(1±τ) and vanishes sharply outside of the ellipse (35, 39, 41). The corresponding theorem is known as the Elliptic Law, and its validity extends beyond the Gaussian matrix distributions (42, 43). However, in the context of our investigation, it is the density of real eigenvalues of X that appears to be most relevant.

Denote by ρN(r)(x) the density of real eigenvalues of N×N matrices X [9] averaged over all realizations of X. It is convenient to normalize ρN(r)(x) in such a way that ∫αβρN(x) dx gives the average number of real eigenvalues of X in the interval [α,β]. A crucial observation is that ρN(r)(x) is directly related to the averaged value of the modulus of the determinant that appears in Eq. 10. Namely,〈|det(xδij−Xij)|〉XN=CN(τ)ex22(1+τ)ρN+1(r)(x),[11]where CN(τ)=21+τ (N−1)!/(N−2)!! and ρN+1(r)(x) is the average density of real eigenvalues of matrices X of size (N+1)×(N+1). For the limiting case τ=0, this relation appeared originally in ref. 44, and it can be extended to any τ∈[0,1) without much difficulty (see Supporting Information for a derivation of Eq. 11 following the approach of ref. 45). In the limiting case of real symmetric matrices τ=1, all eigenvalues of X are real and relation [11] is also valid (28).

Combining Eqs. 10 and 11 and changing the variable of integration from t to λ=m+tτ, one can express 〈Ntot〉 for system [2] with N degrees of freedom in terms of the density of real eigenvalues in the elliptic ensemble of random matrices [9] of size (N+1)×(N+1),〈Ntot〉=KN(τ)mN∫−∞∞e−NS(λ)ρN+1(λN) dλ2π,[12]where S(λ)=(λ−m)2/(2τ)−λ2/[2(1+τ)] and KN(τ)=N−N+12CN(τ)/τ. The importance of this relation is due to the fact that ρN(r)(x) is known, in closed form, in terms of Hermite polynomials (39). This allows us to carry out an asymptotic evaluation of the integral in [12] and calculate 〈Ntot〉 in the limit N→∞. The key finding that emerges from this calculation is that 〈Ntot〉 changes drastically around m=1. If m>1, thenlimN→∞〈Ntot〉=1.[13]On the other hand, if 0<m<1, then, to leading order in the limit N→∞,〈Ntot〉=γτeN∑tot(m),[14]where ∑tot(m)=(1/2)(m2−1)−ln⁡m>0 for all 0<m<1. Therefore, if m<1, then 〈Ntot〉 grows exponentially with N. The factor in front of the exponential in Eq. 14 is given by γτ=2(1+τ)/(1−τ) as long as τ<1. The gradient limit τ=1 can be approached by scaling τ with N. Setting τ=1−u2/N,0≤u<∞, one obtains γτ=4N/π∫01−m2e−u2p2dp. This regime describes a weakly nongradient flow. The corresponding regime for ensembles of asymmetric matrices was discovered long ago (46, 47).

Close to the phase transition point m=1, the complexity exponent vanishes quadratically, ∑tot=(1−m)2 as m→1, implying that the width of the transition region around m=1 is 1/N. According to the general lore of phase transitions, for large but finite N, there exists a “critical regime” m=1+κN−1/2 where the number of equilibria changes smoothly between the two phases [13] and [14]. A quick inspection of Eq. 12 shows that the corresponding crossover profile is determined by the profile of ρN(r)(x) in the vicinity of the “spectral edge” x=(1+τ)N (see Materials and Methods). After rescaling λ, λ=1+τ+ζ1−τ2/N, the density ρN(r)(λN) converges to (1/1−τ2)ρedge(r)(ζ) in the limit N→∞, where (39)ρedge(r)(ζ)=122π{erfc(2ζ)+12e−ζ2[1+erf(ζ)]},[15]with erf(x)=1−erfc(x)=(2/π)∫0xe−t2 dt. In terms of ρedge(r)(ζ), the critical crossover profile is given by〈Ntot〉=γτeκ2∫−∞∞e−t2/2 ρedge(r)(cτt+κγτ/2)dt,[16]where cτ=τ/(1−τ). The right-hand side of [16] interpolates smoothly between the two regimes [13] and [14], when parameter κ runs from κ=−∞ to κ=+∞.

Although our investigation is concerned with the ensemble average of the number of equilibria Ntot, we expect that, in the limit N→∞, the deviations of Ntot from its average 〈Ntot〉 are relatively small. This is certainly the case above the critical threshold, for m>1. For, under some additional technical assumptions on the decay of correlations for f(x), system [2] will almost certainly have at least one stationary point (see ref. 48 for the relevant results about the maxima of homogeneous Gaussian fields). Therefore, Ntot≥1, and the established convergence of 〈Ntot〉 to 1 in the limit N→∞ actually implies that the probability for Ntot to take other values than 1 is close to zero for large N. The problem of estimating the deviation of Ntot from its average value in the opposite regime 0<m<1 is much harder and is an open and interesting question.^†

Discussion

In this paper, we introduced a model describing generic large complex systems and examined the dependence of the total number of equilibria in such systems on the system complexity as measured by the number of degrees of freedom and the interaction strength. The inspiration for our work came from May’s pioneering study (1) of the neighborhood stability of large model ecosystems. Our outlook is complementary to that of May’s in that it adopts a global point of view, which is not limited to the neighborhood of the presumed equilibrium.

In the context of model ecosystems, our analysis is applicable to complex multispecies communities in which each kind of species on its own becomes extinct and thus interaction between species is key to persistence of the community. The key feature of our analysis is that, in the presence of interactions, as the complexity increases, there is an abrupt change from a simple set of equilibria (typically, a single equilibrium for large number of species N≫1) to a complex set of equilibria, with their total number growing exponentially with N. In the latter regime, we expect the stable equilibria to be only a tiny proportion of all of the multitude of equilibria (see discussion below), which is indicative of long relaxation times and transient nonequilibrium behavior.

We expect this sharp transition in the phase portrait to be shared by other systems of randomly coupled autonomous ODEs with large numbers of degrees of freedom. To that end, it is appropriate to mention that, very recently, a similar explosion in complexity was reported in a model of a neural network consisting of randomly interconnected neural units (22). The model considered in ref. 22 is essentially of form [2] but with the particular choice of fi=∑jJijS(xj) where S is an odd sigmoid function representing the synaptic nonlinearity and Jij are independent centered Gaussian variables representing the synaptic connectivity between neuron i and j. Although Gaussian, the corresponding (nongradient) vector field is not homogeneous and thus seems rather different from our choice and not easily amenable to a rigorous analysis. Nevertheless, a shrewd semiheuristic analysis of ref. 22 revealed that, close to the critical coupling threshold, the two models actually display very similar behavior, described essentially by the same exponential growth in the total number of equilibria with rate ∑tot(m). This fact points toward considerable universality of the transition from [13] to [14] and suggests that the crossover function [16] may be universal as well.

Our model captures an abrupt change in the dynamics of large complex systems on the macroscopic scale. At the same time, zooming in to classify each and every equilibrium point into locally stable or unstable seems a hard task. For, although linearizing the field f(x) around a given equilibrium is fairly straightforward, with the outcome being the Jacobian matrix [8], conditioning by the positions of equilibria and taking into account all eventualities is a highly nontrivial task. Given the stochastic setup of our model, the question about stability of individual equilibria may even be the wrong question to ask, whereas addressing the statistics of the number of stable equilibria seems very appropriate.

Arguments similar to those in Results yield the ensemble average of the total number of stable equilibria, 〈Nst〉, over all realizations of the vector field f(x) in terms of the random matrix integral (compare with Eq. 10),〈Nst〉=N−N2mN∫−∞∞〈det(x δij−Xij)χx(Xij)〉XNe−Nt22dt2π/N,[17]where χx(X)=1, if all N eigenvalues of matrix X have real parts less than the spectral parameter x=N(m+tτ), and χx(X)=0 otherwise. In the limiting case of a purely gradient dynamics τ=1, the rescaled Jacobian matrix X is real symmetric with all N eigenvalues real. In this case, the above integral can be related to the probability density of the maximal eigenvalue of the GOE matrix (29, 30), with the latter being a well-studied object in the random matrix theory (see, e.g., ref. 50 and references therein). This observation can then be used to evaluate 〈Nst〉 for N≫1. One finds (29) that 〈Nst〉→1 if m>1, whereas, if 0<m<1, then, to leading order in N, 〈Nst〉∝eNΣst, with 0<Σst<Σtot. Thus, in the case of purely gradient dynamics, as the complexity increases, large nonlinear autonomous systems assembled at random undergo an abrupt change from a typical phase portrait with a single stable equilibrium to a phase portrait dominated by an exponential number of unstable equilibria with an admixture of a smaller, but still exponential in N, number of stable equilibria.

It was suggested to us by J.-P. Bouchaud that, in the general case of nongradient dynamics 0≤τ<1, it would be natural to expect a further phase transition in the plane (m,τ) such that, below a certain number τc(m), stable equilibria are no longer exponentially abundant in the limit N→∞ [i.e., ∑st(m,τ)→0], with further implications for the global dynamics. Unfortunately, for a fixed 0≤τ<1, only a vanishing fraction of order N−1/2 of eigenvalues of X remain real, and the relation of the integral in Eq. 17 to statistics of the largest real eigenvalue in the elliptic ensemble seems to be lost. This fact has prevented us, so far, from reliable counting of stable equilibria for the general case of nongradient flows. In principle, for given values of parameters N,τ,m, one may attempt to evaluate the random matrix ensemble average in the integral in [17] numerically, and then evaluate numerically the integral itself. Although such a procedure seems tractable, its actual implementation with sufficient precision is not straightforward, especially in the limit N→∞, due to the exponentially large values involved. Clarification of the status of the picture suggested by J.-P. Bouchaud and related studies remain an important outstanding issue.

Materials and Methods

Kac−Rice Formula.

The expected number 〈Ntot〉 of simultaneous solutions to the system of Eq. 6 in ℝN is given by the formula (see, e.g., ref. 34)〈Ntot〉=∫ℝN〈δ(−μx+f(x))|det(−μδij+Jij(x))|〉 dxN,[18]

where δ(x) is the multivariate Dirac δ-function, dxN is the volume element in ℝN, and Jij(x)=∂fi/∂xj are matrix elements of the Jacobian matrix J=(Jij). By our assumptions, the random field f(x) is homogeneous and isotropic. For such fields, samples of f and J taken in one and the same spatial point x are uncorrelated, 〈fl⋅∂fi/∂xj〉=0 for all i,j,l. This is well known and can be checked by straightforward differentiation. In addition, the field f is Gaussian; hence the f(x) and J(x) are actually statistically independent. This simplifies the evaluation of the integral in [18] considerably. Indeed, the statistical average in [18] factorizes and, because 〈|det(−μδij+Jij(x))|〉 does not vary with x, the integrand effectively reduces to〈δ((−μx+f(x))〉=∫dkN(2π)N e−μk⋅x〈eik⋅f(x)〉.[19]

Furthermore, at every spatial point x, the vector f(x) is Gaussian with uncorrelated and identically distributed components,〈fi(x)fj(x)〉=δijσ2, σ2=2v2|Γ′V(0)|+2a2|Γ′A(0)|N−1N.

Therefore, 〈eik⋅f(x)〉=e−σ2|k|2/2, and, after evaluating the integral on the right-hand side in [19], one arrives at [7].

Real Elliptic Matrices and Asymptotics of 〈N〉tot.

The joint probability density (JPD) function PN(X) of the matrix entries in the elliptic ensemble of real Gaussian random matrices X of size N×N is given byPN(X)=ZN−1⁡exp[−12(1−τ2)Tr(XXT−τX2)],[20]

where ZN is the normalization constant and τ∈[0,1). It is straightforward to verify that the covariance of matrix entries Xij is given by the expression specified in [9]. The mean density of real eigenvalues of ρN(r)(x) in the elliptic ensemble [20] is known, in closed form, in terms of Hermite polynomials (see ref. 39). Assuming, for simplicity, that N+1 is even, one has ρN+1(r)(x)=ρN+1(r),1(x)+ρN+1(r),2(x), whereρN+1(r),1(x)=12π∑k=0N−1 |ψk(τ)(x)|2k![21]

andρN+1(r),2(x)=12π(1+τ)(N−1)!ψN(τ)(x)∫0xψN−1(τ)(u)du.[22]

Here ψk(τ)(x)=e−x2/[2(1+τ)]hk(τ)(x) and hk(τ)(x) are rescaled Hermite polynomials, hk(τ)(x)=(1/π)∫−∞∞e−t2(x+it2τ)k dt. This, together with Eq. 12, allows one to evaluate 〈N〉tot in the limit N→∞. We shall sketch the corresponding evaluation below.

The asymptotics of ρN(r)(x) in the bulk and at the edge of the support of the distribution of real eigenvalues in the real elliptic ensemble were found in ref. 39, and, outside of the support, it can also be readily extracted using [21] and [22]. In particular, in the bulk, i.e., for |x|<(1+τ)N, the contribution of [21] to ρN(r)(x) is dominant, and, to leading order in N,ρN+1(r)(λN)||λ|<1+τ=12π(1−τ2).[23]

At the same time, for |x|>(1+τ)N, both [21] and [22] yield exponentially small contributions to ρN+1(r)(x), with [22] being dominant. Our evaluation yieldsρN+1(r)(λN)|λ>(1+τ)=Q(λ)exp[−NΨ(λ)], whereQ(λ)=[τ2π(1+τ)1λ2−4τ)(λ−λ2−4τ)]1/2,[24]Ψ(λ)=λ22(1+τ)−18τ(λ−λ2−4τ)2−lnλ+λ2−4τ2τ.[25]

The form of [12] suggests the application of the Laplace method. One easily finds that S(λ) has a minimum at λ*=m(1+τ), which belongs to the domain |λ|<1+τ as long as 0<m<1. Thus, applying the Laplace method and taking into account the asymptotic formulaKN(τ)≈2πNe−N/2(1+τ)/τ,[26]

one arrives at the asymptotic expression [14] for 〈N〉tot in the parameter range 0<m<1. For m>1, the saddle point occurs in the domain λ>1+τ so that the analysis requires search for the minimum of S(λ)+Ψ(λ). After straighforward algebra, we find (d/dλ)(S(λ)+Ψ(λ))=(1/2τ)(λ+λ2−4τ)−(m/τ), which is equal to zero at λ=λ*=m+(τ/m)>1+τ. One also verifies that this is a point of minimum for S(λ)+Ψ(λ), and a further simple calculation then yields S(λ*)+Ψ(λ*)=−ln(m/τ). Calculating the saddle point contribution then yields [13].

The above asymptotic analysis assumes that 0≤τ<1. Let us now discuss the modifications required to study the scaling regime of weakly nongradient flow τ→1 for 0<m<1. We only need to evaluate the leading contribution to ρN+1(r) given by [21]. By making use of the above integral representation for the scaled Hermite polynomials hk(τ)(x) and applying the scaling τ=1−u2/N, we can writeρN+1(r),1(λN)=12πN2πτe−Nλ21+τ+Nλ2τIN+1(λ),

with IN(λ) given byIN(λ)=πNe−Nλ22∫−∞∞e−u2p22ΦN−2(N2(p2+λ22)),

where ΦN(a)=e−a∑k=0Nak/k!. Recalling that the limit of ΦN(a) as N→∞ is 1 if 0<a<1 and 0 if a>1, one obtainsρN+1(r)(λN)=1πN∫01−λ2/4e−u2p2dp.

Substituting this expression into the integrand in [12] and evaluating the integral in the limit N→∞ (hence, τ→1) by the Laplace method then yields 〈N〉tot in the weakly nongradient regime.

Finally, our calculation of the profile of 〈N〉tot in the transitional region m=1+κN−1/2 uses the fact that, in such a regime, the main contribution to the integral [12] comes from the neighborhood of the spectral edge, λ=1+τ+(ζ1−τ2)/N, where we have, to the leading order in N,e−N[S(λ)+12]mN=exp[−1−τ2τ(κ2+ζ2)+1−τ2τκζ].

This, together with [26] and [15], converts [12] to [16].

SI Text

In this supporting information, we express the mean density of real eigenvalues in the real elliptic ensemble (see Eq. S9) of N×N matrices XN in terms of the modulus of the spectral determinant of real elliptic matrices of size (N−1)×(N−1). Our derivation is based on the method suggested in ref. 44; however, our Jacobian computation differs from the one given in ref. 44 and may be of interest on its own. For a similar calculation in the context of complex matrices, see ref. 45.

Housholder Reflections and Partial Triangularization of Real Matrices

The key idea of ref. 44 is based on using Householder reflections described by matricesPv=IN−2v⊗vT, vTv=1,[S1]where IN is the N×N identity matrix, and v is a column vector in ℝN of unit length. Its transpose, vT, is a row vector, so that the Kronecker product v⊗vT is a matrix. Pv describes the reflection about the hyperplane with normal v and passing through the origin. Obviously, Pv is symmetric and orthogonal: PvT=Pv and Pv2=IN.

Let e1 be the first vector of the standard basis in ℝN, i.e., e1T=(1,0,…,0). For any unit vector x,|x|=1, definev=x+e12(1+x1)[S2]and consider the Hausholder reflection Pv built from the above v according to Eq. S1. Then it is easy to check that Pvx=−e1. This implies that, for any nonzero vector x, there exists a Housholder reflection such that Pvx=ke1, with k=−|x|.

Let λ be a real eigenvalue of the N×N matrix A(N) with real entries Aij(N), i.e., A(N)x=λx for some column N-vector x of unit length. Our goal is to demonstrate that it is always possible to represent that matrix asA(N)=P(λwT0A(N−1))PT[S3]for some real (N−1)-component vector w and a real matrix A(N−1) of size (N−1)×(N−1). To verify this, take the Housholder reflection acting on the eigenvector x as Pvx=ke1. Obviously, PvA(N)x=λPvx=λke1. In view of Pv2=IN, we can rewrite the left-hand side as PvA(N)PvPvx=kPvA(N)Pve1, which, after denoting A˜=PvA(N)Pv, imples A˜e1=λe1. Using the definition of e1, we see that A˜11=λ and A˜1j=0 for all j=2,…,N. Therefore, A˜ can be written as A˜=(λw0A(N−1)), and, as A(N)= PvA˜Pv, the relation [S3] follows.

Considering now the volume element dA(N)=∏i,jNdAij(N), our next goal is to write it down in terms of N2 independent variables parametrizing the right-hand side of [S3], that is, (N−1)2 variables parametrizing A(N−1), N−1 components of w, one parameter for λ, and the remaining N−1 parameters for representing the matrix P. A convenient parametrization for P comes from using [S2] for the vector v, which shows that the last N−1 components of that vector (v2,…,vN)T≡q can be used as independent variables, whereas normalization fixes the first component. Writing vT=(1−qTq,q) with |q|<1 and using [S1] yields an explicit parametrization,P=(2qTq−1−2qT1−qTq−2q1−qTqIN−1−2q⊗qT).[S4]The problem therefore amounts to calculating the Jacobian of the transformation A(N)→(λ,w,q,A(N−1)). To that end, we start with differentiating [S3], which givesdA(N)=P{[(PdP),(λwT0A(N−1))]+(dλdwT0dA(N−1))}PT,where we made use of dP P=−P dP and also used the notation [A,B]=AB−BA for the matrix commutator. A direct calculation using [S4] shows that the matrix (PdP) can be symbolically written as(PdP)=(0−dbTdbdF),[S5]where the expression for dF is immaterial for our goals, anddb=(−21−qTq+1−2qTq1−qTqq⊗qT)dq+(−41−qTq+1−2qTq1−qTq)(dqTq)q.[S6]Substituting [S5] into the expression for dA(N), we find thatdA(N)=P(dλ−wTdbdwT+dbT(λ−A(N−1))−wT dF(λ−A(N−1)) dbdA(N−1)+db⊗wT+[dF,A(N−1)])PT.From this expression, we easily read off the required Jacobian to be given byJacobian=|det(λIN−1−A(N−1))| |∂b∂q|,[S7]where the last factor symbolically denotes the part of the Jacobian coming from the transformation db→dq described in [S6]. A straightforward calculation shows that|∂b∂q|=2N(1−qTq)N2−1,so that finally we arrive at the change-of-variables formuladA(N)=2N(1−qTq)N2−1|det(λIN−1−A(N−1))|dλ dwN−1 dqN−1 dA(N−1).[S8]

Elliptic Ensemble of Gaussian Random Matrices

The JPD of the elliptic ensemble of Gaussian random matrices XN of size N×N whose entries have covariance 〈XijXnm〉=δinδjm+τδjnδim is given byP(XN)=ZN−1⁡exp−12(1−τ2)[Tr(XNXNT)−τ Tr(XN2)],[S9]where ZN is the corresponding normalization factorZN=2N/2πN(N+1)/2(1+τ)N(N+1)/4(1−τ)N(N−1)/4.Our goal is to find the JPD of variables λ,w,q,XN−1 used to perform the partial triangulation of XN via [S3], with the role of A(N) played now by XN. We haveTr(XNXNT)=λ2+wTw+Tr(XN−1XN−1T),Tr(XN2)=λ2+Tr(XN−12).Taking into account the change-of-variables formula S8, we see that the corresponding JPD can be written asP˜(λ,w,q,XN−1)=2N(1−qTq)N2−1e−1(1−τ2)wTwZNZN−1e−12(1+τ)λ2|det(λIN−1−XN−1)|P(XN−1).By definition, the density of real eigenvalues ρN(r)(λ) is obtained by integrating the above JPD over variables |q|<1,−∞<w<∞ and finally over XN−1. After performing the integrals over q and w, we immediately arrive at the relationρN(r)(λ)=1CN−1(τ)e−12(1+τ)λ2〈|det(λIN−1−XN−1)|〉XN−1,where CN−1(τ) is a certain constant that can be found explicitly. This is precisely equivalent to Eq. 11, with the obvious change of notations: λ→x and N→N+1.

Acknowledgments

The authors thank David Arrowsmith, Matthew Evans, and Vito Latora for comments on various versions of this paper and J.-P. Bouchaud for his comments on the number of stable equilibria in nongradient systems. The research of Y.V.F. was supported by EPSRC Grant EP/J002763/1, “Insights into Disordered Landscapes via Random Matrix Theory and Statistical Mechanics.” B.A.K. thanks the Isaac Newton Institute for Mathematical Sciences for its hospitality during the program Periodic and Ergodic Spectral Problems, supported by EPSRC Grant EP/K032208/1.

Footnotes

↵¹To whom correspondence may be addressed. Email: y.fyodorov{at}qmul.ac.uk or b.khoruzhenko{at}qmul.ac.uk.

Author contributions: Y.V.F. and B.A.K. performed research and wrote the paper.
The authors declare no conflict of interest.
↵^†In this context, we would like to mention the recent work of Subag (33), who proved that the deviations of Ntot from 〈Ntot〉 in the spherical p-spin model are negligible in the limit of large system size. Although that model is different from ours, it is not dissimilar to the gradient limit of τ=1 of our model (32); for instance, the average number of equilibria grows exponentially with N (30). Thus one might hope to adopt the technique of ref. 33 to our model. Another relevant reference is ref. 49.
This article is a PNAS Direct Submission.
This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1601136113/-/DCSupplemental.

References

↵
1. May RM
(1972) Will a large complex system be stable? Nature 238(5364):413–414.
.
OpenUrl CrossRef PubMed
↵
1. Gardner MR,
2. Ashby WR
(1970) Connectance of large dynamic (cybernetic) systems: Critical values for stability. Nature 228(5273):784.
.
OpenUrl CrossRef PubMed
↵
1. Porter CE
(1965) Statistical Theories of Spectra: Fluctuations. A Collection of Reprints and Original Papers (Academic, New York).
.
↵
1. Allesina S,
2. Si T
(2015) The stability complexity relationship at age 40: A random matrix perspective. Popul Ecol 57(1):63–75.
.
OpenUrl CrossRef
↵
1. James A, et al.
(2015) Constructing random matrices to represent real ecosystems. Am Nat 185(5):680–692.
.
OpenUrl CrossRef PubMed
↵
1. McCann KS
(2000) The diversity-stability debate. Nature 405(6783):228–233.
.
OpenUrl CrossRef PubMed
↵
1. Boyd IL
(2012) The art of ecological modeling. Science 337(6092):306–307.
.
OpenUrl Abstract/FREE Full Text
↵
1. Allesina S,
2. Tang S
(2012) Stability criteria for complex ecosystems. Nature 483(7388):205–208.
.
OpenUrl CrossRef PubMed
↵
1. Johnson S,
2. Domínguez-García V,
3. Donetti L,
4. Muñoz MA
(2014) Trophic coherence determines food-web stability. Proc Natl Acad Sci USA 111(50):17923–17928.
.
OpenUrl Abstract/FREE Full Text
↵
1. Yodzis P,
2. Innes S
(1992) Body-size and consumer-resource dynamics. Am Nat 139(6):1151–1173.
.
OpenUrl CrossRef
↵
1. Williams RJ,
2. Martinez ND
(2004) Stabilization of chaotic and non-permanent food-web dynamics. Eur Phys J B 38(2):297–303.
.
OpenUrl CrossRef
↵
1. Belgrano A,
2. Scharler UM,
3. Dunne J,
4. Ulanowicz RE
1. Dunne JA,
2. Brose U,
3. Williams RJ,
4. Martinez ND
(2005) Modeling food–web dynamics: Complexity-stability implications. Aquatic Food Webs: An Ecosystem Approach, eds Belgrano A, Scharler UM, Dunne J, Ulanowicz RE (Oxford Univ Press, Oxford), pp 117–129.
.
↵
1. Rohr RP,
2. Saavedra S,
3. Bascompte J
(2014) Ecological networks. On the structural stability of mutualistic systems. Science 345(6195):1253497.
.
OpenUrl Abstract/FREE Full Text
↵
1. Dai L,
2. Vorselen D,
3. Korolev KS,
4. Gore J
(2012) Generic indicators for loss of resilience before a tipping point leading to population collapse. Science 336(6085):1175–1177.
.
OpenUrl Abstract/FREE Full Text
↵
1. Hastings A,
2. Powell T
(1991) Chaos in a three-species food chain. Ecology 72(3):896–903.
.
OpenUrl CrossRef
↵
1. McCann K,
2. Yodzis P
(1994) Non-linear dynamics and population disappearances. Am Nat 144(5):873–879.
.
OpenUrl CrossRef
↵
1. Galla T,
2. Farmer JD
(2013) Complex dynamics in learning complicated games. Proc Natl Acad Sci USA 110(4):1232–1236.
.
OpenUrl Abstract/FREE Full Text
↵
1. Haldane AG,
2. May RM
(2011) Systemic risk in banking ecosystems. Nature 469(7330):351–355.
.
OpenUrl CrossRef PubMed
↵
1. Farmer JD,
2. Skouras S
(2013) An ecological perspective on the future of computer trading. Quant Finance 13(3):325–346.
.
OpenUrl CrossRef
↵
1. de Jong H
(2002) Modeling and simulation of genetic regulatory systems: A literature review. J Comput Biol 9(1):67–103.
.
OpenUrl CrossRef PubMed
↵
1. Sompolinsky H,
2. Crisanti A,
3. Sommers H-J
(1988) Chaos in random neural networks. Phys Rev Lett 61(3):259–262.
.
OpenUrl CrossRef PubMed
↵
1. Wainrib G,
2. Touboul J
(2013) Topological and dynamical complexity of random neural networks. Phys Rev Lett 110(11):118101.
.
OpenUrl CrossRef PubMed
↵
1. Stadler PF,
2. Fontana W,
3. Miller JH
(1993) Random catalytic reaction networks. Physica D 63(3-4):378–392.
.
OpenUrl CrossRef
↵
1. Angelani L,
2. Parisi G,
3. Ruocco G,
4. Viliani G
(2000) Potential energy landscape and long-time dynamics in a simple model glass. Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics 61(2):1681–1691.
.
OpenUrl PubMed
↵
1. Cugliandolo LF,
2. Kurchan J,
3. Le Doussal P,
4. Peliti L
(1997) Glassy behaviour in disordered systems with nonrelaxational dynamics. Phys Rev Lett 78(2):350–353.
.
OpenUrl CrossRef
↵
1. Zhou JX,
2. Aliyu MD,
3. Aurell E,
4. Huang S
(2012) Quasi-potential landscape in complex multi-stable systems. J R Soc Interface 9(77):3539–3553.
.
OpenUrl Abstract/FREE Full Text
↵
1. Schoenberg IG
(1938) Metric spaces and completely monotone functions. Ann Math 39:811–841.
.
OpenUrl CrossRef
↵
1. Fyodorov YV
(2004) Complexity of random energy landscapes, glass transition, and absolute value of the spectral determinant of random matrices. Phys Rev Lett 92(24):240601.
.
OpenUrl CrossRef PubMed
↵
1. Fyodorov YV,
2. Nadal C
(2012) Critical behavior of the number of minima of a random landscape at the glass transition point and the Tracy-Widom distribution. Phys Rev Lett 109(16):167203.
.
OpenUrl CrossRef PubMed
↵
1. Auffinger A,
2. Ben Arous G,
3. Cerny J
(2013) Random matrices and complexity of spin glasses. Commun Pure Appl Math 66(2):165–201.
.
OpenUrl CrossRef
↵
1. Nicolaescu L
(2014) Complexity of random smooth functions on compact manifolds. Indiana Univ Math J 63(4):1037–1065.
.
OpenUrl CrossRef
↵
1. Fyodorov YV
(2015) High-dimensional random fields and random matrix theory. Markov Processes Related Fields 21(3):483–518.
.
OpenUrl
↵
1. Subag E
(2015) The complexity of spherical p-spin models—A second moment approach. arXiv:1504.0225.
.
↵
1. Azais J-M,
2. Wschebor M
(2009) Level Sets and Extrema of Random Processes and Fields (John Wiley, New York).
.
↵
1. Akemann G,
2. Baik J,
3. Di Francesco P
1. Khoruzhenko BA,
2. Sommers H-J
(2011) Non-Hermitian ensembles. The Oxford Handbook of Random Matrix Theory, eds Akemann G, Baik J, Di Francesco P (Oxford Univ Press, Oxford), pp 376–397.
.
↵
1. Akemann G,
2. Kanzieper E
(2007) Integrable structure of Ginibre’s ensemble of random real matrices. J Stat Phys 129(5):1159–1231.
.
OpenUrl CrossRef
↵
1. Sommers HJ,
2. Wieczorek W
(2008) General eigenvalue correlations for the real Ginibre ensemble. J Phys A Math Theor 41(40):405003.
.
OpenUrl CrossRef
↵
1. Borodin A,
2. Sinclair CD
(2009) The Ginibre ensemble of real random matrices and its scaling limits. Commun Math Phys 291(1):177–224.
.
OpenUrl CrossRef
↵
1. Forrester PJ,
2. Nagao T
(2008) Skew orthogonal polynomials and the partly symmetric real Ginibre ensemble. J Phys A Math Theor 41(37):375003.
.
OpenUrl CrossRef
↵
1. Akemann G,
2. Phillips MJ
(2014) The Interpolating Airy Kernels for the β = 1 and β = 4 elliptic Ginibre ensembles. J Stat Phys 155(3):421–465.
.
OpenUrl CrossRef
↵
1. Sommers H-J,
2. Crisanti A,
3. Sompolinsky H,
4. Stein Y
(1988) Spectrum of large random asymmetric matrices. Phys Rev Lett 60(19):1895–1898.
.
OpenUrl CrossRef PubMed
↵
1. Girko VL
(1995) The Elliptic Law: Ten years later. Random Oper Stochastic Equations 3(3):257–302.
.
OpenUrl
↵
1. Nguyen H,
2. O’Rourke S
(2015) The elliptic law. Int Math Res Not 2015(17):7620–7689.
.
OpenUrl Abstract/FREE Full Text
↵
1. Edelman A,
2. Kostlan E,
3. Shub M
(1994) How many eigenvalues of a random matrix are real? J Am Math Soc 7:247–267.
.
OpenUrl CrossRef
↵
1. Fyodorov YV,
2. Khoruzhenko BA
(2007) On absolute moments of characteristic polynomials of a certain class of complex random matrices. Commun Math Phys 273(3):561–599.
.
OpenUrl CrossRef
↵
1. Fyodorov YV,
2. Khoruzhenko BA,
3. Sommers H-J
(1997) Almost-Hermitian random matrices: Eigenvalue density in the complex plane. Phys Lett A 226(1-2):46–52.
.
OpenUrl CrossRef
↵
1. Efetov KB
(1997) Directed quantum chaos. Phys Rev Lett 79(3):491–494.
.
OpenUrl CrossRef
↵
1. Piterbarg VI
(1996) A͡symptotic Methods in the Theory of Gaussian Processes and Fields (Am Math Soc, Providence, RI).
.
↵
1. Nicolaescu L
(2015) A CLT concerning critical points of random functions on a Euclidean space. arXiv:1509.06200.
.
↵
1. Majumdar SN,
2. Schehr G
(2014) Top eigenvalue of a random matrix: Large deviations and third order phase transition. J Stat Mech 2014:P01012.
.
OpenUrl CrossRef

Article Alerts

Email Article

Citation Tools

Request Permissions

Proceedings of the National Academy of Sciences: 113 (25)

Table of Contents

Submit

Article Classifications

Physical Sciences
Applied Mathematics

[1] ↵
May RM
(1972) Will a large complex system be stable? Nature 238(5364):413–414.
.
OpenUrl CrossRef PubMed

[2] May RM

[3] ↵
Gardner MR,
Ashby WR
(1970) Connectance of large dynamic (cybernetic) systems: Critical values for stability. Nature 228(5273):784.
.
OpenUrl CrossRef PubMed

[4] Gardner MR,

[5] Ashby WR

[6] ↵
Porter CE
(1965) Statistical Theories of Spectra: Fluctuations. A Collection of Reprints and Original Papers (Academic, New York).
.

[7] Porter CE

[8] ↵
Allesina S,
Si T
(2015) The stability complexity relationship at age 40: A random matrix perspective. Popul Ecol 57(1):63–75.
.
OpenUrl CrossRef

[9] Allesina S,

[10] Si T

[11] ↵
James A, et al.
(2015) Constructing random matrices to represent real ecosystems. Am Nat 185(5):680–692.
.
OpenUrl CrossRef PubMed

[12] James A, et al.

[13] ↵
McCann KS
(2000) The diversity-stability debate. Nature 405(6783):228–233.
.
OpenUrl CrossRef PubMed

[14] McCann KS

[15] ↵
Boyd IL
(2012) The art of ecological modeling. Science 337(6092):306–307.
.
OpenUrl Abstract/FREE Full Text

[16] Boyd IL

[17] ↵
Allesina S,
Tang S
(2012) Stability criteria for complex ecosystems. Nature 483(7388):205–208.
.
OpenUrl CrossRef PubMed

[18] Allesina S,

[19] Tang S

[20] ↵
Johnson S,
Domínguez-García V,
Donetti L,
Muñoz MA
(2014) Trophic coherence determines food-web stability. Proc Natl Acad Sci USA 111(50):17923–17928.
.
OpenUrl Abstract/FREE Full Text

[21] Johnson S,

[22] Domínguez-García V,

[23] Donetti L,

[24] Muñoz MA

[25] ↵
Yodzis P,
Innes S
(1992) Body-size and consumer-resource dynamics. Am Nat 139(6):1151–1173.
.
OpenUrl CrossRef

[26] Yodzis P,

[27] Innes S

[28] ↵
Williams RJ,
Martinez ND
(2004) Stabilization of chaotic and non-permanent food-web dynamics. Eur Phys J B 38(2):297–303.
.
OpenUrl CrossRef

[29] Williams RJ,

[30] Martinez ND

[31] ↵
Belgrano A,
Scharler UM,
Dunne J,
Ulanowicz RE
Dunne JA,
Brose U,
Williams RJ,
Martinez ND
(2005) Modeling food–web dynamics: Complexity-stability implications. Aquatic Food Webs: An Ecosystem Approach, eds Belgrano A, Scharler UM, Dunne J, Ulanowicz RE (Oxford Univ Press, Oxford), pp 117–129.
.

[32] Belgrano A,

[33] Scharler UM,

[34] Dunne J,

[35] Ulanowicz RE

[36] Dunne JA,

[37] Brose U,

[38] Williams RJ,

[39] Martinez ND

[40] ↵
Rohr RP,
Saavedra S,
Bascompte J
(2014) Ecological networks. On the structural stability of mutualistic systems. Science 345(6195):1253497.
.
OpenUrl Abstract/FREE Full Text

[41] Rohr RP,

[42] Saavedra S,

[43] Bascompte J

[44] ↵
Dai L,
Vorselen D,
Korolev KS,
Gore J
(2012) Generic indicators for loss of resilience before a tipping point leading to population collapse. Science 336(6085):1175–1177.
.
OpenUrl Abstract/FREE Full Text

[45] Dai L,

[46] Vorselen D,

[47] Korolev KS,

[48] Gore J

[49] ↵
Hastings A,
Powell T
(1991) Chaos in a three-species food chain. Ecology 72(3):896–903.
.
OpenUrl CrossRef

[50] Hastings A,

[51] Powell T

[52] ↵
McCann K,
Yodzis P
(1994) Non-linear dynamics and population disappearances. Am Nat 144(5):873–879.
.
OpenUrl CrossRef

[53] McCann K,

[54] Yodzis P

[55] ↵
Galla T,
Farmer JD
(2013) Complex dynamics in learning complicated games. Proc Natl Acad Sci USA 110(4):1232–1236.
.
OpenUrl Abstract/FREE Full Text

[56] Galla T,

[57] Farmer JD

[58] ↵
Haldane AG,
May RM
(2011) Systemic risk in banking ecosystems. Nature 469(7330):351–355.
.
OpenUrl CrossRef PubMed

[59] Haldane AG,

[60] May RM

[61] ↵
Farmer JD,
Skouras S
(2013) An ecological perspective on the future of computer trading. Quant Finance 13(3):325–346.
.
OpenUrl CrossRef

[62] Farmer JD,

[63] Skouras S

[64] ↵
de Jong H
(2002) Modeling and simulation of genetic regulatory systems: A literature review. J Comput Biol 9(1):67–103.
.
OpenUrl CrossRef PubMed

[65] de Jong H

[66] ↵
Sompolinsky H,
Crisanti A,
Sommers H-J
(1988) Chaos in random neural networks. Phys Rev Lett 61(3):259–262.
.
OpenUrl CrossRef PubMed

[67] Sompolinsky H,

[68] Crisanti A,

[69] Sommers H-J

[70] ↵
Wainrib G,
Touboul J
(2013) Topological and dynamical complexity of random neural networks. Phys Rev Lett 110(11):118101.
.
OpenUrl CrossRef PubMed

[71] Wainrib G,

[72] Touboul J

[73] ↵
Stadler PF,
Fontana W,
Miller JH
(1993) Random catalytic reaction networks. Physica D 63(3-4):378–392.
.
OpenUrl CrossRef

[74] Stadler PF,

[75] Fontana W,

[76] Miller JH

[77] ↵
Angelani L,
Parisi G,
Ruocco G,
Viliani G
(2000) Potential energy landscape and long-time dynamics in a simple model glass. Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics 61(2):1681–1691.
.
OpenUrl PubMed

[78] Angelani L,

[79] Parisi G,

[80] Ruocco G,

[81] Viliani G

[82] ↵
Cugliandolo LF,
Kurchan J,
Le Doussal P,
Peliti L
(1997) Glassy behaviour in disordered systems with nonrelaxational dynamics. Phys Rev Lett 78(2):350–353.
.
OpenUrl CrossRef

[83] Cugliandolo LF,

[84] Kurchan J,

[85] Le Doussal P,

[86] Peliti L

[87] ↵
Zhou JX,
Aliyu MD,
Aurell E,
Huang S
(2012) Quasi-potential landscape in complex multi-stable systems. J R Soc Interface 9(77):3539–3553.
.
OpenUrl Abstract/FREE Full Text

[88] Zhou JX,

[89] Aliyu MD,

[90] Aurell E,

[91] Huang S

[92] ↵
Schoenberg IG
(1938) Metric spaces and completely monotone functions. Ann Math 39:811–841.
.
OpenUrl CrossRef

[93] Schoenberg IG

[94] ↵
Fyodorov YV
(2004) Complexity of random energy landscapes, glass transition, and absolute value of the spectral determinant of random matrices. Phys Rev Lett 92(24):240601.
.
OpenUrl CrossRef PubMed

[95] Fyodorov YV

[96] ↵
Fyodorov YV,
Nadal C
(2012) Critical behavior of the number of minima of a random landscape at the glass transition point and the Tracy-Widom distribution. Phys Rev Lett 109(16):167203.
.
OpenUrl CrossRef PubMed

[97] Fyodorov YV,

[98] Nadal C

[99] ↵
Auffinger A,
Ben Arous G,
Cerny J
(2013) Random matrices and complexity of spin glasses. Commun Pure Appl Math 66(2):165–201.
.
OpenUrl CrossRef

[100] Auffinger A,

[101] Ben Arous G,

[102] Cerny J

[103] ↵
Nicolaescu L
(2014) Complexity of random smooth functions on compact manifolds. Indiana Univ Math J 63(4):1037–1065.
.
OpenUrl CrossRef

[104] Nicolaescu L

[105] ↵
Fyodorov YV
(2015) High-dimensional random fields and random matrix theory. Markov Processes Related Fields 21(3):483–518.
.
OpenUrl

[106] Fyodorov YV

[107] ↵
Subag E
(2015) The complexity of spherical p-spin models—A second moment approach. arXiv:1504.0225.
.

[108] Subag E

[109] ↵
Azais J-M,
Wschebor M
(2009) Level Sets and Extrema of Random Processes and Fields (John Wiley, New York).
.

[110] Azais J-M,

[111] Wschebor M

[112] ↵
Akemann G,
Baik J,
Di Francesco P
Khoruzhenko BA,
Sommers H-J
(2011) Non-Hermitian ensembles. The Oxford Handbook of Random Matrix Theory, eds Akemann G, Baik J, Di Francesco P (Oxford Univ Press, Oxford), pp 376–397.
.

[113] Akemann G,

[114] Baik J,

[115] Di Francesco P

[116] Khoruzhenko BA,

[117] Sommers H-J

[118] ↵
Akemann G,
Kanzieper E
(2007) Integrable structure of Ginibre’s ensemble of random real matrices. J Stat Phys 129(5):1159–1231.
.
OpenUrl CrossRef

[119] Akemann G,

[120] Kanzieper E

[121] ↵
Sommers HJ,
Wieczorek W
(2008) General eigenvalue correlations for the real Ginibre ensemble. J Phys A Math Theor 41(40):405003.
.
OpenUrl CrossRef

[122] Sommers HJ,

[123] Wieczorek W

[124] ↵
Borodin A,
Sinclair CD
(2009) The Ginibre ensemble of real random matrices and its scaling limits. Commun Math Phys 291(1):177–224.
.
OpenUrl CrossRef

[125] Borodin A,

[126] Sinclair CD

[127] ↵
Forrester PJ,
Nagao T
(2008) Skew orthogonal polynomials and the partly symmetric real Ginibre ensemble. J Phys A Math Theor 41(37):375003.
.
OpenUrl CrossRef

[128] Forrester PJ,

[129] Nagao T

[130] ↵
Akemann G,
Phillips MJ
(2014) The Interpolating Airy Kernels for the β = 1 and β = 4 elliptic Ginibre ensembles. J Stat Phys 155(3):421–465.
.
OpenUrl CrossRef

[131] Akemann G,

[132] Phillips MJ

[133] ↵
Sommers H-J,
Crisanti A,
Sompolinsky H,
Stein Y
(1988) Spectrum of large random asymmetric matrices. Phys Rev Lett 60(19):1895–1898.
.
OpenUrl CrossRef PubMed

[134] Sommers H-J,

[135] Crisanti A,

[136] Sompolinsky H,

[137] Stein Y

[138] ↵
Girko VL
(1995) The Elliptic Law: Ten years later. Random Oper Stochastic Equations 3(3):257–302.
.
OpenUrl

[139] Girko VL

[140] ↵
Nguyen H,
O’Rourke S
(2015) The elliptic law. Int Math Res Not 2015(17):7620–7689.
.
OpenUrl Abstract/FREE Full Text

[141] Nguyen H,

[142] O’Rourke S

[143] ↵
Edelman A,
Kostlan E,
Shub M
(1994) How many eigenvalues of a random matrix are real? J Am Math Soc 7:247–267.
.
OpenUrl CrossRef

[144] Edelman A,

[145] Kostlan E,

[146] Shub M

[147] ↵
Fyodorov YV,
Khoruzhenko BA
(2007) On absolute moments of characteristic polynomials of a certain class of complex random matrices. Commun Math Phys 273(3):561–599.
.
OpenUrl CrossRef

[148] Fyodorov YV,

[149] Khoruzhenko BA

[150] ↵
Fyodorov YV,
Khoruzhenko BA,
Sommers H-J
(1997) Almost-Hermitian random matrices: Eigenvalue density in the complex plane. Phys Lett A 226(1-2):46–52.
.
OpenUrl CrossRef

[151] Fyodorov YV,

[152] Khoruzhenko BA,

[153] Sommers H-J

[154] ↵
Efetov KB
(1997) Directed quantum chaos. Phys Rev Lett 79(3):491–494.
.
OpenUrl CrossRef

[155] Efetov KB

[156] ↵
Piterbarg VI
(1996) A͡symptotic Methods in the Theory of Gaussian Processes and Fields (Am Math Soc, Providence, RI).
.

[157] Piterbarg VI

[158] ↵
Nicolaescu L
(2015) A CLT concerning critical points of random functions on a Euclidean space. arXiv:1509.06200.
.

[159] Nicolaescu L

[160] ↵
Majumdar SN,
Schehr G
(2014) Top eigenvalue of a random matrix: Large deviations and third order phase transition. J Stat Mech 2014:P01012.
.
OpenUrl CrossRef

[161] Majumdar SN,

[162] Schehr G

Main menu

User menu

Search

New Research In

Nonlinear analogue of the May−Wigner instability transition

Significance

Abstract

Model

Results

Discussion

Materials and Methods

Kac−Rice Formula.

Real Elliptic Matrices and Asymptotics of 〈N〉tot.

SI Text

Housholder Reflections and Partial Triangularization of Real Matrices

Elliptic Ensemble of Gaussian Random Matrices

Acknowledgments

Footnotes

References

Citation Manager Formats

Article Classifications

You May Also be Interested in

Similar Articles

Articles

PNAS Portals

Information

Main menu

User menu

Search

New Research In

Physical Sciences

Featured Portals

Articles by Topic

Social Sciences

Featured Portals

Articles by Topic

Biological Sciences

Featured Portals

Articles by Topic

Nonlinear analogue of the May−Wigner instability transition

Significance

Abstract

Model

Results

Discussion

Materials and Methods

Kac−Rice Formula.

Real Elliptic Matrices and Asymptotics of 〈N〉tot.

SI Text

Housholder Reflections and Partial Triangularization of Real Matrices

Elliptic Ensemble of Gaussian Random Matrices

Acknowledgments

Footnotes

References

Citation Manager Formats

Sign up for Article Alerts

Article Classifications

Jump to section

You May Also be Interested in

Similar Articles

Articles

PNAS Portals

Information