Aging by Near-Extinctions in Many-Variable Interacting Populations

Interactions between species in ecosystems may lead to large fluctuations in their population sizes. Theoretical models play a central role in understanding these fluctuations in nature and experiments, both for several species [1–4] and for many species [5]. The dynamics of populations that interact and reproduce are often modeled by coupled ordinary differential equations for population sizes {𝑥𝑖}. They are non-negative variables, 𝑥𝑖 ≥0, and must remain so throughout the dynamics. The boundary values 𝑥𝑖⁡(𝑡) =0 represent extinct populations: If a population is extinct at some time 𝑡, it must remain so at all later times. Namely, 𝑥𝑖 =0 is an absorbing value for 𝑥𝑖. These requirements are satisfied by a broad class of differential equations of the form [6]

Examples in this class include the Lotka-Volterra equations for which 𝑔𝑖⁡(→𝑥)=𝐵𝑖−∑𝑗𝐴𝑖⁢𝑗⁢𝑥𝑗, with the matrix 𝐀 encoding the interactions between populations; resource-competition models [7]; and the replicator equations employed in evolution and game theory [6].

It is well known that, depending on the shape of the functions 𝑔𝑖, few variable systems of the form (1) can exhibit different long-time behaviors such as stationarity, periodicity, or chaos [6]. Remarkably, the existence of absorbing hyperplanes has also been shown to lead, in some cases, to robust heteroclinic cycles [8,9]. A classical example is the three-species Lotka-Volterra system with rock-paper-scissors-type interactions [10], where each species hinders the growth of the next. There, trajectories are attracted to a cycle connecting three unstable fixed points, each with a single surviving population; see Fig. 1(a). As time increases, they pass ever closer to these fixed points, resulting in slowdown of the dynamics, with exponentially increasing sojourn times in their vicinity and rapid transitions between them [10,11].

FIG. 1.

In models characterized by a large number 𝑆 of variables, recent works find that an analogous slowdown emerges generically for random interaction coefficients. It is known that aging (a situation in which the system does not asymptotically settle to a fixed point but keeps exploring the phase space with a velocity that nevertheless decays with the elapsed time) can occur in many-variable Lotka-Volterra systems with random asymmetric interactions [13] and replicator equations with nearly antisymmetric random interactions [14]. Here, some populations experience ever-longer periods near extinction ( 𝑥𝑖≃0 and 𝑥𝑖 closer to zero in successive near-extinction periods) before eventually returning to 𝑥𝑖=𝑂⁡(1); see Fig. 1(b). Such dips and “blooms” are documented in experiments and field data (e.g., Refs. [15,16]) and are ecologically significant, as they may lead to extinctions in actual finite populations. The properties of these dynamics have remained elusive, however. The analogy with low-dimensional examples such as in Fig. 1(a) is limited. For one, in the many-variable case, the system does not approach a limit cycle (at least if the limit 𝑆 →∞ is taken before 𝑡 →∞). Second, large dynamical systems of the form (1) may possess many fixed points with different properties (e.g., the fraction of variables for which 𝑥𝑖 =0 or their instability index), and linking the characteristics of fixed points to the dynamics remains an open problem.

In this Letter, we propose a high-dimensional model of the form (1) that provides insights into the connection between aging and absorbing values, by bypassing some of the difficulties inherent to the many-variable Lotka-Volterra and replicator equations. Fixed points of Eq. (1) satisfy either 𝑥𝑖 =0 or 𝑔𝑖⁡(→𝑥)=0 for every 𝑖. Since the unique aging behavior of these systems is tied to the existence of absorbing values ( 𝑥𝑖 =0), we introduce a model with two absorbing values for each variable, which we refer to as the mirrored-extinction model. Specifically, we consider the evolution of 𝑆 degrees of freedom {𝑥𝑖}𝑖=1,…,𝑆, with 0 ≤𝑥𝑖 ≤1 for all 𝑖:

where 𝜶 is a zero-mean Gaussian random matrix with independent and identically distributed entries (referred to as asymmetric interactions). We take 𝔼⁡[𝛼2𝑖⁢𝑗]=1/𝑆, which sets the units of time. From an ecological perspective, the interactions in Eq. (2) affect the growth rates of populations but not their maximal size, which might be limited by other factors; see, for example, Refs. [17,18]. Equation (2) can be extended by adding a species-dependent growth rate 𝑔𝑖 to the sum: ∑𝑗𝛼𝑖⁢𝑗⁢𝑥𝑗→𝑔𝑖+∑𝑗𝛼𝑖⁢𝑗⁢𝑥𝑗 (so that when a species is alone it undergoes simple logistic growth with a growth rate 𝑔𝑖, similarly to the Lotka-Volterra equations) and is solvable just as described below and with the same qualitative outcomes [12].

The resulting dynamical system has many fixed points where all degrees of freedom are at their absorbing values, either 𝑥𝑖 =0 or 𝑥𝑖 =1, allowing us to focus on the effects of these absorbing boundaries. It displays aging, similarly to the Lotka-Volterra case, but with the 𝑥𝑖 spending ever-longer times close to either 𝑥𝑖 =0 or 1 with rapid transitions between these two values; see Fig. 1(b). Importantly, the model in Eq. (2) is exactly solvable in high dimension, allowing us to obtain detailed information on the link between near-extinction processes and aging, beyond other models that also exhibit both phenomena [13,14].

The mechanism for aging found here is drastically different from that at play in aging of usual spin glasses following a quench, where the system’s energy is reduced until it reaches an energy surface dominated by marginally stable fixed points and spends its time there [19–21]. This includes Lotka-Volterra dynamics with symmetric interaction matrices 𝛼𝑖⁢𝑗 [22,23], where →𝑔⁡(→𝑥) is the gradient of a potential, thus permitting a mapping to a spin-glass phase. This form of aging is known to disappear when asymmetry is introduced to the interaction coefficients [24,25].

In contrast, here we show that aging happens in Eq. (2), as variables are driven close to their absorbing values: The probability 𝑃⁡(𝑥𝑖) at long times concentrates about {0,1}, as shown below in Eq. (13). Near fixed points, the dynamics slow down, as manifested in the autocorrelation 𝐶⁡(𝑡′+𝜏,𝑡′) of 𝑥𝑖⁡(𝑡), which, as 𝑡′ grows, relaxes more slowly with 𝜏, as shown below in Eq. (9). Similarly to the three-variable example in Fig. 1(a), typical systems go very close to fixed points which are therefore long lived; see Figs. 1(b) and 1(c). This happens despite these fixed points being unstable, which we show by calculating the spectrum [Eq. (17)] of the linearized dynamics around the fixed points approached at long times. This provides a mechanism for aging in the absence of an underlying energy function. We find that, in the long-time limit, the system moves between infinitely many unstable fixed points that all have the same finite fraction of unstable directions and the same stability spectrum. They are neither the most stable nor the most abundant fixed points.

Dynamical mean field theory.— To analyze the many-variable dynamics (2), we use dynamical mean field theory (DMFT) [26,27]. In the limit 𝑆 →∞ and for 𝑥𝑖 sampled independently at the initial time, the dynamics of a single degree of freedom 𝑥⁡(𝑡) are exactly described by a stochastic differential equation:

with 𝜉⁡(𝑡) a zero-mean Gaussian process. This stems from the fact that the term 𝜉𝑖⁡(𝑡)≡∑𝑗𝛼𝑖⁢𝑗⁢𝑥𝑗⁡(𝑡) appearing in Eq. (2) is the sum of many weakly correlated contributions. As is usual in DMFT, this expression for 𝜉𝑖⁡(𝑡) yields a self-consistent closure relation that reads 𝐶⁡(𝑡,𝑡′)≡⟨𝜉⁡(𝑡)⁢𝜉⁡(𝑡′)⟩=⟨𝑥⁡(𝑡)⁢𝑥⁡(𝑡′)⟩. Here, the angular brackets ⟨.⟩ denote an average over the initial conditions 𝑥⁡(0) and realizations of the noise 𝜉⁡(𝑡). The derivation of Eq. (3) follows a standard procedure [13,26–29] and is detailed in Ref. [12]. To analyze the dynamics, it is, therefore, very helpful to solve for the autocorrelation function 𝐶⁡(𝑡,𝑡′).

To proceed, we introduce the transformation 𝑢⁡(𝑡)=ln⁡[𝑥⁡(𝑡)/(1−𝑥⁡(𝑡))] that sends the boundaries of the domain [0, 1] to (−∞,+∞) and for which Eq. (3) becomes

with the closure relation

where

Aging and the autocorrelation function.— We start by showing that the mean-square displacement of 𝑢⁡(𝑡) is ballistic. Denote the autocorrelation 𝐺⁡(𝑡,𝑡′)≡⟨𝑢⁡(𝑡)⁢𝑢⁡(𝑡′)⟩, which by Eq. (4) is related to 𝐶⁡(𝑡,𝑡′) by 𝐶⁡(𝑡,𝑡′)=∂𝑡∂𝑡′𝐺⁡(𝑡,𝑡′). We take initial conditions such that 𝑢⁡(0) =0 or, equivalently, 𝑥⁡(0) =1/2; the long-time behavior of the correlation function is insensitive to this choice. The closure equation in Eq. (5) can then be written as

𝑢⁡(𝑡) and 𝑢⁡(𝑡′) are jointly Gaussian with zero means, from which it follows that 1/16≤⟨𝑓⁡(𝑢⁡(𝑡))⁢𝑓⁡(𝑢⁡(𝑡′))⟩≤1; see Ref. [12]. Therefore, 𝑡⁢𝑡′/16<𝐺⁡(𝑡,𝑡′)<𝑡⁢𝑡′, so ⟨𝑢⁢(𝑡)2⟩=𝐺⁡(𝑡,𝑡)∼𝑡2, corresponding to ballistic growth of 𝑢⁡(𝑡). We show below that 𝑢⁡(𝑡) nonetheless repeatedly crosses the origin at arbitrarily long times.

The long-time expression for 𝐺⁡(𝑡,𝑡′) can be worked out from Eq. (6). Here, we present a different but equivalent derivation, which makes explicit the aging properties of the model. Motivated by the ballistic growth of 𝑢⁡(𝑡), we introduce 𝑧⁡(𝑡)≡𝑢⁡(𝑡)/𝑡, and we rescale time though 𝑠 ≡ln⁡(𝑡). The resulting dynamics read

together with the closure relation [from Eq. (5)]

Because 𝑧⁡(𝑠) is a Gaussian process with finite 𝑂⁡(1) variance as 𝑠 →∞, in the long-time limit this equation reads

Equations (7) and (8) map the original many-body dynamics of Eq. (2), in the long-time limit, to chaotic dynamics of random neural networks of the form discussed in Ref. [30]. As in Ref. [30], at large 𝑠, we expect the process in Eq. (7) to reach a time-translation-invariant chaotic state characterized by

In the original timescale 𝑡 =𝑒𝑠, this corresponds to autocorrelation of the form

at fixed 𝛽 ≡𝜏/𝑡′. 𝐶⁡(𝑡′+𝜏,𝑡′) thus relaxes more slowly with 𝜏 as 𝑡′ grows, a hallmark of aging, here with correlation time growing linearly with the elapsed time. Accordingly, from Eq. (7), the 𝑧⁡(𝑠) autocorrelation function also admits a time-translation-invariant form at large times:

which is related to 𝐺⁡(𝑡,𝑡′) through lim𝑡′→∞⁡𝐺⁡(𝑡′+𝜏,𝑡′)/𝑡′⁢(𝑡′+𝜏)=^Δ⁡(ln⁡(1+𝛽)). We now sketch the derivation of ^Δ. Following Ref. [30], ^Δ and ^𝐶 are related by ^𝐶⁡(𝑠)=−^Δ′′⁡(𝑠)+^Δ⁡(𝑠), which, together with Eq. (8), implies that ^Δ⁡(𝑠) satisfies an equation for the motion of a classical particle in a potential 𝑉:

where the potential depends parametrically on the initial condition Δ0 ≡^Δ⁡(0) and reads

The condition ^Δ⁡(𝑠) =^Δ⁡(−𝑠) implies ^Δ′⁡(0) =0 so that the ^Δ⁡(𝑠) trajectory has zero initial kinetic energy. The only physically relevant trajectory is, therefore, the one converging to the unstable fixed point Δ* with the same potential energy as the initial condition and related to Δ0 by 𝑉′⁡(Δ*,Δ0)=0 and 𝑉⁡(Δ*,Δ0) =𝑉⁡(Δ0,Δ0). This gives Δ0 ≃0.476 and Δ* ≃0.427.

The correlation 𝐶⁡(𝑡′+𝜏,𝑡′)=(1/𝑆)⁢∑𝑖𝑥𝑖⁡(𝑡′+𝜏)⁢𝑥𝑖⁡(𝑡′), obtained by running the dynamics (2), is thus expected by Eq. (9) to collapse when plotted against 𝜏/𝑡′, as indeed seen in Fig. 1(d), and it matches the correlation function ^𝐶⁡(𝑠) obtained by numerically solving Eq. (10) with the appropriate initial conditions. Note that Δ0 is linked to the long-time growth of ⟨𝑢⁢(𝑡)2⟩, as 𝐺⁡(𝑡,𝑡)/𝑡2→𝑡→∞Δ0. Additionally, the autocorrelation satisfies

so that the system continues to evolve, as the correlation with the state at any time is later partially lost. Equation (10) implies a power law relaxation of 𝐶⁡(𝑡,𝑡′) in the aging regime to its plateau value Δ*:

with 𝑘=√|𝑉′′⁡(Δ*,Δ0)|≃0.492.

Single-variable dynamics.— The dynamics (2) pass very close to fixed points at long times. To see this, we calculate the probability distribution of 𝑥 at time 𝑡, 𝑃𝑡⁡(𝑥), taken over many variables in Eq. (2) or, equivalently, over different realizations of Eq. (3). Using the fact that 𝑢⁡(𝑡) is Gaussian and that 𝑥⁡(𝑡)=𝑓⁡(𝑢⁡(𝑡)), it reads

In particular, this implies

This shows that the system (2) asymptotically approaches fixed points of the dynamics, where all 𝑥𝑖 ∈{0,1}. Furthermore, at large but finite times, the probability to find 𝑥⁡(𝑡) away from the boundaries of [0, 1] decays as 1/𝑡, with Eq. (12) giving

for any fixed 𝜀 ∈[0,1/2]. The probability is, thus, concentrated exponentially close in time to 0 and 1. Yet the system continues to evolve [see Eq. (11)], so that none of these fixed points are stable: At long times, the system transitions between unstable fixed points, spending ever-longer times in their vicinity with fast transitions between them.

In the long-time limit, since 𝑥⁡(𝑠)=Θ⁡(𝑧⁡(𝑠)) for 𝑠 →∞, 𝑥⁡(𝑠) asymptotically approaches a time-translation-invariant two-state process. This is illustrated in Fig. 1(c). Note that, as Δ* >0 in Eq. (11), Eq. (13) is not the ergodic measure (in log time) of a single variable 𝑥𝑖⁡(𝑡). In Ref. [12], we show that for a given degree of freedom the log-time ergodic measure is given by

with 𝑝=[1+Erf⁢(¯𝜉/√2⁢𝜒)]/2, where ¯𝜉 is a zero-mean Gaussian random variable with variance ⟨¯𝜉2⟩ =Δ* and 𝜒=∫∞0𝑑𝑠⁢𝑒−𝑠⁢[^𝐶⁡(𝑠)−Δ*]. So, in a given realization of Eq. (2), each variable has an “identity” expressed in the fraction of time (in log time) it spends near 0 and 1.

Stability of visited fixed points.— We found above that at long times the system approaches fixed points but eventually leaves their vicinity, signaling that they are unstable. We now calculate their entire stability spectrum. The linearized dynamics close to a fixed point 𝒙* are ˙𝛿⁢𝑥𝑖=𝐽𝑖⁢𝑗⁢𝛿⁢𝑥𝑖 with a diagonal matrix 𝐽𝑖⁢𝑗=𝛿𝑖⁢𝑗⁢𝜆*𝑖. The growth rates 𝜆*𝑖, positive when growing in the direction away from the boundaries, are given by

The stability spectrum of the visited fixed points is therefore equal, at long times, to the empirical distribution in the many-variable dynamics (2) of 𝜆𝑖⁡(𝑡)≡[1−2⁢𝑥𝑖⁡(𝑡)]⁢𝜉𝑖⁡(𝑡) for 𝑖 =1…𝑆. In the 𝑆 →∞ limit, the stability spectrum is thus equal to the distribution of 𝜆⁡(𝑡) =[1−2⁢𝑥⁡(𝑡)]⁢𝜉⁡(𝑡) in the DMFT framework. It can also be shown that the 𝜆𝑖⁡(𝑡) are independent and identically distributed random variables [12]; therefore, the spectrum is self-averaging.

The joint distribution of 𝜉⁡(𝑡) and 𝑢⁡(𝑡) is Gaussian, with correlations ⟨𝑢⁢(𝑡)2⟩=𝐺⁡(𝑡,𝑡), ⟨𝜉⁢(𝑡)2⟩ =𝐶⁡(𝑡,𝑡), and cross-correlation ⟨𝑢⁡(𝑡)⁢𝜉⁡(𝑡)⟩. Changing variables from (𝑢,𝜉) to (𝑢,𝜆) and integrating over 𝑢, we obtain an expression for the distribution of 𝜆⁡(𝑡), reproduced in Ref. [12]. Taking its long-time limit, we find that the dynamics (2) transition between fixed points which all have the same stability spectrum:

with 𝜅∞=√1/(2⁢Δ0)−1≃0.224; see Fig. 2. This distribution has a finite fraction of unstable directions, given by

Thus, the system approaches unstable fixed points. This can be compared with the statistics of the full distribution of fixed points of Eq. (2) with all 𝑥𝑖 ∈{0,1}. There are 2𝑆 of them, and the average number of those with 𝛼⁢𝑆 unstable directions is given by the binomial law ⟨𝒩𝛼⟩∼exp⁡(𝑆⁢𝑔⁡(𝛼)), with 𝑔⁡(𝛼)=−𝛼⁢ln⁡𝛼−(1−𝛼)⁢ln⁡(1−𝛼). Therefore, in typical fixed points, half of the directions are unstable: 𝛼 =1/2. The dynamics, therefore, selects in the long-time limit fixed points that are exponentially rare (compared to the typical ones) but that are not the most stable ones existing in the phase space, which are marginal ( 𝛼 =0).

FIG. 2.

In conclusion, we propose an exactly solvable many-variable model for the dynamics of interacting populations with absorbing boundary values. Its dynamics slow down with a correlation time that grows as the age of the system; see Eq. (9). The system evolves in the vicinity of fixed points: In the long-time limit, all variables are found exponentially close in time to absorbing values; see Eq. (14). The time it takes for a variable to leave the vicinity of one absorbing value to visit the vicinity of the other is, therefore, proportional to the age of the system. This explains the scaling of the aging [Eq. (9)]. All these fixed points are unstable, as shown in Eq. (17), in contrast with marginal fixed points reached in usual glassy dynamics [19–21]. In the future, it would be interesting to understand how this scenario is adapted to other many-variable interacting population dynamics, such as the Lotka-Volterra model, where fixed points have degrees of freedom that are not at absorbing values. Fingerprints of these phenomena might be observed, as an increase in correlation time combined with population blooms, in experiments that follow interacting species starting from similar population sizes.