domain: math sdv

text:
Maxime Clenet

François Massol

Jamal Najim

Impact of a block structure on the Lotka-Volterra model

Keywords: Lotka-Volterra model, Block structure, Linear Complementarity Problems, Large Random Matrices, Stability of food webs

The Lotka-Volterra (LV) model is a simple, robust, and versatile model used to describe large interacting systems such as food webs or microbiomes. The model consists of n coupled differential equations linking the abundances of n different species. We consider a large random interaction matrix with independent entries and a block variance profile. The ith diagonal block represents the intra-community interaction in community i, while the off-diagonal blocks represent the inter-community interactions. The variance remains constant within each block, but may vary across blocks.

We investigate the important case of two communities of interacting species, study how interactions affect their respective equilibrium. We also describe equilibrium with feasibility (i.e., whether there exists an equilibrium with all species at non-zero abundances) and the existence of an attrition phenomenon (some species may vanish) within each community.

Information about the general case of b communities (b > 2) is provided in the appendix.

Introduction

Understanding large ecosystems and the underlying mechanisms that support high species diversity is a major challenge in theoretical ecology. Motivated by the seminal work of May [START_REF] May | Will a Large Complex System be Stable?[END_REF], the introduction of random matrices has been a key mathematical step in modeling high-dimensional ecosystems. These tools have expanded our ability to understand the nature of interactions and how food webs can resist to small perturbations (stability) [START_REF] Allesina | Stability criteria for complex ecosystems[END_REF][START_REF] Tang | Correlation between interaction strengths drives stability in large ecological networks[END_REF].

Differential equations are frequently used in ecology to describe a system of interacting species. One of the most common models is the Lotka-Volterra (LV) model [START_REF] Lotka | Elements of Physical Biology[END_REF][START_REF] Volterra | Fluctuations in the Abundance of a Species considered Mathematically[END_REF], which has been the subject of research in both ecology [START_REF] Wangersky | Lotka-Volterra Population Models[END_REF][START_REF] Jansen | A permanence theorem for replicator and Lotka-Volterra systems[END_REF][START_REF] Law | Self-Assembling Food Webs: A Global Viewpoint of Coexistence of Species in Lotka-Volterra Communities[END_REF] and mathematics [START_REF] Goh | Feasibility and stability in randomly assembled Lotka-Volterra models[END_REF][START_REF] Goh | Global Stability in Many-Species Systems[END_REF][START_REF] Taylor | Consistent scaling and parameter choice for linear and Generalized Lotka-Volterra models used in community ecology[END_REF][START_REF] Hofbauer | Evolutionary Games and Population Dynamics[END_REF][START_REF] Takeuchi | Global dynamical properties of Lotka-Volterra systems[END_REF]. Certain properties of this model, such as its stability [START_REF] Gibbs | Effect of population abundances on the stability of large random ecosystems[END_REF], have raised much interest. The conditions under which all species survive, referred to as feasibility, have also motivated many works [BN21, GAS + 17, Sto18].

In nature, ecological networks are rather structured, and many studies have investigated the network structures that contribute to the stability of a given community [TF10, AGB + 15]. One common network structure is food web compartmentalization, also known as modularity. The underlying concept is that the network is structured in the form of groups of nodes that interact more strongly within their group and more weakly between groups. A mathematical formulation of modularity was defined by Newman [START_REF] Newman | Modularity and community structure in networks[END_REF]. Subsequently, modularity has been of great importance in ecology [GSSP + 10], in complex networks [START_REF] Variano | Networks, Dynamics, and Modularity[END_REF], and in community detection (for a complete review, see Fortunato [START_REF] Fortunato | Community detection in graphs[END_REF]).

May had already mentioned that a multi-community structure should improve stability [START_REF] May | Will a Large Complex System be Stable?[END_REF], a hypothesis later investigated by Pimm [START_REF] Pimm | The structure of food webs[END_REF]. In the same framework as May, Grilli et al. [START_REF] Grilli | Modularity and stability in ecological communities[END_REF] studied the effect of modularity on the stability of the Jacobian of a system, the so-called "community matrix". However, studies show that modularity improves the persistence (:= non-extinction of species, generally related to their resistance to external perturbations) of species in the dynamical system [START_REF] Stouffer | Compartmentalization increases food-web persistence[END_REF].

In this article, we study the Lotka-Volterra model where we consider a block structure network representing the inter-and intra-community interactions. Of particular interest are the interactions between the communities that affect their respective equilibrium and stability.

Each block is identified by its interaction strength, which is the standard deviation of the random part of the interactions. The idea that interaction strength plays a key role in the stability of ecosystems was introduced by May [START_REF] May | Will a Large Complex System be Stable?[END_REF]. For the sake of mathematical simplicity, we limit our model to two communities, although we can extend the model to more complex food webs and multi-community frameworks.

We study the existence and stability of an equilibrium, together with its properties. When an attrition phenomenon occurs (some species may vanish), we describe the proportion of surviving species and their distribution. We also provide conditions for which the equilibrium is feasible (i.e., whether there exists an equilibrium with all species at non-zero abundances).

Model and assumptions. The LV model is a standard model in ecology to study the dynamics of a community of species over time. It is defined by a system of n differential equations

dx k dt (t) = x k (t)   r k -x k (t) + ℓ∈[n] B kℓ x ℓ (t)   , k ∈ [n] = {1, • • • , n} . (1)
The abundance of species k at time t is represented by x k (t) with x = (x 1 , • • • , x n ) the vector of abundances. Parameter r k corresponds to the growth rate of species k. The coefficient B kℓ represents the impact of species ℓ on species k. The n × n matrix B = (B kℓ ), which represents the interaction network, is decomposed into a block structure. This structure differentiates various groups of species in the form of communities that interact with each other. On the one hand, the diagonal blocks of B correspond to interactions within each community, each with its own interaction strength. On the other hand, the off-diagonal blocks correspond to the impact of the communities on each other. Analytically and within the framework of two communities, the matrix B = (B kℓ ) n,n is defined by blocks using random matrices A ij and interaction strengths s ij :

B = 1 √ n s 11 A 11 s 12 A 12 s 21 A 21 s 22 A 22 , (2)
where I 1 = [n 1 ] (resp I 2 = {n 1 +1, • • • , n}), the subset of [n] of size |I 1 | = n 1 (resp |I 2 | = n 2 -here and below n = n 1 + n 2 ) matching the index of species belonging to community 1 (resp community 2) and β = (β 1 , β 2 ) with β 1 = n 1 /n , β 2 = 1 -β 1 = n 2 /n. The random matrix A ij is non-Hermitian of size |I i | × |I j | with standard Gaussian entries i.e. N (0, 1). The Gaussianity assumption clarifies the explanations, but can be relaxed under certain circumstances (see the corresponding sections for details). Notice a normalization parameter 1/ √ n in the matrix B. This enables the interaction matrix B to have a macroscopic effect on system (1): From an ecological perspective, an increase in the number of species may not necessarily lead to a corresponding increase in the overall strength of interactions between one species and all others.

E   ℓ∈[n]
The relative strength of interactions within and between blocks is controlled by the four s ij coefficients, which can be grouped together in a matrix s:

s = s 11 s 12 s 21 s 22 .
The diagonal terms (s 11 , s 22 ) represent the interaction strength in each community. The off-diagonal term s 12 (resp. s 21 ) represents the interaction strength of the impact of community 2 on community 1 (resp. community 1 on community 2). The lower the value of s ij , the lower the rates of interaction between species. Note that in the case of a unique community, s is the interaction strength coefficient, i.e. the standard deviation of the interspecific coefficients of the LV model.

Remark 1. For the sake of simplicity, the results are presented in the case of two interacting communities but can be extended to the case of b communities, see Appendix E.

There are two scenarios of interest: Let us consider two separate groups of species that follow the dynamics described in model (1). The matrices A ij are each sampled once. In the first scenario, we consider a very weak interaction between the two communities

s = 1/2 ε ε 1/2 , ε > 0 ,
and the intra-communities interactions are small enough (Fig 1a). We observe that both communities dynamics converge to a feasible equilibrium in the sense that all species survive (Fig. 1b). In the second scenario, we increase the interactions between the communities, i.e. the standard deviation matrix is defined by

s = 1/2 1 1 1/2 .
It is no longer possible for both communities to maintain the feasibility of all species. Some species are likely to disappear (Fig. 2).

Properties of the dynamical system. We are interested in the effect of a block structure on the food web, limit our study to the 2-blocks case (2) and focus on the model with constant growth rate1 r k = 1:

dx k dt = x k (1 -x k + (Bx) k ) , k ∈ [n] . (3)
Of major interest is the existence and uniqueness of an equilibrium

x * = (x * k ) k∈[n]
. The LV system is an autonomous differential system. If the initial conditions are positive i.e. x(0) > 0 (componentwise), it implies x(t) > 0 for every t > 0. However, some of the components x k (t) may converge to zero if the equilibrium x * has components equal to zero. An equilibrium to the LV system should hence satisfy the following set of constraints:

x

* k (1 -x * k + (Bx * ) k ) = 0 , ∀k ∈ [n] , x * k ≥ 0 . (4)
Two substantially different situations arise, that we will study hereafter. First, if x * has vanishing components, the equilibrium equations are cast into a nonlinear optimization problem, which has been studied by Clenet et al. [START_REF] Clenet | Equilibrium and surviving species in a large Lotka-Volterra system of differential equations[END_REF] in the case of a single community.

If the equilibrium is feasible, that is x * > 0, then the equilibrium set of equations becomes a linear equation:

x

* = 1 + Bx * . (5)
In the context of a single community, the existence of a positive solution has been studied by Bizeul and Najim [START_REF] Bizeul | Positive solutions for large random linear systems[END_REF] and extended for more complex food webs in [START_REF] Akjouj | Feasibility of sparse large Lotka-Volterra ecosystems[END_REF][START_REF] Clenet | Equilibrium in a large Lotka-Volterra system with pairwise correlated interactions[END_REF][START_REF] Liu | Feasibility and stability in large Lotka Volterra systems with interaction structure[END_REF].

A further consideration which will be addressed is whether the equilibrium x * is asymptotically globally stable, i.e. if for every initial vector x 0 > 0 the solution of (3), which starts at x(0) = x 0 , satisfies

x(t) ---→ t→∞ x * .
In the sequel, the term "stability" will refer to "asymptotic stability".

Outline of the article. In Section 2, we describe sufficient conditions for the existence and uniqueness of a stable equilibrium in the model (3), see Theorem 3. Section 3 is devoted to the study of the properties of the species that survive in each of the communities. Finally, in Section 4 we provide conditions under which the equilibrium is feasible, see Theorem 4.

Existence of a unique equilibrium

In Figures 1 and2, we notice that for different interaction coefficients s, the system converges to an equilibrium (with or without vanishing species). Theorem 3 below will provide the adequate theoretical framework.

Theoretical background

Non-invadability condition. The research of equilibrium points of (3) is equivalent to the identification of solutions of system (4). However, the number of potential solutions can be extremely large. In order for the equilibrium x * to be stable, there exists a necessary condition, known in ecology as the non-invadability condition [START_REF] Law | Permanence and the Assembly of Ecological Communities[END_REF], namely that 1

-x * k + (Bx * ) k ≤ 0 . (6)
In model (3), the non-invadability condition for a given species k ∈ [n] is equivalent to

1 x k dx k dt x k →0 + ≤ 0 . (7)
Condition (7) describes the fact that if we add a species to the system at a very low abundance, it will not be able to invade the system. As a consequence, the number of possible solutions should solve the following set of constraints:

   x * k (1 -x * k + (Bx * ) k ) = 0 for k ∈ [n] , 1 -x * k + (Bx * ) k ≤ 0 for k ∈ [n] , x * ≥ 0 componentwise . (8)
This casts the search of a nonnegative equilibrium problem into the class of linear complementarity problems (LCP). For a reminder of the definition of an LCP problem, see for instance [START_REF] Clenet | Equilibrium and surviving species in a large Lotka-Volterra system of differential equations[END_REF]. In the following, we recall the main Theorem for proving the existence and uniqueness of a single equilibrium.

The equilibrium x * and its stability. Let X ⊤ be the transpose of the matrix X.

Definition 1 (Lyapunov diagonal stability). A matrix M is called Lyapunov diagonally stable, denoted by M ∈ D, if and only if there exists a diagonal matrix D with positive diagonal elements such that DM +M ⊤ D is negative definite, i.e. all eigenvalues are negative. This class of matrix was already mentioned in Volterra's historical paper [START_REF] Volterra | Leçons sur la théorie mathématique de la lutte pour la vie[END_REF] and in Logofet's book [START_REF] Logofet | Matrices and graphs: stability problems in mathematical ecology[END_REF]Chap. 4], in relation with the stability of LV models.

Proposition 1 (Takeuchi et al. [START_REF] Takeuchi | Global stability of ecosystems of the generalized volterra type[END_REF]). If M ∈ D then -M is a P-matrix i.e. all its principal minors (sub-determinants) are strictly positive:

det(M I ) > 0 , ∀ I ⊂ [n] , M I = (M kℓ ) k,ℓ∈I .
Recall System (1) with different growth rates for each species and consider matrix B is arbitrary,

d y k dt = y k (r k + ((-I + B)y) k ) , k ∈ [n] . (9)
The LCP associated with (9) is as follows

   y * k (r k -y * k + (By * ) k ) = 0 for k ∈ [n] , r k -y * k + (By * ) k ≤ 0 for k ∈ [n] , y * ≥ 0 componentwise . (10)
Takeuchi and Adachi [Tak96, Th. 3.2.1] establish a criterion for the existence of a unique globally stable equilibrium y * of (9).

Theorem 2 (Takeuchi and Adachi [START_REF] Takeuchi | The existence of globally stable equilibria of ecosystems of the generalized Volterra type[END_REF]). If -I + B ∈ D, then System (10) admits a unique solution. In particular, for every r ∈ R n , there is a unique equilibrium y * to (9), which is globally stable in the sense that for every y 0 > 0, the solution to (9) which starts at y(0) = y 0 satisfies: y(t) ---→ t→∞ y * .

Sufficient condition in the block model (3)

We recall the definition of the Stieltjes transform and some of its properties in Appendix A (for more details, see [START_REF] Bai | Spectral analysis of large dimensional random matrices[END_REF]). For a wide range of parameters (β, s), we aim to ensure the existence of a globally stable equilibrium x * of (3) associated to LCP (8). Denote by ∥x∥ ∞ the sup norm of a vector and by ∥X∥ ∞ its induced operator norm, i.e.

∥x∥ ∞ = max k∈[n] |x k | and ∥X∥ ∞ = max k∈[n] n ℓ=1 |X kℓ | .
Let X, Y be matrices of the same size, then

X • Y is their Hadamard product i.e. (X • Y ) ij = X ij Y ij and consider s • s = s 2 11 s 2 12 s 2 21 s 2 22 . Theorem 3. Assume that diag(β) (s • s) + (s • s) ⊤ ∞ < 1 , then a.s. matrix (I -B) + (I -B) ⊤ is
eventually positive definite: with probability one, there exists N depending on matrix B's realization such that for n ≥ N , (I -B) + (I -B) ⊤ is positive definite. In particular, I -B ∈ D. There exists a unique vector solution to the LCP (8). This vector x * is the unique (random) globally stable equilibrium of (3).

Remark 2. Theorem 3 can be extended in two directions. The Gaussianity assumption can be relaxed to any reasonable distribution with finite second moment, and growth rates different from one i.e. r k ̸ = 1 can be considered (see for instance [BN21, Sections 4.2 and 4.3] in the context of a single community).

Sketch of proof.

From Theorem 2, we need to verify the Lyapunov diagonally stable condition of the matrix (-I + B) ∈ D by analyzing its largest eigenvalue

(-I + B) + (-I + B ⊤ ) = -2I + 1 √
). Given m(z) = (m 1 (z), • • • , m n (z)
), the QVE associated to the matrix H is decomposed as

k ∈ I 1 , - 1 m k (z) = z + ℓ∈I1 2s 2 11 n m ℓ (z) + ℓ∈I2 1 n s 2 12 + s 2 21 m ℓ (z) , k ∈ I 2 , - 1 m k (z) = z + ℓ∈I1 1 n s 2 12 + s 2 21 m ℓ (z) + ℓ∈I2 2s 2 22 n m ℓ (z) . Denote by 1/m(z) = (1/m 1 (z), • • • , 1/m n (z)) ⊤ , 1 Ii a vector whose entries are 1's of size |I i | and S = 1 n 2s 2 11 1 I1 1 ⊤ I1 (s 2 12 + s 2 21 )1 I1 1 ⊤ I2 (s 2 12 + s 2 21 )1 I2 1 ⊤ I1 2s 2 22 1 I2 1 ⊤ I2 ,
the QVE can be written in the standard form ∞ . This information gives an asymptotic bound on the support of the matrix H associated with (11), i.e. asymptotically ∀ ε > 0 there exists N depending on matrix B's realization, such that for n ≥ N

- 1 m(z) = z + Sm(z) . (11
λ max (H) ≤ 2 ∥S∥ 1/2 ∞ + ε . Recall that -2I + H is negative definite iff λ max (H) < 2. This condition is fulfilled if 2 ∥S∥ 1/2 ∞ < 2 , or equivalently ∥S∥ ∞ < 1 .
Note that this condition is sufficient but not necessary. Given the particular shape of the matrix S, computing its norm is equivalent to computing the norm of a matrix of size 2 × 2

∥S∥ ∞ = diag(β) (s • s) + (s • s) ⊤ ∞ = β 1 0 0 β 2 2s 2 11 s 2 12 + s 2 21 s 2 12 + s 2 21 2s 2 22 ∞
, which completes the proof. We can then rely on Theorem 2 to conclude.

Remark 3.

1. In the context of a unique community, suppose that s = s11 ⊤ , then the previous condition takes the simpler form s < 1/ √ 2, a condition already mentionned in [START_REF] Clenet | Equilibrium and surviving species in a large Lotka-Volterra system of differential equations[END_REF]. Indeed, starting from the condition of Theorem 3, the condition on the matrix is

1/2 0 0 1/2 1/2 2s 2 2s 2 2s 2 2s 2 1/2 0 0 1/2 1/2 ∞ = 1s 2 1s 2 1s 2 1s 2 ∞ = 2s 2 ,
The same sufficient condition is obtained 2s

2 < 1 ⇔ s < 1/ √ 2.
2. The condition given in Theorem 3 is sufficient to guarantee a stable unique solution to LCP (8) but not necessary. Even in the single community case, finding optimal thresholds remains an open question (see [ABC + 22, Table 1]).

Surviving species

In Section 2, we have given conditions on matrix s and on β = (β 1 , β 2 ) for the existence of a globally stable equilibrium x * to (3) under the non-invadability condition. The equilibrium vector x * is random and depends on the realization of matrix B. Moreover since s has fixed components and does not depend on n, the equilibrium x * will feature vanishing components (see the original argument for a unique community in [DVR + 18] and the discussion in [START_REF] Bizeul | Positive solutions for large random linear systems[END_REF]). In an ecological context, we differentiate two kind of components in vector x * , the non-vanishing components x * k > 0 corresponding to surviving species and the vanishing ones x * k = 0 corresponding to the species going to extinction:

x k (t) ---→ t→∞ 0 .
Hereafter, we describe statistical properties of x * : the proportion of surviving species in each community, the distribution of the corresponding abundances, which turns out to be a truncated Gaussian, etc.

Heuristics for the properties of surviving species

Starting from the model (3), the set of surviving species in community i ∈ {1, 2} is defined as

S 1 = {k ∈ I 1 , x * k > 0} ; I 1 = [1, β 1 n] , S 2 = {k ∈ I 2 , x * k > 0} ; I 2 = [β 1 n + 1, n] .
Given the random equilibrium x * , we introduce the following quantities for each community

i ∈ {1, 2} pi = |S i | |I i | , mi = 1 |S i | k∈Ii x * k , σ2 i = 1 |S i | k∈Ii (x * k ) 2 .
Quantity pi represents the proportion of surviving species in community i, mi the empirical mean of the abundances of the surviving species in community i and σ2 i , the empirical mean square of the surviving species in community i.

Denote by Z ∼ N (0, 1) a standard Gaussian random variable and by Φ the cumulative Gaussian distribution function:

Φ(x) = x -∞ e -u 2 2 √ 2π du .
Heuristics 1. Let s be the 2 × 2 matrix of interaction strengths and assume that the condition of Theorem 3 holds, then the following system of four equations and four unknowns (p 1 , p 2 , σ 1 , σ 2 )

p 1 = 1 -Φ(δ 1 ) , (12)
p 2 = 1 -Φ(δ 2 ) , (13)
(σ 1 ) 2 = 1 + 2∆ 1 E(Z|Z > δ 1 ) + ∆ 2 1 E(Z 2 |Z > δ 1 ) , (14)
(σ 2 ) 2 = 1 + 2∆ 2 E(Z|Z > δ 2 ) + ∆ 2 2 E(Z 2 |Z > δ 2 ) , (15)
where ----→ n→∞ σ * i .

∆ i = p 1 (σ 1 ) 2 β 1 s 2 i1 + p 2 (σ 2 ) 2 β 2 s 2 i2 and δ i = -1 ∆ i , (16)
In order to simplify the following calculations, we denote by

∆ * i = p * 1 (σ * 1 ) 2 β 1 s 2 i1 + p * 2 (σ * 2 ) 2 β 2 s 2 i2 and δ * i = -1 ∆ * i .
There is a strong matching between the solutions obtained by solving ( 12)-( 15) and their empirical counterparts obtained by Monte-Carlo simulations. This is illustrated in Fig. 3.

s 12 s 21 1/ √ 2 , β = 1 2 , 1 2 .
When off-diagonal block interactions s 12 , s 21 increase, the proportion of surviving species p * decreases but the root mean square of their equilibrium abundances σ * increases.

Construction of the heuristics

Obtaining information about the fixed point is equivalent to solving the LCP problem

x * k 1 -x * k + n ℓ=1 B kℓ x * ℓ = 0 , ∀k ∈ [n] .
Consider the random variables:

∀k ∈ [n], Žk = ℓ∈S1∪S2 B kℓ x * ℓ .
We assume that asymptotically the x * ℓ 's are independent from the B kℓ 's, an assumption supported by the chaos hypothesis, see for instance Geman and Hwang [START_REF] Geman | A chaos hypothesis for some large systems of random equations[END_REF]. Denote by E x * = E( • | x * ) the conditional expectation with respect to x * . Notice that conditionally to x * , the Žk 's are independent Gaussian random variables, whose first two moments can easily be computed, see Appendix C for the details:

∀k ∈ I i , Var x * ( Žk ) = p1 σ2 1 β 1 s 2 i1 + p2 σ2 2 β 2 s 2 i2 .
Notice that the fact that Var x * ( Žk ) only depends on p1 , p2 , σ1 , σ2 (which are converging quantities when n → ∞) supports the idea that Žk is unconditionally a Gaussian random variable with second moment:

Var( Žk ) = p * 1 (σ * 1 ) 2 β 1 s 2 i1 + p * 2 (σ * 2 ) 2 β 2 s 2 i2 , where p * 1 , p * 2 , σ * 1 , σ * 2 are resp.
the limits of p1 , p2 , σ1 , σ2 . We can introduce two families of standard Gaussian random variables (Z k ) k∈I1 and (Z k ) k∈I2 :

∀k ∈ I i , Z k = Žk var( Žk ) = Žk p * 1 (σ * 1 ) 2 β 1 s 2 i1 + p * 2 (σ * 2 ) 2 β 2 s 2 i2 . Consider the equilibrium x * = (x * k ) k∈[n]
, the definition of the LCP equilibrium implies if k ∈ S 1 ∪ S 2 :

x * k (1 -x * k + (Bx * ) k ) = 0 and 1 + (Bx * ) k = 1 + Žk > 0 .
We finally obtain the following relationship for the surviving species:

x * k = 1 + ∆ * i Z k for k ∈ S i . (17)
Note that ∆ * i corresponds to the average variance of the interactions on community i.

Heuristics (12)-(13). We can write the first two equations:

P(x * k ≥ 0 | k ∈ S 1 ) = P(Z k < δ * 1 | k ∈ S 1 ) = 1 -Φ(δ * 1 ) , P(x * k ≥ 0 | k ∈ S 2 ) = P(Z k < δ * 2 | k ∈ S 2 ) = 1 -Φ(δ *
2 ) , and finally obtain (12) and (13):

p * 1 = 1 -Φ(δ * 1 ) and p * 2 = 1 -Φ(δ * 2 ) .
Heuristics (14)-(15). Our starting point is the following generic representation of an abundance at equilibrium (either of a surviving or vanishing species) in the case k ∈ S i :

x * k = (1 + ∆ * i Z k ) 1 {Z k >δ * i } = 1 {Z k >δ * i } + (∆ * i Z k ) 1 {Z k >δ * i } .
Taking the square, we get:

(x * k ) 2 = (1 + ∆ * i Z k ) 2 1 {Z k >δ * i } = 1 {Z k >δ * i } + 2∆ * i Z k 1 {Z k >δ * i } + (∆ * i ) 2 Z 2 k 1 {Z k >δ * i } .
Summing over S i and normalizing, we get

1 |S i | k∈Si (x * k ) 2 = 1 |S i | k∈Si 1 {Z k >-δ * i } + 2∆ * i 1 |S i | k∈Si Z k 1 {Z k >-δ * i } + (∆ * i ) 2 1 |S i | k∈Si Z 2 k 1 {Z k >δ * i } , σ2 i (a) = 1 + 2∆ * i |I i | |S i | 1 |I i | k∈Ii Z k 1 {Z k >δ * i } + (∆ * i ) 2 |I i | |S i | 1 |I i | k∈Ii Z 2 k 1 {Z k >δ * i } , σ2 i (b) ≃ 1 + 2∆ * i 1 P(Z > δ * i ) E(Z1 {Z>δ * i } ) + (∆ * i ) 2 1 P(Z > δ * i ) E(Z 2 1 {Z>δ * i } ), σ2 i ≃ 1 + 2∆ * i E(Z | Z > δ * i ) + (∆ * i ) 2 E(Z 2 | Z > δ * i ).
where (a) follows from the fact that

|S i | = k∈Si 1 {Z k >δ * i } (by definition of S i ), (b) from the law of large numbers 1 |Ii| k∈Ii Z j k 1 {Z k >δi} ----→ n→∞ EZ j 1 {Z>δ * i } , j = 1, 2 and |Si| |Ii| ----→ n→∞ P(Z > δ * i )
with Z ∼ N (0, 1). It remains to replace σi by its limit σ * i to obtain ( 14)-(15). We finally obtain the third and fourth equations:

(σ * 1 ) 2 = (1 + λ * 1 ) 2 + 2(1 + λ * 1 )∆ * 1 E(Z | Z > δ * 1 ) + (∆ * 1 ) 2 E(Z 2 | Z > δ * 1 ) , (σ * 2 ) 2 = (1 + λ * 2 ) 2 + 2(1 + λ * 2 )∆ * 2 E(Z | Z > δ * 2 ) + (∆ * 2 ) 2 E(Z 2 | Z > δ * 2 ) .

General properties of the ecosystem

The properties at equilibrium, such as the proportion and mean square of the abundance of surviving species, can be computed for each community by solving the system of equations in Heuristics 1. An additional property, the mean abundance of the surviving species at equilibrium for each community (m * 1 , m * 2 ), can be calculated using a method similar to the mean square of the abundances (see Appendix C.2 for the details of the computations).

m * 1 = 1 + ∆ * 1 E(Z|Z > δ * 1 ) , (18)
m * 2 = 1 + ∆ * 2 E(Z|Z > δ * 2 ) . ( 19
)
The two equations are not necessary for solving Heuristics 1, but they provide new information. In particular, strong inter-or intra-community interactions increase the mean abundance of the surviving species (see Fig. 4). Conditional on each community, one can easily extend the properties of each community to the whole ecosystem. We denote by p * the proportion, (σ * ) 2 the mean square and m * the mean of surviving species. We observe the linear effect of community size β on general properties:

1. Proportion of surviving species.

P(x * k ≥ 0) = P(x * k ≥ 0 | k ∈ I 1 )P(k ∈ I 1 ) + P(x * k ≥ 0 | k ∈ I 2 )P(k ∈ I 2 ) , p * = p * 1 β 1 + p * 2 β 2 .
2. Mean square of the abundance of the surviving species.

E((x * k ) 2 ) = E((x * k ) 2 | k ∈ I 1 )P(k ∈ I 1 ) + E((x * k ) 2 | k ∈ I 2 )P(k ∈ I 2 ) , (σ * ) 2 = (σ * 1 ) 2 β 1 + (σ * 2 ) 2 β 2 .
3. Mean of the abundance of the surviving species

E(x * k ) = E(x * k | k ∈ I 1 )P(k ∈ I 1 ) + E(x * k | k ∈ I 2 )P(k ∈ I 2 ) , m * = m * 1 β 1 + m * 2 β 2 .

Distribution of the surviving species

We may recall the following representation of the abundance x * k of a surviving species when k ∈ S i :

x * k = 1 + ∆ * i Z k if k ∈ S i ,
where Z k ∼ N (0, 1) and

Z k > δ * i = δ i (p * i , σ * i ) defined in (17
). This representation allows to characterize the distribution of x * k of each community. It turns out that the surviving species of each community follow a truncated Gaussian distribution.

Heuristics 2. Let s be the 2 × 2 matrix of interaction strengths and assume that the condition of Theorem 3 holds and let (p * 1 , p * 2 , σ * 1 , σ * 2 ) be the solution of the system (12)-( 15). Recall the definition (16) of ∆ i and δ i and denote by

δ * i = δ i (p * i , σ * i ). Let x * k > 0 be a positive component of x * belonging to the community i, then the law of x * k is L(x * k ) ----→ n→∞ L (1 + ∆ * i Z | Z > δ * i ) ,
where Z ∼ N (0, 1). Otherwise stated, asymptotically for k ∈ S i , x * k admits the following density

f k (y) = 1 {y>0} Φ(-δ * i ) 1 ∆ * i √ 2π exp - 1 2 y ∆ * i + δ * i 2 . ( 20
)
The heuristics simply follows from the fact that if x * k is a surviving species and k ∈ S i then

x * k = 1 + ∆ * i Z k ,
conditionally on the fact that the right-hand side of the equation is positive, that is Z k > δ * i . A simple change of variable yields the density -details are provided in Appendix C.

Fig. 4 illustrates the matching between the theoretical distribution obtained by equation ( 20) and a histogram obtained by generating the interaction matrix for 2 communities. In Fig. 5, the validity of heuristics in the case of non-Gaussian entries is illustrated.

Remark 4. The proof relies on the Gaussianity assumption, but we are convinced that it could be extended beyond. In particular, in Figure 5, non-Gaussian entries centered E(B kℓ ) = 0 with variance one E(|B kℓ | 2 ) = 1 are considered. The distribution of surviving species still fits the truncated Gaussian in this case. and s = 1/2 2/3 1/3 1/4 . Notice in particular that the histogram is well predicted by the theoretical distributions even if the entries are non-Gaussian.

Feasibility

Recall s the interaction parameter in the case of a unique community. According to the work of Dougoud et al.. [DVR + 18], if s is fixed (i.e. does not depend on n) then there can be no feasible equilibrium at large n. Following this work, Bizeul and Najim [START_REF] Bizeul | Positive solutions for large random linear systems[END_REF] provided the appropriate normalization of s to have a feasible equilibrium. The threshold corresponds to s ∼ 1/ 2 log(n). The equilibrium is feasible almost surely when s is less than this threshold value, i.e. when elements of random matrix B are divided by 2n log(n) or a larger factor. Some extensions of these results have been made in the sparse case [START_REF] Akjouj | Feasibility of sparse large Lotka-Volterra ecosystems[END_REF] and with a mean and pairwise correlated entries [START_REF] Clenet | Equilibrium in a large Lotka-Volterra system with pairwise correlated interactions[END_REF]. In this section, conditions are given on the matrices s to get a feasible equilibrium in each community, called co-feasibility. We then provide some ecological interpretations.

Theoretical analysis of the threshold

Recall the notation x = (x k ) k∈[n] and denote by ∥x∥ ∞ = max k∈[n]
|x k |. We are interested in the existence of a feasible solution of the fixed point problem associated with the model (3). To consider this problem, we extend the computations of Bizeul and Najim in the framework of a block structure network. Consider s such that I -B is invertible. The problem is defined by

x * = 1 + Bx * ⇔ x * = (I -B) -1 1 , (21)
The problem (21) admits a unique solution. We consider a matrix s which depends on n, i.e. s = s n such that:

s n ----→ n→∞ 0 ⇔ ∀ i, j ∈ {1, 2} , s ij ----→ n→∞ 0 .
Note that for sufficiently large n, the problem satisfies the sufficient condition of Theorem 3 to have a unique globally stable equilibrium.

Let matrix B n depending on the interaction matrix s n defined by

B n = V s n V ⊤ • 1 √ n A 11 A 12 A 21 A 22 , (22)
where

V ∈ M n×2 , V = 1 I1 0 0 1 I2 .
The spectral radius of 1

√ n

A 11 A 12 A 21 A 22 a.s. converges to 1 due to the circular law [START_REF] Tao | Random matrices: Universality of ESDs and the circular law[END_REF]. So as long as s n is close to zero, the matrix I -B n is eventually (for large enough values of n) invertible.

Theorem 4 (Co-feasibility for the 2-blocks model). Assume that matrix B n is defined by the 2-blocks model (22). 1. If there exists ε > 0 such that eventually

Let β = (β 1 , β 2 ), β 1 + β 2 =
(s n • s n )β ⊤ ∞ ≥ (1 + ε)(s * n ) 2 then P min k∈[n]
x k > 0 ----→ n→∞ 0 .

2. If there exists ε > 0 such that eventually

(s n • s n )β ⊤ ∞ ≤ (1 -ε)(s * n ) 2 then P min k∈[n] x k > 0 ----→ n→∞ 1 .
A sketch of proof is postponed to Appendix D, the extension to the b-blocks case can be found in Appendix E.

Proof of Theorem 4 strongly depends on the assumption of Gaussianity and equal growth rates of each species. However, according to the approach of Bizeul et al. [START_REF] Bizeul | Positive solutions for large random linear systems[END_REF], these assumptions could be relaxed. In particular, the phenomenon seems to be universal, i.e. the feasibility threshold works for a wide range of distribution choices.

n the critical regime s ∝ 1/ log(n) or equivalently s -1 ∝ log(n). We thus introduce matrix κ defined by

κ = 1 log(n) s -1 11 s -1 12 s -1 21 s -1

22

.

Notice that at criticality κ will be of order O(1). This will be convenient for ecological interpretations. Using the inequality of Theorem 4, the co-feasibility condition on κ writes

(s n • s n )β ⊤ ∞ < (s * n ) 2 ⇔ max 2β 1 κ 2 11 + 2β 2 κ 2 12 , 2β 1 κ 2 21 + 2β 2 κ 2 22 < 1 . ( 23
)
If for i = 1, 2, β i = 1 2 and the entry of the matrix κ are equal, then condition (23) gives the threshold κ ij > √ 2, and we recover the same critical threshold 2 log(n) as in [START_REF] Bizeul | Positive solutions for large random linear systems[END_REF].

Remark 5. Assume κ 11 = κ 22 = ν 1 and κ 12 = κ 21 = ν 2 , condition (23) is reformulated as:

max 2β 1 ν 2 1 + 2β 2 ν 2 2 , 2β 1 ν 2 2 + 2β 2 ν 2 1 < 1 .
If β 1 , β 2 and ν 2 are fixed, then the phase transition on the intra-community interactions occurs at

ν 1 > min   β 1 1 2 -β2 ν 2 2 , β 2 1 2 -β1 ν 2 2   .
In Fig. 6, the phase transition is represented for a selected set of parameters. Note that the transition is rather smooth. The threshold depends on ν 2 . Increasing ν 2 (decreasing the inter-block interactions) lowers the co-feasibility threshold to at least 1 (for communities of the same size).

Preservation of co-feasibility

Equation (23) defines a "co-feasibility domain" and gives a constraint in five dimensions. The two communities of species can be studied independently i.e. the two components of equation ( 23) respectively give the feasibility condition for each community:

If 2β1 κ 2 11 + 2β2 κ 2 12 < 1, then community 1 is feasible. If 2β1 κ 2 21 + 2β2 κ 2 22
< 1, then community 2 is feasible.

The first community (resp. the second one) will be affected by changing κ 11 , κ 12 (resp. κ 21 , κ 22 ). In general, increasing the inter-or intra-interaction strength will decrease the probability of having a co-feasible equilibrium.

If κ 12 = κ 21 = ∞, then condition (23) gives the co-feasibility conditions for each community:

s 11 < 1 2β 1 log(n) and s 22 < 1 2β 2 log(n) .
For the same s, it means s < 1 2 log(n) max(β 1 , β 2 ) .

As an example of application, suppose we start with co-feasible communities of equal size (β 1 = β 2 = 0.5) and add interactions between these two groups, co-feasibility may be dropped (see Fig. 2). The co-feasibility domain is illustrated in Fig. 7. It shows a threshold where the co-feasibility property is satisfied above the curve. This means that the lower the values of κ 11 and κ 22 , i.e. the stronger the interactions within the groups, the more likely the co-feasibility property is lost. We can conclude that an independent group structure is more likely to be co-feasible and therefore stable, which supports previous work on compartmentalization models [START_REF] Stouffer | Compartmentalization increases food-web persistence[END_REF].

& R P P X Q L W \ ( 1 1 ) & R P P X Q LW \ ( 2 2)
,QWHU FRP PXQ LWLHV ( The z-axis (resp x-axis) is the strength of interaction within community 1 -κ 11 (resp community 2 -κ 22 ). The y-axis is the inter-community interactions κ 12 = κ 21 . The colored area illustrates the threshold between the co-feasible and non-co-feasible domains in the system (3).

Impact of the community size

For a fixed matrix κ, the condition to have a co-feasible fixed point can be computed as a function of the size of each community i.e. the pair β = (β 1 , β 2 ). Starting from the co-feasibility inequality (23):

max 2β 1 κ 2 11 + 2β 2 κ 2 12 , 2β 1 κ 2 21 + 2β 2 κ 2 22 < 1 ,
the two components are studied independently,

2β 1 κ 2 11 + 2(1 -β 1 ) κ 2 12 < 1 ⇔ β 1 2 κ 2 11 - 2 κ 2 12 < 1 - 2 κ 2 12 ⇒ β 1 < 1 -2 κ 2 12 2 κ 2 11 -2 κ 2 12 if κ 11 < κ 12 .
Similarly, one has

2β 1 κ 2 21 + 2(1 -β 1 ) κ 2 22 < 1 ⇒ β 1 > 1 -2 κ 2 22 2 κ 2 21 -2 κ 2 22 if κ 22 < κ 21 .
In the case where the intra-community interactions (κ 11 , κ 22 ) are smaller than the inter-community interactions (κ 12 , κ 21 ), we obtain an upper and a lower bound for the admissible size of each community β 1 , β 2 to have a co-feasible equilibrium. In Fig. 8, different cases of the co-feasibility zone are represented according to the inter-community interactions (κ 12 , κ 21 ). If the intra-community interactions are different, the community with the lowest interaction κ ii is advantaged i.e. the size of the community can be larger.

Connection increases co-feasibility

In Section 4.2, we analyzed the co-feasibility condition for a scenario involving two communities. We presented a co-feasibility domain defined by

(s n • s n )β ⊤ ∞ < 1 2 log(n) := (s * n ) 2 , ⇔ max β 1 s 2 11 + β 2 s 2 12 , β 1 s 2 21 + β 2 s 2 22 < 1 2 log(n) , ⇔ β 1 s 2 11 + β 2 s 2 12 < 1 2 log(n) and β 1 s 2 21 + β 2 s 2 22 < 1 2 log(n) . ( 24
)
These two distinct conditions within the communities have led us to the conclusion that community isolation is beneficial for coexistence. However, general constraints that affect all interactions could be considered at the ecosystem scale. To this end, we introduce a complementary condition: the global variance of the interaction coefficients in the ecosystem remains invariant, i.e.

β 1 s 2 11 + β 2 s 2 12 + β 1 s 2 21 + β 2 s 2 22 = Γ , with Γ > 0.
In order to remove the n dependency, we define γ 2 ij = s 2 ij log(n) and Γ n = Γ log(n) (note that the γ coefficients are defined differently from the κ coefficients). Combined with condition (24), we get the following system of equations:

     β 1 γ 2 11 + β 2 γ 2 12 + β 1 γ 2 21 + β 2 γ 2 22 = Γ n , 2β 1 γ 2 11 + 2β 2 γ 2
12 < 1, (feasibility condition for community 1) , 2β 1 γ 2 21 + 2β 2 γ 2 22 < 1, (feasibility condition for community 2) .

(

) 25
Assuming the fixed intra-community variances γ 11 and γ 22 , we seek to determine the co-feasibility conditions for the inter-community variances γ 12 and γ 21 . In this case, the constraint on the total variance corresponds to the equation of an ellipse in the (γ 12 , γ 21 ) plane:

β 1 γ 2 11 + β 2 γ 2 12 + β 1 γ 2 21 + β 2 γ 2 22 = Γ n , ⇔ β 2 Γ n -β 1 γ 2 11 -β 2 γ 2 22 γ 2 12 + β 1 Γ n -β 1 γ 2 11 -β 2 γ 2 22 γ 2 21 = 1 ,
The equations below provide the values for the semi-major axis a and the semi-minor axis b of the ellipse:

a = Γ n -β 1 γ 2 11 -β 2 γ 2 22 β 2 , b = Γ n -β 1 γ 2 11 -β 2 γ 2 22 β 1 .
Note that if both communities are of equal size (β 1 = β 2 ), a circle with radius a is obtained.

From a visual standpoint, the conditions (25) are depicted in Figure 9. Since the coefficients (γ 12 , γ 21 ) are non-negative, we are only interested in the positive orthant. The feasibility condition for community 1 is given by the horizontal axis defined by

2β 1 γ 2 11 + 2β 2 γ 2 12 < 1 ⇔ γ 12 < 1 -2β 1 γ 2 11 2β 2 ,
and the one of community 2 is given by the vertical axis defined by

2β 1 γ 2 21 + 2β 2 γ 2 22 < 1 ⇔ γ 21 < 1 -2β 2 γ 2 22 2β 1 .
The intersection between the vertical (resp. horizontal) line and the ellipse occurs when the semi-major axis (resp. semi-minor axis) exceeds the vertical condition a >

(1

-2β 1 γ 2 11 )/(2β 2 ) (resp. horizontal condition b > (1 -2β 2 γ 2 22 )/(2β 1 )
). Remark 6. By replacing the feasibility condition of γ 21 in the equation of the ellipse, we can derive the intersection between the vertical axis and the ellipse as follows:

γ 2 21 = Γ n -β 2 γ 2 22 -1/2 β 1 ,
equivalent to the feasibility condition of γ 12 (by replacing Γ n ):

γ 12 < 1 -2β 1 γ 2 11 2β 2 .
We identify the range of co-feasibility between the two groups for γ 21

Γ n -β 2 γ 2 22 -1/2 β 1 < γ 2 21 < 1 -2β 2 γ 2 22 2β 1 .
This simple framework allows for testing different scenarios. Figure 9a is the reference figure. It represents a situation where all the potential interactions between communities lead to co-feasibility. In Figure 9b, the total system variance is increased, which results in a reduction of co-feasibility options: the interactions between the two communities must be high γ 12 , γ 21 ≫ 0. When the intra-community variances are increased as shown in Figure 9c, the ellipse shrinks, and the set of inter-community interaction variances is reduced. This observation reinforces the findings of Section 4.2 where large interaction between communities enables co-feasibility through community isolation. In the concluding example in Figure 9d, we reduce only the γ 11 interaction in community 1. We observe that the impact of community 1 on 2 γ 21 must be weaker, but the impact of community 2 on 1 γ 12 can no longer be weak. Weaker interactions within community 1 imply a stronger connection between the communities for co-feasibility.

Discussion

In this paper, we described a model of the dynamics of species abundances when the interaction among species is structured in multiple communities. The main interest is to outline the effect of a block structure on the stability and persistence of species. We defined an interaction matrix per block which has several characteristics such as the strength of the interactions s and the size of the community β. Specifically, we described the dynamics and properties of each community in the system (feasibility, proportion of surviving species, mean and root mean square of the abundances of surviving species) and their effect on each other.

In this context, we focused most of our analysis on the case of two interacting communities. However, our results can be extended to more than 2 communities (see Appendix E). First, theoretical conditions were given for a unique globally stable equilibrium in the model (3) with surviving and vanishing species. This follows from Lyapunov conditions related to a result of Takeuchi and Adachi [START_REF] Takeuchi | The existence of globally stable equilibria of ecosystems of the generalized Volterra type[END_REF] and random matrix theory. These stability results had been found in the case of a single community by Clenet et al. [START_REF] Clenet | Equilibrium and surviving species in a large Lotka-Volterra system of differential equations[END_REF]. This complements the stability properties in the Lotka-Volterra system studied by Stone [START_REF] Stone | The feasibility and stability of large complex biological networks: a random matrix approach[END_REF] and Gibbs et al. [START_REF] Gibbs | Effect of population abundances on the stability of large random ecosystems[END_REF]. Recent random matrix methods allow us to describe the spectrum of a block matrix and plot it numerically. For a detailed discussion of random matrices in the Lotka-Volterra model, see Akjouj et al. [ABC + 22].

Subsequently, we gave heuristics on the surviving species (proportion, mean and root mean square of their abundances). These heuristics have also been found in the case of a single community by Clenet et al. [START_REF] Clenet | Equilibrium and surviving species in a large Lotka-Volterra system of differential equations[END_REF]. From a physicist's point of view and using the methods of Bunin [START_REF] Bunin | Ecological communities with Lotka-Volterra dynamics[END_REF] and Galla [START_REF] Galla | Dynamically evolved community size and stability of random Lotka-Volterra ecosystems[END_REF] where they consider a different growth rate for each species. The study of the stability and properties of surviving species in the LV system has also been carried out by Pettersson et al. [START_REF] Pettersson | Predicting collapse of complex ecological systems: quantifying the stability-complexity continuum[END_REF][START_REF] Pettersson | Stability of ecosystems enhanced by speciesinteraction constraints[END_REF]. From an ecological point of view, heuristics are derived from the properties of interactions between multiple communities.

In a third part, we studied the condition under which the feasibility threshold exists where all species coexist. We extend the feasibility results found by Bizeul and Najim [START_REF] Bizeul | Positive solutions for large random linear systems[END_REF] in the case of a block structure. A co-feasibility threshold was found in the form of an inequality that must be verified to have a feasible community set. This complements the recent results on interactions with a sparse structure [START_REF] Akjouj | Feasibility of sparse large Lotka-Volterra ecosystems[END_REF] and interactions with a correlation profile [START_REF] Clenet | Equilibrium in a large Lotka-Volterra system with pairwise correlated interactions[END_REF]. We notice that to maximize the probability of co-feasibility, we need to minimize the interactions between the communities. Additionally, a community with weaker interactions can exhibit a larger total abundance in the ecosystem while maintaining the co-feasibility threshold. At the ecosystem level, when a generic constraint that affects all interactions is added, weaker interactions within one of the communities suggest a stronger connection between the communities for co-feasibility.

There are still many mathematical and ecological questions that remain unanswered in this type of model. First, a rigorous mathematical proof of the heuristics presented here would be of interest, although the LCP procedure induces an a priori statistical bias that is difficult to handle. This issue is still pending in the single community case [START_REF] Clenet | Equilibrium and surviving species in a large Lotka-Volterra system of differential equations[END_REF] and appears to be challenging to address. Recently, Akjouj et al. [START_REF] Akjouj | Equilibria of large random Lotka-Volterra systems with vanishing species: a mathematical approach[END_REF] provided a rigorous proof using an approximate message passing (AMP) approach in the single community model with an interaction matrix taken from the Gaussian Orthogonal Ensemble (GOE). Their approach was based on work by Hachem [START_REF] Hachem | Approximate Message Passing for sparse matrices with application to the equilibria of large ecological Lotka-Volterra systems[END_REF].

Second, we could extend the heuristics for two different scenarios. On the one hand, it would be interesting to add pairwise correlations between species coefficients A ij . This has already been done by physicists, see [START_REF] Barbier | Generic assembly patterns in complex ecological communities[END_REF][START_REF] Poley | Generalized Lotka-Volterra model with hierarchical interactions[END_REF]. In the study of feasibility, it was shown that a correlation profile does not change the feasibility threshold [START_REF] Clenet | Equilibrium in a large Lotka-Volterra system with pairwise correlated interactions[END_REF]. On the other hand, for the sake of simplicity, we have chosen to set the growth rates equal to the same value r k = 1 for k ∈ [n]. It would be relevant to control the distribution of the growth rate as in [SCG + 18] or to consider structural stability as in Saavedra et al. [SRB + 17], i.e. how much can the growth rates be perturbed (initially all equal to 1) without changing the type of equilibrium x * obtained.

There are many applications of this kind of models in ecology. We could consider a spatial structure that accounts for spatial proximity in the sense that two nearby communities tend to be more strongly connected. For example, in an aquatic environment, we could imagine the existence of an up/down gradient in a water column. Fig. 10 illustrates a situation where three communities are involved. Originally introduced by R.T. Paine [START_REF] Paine | Food Web Complexity and Species Diversity[END_REF][START_REF] Paine | The Pisaster-Tegula Interaction: Prey Patches, Predator Food Preference, and Intertidal Community Structure[END_REF], the concept of keystone species is widely used in ecology i.e. one species controls the coexistence of the others and species are lost after the removal of this keystone species. Mouquet et al. [START_REF] Mouquet | Extending the concept of keystone species to communities and ecosystems[END_REF] suggested extending the concept of keystone species to communities. In the block system, one could analyze the existence of a keystone community that would have disproportionately large effect on other communities. In a metacommunity dynamic, Resetarits et al. [START_REF] Resetarits | Testing the keystone community concept: effects of landscape, patch removal, and environment on metacommunity structure[END_REF] have explored the concept of keystone communities, where some patches have stronger effect on others.

One could imagine that the same species is present several times in the system, but in different blocks, see Gravel et al. [START_REF] Gravel | Stability and complexity in model meta-ecosystems[END_REF]. In this case, the inter-blocks represent interactions between spatially isolated communities (so should be less strong). If each diagonal or non-diagonal block is a copy of the same interaction pattern (possibly slightly perturbed) and we can add linear effects to the system to represent emigration and immigration, then we could study the feasibility properties of this system. In [START_REF] Gravel | Stability and complexity in model meta-ecosystems[END_REF], they found that stability is most likely when dispersal (which controls off-diagonal blocks) is intermediate.

Last but not least, it would be relevant to compare the patterns obtained with data in ecology, as in the recent article by Hu et al. [HAB + 22] in the case of a single community.

Starting from the QVE (11) associated to the matrix H, the system takes the simpler form

-1 µ(z) = z + 2β 1 s 2 11 µ(z) + β 2 (s 2 12 + s 2 21 )ν(z) -1 ν(z) = z + β 1 (s 2 12 + s 2 21 )µ(z) + 2β 2 s 2 22 ν(z)
, where ∀k ∈ I 1 , m k (z) = µ(z) and ∀k ∈ I 2 , m k (z) = ν(z). All the knowledge of equation ( 11) relies on the functions µ(z) and ν(z). Then, using RMT theory, the resolvent G of the symmetric matrix H can be approximated by

G(z) = (H -zI) -1 ≃ diag(µ(z)1 ⊤ I1 , ν(z)1 ⊤ I2
) . From Remark 7, the trace of the resolvent is equal to the Stieltjes transform

g(z) = 1 n Tr(G) ≃ β 1 µ(z) + β 2 ν(z)
of the spectral measure. Finally, the spectral density can be obtained using a the Stieltjes inversion (Prop. 5). The spectral density of the matrix H can be computed numerically by an iterative scheme. The initial condition of the two measurements (µ, ν) is µ 0 = ν 0 = -1 z . Then, the iterative scheme Remark 8. To handle the Stieltjes inversion (Prop.5) numerically, it is similar as starting with z = x+ϵi, ϵ ≈ 10 -3 .

-1 µp = z + 2β 1 s 2 11 µ p-1 + β 2 (s 2 12 + s 2 21 )ν p-1 -1 νp = z + β 1 (s 2 12 + s 2 21 )µ p-1 + 2β 2 s 2 22 ν p-1 , converge to µ ∞ = lim
In Fig. 11, we present the numerical estimation of the spectral density for different types of interactions of the matrix H.

C Remaining computations C.1 Moments of Žk

We compute hereafter the conditional variance of Žk = (Bx * ) k with respect to x * . We rely on the following identities ∀k ∈ I i , ∀ℓ ∈ I j and ∀o ∈ I q :

EB kℓ = 0 , E(B 2 kℓ ) = s 2 ij n , EB kℓ B ko = 0 (ℓ ̸ = o) .
We first compute the conditional mean:

∀k ∈ I i , E x * ( Žk ) = ℓ∈[n] E(B kℓ )x * ℓ = ℓ∈S1 E(B kℓ )x * ℓ + ℓ∈S2 E(B kℓ )x * ℓ = 0 .
We now compute the second moment:

∀k ∈ I i , E x * ( Ž2 k ) = E x *   ℓ∈[n] B kℓ x * ℓ   2 = E x * ℓ,o∈[n] B kℓ B ko x * ℓ x * o , = ℓ∈S1∪S2 E(B 2 kℓ )x * 2 ℓ + ℓ̸ =o E(B kℓ B ko )x * ℓ x * o , = s 2 i1 n ℓ∈S1 x * 2 ℓ + s 2 i2 n ℓ∈S2 x * 2 ℓ , (a)
≃ where the approximation in (a) follows from the fact that

β 1 p1 σ2 1 s 2 i1 + β 2 p2 σ2 2 s 2 i2 , (a) β = [1/2, 1/2], s = 1/ √ 2 1/ √ 2 1/ √ 2 1/ √ 2 (b) β = [1/2, 1/2], s = 1/2 1 1 1/8 (c) β = [1/2, 1/2], s = 1 1/2 1/8 1 (d) β = [3/4, 1/4], s = 1/3 1/5 1 1/2
1 n 2 ℓ,o∈[n] x * ℓ x * o = 1 n 2 ℓ̸ =o x * ℓ x * o + O 1 n .
We can now compute the variance:

∀ k ∈ I i , Var x * Žk = E x * Ž2 k -E x * Žk 2 = β 1 p1 σ2 1 s 2 i1 + β 2 p2 σ2 2 s 2 i2 .

C.2 Details of heuristics of the mean

Our starting point is the following generic representation of an abundance at equilibrium (either of a surviving or vanishing species) in the case k ∈ S i :

x

* k = (1 + ∆ * i Z k ) 1 {Z k >δ * i } = 1 {Z k >δ * i } + (∆ * i Z k ) 1 {Z k >δ * i } .
Summing over S i and normalizing,

1 |S i | k∈Si x * k = 1 |S i | k∈Si 1 {Z k >δ * i } + ∆ * i 1 |S i | k∈Si Z k 1 {Z k >δ * i } , mi (a) = 1 + ∆ * i |I i | |S i | 1 |I i | k∈Ii Z k 1 {Z k >δ * i } , mi (b) ≃ 1 + ∆ * i 1 P(Z > δ * i ) E(Z1 {Z>δ * i } ), mi ≃ 1 + ∆ * i E(Z | Z > δ * i )
. where (a) follows from the fact that

|S i | = k∈Si 1 {Z k >δ * i (by definition of S i ), (b) from the law of large numbers 1 |Ii| k∈Ii Z k 1 {Z k >δi} ----→ n→∞ EZ1 {Z>δ * i } and |Si| |Ii| ----→ n→∞ P(Z > δ * i )
with Z ∼ N (0, 1). It remains to replace mi by its limit m * i to obtain the heuristics of the mean:

m * 1 = 1 + ∆ * 1 E(Z | Z > δ * 1 ) , m * 2 = 1 + ∆ * 2 E(Z | Z > δ * 2 ) .
C.3 Density of the distribution of the surviving species.

Assume that x * > 0, and let f = R → R be a bounded continuous test function, then ∀k ∈ S

i Ef (x * k ) = E f (1 + λ * i + ∆ * i Z k ) Z k > δ * i , = ∞ -∞ f (1 + λ * i + ∆ * i u) 1 {u>δ * i } 1 -Φ(δ * i ) e -u 2 2 √ 2π du , = ∞ 0 f (y)e -1 2 y ∆ * i +δ * i 2 1 √ 2πΦ(-δ * i )∆ * i dy ,
hence the density of x * k , ∀k ∈ S i .

D Sketch of proof of Theorem 4

The first step consists in decomposing the equilibrium x * :

x * k = e ⊤ k x * = e ⊤ k (I -B) -1 1 = ∞ ℓ=0 e ⊤ k B ℓ 1 = 1 + e ⊤ k B1 + e ⊤ k B 2 (I -B) -1 1 , = 1 + Z k + R k , where Z k = n ℓ=1 B kℓ , ∀k ∈ [n]. One can prove that ∀k ∈ [n], R k is a negligible term if n
is sufficiently large. From a technical point, it relies on Gaussian concentration of Lipschitz functionnals and we are confident that the techniques applied in [START_REF] Bizeul | Positive solutions for large random linear systems[END_REF] will succeed in handling R k . However, this part of the proof is not been treated here since we want to stick a concise argumentation of the proof which gives the reader information about the critical bound of the feasibility threshold.

The feasibility of the two communities is studied independently. Using Gaussian addition properties, a simpler form of Z k is first deduced. Consider a family ( Žk ) k∈[n] of i.i.d. random variables N (0, 1).

If k ∈ I 1 , Z k = ℓ∈I1 B kℓ + ℓ∈I2 B kℓ , ∼ N 0 , β 1 s 2 11 + N 0 , β 2 s 2 12 , ∼ β 1 s 2 11 + β 2 s 2 12 Žk . Similarly if k ∈ I 2 , Z k ∼ β 1 s 2 21 + β 2 s 2 22
Žk . Given β = (β 1 , β 2 ), conditions on the matrix s are inferred to have:

P min k∈[n] x k > 0 = 1 ⇔ P min k∈[n] Z k > -1 = 1 . (26)
In order to compute a tractable form of min

k∈[n] Z k , an additional approximation is made, if n is large enough ∀i ∈ {1, 2}, min k∈Ii Žk ∼ -2 log(β i n) ≈ -2 log(n) (27) min k∈[n] Z k = min β 1 s 2 11 + β 2 2 12 min k∈I1 Žk , β 1 s 2 21 + β 2 s 2 22 min k∈I2 Žk , (a)
≃ min

β 1 s 2 11 + β 2 s 2 12 -2 log(n) , β 1 s 2 21 + β 2 s 2 22 -2 log(n) , = min -2β 1 s 2 11 log(n) + 2β 2 s 2 12 log(n), -2β 1 s 2 21 log(n) + 2β 2 s 2 22 log(n) , = -max 2β 1 s 2 11 log(n) + 2β 2 s 2 12 log(n), 2β 1 s 2 21 log(n) + 2β 2 s 2 22 log(n) .
where (a) comes from the approximation (27). The condition min

k∈[n] Z k > -1 asymptotically boils down to max 2β 1 s 2 11 log(n) + 2β 2 s 2 12 log(n), 2β 1 s 2 21 log(n) + 2β 2 s 2 22 log(n) < 1 , ⇔ max 2β 1 s 2 11 log(n) + 2β 2 s 2 12 log(n), 2β 1 s 2 21 log(n) + 2β 2 s 2 22 log(n) < 1 , ⇔ max β 1 s 2 11 + β 2 s 2 12 , β 1 s 2 21 + β 2 s 2 22 < 1 2 log(n) , ⇔ (s n • s n )β ⊤ ∞ < 1 2 log(n) := (s * n ) 2 .

E Extension to the b-blocks model E.1 Interaction matrix with b communities

Within the framework of b communities, the matrix B = (B kℓ ) n,n is defined as

B = V sV ⊤ • 1 √ n A , (28)
where

V ∈ M n×b , V =      1 I1 0 • • • 0 0 1 I2 • • • 0 . . . . . . . . . . . .
Following the same arguments as in Theorem 3 and relying on Theorem 2, given the particular shape of the matrix S, computing its norm is equivalent to computing the norm of a matrix of size b ----→ n→∞ σ * i .

∥S∥ ∞ = diag(β) (s • s) + (s • s) ⊤ ∞ .

E.3 Surviving species

E.4 Distribution of the surviving species

Let s be the b × b matrix of interaction strengths and assume that the condition of Theorem 3 holds. Let x * the solution of (8) and (p * , σ * ) the solution of the heuristic 3. Recall the definition of ∆ i , δ i and denote by δ * i = δ i (p * i , σ * i ). Let x * k > 0 a positive component of x * belonging to the community i, then:

L(x * k ) ----→ n→∞ L 1 + ∆ * i Z Z > δ * i ,
where Z ∼ N (0, 1). Otherwise stated, asymptotically ∀k ∈ S i , x * k admits the following density

f k (y) = 1 {y>0} Φ(-δ * i ) 1 ∆ * i √ 2π exp - 1 2 y ∆ * i + δ * i 2 . (30)

E.5 Feasibility

We consider a growing scaling matrix s n ----→ n→∞ 0 ⇔ ∀i, j ∈ {1, b}, s ij ----→ n→∞ 0 .

Let B n a matrix defined by

B n = V s n V ⊤ • 1 √ n A . (31)
The spectral radius of 1 √ n A a.s. converges to 1 (circular law). So as long as s n is close to zero, the matrix I -B n is eventually invertible.

Recall the problem which admits a unique solution defined by x k > 0 ----→ n→∞ 0 .

x * = 1 + Bx * ⇔ x * = (I -B) -1 1 , (32)
2. If there exists ε > 0 such that eventually

(s n • s n )β ⊤ ∞ ≤ (1 -ε)(s * n ) 2 then P min k∈[n]
x k > 0 ----→ n→∞ 1 .

Sketch of proof.

Starting from the decomposition, the equilibrium x * :

x * k = 1 + Z k + R k ,
where Z k = n ℓ=1 B kℓ , ∀ k ∈ [n] and we assume that ∀ k ∈ [n] , R k is a negligible term if n is sufficiently large.

The feasibility of the b communities is studied independently. Using Gaussian addition properties, a simpler form of Z k is derived. Consider a family ( Žk ) k∈[n] of i.i.d. random variables N (0, 1). x k > 0) = 1 ⇔ P( min

k∈[n]
Z k > -1) = 1 .

In order to compute a tractable form of min

⇔ (s n • s n ) 2 β ⊤ ∞ < 1 2 log(n) := (s * n ) 2 .

B

  kℓ x ℓ (t)

Figure 1 :

 1 Figure 1: Dynamics of the model (1) with 2 distinct communities of 5 species each and interaction matrix given by (2). The two communities converge to their feasible equilibrium point and do not interact. In Fig. (a), a heat map illustrates the interaction matrix (2).Figure (b) shows the dynamics where the species of each community reach a feasible equilibrium.

  Abundance dynamics of two communities of five species

Figure 2 :

 2 Figure 2: Dynamics of Model (1) with 2 distinct communities of 5 species each with interaction matrix given by (2). In Fig. (a), representation of the interaction matrix (2) when the interactions between the communities are strong. Fig. (b) shows the dynamics where the species of each community reach an equilibrium. Notice that there are vanishing species in each community.

)Following

  Theorem 2.1 in Ajanki et al. [AEK19], ∀z ∈ C + , Equation (11) has a unique solution m = m(z) where n -1 m i is the Stieltjes transform of a probability measure and the support of the associated measure is included in [-Σ, Σ], where Σ = 2 ∥S∥ 1/2

  admits a unique solution (p * 1 , p * 2 , σ * 1 , σ * 2 ) and for i ∈ {1, 2} pi a.s.

  (a) Parameters (p * 1 , σ * 1 ) versus (s21, s12). (b) Parameters (p * 2 , σ * 2 ) versus (s21, s12).

Figure 3 :

 3 Figure 3: Comparison between the theoretical solutions (p * 1 , p * 2 , σ * 1 , σ * 2 ) of (12)-(15) and their empirical Monte Carlo counterpart (the star marker) as functions of the off-diagonal block interaction strength (s 12 , s 21 ). The left column is associated to the properties of community 1. The right column is associated to the properties of community 2. Matrix B has size n = 500 and the number of Monte Carlo experiments is 500. The parameters are s = 1/3s 12 s 21 1/ √ 2 , β = 1 2 , 1 2 .When off-diagonal block interactions s 12 , s 21 increase, the proportion of surviving species p * decreases but the root mean square of their equilibrium abundances σ * increases.

Figure 4 :

 4 Figure 4: Distribution of surviving species in each community. The x-axis represents the value of the abundances and the histogram is built upon the positive components of equilibrium x * associated to each community. The blue-solid line (resp. red-solid line) represents the theoretical distribution of community 1 (resp. community 2) for parameters s as given by Heuristics 2. The blue-dashed vertical line (resp. reddashed vertical line) corresponds to the mean abundance of community 1 (resp. community 2) given by equations (18)-(19). All th ematrix entries A ij 's are independent Gaussian N (0, 1); the parameters are set to n = 2000, β = (0.75, 0.25) and s = 1/2 1/ √ 2 1/5 1/9 .

Figure 5 :

 5 Figure 5: Distribution of surviving species in each community with non-Gaussian entries. The x-axis represents the value of the abundances. The histogram is built upon the positive components of equilibrium x * associated to each community. (blue for community 1, red for community 2). The solid lines represent the theoretical distributions associated to parameters s as given by Heuristics 2. The entries of the A ij matrices are uniform U(-√ 3, √ 3) with variance 1 and the parameters are set to n = 2000, β = (0.5, 0.5)

  1 represents the proportion of each community. Let s n ----→ n→∞ 0 and denote by s * n = 1/ √ 2 log n the critical threshold. Let x n = (x k ) k∈[n] be the solution of (21).

Figure 6 :

 6 Figure 6: Transition towards co-feasibility for the 2-blocks model (2). For each value ν 1 on the x-axis, we simulate 500 matrices B of size n = 5000 with two communities of the same size (β 1 = β 2 = 0.5) with the inter-block interactions fixed at s 21 (ν 2 ) = s 21 (ν 2 ) = 1/2 log(n) and compute the solution x of Theorem 4 at the scaling for the intra-block interactions s 11 (ν 1 ) = s 22 (ν 1 ) = 1/ν 1 log(n). The curve represents the proportion of feasible solutions x obtained for the 500 simulations. The dotdashed vertical line corresponds to ν 1 = β1

Figure 7 :

 7 Figure 7: Representation of the co-feasibility phase diagram. The co-feasible domain is above the surface.The z-axis (resp x-axis) is the strength of interaction within community 1 -κ 11 (resp community 2 -κ 22 ). The y-axis is the inter-community interactions κ 12 = κ 21 . The colored area illustrates the threshold between the co-feasible and non-co-feasible domains in the system (3).

Figure 8 :

 8 Figure 8: Representation of the co-feasibility domain depending on the fixed intra-community interaction.In (a), (b), (c), a different scenario of intra-community interaction is presented. Each panel represents the upper-bound (blue curve) and the lower-bound (red curve) of the size of community 1 as a function of the interaction between two communities (κ 12 , κ 21 ). The blue area is the admissible zone to have a co-feasible fixed point in (3). The size of community 2 is equal to β 2 = 1 -β 1 .

  (a) γ11 = γ22 = 0.5, Γn = 0.6. (b) γ11 = γ22 = 0.5, Γn = 0.8. (c) γ11 = γ22 = 0.7, Γn = 0.6. (d) γ11 = 0.2, γ22 = 0.5, Γn = 0.6.

Figure 9 :

 9 Figure 9: Graphical representation of equation system (25). The solid line circle constrains the total variance of the system (equation 1 in (25)), while the dashed and dotted lines correspond to feasibility conditions for community 1 (equation 2 in (25)) and 2 (equation 3 in (25)), respectively. In each figure, two communities of equal size are considered with β 1 = β 2 = 1/2, creating the solid line circle (rather than an ellipse) in this particular case. Each figure illustrates a distinct interaction situation between two communities outlined in the caption. Figure (a) serves as the reference against which the other figures (b)-(c)-(d) are compared. The co-feasibility arises only for the values of (γ 12 , γ 21 ) on the line found within the inner square at the bottom left-hand corner of the figure. When values of (γ 12 , γ 21 ) on the line are outside this inner square, then one or the other community is not feasible anymore.

Figure 10 :

 10 Figure10: In (a), a representation of the gradient of interaction between three communities in a water column is represented. The blue arrows correspond to strong interaction strength due to their spatial proximity. On the opposite, the communities 1 and 3 are separated, the green arrow represents a weaker interaction. In (b), the block matrix associated with this type of model is displayed. The colors of the blocks corresponds to the colors of the arrows. The red colored block corresponds to intra-community interactions.

  p→+∞ µ p and ν ∞ = lim p→+∞ ν p . The last step consist of using the property of the Stieltjes inversion (Prop.5).

Figure 11 :

 11 Figure 11: Spectrum (histogram) of the Hermitian random matrix H (n = 1000), conditions on (β, s) are given in each panel. The numerical approach is used to compute the solid line spectrum distribution. An upper bound for the largest eigenvalue of H, given by 2 ∥S∥ 1/2 ∞ , is denoted by the dashed vertical line.

  I 1 = [n 1 ], I 2 = {n 1 + 1, • • • , n 1 + n 2 },. . . , I b = {n 1 + • • • + n b-1 + 1, • • • , n} the subset of [n] of size |I i | := n imatching the index of species belonging to community i and β = (β 1 , β 2 , .., β b ) with ∀i ∈ [b], β i = n i /n and b i=1 β i = 1. The random matrix A ij is non-Hermitian of size |I i | × |I j | with standard Gaussian entries i.e. N (0, 1). Recall that 1 Ii be the element-wise vector of 1 with size |I i |.

E. 2

 2 Existence of a unique equilibriumLet H be the symmetric matrixH = B + B ⊤ ∀ i, j ∈ [b], H ij is a matrix of size |I i | × |I j | and each off-diagonal entries follow a Gaussian distribution N (0, s 2 ij + s 2 ji ).The QVE associated to the matrix H is decomposed ask ∈ I i , + s 2 ji m ℓ . Given m(z) = (m 1 (z), • • • , m n (z)), denote by 1/m(z) = (1/m 1 (z), • • • , 1/m n (z)) and S = 1 n V (s + s ⊤ )V ⊤ the QVE can be written in the standard form -1 m(z)= z + Sm(z) .

Heuristics 3 .

 3 Let s be the b × b matrix of interaction strengths and assume that the condition of Theorem 3 holds, then the following system of 2b equations and 2b unknowns p = (p 1 , p 2 , .., p b ),σ = (σ 1 , σ 2 , .., σ b ) ∀i ∈ [b] , p i = 1 -Φ(δ i ) , ∀i ∈ [b] , (σ i ) 2 = 1 + 2∆ i E(Z|Z > δ i ) + ∆ 2 i E(Z 2 |Z > δ i ) ,where∆ i = b j=1 p j (σ j ) 2 β j s 2 ij ; δ i = -1 ∆ i ,admits a unique solution (p * , σ * ) and ∀ i ∈ [b]

Theorem 6 (

 6 Co-feasibility for the b-blocks model). Assume that matrix B n is defined by the b-blocks model (31). Let β = (β 1 , β 2 , .., β b ), b i=1 β i = 1 represents the proportion of each community. Let s n ----→ n→∞ 0 and denote by s* n = 1/ √ 2 log n. Let x n = (x k ) k∈[n]be the solution of (32).1. If there exists ε > 0 such that eventually(s n • s n )β ⊤ ∞ ≥ (1 + ε)(s * n ) 2 then P min k∈[n]

  If k ∈ I i , Z k = Given β = (β 1 , β 2 , .., β b ), conditions on the matrix s are inferred to have P( min k∈[n]

  k∈[n]Z k , an additional approximation is made, if n is large enoughmin k∈Ii Žk ∼ -2 log(β i n) ≃ -2 log(n) .Following the approximation (33), the condition mink∈[n]Z k > -1 asymptotically boils down to max

  , Barbier et al. [BABL18] and Poley et al. [PBG23] have extended the heuristics in the block and cascade model. Previously, obtaining properties on surviving species in the LV model (not normalized by √ n) was already done by Servan et al. [SCG + 18]

The simplifying assumption r k = 1 allows tractable computations and could be extended to r k = c with c > 0. However, if the growth rate is different for each species, the mathematical development and result may be strongly affected and will be discussed in each section.

* Supported by CNRS Project 80 Prime -KARATE.

A Stieljes transform

We provide some reminders regarding Stieltjes transforms, a central element of proofs in random matrix theory. We denote by C + := {z ∈ C : Im(z) > 0} the upper half of the complex plane.

Definition 2 (Stieltjes transform). Let ν ∈ P(R) be a probability measure. The Stieltjes transform of ν, denoted by g ν : C + → C, is defined by

Remark 7. Let ν X be the empirical measure of the eigenvalues λ 1 (X), • • • , λ n (X) of the symmetric matrix X ∈ M n (C) define by

then the associated Stieltjes transform is given by

where Q = (X -zI) -1 is the resolvent of the matrix X and Tr(Q) is the trace of matrix Q.

B Numerical methods

All figures and code are available on Github [START_REF] Clenet | Impact of a block structure on the lotka-volterra model[END_REF]. The code is written in Python. To verify the system of equations of heuristics 1, the simulations on the properties of surviving species are performed in two distinct methods (see Fig. 3). On the one hand, we use a standard solver (cf. scipy.optimize) to find the theoretical solutions by finding a local minimum of the system of equations (a modification of the Powell hybrid method). On the other hand, we simulate a large number of matrix B, each corresponding to an experiment, and we resolve the associated LCP problem using the Lemke's algorithm (see the lemkelcp package [START_REF] Lamperski | Lemke's algorithm for linear complementarity problems[END_REF]). The empirical solutions are computed using a Monte Carlo experiment, i.e. we use the LCP solution to compute the properties of the surviving species and we make an average over the ensemble of experiments. As a baseline, the dynamics of Lotka-Volterra are achieved by a Runge-Kutta method of order 4 (RK4) implemented in the code. Spectrum: a computer based approach. Theorem 3 only provides sufficient conditions for the existence of a unique stable equilibrium and is based on the rough asymptotic upper bound estimation Σ = 2 ∥S∥ 1/2 ∞ . We can assess the sharpness of this bound by comparing it to the limiting spectrum of matrix H, which can be plotted via numerical simulations. An efficient way to compute numerically the spectrum of the matrix H comes from the system of non linear equations (11).