Generalized Lotka-Volterra Systems with Time Correlated Stochastic Interactions

Understanding the mechanisms driving the biodiversity patterns observed across various ecosystems has long been a central challenge in community ecology [1,2]. Traditionally, ecologists believed that more complex ecosystems, containing a greater number of species and interactions, should exhibit greater stability [3]. This notion was fundamentally challenged by May, who introduced the concept of stability bound for randomly assembled ecosystems [4], showing that the larger the number of interactions and their variability, the less stable the system is. This result is known as “diversity-stability paradox” [5]. One key area of research that has gained momentum since May’s pivotal work focuses on understanding the role of interaction networks on the stability and species coexistence within large communities [2,6–10]. Given the inherent unknowns in species interactions, several recent works have proposed modeling the dynamics of interacting species through generalized Lotka-Volterra (GLV) equations with quenched random disorder (QGLV), where the underlying interaction network is fully connected, leading to a number of interesting results [11–15]. The phase diagram of these models is essentially divided into three regions of qualitatively different behaviors: a system may converge to a fixed point, reach a multiple-attractors state, or have populations which grow indefinitely [13,14,16]. Furthermore, the addition of demographic noise to the QGLV leads to new phases such as a Gardner phase [17].

The QGLV model assumes that species interactions remain constant over time. However, empirical ecological systems are characterized by temporal fluctuations in species interactions, influenced by variations in environmental conditions, resource availability, and other factors that operate on a timescale comparable to population dynamics [18–22]. In the single fixed point phase, the stability of the QGLV model decreases as the fraction of nonzero interactions and the heterogeneity of the interaction strengths increases [13,23]. In particular, such ecological communities dynamically mitigate instability by reducing species diversity, eventually achieving a marginally stable state [13], consistent with previously obtained theoretical bounds [5]. Moreover, the distribution of the stationary populations within the ecological community (known as species abundance distribution, or SAD), as obtained in the unique equilibrium phase and in the limit of a large number of species within the dynamical mean field theory (DMFT), is a truncated Gaussian [14,16]. This distribution is very different from the heavy tail SAD observed in empirical microbial [24], plankton [25], or forest [26] communities.

In the present Letter, within the established framework of the GLV model featuring a fully connected random interaction network, we consider time-dependent species interactions and a Monod functional response, commonly used for modeling the growth of microorganisms [27].

Specifically, we adopt the hypothesis that any interaction between an arbitrary pair of species can be modeled as stochastic colored noise, which we call annealed GLV (AGLV). The introduction of temporal stochastic fluctuations in the strengths of species interactions yields results that fill some of the aforementioned gaps for the QGLV. In fact, unlike previous models where environmental noise was introduced externally with predefined statistical properties [24,28], our model incorporates an effective random environment whose characteristics are determined self-consistently. In this way, every species experiences an environment shaped by the collective presence of all other species. We find that temporal fluctuations in the species interactions promote species diversity and generate SAD that aligns with data, where the majority of species are rare, and most individuals belong to a small fraction of all species [29].

Let us consider 𝑥𝑖⁡(𝑡), the population density at time 𝑡 of the species 𝑖. Then the dynamics of the AGLV system with colored noise for 𝑆 interacting species are given by

with 𝑖 =1,…,𝑆, and where 𝛼𝑖⁢𝑗⁡(𝑡)=𝜇/𝑆+𝜎⁢𝑧𝑖⁢𝑗⁡(𝑡)/√𝑆 for 𝑖 ≠𝑗 and {𝑧𝑖⁢𝑗⁡(𝑡):  𝑡>0} are independent Gaussian random variables with ‾𝑧𝑖⁢𝑗⁡(𝑡)=0, ‾𝑧𝑖⁢𝑗⁡(𝑡)⁢𝑧𝑘⁢𝑙⁡(𝑡′)=𝛿𝑖,𝑘⁢𝛿𝑗,𝑙⁢𝑄⁡(Δ⁢𝑡|𝜏)=𝛿𝑖,𝑘⁢𝛿𝑗,𝑙⁢[(1+2⁢𝜏/𝜏0)/2⁢𝜏]⁢𝑒−Δ⁢𝑡/𝜏, where 𝛿𝑥,𝑦 is the Kronecker 𝛿, and Δ⁢𝑡 =|𝑡 −𝑡′|. 𝐽⁡(𝑥) is a generic function of 𝑥, and in particular we will consider two cases: 𝐽⁡(𝑥) =𝑥, to make direct comparison with QGLV, and 𝐽⁡(𝑥)=[𝑥/(1+𝑐⁢𝑥)], where 𝑐 is known as handling time [30]. ℎ𝑖⁡(𝑡) is a possible time-dependent external field. The amplitude of the noise 𝑄⁡(Δ⁢𝑡|𝜏) is chosen so that in the limit 𝜏 →∞ we recover the case of quenched disorder. For simplicity, we set 𝜏0,𝑟𝑖,𝐾𝑖 =1 for 𝑖 =1,…,𝑆 and work with dimensionless variables and parameters. From this general annealed formulation with colored noise, the limit 𝜏 →0 corresponds to annealed white-noise (AWN) dynamics. In Fig. 1 we show the effect of time correlated noise on the species abundance evolution for 𝐽⁡[𝑥𝑗⁡(𝑡)] =𝑥𝑗⁡(𝑡). In our model fluctuations scale like the size of the population, generating a time-dependent multiplicative noise that keeps the trajectories away from zero, thus preventing their extinction for any finite 𝜏. This behavior is substantially different from the one of models with quenched noise, where all species’ populations with a negative growth rate are doomed to extinction in the infinite time limit. Species populations undergo fluctuations between high and low abundances, with an average frequency which depends on the value of 𝜏. This fluctuating behavior is instrumental in promoting the coexistence of multiple species within the ecosystem and is present for all ranges of finite 𝜏, including the limit 𝜏 →0. As we show, the same results hold for 𝐽⁡[𝑥𝑗⁡(𝑡)]={𝑥𝑗⁡(𝑡)/[1+𝑐⁢𝑥𝑗⁡(𝑡)]}. Facilitating species coexistence through fluctuations is a mechanism that has also been observed in the chaotic phase of the QGLV [15,31,32] or similar models [29], leading to a large pool of very rare species in the systems, supporting the insurance theory of biodiversity [26,29].

FIG. 1.

The DMFT for the general AGLV Eq. (1) is given by [see Sec. 1 in Supplemental Material (SM) [33] ]

where 𝑀⁡(𝑡)=E⁡[𝑥⁡(𝑡)], E[⋅] indicates the expected value, and in the following we set ℎ⁡(𝑡) =0. The self-consistent Gaussian noise 𝜂⁡(𝑡) is such that E⁡[𝜂⁡(𝑡)] =0 and E⁡[𝜂⁡(𝑡)⁢𝜂⁡(𝑡′)]=𝑄⁡(Δ⁢𝑡|𝜏)⁢E⁡{𝐽⁡[𝑥⁡(𝑡)]⁢𝐽⁡[𝑥⁡(𝑡′)]}.

From Fig. 5 in SM [33], we can see that at stationarity, the (connected) autocorrelation function of 𝐽⁡[𝑥⁡(𝑡)] has an exponentiallike decay:

and exploiting Eq. (3) we can simplify the self-consistency for 𝜂 as E⁡[𝜂⁡(𝑡)⁢𝜂⁡(𝑡′)]=𝑄⁡(Δ⁢𝑡|‾𝜏)⁢E⁡[𝐽⁡(𝑥)2], at least in the relevant regime |Δ⁢𝑡| =|𝑡 −𝑡′| ≪𝜏𝑥 in which the connected autocorrelation function can be approximated by the nonconnected autocorrelation function, with the new effective timescale ¯𝜏=1/[1/𝜏+(1−E2⁡[𝐽⁡(𝑥)]/E⁡[𝐽⁡(𝑥)2])/𝜏𝑥]. (see Sec. 2 in SM [33] for further details).

With this simplification we can now use the unified colored noise approximation (UCNA) [34,35] on Eq. (2), which, for both cases of 𝐽⁡(𝑥) considered here, leads to the same stationary SAD:

where 𝑥 >0, 𝑍⁡(𝛿𝜏,𝐷,¯𝜏) is the normalization constant, that can be computed analytically, 𝛿𝜏=(1+𝜇⁢𝑀*)/𝐷⁡(𝜏), 𝐷⁡(𝜏)=𝜎2⁢E*⁡[𝐽⁢(𝑥)2]⁢(1+2⁢𝜏)⁢¯𝜏⁡(𝜏)/2⁢𝜏, and ¯𝑥 =1 +𝜇⁢𝑀* [ E*⁡(·) denotes the average with the distribution 𝑃*𝜏 and 𝑀*=E*⁡[𝐽⁡(𝑥)]]. Notice that 𝑃*𝜏⁡(𝑥) is essentially an interpolation between a truncated Gaussian and a Gamma distribution. The former is known to be the solution for the SAD of the DMFT in the case of random quenched interactions in the single equilibrium phase [14,16], while the latter is shown below to be the exact solution of the AWN case, corresponding to the limit ¯𝜏 ∼𝜏 →0 of Eq. (2).

UCNA is recognized for its exactness in both 𝜏 =0 and 𝜏 =+∞, proving to be a reliable interpolation method for intermediate values of 𝜏 [34,35]. This is explicitly illustrated in Fig. 2, where the analytical solution of the SADs for various 𝜏 values, obtained through both the DMFT and UCNA, is compared with numerical solutions derived from integrating Eq. (1) for a system comprising 𝑆 =100 species with 𝐽⁡(𝑥) =𝑥. In Sec. 2.2 of SM [33], we show that the same result also holds for the case with the Monod functional response. We have also performed a sensitivity analysis, showing that the analytical solutions match the numerical ones very well for a wide range of parameters.

FIG. 2.

As demonstrated in both Figs. 1 and 2, the introduction of time-dependent fluctuations in interactions promotes species coexistence. This is attributed to the induced fluctuating behavior, causing the species’ growth rate to transition from negative to positive, preventing extinctions, as also verified numerically (see Fig. 2). Only when 𝜏 =+∞ a Dirac’s 𝛿 develops at 𝑥 =0, whose mass represents the fraction of extinct species. However, in practical applications, distinguishing rare species from extinct ones in samples with a finite number of species, 𝑆, requires consideration. This can be addressed by introducing a threshold value, 𝑥𝑡⁢ℎ. As shown in Fig. 3, we find that the fraction of extinct species increases with 𝜏. This finding indicates that rapid stochastic variability in species interactions promotes coexistence, while slow or static interactions are detrimental. Each point in Fig. 3 is given by ∫𝑥𝑡⁢ℎ0𝑃*𝜏⁡(𝑥)𝑑𝑥 with 𝑃*𝜏⁡(𝑥) obtained as in Fig. 2. The horizontal dashed lines in Fig. 3 represent the number of extinct species in the quenched disorder limit [14,16]. In Sec. 2.3 of SM [33], we also show that this result holds for different values of the extinction threshold.

FIG. 3.

In the white-noise 𝜏 →0 limit the DMFT equation is the same of Eq. (2), but in this case with E⁡[𝜂⁡(𝑡)⁢𝜂⁡(𝑡′)]=Σ2⁡(𝑡)⁢𝛿⁢(𝑡−𝑡′), Σ2⁡(𝑡)=E⁡{𝐽⁢[𝑥⁡(𝑡)]2}, and the multiplicative noise term 𝑥⁡(𝑡)⁢𝜂⁡(𝑡) should be interpreted in the Stratonovich sense [36]. At stationarity, the self-consistency imposes 𝑀*=E*⁡[𝐽⁡(𝑥)] and Σ*2=E*⁡[𝐽⁢(𝑥)2]. The exact stationary distribution 𝑃*0 can be derived from the Fokker-Planck equation corresponding to Eq. (2), and it reads [for any 𝐽⁡(𝑥)]:

and it coincides with the limits of 𝑃*𝜏⁡(𝑥) when 𝜏 →0 as discussed in detail in the section “White-Noise Limit” in SM [33]. We also have lim𝜏→0⁡𝛿𝜏=2⁢(1+𝜇⁢𝑀*)/(𝜎2⁢Σ*1)≡𝛿 with 𝑀*=1/(1−𝜇) and Σ*2=𝛿⁢(𝛿+1)/𝛽2 for 𝐽⁡(𝑥) =𝑥. Similar results are obtained for the case with the functional response (see SM [33], Sec. 2.2).

For 𝐽⁡(𝑥) =𝑥 we can write explicitly the SAD’s parameters as a function of 𝜇 and 𝜎 as (see Sec. 2.4 in SM [33])

while for the case with functional response, we can simply solve numerically the equations for 𝑀* and Σ* as a function of 𝛽 and 𝛿.

The predicted SADs by Eq. (4) through the DMFT and UCNA are plotted as continuous lines in Fig. 2. The parameters are obtained by first fitting the autocorrelation, checking the agreement with the empirical parameters (error below 5%; see Sec. 2 in SM [33]). In panel (a) we also plot, as a dark blue dashed line, the AWN solution 𝑃*0⁡(𝑥) given by Eq. (5). In this case, the distribution parameters are directly calculated from Eq. (6) as a function of 𝜇 and 𝜎.

To address the empirical relevance of our theoretical findings, we have analyzed two datasets: trees from a tropical forest [37] and coral reefs [38]. Figure 4 shows the SADs using Preston (log2) scale (as usually done in this context [1]). The empirical SADs were then compared with the truncated Gaussian and Gamma distributions that can be derived from the QGLV and AGLV, respectively. We find that the empirical SADs are accurately described by the AGLV model, while the QGLV model fails to reproduce these patterns. This comparison proves the empirical relevance of our model, demonstrating its applicability to real-world ecological data.

FIG. 4.

As is well-known in the literature, the QGLV model exhibits certain pathologies that lead to the existence of a region in the phase diagram where densities diverge [39]. These divergences also persist in the AGLV model. Actually, the temporal fluctuations in interactions expand the region of the phase space where the model does not converge (see Sec. 4 in SM [33] for further details). However, when the Monod functional response is introduced through a bounded 𝐽⁡(𝑥), we prevent indefinite population growth, eliminating any divergences and dependencies on the initial conditions. Therefore, the phase diagram displays only the stable stationary state for any values of 𝜇 and 𝜎. The stationary solution, Eq. (4), does not depend on the specific form of 𝐽⁡(𝑥). The specific choice of 𝐽⁡(𝑥) enters only through self-consistencies, i.e., 𝐽⁡(𝑥) determines the random environment experienced by the species.

In this Letter, we have undertaken an investigation into the GLV equations with annealed disorder, incorporating finite correlation time and simple functional responses. We have determined the corresponding dynamical mean-field equations for a large number of species, which do not depend on the specific form of 𝐽⁡(𝑥). The inclusion of temporal stochastic fluctuations in the strengths of species interactions has resulted in a diverse range of ecologically significant outcomes. First, the introduction of annealed disorder in the GLV equations, for any finite correlation time, has exerted a substantial positive influence on the biodiversity of the system. Specifically, when the dynamics of the system converges to the stationary distribution, we observe fluctuations in the dynamics of species populations, where species abundances alternate between high and low values, favoring the coexistence of all species and supporting the insurance theory of biodiversity [26,29]. This is, in fact, a similar outcome to what QGLV models found in the chaotic phase [15,32] when introducing an immigration rate 𝜆, which explicitly prevents extinctions. Second, in the white-noise limit, the DMFT leads to the stochastic logistic model, a phenomenological model that proved to be consistent with several macroecological laws in microbial ecosystems [24,28]. In particular, the analytical species abundance distribution derived from the DMFT follows the Gamma distribution, a widely utilized probability distribution in macroecology [1,37]. Again, similar truncated fat-tailed distribution has recently been numerically found in the chaotic phase of the QGLV [32] and in the GLV with strong interactions [29,40,41] with immigration. Eventually, we have presented a refinement of the model, through the inclusion of a simple functional response. We have shown that it maintains not only the core phenomenology described above but also rectifies any nonphysical divergences observed in the classic QGLV and in the AGLV. Thus, it enhances the model’s realism and applicability without sacrificing its fundamental characteristics and predictive capabilities. This work opens various avenues of research, including the integration of quenched and annealed disorder and the correlations between pairs of interacting species [14,16] or more complex hierarchical correlation structures [42]. More generally, the methodology presented here can be exploited to study the effect of annealed disorder also in other ecological dynamics. Moreover, to further improve the AGLV model, it would be valuable to explore sparse interaction networks instead of the fully connected ones examined here and in previous works [12–17] for the quenched version. However, unlocking the dynamics of this intriguing scenario requires a more comprehensive generalization of the DMFT approach [43]. The exploration of such directions holds significant promise for advancing the modeling of large-scale ecosystem dynamics, understanding emergent macroecological patterns observed in empirical data, and investigating the influence of environmental fluctuations on species coexistence.