Imane Akjouj

Matthieu Barbier

Maxime Clenet

Walid Hachem

Mylène Maïda

François Massol

Jamal Najim
email: jamal.najim@univ-eiffel.fr

Viet Chi Tran

Complex systems in Ecology: a guided tour with large Lotka-Volterra models and random matrices

Keywords: 92D40, 92D25, 05C90, 34C35, 34D05, 60J80, 60F17 random matrices, dynamical systems, faisability, stability, structured networks, population dynamics

. We will discuss results and questions arising when we allow a combination of deterministic structure and randomness in Lotka-Volterra systems, see Sections 4 and 6(a).

Ecosystems represent archetypal complex dynamical systems, often modelled by coupled differential equations of the form

dx i dt = x i φ i (x 1 , • • • , x N ) ,
where N represents the number of species and x i , the abundance of species i. Among these families of coupled differential equations, Lotka-Volterra (LV) equations, corresponding to

φ i (x 1 , • • • , x N ) = r i -x i + (Γ x) i ,
play a privileged role, as the LV model represents an acceptable trade-off between complexity and tractability. Here, r i is the intrinsic growth of species i and Γ stands for the interaction matrix: Γ ij represents the effect of species j over species i. For large N , estimating matrix Γ is often an overwhelming task and an alternative is to draw Γ at random, parameterizing its statistical distribution by a limited number of model features. Dealing with large random matrices, we naturally rely on Random Matrix Theory (RMT).

Introduction: Complex networks, randomness and large dimension (a) Goals of this study

Ecosystems can be seen as archetypal complex dynamical systems, as they usually consist of a large number of interacting components with heterogeneous properties. In the present article, these components are species or sub-populations that evolve according to their own demographic dynamics and interact through various mechanisms (such as competition, predation or facilitation). These dynamics generally take the form:

dx i dt = x i φ i (x 1 , • • • , x N ) (1.1)
where x i (t) represents the (dimensionless) abundance or density of species i population, the function φ i encapsulates all the sources of growth and mortality, and N is the (large) number of interacting species. We focus in this article on the Lotka-Volterra (LV) model [START_REF] Lotka | Elements of Physical Biology[END_REF][START_REF] Volterra | Variations and fluctuations of the number of individuals in animal species living together[END_REF] which writes:

φ i (x) := φ LV i (x) = r i -x i + (Γ x) i , (1.2)
where r i represents the intrinsic growth rate of species i and Γ = (Γ ij ) is a N × N matrix with Γ ij representing the effect of species j on species i, see Section 2 for a detailed introduction of the LV model. Historically, LV model was designed to understand the interaction between two species (N = 2), which abundances often show an oscillating behavior (think of the celebrated Hare-Lynx example [START_REF] Odum | Fundamentals of ecology[END_REF]). In this study, we will focus on large foodwebs (N 1), equilibria and their stability. Oscillating behaviors in large dimension are beyond the scope of this survey but are discussed in the "open problems" section (Section 6(b)). In the Supplementary Material File (referred to as Supp-Mat file in the sequel) [4, Section 2], we provide a brief presentation of the N = 2 case.

Notice that φ i , often referred to as the net growth rate:

φ i = 1 x i dx i dt
represents the growth rate of abundance x i , i.e. its per capita rate of abundance change. Equation (1.1) incorporates two essential properties of biological dynamics: first, a species can be extinct (equilibrium at x i = 0), and second, a small population x i → 0 displays an exponential growth or decay with rate φ i | xi=0 referred to as the invasion growth rate (e.g. [START_REF] Marrow | The coevolution of predator-prey interactions : ESSS and Red Queen dynamics[END_REF][START_REF] Champagnat | Unifying evolutionary dynamics: from individual stochastic processes to macroscopic models via timescale separation[END_REF], see also the Supp-Mat file [4, part 4

(d)]).

To analyze and model such complex dynamical systems, a set of formal tools, including network theory, dynamical systems and Random Matrix Theory (RMT), have proven successful across a variety of scientific disciplines. This review is meant to offer a constructive viewpoint on the connection between these mathematical tools, especially RMT, and ecological systems, addressed to an interdisciplinary audience spanning theoretical ecology, complex systems, statistical physics and mathematical biology. We propose to focus on formal problems inspired by complex ecosystems, with the ultimate aim of answering ecological questions regarding the conditions of species coexistence, community diversity and ecosystem stability.

Through the prism of the Lotka-Volterra model with random interactions, our aim here is to give a diverse overview of concepts and questions that have proven fruitful in that line of work. In particular, we wish to point out results that we conjecture may also hold qualitatively, or even quantitatively, beyond this particular model, in the hope that such "universal" 1 behaviors might be shown in the future to capture some aspects of real ecosystems.

The connection of this model to empirical data is an important question which is beyond the scope of this review. We discuss this aspect in the discussion section (Section 6).

(b) Historical context

The first eminent proponent of applying results from complex dynamical systems to ecosystems was Robert May [START_REF] May | Will a large complex system be stable?[END_REF]. Inspired by the success of RMT for modelling unknown interactions in complex physical systems, such as large atomic nuclei [START_REF] Wigner | Random matrices in physics[END_REF], he argued as follows:

Obs Empirical observation: ecosystems with a large number of species appear to exist and persist for long times. Hyp 1 If we assume that species abundances are poised at some dynamical equilibrium

x * = (x * 1 , • • • , x * N ) with all x * i > 0 dx i dt = x * i φ i (x * ) = 0 (1.3)
Hyp 2 If we assume that their interactions are sufficiently complex to be modelled as random, i.e. more precisely, that the Jacobian matrix governing the linearized dynamics around the equilibrium

J ij = ∂(x i φ i ) ∂x j (x * ) (1.4)
is modelled as

J ii = -1 + M ii C ii J ij = M ij C ij (i = j) i.e. J = -I + C • M
where the random variables M ij are independent identically distributed (i.i.d.) centered random variables of variance V , C ij are i.i.d. with Bernoulli distribution with parameter C, called the connectance (such that a species has N × C links on average), and • denote the Hadamard product of matrices (pointwise multiplication of entries),

Res Then the equilibrium is stable for N CV < 1 and unstable for N CV > 1.

A formal presentation of May's result, in particular at a probabilistic level and emphasising the role of the large dimension of N , is provided in Supp-Mat file [START_REF] Akjouj | Supplementary Material / Complex systems in Ecology: a guided tour with large Lotka-Volterra models and random matrices[END_REF], (Section 1, Prop. 1.2).

May presented his conclusion as a paradox: earlier ecologists had imagined the observed complexity of natural ecosystems as a positive feature favoring their persistence (i.e. the more connected and/or strongly interacting species were, the more likely they would coexist), which could be contradicted by the loss of stability predicted in May's simplistic model. May's result, however, did not spell out where in the above argument lies the "paradox", i.e. which part of the argument should be reevaluated in the light of empirical evidence, or instead, should lead us to reevaluate earlier assumptions about ecosystems.

This result had a lasting influence on the theoretical ecology literature, yet it never received clear support from ecological data [START_REF] Tang | Correlation between interaction strengths drives stability in large ecological networks[END_REF][START_REF] Jacquet | No complexity-stability relationship in empirical ecosystems[END_REF][START_REF] Neutel | Reconciling complexity with stability in naturally assembling food webs[END_REF][START_REF] James | Constructing Random Matrices to Represent Real Ecosystems[END_REF]. After that, theoretical and empirical work on this topic has been diverse in scope and focus, but can be interpreted as questioning each of the four points summarized above 2 . In response to Obs, ecosystems might be less diverse and complex than they appear [START_REF] Hubbell | The unified neutral theory of biodiversity and biogeography[END_REF] (e.g. few species may coexist at any location, or species may exist in large number N but with low interaction variance V ). Contra Hyp 1, real ecosystem dynamics might be far from Robert [START_REF] May | Will a large complex system be stable?[END_REF] used a mathematical result to anchor a broader qualitative intuition: assuming a system with random features at dynamical equilibrium, there is an axis of effective diversity (equal to the product N CV in his model), combining system size (species number) and statistics (notably variance) of interaction, along which a qualitative transition from stability to instability is expected. Yet how likely is such an equilibrium to arise in the first place, and are there other qualitative transitions along this axis (or other axes that are as broadly relevant)?

equilibrium, e.g. dominated by transient excursions [START_REF] Hastings | Transient phenomena in ecology[END_REF]. Contra Hyp 2, ecological interactions might be structured so that the resulting dynamics deviate importantly from predictions from full randomness [START_REF] Allesina | Stability criteria for complex ecosystems[END_REF] (suggesting a greater role of non-random structure in ecology). Finally, even admitting Res, the loss of stability when N CV > 1 does not imply the extinction of species: the ecosystem may still persist in a steady out-of-equilibrium state [START_REF] Roy | Complex interactions can create persistent fluctuations in high-diversity ecosystems[END_REF], and thus this mathematical statement may not be interpreted ecologically as a limit on a realistic number of species.

The discussion that May's result generated [START_REF] Lehman | Biodiversity, stability, and productivity in competitive communities[END_REF][START_REF] Kokkoris | Variability in interaction strength and implications for biodiversity[END_REF] also highlights the issues arising from misunderstanding mathematical notions of equilibrium existence, equilibrium stability, and system persistence, as well as the gap between measurable quantities in ecology and those relevant to the analysis of ODEs [START_REF] Arnoldi | Unifying dynamical and structural stability of equilibria[END_REF][START_REF] Loreau | Biodiversity and ecosystem stability: a synthesis of underlying mechanisms[END_REF][START_REF] Wang | Ecosystem stability in space: α, β and γ variability[END_REF]. However, the "May-Wigner transition" was an influential result because, while it was proven for a particular case (a random linearized dynamical system), there were many reasons -both mathematical and empirical, from successes of RMT in other fields -to conjecture that it could have far broader generalizations.

(c) Qualitative questions

May's original work did not actually state what happens in arbitrary nonlinear ecological dynamics, besides the fact that it may exhibit a transition to instability (cf. Fig. 1). For instance, can some species go extinct without a loss of stability? [START_REF] Gardner | Connectance of large dynamic (cybernetic) systems: critical values for stability[END_REF][START_REF] Roberts | The stability of a feasible random ecosystem[END_REF][START_REF] Gilpin | Stability of feasible predator-prey systems[END_REF].

One aim of this review is to show that such questions, which cannot be addressed in May's setting, can be and have been fruitfully studied in subsequent decades using the LV model (2.1). A model that is simple, can be derived as an approximation of individual-based models (see Supp-Mat file [4, Section 3]), and is rich enough to demonstrate various dynamical phenomena.

Another important line of inquiry that has grown in response to May's work on stability has focused on imposed structure, searching for various deviations from full randomness that might possibly allow to restore feasibility, stability, or both. These structural features can be motivated by empirical networks, such as reproducing in broad strokes the hierarchical or group structure of a foodweb [START_REF] Allesina | The stability-complexity relationship at age 40: a random matrix perspective[END_REF], the size-dependence of trophic interactions [START_REF] Otto | Allometric degree distributions facilitate food-web stability[END_REF], or accounting for the spatial structure of ecosystems [START_REF] Gravel | Stability and complexity in model meta-ecosystems[END_REF][START_REF] Mougi | Food-web complexity, meta-community complexity and community stability[END_REF]. They can also arise from theoretical mechanisms, such as

(d) Outline of the article

In Section 2, we introduce the LV system of coupled differential equations. We discuss and motivate the large random interactions model, and introduce two such models, with independent interactions (i.i.d. model), and correlated reciprocal interactions (elliptic model). In Section 3, we present results for the two reference models (i.i.d. and elliptic) and address the question of the existence of a unique equilibrium x * = (x * i ) for a LV system with random interactions and its stability, that is the conditions for which x(t) → x * . We also focus on the feasibility of the equilibrium, that is the conditions for which no species vanishes at the equilibrium: x * i > 0 for all i. If feasibility is not reached, we describe the composition of the equilibrium in terms of surviving species x * i > 0 and vanishing species x * i = 0, together with various properties of the equilibrium. Beyond i.i.d. and elliptic models, there is a need to consider more realistic structures of interaction matrices which take into account important features of real foodwebs. In Section 4, we focus on sparse models and kernel matrices. We present a quick survey of the related literature in theoretical ecology and recall associated basic RMT results. In Section 5, we extend Lotka-Volterra models by adding extra randomness yielding to stochastic differential equations (SDEs).

Finally, we summarize in Section 6 what has been covered in this review and what has not. We also present some directions of investigation and open problems of interest.

We briefly describe the Supplementary Material File [START_REF] Akjouj | Supplementary Material / Complex systems in Ecology: a guided tour with large Lotka-Volterra models and random matrices[END_REF] which contents complete the main exposition. In [4, Section 1], we present May's model. In [4, Section 2], we present the deterministic LV model of size 2, for pedagogical purpose. [START_REF] Akjouj | Supplementary Material / Complex systems in Ecology: a guided tour with large Lotka-Volterra models and random matrices[END_REF]Section 3] establishes the connexion between individual-based models and LV model, while [START_REF] Akjouj | Supplementary Material / Complex systems in Ecology: a guided tour with large Lotka-Volterra models and random matrices[END_REF]Section 4] shows how individual-based models studied at different asymptotics yield community models with noise. In [4, Section 5], we provide mathematical details related to Section 4, in particular a certain type of deterministic networks and kernel matrices. More precisely, concerning the latter, we emphasize the connexion between general kernel matrices and Marchenko-Pastur distribution associated to large random covariance matrices.

The Lotka-Volterra model (a) The Lotka-Volterra system of differential equations

We are interested in many-species ecological dynamics and our main object of study will be the following LV system of differential equations:

dx i dt = x i (r i -x i + (Γ x) i ) , (2.1)
where

i ∈ [N ] := {1, • • • , N } and x = (x 1 , • • • , x N ).
The parameter N represents the number of species, supposed large, x i = x i (t) is a dimensionless quantity evolving in time t, in relation with the abundance of species i, r i represents the intrinsic growth of species i and Γ = (Γ ij ) is a N × N matrix reflecting the interaction effect of species j on the growth of species i:

(Γ x) i = j Γ ij x j .
Notice that if there are no interactions, i.e. Γ = 0, we recover for each species the standard logistic (Verhulst) differential equation which in particular reflects the limited amount of available resources.

royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 representing the net growth rate of species i. As mentioned above, we detail in [4, Section 3] how this LV model naturally emerges from an individual-based model when we consider a certain asymptotics of the birth, death and interaction rates. Let us also formalise in Remark 2.1 below how one can obtain the dimensionless, normalized system of equations (2.1), where each species has a comparable importance, from a similar system dealing with the actual abundances of species.

Remark 2.1 (Obtaining the dimensionless LV model (2.1)). In ecology, a natural formalization of the LV model involving the actual abundances or densities X i of the modelled species reads as:

dX i dt = X i r i -D i X i + (M X) i , (2.2)
where r i is species i intrinsic growth rate, D i is species i density-dependent term and M is the matrix of interaction coefficients. In order to obtain the dimensionless version of the Lotka-Volterra model presented in Equation (2.1), the following changes of variables are needed:

x i := D i X i , Γ ij := M ij /D j . (2.3)
Plugging these new variables into Eq. (2.2), one obtains Eq. (2.1).

In a large ecosystem consisting of N species (N 1), the precise knowledge of the interaction matrix Γ = (Γ ij ) among these species is often out of reach. An interesting alternative is to model the N × N matrix Γ with random entries and to rely on RMT. The statistical properties of the entries may then reflect a partial knowledge on the ecological interaction network.

(b) Two random interactions models

In this section, we precisely describe two models for the matrix Γ with random interactions. The first is the simplest theoretical baseline of i.i.d. entries with zero mean, while the second model is a natural extension that allows to represent more types of ecological interactions [START_REF] Tang | Reactivity and stability of large ecosystems[END_REF]. Other (more involved) models will be discussed in Section 4.

Independent and Identically Distributed entries.

In this model, the entries Γ ij are i.i.d., with a N -dependent common distribution and can be:

(i) Γ ij = A ij √ N or (ii) Γ ij = A ij α N √ N (2.4)
where the A ij 's are i.i.d., E A ij = 0, E A 2 ij = 1 and a distribution independent from N . In the case (i), the N -1/2 -normalization casts matrix Γ = A √ N into the framework of RMT, where the limiting properties of the spectrum and eigenvectors of matrix Γ are well described.

The circular law (cf. [START_REF] Bordenave | Around the circular law[END_REF]) asserts that the spectrum of Γ converges toward the uniform distribution on the disk of radius 1, see Fig. 2-(a). Moreover, denote by ρ(Γ ) Γ 's spectral radius:

ρ(Γ ) := max{|λ(Γ )| , λ(Γ ) ∈ C eigenvalue of Γ },
then its asymptotic behaviour is well-known: ρ(Γ )

a.s ----→ N →∞
1 , where a.s --→ stands for the almost sure convergence. As a consequence, matrix Γ has a non-vanishing macroscopic effect on the dynamical system (2.1) even for large N .

In the case (ii), there is an extra normalization term α N which may or may not depend on N . If it does not depend on N , then it is simply a parameter that allows to tune the variance of the entries since var(Γ ij ) = 1 α 2 N . If α N grows to infinity as N → ∞, it has the effect of asymptotically squeezing to zero the contribution of the interaction matrix Γ as

ρ(Γ ) = ρ(A/ √ N ) α N ----→ N →∞ 0 .
Despite this, we will see in Section 3(d) that α * N = 2 log(N ) is the threshold to reach feasibility (no vanishing species -see the formal definition in Section 3(a)-(ii)).

The elliptic model.

Two assumptions of the i.i.d. model are commonly relaxed to describe a wider range of ecological scenarios. First, while the i.i.d. model assumes that interactions have zero mean, ecological networks often contain interactions with a prescribed sign. Second, the i.i.d. model enforces the reciprocal interactions Γ ij and Γ ji to be uncorrelated. However, a large literature, for random as well as for deterministic interactions, deals with symmetric matrices Γ ij = Γ ji , which can for instance arise in ecology in the case of competition [START_REF] Macarthur | Species packing, and what competition minimizes[END_REF], or skew-symmetric matrices Γ ij = -Γ ji [START_REF] Redheffer | Volterra multipliers I[END_REF], which were originally proposed for predator-prey interactions [START_REF] Volterra | Variations and fluctuations of the number of individuals in animal species living together[END_REF]. The random symmetric case is well-known in RMT and is referred to as the Wigner model, see for instance [START_REF] Bai | Spectral analysis of large dimensional random matrices[END_REF]Chap. 2].

These various cases can be unified into the elliptic model [START_REF] Girko | Elliptic law[END_REF][START_REF] Sommers | Spectrum of large random asymmetric matrices[END_REF][START_REF] O'rourke | a Low rank perturbations of large elliptic random matrices[END_REF], which displays a richer statistical structure than the i.i.d. model. The entries of matrix Γ write

Γ ij = A ij √ N + µ N , (2.5)
where µ is a deterministic quantity, A ij is centered with variance equal to 1. In particular, µ N stands for Γ ij 's expectation. The second feature of the model is the existence of a correlation between the random variables A ij and A ji (i = j):

corr(A ij , A ji ) = ξ ∈ [-1, 1] (2.6)
while A kk and {A ij , A ji } are independent for k, i, j and i < j.

Under the elliptic model and in the case where µ = 0 and |ξ| = 1, the spectrum of Γ converges toward the uniform law on the ellipse see Figure 2, (b) and (c) (hence the name). If µ > 1 and does not belong to E ξ , then it has been shown in [START_REF] O'rourke | Low rank perturbations of large elliptic random matrices[END_REF] that we witness an extra outlier located near µ + ξ µ (i.e. a single random eigenvalue of Γ will converge to µ + ξ µ as N → ∞), see Figure 3. We may also consider this model with an extra normalization, as in (2.4)-(ii):

E ξ = z = x + iy, x 2 (1 + ξ) 2 + y 2 (1 -ξ) 2 ≤ 1 ,
Γ ij = A ij α N √ N + µ N . (2.7)
Among the key issues we will address is the study of the various regimes of the dynamical system (2.1) with random interactions: existence of an equilibrium, study of its stability, etc. Remark 2.2. Using random interaction models means that we hope to predict essential aspects of an ecosystem only on the basis of statistical features of the network, rather than its detailed structure. In May's framework (see Sec. 1(b)), the main parameter turned out to be the complexity N CV, combining the number N of species, the connectance C and the variance V of the interactions. The system's equilibrium is predicted to go from stable to unstable as its complexity parameter N CV crosses a threshold.

Our LV models are parametrized by the variance 1 α 2 of the interactions, their bias µ and their correlations ξ etc. One of our goal in the sequel will be to investigate how these parameters can be combined to characterize the behavior of the corresponding dynamical system.

Equilibrium, coexistence and stability

Due to the form of the dynamical system (2.1), standard ODE results yield that if

x(t = 0) = x 0 > 0 (resp. x(t = 0) = x 0 ≥ 0) componentwise, then x(t) > 0 (resp. x(t) ≥ 0) for all t > 0. We are thus interested in equilibrium points x * = (x * i ) i∈[N ] satisfying x * i φ LV i (x * ) = x * i (r i -x * i + (Γ x * ) i ) = 0 , i ∈ [N ] and x * i ≥ 0 .
Let us present the main properties of the equilibria that have been investigated in the literature.

Outline of the section

In Section 3(a), Paragraph (iii), we list the main questions we want to address. In Section 3(c), we provide general criteria to assess the existence of a unique stable equilibrium. In Section 3(d), we address the question of feasibility. In Section 3(e), we provide heuristics on the number of surviving species in the case where the equilibrium is not feasible. In Section 3(f), we present in Table 1 an overview of the main results and remaining open questions. The phase diagram in Fig. 9 may also help the reader to navigate among the various results.

(a) Terminology and main questions

Consider the general system (1.1). If needed, we will use φ LV i instead of φ i .

(i) Various notions of stability

The most common notion of stability in dynamical systems is the so-called Lyapounov stability:

The equilibrium x * is Lyapounov stable if for any neighborhood U of x * , there exists a neighborhood W of x * such that

x(0) ∈ W =⇒ x(t) ∈ U for all t ≥ 0 .
From now on, we will simply refer to stability instead of Lyapounov stability.

The equilibrium x * is said to be asymptotically stable, if and only if it is stable and the neighborhood W can be chosen so that

x(0) ∈ W =⇒ x(t) ----→ t→∞ x * .
One can also get interested in global stability, in the sense that x * is stable and

∀ x(0) ∈ (0, ∞) N , x(t) ----→ t→∞ x * .
The following theorem provides a necessary condition for stability to hold (notice that the result below holds for general functions φ i ).

Theorem 3.1 (see Takeuchi [START_REF] Takeuchi | Global dynamical properties of Lotka-Volterra systems[END_REF], Theorem 3.2.5). Consider the system

dx i dt = x i φ i (x 1 , • • • , x N ), i ∈ [N ] , (3.1)
with all φ i continuous. If an equilibrium point x * ≥ 0 of (3.1) is stable, then φ i (x * 1 , • • • , x * N ) ≤ 0, i ∈ [N ].
We will explain below how to relate Γ 's properties and the stability of the equilibrium.

(ii) Invasion, extinction, feasability and uninvadability

An important feature of Lotka-Volterra dynamics (2.1) is the fact that species can be extinct: x i = 0 entails dx i /dt = 0 and is always a possible equilibrium value for species i. On the other hand, if at some time x i (t) > 0, then positivity is maintained at all later times [START_REF] Hofbauer | The Theory of Evolution and Dynamical Systems: Mathematical Aspects of Selection[END_REF]. However, the solution can tend to zero asymptotically (see also Section 5 for finite populations). The fact that, for all species i, x i = 0 is always associated to some possible equilibria has motivated the notion of species invasion where one takes a system where x i (0) = 0, and perturbs it at time t by setting x i to a nonzero value. In particular, we often ask whether a species can invade from rare, i.e. what is the asymptotic behavior when x i (t) is set to a small positive value? This can be answered by considering the net growth rate φ LV i , defined for species i in Eq. (1.2) and representing its per capita rate of abundance change. The constant

φ LV i | xi=0 = φ LV i (x 1 , • • • , x i-1 , 0, x i+1 , • • • , x N ) = r i - j =i Γ ij x j (3.2)
defines the rate of exponential growth or decay of a small population x i ≈ 0 in an environment where the other species start with sizes given by x. This quantity is called the invasion growth rate in the community ecology literature, or invasion fitness (or just fitness) in evolution. The net growth rate φ LV i also has a probabilistic individual-based interpretation since it is related to the survival probability of a birth and death process that approximates locally x i (t) when x i is small and considering that the other sizes also do not vary much (see [4, Section 4(c)]). The species is said to be deterministically (or asymptotically) extinct if x i → 0, and permanent otherwise [START_REF] Hofbauer | Evolutionary games and population dynamics[END_REF][START_REF] Law | Permanence and the assembly of ecological communities[END_REF][START_REF] Metz | Adaptative Dynamics, a Geometrical Study of the Consequences of Nearly Faithful Reproduction[END_REF].

Ecologists have investigated invasion sequences, e.g. whether different dynamics and attractors are reached depending on the order in which species are introduced [START_REF] Robinson | Does invasion sequence affect community structure?[END_REF]. The same notions arise in an evolutionary context, where mutant types can be modelled as initially absent species. Champagnat et al. [START_REF] Champagnat | Polymorphic evolution sequence and evolutionary branching[END_REF] construct the polymorphic evolution sequence that alternates phases where the dynamics is described by the LV system (2.1) and phases of invasion of new arriving species, generalizing the trait substitution sequence process introduced in [START_REF] Metz | Adaptative Dynamics, a Geometrical Study of the Consequences of Nearly Faithful Reproduction[END_REF]. Merging coevolution and invasion sequences in simulation studies has also been a long-standing endeavour among ecologists [START_REF] Rummel | A theory of faunal buildup for competition communities[END_REF][START_REF] Taper | Models of character displacement and the theoretical robustness of taxon cycles[END_REF][START_REF] Calcagno | Diversity spurs diversification in ecological communities[END_REF][START_REF] Romanuk | Simulated evolution assembles more realistic food webs with more functionally similar species than invasion[END_REF][START_REF] Case | Invasion resistance arises in strongly interacting species-rich model competition communities[END_REF].

Beyond deterministic extinctions and invasions, we explain in Sec. 5 and [4, Section 3] how the expression for fitness arises in a probabilistic individual-based description and connects to survival probability in such a context.

The notion of feasibility will refer to an equilibrium x * = (x * i ) of (1.1) where all species coexist:

x * i > 0 for i ∈ [N ] .
Such an equilibrium will be called feasible.

The notion of uninvadability will refer to an equilibrium x * = (x * i ) of (1.1) satisfying:

x * i φ i (x * ) = 0 for i ∈ [N ]
where either (a) φ i (x * ) = 0 and x * i > 0, in which case the species x * i is said to survive or (b) x * i = 0 and φ i (x * ) ≤ 0, in which case the species is said to vanish (or to be extinct).

Such an equilibrium will be called uninvadable. Notice that it is stated in Theorem 3.1 that an equilibrium, to be stable, must be uninvadable.

(iii) Main questions we want to address

This raises a number of fundamental questions, in particular:

(i) Given parameter r and a random model for matrix Γ , is there a unique equilibrium? Is it locally stable? globally stable ? What is the proportion of surviving species and what are their statistical features? (ii) Given parameter r and a random model for matrix Γ , what are the conditions to get a feasible equilibrium? Is it unique? stable? (iii) Given an interaction matrix Γ , can we characterize the domain of growth rates r that allow feasibility?

We focus hereafter on questions (i) and (ii). We will not develop answers to question (iii), which assumes that r is likely to vary with environmental conditions while Γ is biologically fixed, since the main object of this study is LV systems when N is large and the interaction matrix Γ contains a random component. However, extensive works on question 3 can be found elsewhere, in the context of structural stability of LV equilibria [START_REF] Bastolla | Biodiversity in model ecosystems, I: coexistence conditions for competing species[END_REF][START_REF] Saavedra | A structural approach for understanding multispecies coexistence[END_REF].

There are several levels of answers to these questions, depending of the expected level of rigor. As we will see, we can provide fairly precise mathematical answers for a narrow range of assumptions (tight assumptions on the interactions Γ ). Theoretical physics tools and computer experiments will substantially relax these assumptions and widen our understanding.

(b) Linear Complementarity Problem: an important concept to study equilibrium dynamics

In this section, we provide the definition of the Linear Complementarity Problem (LCP), which is part of the theory of mathematical programming (see [START_REF] Murty | Linear Complementarity, Linear and Nonlinear Programming[END_REF][START_REF] Cottle | The linear complementarity problem[END_REF] for standard references). LCP has already been used in ecological contexts in [START_REF] Takeuchi | Global dynamical properties of Lotka-Volterra systems[END_REF]Chapter 3]. Given a N × N matrix M and a N × 1 vector q, we say that the LCP (M, q) admits a solution (z, w) ∈ R N × R N if there exist two such vectors satisfying the following set of constraints:

     w = M z + q ≥ 0 , z ≥ 0 , w T z = 0 .
In this case, we simply write z ∈ LCP (M, q) since w = M z + q can be inferred from z.

Consider the LV dynamics (2.1). An uninvadable equilibrium

x * = (x * i ) (if it exists) will satisfy      x * i φ LV i (x * ) = 0 , x * i ≥ 0 , φ LV i (x * ) ≤ 0 for all i ∈ [N ] . (3.3)
With the explicit form of φ LV i (x * ), see (1.2), this exactly means that x * ∈ LCP (I -Γ, -r).

(c) Criteria for existence and uniqueness of a globally stable equilibrium.

As mentioned above, an equilibrium point x * = (x * i ), if it exists, should satisfy

x * i (r i -x * i + (Γ x * ) i ) = 0, i ∈ [N ] . (3.4)
Hence, either x * i = 0, and the species i vanishes at equilibrium or r i -x * i + (Γ x * ) i = 0. We a priori do not know beforehand which species vanish and which ones remain. Moreover, uniqueness of the equilibrium may not be guaranteed.

A systematic way to find all the solutions is to arbitrarily fix the vanishing species x * i = 0 for i ∈ A, A being any subset of [N ] (2 N possibilities), then to solve the remaining set of equations

r j -x * j + (Γ x * ) j = 0, j ∈ A c := [N ] \ A . (3.5)
If the obtained x * j are positive, then (0, i ∈ A; x * j , j ∈ A c ) is a possible solution. Adding an uninvadability condition may considerably reduce the number of solutions, which can then be analyzed in the LCP framework.

The following result due to Takeuchi and Adachi [START_REF] Takeuchi | Global dynamical properties of Lotka-Volterra systems[END_REF] provides a sufficient condition for the existence of a unique equilibrium and the global stability of the LV system. It is the cornerstone to establish single equilibrium/stability conditions for interaction matrices with random entries and is based on the explicit construction of a Lyapunov function. Theorem 3.2 (Takeuchi and Adachi 1980, see [START_REF] Takeuchi | Global dynamical properties of Lotka-Volterra systems[END_REF] Th. 3.2.1). Consider the system (2.1):

dx i dt = x i (r i -x i + (Γ x) i ) = x i (r i + [(-I + Γ )x] i )
and assume that there exists a diagonal matrix Υ with positive diagonal elements such that matrix Υ (I -Γ ) + (I -Γ ) * Υ is positive definite. Then there exists a unique equilibrium point x * solution of the LCP (I -Γ, -r) and this equilibrium is globally stable in the sense that it is stable and

∀ x(0) ∈ (0, ∞) N , x(t) ----→ t→∞ x * .
In Takeuchi's language, -I + Γ is said to be Volterra-Lyapounov stable (VL-stable). As explained in [START_REF] Takeuchi | Global dynamical properties of Lotka-Volterra systems[END_REF], it is in general difficult to determine whether a given matrix is VL-stable. But one can compare it to other better known set of matrices : any negative diagonal dominant matrix, or negative semidefinite matrix is VL-stable and conversely, if A is VL-stable then A has all its principal minors positive. It is possible to be more precise when all off-diagonal elements have the same sign but not in general.

Combining Theorem 3.2 and standard results in RMT on the limit of the largest eigenvalue of a Wigner matrix, one can determine a set of parameters for which our random interaction models are VL-stable. We prove the following result:

(i) Matrix Γ is given by model (2.4)-(ii) and α > √ 2 , (ii) Matrix Γ is given by model (2.7) and the parameters (α, µ, ξ) ∈ R + × R × [-1, 1] satisfy α > 2(1 + ξ) and µ < 1 2 + 1 2 1 - 2(1 + ξ) α 2 .
Then almost surely, there exists N large enough such that the system (3.3) admits a unique solution x * = (x * i ), x * i ≥ 0. Moreover, this equilibrium is globally stable (see Section 3(a)-(i) for the definition).

We provide the proof under Assumption (i). For the general proof, see [57, Prop. 2.6].

Elements of proof. Assume (i) and let Υ in Th. 3.2 be given by Υ = I. Compute

I -Γ + I -Γ * = 2I -(Γ + Γ * ) = 2I - √ 2 α A ij + A ji √ 2N . (3.6) Notice that W = Aij +Aji √ 2N
is a symmetric matrix with independent centered entries on and above the diagonal with variance var

Aij +Aji √ 2 = 1 + δ ij
, where δ ij is the Kronecker symbol with value 1 if i = j and zero else. In RMT, W is referred to as a Wigner matrix, and its properties are well-studied. In particular, it is known that its largest eigenvalue λmax(W ) behaves as follows:

λmax(W ) a.s. ----→ N →∞ 2 ,
see for instance [START_REF] Bai | Spectral analysis of large dimensional random matrices[END_REF]Theorem 5.1]. Going back to (3.6), we get:

I -Γ + I -Γ * = 2I - √ 2 α W
and the smallest eigenvalue of this matrix is

λ min 2I - √ 2 α W = 2 - √ 2 α λmax(W ) a.s. ----→ N →∞ 2 1 - √ 2 α .
Taking α > √ 2 yields the desired result.

It is worth noticing that the (i) and (ii) respectively are sufficient conditions to get a unique globally stable equilibrium but simulations and heuristic arguments (see e.g. [START_REF] Bunin | Ecological communities with Lotka-Volterra dynamics[END_REF]) indicate that these conditions are not tight and one could observe a unique equilibrium for smaller α's.

Beyond Theorem 3.2, there exists other stability criteria of a similar flavour that may be useful. One can cite for example Goh [START_REF] Goh | Global stability in many-species systems[END_REF] that provides stability results in the case of the existence of a feasible equilibrium, or Champagnat et al. [START_REF] Champagnat | Convergence to equilibrium in competitive Lotka-Volterra and chemostat systems[END_REF] that provide a variant to Theorem 3.2. Those results could have similar RMT interpretations as Theorem 3.3.

(d) Unique feasible equilibrium.

Once the conditions are met so that the LV system of coupled equations admits a unique equilibrium, a natural question arises: is this equilibrium feasible, in the sense that x * > 0 componentwise? A negative answer has been brought by Dougoud et al. [START_REF] Dougoud | The feasibility of equilibria in large ecosystems: A primary but neglected concept in the complexity-stability debate[END_REF] in the case where α is fixed. As a first conclusion, getting feasibility requires α = α N ∞. This implies a qualitative change in the nature of the interactions since the random part of the interaction matrix would

ρ(Γ ) = O 1 α N → 0 (i.i.d. model) , ρ(Γ ) = | µ | + O 1 α N (elliptical model) ,
where O stands for the standard big O notation. In Bizeul and Najim [START_REF] Bizeul | Positive solutions for large random linear systems[END_REF], the feasibility threshold α N ∼ 2 log(N ) is established. We start with an argument of Dougoud et al. [START_REF] Dougoud | The feasibility of equilibria in large ecosystems: A primary but neglected concept in the complexity-stability debate[END_REF].

No feasibility if α is fixed.

Assume that the equilibrium point x * is feasible, then x * i > 0 for i ∈ [N ] and the equations (3.4) are equivalent to the linear system

(I -Γ )x * = r.
Let Γ be given by the i.i.d. model (2.4). It is well-known in RMT that the spectral radius ρ(A/ √ N ) almost surely converges to 1 (see e.g. [START_REF] Geman | The spectral radius of large random matrices[END_REF]). As a consequence, for every α > 1, ρ(Γ ) < 1 eventually (i.e. almost surely for large N ) and matrix (I -Γ ) is almost surely invertible for large N . Hence, the following algebraic representation of the equilibrium:

x * = I - A α √ N -1
r .

(3.7)

In the simpler case where r = 1 N , the N × 1 vector of ones, Geman and Hwang [START_REF] Geman | A chaos hypothesis for some large systems of random equations[END_REF] have proved that asymptotically, for every finite M ,

(x * 1 , • • • , x * M ) L ----→ N →∞ N M (1 M , σ 2 α I M ) ,
where

σ 2 α = 1 4α 2 -1 depends on α, N M (a, C) N M (a, C
) is the multivariate normal distribution with mean a and covariance matrix C, and L -→ stands for the convergence in distribution. As a consequence, Dougoud et al. [START_REF] Dougoud | The feasibility of equilibria in large ecosystems: A primary but neglected concept in the complexity-stability debate[END_REF] argued that under this interaction regime (fixed α) observing a feasible equilibrium was unlikely. In fact, the theoretical result by Geman and Hwang [START_REF] Geman | A chaos hypothesis for some large systems of random equations[END_REF] asserts that each component x * i of the equilibrium asymptotically behaves as an independent Gaussian random variable centered at 1, with a variance independent from N , hence the heuristics

P( min i∈[N ] x * i > 0) i∈[N ] P(x * i > 0) ----→ N →∞ 0 .
Otherwise stated, the initial assumption that x * > 0 is very unlikely to happen and is asymptotically a large deviation. This a priori analysis motivates the study of a feasible equilibrium under the regime

α = α N ----→ N →∞ ∞.
Feasibility when α N grows to infinity.

In the case where r = 1, there is a sharp phase transition around the threshold value α * N ∼ 2 log(N ) for both models (2.4)-(ii) and (2.7). Below the threshold, there is no feasibility with very high probability while above the threshold, feasibility occurs with probability growing to 1.

(i) If ∃ ε > 0 such that α N ≤ (1 -ε)α * N then P(min i∈[N ] x * i > 0) ----→ N →∞ 0 , (ii) If ∃ ε > 0 such that α N ≥ (1 + ε)α * N then P(min i∈[N ] x * i > 0) ----→ N →∞ 1 .
Although the full proofs of the theorem are involved, simple heuristics capture the phase transition in the i.i.d. case and is presented hereafter.

Remarks

• These results are established in [START_REF] Bizeul | Positive solutions for large random linear systems[END_REF] and [START_REF] Clenet | Equilibrium in a large Lotka-Volterra system with pairwise correlated interactions[END_REF] for the i.i.d. case and the elliptic case. • In the case where r = 1 is still positive componentwise, there is not a sharp threshold at α * N = 2 log(N ) but rather a transition buffer [α * min,N , α * max,N ] from non-feasibility (α N < α * min,N ) to feasibility (α N > α * max,N ). Details can be found in [START_REF] Bizeul | Positive solutions for large random linear systems[END_REF]Section 4.2].

Heuristics to understand the phase transition in Theorem 3.4. Assume that Γ is given by model (2.4)-(ii). In the representation (3.7) of the equilibrium, expand the inverse matrix as a Neumann series and only consider the first order expansion:

x * = I - A α N √ N -1 1 = 1 + A α N √ N 1 + • • • Every component x * k of x * writes x * k = 1 + Z k α N + • • • where Z k = j∈[N ] A kj √ N .
Notice that the Z k 's are i.i.d. N (0, 1). Pushing the approximation and taking the minimum yields

min k∈[N ] x * k 1 + min k∈[N ] Z k α N .
Now standard results from extreme value theory yield

min k∈[N ] Z k ∼ -2 log(N ), hence min k∈[N ] x * k 1 - 2 log(N ) α N .
The relative position of α N with respect to α * N = 2 log(N ) yields the desired result.

(e) Unique equilibrium with vanishing species

In the case of a unique equilibrium with species vanishing when t → +∞, it is interesting to understand some properties of the survivors such as the individual distribution of the abundance of a given species, the number of vanishing species, etc. Various techniques (such as the replica method from theoretical physics) yield quantitative heuristics validated by simulations. A full mathematical analysis remains currently out of reach.

(i) Number of vanishing species.

We mentioned earlier that should the parameter α (related to the strength of the interaction) be constant or less than 2 log(N ), the equilibrium x * will feature vanishing components x * i = 0 representing disappearing species.

In this section, we address the question of estimating the proportion of surviving species p = p(α) as a function of parameter α. In [START_REF] Bunin | Ecological communities with Lotka-Volterra dynamics[END_REF], Bunin provides heuristics based on the cavity method to address this question, while in [START_REF] Galla | Dynamically evolved community size and stability of random Lotka-Volterra ecosystems (a)[END_REF], Galla establishes equations comparable to those of Bunin based on dynamical generating functionals techniques. Both heuristics apply for the elliptical model (2.7). For the i.i.d. model, a simple order statistics argument can be found in [START_REF] Clénet | Equilibrium and surviving species in a large Lotka-Volterra system of differential equations[END_REF][START_REF] Clenet | Surviving species in a Large Lotka-Volterra system of differential equations[END_REF]. These a priori different methods yield the same equations from which one can extract p(α).

Given the random equilibrium x * , we introduce the following quantities:

S = {i ∈ [N ], x * i > 0} , p = |S| N , m2 = 1 |S| i∈[N ] (x * i ) 2 .
Denote by Z ∼ N (0, 1) and by Φ the cumulative Gaussian distribution function:

Φ(x) = x -∞ e -u
m √ p Φ -1 (1 -p) + α = 0 , (3.8) 1 + 2m √ p α E Z | Z > - α m √ p + m 2 p α 2 E Z 2 | Z > - α m √ p = m 2 (3.9)
admits a unique solution (p * , m * ) and the following convergence holds p a.s.

----→ n→∞ p * and m a.s.

----→ n→∞ m * .

Remark 3.6. Notice that the condition α > √ 2 guarantees by Theorem 3.3 that a.s. eventually there exists a unique equilibrium. This condition is sufficient but might not be necessary. In the simulations hereafter, we also test the case where α ∈ (1, √ 2] and observe that with high probability, there exists a unique equilibrium and a good matching with equations (3.8)-(3.9).

Simulations. We describe the setup of simulations to be compared with the theoretical values announced in Conjecture 3.5. Similar simulations can be found in [START_REF] Bunin | Ecological communities with Lotka-Volterra dynamics[END_REF].

We fix N = 1000 and draw L independent realizations of matrices A (i) . We then compute the corresponding equilibria x * (i) (α) and their related quantities (p (i) (α), m(i) (α)) for a given α > 0. We finally compare the empirical Monte Carlo averages:

pL (α) = 1 L L i=1 p(i) (α) and mL (α) = 1 L L i=1 m(i) (α)
to their theoretical counterparts p * (α), m * (α), solutions of (3.8) and (3.9). For α ∈ ( √ 2, 2 log(1000)), we consider L = 500 repeated samplings and for α ∈ (1, √ 2], L = 100. As shown in Figure 4, the matching is remarkable, even for α below √ 2. (ii) Single species distribution.

The previous conjecture provides an estimation of the proportion of surviving species p * (α). We go here one step further and describe the distribution of a given abundance x * i where index i corresponds to a surviving species. Conjecture 3.7. Let α ∈ √ 2, 2 log(N ) and let i ∈ S, i.e. i corresponds to a surviving species. Let p * , m * be the solutions of (3.8) and (3.9) and Z ∼ N (0, 1). Then the distribution of x * i is a truncated Gaussian:

L(x * i ) = L 1 + m * √ p * α Z Z > - α m * √ p * .
Otherwise stated, x * i admits the following density:

f * (v) = 1 (v>0) Φ(-δ) δ √ 2π exp - δ 2 (v -1) 2 2 where δ = α m * √ p *
and Φ stands for the cumulative Gaussian distribution.

The matching between the theoretical density f * given in Heuristics 3.7 and a histogram of a given equilibrium x * is illustrated in Figure 5. In particular, the theoretical distribution matches, even with non-Gaussian entries (see Fig. 5b).

(iii) Interactions between survivors

When only a fraction of species survive in the unique LCP equilibrium, one can also ask how the interactions restricted to the survivors are modified. Mathematically, this boils down to consider the submatrix (Γ ij ) i,j∈S . Of course, the lines and columns that are selected depend on the initial realisation of the matrix Γ and it is not an easy task to predict the new statistical features of the entries. Nevertheless, heuristics for these quantities have been given in [START_REF] Bunin | Interaction patterns and diversity in assembled ecological communities[END_REF], using the cavity method. These authors have obtained general formulas for the model (2.5) that can be found in [START_REF] Bunin | Interaction patterns and diversity in assembled ecological communities[END_REF], but for the sake of simplicity, we present it here in the case ξ = 0. surviving species, the entry Γ ij is still Gaussian but with the following bias and correlation:

E[Γ ij | x * ] - µ N = - x * i -1 + µx * N × (x * ) 2 x * j and Corr Γ ij , Γ ik | x * = - x * j x * k N × (x * ) 2 , ( 3.10)
where, for any vector

x = (x i , i ∈ [N ]), x 2 = (x 2 i , i ∈ [N ]) and x := 1 N N i=1 x i .
In [START_REF] Barbier | Fingerprints of high-dimensional coexistence in complex ecosystems[END_REF], a slightly different maximum likelihood viewpoint is adopted, where the elliptic model (2.5) is considered as a prior distribution and the corresponding posterior distribution is computed knowing the equilibrium, formulated as a linear constraint as in (3.5). Formulas (3.10) are recovered in this context. It is possible to study how these correlations impact spectral properties of the restricted matrix [START_REF] Barbier | Fingerprints of high-dimensional coexistence in complex ecosystems[END_REF][START_REF] Baron | Non-Gaussian random matrices determine the stability of Lotka-Volterra communities[END_REF]. Rigorous proofs remain out of reach.

(f) An overview of the results and some open questions

In the next table, we summarize some of the results presented here and list a few open questions.

α α fixed α = α N ∞ value √ 2

Structured models

In the previous section, we have studied at large the LV system (2.1) where the matrix Γ , supposed random, has either i.i.d. entries or follows the elliptic model, see Section 2(b). Another line of research focused on a model of the Jacobian of the system near equilibrium:

J = -I + M (4.1)
where M is random, the question being then to understand the relative localisation of the spectrum of M with respect to -1 to conclude on the stability of the system. This approach is historically the first with May's paper [START_REF] May | Will a large complex system be stable?[END_REF], a presentation of which is provided in [4, Section 1].

To progress toward a more realistic description of the reality, one is tempted to consider more involved models of random matrices to take into account more properties of the complex systems such as sparsity, existence of underlying structures, randomness beyond independence, etc.

For instance, the question of the effect of the structure of the ecological network on its feasibility and stability already appeared in [START_REF] Pimm | The structure of food webs[END_REF], where Pimm argued that connectance is not the only parameter that can influence the feasibility and stability of the networks, and started a theoretical study of structured (both deterministic and random) networks.

In this section, we present a variety of random matrix models beyond the i.i.d. and elliptic ones, highlighting their use in theoretical ecology and listing mathematical results and questions of interest. Often, mathematical results are not directly available and massive simulations remain the main approach to exploit the potentialities of such models. We start this collection of structured models with a quick reminder of what ecologists understand by "sparsity", since most of the model structures revolve around the idea of structuring where the zeros are in matrix Γ .

Sparsity in ecological networks. Empirically, in an ecosystem with N species, even if the maximal number of interactions is N 2 , the real number L of nonzero interactions is often much smaller. We define the connectance as

C = L N 2 . (4.2)
For a LV system, this means that matrix Γ has N 2 -L null entries. Although the interpretation is less obvious for random jacobian matrices (4.1), the notion of connectance already appeared in [START_REF] May | Will a large complex system be stable?[END_REF][START_REF] Gardner | Connectance of large dynamic (cybernetic) systems: critical values for stability[END_REF] and is important to model sparsity. Recently, Grilli et al. [START_REF] Grilli | Feasibility and coexistence of large ecological communities[END_REF] worked explicitly on the interaction matrix of a LV system and studied the stability and feasibility of the equilibrium as a function of various parameters, among which the connectance (see also [START_REF] Marcus | Local and collective transitions in sparsely-interacting ecological communities[END_REF][START_REF] Dunne | Food-Web Structure and Network Theory: The Role of Connectance and Size[END_REF]). Based on empirical evidence, [START_REF] Busiello | Explorability and the origin of network sparsity in living systems[END_REF] suggested that food webs can have very low connectance.

The notion of sparsity used thereafter is different from the one usually employed in mathematics, where matrices or networks are said to be sparse when the ratio L/N 2 goes to zero with N and dense when C is of the order of a positive constant for large N , the graph is said to be dense. In ecology, "sparse networks" are understood more broadly as networks in which not all interactions exist (i.e. matrices with zeros). In the following, we will consider sparse networks, in the ecological sense, with a wider range of connectance including the regime when C = O(1), which is called dense in mathematics.

Beyond the connectance, it is possible to take into account the structure of the network by setting some of the entries to zero, thus enforcing ecological sparsity, but in a structured manner. For this, we may use a matrix ∆ = (∆ ij ) where ∆ ij equals 1 if species j has an effect on species i and 0 otherwise. If one draws the system interactions as a graph, then ∆ can be interpreted as the adjacency matrix of this graph and the interaction matrix Γ or the community matrix M can then be represented as proportional to ∆ • A, with (∆ • A) ij = ∆ ij A ij and A is random either i.i.d. or elliptic. In such a model, ∆ represents the structure of the system summarizing ecological sparsity and A the (random) intensity of the interactions.

(a) The simplest model for sparsity: Erd ős-Rényi graphs

When all species play the same role in the foodweb and the only parameter of interest is the average number of interactions for a given species, it is natural to choose ∆ as the adjacency matrix of an Erdős-Rényi (ER) graph of size N : each coefficient of the random matrix ∆ has probability p to be nonzero, equal to 1, and probability 1 -p to be put to zero, independently of the others. The average number of edges in the graph is pN 2 , hence the connectance C equal to p.

ER in the mathematical literature

ER graphs are reference models in random graphs and their geometric properties are well-known (see e.g. [START_REF] Bollobás | Random graphs[END_REF][START_REF] Durrett | Random graph dynamics[END_REF][START_REF] Der Hofstad | Random Graphs and Complex Networks[END_REF]), together with the spectral properties of their adjacency matrices. In the regime when C = O(1), which is called dense by mathematicians, the ER matrix is a rank one deformation

√ N C ,
where • refers to the spectral norm when applied to a matrix. Notice that the precise understanding of the extreme eigenvalues of ∆ √ CN in sparse or very sparse regimes is still an active subject in RMT. A concise overview can be found in the introduction of [START_REF] Alt | Poisson statistics and localization at the spectral edge of sparse Erdős-Rényi graphs[END_REF].

ER in the ecological literature : sparsity increases stability

As developed in [4, Section 1], the case when M ij = ∆ ij A ij , with ∆ the adjacency matrix of a dense ER graph and A has i.i.d. centered entries with variance V has been already considered by May. This model is equivalent to the full model, where the entries have variance CV . In this case sparsity increases stability: the stability condition N V C < 1 is easily satisfied for small C.

The case when ∆ is the adjacency matrix of an ER graph but the model for the matrix A is more involved has been studied in particular in [START_REF] Allesina | Stability criteria for complex ecosystems[END_REF]. They use models for A that are of the same flavour as the elliptic model -for example, (A ij , A ji ) i<j both positive to model mutualistic systems or with opposite sign to model a prey-predator situation. As illustrated in Figure 6, in the mutualistic case, outliers with a large real part may strongly affect the stability. In [START_REF] Allesina | Stability criteria for complex ecosystems[END_REF], the authors also establish an explicit stability criterion adapted to each case, generalizing May's criterion and emphasize again that sparsity increases stability.

Γ = 1 α ∆ • A √ d = 1 α ∆ ij A ij √ d ij ,
where the A ij 's are i.i.d. and α is an extra normalization which may depend on N . Notice that the normalization is no longer √ N but √ d accounting for the fact that there are exactly d non-null entries per row. For such a model the connectance C equals:

C = d N
and the interest lies in "small" values of d.

(c) Introducing modularity through Stochastic Block Model (SBM)

Beyond the ER case, when every species equally interacts with another, it is often more realistic to consider communities within the ecosystem (also called modules), that is groups of species sharing the same connexion patterns. This leads to the celebrated Stochastic Block Model (SBM), introduced in [START_REF] Holland | Stochastic blockmodels: some first steps[END_REF], (see also [START_REF]Community detection and stochastic block models: recent development[END_REF][START_REF] Lee | A review of stochastic block models and extensions for graph clustering[END_REF] for reviews). Let r ∈ N be the number of communities. Given

• a vector of positive real numbers (π 1 , • • • , πr) such that r i=1 π i = 1, • an r × r matrix P, the corresponding SBM is a random graph whose vertices are partitioned into r communities C 1 , • • • , Cr, where each node belongs to the community C i , i ∈ {1, • • • , r}, with probability π i . Then, an edge between a vertex u ∈ C i and a vertex v ∈ C j exists with probability p ij , independently of all other edges.

SBM in the mathematical literature

There exists a huge literature on the SBM, initially introduced to analyze social networks, and extensively used in machine learning for modelling complex networks and address the community detection problem. The goal there is to design algorithms to cluster the different communities and estimate accurately matrix P , see for example [START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF][START_REF] Baskerville | Spatial Guilds in the Serengeti Food Web Revealed by a Bayesian Group Model[END_REF].

Again using the Hadamard product ∆ • A, the spectrum of the adjacency matrix ∆ associated to a SBM can be described, at least in simple cases. Consider for example a SBM with r = 2 communities of equal size (π 1 = π 2 = 1/2) and let p q q p , with p and q of order 0(1) (dense case). Then ∆ is a rank-two perturbation of a matrix with centered independent entries. Depending on the values of p+q 2 and p-q 2 , there can be up to two outliers in its spectrum. As in the ER case, sparse cases has also been recently considered [START_REF] Benaych-Georges | Spectral radii of sparse random matrices[END_REF].

SBM in the ecological literature : modularity increases stability

In the seventies, May and Pimm already considered rudimental forms of the SBM into the framework of the Jacobian model (4.1), to take into account some features of ecosystems such as modularity and compartmental models.

In [START_REF] May | Will a large complex system be stable?[END_REF], May presents a simple occurrence of SBM. He considers a SBM with r modules and a probability vector (c 1 , . . . , cr). This SBM corresponds to modules with no interactions, while within the ith block made of d i species, the interactions behave like an ER graph with connectance c i and variance V i . The stability condition reads :

max i∈[r] c i V i d i < 1,
hence modularity increases stability. This phenomenon is illustrated in Figure [START_REF] Wigner | Random matrices in physics[END_REF].

In [START_REF] Pimm | The structure of food webs[END_REF], Pimm addresses the following question "should model systems be organized into compartments of species characterized by strong interactions within compartments, but weak interactions among the compartments ?" A random version of such a model would correspond to a SBM with a matrix P having large diagonal coefficients and small off-diagonal ones.

(a) (b)

Figure 8: SBM with two communities of 500 species, N = 1000. The probability of interactions inside the first (resp. second) community is c 1 (resp. c 2 ) and the probability of interaction with species of the other community is equal to q. (a) Spectrum of the interaction matrix normalized by

1 √ (c1+q)N
, where ∆ is a SBM with c 1 = 0.5 > c 2 = 0.2 > q = 0.02 and A has i.i.d. Gaussian N (0, 1) entries. (b) Spectrum of the interaction matrix normalized by 1 √ qN where ∆ is a bipartite symmetric model, that is a SBM with q = 0.5 and c 1 = c 2 = 0 and A has i.i.d |N (0, 1)| entries.

More recently, the effect of modularity on the stability of the networks is extensively explored in [START_REF] Thébault | Stability of ecological communities and the architecture of mutualistic and trophic networks[END_REF] in the framework of a tamed version of LV equations. They evaluate modularity through an index, introduced by [START_REF] Barber | Modularity and community detection in bipartite networks[END_REF] and, through simulations, illustrate that persistence, that is the number of surviving species, increases with modularity in trophic networks (see also [START_REF] Stouffer | Compartmentalization increases food-web persistence[END_REF]) but decreases with modularity in mutualistic networks. It would be interesting to investigate whether mathematical results on SBM could help to understand their observations. The question of stability can be important for plant-pollinator ecosystems, corresponding to bipartite mutualistic networks (see e.g. [START_REF] Thébault | Stability of ecological communities and the architecture of mutualistic and trophic networks[END_REF][START_REF] Billiard | Continuous limits of large plant-pollinator random networks and some applications[END_REF]). In [START_REF] Billiard | Continuous limits of large plant-pollinator random networks and some applications[END_REF], the evolution of abundances can be approximated, when the number of species N tends to infinity, by a kinetic integrodifferential equation where the dense graphs are replaced by graphons. The theory of graphons is mathematically well developed but beyond the scope of this review (see for example [START_REF] Lovász | Large networks and graph limits[END_REF]). (d) Nested models : a few generalist and many specialist species.

In the Erdős-Rényi model or in SBM, the network is determined by considering the absence or presence of edges for each pair of vertices independently of the others. Other models of random graphs are defined by specifying the degree distributions. For example in the configuration models (also known as Molloy-Reed-Bollobás, see e.g. [START_REF] Bollobás | Random graphs[END_REF][START_REF] Der Hofstad | Random Graphs and Complex Networks[END_REF][START_REF] Molloy | A critical point for random graphs with a given degree sequence[END_REF]), independent random variables distributed with the target degree distribution are associated to each vertex and edges are formed by pairing at random the half-edges.

By choosing heavy-tailed degree distributions, one can create a few vertices with very high degrees (corresponding to generalist species) and a majority of vertices with low degree (corresponding to specialists). Such ecosystems are called nested. They have been modelled and studied, at least through simulations, see for instance [START_REF] Thébault | Stability of ecological communities and the architecture of mutualistic and trophic networks[END_REF] following [START_REF] Okuyama | Network structural properties mediate the stability of mutualistic communities[END_REF].

Nested ecosystems can also be described through random graphs with given expected degrees. This model is known as the Chung-Lu model : take a deterministic sequence w = (w 1 , . . . , wn), that will correspond to the expected degrees and draw an edge between vertex i and vertex j with probability w i w j / n i=1 w i independently of all other edges. If we choose all the weights to be equal to pn, we are obviously back to the ER model with connectance p but nested ecosystems can be modeled by choosing a power-law distribution for the weights, that is

w i = c -1 β-1 i
, for i greater or equal to some i 0 . In this case, the number of species interacting with k others is proportional to k -β . The spectrum of the adjacency matrix of such a graph has been studied in [START_REF] Chung | Spectra of Random Graphs with Given Expected Degrees[END_REF] where a phase transition is shown to occur at β = 2.5 : for β > 2.5, the largest eigenvalue behaves like √ m, the maximal degree in the graph, whereas for β < 2.5, the largest eigenvalue behaves like d, the weighted average of the square of the expected degrees. It would be interesting to investigate whether these mathematical results could be effectively used in the study of nested ecosystems.

(e) Kernel matrices

Definition of the model. Part of the literature on ecological networks considers that the interactions between two species depends on a distance between their respective values of some functional traits. The examples that we will present below, fit into the mathematical framework of kernel matrices. We have

Γ ij or M ij = f (g(x i , x j )), (4.3)
where x i is a vector modelling the traits of species i, g : R p × R p -→ R is a symmetric function, denoted as the kernel (corresponding to the measure of the distance), and f : R -→ R a function called the envelope. Examples are g(x, y) = x T y, or xy 2 and f (x) = exp(ex) or (1 + x) a etc.

Kernel matrices in the mathematical literature. Among these models, the first interesting and well studied case is the so-called Wishart case3 , when g(x, y) = x T y and f (x) = x. If the entries of the vectors are i.i.d. centered and normalized, then it is well known (see [START_REF] Marčenko | The spectrum of random matrices[END_REF]) that, if the ratio p/N of the number of traits over the number of species converges to τ, the spectrum converges almost surely to the Marcenko-Pastur distribution. We will develop the mathematical theory of kernel matrices in [START_REF] Akjouj | Supplementary Material / Complex systems in Ecology: a guided tour with large Lotka-Volterra models and random matrices[END_REF]Section 5(b)]. The main message is that in the RMT regime, that is when the number of traits is large and of the same order as the number of species, and if g and f are reasonable, they have almost no influence on the spectrum, meaning that, in the models, "any" kernel matrix could be without harm replaced by a Wishart matrix or a Gaussian kernel matrix.

Kernel matrices for ecological networks. Dieckmann and Doebeli consider a simple coevolutionary model [START_REF] Dieckmann | On the origin of species by sympatric speciation[END_REF] (see also [START_REF] Champagnat | Unifying evolutionary dynamics: from individual stochastic processes to macroscopic models via timescale separation[END_REF][START_REF] Roughgarden | Theory of population genetics and evolutionary ecology: an introduction[END_REF]) where interactions are of the form (4.3) with f (x) = ηc exp(-x 2 /2σ 2 c ) and g(x, y) = xy 2 , ηc ∈ R being the strength of the competition of mutualism and σc > 0 being the width of the ecological niches. The function f corresponds to a Gaussian kernel (also used in [START_REF] Meszéna | Competitive exclusion and limiting similarity: A unified theory[END_REF]). In [START_REF] Nuismer | Coevolution and the architecture of mutualistic networks[END_REF], they develop a phenotype matching model, where the interaction is stronger when the traits of two species are close but also phenotype difference (or threshold) model, in which successful interaction depends on the degree to which the trait of the second species surpasses the trait of the first one 4 (see also [100,101]). Other models involving thresholds can be found in [102] and in [103], the models involving a combination of characteristics of the species taken separately and a measure of the similarity between the traits.

Let us finally mention the work of [104], which lies in the LV framework, with a kernel matrix Γ for which the distance between the traits, is determined through a distance between species in their phylogenetic tree. The author addresses the questions of feasability and stability as we have detailed it for the elliptic case and explicitly uses the link with Wishart matrices mentioned above. It would be interesting to investigate how this point of view applies to other contexts.

Other community models with noise

The models presented in the previous sections can be complexified. For instance, in [105,106], a dynamical system allowing to incorporate homogeneous Gaussian vector field with gradient and solenoidal contributions (resuming to May's equation when linearizing around critical points) is presented. It is shown that in this case there can be a exponential number of equilibria, the majority being unstable. Also, in [107], interactions corresponding to locally tree-like networks are modeled, with a focus on transient early phases in Equations 2.1.

Stability questions as well as out of equilibrium studies have also been investigated in other kinds of stochastic Lotka-Volterra systems, when not only the interaction matrix is random but when also noise is added to the dynamics. We review different models with noise that are obtained from individual-based stochastic models by various rescalings. The evolution is described by stochastic differential equations (SDEs), of Feller or Ornstein-Uhlenbeck types, see [108] for a review. Other possible equations are presented in [4, Section 3] with time-scales separations and [109] for discussions in view of ecological systems. For such systems of stochastic differential equations (such as the ones in Prop. 5.2 and 5.3), the notion of stability for ODEs is replaced by the notion of equilibrium distribution -the only stationary state being often zero, corresponding to the extinction of all species. More precisely, the solution X(t) of a stochastic differential equation converges in distribution to the stationary probability measure µ on R N if for every continuous bounded function f , lim t→+∞ E(f (X(t)) = R N f (x)dµ(x) (see [110]). Long time behaviour of SDEs is the subject of an abundant literature in mathematics [111,112].

The LV eq. (2.1) can arise as limits of Individual Based Models (IBM), see [START_REF] Champagnat | Unifying evolutionary dynamics: from individual stochastic processes to macroscopic models via timescale separation[END_REF], when population sizes are large. We introduce a scaling parameter K referred to as a carrying capacity and assume that the N species have initially sizes of order K. The individuals can give birth to new individuals of the same species or die. More precisely, the natural birth and death rates of an individual of species i ∈

[N ] are b K i and d K i such that b K i -d K i = r K i .
The competition pressure (or extra death rate) exerted by an individual of species j on an individual of species i is Γ K ij , if the latter is non positive. In case Γ K ij is positive, it can be considered as an extra birth rate due to cooperation between the species i and j. Additionally an individual of species i experiences an intra-specific death rate proportional to the size of its species (each individual experiences an extra death term due to the logistic competition and equal to the size of the population i over K). The natural IBMs associated with (2.1) have two levels of stochasticity: (i) matrix Γ K is random, (ii) the occurrence of birth/death events is random. Hereafter, we work conditionally to Γ K .

Individual-based models have long been used for simulations in Ecology [113-116], (see [117,118] for software performing IBM simulations). By presenting the fluctuations arising in the convergence of the IBM abundances to LV abundances, we provide a link between the equations considered in this review and these algorithms (see [START_REF] Akjouj | Supplementary Material / Complex systems in Ecology: a guided tour with large Lotka-Volterra models and random matrices[END_REF]Part 3] for more mathematical details).

(t) = (Y K 1 (t), • • • , Y K N (t))
the vector of all the species' sizes, with forall i ∈ [N ], the following convergence in probability

lim K→+∞ Y K i (0) K = x i (0). (5.1)
The stochasticity of the birth and death events gives rise to an additional noise process compared to (2.1). More precisely, Y K (t) now satisfies the following SDE: for all i ∈ [N ],

dY K i (t) = r K i -Y K i (t) -(Γ K Y K (t)) i Y K i (t) dt + dM K i (t) , (5.2)
where

M K i is a martingale process, with E(M K i (t)) = 0, Cov(M K i (t), M K j (t)) = 0 and Var(M K i (t)) = 1 K E t 0 b K i + d K i + Y K i (s) + (Γ K Y K (s)) i Y K i (s) ds , (5.3)
(this is the analogue of Eq. (3.4) derived for N = 1 in [4, Part 3]).

We now detail two different limits that can be derived from the IBM depending on the chosen parameters and rescalings: (i) the Lotka-Volterra ODEs (2.1) with a fluctuation process of Ornstein-Uhlenbeck type when K → +∞ without time rescaling; (ii) the Feller-type diffusions when the birth-death dynamics is nearly critical and when time is also rescaled by K, that is when considering the process at times Kt, t ∈ [0, T ]. (a) Large population limit and fluctuation around the ODE (2.1) First, we let K → +∞ without rescaling time and with a fixed number N of species. We consider the rescaled process:

X K (t) := Y K (t) K .
Here, for all i, j ∈

[N ], b K i = b i , d K i = d i , Γ K ij = Γ ij K , (5.4)
where the quantities b i , d i and Γ ij do not depend on K and r i := b i -d i . Notice that the competition term Γ K ij can be understood as the extra death rate exerted by an individual of the species j on an individual of the species i. When the population is large and of order K, the interaction between the individuals of a given pair are weaker and therefore the competition term is rescaled in 1/K.

An averaging phenomenon appears (similarly to the law of large numbers): from (5.3), we can see heuristically that the noise disappear and the evolution equations (5.2) can be approximated by Lotka-Volterra ODEs (see e.g. [119, Theorem 2.1 p.456], or for more generalizations to measurevalued processes in [START_REF] Champagnat | Unifying evolutionary dynamics: from individual stochastic processes to macroscopic models via timescale separation[END_REF]108,[120][121][122] for rigorous proofs). Indeed, the variance of the martingale part in X K i is of order 1/K and converges to zero when K → +∞.

Proposition 5.1. Assume (5.1), (5.4) and assume that for all i ∈

[N ], sup K∈N E X K i (0) 3 < +∞, (5.5)
then, when K → +∞, the abundance processes (X K (t)) t≥0 converges uniformly on every compact time interval [0, T ] with T > 0 (as R N -valued processes) and in probability to population densities (x(t))

dx i dt = x i (r i -x i + (Γ x) i ) , i ∈ [N ]. (5.6)
Put formally, this means that

∀ε, T > 0 , P sup t≤T X K (t) -x(t) > ε ----→ K→∞ 0 .
We now consider the fluctuation process associated with this convergence:

η K (t) = √ K X K (t) -x(t) . (5.7) It is a R N -valued vector process whose ith coordinate is √ K(X K i (t) -x i (t)).
Another reformulation is that the stochastic process can be expressed as:

X K (t) = x(t) + η K (t) √ K .
Applying [119, Theorem 2.3, Chapter 11], we obtain that: Proposition 5.2. Under the assumption of Proposition 5.1 and assuming that, in distribution,

lim K→+∞ η K 0 = η 0 ∈ R N , (5.8)
then, when K → +∞, the process (η K (t)) t≥0 converges in distribution, and for the topology of uniform convergence on [0, T ] for T > 0, to the solution of the Ornstein-Uhlenbeck SDE:

dη(t) = r -2x(t) -Γ x(t) • η(t) dt + x(t) • Γ η(t) dt + diag (b + d + x(t) + Γ x(t)) • x(t) dB(t) , (5.9)
with initial condition η 0 defined in (5.8), B a standard N -dimensional Brownian motion, • the Hadamard product and diag(•) the diagonal matrix with diagonal entries the components of the vector in the argument. Equivalently, the componentwise definition of η = (η i ) is given for all i ∈ [N ] by:

dη i (t) = r i -2x i (t) -(Γ x) i (t) η i (t) dt + x i (t) Γ η(t) i dt + b i + d i + x i (t) + (Γ x) i (t) x i (t) dB i (t) .
Conditionally on matrix Γ , the solution x(t) of (2.1) is deterministic and hence the stochastic differential equation (5.9) is of Ornstein-Uhlenbeck type, i.e.

dη(t) = A(t)η(t)dt + Σ(t)dB(t)

where

A(t) = diag r -2x(t) -Γ x(t) + diag x(t) Γ, Σ(t) = diag b + d + x(t) + Γ x(t) • x(t) .
The solution is a centered Gaussian process with covariance function:

Cov η(t), η(s) = t∧s 0 e -v 0 A(u)du Σ 2 (v)e -v 0 A(u)du dv ,
a N × N matrix-valued function. This central limit theorem quantifies the convergence rate in Proposition 5.1 and allows for example to compute confidence intervals.

(b) Noisy versions of Lotka-Volterra equations (i) Limiting Feller diffusion for large population with time rescaling

Another way to exhibit SDEs is to consider diffusive time rescaling in the almost-critical case (e.g. [6, Section 4.2], when the growth rates are close to zero). More precisely, if i.e. Σ = diag(σ i ). The fact that the birth and death rates are of order K corresponds to accelerating the time proportionally to the factor K that also rescaling the population size. That the species i is close to criticality appears in the fact that both the birth and the death rate have the same leading term in Kσ i with the same coefficient σ i . This coefficient can depend on the species.

b K i = Kσ i + b i , d K i = Kσ i + d i , and Γ K ij = Γ ij K , ( 5
Proposition 5.3. Assume (5.1), (5.5) and the rates (5.10). Then, when K → +∞, the process (X K (t)) t≥0 converges uniformly on every compact time interval [0, T ] with T > 0 and in distribution to the solution of the following Feller stochastic differential equation:

dX(t) = (r -X(t) + Γ X(t)) • X(t) dt + 2ΣX(t) • dB t (5.11)
where B is a N -dimensional standard Brownian motion independent of Γ , and where the function x → √ x is applied elementwise to the vector ΣX(t).

The random noise appearing in (5.11) comes from the rapid successions of birth and death events in this accelerated time-scale. See where B(t) is a N -dimensional standard Brownian motion independent of Γ , and where the function f : R + → R + is applied elementwise to any vector x = (x i ), i.e. f (x) = (f (x i )).

In the framework of a symmetric interaction matrix, Biroli et al. [123] unveil, by the replica method, the large-N system behavior of X, recovering the parameter regions for which the system has a single equilibrium or multiple equilibria. Another generalization was done by Roy et al. in [124], where Γ can follow the general elliptic model. These authors study the large-N limit of the SDE (5.12) by relying on a dynamical mean field approach based on the dynamical cavity method, as detailed in [125].

A better understanding of these results from a mathematical perspective, as well as their generalization to more sophisticated models than the elliptical model, are interesting and useful research directions which have not been undertaken so far to the best of our knowledge. One of the first mathematical formalizations of this class of problems goes back to the work of Ben Arous and Guionnet [126], who were interested in the dynamics of mean-field spin glasses. In that setting, their analogue of our SDE (5.12) is a Langevin version of the so-called Sherrington-Kirkpatrick model for the spin glass dynamics. In the same line of thought, Faugeras et al. [127,128]

Discussion

Our guided tour through large LV models has highlighted the importance of some features of matrix Γ to understand the dynamics of complex ecological systems. Among these features, the first two statistical moments of the interactions, µ/N and σ 2 /(N α 2 ), are of paramount importance, as well as the scaling of the normalisation applied to matrix Γ (1/N for the deterministic part of the matrix, 1/ √ N or 1/α N √ N for its random part). Our review covers some key topics on the analysis of large LV models, i.e. conditions leading to a feasible equilibrium, to a unique maximal equilibrium and the link between the interaction matrix and equilibrium stability (Sec. 3). In Sec. 4 and [4, Section 5], we mention some sophistications (block, trait-based, sparse matrices) of the model that can make it more realistic. Theoretical results obtained on kernel matrices (many traits, many species) [START_REF] Akjouj | Supplementary Material / Complex systems in Ecology: a guided tour with large Lotka-Volterra models and random matrices[END_REF]Thm. 5.3] suggest that the study of LV models obtained from trait-based interaction matrices could be simplified by assuming that many interaction matrix spectra will resemble that of a Wishart matrix. However, this result does not yield any direct conclusion on the feasibility of an equilibrium of the LV model, nor does it help generalize the heuristics 3.5 to such structured models.

We have not considered models based on fully organized interaction matrices, which are not amenable to an RMT analysis but sometimes allow some direct analysis through classic analytical tools (e.g. Lyapunov functions, see also [129]). For instance, modelling food webs as strictly organized by trophic levels (i.e. species from trophic level k can only be positively affected by some species of level k -1 and negatively affected by some species of level k + 1), it is possible to express coexistence conditions as mathematical conditions on the covering of the food web by pairs of interacting species [130,131]. This finding, which highlights the difficulty of having a feasible equilibrium in level-organized food webs, hints at the potential importance of omnivory (i.e. that interactions are not strictly organized by levels, such that a predator can also feed on the prey of its prey species) to explain the stability of real food webs. An interesting endeavour could be to formalize omnivory in the context of random interaction matrices in order to tease apart the effect of omnivory from that of the food web being acyclic (which would also break the trophic level-based nature of the interaction matrix used by [130,131]).

(a) Empirical data and real interaction networks

There is a huge literature on statistical models for ecological networks, see [132][133][134][135][136][137][138][139]. The variety of data available and the specificity of each ecosystem explain this vast corpus, which in itself would deserve an entire review.

Many papers consider variations of the LV system, for instance integrating functional responses [140] (and references therein). The inference and empirical testing of models is limited by the type and volume of data that can be obtained through the observation of natural systems or controlled experiments. For most systems, we can only access abundances x i (often in a single snapshot that may not be at equilibrium, less frequently in time series). For some systems, especially with predator-prey or pollinator-plant interactions, we can also know which interactions Γ ij are zero or nonzero, as species may be limited in which partners they can interact with, but quantitative estimates of nonzero interactions remain out of reach. Therefore, exhaustive information on parameters r i and Γ ij is not accessible without additional assumption [START_REF] Jacquet | No complexity-stability relationship in empirical ecosystems[END_REF]: the main method for estimating them reliably is to observe the growth rates and/or equilibrium abundances of each species in isolation and in combinations, which is only possible in small-N experiments. Abundance distributions have been exploited to fit many other models, e.g neutral theory [START_REF] Hubbell | The unified neutral theory of biodiversity and biogeography[END_REF]. Additional data may help estimation and model validation [START_REF] Barbier | Fingerprints of high-dimensional coexistence in complex ecosystems[END_REF]141,142]. Finally, notice that although the models discussed above have a high number of parameters, their distributions have low degrees of freedom, see for instance the elliptic model.

It is important to keep in mind that real interaction networks such as food webs are likely not randomly structured. Interactions between species play a crucial role in their respective fitness, so that when given sufficient variation in traits governing interaction strength, coevolution of species will lead to particular network structures [START_REF] Allhoff | Evolutionary food web model based on body masses gives realistic networks with permanent species turnover[END_REF][START_REF] Låstad | The geometry of evolved community matrix spectra[END_REF]101,143]. In the sequential invasion simulations of food webs [START_REF] Låstad | The geometry of evolved community matrix spectra[END_REF], the spectrum of the Jacobian matrix of the underlying LV system does not converge towards an ellipsis as the number of invasion attempts increases, but rather becomes bimodal, thus highlighting a possible structuring effect of evolution sensu lato. In another

Figure 1 :

 1 Figure1: Schematic cartoon of early arguments from randomness in complex ecological systems. Robert[START_REF] May | Will a large complex system be stable?[END_REF] used a mathematical result to anchor a broader qualitative intuition: assuming a system with random features at dynamical equilibrium, there is an axis of effective diversity (equal to the product N CV in his model), combining system size (species number) and statistics (notably variance) of interaction, along which a qualitative transition from stability to instability is expected. Yet how likely is such an equilibrium to arise in the first place, and are there other qualitative transitions along this axis (or other axes that are as broadly relevant)?

5 Figure 2 :

 52 Figure 2: Spectrum of 500 × 500 centered random matrices with (a) N -1/2 -normalized i.i.d. standard Gaussian entries, (b) elliptic distribution ξ = -0.5, (c) elliptic distribution ξ = 0.5. The points represent the eigenvalues. The solid curve represents the boundary of the support of (a) the circular law (uniform law on the disc), (b) and (c) the elliptic distribution with parameter ξ.

  royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3 :

 3 Figure 3: Spectrum of a deformed 500 × 500 elliptic matrix with parameters ξ = 0.5, µ = 2. The solid line represents ellipse E ξ , the boundary of the support of the limiting spectral distribution for an elliptic model. Notice the extra outlier near 2 + ρ 2 as expected.

  royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Theorem 3.3 (Unique equilibrium and stability under uninvadability condition). Let α N = α be fixed. If one of the following conditions is satisfied:

  royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . have a (macroscopic) vanishing effect:

Theorem 3 . 4 .

 34 Let Γ be given by model (2.4)-(ii) or (2.7); in the latter case assume that µ < 1. Assume α N ----→ N →∞ ∞ and denote by α * N = 2 log(N ). Let r = 1. Then x * given by (3.7) is well-defined and

  royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 4 :

 4 Figure 4: The plot represents a comparison between the theoretical proportion of surviving species p * (α) (left) and second moment m * (α) (right) computed as solutions of (3.8)-(3.9), and their empirical Monte Carlo counterpart (p L (α), mL (α)). The parameter α on the x-axis ranges from 1 to 2 log(N ) 3.71. The value α = √ 2 corresponds to the theoretical lower bound provided by Theorem 3.3, which guarantees a stable equilibrium; the value α = 2 log(N ) is the upper-limit above which we have no extinction (p * = 1). Notice that for α ∈ (1, √2], the simulations show a remarkable matching with the heuristics despite no theoretical guarantee.

  royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (a) Gaussian entries. (b) Uniform entries.

Figure 5 : 3 .

 53 Figure 5: Distribution of surviving species. The solid line represents the theoretical distribution f * as given by Conjecture 3.7. The histogram is built by solving the LCP problem with an interaction matrix of size N = 2000 and parameter set to α = √ 3. In plot (A), the entries are Gaussian N (0, 1) and the interaction strength is fixed to α = 2. In plot (B), the entries are uniform U([-√ 3, √3]) which implies that the entries are centered with variance one.

Conjecture 3 . 8 .

 38 Let Γ be a non-centered Ginibre matrix, that is obeying model (2.5) with ξ = 0. Assume that α ∈ √ 2, 2 log(N ) and, for any i ∈ [N ], r i = 1. Then, for any i, j ∈ S, i.e. i, j correspond to royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 6 :

 6 Figure 6: Spectrum of the interaction matrix for N = 1000 species. ∆ is ER symmetric with C = 0.1. For the competitive (resp. mutualistic) model, A ij = A ji with distribution -|N (0, 1)| (resp. |N (0, 1)|) variables. For predator-prey A ij = -A ji , distribution |N (0, 1)|.

Figure 7 :

 7 Figure 7: Deterministic model with N = 1000, d = 10. Each species interacts with d species. Nonnull interactions are i.i.d. N (0, 1). (a) Spectrum of the 1 √ d -normalized interaction matrix. (b) Proportion of feasible equilibrium: each point represents the proportion of feasible solutions x over 200 realizations of random matrices Γ N for different values of κ, with α N = κ log(N ).

  royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. 10 )

 10 where σ i , b i , d i and Γ ij do not depend on K. As usual, we denote by b = (b i ), d = (d i ), r = (b i -d i ) N × 1 vectors and define the matrix Σ as the N × N diagonal matrix with entries (σ i ),

[ 4 ,

 4 Part 4(a)] for details. (ii) Variations around the Feller equations In [123], Biroli et al. added an immigration factor λ > 0. The SDEs they consider are:dX(t) = X(t) • (1 N -X(t) + Γ X(t)) dt + λ1 N dt + f (X(t)) • dB(t),(5.12)

  used the approach of Ben Arous and Guionnet to study a diffusion version of the so-called Hopfield model for biological neural networks. Details are given in [4, Part 4(b)].

  royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	May's intuition: assuming an equilibrium with random-like features, it should generically be unstable for high diversity
		stable	unstable	of interactions and species) (number and variance effective diversity
		May threshold
	What May did not discuss:	before Dynamics leading up to equilibrium: -other qualitative transitions? -do we lose coexistence? -does structure emerge?	after -do species go extinct, reducing diversity? -do species survive, but in a fluctuating state?

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	Using φ LV i (x) := r i -x i + (Γ x) i introduced in (1.2), Eq. (2.1) follows the generic form
	dx i dt	= x i φ LV i (x) ,
	introduced in (1.1) with φ LV i

  The following system of two equations and two unknowns (p, m)

		2, 2 log(N ) and assume
	that Γ follows Model (2.4)-(ii).
	2
	2 2π √	du .

Table 1 :

 1 The table above summarizes the different regimes (depending on α) of the LV system (2.1) where r i = 1 and Γ follows i.i.d. model (2.4)-(ii) with standard Gaussian entries. By "OQ?", we mean open question; in typewriter font, the result is proven mathematically; in italics, the result is established at a physical level of rigor.

	2 log(n)

  royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Denote by Y K i (t) the size of species i at time t and by Y K

  t≥0 , royalsocietypublishing.org/journal/rspa Proc R Soc A 0000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . for i ∈ [N ], whose evolution is described by the system of ODEs (2.1):

Beware that in RMT, universality has a different meaning: for a given random matrix, a result is universal if it does not depend on the distribution of the underlying random variables but only on its first moments.

May's own interpretation was probably focused on the third (Nature had to employ "devious strategies", in his words, to allow high-diversity ecosystems to exist in a stable equilibrium).

That is the empirical covariance matrix of the vectors X1, . . . , X N .

One can think of fruit and beak sizes in a plants-birds interaction network.

Unique stable equilibrium

Acknowledgements. All authors are supported by the CNRS 80 prime project KARATE. I.A., M.C., W.H., M.M., J.N. are supported by GdR MEGA. I.A. and M.M. are funded by Labex CEMPI (ANR-11-LABX-0007-01). M.C., W.H., J.N. and V.C.T. are supported by Labex Bézout (ANR-10-LABX-0058). F.M. and V.C.T. are funded by ANR EcoNet (ANR-18-CE02-0010). I.A. and V.C.T. acknowledges support of the "Chaire Modélisation Mathématique et Biodiversité (MMB)" of Veolia Environnement-Ecole Polytechnique-Museum National d'Histoire Naturelle-Fondation X. V.C.T. is partly funded by the European Union (ERC-AdG SINGER-101054787). M.M. and V.C.T. thank for its hospitality the CRM Montréal (IRL CNRS 3457) where part of this work was completed.

sequential invasion model, Bastolla et al. [143] found that the evolution of food webs enforces constraints on species diversity through competition among species preying on the same prey species and through an increase in competition strength at higher trophic levels due to the propagation of fluctuations and energy dissipation. In food web models incorporating realistic, non-invasion-based evolution, species tend to naturally group into trophic levels [START_REF] Allhoff | Evolutionary food web model based on body masses gives realistic networks with permanent species turnover[END_REF]144,145], an absent feature in most random matrix models. The various structures described in section 4 can be used to investigate their effects on equilibria and their stability, but it is still a long way from assessing the effects of evolutionary processes on these properties, e.g. the effect of the speciation-to-invasion ratio on food web susceptibility to further invasions [146].

(b) Open mathematical and modelling problems Multiple equilibria and non-equilibrium attractors

While our review focuses on conditions under which there exists a unique stable equilibrium (Sec. 3), other regimes including multiple equilibria or non-equilibrium attractors (e.g. chaos or cycles) have been studied using physics tools [START_REF] Roy | Complex interactions can create persistent fluctuations in high-diversity ecosystems[END_REF]123,147], see for instance Fig. 9. A mathematical understanding of these regimes and the associated thresholds is challenging and would represent an important step in the understanding of LV systems.

Quantitative metrics for stability

As hinted in Sec. 5 and [4, Sec. 3], other dynamical properties of LV systems (e.g. dynamics of fluctuations around equilibria), can also be studied using stochastic differential equation analogues to the LV ODEs. More quantitative metrics of stability are also of ecological interest, e.g. quantifying how much attractors change in response to perturbations of model parameters [149], assessing how much model parameters can change without changing the dynamical attractor (i.e. structural stability, [START_REF] Saavedra | A structural approach for understanding multispecies coexistence[END_REF]), or characterizing the time and trajectory of return to an attractor after a perturbation in the abundances (e.g. chap. 58 in [150]). When perturbations are small, one can linearize the dynamics around an attractor, and compute stability metrics for the i.i.d. or elliptic model [START_REF] Arnoldi | Unifying dynamical and structural stability of equilibria[END_REF]151,152]. Beyond the linearized regime, tools such as Freidlin-Wentzell theory can help quantify the basins of dynamical attractors and transition times between them [153].

Counting the number of coexisting species

In Section 3(e), we provide heuristics to compute the proportion of vanishing species for a given equilibrium and refer to articles relying on techniques from Physics. A mathematical computation still remains out of reach and could be phrased as understanding the properties of solutions of Linear Complementarity Problems with random matrices as input. Part of the challenge lies on the fact that the properties of interest of the solution of a LCP (proportion of non-vanishing components) do not only rely on the spectral properties of the random matrices at hand. Table 1 in Section 3 lists a number of open questions when studying a specific LV model.

A mathematical understanding of the cavity method

The cavity method, originally introduced in the context of spin glasses [125], has become extremely versatile and in general yields to closed-form expressions accurately matching simulation results. There are quite few mathematically rigorous counterparts of the results obtained through the cavity method (one can cite e.g. [154] for sparse graphs). In the context of ecosystems, we have worked in this direction in [155,156], using mathematical results on the Approximate Message Passing algorithm to recover rigorously some equations obtained in [START_REF] Bunin | Ecological communities with Lotka-Volterra dynamics[END_REF] through the cavity method. There is still much to do to get a wider mathematical formalization of the cavity method for the analysis of non-spectral properties of large random matrix observables.