Competitive exclusion and Hebbian couplings in random generalised Lotka-Volterra systems

We study communities emerging from generalised random Lotka--Volterra dynamics with a large number of species with interactions determined by the degree of niche overlap. Each species is endowed with a number of traits, and competition between pairs of species increases with their similarity in trait space. This leads to a model with random Hopfield-like interactions. We use tools from the theory of disordered systems, notably dynamic mean field theory, to characterise the statistics of the resulting communities at stable fixed points and determine analytically when stability breaks down. Two distinct types of transition are identified in this way, both marked by diverging abundances, but differing in the behaviour of the integrated response function. At fixed points only a fraction of the initial pool of species survives. We numerically study the eigenvalue spectra of the interaction matrix between extant species. We find evidence that the two types of dynamical transition are, respectively, associated with the bulk spectrum or an outlier eigenvalue crossing into the right half of the complex plane.


I. INTRODUCTION
The foundations of the theory of disordered systems date back close to 50 years [1].Initially, the aim was to understand certain magnetic states in condensed matter physics ('spin glasses') [2].However, it became clear that applications of the tools developed for disordered systems had a reach far beyond the boundaries of physics.Methods such as replica theory or dynamic generating functionals were quickly adapted and used to answer questions in neural networks [3][4][5], to study the Minority Game [6] (sometimes presented as a simple model of a financial market), or indeed evolutionary bi-matrix games and socalled Nash equilibria [7].
The defining feature of disordered systems is the presence of quenched disorder.That is, the system is made up of many constituents, and the interactions between these are determined by coefficients that are drawn at random at the beginning, but then remain fixed as the dynamics of the system unfolds.The disorder leads to complicated energy landscapes.The number of local minima can grow exponentially in the size of the system, and is often organised in a hierarchical manner.Dynamic phenomena in disordered systems include ergodicity breaking and so-called ageing [2,8].
Ideas and methods from the physics of disordered systems have also been used to study complex ecosys-tems [9][10][11][12][13][14][15][16][17].The word 'complex' in this context indicates that the ecosystem is composed of a large number of species, and that these species are subject to randomly drawn interaction coefficients.In this paper we continue this line of work, and focus on a Lotka-Volterra system with Hopfield-like interactions [3][4][5].More specifically, we are interested in a set of N species (N ≫ 1), whose abundances develop in time following a generalised Lotka-Volterra equation (details will follow in Sec.II).This involves an N × N matrix a ij of interaction coefficients.Existing work on the statistical physics of complex ecosystems has mostly focused on the case in which the interaction matrix is drawn from distributions with either no correlations between different matrix elements, or only correlations between diagonally opposed entries [9,11,[14][15][16][17][18][19].There is also work on cases in which the matrix is composed of blocks, and where the elements in different blocks have different statistics [20].One common element shared by many existing random Lotka-Volterra models is that the finest level of modelling is set by the interaction coefficients.No further assumptions are made about the properties of the species, and how the species interactions come about from these properties.
The Hopfield model is inspired by structures first used in neural networks [5,21,22].Translated into the language of ecology, the starting point is now a set of species and a set of traits.Each species can either possess or not possess a given trait.This assignment of traits to species, in turn, determines how species will interact.Broadly speaking, the interaction between two species will be more competitive the more traits they share (i.e. the more similar the two species are).This type of interaction structure has also been studied in models combining resources and consumers, both in economics and in ecology [12,13,[23][24][25].A particularly notable model is that by MacArthur and collaborators [26][27][28].Analyses of random replicator systems with 'Hebbian' interactions [29] have shown interesting statistical mechanics, and in particular types of phase transition that are different from those seen in replicator systems with Gaussian couplings.
In this paper, we set out to characterise the behaviour of a Lotka-Volterra system with Hopfield-like interactions, where we allow for a degree of 'mild' dilution (the system is not fully connected, but each species still interacts with an extensive number of other species).A system of replicator equations with such interactions was studied in [29].Our aim is to calculate the statistics of fixed points in the phase where such fixed points are attained and identify the onset of instability.As in the system with Gaussian interactions, we find that only a proportion of the initial species survive at stable fixed points.Recent work [30] on Gaussian systems has shown that the reduced interaction matrix (the matrix of interaction coefficients among the surviving species) has intricate statistics.Specifically, its bulk and outlier eigenvalues can be related to different types of dynamic phase transitions.As we will show, the types of phase transition seen in our model differ from those in the Gaussian model.One aim of the current paper is therefore to establish (in simulations) how these transitions relate to the spectra of the interaction matrix of the extant species.
The remainder of the paper is organised as follows.In Sec.II we define the model and introduce the necessary notation.Sec.III then contains the mathematical analysis.This is based on so-called 'generating functionals' and dynamic mean field theory.The phase diagram and further behaviour of the model are then discussed in Sec.IV.In Sec.V we finally turn to a study of the spectra of the reduced interaction matrix and their relation to the phase diagram.We conclude the paper with a discussion and an outlook in Sec.VI.

II. MODEL DEFINITIONS
We will study the following generalised Lotka-Volterra equation (gLVE) where the x i ≥ 0 represent the abundances (or population densities) of different species, i = 1, . . .N .We always assume initial conditions for which all x i are strictly positive.
The quantities u i > 0 denote the strength of intraspecific competition, and the a ij = c ij J ij represent the interspecific interactions.The K i (together with the u i ) set the carrying capacities of the species in the absence of interactions between different species (x i then tends to K i /u in the long run).We focus on the case u i ≡ u for all i, noting that u controls the time scale on which the non-interacting system approaches the fixed point x i ≡ K i /u.We allow for general positive values of u throughout our analysis, but in an effort to keep the number of parameters manageable we set K i ≡ 1.
The dilution variables c ij ∈ {0, 1} (i ̸ = j) determine which species interact with one another, i.e. they set the topology of the interaction network.For each pair i < j, the coefficients c ij and c ji are chosen from a Bernoulli distribution with We thus have P (c ij = 1) = c for all i ̸ = j, i.e. c is the analog of what May called 'connectance' [9].The parameter Γ is restricted to the range from −1 to 1 by construction, but we note that not all choices of pairs (c, Γ) are possible (see Supplemental Material [31] for details).We note that the choice of the diagonal coefficients c ii is irrelevant as we set J ii = 0 below Throughout our paper, c is chosen not to scale with N [c = O(N 0 )].This means that each species interacts with an O(N ) number of other species.We are therefore not studying a 'dilute' system in the sense of random matrix theory.The extensive connectivity allows us to use established methods from dynamic mean-field theory.Truly dilute systems with c = O(1/N ) (and where consequently each species only interacts with a finite number of other species) can be expected to behave very differently, see e.g.[32], and a theoretical analysis would be much more intricate.
Interaction links in our system are directed, that is, an effect of the presence of species j on the dynamics of i does not necessarily imply the reverse.The parameter Γ measures the correlations between c ij and c ji .A choice of Γ = −1 implies c ij = 1 − c ji with probability one, and Γ = 1 means that c ij = c ji with probability one.
The matrix J ij determines the strength of the effects of the presence of species j on the dynamics of the abundance of species i. Positive values of J ij imply that the population of species j is beneficial to the growth of species i, while negative values imply a detrimental interaction.In ecological terms the signs of the pair of interactions (sgn J ij , sgn J ji ) determine whether two species are in a mutualistic relation (+, +), whether they compete with one another (−, −), or whether there is an antagonistic predatorprey relation between i and j (±, ∓) [10].
In this work, the interspecific interaction is chosen according to a niche overlap heuristic (see e.g.[33]).
We assume each species is described by a set of binary traits, labelled µ = 1, . . ., P , and that a pair of species will compete in proportion to the similarity between the two species (i.e, the number of traits which both or neither species possess).We write ξ µ i = +1, if species i has trait µ, and ξ µ i = −1 if the species does not possess the trait.Interactions are then assumed to be of the form We have here set P = αcN , with α > 0 a model parameter (in simulations P is restricted to integer values.)That is to say, we assume that the number of traits is proportional to the number of species in the system.The interaction in Eq. ( 3) is reminiscent of the Hopfield model, used in the context of neural networks [4,21].This suggests interesting phase behaviour when α = O(N 0 ), which is the regime we focus on.We have normalised the interaction strength by cN , the mean number of species that any one species will interact with.We will refer to the random variables c ij and ξ µ i as the disorder of the system.
The traits ξ µ i are chosen to be ξ µ i = ±1 with equal probability, and there is no correlation between the different ξ µ i .This implies that the distribution of the J ij approaches a Gaussian as N → ∞, reminiscent of the model studied for example in [14,15].We note, however, that the Hopfield structure introduces correlations between the different J ij , which are different from the correlations studied in the earlier literature.We highlight again the structural similarity to MacArthur's consumer-resource model [26][27][28] (see also [24,25] for statistical physics studies), noting though that the latter model is more sophisticated, with dynamical equations both for consumers and resources.

III. GENERATING FUNCTIONAL ANALYSIS AND STABILITY
A. Generating functional and effective process We analyse the system in Eq. (1) using dynamic generating functionals, an established method in the theory of disordered systems [4,34].This leads to an effective 'dynamic mean field theory'.Similar approaches have been used to study Lotka-Volterra systems with Gaussian random couplings [11,15,35].We note that an alternative approach is based on the so-called cavity method [14,18,36].We also add that the dynamics admit a Liapunov function when c = 1 (leading to a symmetric interaction matrix).Methods from the equilibrium statistical physics of disordered systems such as the replica approach can be used in this special case (for examples see [16,17]).
The outcome of the application of these techniques is an effective stochastic process for a 'representative species'.The ensemble of realisations of stochastic processes is statistically equivalent to the set of single-species trajectories x i (t) of the disordered dynamical system in Eq. ( 1).The dynamic mean-field description becomes exact in the thermodynamic limit (N → ∞).Overall, in this limit, the infinite-dimensional deterministic dynamical system in Eq. ( 1) is traded for an effective single-species process which is non-local in time (it involves retarded self-interaction) and contains coloured noise.
The generating functional analysis begins from where we have introduced the perturbation fields h i (t) in order to calculate linear response functions.These fields are not actually part of the model and are set to zero at the end of the calculation, as well as in all simulations shown in the paper.For more details see the Supplementary Material (SM).The generating functional of this dynamical system is given by where the average is over paths [x 1 (t), . . ., x N (t)] of the dynamics in Eq. ( 4).The ψ i (t) constitute a source field.The generating functional in Eq. ( 5) is the Fourier transform of the probability measure in the space of paths generated by Eq. ( 4).
The final outcome of the generating-functional analysis is a set of equations for the dynamic macroscopic order parameters of the problem.For the Lotka-Volterra model these are where δ denotes a functional derivative and ⟨• • • ⟩ 0 stands for an average over random initial conditions.The overbar • • • represents the average over the disorder, i.e. over the c ij and ξ µ i .The order parameters can be obtained from the disorder-averaged generating functional as derivatives with respect to the fields ψ i (t) and/or h i (t), evaluated at ψ i (t) ≡ 0 and h i (t) ≡ 0. The order parameters in Eqs. ( 6) are determined self-consistently from an effective process for a single representative ('mean field') species.The procedure to derive the effective equations is well-documented [6,11,37,38], therefore, we only report the final result (a more detailed derivation can be found in the SM).The effective single-species process for the model is given by where I is the identity operator and η(t) is coloured Gaussian noise with zero mean and correlations in time, given by The order parameters in Eqs. ( 6) are to be obtained self-consistently from the following expressions, where the average ⟨• • • ⟩ * is performed over realisations of the process in Eq. (7).Eqs. ( 7)-( 9) form a closed system and have to be solved self-consistently.

B. Fixed point analysis
There is no realistic prospect for a general analytical solution of the effective dynamics in Eq. ( 7).One alternative is to use Monte-Carlo methods to construct sample paths for the effective process and solutions for the dynamic order parameters.For example via the Eissfeller-Opper procedure [39], or using the more recent approach in [18].The latter reference explicitly discusses applications to random Lotka-Volterra systems with Gaussian disorder.
Here we will instead follow [11,15,35] and focus on analytical solutions in the parameter regime in which the dynamics approach stable fixed points.This is motivated by observations from the numerical integration of Eq. ( 1).We find that, for certain parameters, the system tends to a unique fixed point, which is independent of initial conditions.Fig. 1 shows examples of parameter regions in which this is the case.Broadly speaking, we observe two different types of behaviour: (i) the population densities converge to a fixed point, or (ii) they diverge.These types of behaviour occur in different regions of parameter space (Fig. 1).There is a thin boundary between the two regions where other behaviour (e.g.periodic behaviour or persistent irregular motion) can appear, as evidenced by the occasional green or light blue pixel in Fig. 1.We attribute this to the fact that the system size N is necessarily finite in numerical experiments, and we expect that this behaviour will become increasingly rare as N → ∞.
We will thus assume that each path in the ensemble of trajectories of the effective process eventually arrives at a unique fixed point, x = lim t→∞ x(t).Each realisation of the noise variable η(t) in Eq. ( 9) also approaches a stationary value η.We note that x and η will be random variables, differing across realisations of the effective dynamics.We can then write These relations can be understood as follows: if all realisations of the effective dynamics approach stationary values then M (t) will approach a constant, given by ⟨x⟩ * .Furthermore, we assume that the response function G(t, t ′ ) becomes time-translation invariant for large t, i.e.G(t, t ′ ) = G(t − t ′ ).Causality implies that G(t − t ′ ) = 0 for t < t ′ .Finally, given that all trajectories of the effective dynamics approach fixed points, the correlation function C loses all time dependence and so we have written C(t, t ′ ) = q.This is consistent with Eq. ( 8), the noise variables η(t) also approach a random but timeindependent value for all realisations.The mean of the random variable η is zero and using Eq. ( 8), its variance is given by where From now on, we will write η = √ q Σz, where z is a standard Gaussian random variable, and Setting the time derivative on the left-hand side of Eq. ( 7) to zero, and using Eqs.( 10)-( 12) we find For a given value of z this is to be solved for x, subject to the constraint that abundances are nonnegative, i.e, x(z) ≥ 0. Irrespective of the value of z, Eq. ( 15) always has the solution x(z) = 0. Additionally, a second non-negative solution is possible for some values of z.As we will confirm in simulations, the physically meaningful solution is given by For given order parameters q and χ, the function x(z) in Eq. ( 15) is therefore piecewise linear, with one piece equal to zero.The denominator in Eq. ( 15) always comes out positive.Therefore, we have the solution x(z) > 0 when z < ∆, and x(z) = 0 when z ≥ ∆, with Given that z is a Gaussian random variable, the abundances of extant species at the fixed point follow a clipped Gaussian distribution.This is similar to what was reported in other random Lotka-Volterra models, see e.g.[14].An explicit example of a species abundance distribution can be found for instance in [15].
Using this fixed point ansatz, the relations for the order parameters in Eq. ( 10) can be written in the following form [14,15] where It is now convenient to introduce the following functions for n = 0, 1, 2. We then find from Eqs. ( 17) Eqs. (19) together with Eq. ( 16) form a closed set for the set of unknowns q, χ, M and ∆, which is to be solved as a function of the model parameters u, c, Γ and α.
Recalling that x(z) > 0 if, and only if, z < ∆ we identify f 0 (∆) as the fraction of surviving species, Eqs. ( 19) can be solved parametrically.We fix u, c, Γ and ∆ and then solve for the set of χ, q, M and α.
In detail, we find the following cubic equation for χ, valid for c < 1, Further, we have from Eqs. (19), where the f n are to be evaluated at ∆.The relations in Eqs. ( 21) and ( 22) are also valid for c = 1, and can then be simplified as outlined in Appendix B.
The validity of the predictions from Eqs. ( 21) and ( 22) is confirmed by direct numerical integration of the gLVE in Fig. 2.

Diverging abundance
Model with c < 1.The first and second relations in Eqs.(22) indicate that the order parameters M and q both diverge in the system with c < 1 when f 0 (∆) = f 2 (∆).The latter implies ∆ = 0.The value of α for which this occurs can (for a given choice of c, u and Γ) be obtained from the third relation in Eqs.(22), with χ being the relevant root of Eq. ( 21).Using Eq. ( 21) the susceptibility χ is found to remain finite at the transition.We note that f 0 (∆) > 0 for all relevant values of ∆.
Model with c = 1.The fully connected system also shows two types of divergences: (i) The quantities M and q both diverge when f 0 (∆) = f 2 (∆), see Eqs. (B1).The susceptibility then remains finite; (ii) Eqs.(B1) further indicate, that M and q also diverge in the model with c = 1 when f 2 0 (∆) = uf 2 (∆).This latter condition results in α = u.From Eqs. (B1) the susceptibility χ is then seen to diverge as well (the divergences of M , q, and χ take place simultaneously).
We note that the divergencies resulting from f 0 = f 2 and f 2 0 = uf 2 can take place at different locations in parameter space for the model with c = 1.If this is the case, and starting in the stable phase, the divergence that occurs first will determine the loss of stability in the fully connected system.For u < 1/2 the transition of type (ii) takes place first as α is increased (M , q, and χ diverge), and for u > 1/2 the transition of type (i) is instead observed (q, M diverge, χ remains finite).At present we do not have any further intuition regarding any significance or special role of the value u = 1/2.

Linear instability
The system also shows a linear instability which can be identified using the procedure established in [11,40].We write x(t) = x + y(t) and η(t) = √ q Σz + ζ(t), where y(t) and ζ(t) are small perturbations about the fixed point of the trajectories of the effective process in Eq. (7).Expanding to first order in these perturbations we find that with We also have the self-consistency relation where D (t, t ′ ) = ⟨y(t)y (t ′ )⟩ * .When x = 0, Eq. ( 23) becomes Eq. ( 15), together with the observation that the denominator in this equation is strictly positive, implies that 1 − √ q Σz < 0 when x = 0.This allows us to conclude that perturbations on extinct species decay, and do not contribute to any linear instability.
For fixed points x > 0 we find from Eqs. ( 15) and ( 23) that To identify the onset of linear instability we follow [11,40].We move to Fourier space, writing ω for the variable conjugate to time t, and using tildes to indicate Fourier transforms.
Focusing on the mode with ω = 0 and following steps similar to those in [11,15,40] we then find from Eq. ( 26) The left-hand side is manifestly non-negative, so a change of sign of the expression inside the square bracket on the right-hand side indicates an inconsistency (and divergence of |ỹ(0)| 2 ).Using Eqs.(19) this is shown to occur when For c < 1 the expression in the square brackets is never zero.This leaves us with the condition f 0 = f 2 , which is the same as we obtained for the divergence of M and q.If c = 1, the term in the square bracket is zero if χ → ∞, which using Eq.(B1) we can write as f 2 0 − uf 2 = 0. From this, we conclude that in our model the linear instability is always accompanied by the instability with diverging mean abundance.This is markedly different from the behaviour of the gLVE model with Gaussian random interactions.In this Gaussian model there are instances where the linear instability sets in as the variance of interactions is increased, but where abundances remain finite and the divergence only occurs at a later point at even higher variance of the interactions.This leads to a phase with multiple attractors between the two transitions [14,16,36].Our analysis indicates that the model with Hopfield-like couplings does not have such a multiple-attractor regime.

IV. PHASE DIAGRAM AND FURTHER BEHAVIOUR OF THE MODEL
A. Phase diagram for the fully connected system (c = 1) The phase diagram of the fully connected model is shown in Fig. 3(a).We recall that, for c = 1, the only model parameters are the self-interaction coefficient u and the ratio of the number of traits to the number of interspecies interactions in the original pool (α = P/cN ).For a fixed value of α, the system shows a unique stable fixed point for u > u c (α), where u c (α) marks the onset of instability.The line in Fig. 3(a), obtained from Eq. (B1), shows the phase boundary between the stable and unstable regions.At this boundary M and q diverge, and if u c < 1/2 we also observe a divergence of χ.
The two types of trajectory in the stable and divergent phases are illustrated in the right panel of Fig. 3(b).In the stable phase the system reaches a fixed point, for any one realisation of the interaction matrix (two examples are illustrated in green and red respectively).
Fig. 3(b) also shows two examples in which the species abundances diverge (blue and orange).The divergence occurs at a finite time.We will discuss this further in Sec.IV C.

B. Phase diagram for connectivity c < 1
Fig. 4 shows how the phase diagram for the system with c < 1 depends on the connectivity c and the symmetry parameter Γ.In all cases there is a single phase boundary, where the divergence of M and q and the onset of linear instabilities coincide.This phase boundary separates a region where trajectories converge to a single globally stable fixed point (phase to the right of the line), from a region where trajectories are unbounded and diverge in finite time (phase to the left).
The phase diagrams in Figs. 3 and 4 show that the system is in the stable phase for small values of α (i.e. a small number of traits relative to the number of species in the initial pool), or large values of u (i.e.large negative self-interaction).This is the consequence of two competing effects, the selfinteraction (parametrised by u) which stabilises the system, and the interaction between species (induced  The system is stable to the right of the line, abundances diverge on the left. by competition of similar species) which promotes instability.When u is large and/or α is small, the stabilizing effect of the intra-species interaction dominates over the interactions across species.In the extreme limit α = 0 (no interaction between different species), each abundance follows a separate logistic equation, ẋi = x i (1 − ux i ), and converges to x i = 1/u.When α is small but non-zero, the system consists of weakly interacting species.The effect of the interactions between species is then a small perturbation to the logistic behaviour of individual species, and does not change the convergence to a fixed point.This can be confirmed from Eqs. ( 21) and ( 22) by taking the limit α → 0, which results in all species surviving with fixed point abundance x * i = 1/u (ϕ → 1, M → 1 u , and Var[x] → 0).A similar result is obtained for u ≫ 1 at fixed value of α.
Conversely, for low values of u or large values of α the system is unstable.In this situation, the stabilising self-interaction is not sufficient to overcome the destabilising effect of the random interactions be-tween species.
The most interesting behaviour takes place at the phase boundary, where the effect of the intraspecific and interspecific nonlinearities are of comparable magnitude.From Eqs. ( 21) and ( 22) we can conclude that ϕ → 1/2, and M ∼ (α c − α) −1 as the system approaches the instability (from the stable phase).Further details can be found in Appendix C.
We further note that decreasing the value of the symmetry parameter Γ, increases the range of the stable region in the phase diagrams in Fig. 4.This is similar to the effect of increasing the fraction of predator-prey interactions in Lotka-Volterra models with Gaussian interactions [10,14,15].Indeed, the effect of a reduction of Γ is to increase the fraction of species pairs i, j with c ij = 1 and c ji = 0, that is the proportion of uni-directional interactions.
Interestingly, the effect of varying the 'connectance' c is not straightforward.As can be seen in Fig. 4 an increased connectivity can, depending on the other model parameters, turn a previously stable system into an unstable one, or vice versa, stabilise a previously unstable system.

C. Finite-time divergence of the mean abundance
As mentioned earlier, the divergence of the abundances in the divergent phase occurs at finite time.This has previously been reported in the model with Gaussian interactions [18], and can be justified heuristically from the Lotka-Volterra equations.Indeed, Eq. ( 1) has a second-order non-linearity in the abundances x i .This can lead to dynamics of the form ẋ ∼ x 2 , which in turn implies a solution of the form x(t) = (c − t) −1 , where c is an integration constant.This results in a divergence at finite time.Fig. 5 shows the time, t div , at which the divergence occurs for different choices of the model parameters.This time grows as one approaches the stability line (from inside the unstable phase).When the stability line is crossed (into the stable phase), the time-todivergence diverges itself (t div → ∞), i.e. the divergence no longer occurs.Results from the numerical integration of the gLVE suggest that the divergence of the abundances is of the form M ∼ (t div − t) ν , where ν ≈ 1, as shown in Fig. 6.This behaviour appears to be independent of initial conditions, the values of the parameters u and α, and the initial number of species N .

V. REDUCED INTERACTION MATRIX AND ITS EIGENVALUE SPECTRUM
Ref. [30] recently established a close connection between different instabilities in the Gaussian random Lotka-Volterra model and the eigenvalue spectrum of the interaction matrix of the surviving species.More specifically, the spectrum of this reduced interaction matrix is composed of a bulk region and a potential outlier eigenvalue.As parameters are changed (starting from within the stable phase) either the bulk spectrum or the outlier eigenvalue can cross into the right half of the complex plane.In the Gaussian model, the crossing of the outlier is associated with a transition marked by the divergence of abundances, and a crossing of the bulk with a linear instability.
In this section, we explore in numerical simulations how the different transitions in the gLVE model with Hopfield-like interactions relate to the eigenvalue spectrum of the matrix of interactions between surviving species.

A. Spectrum of the original interaction matrix
Before we discuss the spectra of the reduced interaction matrix, we make a few remarks on the initial interaction matrix α ij = c ij J ij among all species.Throughout this section we set the diagonal elements of this matrix to zero, the only effect of self-interaction (the term −ux i ) is a simple shift of this spectrum.In the large-N limit the central limit theorem applies to J ij = − 1 cN αcN µ=1 ξ µ i ξ µ j , so each off-diagonal entry α ij of the interaction matrix is either a Gaussian random variable (if Calculating the correlations between pairs of elements we obtain where we have used Eq. ( 2) and the fact that J ij is symmetric.This means that only diagonally opposed pairs of elements are correlated, and that their correlation is determined by both, Γ and c.
Based on a theory that only takes into account correlations between diagonally opposed matrix entries, one might then expect an elliptic spectrum [41], with support given by the ellipse with τ = Γ(1−c)+c.However, as illustrated in Fig. 7, this is an approximation to the true spectrum at best for large values of α.For intermediate values of α (an example is shown in orange in the figure), the eigenvalue spectrum appears to have a triangular shape, and for small values of α (shown in green), the spectrum becomes even more skewed and eventually appears to consist of two separate components (example shown in red).While we cannot fully exclude finitesize effects (the spectra in Fig. 7 are for N = 5000), we believe that the deviations from an elliptical spectrum in Eq. ( 31) are due to higher-order correlations between entries of the interaction matrix.For example, it has been shown in Ref. [42] that cyclic correlations can result in eigenvalue spectra with shapes similar to the ones in Fig. 7.For c = 1 the interaction matrix is a (scaled and shifted) Wishart matrix, so its spectrum follows the Marcenko-Pastur law.

B. Eigenvalues of the reduced interaction matrix
We now conclude the analysis of the model with a numerical study of the spectra of the reduced interaction matrix, that is, the interaction matrix between species that survive at the fixed points of the gLVE.Fig. 8 shows the spectra of this matrix for the case c < 1, and for a choice of Γ less than one.This means that the initial interaction matrix is not symmetric.The reduced matrix is not symmetric either, and as a consequence its eigenvalues will generally be complex.As seen in Fig. 8, the spectrum is not elliptic, and we have found no evidence of an outlier eigenvalue in this scenario (c < 1).In the figure we have fixed u, and varied α.The data suggests that the phase transition at α = α c (u) coincides with the point at which the right-most bulk eigenvalue crosses the imaginary axis into the right half-plane.
In Fig. 9 we study the fully connected system for two different values of the self-interaction strength u.The original interaction matrix in the fully connected model is symmetric by construction, and so is the reduced interaction matrix.As a consequence all eigenvalues are real.
Panel (a) focuses on the case u > 1/2.We find no signs of outlier eigenvalues, and again the data indicates that the transition to instability occurs when the leading bulk eigenvalue crosses into the positive half of the real axis.
Panel (b) shows a scenario in which u < 1/2.In contrast with the situation in (a), an outlier eigenvalue now becomes apparent, and the transition to instability in the gLVE at α = α c (u) now appears to coincide with the point at which the outlier becomes positive.
The connection between the different types of transition and the behaviour of the spectrum of the reduced interaction matrix is shown in Table I.We recall that the mean abundance and the second moment of the abundances diverge at all transitions, and that the onset of the linear instability always coincides with the point of diverging abundances.There are thus only two types of transition, one in which the susceptibility remains finite (χ < ∞), and another for which it diverges (χ → ∞).The table indicates that the former transition (χ finite) appears to coincide with the bulk spectrum of the reduced matrix crossing into the right half of the complex plane.The transition at which χ → ∞ (along with the divergences of M and q) on the other hand seems to be seen when the outlier eigenvalue of the reduced matrix in the fully connected system reaches the origin.
We stress that these are numerical observations, and that these findings should therefore be seen mostly as conjectures at this point.In principle, the spectrum of the reduced interaction matrix can likely be calculated in our model, adapting the method used in Ref. [30].However, this involves a substantial calculation and is beyond the scope of the current paper.

VI. DISCUSSION
To summarise, we have carried out a generating functional analysis of a random generalised Lotka-Volterra system with interactions determined by niche overlap.Species interactions in the model are governed by Hopfield-like couplings subject to mild dilution (the remaining connectivity is still extensive).We have computed the statistics of surviving species in the stable fixed point phase, and we have analytically determined the onset of instability.Similar to the gLVE with Gaussian interactions, asymmetry in the connectivity matrix promotes stability.c < 1 q, M diverge bulk spectrum χ remains finite crosses axis c = 1 u > 1/2 q, M diverge bulk spectrum χ remains finite crosses axis u < 1/2 q, M, χ outlier eigenvalue all diverge crosses axis (at u = α) That is to say, the system becomes more stable when there is a larger fraction of unidirectional interactions (c ij = 1, but c ji = 0).In contrast with the Gaussian model, the linear instability against small perturbations cannot be separated from an instability at which species abundances diverge.As a consequence, there is no phase with multiple stable fixed points our model.Despite some common features, the statistical mechanics of the Gaussian and Hopfield-like models are therefore rather distinct.
Our analysis shows two types of transitions to divergent abundances, one in which the integrated response χ remains finite, and another in which χ diverges.This raises interesting questions about the exact nature of memory onset in the system (a diverging integrated response indicates persistent memory of perturbations).Future work could focus on the precise shape of the response function, where the numerical methods in [18] might prove particularly useful.Given that the fully connected system has symmetric couplings it would also be interesting to see how crossing each of the different types of transition affects the energy landscape.A natural approach here might be the replica method and suitable levels of replica symmetry breaking [16,19].
Numerical simulations provide evidence that the transition at which the integrated response remains finite (χ < ∞) is associated with the bulk spectrum of the reduced interaction matrix (the matrix of interactions between extant species) crossing the axis.The transition at which χ diverges on the other hand appears to be signalled by an outlier eigenvalue crossing the imaginary axis.
These findings in simulations reinforce the intriguing analytical result obtained recently in [30].Namely, the eigenvalues of the interaction matrix in the community of surviving species can be used to decide the stability of feasible equilibria, that is to say, fixed points with non-negative species abundances.In the traditional approach to ecosystem stability by Robert May [9], based on the eigenvalue spectra of random matrices, no actual dynamics are specified, and the feasibility of the assumed equilibria remains unclear.Any fixed point of the generalised Lotka-Volterra model on the contrary is feasible by construction.The study of the spectra of reduced interaction matrices resulting from Lotka-Volterra dynamics can therefore contribute to establishing how May's approach can be adapted to include feasible equilibria.Noting that previous work [14,30] has shown that the statistics of the reduced interaction matrix in random Lotka-Volterra models can be quite different from those of the original interaction matrix, it would be interesting to study the statistics of the J ij among survivors in the present model in more detail.In particular, the niche overlap between surviving species.
On a broader level, our study highlights two common facets of work on the statistical physics of complex systems, which were also seen for example 30-40 years ago when physicists studied neural networks, or 15-20 years ago when a number of physicists worked on the Minority Game.On the one hand, tools from physics can make a difference in problems from other disciplines.In our system (and other models of complex ecosystems more generally) this is the study of feasible equilibria with methods from spin glass physics.At the same time, studying problems arising in other areas can reveal new types of physics and complexity, which one would perhaps not find within the strict boundaries of traditional physics.In our case, these are the different types of phase transition in the generalised Lotka-Volterra model.We think that this mutually beneficial relation of physics and adjacent disciplines is what makes the field of complex systems particularly attractive.
Next, we compute the value of χ.Since only f 0 and f 2 are present in Eq. ( 21), and only f 2 is divergent, we group in terms proportional to f 2 to obtain 0 which we can check always has χ = −1/u as its negative solution.Finally, from Eq. ( 22) we obtain M = 1/u and Var[x] = 0.
As expected, these values are independent of c and Γ, since in the limit of absent interactions Eq. ( 4) becomes a set of independent logistic maps.In this case, all species survive with abundance 1/u, which is what we obtain.

Limit α → αc
There are two different scenarios for the limit α → α c (where α c is the location of the phase transition).

C. Disorder average
We follow mainly [3,4].The only term in Eq. (S5) containing the disorder variables {c ij } and {ξ To simplify the calculations we will only keeps the terms that are leading order in N , since we eventually intend to take the thermodynamic limit N → ∞.To estimate the order of each term we use their variance, thus and implying that We first perform the average of (S16) over the {c ij }.We begin by Taylor expanding the exponential where we have used J ij = J ji , which follows directly from Eq. ( 3), and we have introduced the abbreviation b ij = dt xi x j .The different factors in the product are independent (since i < j) so we can average over {c ij } directly to obtain where we have re-exponentiated the expansion with the understanding that the result is only valid in the thermodynamic limit, and used the notation Now we calculate the average with respect to the variables {ξ µ i }.To do this we rewrite the exponentials in Eq. (S21) as were we have set J 2 ij = ⟨J 2 ij ⟩ in the second term (keeping in mind that we are only interested in leading-order terms in the thermodynamic limit).
Again, retaining only leading-order terms, we can express the different combinations of b ij using the order parameters and find  The sum over µ in the penultimate line produces αcN identical factors, as indicated in the last line.At this point, we have 2c dt dt ′ (v 2 L(t,t ′ )Q(t,t ′ )+w 2 K(t,t ′ )K(t ′ ,t)) × e −iαcN dt K(t,t) exp i 1 N dt = D[z, w, ẑ, ŵ]e i dt [ẑ(t)z(t)+ ŵ(t)w−z(t)w(t)]− 1 2N i ( dt [ẑ(t)x i (t)− ŵ(t)x i (t)]) 2 +O(N −2 ) = D[z, w, ẑ, ŵ]e i dt [ẑ(t)z(t)+ ŵ(t)w−z(t)w(t)]− where in the last steps we have used the Taylor expansion log[cos(x)] = −x 2 /2 + O(x 4 ), and expanded the resulting square to introduce the order parameters.This expression can be further simplified by carrying out the integration in the z, w variables as follows = dx e ixx f (x, ŷ)δ(ŷ − x) = f (x, ŷ)e ixŷ . (S29) Applying the identity in Eq. (S29) to Eq. (S28) we obtain  Notice that the perturbation field h(t) and the dynamical noise η(t) appear in the same way in Eq. (S53).We can simplify the equations by removing h(t) and writing the response function as

S3. FIXED-POINT ANALYSIS
A. Fixed-point ansatz In addition to focusing on fixed points we make three main assumptions: time translation invariance (TTI), finite integrated response (FIR), and weak long-term memory (WLTM).We will follow mainly section 4.5 of [6] to derive the consequences of these assumptions.

(S56)
As a consequence of TTI, we also have G(t, t ′ ) = G(t − t ′ ) = G(τ ), and by causality this response function can only be non-zero for t > t ′ .

FIG. 1 .
FIG. 1.Fixed points and divergences of the dynamics in Eq. (1).Each panel illustrates the behaviour of the model for different choices of c and Γ.The heatmap indicates the fraction of samples that converge to a fixed point after numerical integration of the gLVE.The criteria for the identification of convergence or divergence are described in Appendix A. The dashed lines are theoretical predictions for the onset of divergence (see Sec. III C).

FIG. 2 .
FIG.2.Test of analytical predictions for the order parameters against numerical simulations.The figure shows the fraction of surviving species ϕ, the mean abundance M , and the variance of abundances (q − M 2 ) as a function of the model parameter α (where P = αcN is the number of traits each species in the original system possesses).Lines are from the theory, derived in Eqs.(21) and(22), markers from numerical integration of the gLVE (N = 1000, tmax = 30, averaged over 10 realisations of the disorder).The remaining model parameters are c = 0.5, and Γ = 0.3.Vertical dashed lines indicate the onset of divergence as determined from the theory in Sec.III C.

FIG. 3 .
FIG. 3. Phase behaviour of the fully connected model (c = 1).Panel (a): Phase diagram for the model with c = 1, the only model parameters are then u and α.The system is stable to the right of the lines.At the dot-dashed line (u < 1/2) q, M and χ all diverge, and at the dashed line M and q diverge, but χ remains finite.Panel (b): Illustration of the behaviour of the abundances of individual species in the two different phases (convergence to a fixed point shown in green and red, diverging abundances in orange and blue).

9 FIG. 4 .
FIG. 4. Phase diagram for different choices of the connectivity c, and the symmetry parameter Γ.The coloured lines in each panel indicate where the linear instability occurs.The instability coincides with the divergence of M and q.The system is stable to the right of the line, abundances diverge on the left.

FIG. 5 .
FIG.5.Finite-time divergence of abundances.The heatmaps indicate the time, t div , at which abundances diverge, for initial conditions xi(0) = 1.Data is obtained from numerical integration of the gLVE.The dotted line is the phase boundary predicted by the theory.To the right of the phase boundary the system is in the stable phase, so that no divergence occurs.

20 FIG. 8 .
FIG. 8. Eigenvalue spectrum of the reduced interaction matrix for u = 5, c = 0.5, Γ = c/(c − 1), for different choices of the model parameter α.The vertical dashed lines indicate the real part of the right-most eigenvalue.

FIG. 9 .
FIG. 9. Eigenvalue spectrum of the reduced interaction matrix in the fully connected system.The matrices are symmetric, and their eigenvalues are therefore real-valued.Panel (a) is for u = 0.7, panel (b) for u = 0.3.

TABLE I .
Types of phase transition in the gLVE model with Hopfield-like interactions.The table summarises the different transitions, giving details about the nature of the divergence at the transition, and the associated behaviour of the spectrum of the reduced interaction matrix.