# Will a large complex model ecosystem be viable? The essential role of positive interactions

Rudolf P Rohr, Louis-félix Bersier, Roger Arditi

► To cite this version:

Rudolf P Rohr, Louis-félix Bersier, Roger Arditi. Will a large complex model ecosystem be viable? The essential role of positive interactions. Ecology, 2025, 106 (3), pp.e70064. 10.1002/ecy.70064. hal-04997006

HAL Id: hal-04997006

https://hal.sorbonne-universite.fr/hal-04997006v1

Submitted on 19 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Will a large complex model ecosystem be viable?
The essential role of positive interactions

Rudolf P. Rohr<fnref n="1" /> ID

Louis-Félix Bersier<fnref n="1" /> ID

Roger Arditi<fnref n="2" /> ID

**Correspondence**

Roger Arditi
Email: roger.arditi@gmail.com

**Funding information**

Swiss National Science Foundation,
Grant/Award Numbers: 320030L-227556,
31003A_182386, CRSII5_202290

**Handling Editor:** Tadashi Fukami

**Abstract**

Ecologists have documented many characteristics of natural systems that foster ecosystem persistence, and it might be deduced that such strategies are essential for counteracting the negative effect of complexity on local stability that was suggested by R.M. May in his influential work of the 1970s. However, we show that the loss of local stability does not necessarily imply total ecosystem extinction. A more general criterion of ecosystem viability is the long-term persistence of any number of surviving species—not necessarily all of them. With this approach, we show that persistence increases with complexity, contrary to previous theoretical findings. In particular, positive interactions (mutualistic or prey-to-predator) play a crucial role in creating ecological niches, which sustain biodiversity with increasing complexity.

**KEYWORDS**

complexity, model ecological communities, mutualistic interactions, persistence, positive interactions, prey-predator interactions, viability

## INTRODUCTION

The relationship between complexity and stability has long been debated in ecology (Allesina & Tang, 2012; Cohen & Newman, 1985; Elton, 1958; Hatton et al., 2024; MacArthur, 1955; May, 1972, 1973; McCann, 2000; Mougi & Kondoh, 2012, 2014; Odum, 1953; Pimm, 1984; Tilman & Downing, 1994; Tregonning & Roberts, 1979). Laboratory and field experiments (Hooper & Vitousek, 1997; McCann, 2000; McGrady-Steed et al., 1997; Naeem et al., 1994, 1995, 1996; Naeem & Li, 1997; Tilman et al., 1996) tend to indicate that more diverse ecological communities are more productive and less variable. These findings are in line with the conventional wisdom of ecologists that greater biodiversity favors ecosystem health (Elton, 1958; MacArthur, 1955; Odum, 1953). However, these findings also seem to contradict the

pioneering theoretical studies (Gardner & Ashby, 1970; May, 1972, 1973) that demonstrated that randomly built complex mathematical models tend to be unstable and, thus, that the corresponding ecosystems are unlikely to exist in nature. The rationale originally used by May (1972) to prove the nonviability of large complex systems was as follows: take *S* dynamic variables with random interactions, assume that they are in equilibrium, and examine the local (mathematical) stability of this equilibrium. Pairwise interactions can be of all types: competition (− −), antagonism (− +), mutualism (+ +), commensalism (+ 0), and amensalism (− 0). An equilibrium is locally stable if, after a small disturbance, the system returns to it. The authors who used this approach (Gardner & Ashby, 1970; May, 1972, 1973) observed that the probability of being in a stable local equilibrium decreases when the number of species *S* increases and/or

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2025 The Author(s). Ecology published by Wiley Periodicals LLC on behalf of The Ecological Society of America.

<footnote n="1">¹Department of Biology, University of Fribourg, Fribourg, Switzerland</footnote>

<footnote n="2">²Institute of Ecology and Environmental Sciences (iEES-Paris), Sorbonne Université, Paris, France</footnote>
the connectance c (a measure of the interaction density) increases.

Most importantly, these authors did not investigate the nonequilibrium properties of unstable systems because they did not follow the time course of population trajectories. The limited power of computers at that time only permitted determining whether equilibria were stable or unstable. Resting on the assumption that existing ecological systems in nature are necessarily stable and equating viability with local stability, the conclusion was that complexity could not be an explanation for ecological viability.

Here, we argue that it is incorrect to dismiss all unstable systems: some unstable systems can be viable. Taking the case of a single population as the simplest possible example, the basic model of population biology is the Malthus equation.

$$ \frac{dN}{Ndt} = \lambda, $$

which predicts either unlimited growth (if $\lambda > 0$) or a decline to extinction (if $\lambda < 0$) but never nonzero stable equilibrium. With the help of a model, every conservation ecologist would be happy to find that some endangered species has a growth rate $\lambda > 0$, suggesting that the species is not on the way to extinction. However, this also means, strictly speaking, that the population is on the way to infinity and, thus, is unstable. In this case, instability does not equate with extinction, on the contrary! Of course, the population must ultimately encounter some limitation to growth, which is not accounted for in the basic Malthus equation. Simply put, every model has a realm of applicability and, in a conservation context, when populations are low, the Malthus equation is sufficient.

How can this argument for a single population be extended to a multispecies community? Will an unstable system result in total extinction, or will some species grow indefinitely and thus be viable? If the latter is true, how many species are involved in the collective growth? What are the factors that determine the probability of this growth? What happens when growing populations ultimately encounter unavoidable limitations?

In this paper, we argue that many of the theoretical systems that were considered nonviable because of their instability actually correspond to biologically viable systems: those systems that “escape” from the local (mathematically unstable) equilibrium, with one or more species growing indefinitely. Of course, exploding species must ultimately reach bounds set by limiting resources that were not accounted for by the theoretical models. The only nonviable systems are those in which all populations decline to extinction.

## METHODS

A first building block in the mathematical description of a multispecies situation is the well-known Lotka-Volterra family of models for two species. Depending on the signs of the parameters that quantify the interspecific interactions, these models can describe all pairwise interactions mentioned in the Introduction (Gause & Witt, 1935; May, 1973, 1981; Odum, 1953). The dynamic outcomes of these systems are well known. Antagonism (e.g., predation) can lead to stable coexistence, cycling coexistence, or extinction. Situations with negative interactions only (competition, amensalism) lead either to coexistence or to alternative stable states with extinctions (Gause & Witt, 1935). Particularly interesting are the situations with positive interactions only (mutualism, commensalism): these interactions can lead to coexistence, to extinction, or to population explosion (Gause & Witt, 1935; Goh, 1979).

In order to describe the dynamics of larger communities with more than two species, we will use the classical extension known as the generalized Lotka-Volterra model (GLV; see below). Note that we assume that any pair of species is linked by a single type of interaction (competition, antagonism, mutualism, ...). We do not consider dual interactions, as defined by Mougi (2024).

We start by giving a general outline of the methodological rationale. First, given a number of species S, we sample a per capita interaction strength matrix A (of size $S \times S$) and a vector of equilibrium positive abundances $N^*$ (of size S). Based on this interaction matrix and this equilibrium vector, a unique GLV model is built in such way that $N^*$ is a positive equilibrium of this model. Then, the eigenvalues of the Jacobian matrix evaluated at this point are calculated. Given these eigenvalues, the local stability of the equilibrium point is determined by applying the Routh-Hurwitz criterion.

Second, we generate at random a perturbed vector $N_p$ that is very close to the equilibrium vector $N^*$. This perturbed vector is used as the initial condition for simulating numerically the GLV model parameterized as above. If the equilibrium point $N^*$ is locally stable, the dynamics converge back, of course, to $N^*$; but, if it is unstable, two possibilities can occur:

1. The dynamics converge to another equilibrium point, possibly with some species becoming extinct, or even all of them (system extinction).

2. One or more species can blow up with no limit. In this case, we stop the simulation when any species reaches an abundance of 1000 (chosen as an arbitrarily large number).
At the end of the simulation, each species with an abundance greater than $10^{-4}$ is considered alive, while if its abundance is lower, it is considered to be extinct.

Third, we repeat these two steps for a gradient of initial community richness S ranging from 10 to 250 (or 350) species in increments of 10 and simulate 1000 random communities for each richness S.

In order to avoid unrealistic abundance growth to infinity, which can occur in the second step above, we also devised a modification of the GLV model in such a way that populations cannot grow uncontrolled. The assumption is that any species must ultimately encounter limits to growth. In this modified model, each species is bounded to an abundance $N_{max}$, with its per capita growth rate declining asymptotically to zero when approaching $N_{max}$. This is not mechanistic modeling but the simplest phenomenological model that avoids unlimited growth. The dynamics of this modified model always converge to an equilibrium point. Note that the modified model is built in such a way that the sampled equilibrium point $N^*$ at the first step is still an equilibrium point of the modified model.

At the end of the simulation, we count the number of species still alive as well as the number of species that reached the maximum abundance $N_{max}$. In our numerical simulations, we set $N_{max} = 1000$. We repeat the simulations as in the third step.

Finally, in order to understand the effects of positive and negative interactions, we repeat the modified GLV model simulations: first, with interaction matrices containing only negative or null interactions, and second, with only positive or null interspecific interactions and negative intraspecific interactions (intraspecific competition).

## GLV model

The dynamic equations of the GLV model are

$$\frac{dN_i}{dt} = N_i \times \left( r_i + \sum_{j=1}^{S} a_{ij} \times N_j \right).$$

Given a sampled vector $N^*$ and a sampled interaction matrix $A$, $N^*$ must be a positive equilibrium point of the GLV model. This is ensured by setting the vector of intrinsic growth rates $r$ as

$$r = -AN^*.$$

Hofbauer and Sigmund (1998, p. 4, Theorem 5.2.1) proved mathematically that, in the GLV model, the

existence of a strictly positive trajectory requires the existence of a strictly positive equilibrium point (possibly unstable). In other terms, no trajectory can exist in the absence of an equilibrium. Thus, the existence of an equilibrium is not a restrictive assumption.

The model dynamics are simulated using the numerical integrator lsodsa of the library deSolve (Soetaert et al., 2010) with the R software (R Core Team, 2024). The R code that we used for performing the simulations and for drawing the figures is publicly available (Rohr, 2025). The simulation is stopped when any species reaches an abundance of 1000 or when all species have converged to an equilibrium (i.e., $|dN_i/dt| < 10^{-5}, \forall i$). In theory, it is also possible that one or more species settle on a never-ending cyclic (or chaotic) regime. For this reason, we set a simulation time limit of 1000 steps. However, this situation never occurred in our simulations (see Discussion).

## Sampling the interaction matrices and the equilibrium abundances

The per capita interaction strength matrix $A = [a_{ij}]$ is generated as follows:

1. We set the diagonal elements, that is, the intraspecific competition strength, to $-1$; that is, $a_{ii} = -1$.

2. The matrix $A$ contains $c S (S-1)$ nonzero off-diagonal elements $a_{ij}$. The parameter $c$ is the connectance, that is, the proportion of nonzero interspecific interaction coefficients. The magnitude of these interaction coefficients is sampled at random in a Gaussian distribution of mean zero and of variance $\sigma$, that is, $a_{ij} \sim N(0, \sigma)$.

3. The remaining interspecific interactions are set to zero.

The elements of the equilibrium abundance vector $N^* = [N_i^*]$ are sampled following a log-Gaussian distribution with a mean of zero and a SD of 1, that is, $\ln(N_i^*) \sim N(0, 1)$. Consequently, the equilibrium abundances are, on average, equal to $e^{1/2}$.

## Sampling the perturbations

The initial conditions for simulating the GLV model are set by a perturbation of the equilibrium abundances $N^*$. The elements $N_{p,i}$ of the vector of perturbed abundances $N_p$ are obtained by adding a normally distributed random number of mean zero and SD 0.02 to the equilibrium abundances $N_i^*$, that is, $N_{p,i} \sim N(N_i^*, 0.02)$. Note that, in order to avoid negative abundances in the very rare case
in which the sampled perturbed elements are negative (i.e., $N_{p,i} < 0$), we take their absolute values.

## Modified GLV

In order to avoid unlimited growth, we use the following modification of the GLV model:

$$ \frac{dN_i}{dt} = N_i \times \left( r_i + \sum_j a_{ij} \times N_j \right) \times \left( \frac{N_{\max} - N_i}{N_{\max}} \right). $$

Compared with the original GLV model, we apply the multiplicative term $(\frac{N_{\max} - N_i}{N_{\max}})$. This term avoids unlimited growth by bounding, from above, the abundance of each species to the value $N_{\max}$. When the abundance is below this threshold, the sign of the whole expression is given by the sign of the first bracket, that is, the original GLV model. Moreover, an equilibrium point $N^*$ of the original GLV model remains an equilibrium point of the modified model. When a population converges to the abundance $N_{\max}$, its per capita growth rate declines and converges to zero. Thus, this multiplicative term sets the upper bound $N_{\max}$ to every species abundance. In the numerical simulations, we set $N_{\max} = 1000$. With no loss of generality, we can assume that all species are renormalized in such a way that they all have the same $N_{\max}$.

In order to further understand the properties of the modified model, we now perform the phase-plane analysis of a system of two mutualistic species (Figure 1). The original Lotka–Volterra model is given by the following system of two differential equations:

$$ \frac{dN_1}{dt} = N_1 \times (r_1 + a_{11} \times N_1 + a_{12} \times N_2), $$

$$ \frac{dN_2}{dt} = N_2 \times (r_2 + a_{21} \times N_1 + a_{22} \times N_2), $$

while the corresponding modified Lotka–Volterra model is

$$ \frac{dN_1}{dt} = N_1 \times (r_1 + a_{11} \times N_1 + a_{12} \times N_2) \times \left( \frac{N_{\max} - N_1}{N_{\max}} \right), $$

$$ \frac{dN_2}{dt} = N_2 \times (r_2 + a_{21} \times N_1 + a_{22} \times N_2) \times \left( \frac{N_{\max} - N_2}{N_{\max}} \right). $$

In panels A and B of Figure 1, the blue and red dashed lines represent the nontrivial zero-growth isoclines for species 1 and 2, respectively. We assume the existence of a positive equilibrium, visualized by the yellow dot at the intersection of the isoclines, and we also assume that this equilibrium is unstable (a saddle point), that is, we

<figure><img src="image_2.png" /><figcaption>FIGURE 1 Phase-plane comparison between the original and the modified Lotka–Volterra models for two mutualistic species. Panel (A) shows, in the original Lotka–Volterra model, the typical pattern with an unstable nontrivial equilibrium point (yellow dot) leading either to population blow-up at the upper right or toward extinction (green dot). Panel (B) shows the changes brought about by the modified version of the Lotka–Volterra model. The panel clearly shows that the populations are now bounded by the new stable equilibrium point at the upper right (purple dot). Detailed explanations are given in the Methods section.</figcaption></figure>
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

assume $a_{12} \times a_{21} > a_{11} \times a_{22}$. The green dot is the trivial equilibrium (abundances equal to zero), which, in this case, is always locally stable. The arrows represent the vector field of the differential equations. In panel A, we see that, depending on the perturbation of the positive equilibrium (yellow dot), the system can either converge to extinction (green dot) or blow up to infinity. Panel B shows the effect of the extra term $(\frac{N_{max} - N_i}{N_{max}})$ on the behavior of the modified model. First, there are two new nontrivial zero-growth isoclines (blue horizontal and red vertical dashed lines) and their intersection defines a new equilibrium point represented by the purple dot. Second, examining the modification of the vector field, we can see that the system is prevented from blowing up, with the maximum abundance that can be reached being $N_{max}$. Third, depending on the direction of the perturbation of the unstable equilibrium (yellow dot), the system either converges to the trivial equilibrium (green dot), where both species become extinct, or to the new positive equilibrium (purple dot), where both species are alive.

## RESULTS

Our approach attempts to remain as close as possible to the original May approach with the crucial addition that, with numerical simulations of the GLV model, we are able to follow the time course of population abundances, particularly in the cases of unstable equilibria. The main results are presented in Figures 2–4.

In Figure 2, we represent the proportion of communities with a nontrivial stable equilibrium point (Figure 2A), the proportion of communities with at least one species exploding to 1000 (Figure 2B), and the number of species alive at the end of the simulation (Figure 2C), as functions of the initial community

**FIGURE 2** Effect of initial community richness on the outcome of the generalized Lotka–Volterra dynamics. Panel (A) represents the proportion of communities for which the equilibrium point is locally stable, that is, for which a small perturbation in the equilibrium abundance results in dynamics converging back to this equilibrium. Panel (B) shows the proportion of communities in which at least one species displays unlimited growth. Panel (C) shows the number of species alive at the end of the simulations. In panels (B and C), the simulations are stopped either if the system converges to a new equilibrium point or if at least one species reaches an abundance of 1000. The parameters for sampling the interaction matrices are $c = 0.7$ and $\sigma = 0.2$.
<figure><img src="image_6.png" /><figcaption>**FIGURE 3** Effect of the initial community richness on the outcome of the modified generalized Lotka–Volterra dynamics. Panel (A) shows, as expected, that the number of persistent species is lower than the number of initial species once the local stability of the initial equilibrium point has been lost. However, most importantly, the number of persistent species still increases when increasing the initial community richness. Panel (B) shows that, in the final community made up of the persistent species, the per capita interaction strength is, on average, greater than that in the initial community, demonstrating that positive interactions are ecologically selected. The parameters for sampling the interaction matrix are the same as those in Figure 2.</figcaption></figure>

richness *S*. Figure 2A confirms the original May results. As *S* increases, the probability that the equilibrium point is stable is 1 for lower richness (*S* < 50 in our example), but it declines very fast for higher richness, even declining to zero for sufficiently large *S*. Indeed, complexity begets instability. However, do unstable systems become extinct? In Figure 2B, we can see that this is clearly not the case. This panel shows that, in rich unstable communities, at least one species exhibits unlimited growth. Moreover, Figure 2C shows that species-rich communities can be obtained: this panel shows the number of species still alive at the end of the simulation. We observe that a number of species go extinct (because the red curve is lower than the dashed diagonal). However, after a transition, the number of species remaining alive increases linearly with increasing *S*, meaning that the proportion of lost species decreases. Very large viable communities can be built, and the greater the initial richness is, the greater the viable community richness.

In the above simulations, populations can grow with no limit: any population reaching the value of 1000 is considered on its way to infinity. For more realistic modeling, the GLV model is modified with an attenuation factor that sets $N_{max}$ as the asymptotic limit that no

population can exceed (see Methods). In Figure 3, we represent the number of live species at equilibrium (i.e., the number of persistent species) as a function of the initial richness *S*. We also calculate the change in the mean interaction strength, that is, the difference between the mean interaction strength among the species alive at equilibrium and the mean interaction strength among all initial species. A deviation from zero indicates ecological selection of interactions, which can be either positive (tendency to select positive interactions) or negative (tendency to select negative interactions).

Figure 3A shows that the modified GLV model does not alter qualitatively the main result of Figure 2C: starting from an initial community of richness *S*, a number of species are lost, but arbitrarily rich viable communities can still be built, which is shown by the fact that the average number of persisting species increases with no limit (red curve). Figure 3B shows that there is ecological selection of positive interactions: the mean interaction strength in the viable communities is greater than that in the initial pool.

The crucial role of the sign of mutual interactions is illustrated by Figure 4, in which we represent the proportion of communities for which no species becomes
<figure><img src="image_7.png" /><figcaption>**FIGURE 4** Effects on the dynamics outcome when the interspecific interactions are either all negative (competition and amensalism; panels A, C, E) or all positive (mutualism and commensalism; panels B, D, F). All intraspecific interactions are negative, and the parameters for sampling the interaction matrix are the same as those in Figure 2. Local stability is lost very rapidly (panels A and B). Panel D displays contrasting outcomes for the initially unstable equilibria with positive interactions: either the system converges to a new equilibrium without extinction (points on the main diagonal), or almost all species become extinct (points close to the abscissa). The probability of either of the outcomes occurring is approximately 0.5 (panel F).</figcaption></figure>
extinct. Species richness is still detrimental to stability, whether the interactions are all negative (Figure 4A) or all positive (Figure 4B). If they are all negative, persistent systems are readily obtained, but the number of persistent species cannot exceed a certain number (approximately 60–70 in the example of Figure 4C). If the interactions are all positive, contrasting outcomes are obtained at random (Figure 4D): either the system converges to a new equilibrium with all species alive (points on the diagonal) or almost all species go extinct (lower points, close to the x-axis). These two outcomes occur with an approximate 50–50 probability (Figure 4F).

Finally, we checked with additional simulations that the patterns illustrated by Figures 2 to 4 are not qualitatively sensitive to connectance *c* and/or to interaction strength variability *σ* (see Appendix S1: Figures S1–S11). The effects are quantitative only, with no remarkably qualitatively new pattern.

As previously mentioned, none of the numerous simulations had to be stopped with the time limit criterion. In the unmodified GLV model, in all simulations, each population either reached a nonzero equilibrium, went extinct, or exploded. In the modified GLV model, all populations converged to an equilibrium (either zero or nonzero). This means that never-ending trajectories like cyclic or chaotic regimes never occurred. We observed that alternative steady states do exist: when varying the perturbation only, one same system, with given parameters, can reach different equilibrium points, with different species surviving (see Appendix S1: Figure S12).

DISCUSSION

The well-known theoretical studies by R.M. May and followers seemed to contradict the fact that many complex viable ecosystems exist in nature. Our approach proposes a simple reconciliation. We argue that locally unstable systems must be considered viable when they escape some local equilibrium and grow to some other state (possibly losing a number of species but not all of them). Very species-rich, complex viable systems can readily be built. Richness and/or connectance are not obstacles to system viability; on the contrary: the richer the species pool is, and the more complex interspecific interactions are, the more easily robust, highly efficient communities can be formed. These properties cannot be observed when using local stability as a proxy of viability (Allesina & Tang, 2012; Hatton et al., 2024; May, 1972, 1973; Mougi & Kondoh, 2012, 2014).

This is entirely due to the effect of positive interactions, not only in mutualistic and commensalistic situations, but also in antagonistic situations. In their absence,

that is, in purely competitive systems, Figure 4C shows that species richness cannot exceed a certain limit. This limit depends on the connectance *c* and on the interaction strength variability *σ* (see Appendix S1: Figure S8): the higher *σ* and/or the higher *c*, the higher the interaction coefficients and the lower the numbers of coexisting species. In the language of niche theory, species packing limits the number of competing species that can coexist (MacArthur & Wilson, 1967). This result fully echoes the limiting similarity theory of competing species developed by MacArthur and Levins (1967) and its extension to species-rich competition communities by Vandermeer (1970). These authors demonstrated theoretically that there is a maximum number of species that can coexist and that this number is a function of the mean interaction strength, which is in line with our results displayed in Figure 4C and in Appendix S1: Figure S8.

When only positive interactions are allowed, that is, in purely mutualistic systems, it is well known that this is a highly destabilizing situation (Allesina & Tang, 2012; Goh, 1979). In approximately half of the cases, all populations go extinct, as shown in Figure 4D. However, the same figure shows that, in the other half of the cases, all species survive and grow to their ceiling populations. Thus, although mutualism is locally destabilizing, it is a factor that helps persistence, not the opposite.

Positive interactions create ecological niches and increase the number of species that can coexist. This has recently been theoretically supported by Koffel et al. (2021) by extending the classical niche theory, which considers only competition for common resources, to include positive feedback loops that permit niche extension and, consequently, increase the number of coexisting species. Empirically, positive interactions have also been widely studied in plant communities (Callaway, 2007), such as nitrogen fixation, shade for seedlings, or pollen and seed dispersal by pollinators and other animals (Bascompte & Jordano, 2014). Our study shows the importance of positive interactions for maintaining and extending niches, allowing an increasing number of coexisting species.

Methodologically, the major innovation of our work is that, using numerical simulation of dynamic systems, we have been able to follow the time course of populations, thus escaping from the limited scope of local stability analyses in ecology. This answers the complaint of Allesina and Tang (2012) that, although natural systems generally operate far from a steady state, the theoretical study of large systems is still based on local stability for feasibility reasons. With the latter approach, the essential role of positive interactions is condemned to silence.

Choosing the GLV model to represent community dynamics can be questioned. Indeed, alternate models
exist for pairwise interactions. Regarding pairwise competition, the original Lotka–Volterra model is generally accepted with no modification. Regarding pairwise mutualism, authors commonly introduce a saturation factor in order to avoid explosions (e.g., Mougi & Kondoh, 2012). Regarding pairwise predation, the Lotka–Volterra model of predator–prey interactions has been criticized, particularly by one of the present authors (Arditi & Ginzburg, 2012), because the model ignores predator saturation and predator interference. A large number of alternate predator–prey models have been proposed by many authors (e.g., Arditi & Ginzburg, 1989). However, it has proven impossible to generalize these alternate models to situations with more than one prey species and more than one predator species: all attempts failed to satisfy fundamental logical criteria (Arditi & Michalski, 1996). The GLV model remains irreplaceable in abstract theoretical studies such as the present article.

## AUTHOR CONTRIBUTIONS

**Concepts:** Roger Arditi, Rudolf P. Rohr, and Louis-Félix Bersier. **Methods:** Rudolf P. Rohr and Roger Arditi. **Investigation:** Rudolf P. Rohr, Roger Arditi, and Louis-Félix Bersier. **Figures:** Rudolf P. Rohr and Roger Arditi. **Writing—first draft:** Roger Arditi. **Writing—further versions:** Roger Arditi, Louis-Félix Bersier, and Rudolf P. Rohr.

## ACKNOWLEDGMENTS

This study was supported by Swiss National Science Foundation (SNSF) Project grant 31003A_182386 to Rudolf P. Rohr, Sinergia grant CRSII5_202290 to Rudolf P. Rohr, and Lead Agency grant 320030L-227556 to Rudolf P. Rohr.

## CONFLICT OF INTEREST STATEMENT

The authors declare no conflicts of interest.

## DATA AVAILABILITY STATEMENT

The R code (Rohr, 2025) for performing the simulations and for drawing the figures is available in Figshare at https://doi.org/10.6084/m9.fishshare.28378022.v1.

## ORCID

Rudolf P. Rohr https://orcid.org/0000-0002-6440-2696

Louis-Félix Bersier https://orcid.org/0000-0001-9552-8032

Roger Arditi https://orcid.org/0000-0002-5209-8754

## REFERENCES

Allesina, S., and S. Tang. 2012. "Stability Criteria for Complex Ecosystems." *Nature* 483: 205–8.

Arditi, R., and L. R. Ginzburg. 1989. "Coupling in Predator-Prey Dynamics: Ratio-Dependence." *Journal of Theoretical Biology* 139: 311–326.

Arditi, R., and L. R. Ginzburg. 2012. *How Species Interact: Altering the Standard View on Trophic Ecology*. New York, NY: Oxford University Press.

Arditi, R., and J. Michalski. 1996. "Nonlinear Food Web Models and Their Responses to Increased Basal Productivity." In *Food Webs: Integration of Patterns and Dynamics*, edited by G. Polis and K. Winemiller, 122–133. New York, NY: Chapman and Hall.

Bascompte, J., and P. Jordano. 2014. *Mutualistic Networks*. Princeton, NJ: Princeton University Press.

Callaway, R. M. 2007. *Positive Interactions and Interdependence in Plant Communities*. Dordrecht: Springer Netherlands.

Cohen, J. E., and C. M. Newman. 1985. "When Will a Large Complex System be Stable?" *Journal of Theoretical Biology* 113: 153–56.

Elton, C. S. 1958. *The Ecology of Invasions by Animals and Plants*. London: Chapman and Hall.

Gardner, M. R., and W. R. Ashby. 1970. "Connectance of Large Dynamic (Cybernetic) Systems: Critical Values for Stability." *Nature* 228: 784.

Gause, G. F., and A. A. Witt. 1935. "Behavior of Mixed Populations and the Problem of Natural Selection." *The American Naturalist* 69: 596–609.

Goh, B. S. 1979. "Stability in Models of Mutualism." *The American Naturalist* 113: 261–275.

Hatton, I. A., O. Mazzarisi, A. Altieri, and M. Smerlak. 2024. "Diversity Begets Stability: Sublinear Growth and Competitive Coexistence across Ecosystems." *Science* 383: eadg8488.

Hofbauer, J., and K. Sigmund. 1998. *Evolutionary Games and Population Dynamics*. Cambridge: Cambridge University Press.

Hooper, D. U., and P. M. Vitousek. 1997. "The Effects of Plant Composition and Diversity on Ecosystem Processes." *Science* 277: 1302–5.

Koffel, T., T. Daufresne, and C. A. Klausmeier. 2021. "From Competition to Facilitation and Mutualism: A General Theory of the Niche." *Ecological Monographs* 91: e01458.

MacArthur, R. H. 1955. "Fluctuations of Animal Populations and a Measure of Community Stability." *Ecology* 36: 533–36.

MacArthur, R. H., and R. Levins. 1967. "The Limiting Similarity, Convergence, and Divergence of Coexisting Species." *The American Naturalist* 101: 377–385.

MacArthur, R. H., and E. O. Wilson. 1967. *The Theory of Island Biogeography*. Princeton, NJ: Princeton University Press.

May, R. M. 1972. "Will a Large Complex System be Stable?" *Nature* 238: 413–14.

May, R. M. 1973. *Stability and Complexity in Model Ecosystems*. Princeton, NJ: Princeton University Press.

May, R. M. 1981. "Models of Two Interacting Populations." In *Theoretical Ecology: Principles and Applications*, 2nd ed. 78–104. Oxford: Blackwell.

McCann, K. S. 2000. "The Diversity–Stability Debate." *Nature* 405: 228–233.

McGrady-Steed, J., P. M. Harris, and P. J. Morin. 1997. "Biodiversity Regulates Ecosystem Predictability." *Nature* 390: 162–65.

Mougi, A. 2024. "Dual Species Interaction and Ecological Community Stability." *Ecology* 105: e4251.
Mougi, A., and M. Kondoh. 2012. "Diversity of Interaction Types and Ecological Community Stability." *Science* 337: 349-351.

Mougi, A., and M. Kondoh. 2014. "Stability of Competition-Antagonism-Mutualism Hybrid Community and the Role of Community Network Structure." *Journal of Theoretical Biology* 360: 54-58.

Naeem, S., K. Håkansson, J. H. Lawton, M. J. Crawley, and L. J. Thompson. 1996. "Biodiversity and Plant Productivity in a Model Assemblage of Plant Species." *Oikos* 76: 259.

Naeem, S., and S. Li. 1997. "Biodiversity Enhances Ecosystem Reliability." *Nature* 390: 507-9.

Naeem, S., L. J. Thompson, S. P. Lawler, J. H. Lawton, and R. M. Woodfin. 1994. "Declining Biodiversity Can Alter the Performance of Ecosystems." *Nature* 368: 734-37.

Naeem, S., L. J. Thompson, S. P. Lawler, J. H. Lawton, and R. M. Woodfin. 1995. "Empirical Evidence that Declining Species Diversity May Alter the Performance of Terrestrial Ecosystems." *Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences* 347: 249-262.

Odum, E. P. 1953. *Fundamentals of Ecology*. Philadelphia, PA: Saunders.

Pimm, S. L. 1984. "The Complexity and Stability of Ecosystems." *Nature* 307: 321-26.

R Core Team. 2024. *R: A Language and Environment for Statistical Computing*. Vienna: R Foundation for Statistical Computing.

Rohr, R. P. 2025. R Code for "Will a Large Complex Model Ecosystem Be Viable? The Essential Role of Positive Interactions." Figshare. https://doi.org/10.6084/m9 figshare.28378022.v1.

Soetaert, K., T. Petzoldt, and R. W. Setzer. 2010. "Solving Differential Equations in R: Package deSolve." *Journal of Statistical Software* 33: 1-25.

Tilman, D., and J. A. Downing. 1994. "Biodiversity and Stability in Grasslands." *Nature* 367: 363-65.

Tilman, D., D. Wedin, and J. Knops. 1996. "Productivity and Sustainability Influenced by Biodiversity in Grassland Ecosystems." *Nature* 379: 719-720.

Tregonning, K., and A. Roberts. 1979. "Complex Systems which Evolve towards Homeostasis." *Nature* 281: 563-64.

Vandermeer, J. H. 1970. "The Community Matrix and the Number of Species in a Community." *The American Naturalist* 104: 73-83.

## SUPPORTING INFORMATION

Additional supporting information can be found online in the Supporting Information section at the end of this article.

**How to cite this article:** Rohr, Rudolf P., Louis-Félix Bersier, and Roger Arditi. 2025. "Will a Large Complex Model Ecosystem Be Viable? The Essential Role of Positive Interactions." *Ecology* 106(3): e70064. https://doi.org/10.1002/ecy.70064
# Appendix S1

to the article

## Will a large complex model ecosystem be viable? The essential role of positive interactions

by

Rudolf P. Rohr, Louis-Félix Bersier, Roger Arditi

In this appendix, the supplementary figures S1 to S11 are provided to demonstrate the effects of the parameters σ and c on the patterns illustrated by Figs. 2 to 4 of the main document. The figures show that the effects are quantitative only, with no remarkable qualitatively new pattern.

Figure S12 supports the statement in the last sentence of the Results section in the main document: the figure shows that alternative steady states exist and that they do not change our findings.
<figure><img src="image_8.png" /></figure>

$\sigma = 0.1$ and $c = 0.7$

<figure><img src="image_9.png" /></figure>

$\sigma = 0.2$ and $c = 0.3$

<figure><img src="image_10.png" /></figure>

$\sigma = 0.15$ and $c = 0.7$

<figure><img src="image_11.png" /></figure>

$\sigma = 0.2$ and $c = 0.5$

<figure><img src="image_12.png" /></figure>

$\sigma = 0.3$ and $c = 0.7$

<figure><img src="image_13.png" /></figure>

$\sigma = 0.2$ and $c = 0.9$

<caption>Figure S1. Analogous to Fig. 2A in the main text, with various values of $\sigma$ and $c$.</caption>
<figure><img src="image_14.png" /><figcaption>Figure S2. Analogous to Fig. 2B in the main text, with various values of σ and c.</figcaption></figure>
<figure><img src="image_15.png" /></figure>

$\sigma = 0.1$ and $c = 0.7$

<figure><img src="image_16.png" /></figure>

$\sigma = 0.2$ and $c = 0.3$

<figure><img src="image_17.png" /></figure>

$\sigma = 0.15$ and $c = 0.7$

<figure><img src="image_18.png" /></figure>

$\sigma = 0.2$ and $c = 0.5$

<figure><img src="image_19.png" /></figure>

$\sigma = 0.3$ and $c = 0.7$

<figure><img src="image_20.png" /></figure>

$\sigma = 0.2$ and $c = 0.9$

<caption>Figure S3. Analogous to Fig. 2C in the main text, with various values of $\sigma$ and $c$.</caption>
<figure><img src="image_21.png" /><figcaption>Figure S4. Analogous to Fig. 3A in the main text, with various values of σ and c.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure S5. Analogous to Fig. 3B in the main text, with various values of σ and c.</figcaption></figure>
<figure><img src="image_23.png" /></figure>

$\sigma = 0.1$ and $c = 0.7$

<figure><img src="image_24.png" /></figure>

$\sigma = 0.2$ and $c = 0.3$

<figure><img src="image_25.png" /></figure>

$\sigma = 0.15$ and $c = 0.7$

<figure><img src="image_26.png" /></figure>

$\sigma = 0.2$ and $c = 0.5$

<figure><img src="image_27.png" /></figure>

$\sigma = 0.3$ and $c = 0.7$

<figure><img src="image_28.png" /></figure>

$\sigma = 0.2$ and $c = 0.9$

<caption>Figure S6. Analogous to Fig. 4A in the main text, with various values of $\sigma$ and $c$.</caption>
<figure><img src="image_29.png" /></figure>

$\sigma = 0.1$ and $c = 0.7$

<figure><img src="image_30.png" /></figure>

$\sigma = 0.2$ and $c = 0.3$

<figure><img src="image_31.png" /></figure>

$\sigma = 0.15$ and $c = 0.7$

<figure><img src="image_32.png" /></figure>

$\sigma = 0.2$ and $c = 0.5$

<figure><img src="image_33.png" /></figure>

$\sigma = 0.3$ and $c = 0.7$

<figure><img src="image_34.png" /></figure>

$\sigma = 0.2$ and $c = 0.9$

<caption>Figure S7. Analogous to Fig. 4B in the main text, with various values of $\sigma$ and $c$.</caption>
<figure><img src="image_35.png" /></figure>

$\sigma = 0.1$ and $c = 0.7$

<figure><img src="image_36.png" /></figure>

$\sigma = 0.2$ and $c = 0.3$

<figure><img src="image_37.png" /></figure>

$\sigma = 0.15$ and $c = 0.7$

<figure><img src="image_38.png" /></figure>

$\sigma = 0.2$ and $c = 0.5$

<figure><img src="image_39.png" /></figure>

$\sigma = 0.3$ and $c = 0.7$

<figure><img src="image_40.png" /></figure>

$\sigma = 0.2$ and $c = 0.9$

Figure S8. Analogous to Fig. 4C in the main text, with various values of $\sigma$ and $c$.
<figure><img src="image_41.png" /></figure>

$\sigma = 0.1$ and $c = 0.7$

<figure><img src="image_42.png" /></figure>

$\sigma = 0.2$ and $c = 0.3$

<figure><img src="image_43.png" /></figure>

$\sigma = 0.15$ and $c = 0.7$

<figure><img src="image_44.png" /></figure>

$\sigma = 0.2$ and $c = 0.5$

<figure><img src="image_45.png" /></figure>

$\sigma = 0.3$ and $c = 0.7$

<figure><img src="image_46.png" /></figure>

$\sigma = 0.2$ and $c = 0.9$

<caption>Figure S9. Analogous to Fig. 4D in the main text, with various values of $\sigma$ and $c$.</caption>
<figure><img src="image_47.png" /></figure>

<figure><img src="image_48.png" /></figure>

<figure><img src="image_49.png" /></figure>

<figure><img src="image_50.png" /></figure>

<figure><img src="image_51.png" /></figure>

<figure><img src="image_52.png" /><figcaption>Figure S10. Analogous to Fig. 4E in the main text, with various values of σ and c.</figcaption></figure>
<figure><img src="image_53.png" /></figure>

$\sigma = 0.1$ and $c = 0.7$

<figure><img src="image_54.png" /></figure>

$\sigma = 0.2$ and $c = 0.3$

<figure><img src="image_55.png" /></figure>

$\sigma = 0.15$ and $c = 0.7$

<figure><img src="image_56.png" /></figure>

$\sigma = 0.2$ and $c = 0.5$

<figure><img src="image_57.png" /></figure>

$\sigma = 0.3$ and $c = 0.7$

<figure><img src="image_58.png" /></figure>

$\sigma = 0.2$ and $c = 0.9$

<caption>Figure S11. Analogous to Fig. 4F of the main text, with various values of $\sigma$ and $c$.</caption>
<figure><img src="image_59.png" /><figcaption>Figure S12. This figure shows that alternative steady states exist and that they do not change our conclusions. The parameters for sampling the interaction matrices are $c=0.7$ and $\sigma=0.2$. For any given value of $S$, the various vertically aligned points result from different perturbations of the same initial equilibrium and the same interaction matrix. Since the points are spread vertically, this means that the various simulations reach different final states, with different species surviving.</figcaption></figure>