Numerical implementation of dynamical mean field theory for disordered systems: application to the Lotka-Volterra model of ecosystems

TL;DR

This paper develops a numerical implementation of dynamical mean field theory (DMFT) for disordered systems, specifically applying it to the Lotka-Volterra ecosystem model to analyze complex dynamics like chaos and aging.

Contribution

The paper derives DMFT via the dynamical cavity method and introduces a versatile numerical solution approach applicable to various systems.

Findings

Successfully captures chaotic and aging dynamics in the Lotka-Volterra model

Provides a general numerical framework for DMFT in disordered systems

Demonstrates the method's effectiveness through simulations

Abstract

Dynamical mean field theory (DMFT) is a tool that allows to analyze the stochastic dynamics of $N$ interacting degrees of freedom in terms of a self-consistent $1$-body problem. In this work, focusing on models of ecosystems, we present the derivation of DMFT through the dynamical cavity method, and we develop a method for solving it numerically. Our numerical procedure can be applied to a large variety of systems for which DMFT holds. We implement and test it for the generalized random Lotka-Volterra model, and show that complex dynamical regimes characterized by chaos and aging can be captured and studied by this framework.

Figures

Equations

$$\forall i=\mathrm{1..}S,\frac{d{N}_i}{dt}=\frac{r_i}{K_i}{N}_i(K_i-N_i)-N_i\sum _{j\ne i}{\alpha }_{ij}{N}_j+{\lambda }_i$$

$$\overline{{\alpha }_{ij}}=\mu /S,\overline{({\alpha }_{ij}-\overline{{\alpha }_{ij}}{)}^2}={\sigma }^2/S,\overline{({\alpha }_{ij}-\overline{{\alpha }_{ij}})({\alpha }_{ji}-\overline{{\alpha }_{ji}})}=\gamma {\sigma }^2/S$$

$$\forall i=\mathrm{1..}S,\dot{N_i}=N_i(1-N_i-\sum _{j\ne i}{\alpha }_{ij}{N}_j+h_i(t))$$

$${\tilde{N}}_i(t)=N_i(t)-\sum _{j=\mathrm{1..}S}{\int }_0^t{\frac{\delta N_i(t)}{\delta h_j(s)}|}_{h=0}{\alpha }_{j0}{N}_0(s)ds$$

$$\dot{N_0}=N_0(1-N_0-\sum _{i\ne 0}{\alpha }_{0i}{\tilde{N}}_i+h_0(t))$$

$$\sum _j{\alpha }_{0j}{\tilde{N}}_j=\frac{\mu }{S}\sum _i{N}_i(t)-\frac{\mu }{S}\sum _{ij}{\int }_0^t{\chi }_{ij}(t,s)\left(\frac{\mu }{S}+\sigma a_{j0}\right)N_0(s)ds$$

$$+\sigma \sum _i{a}_{0i}{N}_i(t)-\sigma \sum _{ij}{a}_{0i}{\int }_0^t{\chi }_{ij}(t,s)\left(\frac{\mu }{S}+\sigma a_{j0}\right)N_0(s)ds$$

$$S⟨{\chi }_{ii}{a}_{i0}{a}_{0i}⟩=S⟨{\chi }_{ii}⟩⟨a_{i0}{a}_{0i}⟩=\gamma ⟨{\chi }_{ii}⟩$$

$$S(S-1)⟨{\chi }_{ij}{⟩}_{i\ne j}⟨a_{j0}{a}_{0i}{⟩}_{i\ne j}+\sqrt{S(S-1)}\sqrt{⟨{\chi }_{ij}^2{⟩}_{i\ne j}}\sqrt{⟨a_{j0}^2{a}_{0i}^2{⟩}_{i\ne j}}Z$$

$$\sim 0+S\frac{1}{\sqrt{S}}\frac{1}{S}Z$$

$$\begin{gathered} \dot{N_0}=N_0\{1-N_0-\mu ⟨N_i(t)⟩-\sigma \eta (t) \\ +\gamma {\sigma }^2{\int }_0^t⟨{\chi }_{ii}(t,s)⟩N_0(s)ds+h_0(t)\} \end{gathered}$$

$$\dot{N}=N\{1-N-\mu m(t)-\sigma \eta (t)+\gamma {\sigma }^2{\int }_0^t\chi (t,s)N(s)ds+h(t)\}$$

$$\left\{\eqalign{m(t)&=\mathbb{E}[N(t)]\cr C(t,s)&=\mathbb{E}[N(t)N(s)]\cr\chi(t,s)&=\mathbb{E}[\left.\frac{\delta N(t)}{\delta h(s)}\right|_{h=0}]}\right.$$

$$\dot{N_i}=R_i(N_i)+I_i(N_i)\left(\sum _{1\le j_2<,..<j_p\le S}{\alpha }_{j_2,..j_p}^{i}J(N_{j_2})..J(N_{j_p})+h_i(t)\right)+f_i(N_i){\xi }_i(t)$$

$$\overline{{\alpha }_{j_1,..j_p}^i}=\mu \frac{p!}{2S^{p-1}}{\overline{({\alpha }_{j_1,..j_p}^i{)}^2}}_{con}={\sigma }^2\frac{p!}{2S^{p-1}}{\overline{{\alpha }_{j_1,..j_p}^i{\alpha }_{j_1..i..j_{k-1}{j}_{k+1}..j_p}^{j_k}}}_{con}=\gamma {\sigma }^2\frac{p!}{2S^{p-1}}$$

$$\dot{N}=R_i(N)+I_i(N)\left(\mu m+\sigma \eta +\gamma {\sigma }^2\frac{p(p-1)}{2}{\int }_0^t\chi (t,s)C(t,s{)}^{p-2}J(N(s))ds+h\right)$$

$$+f_i(N)\xi$$

$$\left\{\eqalign{m(t)&=\mathbb{E}[J(N_{t})]^{p-1}\\ C(t,s)&=\mathbb{E}[J(N_{t})J(N_{s})]\\ \chi(t,s)&=\mathbb{E}[J^{\prime}(N_{t})\left.\frac{\delta N_{t}}{\delta h_{s}}\right|_{h=0}]}\right.$$

$$\left\{\eqalign{m^{new}(t)&=\mathbb{E}_{paths}\left[J(N_{t})\right]^{p-1}\\ C^{new}(t,t^{\prime})&=\mathbb{E}_{paths}\left[J(N_{t})J(N_{t}^{\prime})\right]\\ \chi^{new}(t,t^{\prime})&=\mathbb{E}_{paths}\left[J(N_{t})\int ds\;C^{-1}(t^{\prime},s)\;\eta(s)\right]\mathbf{or}\;\mathbb{E}_{paths}\left[\chi_{i}(t,t^{\prime})\right]}\right.$$

$$\chi (t,s)={\sigma }^{-1}\mathbb{E}[J(N_t)\int dx{C}^{-1}(s,x)\eta (x)]$$

$$\left\{\eqalign{R_{i}(x)&=x(1-x)+\lambda\\ I_{i}(x)&=-x\\ J(x)&=x\\ f_{i}(x)&=0}\right.$$

$$\dot{N}=N\{1-N-\mu m(t)-\sigma \eta (t)+\gamma {\sigma }^2{\int }_0^t\chi (t,s)N(s)ds+h(t)\}+\lambda$$

$${\int }_0^t\chi (t,s)N(s)ds={\int }_0^t\chi (u)N(t-u)du{\to }_{t\to \infty }{\int }_0^{\infty }\chi (u)duN(\infty )$$

$$0=N_{\infty }\{1-N_{\infty }-\mu m_{\infty }-\sigma {\eta }_{\infty }+\gamma {\sigma }^2{\chi }_{int}{N}_{\infty }+h\}$$

$$\left\{\eqalign{m_{\infty}&=\mathbb{E}[N_{\infty}]\\ \chi_{int}&=\mathbb{E}[\frac{\delta N_{\infty}}{\delta h}]\\ \mathbb{E}[\eta_{\infty}^{2}]&=\mathbb{E}[N_{\infty}^{2}]}\right.$$

$$N_{\infty }=\max \left(0,\frac{1-\mu m_{\infty }-\sigma {\eta }_{\infty }}{1-\gamma {\sigma }^2{\chi }_{int}}\right)$$

$$\widetilde{C_c}(\omega )=\frac{{\sigma }_h^2}{\frac{(1-\gamma {\sigma }^2{\chi }_{int}{)}^2}{\varphi }-{\sigma }^2+|\omega |\frac{|{\chi }_{int}|}{{\varphi }^2}\pi p_{+}(0)}$$

$$\tilde{\chi }(\omega )=-\frac{\varphi }{1-\gamma {\sigma }^2{\chi }_{int}}+i|\omega |\log |\omega |\frac{p_{+}(0)}{(1-\gamma {\sigma }^2{\chi }_{int}{)}^2-\gamma {\sigma }^2\varphi }$$

$$\forall i=\mathrm{1..}S,\dot{N_i}=N_i(1-N_i-\sum _{j\ne i}{\alpha }_{ij}{N}_j+h_i(t))$$

$$\frac{\partial }{\partial t}{\chi }_{il}(t,t^{\prime })={\chi }_{il}(t,t^{\prime })\left[\frac{d}{dt}\log N_i(t)-N_i(t)\right]+N_i(t)\left[\delta (t-t^{\prime }){\delta }_{il}-\sum _j{\alpha }_{ij}{\chi }_{jl}(t,t^{\prime })\right]$$

$$\frac{\partial }{\partial t}{\chi }_{i\ne l}(t,t^{\prime })\sim -N_i(t){\alpha }_{il}{\chi }_{ll}(t,t^{\prime })\to {\chi }_{i\ne l}\sim {\alpha }_{il}\sim S^{-1/2}.$$

$$\dot{N_i}=f(N_i)-N_i\sum _j{\alpha }_{ij}{N}_j$$

Full text

††: \jpa

1 Introduction

A growing body of work has demonstrated that the properties of communities of interacting species can be studied using tools of statistical mechanics, with the role of the thermodynamic limit being played by the large number of species. Such high-diversity communities, with tens to thousands of species, are ubiquitous and can be found anywhere from microbes in the gut to plants in a rain forest [1]. Most of the works have focused on the properties of fixed-points of the dynamics, and much less is known about the dynamics themselves, in particular when they never reach a fixed-point.

Dynamical mean-field theory (DMFT) is a useful theoretical framework which has been often used in the past to study complex stochastic dynamics of interacting degrees of freedom (spins, agents, neurons, …) [2, 3, 4, 5, 6]. In this work we develop DMFT for models of ecosystems formed by a large number of interacting species [7, 8, 9]. In the limit of large ecosystems, interactions between different species are commonly modeled by taking random interaction strengths [10, 11]. The resulting model consists in generalized Lotka-Volterra equations with random couplings. This leads to interesting problems of statistical physics, similar to ones encountered in the theory of disordered systems. Yet, there are a number of crucial differences; in particular an ecosystem is driven by non-conservative forces, hence its dynamics cannot be mapped in general to the one of a physical system in thermal equilibrium. This leads to complex dynamical regimes which have been discussed in other fields before, mainly in neural networks [6] and game theory [4, 5] (see also [12]).

Solving numerically the equations corresponding to DMFT represents a major difficulty due to the retarded friction kernel and the non-linearity of the equations. In the past this obstacle has been solved—actually circumvented—only for simplified (spherical or truncated) spin-glass models for which DMFT equations greatly simplify and reduce to closed integro-differential equations on correlation and response functions of local degrees of freedom [13]. To the best of our knowledge, a procedure to numerically integrate DMFT is still missing (with the exception of [14] that was restricted to the case of stationary dynamics and Ising spins), especially one able to analyze the complex dynamics relevant for ecosystems. In this work, following ideas developed for DMFT of strongly correlated quantum systems [15], we develop a generic numerical scheme to solve DMFT. Our method lays foundations for the study of high-diversity ecological dynamics, but also provides general tools that can be applied to problems beyond ecology, for example in the fields mentioned above 111Our code is available in a public gitHub repository [16]..

We focus on the generalized Lotka-Volterra model of ecosystems. We first present a derivation of DMFT based on the dynamical cavity method [17], which is a more intuitive procedure compared to the usual ones based on generating functional formalism, such as Martin-Siggia-Rose-DeDominicis-Janssen [4, 9, 14]. We then detail our numerical approach for solving DMFT, and show concrete examples of its implementation. This allows us to test the method, and illustrate its ability to describe and characterize complex dynamics involving chaos and aging. We finally conclude by discussing further directions and possible future applications.

2 The random Lotka-Volterra model

In this section, we introduce the random Lotka-Volterra (rLV) model [8], which describes the dynamics of interacting species. We present its phase portrait in the limit of a large number of species.

2.1 Definition and notations

The ecosystem consists of $S$ species. Each species $i$ is characterized by its population $N_{i}(t)$ which is a positive continuous variable at all times $t$. In the absence of interactions each species may grow until saturation (e.g., due to limitations on resources). The impact of other species is modeled through bilinear interactions. A small immigration rate $\lambda_{i}$ is added so that new individuals arrive to the ecosystem from the outside. The dynamical equations read:

[TABLE]

The different parameters are the intrinsic growth rates $r_{i}$ of the species, their single-species population size (carrying capacities) $K_{i}$ in the environment and the interaction matrix $\alpha$. Within our convention, a positive coefficient $\alpha_{ij}$ indicates that the presence of species $j$ is deleterious to the species $i$, due to predation or competition over resources.

For clarity of presentation, in this article we mainly discuss the case where all $r_{i}$ and $K_{i}$ are set to unity, and the immigration rate $\lambda$ is uniform. The analytical and numerical tools described can be used more generally. Immigration will act as a regularization of the problem. All results will be derived with infinitesimal but finite immigration rate $\lambda>0$. We will separately discuss in the last section the case without immigration. The elements of the interaction matrix $\alpha_{ij}$ are i.i.d. Gaussian random variables with moments:

[TABLE]

The scaling with $S$ ensures a proper large $S$ limit. In this limit, the model becomes characterized by three parameters only: the average strength of interaction $\mu$, the variety of interactions $\sigma$, and their symmetry $\gamma$. More specifically, $\gamma$ ranges from -1 (fully antisymmetric case, where all interactions are of predation-prey type) to 1 (fully symmetric case, where an energy can be defined).

For a real ecosystem with given size $S$, the parameters $\mu$, $\sigma$ and $\gamma$ can be statistically computed from the interaction matrix $\alpha$. Our result then stands for this ecosystem with the relevant parameters value. From numerical simulations, we find that ecosystems with $S>200$ are well described by results obtained in the ”thermodynamic” limit $S\rightarrow\infty$.

2.2 Phase diagram

In the large-$S$ limit, three different dynamical phases are found [8]; see figure 1.

• •

Phase I: Unique Equilibrium. In this regime, corresponding to small $\sigma$, the ecosystem displays only one stable equilibrium. Whatever the initial conditions, each species asymptotically ends up with a given number of individuals which is always the same (it can be zero as some species go extinct). This equilibrium state is stable to local and global perturbations. On figure 2 we display the dynamics of an ecosystem in this phase: each line represents the time evolution of the population of one species.

• •

Phase II: Multiple Attractors. When the variability in the interactions $\sigma$ is increased, the single stable fixed point loses its stability, and the system is left with a huge number of (possibly unstable) equilibria. This phase exhibits a complex dynamics with chaos (or aging dynamics for $\gamma=1$). An example of such dynamics can be seen on figure 3.

• •

Phase III: Unbounded Growth. When the average interaction is negative enough ($\mu<-1$), the interactions are cooperative enough to have a beneficial effect on any given species that overrides the single-species saturation. If we fix a higher $\mu$ and increase the standard deviation $\sigma$, at some point a small community of species will have cooperative interactions stronger than their own saturation and this subgroup of species will thus grow without bound, even though all the other species will die out. This explains the existence of phase III also for $\mu>-1$ for a large-enough $\sigma$. An example of such dynamics is displayed in figure 3. It should be noted that the divergence occurs as a finite time explosion of the ecosystem. The unbounded growth is a pathology of the model that could be cured by a saturation stronger than quadratic.

The borders between phases can be computed analytically: I/II and I/III are exact, but II/III is only approximate [8]. They are shown on figure 1.

The symmetric case $\gamma=1$ is special in the sense that the Multiple Attractor phase is not a chaotic one, but rather a spin glass one [7]: the dynamics gets slower and slower as the system approaches marginally stable states. In this case, a kind of physical energy can be defined and serves as a Lyapunov function. The dynamics corresponds to a gradient descent in a rough energy landscape.

3 The Dynamical Mean Field Theory

In this section, we derive the Dynamical Mean Field Theory (DMFT) using the dynamical cavity method [17]. For simplicity, we first present the derivation in the simplest case of random Lotka-Volterra model, then we extend the result to more general models, and finally we explain a numerical method to solve the DMFT equation. We also checked the relevance of the description by comparing the DMFT results with direct simulations, increasing the size $S$ of the ecosystem.

3.1 Derivation via the cavity method

For simplicity, DMFT is first derived with the simplest random Lotka-Volterra model, presented in equation (1). Our approach holds in more general cases, we will present its generalization in section 3.2. We start from the Lotka-Volterra equations:

[TABLE]

where we have added an external field $h_{i}$ that will be necessary to define the response of the system to a perturbation. The initial conditions are sampled from a product measure: $P\{N_{i}(t=0)\}=\prod_{i=1}^{S}P(N_{i}(t=0))$. For instance, we generally use a uniform distribution in $[0,1]$ for simulation purposes.

The main steps of the derivation are the following:

For given parameters $\mu$, $\sigma$, $\gamma$ and system size $S$, consider a system whose interactions and initial conditions are drawn for the $N_{i=1..S}$ species; 2. 2.

Following the dynamics according to equations (1) defines the trajectories $N_{i=1..S}(t)$; 3. 3.

We add a new species $N_{0}$, and therefore draw its initial condition $N_{0}(0)$ and the interactions $\alpha_{i0}$ and $\alpha_{0i}$ for $i=1..S$; 4. 4.

If $S$ is large enough, the impact of this new species on the previous trajectories is a small perturbation and therefore we only consider linear response for the trajectories $\tilde{N}_{i=1..S}(t)$ in the presence of species ‘0’:

[TABLE]

The partial derivative are to be understood in a functional sense. We introduce the notation $\chi_{ij}(t,s)=\left.\frac{\delta N_{i}(t)}{\delta h_{j}(s)}\right|_{h=0}$ 5. 5.

We plug these new trajectories in the equation for $N_{0}$:

[TABLE]

We introduce the matrix $a_{ij}$: $\alpha_{ij}=\mu/S+\sigma a_{ij}$, so that $a_{ij}$ is a Gaussian with zero mean and $1/S$ variance, verifying in addition $\overline{a_{ij}a_{ji}}=\gamma/S$. All the sums $\sum_{i}$ stand for $\sum_{i=1..S}$, so the interaction term reads:

[TABLE] 6. 6.

We take the large S limit and analyze the statistical properties of all terms. The main idea is that by construction $N_{i=1..S}$ are independent from $\alpha_{i0}$ and $\alpha_{0i}$, therefore one can use central-limit-like arguments. Henceforth, the notation $\langle.\rangle$ refers to the average over the couplings $a_{ij}$ and initial conditions $N_{i}(0)$. We will detail the procedure for the response function term as an example. We start from $\sum_{ij}a_{0i}\chi_{ij}(t,s)a_{j0}$. We consider that the different $\chi_{ij}(t,s)$ are random functions that will depend on the initial conditions $N_{i>0}(0)$ and the interaction matrix $a_{ij>0}$, but are otherwise independent from $a_{j0}$ and $a_{0i}$. We first treat the diagonal part. According to the central limit theorem and up to second order contribution, the term $\sum_{i}a_{0i}\chi_{ii}a_{i0}$ will converge towards its average:

[TABLE]

We now focus on the non-diagonal part. Its average is zero because $\langle a_{0i}a_{j0}\rangle_{i\neq j}=0$. To determine the scaling of its fluctuations we evaluate the variance of its single components obtaining $\langle\chi_{ij}^{2}\rangle_{i\neq j}\langle a_{j0}^{2}a_{0i}^{2}\rangle_{i\neq j}$. It can be shown by perturbation theory in the strength of interactions that $\chi_{ij}$ is of order $S^{-1/2}$ for $i\neq j$ [17] (see A). Regrouping the scalings, we obtain that $\sum_{i\neq j}a_{0i}\chi_{ij}a_{j0}$ behaves as:

[TABLE]

where $Z$ is a centered standard Gaussian. This shows that the non-diagonal term induces corrections of order $S^{-1/2}$ and can therefore be neglected in the large-S limit. After careful evaluation of all terms in equation (2) according to the same procedure, we get:

[TABLE]

where $\eta(t)$ is a Gaussian noise with zero mean and covariance $\langle\eta(t)\eta(s)\rangle_{\eta}=\langle N_{i}(t)N_{i}(s)\rangle$. 7. 7.

Since nothing differentiates $N_{0}$ from any other species, we obtain the self-consistent equation that leads to dynamical mean field theory:

[TABLE]

where $\eta$ is a Gaussian noise with zero mean and correlator $C(t,s)$, and $m(t)$, $C(t,s)$ and $\chi(t,s)$ are given functions. They are self-consistently determined with the relations:

[TABLE]

In these definitions, the averages $\mathbb{E}[.]$ are now taken with respects to the noise trajectories $\eta$ and the initial condition $N(0)$. Therefore, the equation is self-consistent in the three following functions: the average population $m(t)$, the correlator of the noise $C(t,s)$ and the averaged response function $\chi(t,s)$.

To sum up, we started from an $S$-body deterministic system of differential equations, and ended up with a one-body stochastic self-consistent differential equation 222The derivation is similar to the one of the Langevin equation from Newtonian dynamics [18], with the extra-ingredient that the bath corresponds to the rest of the system whose behavior can be self-consistently obtained from the one of $N_{0}$.. It has been mathematically proven [19] for spin glasses that when $S\to\infty$, there is a convergence in law between the statistics of the two descriptions. We expect that this holds true for our class of models as well, due to the similarity of both the equations and the method.

An important additional remark is that the DMFT is valid as long as we consider times that do not diverge with system size $S$. Otherwise, one cannot neglect terms vanishing with $S$ as we did.

3.2 DMFT equation for a general class of models

The derivation above can be performed almost identically in more general cases. The only additional subtlety is that we use the fact that the correlation $\langle N_{i}(t)N_{j}(t)\rangle$ scales as $S^{-1/2}$ for $i\neq j$, as can be shown by perturbation theory in the strength of interactions [17] (see A). Below, we just present the result for a general class of dynamics with a generic and species-dependent response function $R_{i}(N_{i})$, non-linear $p$-body interactions due to $I_{i}(N_{i})$, $J(N_{j})$ and a species scaled thermal noise $f_{i}(N_{i})\xi_{i}(t)$.

[TABLE]

where $\xi_{i}$ is a Gaussian white noise, with variance $2\omega^{2}$. The $i$-dependence of the functions denotes the possible presence of random parameters for each species. For instance, in the general Lotka Volterra case, $R_{i}(N_{i})=r_{i}/K_{i}(K_{i}-N_{i})$ where the $r_{i}$ and $K_{i}$ respectively correspond to species-dependent growth rates and carrying capacities, that we will treat as random variables sampled from given distributions. The coupling tensor satisfies $\alpha^{i}_{j_{1},..j_{p}}=0$ if there exists $k$ such that $i=j_{k}$, so as not to interfere with the self-interaction $R_{i}(N_{i})$. Otherwise, its cumulants are taken as:

[TABLE]

where the notation $\overline{X}_{con}$ denotes the connected average of $X$, i.e. when subtracting their average to the elements. Because of the constraint $1\leq j_{1}<,..<j_{p}\leq S$, when considering the cross correlation, there is only one place for the upper index $i$ to go down.

Within this set-up, the DMFT equation for a given species reads:

[TABLE]

where $\eta$ is a Gaussian noise with zero mean and covariance $\frac{p}{2}C(t,s)^{p-1}$, and $\xi$ is a Gaussian white noise, with variance $2\omega^{2}$. Using subscripts for the different times, we obtain the self-consistent closure:

[TABLE]

where the average $\mathbb{E}[.]$ is now taken with respects to the initial condition distribution, the distribution of species-dependent parameters in the functions $R_{i}$, $I_{i}$ and $f_{i}$, the noise trajectory $\eta$ and the thermal noise $\xi$.

It should be stated that the DMFT we derived with the dynamical cavity technique can also be obtained using generating functional technique of Martin-Siggia-Rose-DeDominicis-Janssen [4, 9, 14].

3.3 Solving numerically the DMFT equation

It is difficult to solve numerically a self-consistent equation where the self-consistency applies to functions. We implemented a strategy which works as pictured in figure 4. In this section, we write down in details the methodology of the algorithm. The different steps of the program are the following:

We start from initial guesses for the correlator matrix $C(t,s)$, the average vector $m(t)$, and the response matrix $\chi(t,s)$. For instance we tried a random average vector for $m$, diagonal or random positive symmetric matrices for $C$, and lower triangular random matrices for $\chi$ (since $\chi$ is a causal function); 2. 2.

Using the correlator, we can sample a Gaussian path as a simple multivariate Gaussian random variable with covariance matrix $p/2C(t,s)^{p-1}$. We draw many ($\#_{traj}$) such Gaussian paths. 3. 3.

For each path, we use our guesses $m(t)$ and $\chi(t,s)$ to numerically integrate the DMFT equation where the initial condition is sampled according to the wanted distribution. We used the uniform measure on $[0,1]$ for example. For each Gaussian path, we get a different population trajectory. 4. 4.

From these trajectories, we compute the updated values of the average population vector, the correlator matrix and the response matrix (see below), using the self-consistent closure:

[TABLE] 5. 5.

We update softly the set of functions: $X^{updated}=(1-a)X+aX^{new}$ with $X$ being respectively $m$, $C$ and $\chi$. The soft reinjection is necessary for the algorithm to converge, and not jump erratically from functions to functions. A reinjection parameter $a=0.3$ seems to be a good choice. 6. 6.

We start a new iteration of the loop, with the updated set of functions.

The convergence of the algorithm is of exponential form in the number of iterations, and is independent of the initial set of functions. On figure 5, we show an example of such a convergence.

Now, let us explain point (iv) in more detail. Obtaining the average population vector, and the correlator matrix from the trajectories is a trivial procedure: one just needs to average. Evaluating the response function $\chi$ is instead more tricky. We studied two different complementary, or alternative, procedures. The first one consists in using Novikov’s theorem [20] (or Stein’s lemma) in the statistical field formulation in order to obtain:

[TABLE]

In this formulation, $C^{-1}$ denotes the matrix inverse of $C$. The detailed derivation is presented in B. This expression is easy to implement, however it is sometimes too greedy for numerics. For instance, in our problem with multiplicative noise, the number of DMFT trajectories to average upon in order to obtain a satisfactory estimate for the response function is too high. We thus derive another relation, by directly applying $\frac{\delta}{\delta h(t^{\prime})}$ to (6). In this way, for each trajectory $i$, we can compute the response function $\chi_{i}(t,t^{\prime})$ via temporal integration, and eventually average over trajectories to obtain $\chi(t,t^{\prime})$. This procedure is less greedy in terms of needed trajectories, however it is of higher numerical complexity in the number of time steps. The details are given in C. In D, we sum up and compare the adequacy of the two methods.

The algorithm we presented here can still be improved in several ways. More specifically, when the response function is needed (when $\gamma\neq 0$), the above algorithm is quite expensive numerically. The complexity of the algorithm might be reduced by proceeding in time slices. Indeed, instead of trying to make the observables converge for the whole time interval $[t_{0},t_{f}]$, it should be numerically more efficient to make them converge on $[t_{0},t_{1}]$, then on $[t_{1},t_{2}]$ using the already converged result of $[t_{0},t_{1}]$, and so on.

The details of the numerical implementation are in G, and a public gitHub repository with the corresponding Python programs is available [16].

3.4 Numerical check of the results

We checked that the numerical solution of DMFT is consistent with the one from direct simulations. More specifically, we sample $\#_{instances}=200$ interaction matrices and initial conditions for an ecosystem of size $S$, run the deterministic dynamics, and aggregate the observables by averaging over the $S$ species and $\#_{instances}$ realizations. This is what we call direct simulations. In the Unique Equilibrium phase, the agreement is excellent. In figure 6, we show the comparison between DMFT and direct numerical simulations in the Multiple Attractors phase. As $S$ increases the direct simulations observables converge at all times to the one from DMFT. It is surprising however that the direct simulations are so different from DMFT for $S=200$. We reckon it is related to the fact that at finite $S$, when sampling the interaction matrix with parameters in the Multiple Attractor phase, there is a non-zero probability to get an interaction matrix that describes an Unbounded Growth ecosystem. This problem makes it difficult to have clean data using direct simulations, and underlines the relevance of DMFT analysis.

4 Application to the random Lotka-Volterra model

From the general case of equation (6), we recover the random Lotka-Volterra DMFT by taking:

[TABLE]

From this we obtain:

[TABLE]

The self-consistent closure is that of Eq. (4). In this section, we show how to get back the stationary results [8] from DMFT, and we study the stability of such stationary solution. All the following results are valid for $\lambda>0$ where the limit $\lambda\to 0$ is subsequently taken.

4.1 How to recover the stationary results

If the ecosystem parameters belong to the Unique Equilibrium phase, each species will eventually reach a final population value and stops changing. We describe this final state using DMFT. The one-species stochastic process becomes time-independent, so the derivative is zero, the average $m(t)$ converges to a number $m(\infty)$, the population $N(t)$ and the Gaussian noise $\eta(t)$ converge to random variables $N(\infty)$ and $\eta(\infty)$. As the process is stationary, we treat the memory kernel as time-translational invariant $\chi(t,s)=\chi(t-s)$ and therefore:

[TABLE]

Introducing the integrated memory kernel $\chi_{int}=\int_{0}^{\infty}du\;\chi(u)$, the DMFT equation (9) finally converges to:

[TABLE]

[TABLE]

From equation (10), $N_{\infty}$ can either be 0 or $(1-\gamma\sigma^{2}\chi_{int})^{-1}(1-\mu m_{\infty}-\sigma\eta_{\infty})$. By a simple linear stability analysis performed on the real system of $S$ species (see E.1), it can be shown that the 0 solution is linearly unstable when the other solution is positive. Therefore, we obtain:

[TABLE]

so the random variable $N_{\infty}$ follows a Gaussian distribution, truncated for negative abundances. We write the closed system of equations in E.2, and show that we end up with the same system as the one from [8]. From it, all observables can be computed numerically as a function of the parameters $(\mu,\sigma,\gamma)$: the fraction of species coexisting at the fixed-point, the mean abundance $N$ of species that survive and the mean response function. Some further analytical results can be derived as well, such as identifying in parameter space the boundary between the Unique Equilibrium phase and the Unbounded Growth phase. However, this analysis is only exact when we are in the Unique Equilibrium phase. It becomes approximate in the Multiple Attractors phase.

4.2 Dynamical stability and the transition line to Multiple Attractors

We now describe the loss of stability of the Unique Equilibrium solution when increasing the variability $\sigma$ of interactions: a dynamical phase transition takes place. The setup follows the one of [5]: starting in the Unique Equilibrium phase, we let the system reach an equilibrium point, then add some small field $h(t)$ which we will take as a Gaussian white noise with covariance $\overline{h(t)h(s)}=\sigma_{h}^{2}\delta(t-s)$, and see how the system responds in perturbation theory. In order to do so, we consider the DMFT equation (9), linearize it around a stationary solution, and perform a Fourier analysis [5]. The detailed calculations are presented in F. Introducing $\tilde{X}$ the Fourier transform of $X$, we obtain the small frequency expansions for both the connected correlator $C_{c}(t,t^{\prime})=\mathbb{E}[\eta(t)\eta(t^{\prime})]-\mathbb{E}[\eta_{\infty}^{2}]$ and the response function:

[TABLE]

where $\phi$, $\chi_{int}$ and $p_{+}(0)$ are properties of the Unique Equilibrium we started from. They correspond respectively to the fraction of surviving species, the integrated response to perturbations, and the value in $0^{+}$ of the surviving species’ distribution. They can be computed using the stationary cavity equations in E.2.

Different things should be noted about these expansions. First, it can be checked that we obtain the same zeroth order condition for the response function as in the stationary cavity study: $\tilde{\chi}(\omega=0)=\chi_{int}$. Secondly, the correlator initially behaves as $(a+|\omega|)^{-1}$ which corresponds to a temporal decay as $1/t^{2}$. But a change of behaviour is displayed when zeroth order term $a$ goes to zero: we observe a $1/f(\sim 1/\omega)$ correlation spectrum, which is an indicator of the chaotic transition [5]. Indeed, with this criterion we find the same transition in parameter space as the one from random matrix theory (the line $\sigma_{c}=\frac{\sqrt{2}}{1+\gamma}$ in the phase portrait in figure 1). Surprisingly, the response function instead does not exhibit a transition at $\sigma_{c}$, except for $\gamma=1$ where the fluctuation-dissipation theorem establishes a direct link between the correlation function and the response function. More complex response functions might be needed to locate the transition in the general case.

5 Numerical solution for the random Lotka-Volterra model

In this section, we present some numerical results for the random Lotka-Volterra DMFT, and show the consistency of both analytics and numerics. The aim is to illustrate the quality of the DMFT results, and present a first description of the dynamical phases (a more complete one will be presented elsewhere).

5.1 Result in the Unique Equilibrium phase

We focus on the correlator $C(t,t^{\prime})=\mathbb{E}[N(t)N(t^{\prime})]$. In the Unique Equilibrium phase, it reaches a plateau as each trajectory converges to a random constant. Moreover, the value of the plateau coincides with the stationary cavity observable $q$. This is indeed the case, as pictured on figure 7. The convergence to the stationary solution is a good check of the validity of our numerical strategy. It is shown more precisely on figure 5.

5.2 Result in the Multiple Attractors phase

In the Multiple Attractor phase we expect a different behaviour. The system does reach a time-translational invariant (TTI) chaotic state. This means that the one-time observables (the mean population $m(t)$, the proportion of alive species $\phi(t)$, or the equal-time correlation $C_{\sigma}(t,t)$) converge to a constant, and the two-time observables become functions of the time difference: $C_{\sigma}(t,t^{\prime})=C_{\sigma}(t-t^{\prime})$. If we focus on large times, we expect a relaxing behaviour for the correlator, as the trajectory decorrelates from itself when it explores the phase space along the chaotic attractor. We observe this phenomenon in the numerical solutions. Moreover, the TTI state depends on how deep in the Multiple Attractors phase the system is. On figure 8, we show the dependence on $\sigma$ of the TTI correlation $C_{\sigma}(t-t^{\prime})$, rescaled as follows. These functions $C_{\sigma}(\tau)$ starts at a TTI value for the equal-time correlation $C_{\sigma}(0)$, then as the trajectories decorrelate from themselves $C_{\sigma}(\tau)$ relaxes towards a TTI final value $C_{\sigma}(\infty)$ over a timescale $a_{\sigma}$. We therefore plot $\frac{C_{\sigma}(t-t’)-C_{\sigma}(\infty)}{C_{\sigma}(0)-C_{\sigma}(\infty)}$. We also denote $Q_{\sigma}=C_{\sigma}(0)-C_{\sigma}(\infty)$ the amplitude of the decorrelation. It is representative of the chaos strength, and this is an order parameter for the chaotic transition. On figure 9 we show the dependence of both the chaos strength $Q_{\sigma}$ and the time scale $a_{\sigma}$ as a function of the chaotic depth $\sigma-\sigma_{c}$. As expected, the chaos strength $Q_{\sigma}$ increases and the chaos time scale $a_{\sigma}$ decreases with the chaotic depth. Our results show that chaos emerges through a second-order out of equilibrium dynamical phase transition. A first attempt to obtain critical exponents is shown in figure 9. A more thorough study will be presented elsewhere.

We recall that we have considered small but finite immigration. Dynamics without immigration is different, as we discuss below.

5.3 Aging dynamics without immigration

We now consider the effect of the absence of immigration on the chaotic dynamics. The main issue is that chaos induces fluctuations that can drive species to extinction in absence of immigration and, hence, potentially kill chaos itself. The sustainability of chaotic dynamics without immigration is therefore far from being granted, actually a very different dynamical behavior can be present when $\lambda=0$. Here we show that this is indeed the case for $\gamma=0$. In figure 10 we compare the correlation functions, normalized by its equal time value, obtained by DMFT for $\gamma=0$ with and without immigration. In the former case (left panel), it is clear that a stationary chaotic state establishes as $C(t,t^{\prime})$ becomes a function of $(t-t^{\prime})$ at large times. On the contrary, without immigration (right panel), $C(t,t^{\prime})$ shows the aging behavior characteristic of glassy system: the correlation function is not a function of $t-t^{\prime}$ and displays a relaxation that is slower the older is the system. This is a nice illustration of how our numerical implementation of DMFT allows to unveil the existence of different and complex dynamical behaviors.

A detailed understanding of the aging chaotic behavior shown in figure 10, its dependence on the degree of asymmetry $\gamma$, and a thorough analysis of how and when chaos fades away is left for a future work.

6 Conclusion

In summary, we have presented a general derivation of DMFT for models of ecosystems based on the dynamical cavity method. We have implemented and tested our numerical method for generalized Lotka-Volterra models of ecosystems and showed that it can capture complex dynamics such as chaos and aging. Future works will be devoted to a thorough analysis of these complex dynamical regimes, and also to improvements of our algorithm along the ways discussed in this paper.

The main contribution of our work is the development of a numerical method to solve DMFT that can be used for many different systems characterized by stochastic dynamics and by a large number of degrees of freedom. One important potential application is to the dynamics of interacting particle glassy systems in the limit of infinite dimensions for which mean-field dynamical equations were derived recently [21, 22].

This work was partially funded by Capital Fund Management - Fondation pour la Recherche, and Simons Foundation collaboration Cracking the Glass Problem (No. 454935 to G. Biroli). G. Bunin acknowledges support by the Israel Science Foundation (ISF) Grant no. 773/18. This research was supported in part by the National Science Foundation under Grant No. PHY-1748958. We thank A. Altieri, J.-P. Bouchaud, D.S. Fisher, T. Galla, A. Manacorda, V. Ros and F. Zamponi for useful interactions on this topic.

Appendix A Scaling of the cross response function and cross correlation

A.1 Cross reponse function

We start from the simplified rLV case:

[TABLE]

We derive this equation with respect to $h_{l}(t^{\prime})$, in the functional sense. We’ll denote $\chi_{il}(t,t^{\prime})=\frac{\delta N_{i}(t)}{\delta N_{l}(t^{\prime})}$. We obtain:

[TABLE]

We know the diagonal response to be of order one at short time, and decaying: $\chi_{ii}\sim 1$. We expect that the response conveyed through correlation loop will be subleading compared to the diagonal response: $\chi_{ii}\gg\chi_{i\neq l}$. Therefore, its scaling can be inferred by considering only the contribution from the diagonal in its time evolution:

[TABLE]

A.2 Cross correlation

If we had the Fluctuation-Dissipation theorem, we would directly get the scaling of the cross correlation $C_{i\neq l}\sim\chi_{i\neq l}\sim S^{-1/2}$. However, this is not the case, but we show that the scaling relation still holds here. We start again from the simplified rLV case:

[TABLE]

where $f(N_{i})=N_{i}(1-N_{i})$, and remember that $\alpha_{ii}=0$. We use perturbation theory in the interaction matrix $\alpha=||\alpha||$. Denoting the solution of the equation $N_{i}^{0}$ when $\alpha=0$ and introducing $N_{i}=N_{i}^{0}+\alpha\delta N_{i}+\alpha^{2}...$, we obtain the first order correction for two different species $i\neq l$ through linear response:

[TABLE]

[TABLE]

From this relation, we compute the connected averages:

[TABLE]

We remark that the last term $\langle N_{j}^{0}(t_{1})N_{j^{\prime}}^{0}(t_{2})\rangle_{con}$ corresponds to the connected correlation $C_{jj^{\prime},con}(t_{1},t_{2})$. Then, from the same argument as the cross response function, we expect the cross correlation to be subleading compared to the diagonal one: $C_{ii,con}\gg C_{i\neq l,con}$. Therefore we only consider the diagonal contributions:

[TABLE]

The last term is a random variable with average $\mu^{2}/S$, and variance $\sigma^{4}/S$, therefore we obtain $C_{i\neq l}\sim S^{-1/2}$.

Appendix B Novikov’s theorem and generating functional formalism

We want to evaluate: $\chi(t,t^{\prime})=\mathbb{E}[J^{\prime}(N_{t})\left.\frac{\delta N_{t}}{\delta h_{t^{\prime}}}\right|_{h=0}]$. For simplicity, we forget about averaging over initial conditions and thermal noises, it does not change the proof. We introduce the distribution of the population trajectories:

[TABLE]

where the brackets denotes functional distributions. $\eta$ is a Gaussian noise, so its probability measure is given up to a normalization factor by:

[TABLE]

Furthermore, the $\mathbb{P}\{N|\eta\}$ distribution is deterministic and follows the DMFT dynamics. It is therefore a Dirac-distribution:

[TABLE]

We now have everything to write down the average:

[TABLE]

When we write the average as an integration over the different paths, they become non-correlated variables. The correlation aspect is taken care of in the distributions. In addition, the distribution $\mathbb{P}\{N|\eta,h\}$ is symmetric in $h(t^{\prime})$ and $\sigma\eta(t^{\prime})$. We also perform an integration by part and find:

[TABLE]

Appendix C Temporal integration of the response function

We remind the DMFT equation (6) for a general class of models here, and we will consider each trajectory (denoted with $i$) simulated through this equation:

[TABLE]

We now apply $\frac{\delta}{\delta h_{i}(t^{\prime})}$. In this way, for each trajectory $i$, we can compute the response function $\chi_{i}(t,t^{\prime})=\frac{\delta N_{i}(t)}{\delta h_{i}(t^{\prime})}$ via temporal integration:

[TABLE]

We thus construct $\chi_{i}$ by temporal integration in $t$ at fixed $t^{\prime}$, using the initial conditions $\chi_{i}(t,t^{\prime})=0$ for $t<t^{\prime}$ from causality. Eventually, we get:

[TABLE]

Appendix D Comparison of the different methods for the response function

[TABLE]

Appendix E Closure in the Unique Equilibrium phase

E.1 Linear stability of dead species

We consider the rLV system of equations (1). For each species, there are two possible equilibria [math] or $N_{i}^{*}=1-\sum_{j\neq i}\alpha_{ij}N_{j}$. In total, assuming the reduced matrix is almost always invertible (which is reasonable), this gives $2^{S}$ possible equilibria for the ecosystem, from which we would have to substract the unreachable ones with negative populations. We can linearize the equation around both possible choices for one species:

[TABLE]

Eventually, we see that the dead species solution is linearly unstable if the alive solution $N_{i}^{*}$ is positive.

E.2 Unique Equilibrium system of equations

We now use the stationary cavity solution from equation (12). The species population distribution is a truncated Gaussian: $p(n)=\phi p_{+}(n)+(1-\phi)\delta(n)$ where $\phi$ is the fraction of surviving species, and $p_{+}$ is a Gaussian distribution whose parameters need to be determined. We inject this form in the closure system (11). We introduce the parameters $q=\mathbb{E}[N_{\infty}^{2}]$, $\Delta=(1-\mu m_{\infty})\sigma^{-1}\sqrt{q}^{-1}$ and the functions $w_{k}(\Delta)=\int_{-\infty}^{\Delta}(\Delta-s)^{k}Ds$ where $Ds$ is the standard Gaussian measure, and eventually obtain:

[TABLE]

It is worth noting that, with the implicit dependence $w_{n}=w_{n}(\Delta)$, the system can be rewritten as:

[TABLE]

Under this form, the first line of system (17) gives $\Delta(\sigma,\gamma)$. Afterwards, we directly have $\chi_{int}(\sigma,\gamma,\Delta)$ and $q(\mu,\sigma,\gamma,\Delta)$.

This system can be numerically solved in the variables $(\chi_{int},\Delta,q)$ as functions of the parameters $(\mu,\sigma,\gamma)$. All observables can then be computed from the solution. For example, the proportion of alive species $\phi=w_{0}(\Delta)$. On figure 11, we detail some analysis on $\phi(\sigma,\gamma)$, and the response to an environmental press $\chi_{int}(\sigma,\gamma)$. See also [8, 9, 23] for an analysis of the system.

Appendix F Linear stability analysis of the One-Equilibrium solution

Starting from the DMFT equation (9), we linearize the system around a Unique Equilibrium. We introduce the relative amplitudes $\delta N(t)=N(t)-N_{\infty}$, $\delta m(t)=m(t)-m_{\infty}$ and $\delta\eta(t)=\eta(t)-\eta_{\infty}$, respectively corresponding to the population, average population and interaction noise. These amplitudes are supposed to go to zero, at least in the Unique Equilibrium phase. The self-consistent relation also holds for these relative amplitudes. Indeed, if we denote $C_{c}(t,s)$ the connected correlator, it verifies:

[TABLE]

From now on, when we write the average $\mathbb{E}$, it will correspond to the average over both the static noise $\eta_{\infty}$ and the dynamical noise $\delta\eta$. The cases when $N_{\infty}=0$ will just give a relaxing exponential and not influence relevantly the correlator nor the response function at large times. Therefore, we will focus on the cases $N_{\infty}>0$, and write the corresponding average $\mathbb{E}_{+}$. The linearization reads:

[TABLE]

We focus on long times, and we assume time-translational invariance for the system: $C_{c}(t,s)=C(t-s)$ and $\chi(t,s)=\chi(t-s)$. This assumption has two consequences. First it transforms the integral term in equation (18) into a convolution product. Secondly, denoting $\tilde{f}(\omega)$ the Fourier transform of $f$, the closure relations become:

[TABLE]

We are computing observables $X$ for alive species only; the global average should be:

[TABLE]

where for relevant observables, the dead species contribution $\mathbb{E}_{dead}[X]$ vanishes in the large-$S$ limit. $\phi$ denotes the fraction of alive species in the Unique Equilibrium around which we are linearizing. It can be computed from the stationary cavity in E.2. Now we send the linearized cavity equation (18) into Fourier space:

[TABLE]

Averaging directly equation (20), and as the perturbation is of zero mean $\langle h\rangle_{h}=0$, we get that the perturbation of the mean population is $\delta m=0$. Finally we apply the relations (19). The three terms $N_{\infty}$, $\tilde{h}$ and $\tilde{\delta\eta}$ verify independence relations for different reasons:

• •

$\tilde{h}$ is independent of $N_{\infty}$ by construction, because we added the perturbation once the steady-state had already been reached;

• •

$\delta\eta$ and $N_{\infty}$ are correlated in the temporal representation, but the Fourier correlation of stochastic variables only stands for the same pulsation. Therefore for any non-zero pulsation $\omega$, $\tilde{\delta\eta}(\omega)$ and $N_{\infty}$ are uncorrelated;

• •

$h$ and $\delta\eta$ are directly uncorrelated, as the noise is sampled from a given covariance $C$.

[TABLE]

Eventually, we end up with the closed forms for the correlator and the response function:

[TABLE]

where the $\phi$ factor appears because we focus on the alive species. Indeed, in order to have normalization, the average $\mathbb{E}_{+}$ corresponds to averaging against the stationary measure which we rescaled: $1/\phi\;p_{+}(N_{\infty})\;H(N_{\infty})\;dN_{\infty}$. Therefore, the absolute contribution of alive species is $\phi\mathbb{E}_{+}$.

Now, as we are interested in the large time behaviour of the system, we perform a small $\omega$ expansion of the equations (21) and (22). This limit needs to be taken carefully because there is a competitive effect in the average between $\omega$ and $\frac{1}{N_{\infty}}$. We finally get the expansions in equations (13) and (14).

Appendix G Details of the numerical strategy for the DMFT solver

We made a gitHub repository with the Python programs we wrote [16]. There is also a runMe.py file which can be directly run in order to produce DMFT solutions and figures such as figure 5 or figure 7. In this section, we write down in details the methodology of the algorithm.

We discretize time in equal units of $dt$ such that $t_{k}=k\;dt$. We also fix the final time we’re interested in as $t_{max}=\#_{time}\;dt$. We usually take $dt=0.1$. The two-dimensional functions then become matrices, and the one-dimensional ones are vectors:

[TABLE]

We will now describe how one iteration of the algorithm is computed numerically. We start from the observables $m_{k}$, $C_{kl}$ and $\chi_{kl}$, and we want to compute the new ones $m^{new}_{k}$, $C^{new}_{kl}$ and $\chi^{new}_{kl}$ after one iteration.

G.1 Sampling of the noise

We will simulate $\#_{traj}$ trajectories that we will refer to as ”species”. Remember that they are independent in DMFT setting. We will then detail the procedure for one species only. For each species, we need a given realization of the gaussian noise at all times $\{\eta_{k=1..\#_{time}}\}=\{\eta(t=0...t_{max})\}$, sampled according to the correlator $C$. Given the discretization, we sample $\{\eta_{k=1..\#_{time}}\}$ as a multivariate gaussian vector with covariance $C_{kl}$. One way to do this is to diagonalize the matrix $C_{kl}$, then in the proper basis all components are independent.

G.2 Numerical integration of the trajectory

For the trajectory of the species, the integration of the differential equation is done with a basic Euler scheme. The Lotka-Volterra system is better simulated in log space. Therefore, if we denote $N_{k}=N(t_{k})$, we implement the scheme:

[TABLE]

with:

[TABLE]

The last $\lambda$-dependent scheme is for numerical stability whenever there is immigration in the system. Using this scheme, we compute the trajectory $N_{k=1..\#_{time}}$.

G.3 Computing the new observables from the trajectories

We sample and integrate the trajectories for $\#_{traj}$ species following the previous procedure. We end up with an array of trajectories:

[TABLE]

From them we can compute the new observables $m$ and $C$ by direct averages:

[TABLE]

As stated in D and section 3.3, the response function $\chi$ is more difficult to compute. The two methods can be used:

[TABLE]

In the last line of the equation $\chi_{kl}^{i}$ is integrated for each species according to C. Eventually, we will start a new iteration of the algorithm, with a soft update:

[TABLE]

After some trials, a reinjection parameter $a=0.3$ is a good value.

G.4 Convergence and the iterative strategy

After one iteration of the algorithm, the new results always present some statistical noise, due to the fact that we average over a finite number $\#_{traj}$ of trajectories. To get a good convergence, we increase the number of trajectories as the iteration goes on. The first iterations are performed with few trajectories; they correspond to rough steps in the configurational space. As the observables get closer to the real solution, we refine the iterations by using more trajectories.

All results are shown with the following scheme: $30$ iterations with $10^{3}$ trajectories each, then $10$ iterations with $10^{4}$ each, and $20$ iterations with $10^{5}$ each. The convergence is considered to be reached when the iteration step becomes lower than a given threshold. More precisely, labeling $C^{i}_{kl}$ the correlator after iteration $i$, we have reached convergence when:

[TABLE]

where $\|M\|_{F}=\#_{time}^{-2}\sum_{kl}M_{kls}^{2}$ is the rescaled Frobenius norm. We use the threshold on the correlator, because we found that it is the most difficult observable to converge. A mixed criterion in all three observables would work as well. On figure 12 we show the convergence in terms of iteration steps.

Appendix H Some examples of numerical solutions

On figure 13, we show an example of a chaotic correlator $C(t,t^{\prime})$. It is to be put in contrast with figure 7, which depicted the Unique Equilibrium plateau type correlator. On figure 14, we show an example of a response numerical solution $\chi(t,t^{\prime})$. The behaviour of $\chi$ does not seem to change drastically between Unique Equilibrium phase and the Multiple Attractors phase.