Generalized Lotka-Volterra Equations with Random, Nonreciprocal Interactions: The Typical Number of Equilibria

Systems of many degrees of freedom with heterogeneous and nonreciprocal (asymmetric) interactions emerge naturally when modeling neural networks , natural ecosystems , economic networks, or agents playing games . The dynamics of these systems are characterized by a large number of attractors such as equilibria, limit cycles, and chaotic attractors. Systems admitting an energy landscape, as it is the case for symmetric interactions, only display equilibria, which are the stationary points of the landscape. A rugged landscape is central in the theory of glassy systems, since local minima are associated with metastable states; as a consequence, in-depth investigations and refined tools for counting and classifying local minima of highly nonconvex landscapes have been developed extensively in the context of glassy physics . Most of these studies focused on systems admitting an energy landscape, though. Recently, the interest in nonconservative systems (devoid of an energy landscape) has grown substantially and pioneering works have shown that such systems can also display many equilibria . Developing a general theory in order to count them and to investigate their stability is a challenging goal, with potentially relevant implications for understanding the dynamics.

Here we address this problem for a prototypical nonconservative dynamical system, the random generalized Lotka-Volterra model (RGLV) that describes the dynamics of population sizes of multiple species with pairwise interactions between them. The RGLV equations are used extensively in theoretical ecology to describe well-mixed ecosystems , and they are related to models used in evolutionary game theory and in economic theory . They are known to admit a multiple equilibria phase when the variability of the random interactions is strong enough ; an interesting feature for theoretical ecology . Our main result is a full characterization of multiple equilibria in terms of average abundance, diversity, and stability as summarized in the phase portrait of Fig. . There is a general expectation that the vast majority (if not all) of the equilibria are linearly unstable when the interactions are asymmetric ; our analysis confirms this surmise, which directly implies a complex dynamical behavior, as the system can never settle in a fixed point, even at long times. In order to properly count the typical number of equilibria, we combine random matrix theory with standard tools in the theory of glasses. We thus go beyond the previous analysis performed for systems with asymmetric interactions , which focused on the average number of equilibria. The latter is in fact much larger than the former and not representative of the typical behavior of the RGLV model, as we shall show below (and as it happens in many other disordered and glassy systems).

FIG. 1.

Quenched complexity $Σ (𝜙, 𝜎)$ of uninvadable equilibria for uncorrelated interactions ( $𝛾 = 0$ ). Black lines correspond to vanishing complexity; the green dotted line to the diversity $𝜙 May (𝜎)$ above which equilibria are linearly unstable (red area); the orange dotted line to the transition between the unique ( $𝜎 < 𝜎 𝑐$ ) and the multiple ( $𝜎 > 𝜎 𝑐$ ) equilibria phases.

The RGLV equations determine the dynamics of a pool of $𝑆 ≫ 1$ species. They read

𝑑 𝑁 𝑖 (𝑡) 𝑑 𝑡 = 𝑁 𝑖 (𝑡) 𝐹 𝑖 (\to 𝑁),

(1)

where $𝑁 𝑖 (𝑡) \geq 0$ is the abundance of species $𝑖$ at time $𝑡$ . The vector $\to 𝐹$ represents the growth rates or forces:

𝐹 𝑖 (\to 𝑁) = 𝜅 𝑖 - 𝑁 𝑖 - 𝜇 𝑆 𝑆 \sum 𝑗 = 1 𝑁 𝑗 - 𝜎 \sqrt 𝑆 𝑆 \sum 𝑗 = 1 𝑎 𝑖 𝑗 𝑁 𝑗 .

(2)

Here $𝜅 𝑖$ are the carrying capacities, $𝜇$ , $𝜎$ are the average interaction strength and the variability, and $𝑎 𝑖 𝑗$ are components of a random matrix encoding the fluctuations in the interactions between the different species . To describe interactions where $𝑎 𝑖 𝑗$ and $𝑎 𝑗 𝑖$ are correlated but not exactly the same, we take them as two variables with a joint Gaussian distribution defined by covariances:

⟨ 𝑎 𝑖 𝑗 𝑎 𝑘 𝑙 ⟩ = 𝛿 𝑖 𝑘 𝛿 𝑗 𝑙 + 𝛾 𝛿 𝑖 𝑙 𝛿 𝑗 𝑘, | 𝛾 | \leq 1

(3)

corresponding to $⟨ 𝑎 2 𝑖 𝑗 ⟩ = ⟨ 𝑎 2 𝑗 𝑖 ⟩ = 1$ and $⟨ 𝑎 𝑖 𝑗 𝑎 𝑗 𝑖 ⟩ = 𝛾$ . In the extreme case $𝛾 = \pm 1$ one obtains perfect correlations $𝑎 𝑖 𝑗 = \pm 𝑎 𝑗 𝑖$ , while for $𝛾 = 0$ the interactions are uncorrelated. We focus on $𝜅 𝑖 = 𝜅$ , but the calculation can be easily generalized to heterogeneous $𝜅 𝑖$ .

Equilibria are configurations $\to 𝑁 *$ satisfying

𝑑 𝑁 * 𝑖 𝑑 𝑡 = 𝑁 * 𝑖 𝐹 𝑖 (\to 𝑁 *) = 0 \forall 𝑖, 𝑁 * 𝑖 \geq 0 .

(4)

Numerical simulations and analytical results reveal two distinct regimes for large $𝑆$ : a unique equilibrium regime in which any arbitrary initialization of the population vector converges to a fixed equilibrium $\to 𝑁 *$ which is globally stable, and a multiple equilibria regime. The transition between the two regimes takes place at $𝜎 𝑐 = \sqrt 2 (1 + 𝛾) - 1$ . Characterizing the multiple equilibria phase when $- 1 < 𝛾 < 1$ is still an open challenge as mappings to physical systems work only for $𝛾 = 1$ and $𝛾 = - 1$ . In the former case the problem is conservative and the force is obtained as the derivative of an energy, $𝐹 𝑖 (\to 𝑁) = - \partial 𝑖 𝐿 (\to 𝑁)$ with $𝐿 (\to 𝑁) = \sum 𝑆 𝑖 = 1 𝑁 𝑖 [(𝑁 𝑖 / 2) - 𝜅 𝑖 + (𝜇 / 2) \sum 𝑆 𝑗 = 1 𝑁 𝑗 + (𝜎 / 2 \sqrt 𝑆) \times \sum 𝑆 𝑗 = 1 𝑎 𝑖 𝑗 𝑁 𝑗]$ . Stable equilibria are identified with metastable states (local minima of the energy). Spin-glass techniques can be used to show that there exist exponentially many (in $𝑆$ ) metastable states, the relevant ones being marginally stable, which makes the system critical and hence very fragile to nonconservative perturbations . This formalism requires the existence of an energy landscape. When $- 1 < 𝛾 < 1$ , dynamical mean field theory has provided information on the dynamics but not directly on the equilibria. Here we tackle this challenge by the Kac-Rice formalism . To study the typical number of equilibria for $𝛾 \neq 1$ we make use of the so-called quenched Kac-Rice formalism introduced in .

There are many equilibria solving , that differ by which species are present. We classify their typical number as a function of their diversity: each equilibrium $\to 𝑁 *$ has a certain number of absent species ( $𝑁 * 𝑖 = 0$ ), and a number $𝑠 (\to 𝑁 *)$ of present species ( $𝑁 * 𝑖 > 0$ ). The diversity is defined as $𝜙 (\to 𝑁 *) = 𝑠 (\to 𝑁 *) / 𝑆 \in [0, 1]$ . This quantity is a central property in ecology, which also sets the stability of the equilibria , as we recall below. Our counting of equilibria at varying $𝜙$ is also motivated by the fact that it is not known a priori which equilibria will affect the systems dynamics (and how), at variance with equilibrium frameworks where the relevant equilibria are marginally stable minima, usually the more numerous ones (see, however, ). Therefore, determining the range of diversities where equilibria are present is crucial. We focus on uninvadable equilibria, such that $𝐹 𝑖 (\to 𝑁 *) < 0$ for any $𝑖$ such that $𝑁 * 𝑖 = 0$ (notice that similar constraints appear naturally in constraint satisfaction problems, too ). These equilibria are relevant as they are stable with respect to small positive fluctuations in the abundance of the absent species. The total number $𝒩 𝑆 (𝜙)$ of uninvadable equilibria with diversity $𝜙$ scales exponentially with $𝑆$ . As known from glassy physics, $𝒩 𝑆 (𝜙)$ is a random variable which in general does not concentrate around its average (it is not self-averaging). In this case, the typical number is obtained by focusing on the large- $𝑆$ limit of its logarithm, which does concentrate around a deterministic value $Σ (𝜙)$ :

l i m 𝑆 \to \infty l o g [𝒩 𝑆 (𝜙)] 𝑆 = l i m 𝑆 \to \infty ⟨ l o g [𝒩 𝑆 (𝜙)] ⟩ 𝑆 \equiv Σ (𝜙) .

(5)

$Σ (𝜙)$ governs the exponential scaling of the typical value of $𝒩 𝑆 (𝜙)$ : borrowing the terminology from glassy physics, we refer to it as the quenched complexity. The computation of the average of the logarithm is done via the replica trick:

⟨ l o g 𝒩 𝑆 (𝜙) ⟩ = l i m 𝑛 \to 0 l o g ⟨ 𝒩 𝑛 𝑆 (𝜙) ⟩ 𝑛 .

(6)

When evaluated at $𝑛 = 1$ the right-hand of side of Eq. gives the annealed complexity associated with the average number of equilibria : $Σ (𝐴) (𝜙) \equiv l i m 𝑆 \to \infty (1 / 𝑆) l o g ⟨ 𝒩 𝑆 (𝜙) ⟩$ . When $𝒩 𝑆 (𝜙)$ is not self averaging, $Σ (𝐴) > Σ$ : the average of $𝒩 𝑆 (𝜙)$ is dominated by exponentially rare ecosystems displaying an unusually large number of equilibria. It is therefore much larger than the typical value, which captures the properties of the ecosystems occurring with probability that is not suppressed exponentially in $𝑆$ .

The main steps of the replicated Kac-Rice computation are explained in the Supplemental Material . The value of $⟨ 𝒩 𝑛 𝑆 (𝜙) ⟩$ can be determined by introducing $𝑛$ copies of the ecosystem and by finding the probability that any $𝑛$ given vectors $\to 𝑁 𝑎$ , $𝑎 = 1, \dots, 𝑛$ satisfy Eq. simultaneously, together with the uninvadability condition. This is a function of order parameters measuring properties of the equilibria, like the amount of correlation between them. The number of equilibria is dominated (according to a large deviation principle) by specific values of these order parameters. The order parameters are the first two empirical moments of the vectors $\to 𝑁 𝑎$ and $\to 𝐹 𝑎$ , i.e., the $2 𝑛$ quantities:

𝑚 𝑎 \equiv l i m 𝑆 \to \infty \sum 𝑆 𝑖 = 1 𝑁 𝑎 𝑖 𝑆, 𝑝 𝑎 \equiv l i m 𝑆 \to \infty \sum 𝑆 𝑖 = 1 𝐹 𝑎 𝑖 𝑆,

(7)

as well as the $𝑛 (𝑛 + 1) + 𝑛 (𝑛 - 1)$ correlations (or overlaps):

𝑞 𝑎 𝑏 \equiv l i m 𝑆 \to \infty \to 𝑁 𝑎 \cdot \to 𝑁 𝑏 𝑆, 𝜉 𝑎 𝑏 \equiv l i m 𝑆 \to \infty \to 𝐹 𝑎 \cdot \to 𝐹 𝑏 𝑆, 𝑧 𝑎 𝑏 \equiv l i m 𝑆 \to \infty \to 𝑁 𝑎 \cdot \to 𝐹 𝑏 𝑆,

(8)

where $𝑧 𝑎 𝑎 = 0$ follows from . These order parameters encode the correlations in the location of the different fixed points in configuration space, which emerge because all the fixed points arise from the same interactions between the species. We consider a symmetric ansatz for the order parameters, i.e., $𝑚 𝑎 = 𝑚$ , $𝑞 𝑎 𝑏 = 𝛿 𝑎 𝑏 𝑞 1 + (1 - 𝛿 𝑎 𝑏) 𝑞 0$ , $𝑝 𝑎 = 𝑝$ , $𝜉 𝑎 𝑏 = 𝛿 𝑎 𝑏 𝜉 1 + (1 - 𝛿 𝑎 𝑏) 𝜉 0 𝑧 𝑎 𝑏 = (1 - 𝛿 𝑎 𝑏) 𝑧$ , which is the simplest approximation that takes such correlations into account. Under this assumption, the moments can be written as an integral over all possible values of the order parameters:

⟨ 𝒩 𝑛 𝑆 (𝜙) ⟩ = \int 𝑑 𝐱 𝑒 𝑆 𝑛 ¯ 𝒜 (𝐱; 𝜙) + 𝑜 (𝑛 𝑆),

(9)

with $𝐱 = (𝑚, 𝑝, 𝑞 1, 𝑞 0, 𝜉 1, 𝜉 0, 𝑧)$ , see the Supplemental Material for details of the calculation of $¯ 𝒜$ and for its explicit expression. The large deviation principle then implies that asymptotically

Σ (𝜙) = ¯ 𝒜 (𝐱 ⋆; 𝜙),

(10)

where $𝐱 ⋆$ is the solution of the saddle-point equations $[𝛿 ¯ 𝒜 (𝐱; 𝜙) / 𝛿 𝐱] | 𝐱 ⋆ = 0$ . This results in self-consistent equations for the typical properties of equilibria at fixed $𝜙$ , such as their typical average abundance $𝑚 *$ or the typical similarity between two equilibria $𝑞 * 0$ .

The Kac-Rice computation allows us to determine the linear stability of the equilibria at each given $𝜙$ with respect to perturbations $𝑁 * 𝑖 \to 𝑁 * 𝑖 + 𝛿 𝑁 * 𝑖$ of the populations of coexisting species. This depends on the spectral properties of the matrix:

𝐻 𝑖 𝑗 (\to 𝑁 *) = (\partial 𝐹 𝑖 (\to 𝑁 *) 𝑑 𝑁 𝑗) 𝑖, 𝑗 : 𝑁 * 𝑖, 𝑁 * 𝑗 > 0 .

(11)

For stable equilibria all the eigenvalues of have negative real part. The asymmetry of the matrix $𝑎 𝑖 𝑗$ implies that are themselves asymmetric random matrices . The typical eigenvalue density (neglecting possible isolated eigenvalues) of $𝐻 𝑖 𝑗$ depends on $\to 𝑁 *$ only through its diversity $𝜙$ . For

𝜙 < 𝜙 May = 1 𝜎 2 (1 + 𝛾) 2

(12)

the density has support on the negative real sector; therefore a typical equilibrium with $𝜙 < 𝜙 May$ (if it exists) is stable. At $𝜙 = 𝜙 May$ , the support of the eigenvalue density touches zero and the corresponding equilibrium is marginally stable; for larger $𝜙$ the equilibrium is unstable. The criterion for linear stability is related to that identified by May in , and we henceforth refer to it as the May stability bound. More details on the Kac-Rice computation, with a thorough discussion of the structure of the equations and their resolution, are given in .

We now present our main results, focusing on the case of uncorrelated interactions $𝛾 = 0$ and setting $𝜅 = 1$ . We find that although the saddle point values $𝐱 *$ depend explicitly on $𝜇$ , the complexity at fixed diversity does not, allowing us to discuss the behavior of $Σ (𝜙)$ as a function of $𝜎$ only. As shown in Fig. , when $𝜎 > 𝜎 𝑐$ there is a range of diversities $𝜙 \in [𝜙 𝑎 (𝜎), 𝜙 𝑏 (𝜎)]$ for which $Σ (𝜙) > 0$ [a negative annealed $Σ (𝜙)$ signifies that no equilibria exist typically ]. The RGLV equations thus admit an exponentially large number of uninvadable equilibria with a continuous distribution of diversities. All the equilibria are unstable, as their diversity exceeds the May stability bound, Eq. . In Fig. we show a cut at fixed $𝜎$ of the plot of Fig. . In addition to the quenched complexity we show the annealed one for comparison. We find that the complexity and the diversity $𝜙 m a x$ associated to the typical, i.e., most numerous equilibria at the given $𝜎$ are overestimated by the annealed calculation. Annealed and quenched complexity only coincide for small $𝜙$ . The point $𝜙 c a v$ where they begin to deviate from one another turns out to coincide with the value of diversity predicted by the cavity method discussed in Refs. . The cavity method assumes the existence of a unique stable equilibrium and allows one to characterize its abundance $𝑚$ and overlap $𝑞 1$ , by imposing consistency relations between the properties of the system with $𝑆 + 1$ and $𝑆$ species. The above result shows that despite being only approximate for $𝜎 > 𝜎 𝑐$ , this method still captures the properties of a given family of equilibria, even though they are exponentially rare with respect to the typical ones at $𝜙 m a x$ .

FIG. 2.

Complexity of equilibria as a function of their diversity, for $𝛾 = 0$ . Main panel: Complexity in the multiple equilibria phase (at $𝜎 = 4$ ). A difference between quenched (magenta) and annealed (blue) is apparent. All the equilibria are unstable ( $𝜙 > 𝜙 May$ ). Inset: Annealed complexity in the unique equilibrium phase (at $𝜎 = 1$ ), negative except at the diversity predicted by the cavity formalism consistent with the existence of a unique equilibrium.

We have studied how the properties of equilibria change as $𝜙$ is increased. Figure shows that imposing a larger diversity leads to less populated (lower average abundance $𝑚 *$ ) equilibria. Similarly, it leads to less correlated (lower overlap $𝑞 * 0$ ) equilibria. Figure shows the $𝜎$ dependence of the special values of $𝜙$ discussed above (it corresponds to Fig. seen from the top). The gray area is the support of the quenched complexity, which increases with $𝜎$ . When $𝜎 \to 𝜎 + 𝑐$ all the special values of $𝜙$ merge together and reach $𝜙 May$ . Correspondingly, the complexity vanishes.

FIG. 3.

Typical averaged population size as a function of diversity $𝜙$ for $𝜎 = 4$ and $𝜇 = 30$ , in the annealed (blue) and quenched (magenta) calculation. More diverse equilibria have a smaller averaged population size $𝑚$ , which for $𝜙 > 𝜙 c a v$ is underestimated by the annealed approximation. The inset is a enlarged plot.

FIG. 4.

Diversity vs variability diagram. The range of possible diversities is indicated by the gray region. Curves of maximal complexity are shown in magenta (quenched) and blue (annealed). The black squares give $𝜙 c a v$ . The orange dashed line corresponds to $𝜙 May$ above which all equilibria are linearly unstable. Inset: Enlargement in the vicinity of $𝜎 𝑐 = \sqrt 2$ .

Just above $𝜎 𝑐$ , where the complexity goes to zero, the quenched and annealed calculations have great discrepancies, see the inset of Fig. , probably due to the larger correlation between equilibria. In fact, the average number of equilibria (annealed calculation) is dominated by equilibria having a diversity $𝜙 a n n m a x$ for which typically there are no equilibria, i.e., the quenched complexity vanishes. This feature had already been identified in Ref. for a slightly different model arising in the context of portfolio optimization (and describing, in its ecological interpretation, species competing for a single common resource).

For larger $𝜎$ the cavity approximation underestimates more strongly the diversity (and thus the instability) with respect to that of typical equilibria at $𝜙 m a x$ . For $𝜎 < 𝜎 𝑐$ , the complexity (annealed and quenched) is non-negative only at $𝜙 = 𝜙 c a v$ , which now correctly describes the diversity of the system as there is a unique equilibrium . The analysis of the multiple equilibria also allows us to characterize thoroughly the transition to an additional phase, the unbounded phase, where some abundances diverge as a function of time; see the Supplemental Material .

Finally, let us focus on the properties of the transition to the unique equilibrium phase at $𝜎 𝑐$ . Following the terminology introduced in , this is a trivialization transition and corresponds to the point at which the total quenched complexity $Σ t o t = Σ (𝜙 M a x)$ first vanishes. The way in which $Σ t o t$ vanishes for $𝜎 \to 𝜎 𝑐$ has been the focus of several works. It has been studied in models with a quadratic single-species confinement potential within an annealed calculation . Importantly, it has also been conjectured to be connected to the emergence of chaos and of a finite Lyapunov exponent . For the RGLV model at $𝛾 = 0$ we find that the complexity grows quadratically with $𝜎$ when entering the multiple equilibria phase, $Σ t o t \sim (𝜎 - 𝜎 𝑐) 2$ as in . As found in models of recurrent neural networks, the emergence of a nonzero complexity is concomitant with the emergence of a complex dynamical behavior, including chaos and aging . We notice that the annealed approximation locates correctly the trivialization transition in this case, and also captures the quadratic increase but with a different prefactor. We do not expect this quadratic behavior to be general, unless the total complexity in the vicinity of $𝜎 𝑐$ is captured by the annealed framework. If this is not the case, our calculation suggests that one should find a different power law for $𝛾 \neq 0$ (see Ref. and the Supplemental Material for more details).

In summary, we have characterized the multiple-equilibria phase of the RGLV equations by computing explicitly the complexity of uninvadable equilibria. On a technical ground our approach, giving access to the quenched complexity, has allowed us to assess when and to what degree the annealed calculation is precise: we have found a transition at the value of diversity $𝜙 c a v$ , below which the annealed calculation is exact and above which the quenched calculation gives a quantitatively different result; the latter regime always includes the maximum of the complexity, which corresponds to the typical equilibria.

We performed the calculation assuming a symmetry of the order parameters with respect to permutations of replicas: we are thus restricting the region of parameter space where to look for solutions of the self-consistent equations obtained from the variation of . For $𝛾 = 1$ it is know that the symmetric assumption is an approximation, as is optimized by parameters that break the symmetry between the replicas. Verifying that replica symmetry breaking is not needed for generic $𝛾$ is a challenge that we leave for further studies.

Our calculations show that for nonreciprocal uncorrelated interactions all the uninvadable equilibria are linearly unstable. This marks a difference with respect to the symmetric case, where marginally stable equilibria are present and correspondingly the dynamics is glassy. With unstable equilibria, a chaotic dynamics is expected in the presence of migration and signatures of it emerge in theoretical models and even in controlled experiments . Similarly to the case of landscape studies, which were instrumental to understand glassy dynamics in terms of local minima and metastable states, it would be very interesting to connect the properties of these unstable equilibria (more generally, of heteroclinic networks formed by them ) to the dynamical behavior. We envisage that invadable equilibria also play a role in the dynamics , and the calculation of their complexity is ongoing, as well as the generalization to inhomogeneous carrying capacities $𝜅 𝑖$ .

Acknowledgements

We thank D. Fisher for discussions on this topic. V. R. also thanks P. Urbani for hints on the calculation of the quenched complexity, and B. Lacroix-A-Chez-Toine, J. Berg, J. Krug, and M. Mungan for suggestions of references. V. R. acknowledges funding by the “Investissements d’Avenir” LabEx PALM (ANR-10-LABX-0039-PALM). G. B. acknowledges support from the Simons Foundation Grant (No. 454935 Giulio Biroli). G. B. acknowledges support from the Israel Science Foundation (ISF) Grant No. 773/18. A. M. T. acknowledges support from the Israeli Science Foundation (ISF) Grant No. 1939/18.

Generalized Lotka-Volterra Equations with Random, Nonreciprocal Interactions: The Typical Number of Equilibria

Abstract

Physics Subject Headings (PhySH)

Article Text

Acknowledgements

Supplemental Material

References (76)