Maxime Clenet

François Massol

Jamal Najim

Equilibrium and surviving species in a large Lotka-Volterra system of differential equations

Lotka-Volterra (LV) equations play a key role in the mathematical modeling of various ecological, biological and chemical systems. When the number of species (or, depending on the viewpoint, chemical components) becomes large, basic but fundamental questions such as computing the number of surviving species still lack theoretical answers. In this paper, we consider a large system of LV equations where the interactions between the various species are a realization of a random matrix. We provide conditions to have a unique equilibrium and present a heuristics to compute the number of surviving species. This heuristics combines arguments from Random Matrix Theory, mathematical optimization (LCP), and standard extreme value theory. Numerical simulations, together with an empirical study where the strength of interactions evolves with time, illustrate the accuracy and scope of the results.

Introduction

Since May's seminal work [START_REF] May | Will a large complex system be stable?[END_REF] and for the past decades, many theoretical studies addressed the issue of the coexistence of species in ecosystems.

Introduced in the 1920s by Lotka [Lot25] and Volterra [Vol26], the Lotka-Volterra (LV) model is a well-known classic in theoretical ecology and mathematics. It represents a first step in our understanding of ecosystems through the variety of its dynamical behaviours (single or multiple equilibria, cycles, chaos), its flexibility (many models can be approximated in the form of a LV model) and its mathematical calculability.

In this article, we consider large LV models with random parameters. Leveraging on the asymptotic understanding of large random matrices which naturally appear enables us to provide insights on equilibria and species coexistence for such models.

Model and assumptions.

Large Lotka-Volterra systems of differential equations arise in various scientific fields such as biology, ecology, chemistry, etc. Although our results are generic in nature and not specific to a given field, we will rely on the ecological terminology in the sequel.

A large system of Lotka-Volterra equations is a system of coupled ordinary differential equations (ODE) that write:

dx k (t) dt = x k (t)   r k -θx k (t) + ∈[n] B k x (t)   , (1)
where k ∈

[n] = {1, • • • , n}.
Here, n represents the number of species in a food web or community, the unknown vector x = (x k ) k∈ [n] is the vector of abundances of the various species and evolves with time t > 0 according to the dynamics (1). Parameter r k represents the intrinsic growth rate of species k, θ is an intraspecific feedback coefficient (most often positive due to competition), and B k is the per capita effect of species on species k.

Remark 1. Notice that without interactions, i.e. B = (B k ) k, ∈[n] = 0 , system (1) is simply a system of uncoupled logistic differential equations.

We shall focus on the model where r k = θ = 1:

dx k dt = x k (1 -x k + (Bx) k ) , k ∈ [n] (2)
with matrix B admitting the following representation:

B = A α √ n + µ n 11 T ,
where A = (A ij ) is a matrix with random standardized (EA ij = 0 and var(A ij ) = 1) independent and identically distributed (i.i.d.) entries with finite fourth moment, α > 0 is an extra parameter reflecting the interaction strength, and µ ∈ R represents an arbitrary trend of the interactions. The n × 1 vector 1 is a vector of ones.

Remark 2. Although matrix B is a complex random object, a result by Tao [Tao13, Theorem 1.7] fully describes its asymptotic spectrum: Assume that |µ| > 1/α, then for any fixed ε > 0, almost surely eventually all the eigenvalues of B but one are in the disk {z ∈ C : |z| ≤ 1/α + ε} while one extra eigenvalue takes the value µ + o(1).

Remark 2 is illustrated in Fig. 1. Figure 1: Spectrum of non-Hermitian matrix B in the complex plan (n = 1000, α = 1). In Fig. 1a, µ = 0 and the solid line circle represents the boundary of the circular law. In Fig. 1b, µ = 2 and there is an eigenvalue in the small dashed circle centered at 2, as predicted by [Tao13, Th. 1.7] -see also Remark 2.

Presentation of the main results.

Unique equilibrium. In the study of the behaviour of x(t) as t → ∞ the existence of an equilibrium x * to Eq. ( 2) is an important prior to any stability property of x(t). By equilibrium, we mean the existence of a vector x * = (x * k ) satisfying

x * k (1 -x * k + (Bx * ) k ) = 0 for k ∈ [n] . General results on LV systems state that x(t) > 0 (componentwise) as long as x(0) > 0 [START_REF] Hofbauer | Evolutionary games and population dynamics[END_REF]. However, a possible equilibrium x * will only verify x * ≥ 0, i.e. some components x * k will take the value zero. In Theorem 2, we provide sufficient conditions on the parameters α and µ to ensure the existence of a unique equilibrium. These conditions rely on the "typical" behaviour of the random matrix B in large dimension n → ∞.

Evaluating the number of surviving species. Given a unique equilibrium x * , an important question is to describe the set of surviving/vanishing species. In this perspective, we introduce the set

S = {i ∈ [n], x * i > 0} (3)
of surviving species. In Section 3, we provide a heuristics to compute asymptotically the ratio |S| n and understand via a system of equations the dependence between parameters α and µ and the number of surviving species. A complementary result addressing the elliptic random matrix model by means of theoretical physics methods can be found in [START_REF] Bunin | Ecological communities with lotka-volterra dynamics[END_REF] (dynamical cavity method) and in [START_REF] Galla | Dynamically evolved community size and stability of random Lotka-Volterra ecosystems[END_REF] (generating functional techniques).

Notice that in [START_REF] Bizeul | Positive solutions for large random linear systems[END_REF], Bizeul and Najim have studied a different normalization for α in the case µ = 0, namely α ∼ 2 log(n), to guarantee the survival of every species (feasibility of the equilibrium). Indeed, a consequence of Dougoud et al.'s results [DVR + 18] is that some species will go to extinction if α > 0 is fixed (i.e. does not increase sufficiently with n).

An empirical study of LV systems with changing interaction strengths. Equipped with results on the existence of a unique equilibrium, one pending question is to understand what happens when the coefficient α varies for the same matrix A and the same parameter µ. In particular, when the value of α increases above a certain critical value, all species will coexist [START_REF] Bizeul | Positive solutions for large random linear systems[END_REF]; conversely, for sufficiently low values of α, the existence of a feasible equilibrium is not warranted anymore, however a unique and stable equilibrium may exist. How species equilibrium abundances change between these two states and how |S| varies will be the focus of Section 4.

Notations

Denote by ρ(C) the spectral radius of matrix C, by C the spectral norm of matrix C, and by u the euclidean norm of vector u. We represent by δ x the Dirac measure at x:

δ x (E) = 1 if x ∈ E 0 else .
We denote by a.s.

--→ the almost sure convergence of random quantities and by weak ---→ the weak convergence of measures. Given a set S, we denote by |S| its cardinality.

Equilibrium and stability results

A primer on Random Matrix Theory. We first recall some results on Random Matrix Theory (RMT), which provides a number of valuable insights to understand the asymptotic behaviour of A. We begin by the almost sure (a.s.) convergence of the spectral radius and the spectral norm:

ρ(A/ √ n) a.s.
----→ n→∞ 1 and A/ √ n a.s.

----→ n→∞ 2 .

We also have the a.s. weak convergence of the spectral measure of A/ √ n to the circular law (see for instance [START_REF] Bordenave | Around the circular law[END_REF]):

(a.s.) 1 n k∈[n] δ λ k (A/ √ n) weak ----→ n→∞ 1 {x 2 +y 2 ≤1} π dx dy , where (λ k (A/ √ n); k ∈ [n]) is the spectrum of A/ √ n.
This convergence is illustrated in Fig. 1a.

The description of the spectral norm of the deterministic part of matrix B is more straightforward:

µ 11 T n = |µ| .
Notice that both the random and deterministic parts of matrix B do not vanish asymptotically and thus have a macroscopic effect on the dynamics of system (1), as recalled in Remark 2 where the asymptotic spectrum of B is described.

The non-invadability condition. A key element to understand the dynamics of the LV system (1) is the existence of an equilibrium

x * = (x * k ) k∈[n] such that x * k (1 -x * k + (Bx * ) k ) = 0 , ∀k ∈ [n] , x * k ≥ 0. ( 4
)
and the study of its stability, that is the convergence of a solution x to the equilibrium x * :

x(t) ---→ t→∞ x * if x(0) is sufficiently close to x * .
It is well known that for LV equations, the fact that x(0) > 0 (componentwise) implies that x(t) > 0 for every t > 0, but one can have some components x k (t) of x(t) vanishing to zero. As a consequence, we will only consider non negative equilibria x * ≥ 0 with possibly vanishing components.

Notice that the situation substantially differs whether x * > 0 or x * has vanishing component. In the former case, the equilibrium set of equations becomes a linear equation:

x * = 1 + Bx * .
In the latter case, the equilibrium equations are no longer linear. In the centered case µ = 0, the existence of a positive solution has been studied in [START_REF] Bizeul | Positive solutions for large random linear systems[END_REF] and requires α 2 log(n) (while we consider α fixed here). A naive and systematic way to solve (4) is to choose a priori a subset I ⊂ [n], to set the corresponding components x I = (x * i ) i∈I to zero, and to solve the remaining linear system:

x I c = 1 |I c | + B I c x I c .
If there exists x I c ≥ 0 that solves the previous equation, then x =

x I x I c satisfies (4) and is a potential equilibrium. The number of subcases I ⊂ [n] is 2 n and in particular grows exponentially as n → ∞.

In order to decrease the number of potential solutions to (4), we first notice that relying on standard properties of dynamical systems, see for instance [Tak96, Theorem 3.2.5], a necessary condition for the equilibrium x * to be stable is that 1

-x * k + (Bx * ) k ≤ 0 . (5)
The condition (5) is better known in ecology as the non-invadability condition [START_REF] Law | Permanence and the Assembly of Ecological Communities[END_REF]. In reference to the ODE (2), the requirement for a given species k ∈ [n] to be non-invasive is equivalent to:

1 x k dx k dt x k →0 + ≤ 0 . (6)
The main interpretation is as follows: if one adds species k with a very low abundance in the system, it will not be able to invade the system as a result of condition (6).

As a consequence, we will now focus on the following set of conditions:

   x * k (1 -x * k + (Bx * ) k ) = 0 for k ∈ [n] , 1 -x * k + (Bx * ) k ≤ 0 for k ∈ [n] , x * ≥ 0 componentwise . (7)
This casts the problem of finding a non negative equilibrium into the class of Linear Complementarity Problems (LCP), which we describe hereafter.

Linear Complementarity Problem (LCP). LCP is a class of problems from mathematical optimization which in particular encompasses linear and quadratic programs; standard references are [START_REF] Murty | Linear complementarity, linear and nonlinear programming[END_REF][START_REF] Cottle | The linear complementarity problem[END_REF]. Given a n × n matrix M and a n × 1 vector q, the associated LCP denoted by LCP (M, q) consists in finding two n × 1 vectors z, w satisfying the following set of constraints:

   z ≥ 0 , w = M z + q ≥ 0 , w T z = 0 ⇔ w k z k = 0 for all k ∈ [n] . (8)
Since w can be inferred from z, we denote z ∈ LCP (M, q) if (w, z) is a solution of (8).

A theorem by Murty [START_REF] Murty | On the number of solutions to the complementarity problem and spanning properties of complementary cones[END_REF] states that the LCP (M, q) has a unique solution (w, z) iff M is a P -matrix, that is:

det(M I ) > 0 , ∀ I ⊂ [n] , M I = (M k ) k, ∈I .
In view of (7), we look for x * ∈ LCP (I -B, -1).

The equilibrium x * and its stability. For a generic LV system

d y k (t) dt = y k (r k + (Cy) k ) , k ∈ [n] , (9)
Takeuchi and Adachi (see for instance [Tak96, Th. 3.2.1]) provide a criterion for the existence of a unique equilibrium y * and the global stability of the LV system.

Theorem 1 (Takeuchi and Adachi [START_REF] Takeuchi | The existence of globally stable equilibria of ecosystems of the generalized volterra type[END_REF]). If there exists a positive diagonal matrix ∆ such that ∆C + C T ∆ is negative definite, then LCP (-C, r) 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 * .
Combining this result (setting C = -(I -B)) with results from RMT, we can guarantee the existence of a globally stable equilibrium x * of (1) for a wide range of the set (α, µ). Denote by

A = (a, m) ∈ R * + × R : a > √ 2, m < 1 2 + 1 2 1 - 2 a 2 (10)
the set of admissible parameters. Proof. We have

I -B + I -B T = 2I -(B + B T ) = 2I - A + A T α √ n + 2µ n 11 * . Notice that 2I -(B + B T ) is positive definite iff the top eigenvalue of B + B T is lower than 2: λ max (B + B T ) < 2 (11)
We first focus on the random part (A+A T )/α which is a symmetric matrix with independent N (0, 2/α 2 ) entries above the diagonal (note that the distribution of the diagonal entries is different from the off-diagonal entries, with no asymptotic effect). In this case, it is well known that the largest eigenvalue of the normalized matrix (or equivalently its spectral norm since the matrix is symmetric) a.s. converges to the right edge of the support of the semi-circle law (see [BS10, Th.

5.2]):

λ max A + A T α √ n a.s. ----→ n→∞ 2 √ 2 α . (12)
In the centered case (µ = 0), condition (11) occurs if α > √ 2. We now consider the general case where µ = 0. Notice that the rank-one perturbation matrix P = 2µ n 11 * admits a unique non zero eigenvalue 2µ. Denote by Ǎ = A+A T α √ n . We are concerned with the top eigenvalue of the symmetric matrix Ǎ + P . Based on a result by Capitaine et al. [CDMF09, Th. 2.1], we have:

λ max ( Ǎ + P ) a.s ----→ n→∞ 2µ + 1 α 2 µ if µ > 1 √ 2α , 2 √ 2 α else.
This result is illustrated in Figure 3.

Assume first that µ ≤ 1 α √ 2 (corresponding to zone C in Fig. 2), then

λ max ( Ǎ + P ) a.s. ----→ n→∞ 2 √ 2 α , which is strictly lower than 2 (cf. condition (11)) if α > √ 2.
Hence λ max ( Ǎ + P ) is eventually strictly lower than 2 under this condition.

Assume now that µ > 1 α √ 2 (corresponding to zone B in Fig. 2), then

λ max ( Ǎ + P ) a.s. ----→ n→∞ 2µ + 1 α 2 µ .
We are interested in the conditions for which 2µ + 1 α 2 µ < 2 or equivalently

2α 2 µ 2 -2α 2 µ + 1 < 0 . (13)
An elementary study of the polynomial ξ

(X) = 2α 2 X 2 -2α 2 X + 1 yields that ξ's discriminant is positive if α > √ 2, ξ(µ ± ) = 0 ⇔ µ ± = 1 2 ± 1 2 1 - 2 α 2 ,
and ξ

1 α √ 2 < 0, so that 1 α √ 2 ∈ (µ -, µ + ). In particular condition (13) is fulfilled if µ ∈ 1 α √ 2 , 1 2 + 1 2 1 - 2 α 2 .
Under this condition, (13) is fulfilled and a.s. lim sup n→∞ λ max ( Ǎ + P ) < 2, which completes the proof: we can then rely on Theorem 1 to conclude. 3 A heuristic approach to the proportion and distribution of the surviving species

(a) No outlier if µ ≤ (α √ 2) -1 . (b) Outlier if µ > (α √ 2) -1 .

Proportion of surviving species

In Section 2, we have presented conditions on parameters α, µ for the existence of a globally stable equilibrium x * to (1) under the non-invadability condition. As x * depends on the realization of matrix B, it is a random vector. Moreover since α > 0 is fixed and does not depend on n, the equilibrium x * will feature vanishing components (see the original argument in [DVR + 18] and the discussion in [START_REF] Bizeul | Positive solutions for large random linear systems[END_REF]). In an ecological context, we shall refer to these non-vanishing components as the surviving species, the vanishing components corresponding to the species going to extinction with x * k = 0 and x k (t) ---→ t→∞ 0 .

In this section, we assume that the A ij 's are N (0, 1)-distributed and describe the proportion of non-vanishing components of the equilibrium x * ; we also describe the distribution of the surviving species x * i > 0 which turns out to be a truncated Gaussian.

Remark 3. The Gaussianity assumption facilitates the explanation of the heuristics but does not seem necessary for the result to hold. In Fig. 6b, the entries are not considered Gaussian but the distribution of the surviving species still matches the truncated Gaussian.

Given the random equilibrium x * , recall the definition of S in (3). We introduce the following quantities:

p = |S| n , m = 1 |S| i∈[n] x * i , σ2 = 1 |S| i∈[n] (x * i ) 2 .
Notice that in the definitions of m and σ * we can replace i∈[n] by i∈S .

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 .
Recall the definition of the set A in (10).

Heuristics 1. Let (α, µ) ∈ A. The following system of three equations and three unknowns (p, m, σ) There is a strong matching between the parameters obtained by solving (14)-( 16) and their empirical counterparts obtained by Monte-Carlo simulations. This is illustrated in Fig. 4. In Fig. 5, we illustrate the sensitivity of σ * to the parameters (α, µ).

σ √ pΦ -1 (1 -p) + α(1 + µ p m) = 0 , ( 14
)
1 + µ p m + σ √ p α E(Z | Z > -δ) = m , (15)
(1 + µ p m) 2 + (1 + µ p m) 2σ √ p α E(Z | Z > -δ) + σ 2 p α 2 E(Z 2 | Z > -δ) = σ 2 (16) where δ = δ(p, m, σ) = α σ √ p (1 + µπm) , (17)
Remark 4. The heuristics above substantially simplifies in the centered model case, where µ = 0 and B = A α √ n . Following (12), assume that α > √ 2. Then the system with two unknowns (p, σ) 14)-( 16) and their empirical Monte Carlo counterpart (the star marker) as functions of the interaction strength α (left) and the interaction drift µ (right). Matrix B has size n = 500 and the number of Monte Carlo experiments is 200. In Column (4a), µ = 0 and α > √ 2 on the x-axis (which guarantees a unique and stable equilibrium x * ). When interaction α -1 increases, the number of surviving species p * decrease but their variance σ * and mean m * increase. In Column (4b), α = 2 and µ ∈ (-0.5, 0.5) on the x-axis. The interaction drift appears to have no impact on the proportion p * of surviving species, whereas it influences their variances and means.

   σ √ pΦ -1 (1 -p) + α = 0 1 + 2σ √ p α E(Z | Z > -δ) + σ 2 p α 2 E(Z 2 | Z > -δ) = σ 2 where δ = α σ √ p

Distribution of surviving species

In the previous section, the proportion p of the surviving species, their mean m and second moment σ2 have been described as empirical counterparts of the solutions p * , m * , (σ * ) 2 of a system of equations. While establishing this system of equations, we will provide the following representation (see (19)) of the abundance x * k of a surviving species:

x * k = 1 + µ p * m * + σ * √ p * α Z k ,
where Z k ∼ N (0, 1) and Z k > -δ * = -δ(p * , m * , σ * ), δ being defined in (17).

We take here advantage of this representation to characterize x k 's distribution, which turns out to be a truncated Gaussian.

Heuristics 2. Let (α, µ) ∈ A, x * the solution of (7) and let (p * , m * , σ * ) the solution of the system (14)-( 16). Recall the definition (17) of δ and denote by

δ * = δ(m * , p * , σ * ). Let x * k > 0 a positive component of x * , then: L(x * k ) ----→ n→∞ L 1 + µp * m * + σ * √ p * α Z Z > -δ * ,
where Z ∼ N (0, 1). Otherwise stated, asymptotically x * k admits the following density The heuristics simply follows from the fact that if x * k is a surviving species then

f (y) = 1 {y>0} Φ(δ * ) α σ * √ 2π p * exp -- 1 2 α σ * √ p * y -δ * 2 . (
x * k = 1 + µ p * m * + σ * √ p * α Z k
conditionally on the fact that the right hand side of the equation is positive, that is Z k > -δ * . A simple change of variable yields the density -details are provided in Appendix B. Fig. 6 illustrates the matching between the theoretical distribution obtained in Heuristics 2 and a histogram obtained by Monte-Carlo simulations. It also illustrates the validity of the heuristics beyond the gaussiannity assumption of the entries.

Construction of the heuristics

We first discuss Heuristics 1 and establish Equations ( 14), ( 15) and (16).

Equation (14). We first recall a result on order statistics of a Gaussian sample. Consider a family (Z k ) k∈[n] of i.i.d. random variables N (0, 1) and the associated order statistics

Z * 1 ≤ Z * 2 ≤ • • • ≤ Z * n . Consider an index nα ∈ [n] where α ∈ (0, 1) is fixed, then the typical location of Z * nα is Φ -1 (α): Z * nα Φ -1 (α) as n → ∞ , (18)
see for instance [START_REF] Smirnov | Limit distributions for the terms of a variational series[END_REF][START_REF] Balkema | Limit distributions for order statistics[END_REF].

Let x * be the equilibrium of (1) and consider the random variable

Žk = i∈S B ki x * i = (Bx * ) k .
We assume that asymptotically the x * i 's are independent from the B ki '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 two first moments can easily be computed, see Appendix B, Section B.1 for the details:

E x * Žk = µ p m and var x * ( Žk ) = pσ 2 α 2 .
Notice that the fact that E x * and var x * ( Žk ) only depend on p, σ and m which are (supposedly) converging quantities supports the idea that Žk is unconditionally a Gaussian random variable with moments:

E Žk = µ p * m * and var( Žk ) = p * (σ * ) 2 α 2 ,
where p * , m * , σ * are resp. the limits of p, m, σ. We now introduce the standard Gaussian random variables (Z k ) k∈ [n] where

Z k = Žk -E Žk var( Žk ) = α Žk -µ p * m * σ * √ p * . Consider the equilibrium x * = (x * k ) i∈[n] . If k ∈ S, that is x * k > 0, we have 1 -x * k + (Bx * ) k = 0 ⇒ 1 + (Bx * ) k > 0 .
This identity has two implications:

x * k = 1 + (Bx * ) k and 1 + (Bx * ) k > 0 if k ∈ S .
Relying on the representation (Bx * ) k = Žk , we obtain the representation

x k = 1 + (Bx * ) k = 1 + µ p * m * + σ * √ p * α Z k if k ∈ S . (19)
and the condition:

1 + (Bx * ) k = 1 + µ p * m * + σ * √ p * α Z k > 0 . If k / ∈ S then 1 + (Bx * ) k = 1 + µ p * m * + σ * √ p * α Z k ≤ 0
by the non invadability condition. Otherwise stated,

Z k ≤ -α(1+µ p * m * ) σ * √ p * if k / ∈ S , Z k > -α(1+µ p * m * ) σ * √ p * if k ∈ S .
Considering the order statistics of the Z k 's we obtain:

Z * 1 ≤ • • • ≤ Z * i ≤ - α(1 + µ p * m * ) σ * √ p * ≤ Z * i+1 ≤ • • • ≤ Z * n .
Now, there are exactly n -|S| = n(1 -p) indices before the threshold corresponding to the components of x * equal to zero. In particular, index i = n(1-p) corresponds to the value

Z * i - α(1 + µ p * m * ) σ * √ p *
Relying on (18), we finally obtain

Φ -1 (1 -p) = - α(1 + µ p * m * ) σ * √ p * .
It remains to replace p by its limit p * to obtain (14).

Equation (15). Our starting point is the following generic representation of an abundance at equilibrium (either of a surviving or vanishing species):

x * k = 1 + µ p * m * + σ * √ p * α Z k 1 {Z k >-δ * } = (1 + µ p * m * ) 1 {Z k >-δ * } + σ * √ p * α Z k 1 {Z k >-δ * }
Summing over S and normalizing,

1 |S| k∈S x * k = (1 + µ p * m * ) 1 |S| k∈S 1 {Z k >-δ * } + σ √ p * α 1 |S| k∈S Z k 1 {Z k >-δ * } , m (a) = (1 + µ p * m * ) + σ √ p * α n |S| 1 n k∈[n] Z k 1 {Z k >-δ * } , m (b) (1 + µ p * m * ) + σ √ p * α 1 P(Z > -δ * ) E(Z1 {Z>-δ * } ), m (1 + µ p * m * ) + σ √ p * α E(Z | Z > -δ * ).
where (a) follows from the fact that |S| = k∈S 1 {Z k >-δ * } (by definition of S), (b) from the law of large numbers

1 n k∈[n] Z k 1 {Z k >-δ} ----→ n→∞ EZ1 {Z>-δ} and |S| n ----→ n→∞ P(Z > -δ *
) with Z ∼ N (0, 1). It remains to replace m by its limit m * to obtain (15).

Eq. ( 16) can be obtained similarly. Details are provided in Appendix B, see Section B.2.

Switching between equilibria: changing interaction strength

In the previous sections, the strength α of the interactions was fixed, cf. equation (2). However, in nature interactions between species are constantly changing due to e.g. abiotic factors such as temperature, which affect the rate at which individuals forage for prey, etc. Our contribution is rooted in the framework of asymptotic dynamics, but many recent ecological studies highlight the importance of taking into account both transient dynamics (out-of-equilibrium abundance values due to frequent perturbations) and shifts between equilibria due to changing environmental conditions [Has01, FN11, NA16]. In the sequel, we discuss a more general framework.

The model and intuition. If we restrict ourselves to a specific environment, a possible ecological interpretation of the fluctuation of interaction strength corresponds to the relationship between the size of the habitat and the probability of contact between individuals from two interacting species (e.g. think of a pool of freshwater containing piscivorous fishes and their prey species -interactions, be them competition or predation, would be potentially more frequent if the volume of water was reduced). In physics, think of particles in motion in a given volume: if the number of particles and the temperature stay constant, reducing the volume should increase the number of interactions between particles. From a model standpoint, let µ be fixed, α = α(t) : R + → ( √ 2, ∞). Consider

dx k dt = x k (1 -x k + (B t x) k ) , k ∈ [n] , (20)
where matrix B t admits the following representation

B t = A α(t) √ n + µ n 11 T and (α(t), µ) ∈ A .
Remark 5. Following Theorem 2, condition (α(t), µ) ∈ A guarantees that there exists a unique equilibrium x * (t) for every t ∈ R + .

We focus on the case of a system that fluctuates between two equilibrium points (Figure 7a). We consider a sudden incident, most often irreversible in the short term, which reduces a portion of habitat, e.g. forest fires. The system transits from a feasible state to a state with vanishing species due to the change of the strength of interactions, modelled by the following step function:

α(t) = α 1 1 [0,t0) + α 2 1 [t0,+∞) , (α 1 , α 2 , t 0 ) ∈ ( √ 2, +∞) 2 × R + (21)
The change of model parameter occurs at t 0 which causes a change in the strength of the interactions going from a value α 1 to α 2 . One may choose α 2 < α 1 and the difference (or ratio) between the two values represents the intensity of the incident.

In large dimension, it is possible to characterize this change by its impact on the number of surviving species in the system (20). At a given time t, the proportion of surviving species p = p(t) ∈ [0, 1] can be computed by resolving the system in Heuristics 1. This function, associated to the step function α given in (21), is represented in Figure 7b. In the case of a sudden incident, the proportion of surviving species predicted by the heuristics has a form similar to α(t) i.e. a step response. In the feasible state, p(t) is close to 1 (i.e. all species coexist); after the transition occurs, some species vanish and here p(t) ≈ 0.87. Beware that the heuristics results follow instantaneously the change of α; however, there is a smoother transition in the dynamics between the two equilibria (respectively corresponding to α 1 and α 2 ) due to the return rate to equilibrium, see for instance [START_REF] Neubert | Alternatives to resilience for measuring the responses of ecological systems to perturbations[END_REF], [START_REF] Arnoldi | How ecosystems recover from pulse perturbations: A theory of shortto long-term responses[END_REF]). This transition is illustrated in Figure 8.

Simulations. We provide hereafter a simulation-based analysis of the impact of the sudden incident on a given ecosystem: Define a ten-species system (20) with a fixed matrix of interactions A with Gaussian entries N (0, 1) and consider α = α(t) as in Fig. 7a. This scenario has a mixed impact on the community, see Fig. 8. While some species benefit from this change through an increase in their abundances, others are severely affected by this shift, some of which become extinct. This phenomenon can be understood as follows: at first (t ≤ t 0 ), the system admits a feasible equilibrium state with α = α 1 = 2.5 > 2 log(10) 2, 14 and the abundances converge to this equilibrium (see Figure 7a). When the transition occurs at t = t 0 , Theorem 2 ensures the convergence to a new equilibrium defined by parameter α 2 . Since α 2 = 1.5 is below the feasibility threshold 2 log(10) 2.14, some species vanish. In other words, this sudden change of model parameter causes an increase of interaction strengths, which has a negative impact on species diversity. Evolution of diversity Finally, we illustrate the evolution of diversity using diversity indicators more suited to the description of changes such as the one represented in Fig. 8 [START_REF] Jost | Entropy and diversity[END_REF]. We introduce here Shannon diversity H , a standard measure of biodiversity in ecology, which is given by

H = - i x i j x j log x i j x j (22)
and ranges from 0 (one species completely dominates the community) to log(n), when each species is equally abundant. When many species become rare while others become more abundant, H drops. Because H varies before species actually vanish, it is a more sensitive index of community diversity than species richness.

The Hill number of order 1, defined as e H , is a diversity measure that takes into account species abundances and varies between 1 and n, i.e. it behaves like an "effective species richness", see e.g. [START_REF] Jost | Partitioning diversity into independent alpha and beta components[END_REF].

In Fig. 9a, we represent the mean of this diversity measure over time for a hundred-species system and observe a negative impact of the variation of the strength of the interactions on diversity. Parameter µ has no impact on diversity at equilibrium (similarly, µ has no impact on the proportion of persistent species), but the lower the value of µ, the slower the transition to a new equilibrium. In other words, the more generally "competitive" the ecosystem is (i.e. very negative values of µ), the longer it takes for transient dynamics to settle near equilibria.

The evolution of the Hill number of order 1 complements the evolution of species richness: when α decreases, the expected number of surviving species decreases (Fig. 7); at the same time, e H decreases even more drastically as the abundance distribution of surviving species gets more heterogeneous. Figure 9b also shows that the variability of the Hill number of order 1 among simulations increases drastically when α decreases. The conclusion is that the more species are lost, the more difficult it is to predict the diversity index as σ * depends on α and strongly influences e H . We define an interaction matrix A and let the dynamics of Lotka-Volterra evolve according to model (20) and we apply the function of variation α(t) of Figure 7a. For each time step, we compute e H . We repeat this scheme a large number of times (here 500), and we average the time series. In (a), we apply this procedure for different values of µ. In (b), we apply this procedure for a fixed µ = 0 and compute the quantiles of the 500 trajectories.

Theoretical estimation of diversity Standard mathematical methods (Taylor's theorem) can be used to obtain a theoretical approximation of the Hill number of order 1 (details are provided in Appendix B.4):

e H ≈ np * 3 2 - 1 2 (σ * ) 2 (m * ) 2 . ( 23
)
This estimator is based on the properties of the persistent species (p * , m * , σ * ) calculated by solving the fixed point equation of the heuristics (1). These three properties depend on the type of the interactions, as indicated by parameters (α, µ) (Figure 4). We compare the accuracy of this estimator with two examples in which the strength of the interactions (α) and the interaction drift (µ) vary (Figure 10). On the one hand, a shift of the interaction drift µ does not affect the proportion of surviving species. Furthermore, the impact of µ on σ * and m * is proportional i.e. σ * m * is equal to a constant. For these reasons, µ does not affect the Hill number (Figure 10b). On the other hand, if α increases, then p * , m * , σ * -----→ α→+∞ 1 and e H → n which is makes sense because when α becomes very large, all abundances converge to 1. If α decreases: p * decreases, and σ * increases faster than m * . This confirms that e H decreases even more drastically as the abundance distribution of surviving species gets more heterogeneous (Figure 10a).

Discussion

In this paper, our main interest was to describe the impact of the strength α and mean µ of interactions in large LV models on the conditions of coexistence of interacting species. Combining results from Takeuchi [Tak96] with standard RMT results, we have provided insights into the study of stability of large random ecosystems -see [START_REF] Stone | The feasibility and stability of large complex biological networks: a random matrix approach[END_REF][START_REF] Gibbs | Effect of population abundances on the stability of large random ecosystems[END_REF].

We have characterized the unique equilibrium properties of the surviving species by resolving a system of equations. From a physicist point of view, similar equations were obtained by Opper and Diederich [START_REF] Opper | Phase transition and 1/f noise in a game dynamical model[END_REF] and studied in a more general framework by Bunin using the dynamical cavity method [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] using generating functional techniques.

The coexistence of many species in random ecosystems was also studied by Servan et al. [SCG + 18] and Pettersson et al. [START_REF] Pettersson | Predicting collapse of complex ecological systems: quantifying the stability-complexity continuum[END_REF], where a more generic case was analyzed with different growth rates. Grilli et al. [GAS + 17] identified the key quantities regulating the parameter space leading to feasible communities. In contrast to previous approaches, an important feature of our model is the monitoring of the impact of interactions by a normalization factor (α √ n). From an ecological point of view, one might expect that the larger the number of species, the weaker the interactions will be due to some dilution of interactions among potential interaction partners, which would justify the use of such normalization factors. From a mathematical viewpoint, the normalizing parameter α captures the range of a unique equilibrium and the threshold for feasibility.

In nature, interactions between species are constantly changing and affected by the environment. Under the assumptions that environmental conditions influence interaction strengths, we have endeavoured in Section 4 to study the consequences of a sudden change of environmental conditions, expressed through a decrease in parameter α. Solving numerically the Lotka-Volterra system confirms the predictions given by heuristics, i.e. that a decrease in α negatively affects equilibrium species richness. A more precise representation of biodiversity dynamics can be obtained through Hill numbers of order 1 which also decreases after the sudden change in α. The dynamics of this diversity measure suggests that the mean of interaction coefficients, µ, affects the duration of transient dynamics, with shorter transient dynamics being associated with more mutualistic interactions (i.e. higher positive values of µ).

Many questions naturally arise as a follow-up. First, a mathematical proof of the heuristics presented here is a challenging prospect because the LCP procedure induces a statistical dependence that is a priori difficult to handle. However, looking for this proof will certainly help extend the results to other underlying assumptions on the parameters of the LV system.

Regarding the extension of the heuristics to other assumptions, two situations could be of particular interest: non-centered elliptical matrix models as in Bunin [START_REF] Bunin | Ecological communities with lotka-volterra dynamics[END_REF] and LV models in which species growth rates are also controlled as in [SCG + 18]. We are confident that such extensions are possible, but might hinge on more sophisticated developments, in particular to include growth rates in the calculus of order statistics.

In this paper we have only considered the case of a full interaction matrix with parameters (α, µ). However, food webs are often structured in compartments [BDB + 11] and/or obey hierarchies (e.g. larger species eat smaller ones) [BAB + 19]. By a numerical analysis of LV systems, one could use the same tools to study more patterned matrices [START_REF] Allesina | Stability criteria for complex ecosystems[END_REF]. Recent studies emphasize the sparsity of real food webs [START_REF] Busiello | Explorability and the origin of network sparsity in living systems[END_REF]. Beyond the feasibility studied by Akjouj and Najim [START_REF] Akjouj | Feasibility of sparse large lotka-volterra ecosystems[END_REF], one could also study the existence and stability of a unique equilibrium in a sparse context.

Finally, variations of the interaction strength highlight the impact of habitat destruction. Many theoretical studies provide mathematical formulas for the return rate to equilibrium [START_REF] Neubert | Alternatives to resilience for measuring the responses of ecological systems to perturbations[END_REF][START_REF] Arnoldi | How ecosystems recover from pulse perturbations: A theory of shortto long-term responses[END_REF]. A further theoretical study of model (20) could provide a more quantitative answer. In this article, we have limited the analysis to the case of a single sudden incident, but other types of fluctuations for the interaction strength could be considered for a better understanding of habitat conservation phenomena. For example, a seasonal model could be appropriate to describe the evolution of the dynamics over the seasons.

[Tak96] Y. Takeuchi. Global dynamical properties of Lotka-Volterra systems. World Scientific, 1996.

[Tao13] T. Tao. Outliers in the spectrum of iid matrices with bounded rank perturbations. Probability Theory and Related Fields, 155(1):231-263, February 2013.

[Vol26] V. Volterra. Fluctuations in the abundance of a species considered mathematically. Nature, 118(2972):558-560, 1926.

A Simulation details

Simulations were performed in Python. All the figures and the code are available on Github [START_REF] Clenet | Equilibrium and surviving species in a large lotka-volterra system of differential equations[END_REF].

Simulations on the properties of surviving species are performed in two different ways. The theoretical solutions are obtained resolving numerically the system of equations of heuristics 1. We use a solver (cf. scipy.optimize) to find a local minimum of the function defined by the system of equations (a modification of the Powell hybrid method). The empirical solutions are computed using a Monte Carlo experiment. We simulate a large number of matrix matrix B, we resolve the associated LCP problem using the Lemke's algorithm. Then, we use the LCP solution to calculate the properties of the surviving species: proportion of survivors, etc. Finally, we make an average on the ensemble of experiments. The Lemke algorithm is implemented in the lemkelcp package and can be found on Github [START_REF] Lamperski | Lemke's algorithm for linear complementarity problems[END_REF]. The dynamics of the Lotka-Volterra are achieved by a Runge-Kutta of order 4 (RK4) implemented in the code.

B Remaining computations

B.1 Details on Equation (14): Moments of Žk .

We compute hereafter the conditional mean and variance of Žk = (Bx * ) k with respect to x * . We rely on the following identities:

EB ki = µ n , E(B ki ) 2 = 1 α 2 n + µ 2 n 2 1 α 2 n , EB ki B kj = µ 2 n 2 (i = j) .
We first compute the conditional mean:

E x * ( Žk ) = i∈[n] E(B ki )x * i = i∈S E(B ki )x * i = µ n i∈S x * i , = µ |S| n 1 |S| i∈S x * i , = µ p m .
We now compute the second moment:

E x * ( Ž2 k ) = E x * i∈S B ki x * i 2 = E x * i,j∈S B ki B kj x * i x * j , = i∈S E(B 2 ki )x * 2 i + i =j E(B ki B kj )x * i x * j , = 1 α 2 n i∈S x * 2 i + µ 2 n 2 i =j x * i x * j , (a)
pσ 2 α 2 + µ 2 p2 |S| 2 i,j∈S

x * i x * j = pσ 2 α 2 + µ 2 p2 |S| 2 i∈S x * i 2 = pσ 2 α 2 + µ 2 p2 m2 ,
where the approximation in (a) follows from the fact that 1 |S| 2 i,j∈S

x * i x * j = 1 |S| 2 i =j x * i x * j + O 1 |S| .
We can now compute the variance:

var x * Žk = E x * Ž2 k -E x * Žk 2 = p σ2 α 2 .

B.2 Details on Equation (16).

As for the proof of (15), we start from the generic representation of x * k :

x

* k = 1 + µ p * m * + σ * √ p * α Z k 1 {Z k >-δ * } = (1 + µ p * m * ) 1 {Z k >-δ * } + σ √ p * α Z k 1 {Z k >-δ * } .
Taking the square, we get

x * 2 k = (1 + µ p * m * ) 2 1 {Z k >-δ} + 2(1 + µ p * m * ) σ * √ p * α Z k 1 {Z k >-δ} + (σ * ) 2 p * α 2 Z 2 k 1 {Z k >-δ * } .
Summing over S and normalizing, we get

1 |S| k∈S (x * k ) 2 = (1 + µ p * m * ) 2 1 |S| k∈S 1 {Z k >-δ * } + 2(1 + µ p * m * ) σ * √ p * α 1 |S| k∈S Z k 1 {Z k >-δ * } + (σ * ) 2 p * α 2 1 |S| k∈S Z 2 k 1 {Z k >-δ * } .
Finally, we conclude by replacing the empirical means by their limits It remains to replace σ by its limit σ * to obtain (16).

1 |S| k∈S Z i k 1 {Z k >-δ * } = E(Z i | Z > -δ * ) , i = 1, 2 .
B.3 Density of the distribution of the persistent species.

Assume that x * > 0, and les f = R → R be a bounded continuous test function.

We have

Ef (x * k ) = E f 1 + σ * √ p * α Z k + µ p * m * Z k > -δ * , = ∞ -∞ f 1 + µ p * m * + σ * √ p * α u 1 {u>-δ * } 1 -Φ(-δ * ) e -

B.4 Theoretical estimation of the diversity index

Recall that |S| = np is the number of surviving species and that

p i = x i j∈S x j
is the frequency of (surviving) species i.

To find a theoretical estimate of Hill number of order 1, we proceed by expansion and set

p i = 1 |S| + δ i , |δ i | 1 |S|
where δ i represents the deviation of species i from the standard frequency if all surviving species have the same frequency. Notice that i∈S δ i = 0.

H = - i∈S p i log(p i ) = - i∈S 1 |S| + δ i log 1 |S| + δ i .
We use the Taylor-Young formula of order 2 to decompose the log:

log 1 |S| + δ i = log 1 |S| + |S|δ i - δ 2 i |S| 2 2 + δ 3 i ε(δ i ) , ≈ log 1 |S| + |S|δ i - δ 2 i |S| 2 2 . H ≈ - i∈S 1 |S| + δ i log 1 |S| + |S|δ i - δ 2 i |S| 2 2 , = - i∈S 1 |S| log 1 |S| + δ i - δ 2 i |S| 2 + δ i log 1 |S| + |S|δ 2 i - δ 3 i |S| 2 2 , = log(|S|) - i∈S δ 2 i |S| 2 + i∈S δ 3 i |S| 2 2 .
Notice that

|S| i=1 δ 3 i |S| 2 2
is negligible since |δ i | |S| -1 . The term 1 corresponds to the maximum value that the Shannon diversity index can take if |S| are present in the system. It remains to develop the second term of the r.h.s.

-

1 2 i∈S δ 2 i |S| = - |S| 2 i∈S x i j∈S x j - 1 |S| 2 , = - |S| 2 i∈S x 2 i ( j∈S x j ) 2 - 2 |S| x i j∈S x j + 1 |S| 2 , = - |S| 2 i∈S x 2 i ( j∈S x j ) 2 + 1 2 , = - |S| 2 i∈S x 2 i |S| 2 ( 1 |S| j∈S x j ) 2 + 1 2 , = - 1 2 1 |S| i∈S x 2 i ( 1 |S| j∈S x j ) 2 + 1 2 , = - 1 2 σ2 ( m) 2 + 1 2 , = - 1 2 σ2 m2 -1 .
Finally the Hill number of order 1 can be computed as:

e H ≈ e log(|S|)-|S| 2 |S| i=1 δ 2 i , ≈ |S|   1 - |S| 2 |S| i=1 δ 2 i   = |S| 1 - 1 2 σ2 ( m) 2 + 1 2 = |S| 2 3 - σ2 ( m) 2 .

( a )

 a No outlier (µ = 0) (b) Presence of an outlier (µ = 2)

Theorem 2 .

 2 Let (α, µ) ∈ A, then a.s. matrix (I -B) + (I -B) T is eventually positive definite: with probability one, for a given realization ω, there exists N (ω) such that for n ≥ N (ω), (I -B ω ) + (I -B ω ) T is positive definite. In particular, there exists a unique (random) globally stable equilibrium x * ∈ LCP (I -B ω , -1) to (7).

Figure 2 :

 2 Figure 2: The shaded area represents the set A given by (10) yielding the existence of a unique (random) globally stable equilibrium x * . Area A is divided into two zones B and C. Both zones correspond to parameters (α, µ) for which matrix 2I -(B + B T ) is definite positive, as stated in Theorem 2. In zone B, λ max (B + B T ) corresponds to a spiked eigenvalue (µ above the critical threshold (α √ 2) -1 ). In zone C, λ max (B + B T ) corresponds to the right edge of the semicircle law. Notice that zone C extends to negative values along the y-axis.

Figure 3 :

 3 Figure 3: Spectrum (histogram) of the Hermitian random matrix B + B T (n = 1000, α = √2). The solid line represents the semi-circular law. In Fig.3a, µ = 0 and there is no oulier. In Fig.3b, µ =1.5 and one can notice the presence of an eigenvalue outside the bulk of the circular law. The dashed line indicates its theoretical value.

  admits a unique solution (p * , m * , σ * ) and p a.s. ----→ n→∞ p * , m a.s. ----→ n→∞ m * and σ a.s. ----→ n→∞ σ * . Associated to this solution (p * , m * , σ * ) is δ * = δ(p * , m * , σ * ).

  admits a unique solution (p * , σ * ). Moreover, p a.s. ----→ n→∞ p * and σ a.s. ----→ n→∞ σ * .

  (a) Parameters (p * , σ * , m * ) versus α. (b) Parameters (p * , σ * , m * ) versus µ.

Figure 4 :

 4 Figure 4: The plots represent a comparison between the theoretical solutions (p * , σ * , m * ) of (14)-(16) and their empirical Monte Carlo counterpart (the star marker) as functions of the interaction strength α (left) and the interaction drift µ (right). Matrix B has size n = 500 and the number of Monte Carlo experiments is 200. In Column (4a), µ = 0 and α > √ 2 on the x-axis (which guarantees a unique and stable equilibrium x * ). When interaction α -1 increases, the number of surviving species p * decrease but their variance σ * and mean m * increase. In Column (4b), α = 2 and µ ∈ (-0.5, 0.5) on the x-axis. The interaction drift appears to have no impact on the proportion p * of surviving species, whereas it influences their variances and means.

Figure 5 :

 5 Figure 5: The 3D plot represents σ * = σ * (α, µ), solution of the system (14)-(16). In contrast to the proportion of surviving species p * , we observe that µ has a major influence over σ * . The graph for the theoretical value of m * has approximately the same behavior with respect to µ and α.

  a) Gaussian entries. (b) Uniform entries.

Figure 6 :

 6 Figure 6: Distribution of surviving species. The x-axis represents the value of the abundances and the histogram is built upon the positive components of equilibrium x * . The solid line represents the theoretical distribution for parameters (α, µ) as given by Heuristics 2. In Fig. (6a), the entries are Gaussian N (0, 1) and the parameters are set to (n = 2000, α = 2, µ = 0.2). In Fig. (6b), the entries are uniform U(-√ 3, √ 3) with variance 1 and the parameters are set to (n = 2000, α = √ 3, µ = 0). Notice in particular that the theoretical distribution matches with non-Gaussian entries.

  Figure 7: (a) Variation of the interaction strength through time, used in the dynamics of a ten-species system depicted in Figure 8 (α 1 = 2.5, α 2 = 1.5). The dotted line represents the feasibility threshold associated to the system. (b) Variation of the proportion of surviving species depending of the variation of α(t) in (a).

Figure 8 :

 8 Figure 8: Abundance dynamics in the case of a community of ten species. The matrix of interactions A and the initial conditions are common and we apply the function of variation α(t) given in Figure 7a. The dashed lines represents species which benefit from habitat variation; solid lines represent species suffering from the change. Dotted lines represent species undergoing extinction.

Figure 9 :

 9 Figure9: Dynamics of the Hill number of order 1 in the case of an ecosystem of a hundred species. The initial conditions are similar for each species. We define an interaction matrix A and let the dynamics of Lotka-Volterra evolve according to model (20) and we apply the function of variation α(t) of Figure7a. For each time step, we compute e H . We repeat this scheme a large number of times (here 500), and we average the time series. In (a), we apply this procedure for different values of µ. In (b), we apply this procedure for a fixed µ = 0 and compute the quantiles of the 500 trajectories.

Figure 10 :

 10 Figure 10: Evolution of the Hill number of order 1 as a function of α (a) and µ (b). The theoretical solutions (solid line) are computed by resolving (1) and integrating the parameter (p * , σ * , m * ) in equation (23). The empirical solutions (star marker) are computed by a Monte-Carlo experiment (100 experiments): we define a matrix B of size 100 × 100, solve the LCP problem and calculate the associated Hill number e H using (22).

  µ p * m * ) 2 + 2(1 + µ p * m * ) σ * √ p * α E(Z | Z > -δ * ) + (σ * ) 2 p * α 2 E(Z 2 | Z > -δ * ) .

  * ) p * σ * dy , hence the density of x * k .

* Supported by CNRS Project 80 Prime -KARATE.

Replacing |S| by np * and σ and m by their limits, we get the desired result:


Guillaume Sarfati

Ananyo Maitra

Raphael Voituriez

Jean-Christophe Galas
email: jean-christophe.galas@sorbonne-universite.fr

Andre Estevez-Torres
email: andre.estevez-torres@sorbonne-universite.fr

André Estevez-Torres

Crosslinking and depletion determine spatial instabilities in cytoskeletal active matter

Introduction

Active matter is composed of agents that transform free energy into mechanical work, thus making non-equilibrium materials with attractive properties such as spatial self-organization and self-propulsion. [START_REF] Sanchez | Spontaneous motion in hierarchically assembled active matter[END_REF] In this context, cytoskeletal active matter, constituted of protein laments set in motion by molecular motors, is particularly interesting because the collective interaction of nanometric subunits creates order at the macroscopic scale. However, the determinants of spatial self-organization in these active systems are not yet well understood.

3-dimensional cytoskeletal active matter is prone to at least four dierent spatial instabilities: i) bending, [START_REF] Chandrakar | Connement Controls the Bend Instability of Three-Dimensional Active Liquid Crystals[END_REF] ii) buckling, 57 iii) global 811 and iv) local contractions. [START_REF] Torisawa | Spontaneous Formation of a Globally Connected Contractile Network in a Microtubule-Motor System[END_REF][START_REF] Alvarado | Molecular motors robustly drive active gels to a critically connected state[END_REF][START_REF] Nédélec | Self-organization of microtubules and motors[END_REF] These instabilities have been observed experimentally but so far in dierent systems, thus precluding a unied analysis of potential transitions between all possible phases. Further, there is currently no clear theoretical understanding of how a gel transitions from one instability mode to another.

Here we consider a mixture composed of microtubules (MTs), kinesin motor clusters that pull on them and a depletant that generates microtubule bundles. [START_REF] Sanchez | Spontaneous motion in hierarchically assembled active matter[END_REF] In the presence of ATP, this mixture generates an active solution in 3-dimensions that we will call `Dogic system' in the following. When the MTs are stabilized with GMPCPP they are shorter and the Dogic system generates continuous chaotic ows [START_REF] Sanchez | Spontaneous motion in hierarchically assembled active matter[END_REF] through an in-plane bending instability both in 2D [START_REF] Martínez-Prat | Selection mechanism at the onset of active turbulence[END_REF] and in 3D. [START_REF] Chandrakar | Connement Controls the Bend Instability of Three-Dimensional Active Liquid Crystals[END_REF] When the MTs are stabilized with taxol they are longer and chaotic ows are still observed but they emerge only after an out-of-plane buckling instability that creates a corrugated sheet. [START_REF] Senoussi | Tunable corrugated patterns in an active nematic sheet[END_REF][START_REF] Strübing | Wrinkling Instability in 3D Active Nematics[END_REF] Bending 15 and buckling [START_REF] Ideses | Spontaneous buckling of contractile poroelastic actomyosin sheets[END_REF] have also been reported in actin-myosin gels.

Besides bending and buckling instabilities, global and local contractions have been re-ported both in actin-myosin gels [START_REF] Bendix | A quantitative analysis of contractility in active cytoskeletal protein networks[END_REF][START_REF] Alvarado | Molecular motors robustly drive active gels to a critically connected state[END_REF] and in MT-motors mixtures. [START_REF] Torisawa | Spontaneous Formation of a Globally Connected Contractile Network in a Microtubule-Motor System[END_REF][START_REF] Foster | Active contraction of microtubule networks[END_REF][START_REF] Nédélec | Self-organization of microtubules and motors[END_REF][START_REF] Surrey | Physical properties determining selforganization of motors and microtubules[END_REF] In the Dogic system, global contractions appear with taxol MTs, [START_REF] Strübing | Wrinkling Instability in 3D Active Nematics[END_REF]1719 while local contractions arise with GMPCPP MTs. [START_REF] Senoussi | Programmed mechano-chemical coupling in reaction-diusion active matter[END_REF][START_REF] Senoussi | Programmed mechano-chemical coupling in reaction-diusion active matter[END_REF] Very recent experiments, however, have observed both types of contractions with GMPCPP MTs. [START_REF] Lemma | Active microphase separation in mixtures of microtubules and tip-accumulating molecular motors[END_REF] Finally, dynamic microtubules in the presence of multiheadedkinesins yield chaotic ows or local contractions in the absence of depletant [START_REF] Roostalu | Determinants of Polar versus Nematic Organization in Networks of Dynamic Microtubules and Mitotic Motors[END_REF] depending on the ratio between MT growth and motor speed, and the ratio between the motor and MT concentration. In summary, the disparity of conditions in which spatial instabilities have been observed in cytoskeletal active matter (Tab. S2) hinders our understanding of these systems. Here, we demonstrate that the four instabilities can be observed in the Dogic system by tuning only two parameters: ATP and depletant concentrations. The transition from bending to buckling and from local to global contractions happens at a critical ATP concentration that is associated to a critical concentration of passive motors that induces a gelation transition between a liquid-like and a solid-like state. [START_REF] Gagnon | Shear-Induced Gelation of Self-Yielding Active Networks[END_REF] A similar transition is observed at constant ATP concentration when the concentration of a passive linker is increased, which supports its generic physical origin, independently of the specic nature of the crosslinking agents.

Results and discussion

Our experiments take place in a long and thin channel of rectangular cross-section with typical dimensions L × W × H = 10 × 1 × 0.1 mm 3 . The active solution is composed of short GMPCPP-stabilized microtubules, clusters of biotinylated kinesin-1 in the presence of streptavidin, the depletant agent pluronic, that forms MTs bundles, and ATP coupled to an ATP-regeneration system (see Fig. 1a and SI Methods). MTs are uorescent and they are imaged using either confocal or epiuorescence microscopy. When the solution is pipetted into the channel, the shear ow generates an initial nematic order along the x-axis.

At high pluronic concentration (5.5 % w/v), the nematic state is deformed in two dierent ways depending on the ATP concentration, before chaotic ow occurs. At 50 µM ATP, the MT bundles deform in the xy-plane forming waves with a typical wavelength of 250 µm along the x-axis that arises within 5 minutes (Fig 1b,c and Movie S1). Confocal images show a minimal deformation of the bundles in the xz-plane. This in-plane bending instability has recently been observed in a similar system both in 3D 4 and in 2D. [START_REF] Martínez-Prat | Selection mechanism at the onset of active turbulence[END_REF] Within 10 -20 min, a superposition of bending instabilities results in a chaotic ow of bundles (Movie S2). At 5 µM ATP, in contrast, the MTs deform in the xz-plane, creating an undulating sheet with a typical wavelength of 200 µm within 150 min (Fig. 1d,e and Movie S3). Buckling happens concomitantly with a 10 % compression of the height of the gel. In the xy-plane, this sheet appears as a periodic structure of black and bright zones, with preserved nematic order of MT bundles along the x-axis. Within 1500 min this instability conduces to an onset of chaotic ows destroying nematic order, although the slow dynamics preclude the observation of a fully developed chaotic state (Movie S4). This out-of-plane buckling instability has recently been reported in a system that is similar but displays two important dierences: longer, taxol-stabilized MTs, and much higher [ATP] (2 mM). [START_REF] Senoussi | Tunable corrugated patterns in an active nematic sheet[END_REF][START_REF] Strübing | Wrinkling Instability in 3D Active Nematics[END_REF] Note however that, to the best our knowledge, the bending to buckling transition has not been observed so far. At a lower pluronic concentration (1.5 %) the behavior of the system changes dramatically (Fig. 2). In these conditions and high [ATP] (100 µM), a chaotic ow is observed during the rst 100 min but afterwards microtubules aggregate, their area growing with an exponential rate of 0.028 min -1 , rst to make 20 µm diameter clumps that coalesce after 1500 min into ∼ 250 µm diameter aggregates separated by a typical distance of ∼ 500 µm (Movie S5). We call this late instability local contractions. At the time resolution of the experiments, the type of instability that yields the initial chaotic ow could not be resolved. In contrast, at low [ATP] (5 µM), a global contraction in the absence of chaotic ow is observed (Movie S6).

The MT network rst contracts along y with an onset rate of 0.015 min -1 until it reaches a steady-state width of 0.3W at t = 1000 min and later contracts along x until contraction stops at 1700 min. This transition from local to global contraction reminds the one reported for actomyosin gels, where a critical concentration of passive linkers induced the percolation of the gel yielding the global contraction state. S3), with buckling being signicantly slower than bending, except at [ATP] c where replicate experiments yield either one or the other instability or a mixture of the two (Fig. S3), with identical time-scales.

The transition between local and global contractions also happens at [ATP] c = 10 µM, suggesting that ATP plays a similar role in both processes. Interestingly, for long, taxolstabilized MTs a similar phase diagram was observed, although with some dierences (Fig. S4).

First, the bending/buckling and local/global transition happened at a slightly higher ATP concentration. Second, the limit between contracting and extensile instabilities occurred at a lower depletant concentration. Finally, in some occasions, buckling and bending ended up producing, respectively, global and local contractions instead of chaotic ows. We hypothesize that these dierences mainly result from the increased length of taxol-stabilized MTs, which would favor gelation at a lower concentration of passive crosslinkers (and thus higher [ATP]) and the formation of bundles at a lower depletant concentration.

In a recent and complementary work, Lemma et al. established the phase diagram of a similar active gel as a function of tubulin, depletant and motor concentration at high ATP concentration. [START_REF] Lemma | Active microphase separation in mixtures of microtubules and tip-accumulating molecular motors[END_REF] They observed the following phases : global and local contractions, active foams, and bending, but not buckling. In particular, Lemma et al. found that increasing depletant promoted the transition from local to global contractions and from these to bending. In our case, depletant concentration induced a transition from global contractions to buckling and from local contractions to bending.

In the following, we will elucidate the eect of [ATP] on the spatial instabilities of the gel by considering two alternative hypothesis. In the kinetic hypothesis, the bending/buckling transition happens when there is a cross-over in the characteristic times of two concurrent processes : the passive contraction of the gel in the z direction promoting buckling, and the activity of the motors, favoring bending, that increases with [ATP]. In the gelation hypothesis, instead, [ATP] inuences the concentration of inactive motors that act as passive crosslinkers of the gel inducing a gelation transition below a critical ATP concentration.

We demonstrate that the gelation hypothesis is the only coherent with our experimental observations.

The bending/buckling transition is not controlled by passive contraction Our rst hypothesis is thus that the reduction in [ATP] reduces the activity and lets the time for the passive contraction to form a lm that buckles. In this kinetic hypothesis two processes are in competition: the passive contraction of MTs into a lm, whose time-scale would be independent of [ATP], and the activity due to the motors, whose dynamics may depend on [ATP] (Fig. 4a). To test it we rst measured the passive contraction dynamics along z in the absence of ATP (Figure S5). In the absence of pluronic but in the presence of motors a slight 5 % expansion of the gel was observed. The strongest contraction, 50 % of the thickness of the gel, was observed in the presence of pluronic alone and it reached 30 % when both kinesin and pluronic were present. The characteristic time of passive contraction in these two cases was in the range 150 -200 min (Fig. S5).

We then repeated the bending/buckling experiment in Fig. 1 but using caged-ATP to let the passive contraction happen in the absence of activity. Initially, the two systems at 5 and 50 µM caged-ATP evolved in the absence of free ATP and after 4 h a light pulse uncaged it (Fig. 4). In both cases, MTs contracted slightly along z during the rst 4 h (Fig. S5b). After uncaging, buckling with preserved nematic order along x was observed at 5 µM ATP within 15 min (Fig. 4b,d). At 50 µM ATP, uncaging resulted in a bending instability destroying nematic order along x (Fig. 4c,e), as in happened in Fig. 1c, although the bending was less clear. Indeed, this is expected because, as discussed in the following, a suspended lm that bends in the xy-plane necessarily deforms in the xz-plane. Importantly, however, pure buckling was not observed at 50 µM ATP even though the uncaging occurred after the MTs lm was formed. This result thus invalidates the kinetic hypothesis.

A passive crosslinker-based theory to explain the transition from one instability mode to another

Gagnon et al. used micro-rheology experiments to show that reducing [ATP] in the Dogic

system results in a gelation transition from a uid to a solid-like state [START_REF] Gagnon | Shear-Induced Gelation of Self-Yielding Active Networks[END_REF] because kinesin motors act as long-lived passive linkers in the absence of ATP. Our second hypothesis is thus that the gelation state of the active uid determines the observed spatial instability, from in-plane bending in the uid phase at high [ATP] to out-of-plane buckling in the solid phase at low [ATP] (Fig. 5a,b). In this subsection, we use a hydrodynamic active matter framework 2326 to show that the observation of dierent instability modes is indeed consistent with a gelation transition controlled by crosslinking.

We model the kinesin-microtubule gel using a consistent description of an orientable and crosslinked active material described in ref. 27 ; this general theory is recapitulated in the SI and summarized here. We assume that at low [ATP] most kinesin motors act as passive crosslinkers, while at high [ATP] most motors are activated and therefore do not act as passive crosslinkers. The concentrations of passive or active motors do not appear directly in our theory; instead we start from a classical model of viscoelasticity [START_REF] Beris | Thermodynamics of Flowing Systems: with Internal Microstructure[END_REF] where crosslinkers eectively control the relaxation rate 1/τ C of a conformation tensor C, which characterizes the polymeric conformation of the microtubule laments. τ C → 0 corresponds to a uid state (at high [ATP]), while τ C → ∞ corresponds to a solid state (at low [ATP]).

In such a permanently crosslinked, solid regime, the deviation of the conformation tensor, δC, from its steady-state value is the classical strain tensor of the gel. [START_REF] Hemingway | Active Viscoelastic Matter: From Bacterial Drag Reduction to Turbulent Solids[END_REF] Next, the degree of alignment of the microtubules is described by the nematic tensor Q, as is usual in the context of liquid crystals. Last, the dynamics of the gel is characterized by the uid ow v.

In the main text, we display the equations of motion corresponding to the situation in which the gel lls the entire channel (which is the case at least at high [ATP]), for simplicity; we average the dynamics of the gel over the thickness of the channel and present the eective two-dimensional dynamics. These two-dimensional hydrodynamic equations governing the dynamics of the coupled variables C, Q, v are As we show in detail in the supplement, this formulation demonstrates that in the uid limit (τ C → 0, associated to high [ATP]), and for high enough active drive ζ, the classical in-plane bending instability characteristic of extensile uids [START_REF] Marchetti | Hydrodynamics of soft active matter[END_REF][START_REF] Voituriez | Spontaneous ow transition in active polar gels[END_REF][START_REF] Aditi Simha | Hydrodynamic Fluctuations and Instabilities in Ordered Suspensions of Self-Propelled Particles[END_REF] is obtained (SI Eq. 24).

∂ t Q = Q • ω -ω • Q -2λA + 1 τ Q H (1) ∂ t C = C • ω -ω • C -2λ C A + 1 τ C B (2) Γv -η∇ 2 ⊥ v = -∇ ⊥ Π + ∇ ⊥ • σ G (3) where A = (1/2)[∇ ⊥ v + (∇ ⊥ v) T ] and ω = (1/2)[∇ ⊥ v -(∇ ⊥ v) T ]. H = -
However, when τ C → ∞ (associated with low [ATP]), i.e., the gel is essentially solid, the bending instability is suppressed when χ > 0 (SI Eq. 63). More quantitatively, in Fourier space, the stability of the gel to a pure bend perturbation with wavenumber q ⊥ is determined by the negativity of the following eigenvalue:

Ξ = - q 2 ⊥ τ Q + 4χ 2 τ C K - τ Q ζ Γ 1 -λ - 2λ C χτ C τ Q (4)
where K is the Frank elasticity of the nematic. Importantly, even when ζ K (i.e., in the highly active regime), if λ C → -1 and χ > 0, as we expect our system parameters to be, this eigenvalue goes from being destabilising (i.e., positive) to stabilising (i.e., negative) as τ C is increased. This demonstrates that greater crosslinking, which leads to a gelation transition, stabilises the gel to bend perturbations.

This thickness averaged description cannot, however, account for the buckling transition since homogeneity along the z direction is assumed. In the supplement, we explicitly demonstrate that when this condition is relaxed and gel disposes of some space to deform in the z direction, it buckles out of plane because the extensile active stress ζQ leads to a negative and destabilising surface tension along the uniaxial direction x, as previously shown 5 (SI Eq. 65). In other words, a buckling gel does not bend in the high τ C solid limit (at least initially). Finally, we also demonstrate that a gel in the uid, low τ C limit that may deform in the z direction, both bends and buckles (SI Eq. 71).

To summarise, the four key results of the theoretical modelling are i. in-plane bend instability at low τ C (uid limit), ii. suppression of this instability at high τ C (solid limit),

iii. buckling instability in the high τ C solid limit, when the gel does not homogeneously ll the channel and iv. both buckling and bending instability in the low τ C uid limit when the gel can deform in the z direction. These theoretical results are compatible with the previous experimental observations if we assume, based on Gagnon et al, [START_REF] Gagnon | Shear-Induced Gelation of Self-Yielding Active Networks[END_REF] that [ATP] drives a gelation transition in the Dogic system. In particular, the third and fourth theoretical results are consistent with buckling being observed at low [ATP] and a mixture of buckling and bending at high [ATP] when the active system may deform in the z axis, as shown in the experiments in Fig. 4.

Note that while we assumed that passive motors act as crosslinkers, the theory is insensitive to the chemical nature of crosslinkers. This suggests that adding passive crosslinkers in a gel, while keeping high ATP and pluronic concentrations, should modify the instability mode of the gel.

Our linear stability theory applies to buckling/bending instabilities that arise close to the initial nematic state but does not inform us about global/local contraction instabilities that happen far from this state. However, the global/local contraction transition has been demonstrated to be associated to a gelation transition in actin gels. [START_REF] Alvarado | Molecular motors robustly drive active gels to a critically connected state[END_REF] These results, together with Gagnon's et al. and Fig. 3 suggest that both the bending/buckling and the global/local contraction transitions at ∼ 10 µM ATP are due to a gelation transition.

The bending/buckling and local/global contraction transitions are controlled by the concentration of passive linkers To experimentally challenge the gelation hypothesis, we introduced passive PRC1 link-ers [START_REF] Chandrakar | Connement Controls the Bend Instability of Three-Dimensional Active Liquid Crystals[END_REF] at constant [ATP]. At high pluronic, we observed that 10 nM PRC1 was sucient to convert a bending into a buckling instability (Fig. 5c,d and Movie S7). At low pluronic, a system displaying local contractions in the absence of PRC1 displayed global ones at 20 nM PRC1 (Fig. 5e,f and Movie S8). In some occasions, we observed that a bending state became global contraction when PRC1 was added (Fig. S6). Importantly, the characteristic times of the bending/buckling instabilities, on the one hand, and of the local/global contraction instabilities, on the other hand, were similar at constant [ATP] and dierent PRC1 (Table S4), in agreement with experiments at [ATP] c (Fig. S3) and in contrast with what was observed at dierent [ATP] in the absence of PRC1 (Table S3).

The gelation transition occurred at [ATP] c ≈ 10 µM in Fig. 3. Using a Michaelis-Menten model, we can estimate the concentration of passive motors as

c p m = c 0 m (1 -[ATP]/(K M + [ATP]))
, where c 0 m is the concentration of all motors and K M = 40 -60 µM is the Michaelis-Menten constant of kinesin-1 motor and ATP. [START_REF] Schnitzer | Kinesin hydrolyses one ATP per 8-nm step[END_REF][START_REF] Hua | Coupling of kinesin steps to ATP hydrolysis[END_REF] Taking c 0 m = 25 nM for our experiments, one nds c p m ≈ 20 nM, which corresponds to a concentration of passive linkers c p l ≈ 5 nM (Fig. S7). This result is consistent with the gelation transition observed at 10 and 20 nM PRC1, respectively for bending/buckling and local/global contractions (Fig. 5c-f).

To further test the gelation hypothesis, we quantied the average velocity ows of the microtubules in the xy-plane during the bending and buckling instabilities at dierent [ATP] in the absence of PRC1 (Fig. 5g,h). During bending, strong ows up to 7 µm/min arise parallel to the y axis as observed previously. [START_REF] Chandrakar | Connement Controls the Bend Instability of Three-Dimensional Active Liquid Crystals[END_REF][START_REF] Martínez-Prat | Selection mechanism at the onset of active turbulence[END_REF] In contrast, during buckling the ows do not show a clear structure and they are signicantly smaller, on the order of 0.4 µm/min, as expected for a solid. The correlation along x of the velocity orientation angles conrms this interpretation (Fig. S8).

These observations clearly show that a large fraction of kinesin motors act as passive crosslinkers at low [ATP] but not at high [ATP], implying that [ATP] controls a transition from a crosslinked gel to a uid that determines the spatial instability displayed by the active gel.

Conclusion

We have observed four important spatial instabilities bending, buckling, local and global contractions in a single type of cytoskeletal active gel and we have shown that they are determined by two factors: the crosslinking of the gel and the concentration of depletion agent. Our theoretical and experimental results demonstrate that the bending-to-buckling and local-to-global contraction transitions are determined by a gelation transition in the active gel, induced either by motors acting as passive crosslinkers when [ATP] decreases, or by passive linkers at a xed [ATP]. However, although the depletant concentration is known to promote lament bundling, [START_REF] Needleman | Synchrotron X-ray Diraction Study of Microtubules Buckling and Bundling under Osmotic Stress: A Probe of Interprotolament Interactions[END_REF] it is not yet clear how this property controls the bending-to-local contraction and buckling-to-global contraction transitions.

Our work generalizes previous results that showed that the local-to-global contraction transition is controlled by the percolation of the gel induced by passive linkers in actin-myosin gels. [START_REF] Alvarado | Molecular motors robustly drive active gels to a critically connected state[END_REF] Building on previous work, [START_REF] Gagnon | Shear-Induced Gelation of Self-Yielding Active Networks[END_REF] our results underline the dual role of motors both as passive and active linkers in cytoskeletal active gels. In particular, we show, by varying [ATP], that motors not only control the degree of active drive, as is generally assumed in active uid theories, [START_REF] Senoussi | Tunable corrugated patterns in an active nematic sheet[END_REF][START_REF] Martínez-Prat | Selection mechanism at the onset of active turbulence[END_REF] but also the mechanical properties of the gel, leading to qualitatively distinct spatial instabilities and gel morphologies. In summary, this work provides a unied view of spatial instabilities in cytoskeletal active gels and provides an experimental and theoretical framework to design active materials with controllable out-of-equilibrium properties. [START_REF] Senoussi | Programmed mechano-chemical coupling in reaction-diusion active matter[END_REF]3436 Materials and methods Active solution: Active linkers are clusters of streptavidin and biotinylated kinesin-1.

The active solution was composed of uorescent GMPCPP-stabilized microtubules (at 0.5-1 mg/mL), kinesin and streptavidin (both at a concentration of 25 nM), ATP and pluronic (from 1 µM to 100 µM and from 1.5% to 5.5%, respectively), an ionic buer, an antioxydant mix and an ATP regeneration system. Bioinspired Soft Materials, DMR-2011846. In the concluding stages of our work, we became aware of a complementary, independent eort by the group of Guillaume Duclos who studied instabilities in kinesin-microtubule active uids. [START_REF] Najma | Dual antagonistic role of motor proteins in uidizing active networks[END_REF]

Table of Contents

Figure 1 :

 1 Figure 1: Bending and buckling instabilities are controlled by ATP concentration at high depletant concentration. a. Scheme of the experimental system with the microtubules in brown, the motor clusters in green and the pluronic depletant in orange, inside a rectangular channel. Cartoons of the microtubule morphology during the bending (b) and buckling (d) instabilities (the MTs appear in black). Time-lapse confocal images of uorescent microtubules showing the bending (c) and buckling (e) from the bottom (xy-plane) and from the side (xz-plane) of the channel (the MTs appear in white). 5.5 % pluronic, [ATP] = 50 µM for bending and 5 µM for buckling. Scale bars are 200 µm.

Figure 2 :

 2 Figure 2: Global and local contractions are controlled by ATP concentration at low depletant concentration. Cartoons of the microtubule morphology during local (a) and global (c) contractions (the MTs appear in black). Time-lapse epiuorescence images for local (b) and global contractions (d) in the yx-plane (the MTs appear in white). 1.5 % pluronic, [ATP] = 5 µM for global and 100 µM for local contractions. Scale bars are 500 µm.

10

 10

Figure 3 :

 3 Figure 3: Phase space of the four spatial instabilities as a function of pluronic and ATP concentrations. The symbols indicate the observed instabilities in a triplicate experiment and the colors are a guide to the eye.

Figure 4 :

 4 Figure 4: The bending/buckling transition is not controlled by passive contraction. a) Scheme of the kinetic hypothesis, where two instabilities are observed depending on ATP concentration. At low [ATP], activity is reduced, which lets the time for the passive contraction of microtubules to form a lm that buckles. At high [ATP], the higher activity induces a bending instability. Confocal images of solutions with respectively 5 µM (b, d) and 50 µM (c, e) caged-ATP. Panels b and c provide time-lapses in the xz-plane with the uncaging event indicated at 3 h 59 min. Panels d and e display images in the xy-plane just after uncaging. Scale bars are 100 µm in xy and 188 µm in z. Pluronic 5.5 % w/v.

  δF/δQ and B = -δF/δC are the molecular elds corresponding to the nematic tensor and the conformation tensor respectively where F is the free energy that would control the dynamics in the absence of activity. Crucially, F includes a coupling between alignment Q and conformation C of the gel with a coecient χ (SI Eq. 8), which is of the form 2χC ij Q ij , where Einstein summation convention has been assumed for the indices. Further, Q and C are also coupled via their respective couplings to the uid ow. Π enforces the two-dimensional incompressibility constrain ∇ ⊥ • v = 0 and σ G is the stress due to the gel which has the form described in ref.27. Activity is introduced in the model via an active stress ζQ which we assume to depend only on the nematic order parameter with ζ > 0 describing a system with an extensile uniaxial active stress. The remaining terms are required to ensure that the correct equilibrium physics is obtained in the absence of activity.

Figure 5 :

 5 Figure 5: The bending/buckling and local/global contraction transitions are controlled by the concentration of passive linkers that induce a gelation transition. (a) Schemes of passive linkers (ATP-depleted kinesin, in red, and PRC1 in black) and active one (ATP-consuming kinesin, in green). (b) Schemes of the active gel at dierent [ATP] and [PRC1] highlighting their liquid/solid-like state and the instabilities observed. Epiuorescence images showing the instabilities observed at high [ATP] (50 µM) and either high (c,d) or low (e,f) pluronic concentration and low (c,e) or high (d,f) [PRC1]. Scale bars are 200 µm for (c,d) and 500 µm for (e,f). g,h) Particle image velocimetry of the bending (g) and buckling (h) instabilities observed at 0 nM PRC1, 5.5 % pluronic and [ATP] as indicated. In panels c-f [MT] = 1 mg/mL.

Figure 6 :

 6 Figure 6: Four spatial instabilities arising in kinesin/microtubule active gels are determined by two factors: depletion and cross-linker activity.

Acknowledgment Z. Gueroui for a kind gift of the K401 plasmid, F. Lam from the microscopy platform at IBPS, L.L. Pontani for providing access to a spinning disk microscope, and A. Senoussi, Y. Vyborna and G. Duclos for insightful discussions. This work has been funded by the European Research Council (ERC) under the European's Union Horizon 2020 programme (grant No 770940, A.E.-T.). A.M. was supported by a TALENT fellowship from CY Cergy Paris Université. For PRC1 synthesis, we acknowledge support from Brandeis NSF MRSEC,

Imaging: The active solution was injected into rectangular cross-section channels (22 × 1.5 × 0.13 mm) made using passivated glass slides separated by paralm. Imaging was done by epiuorescence or confocal microscopy at room temperature.

Analysis: Images were reconstructed with Fiji and analyzed with Python scripts, using openpiv package for PIV analysis.

See Supporting Information for a more detailed description.

Supporting information