# Effects of stage structure on coexistence: mixed benefits

Gaël Bardon  Institute of Mathematics of Bordeaux, University of Bordeaux,
CNRS, Talence, France  Department of Polar Biology, Centre Scientifique de
Monaco, Principality of Monaco  Department of Ecology, Physiology and
Ethology, Institut Pluridisciplinaire Hubert Curien UMR 7178, University of
Strasbourg, CNRS, Strasbourg, France  Frédéric Barraquand  Institute of
Mathematics of Bordeaux, University of Bordeaux, CNRS, Talence, France

###### Abstract

The properties of competition models where all individuals are identical are
relatively well-understood; however, juveniles and adults can experience or
generate competition differently. We study here less well-known structured
competition models in discrete time that allow multiple life history
parameters to depend on adult or juvenile population densities. A numerical
study with Ricker density-dependence suggested that when competition
coefficients acting on juvenile survival and fertility reflect opposite
competitive hierarchies, stage structure could foster coexistence. We revisit
and expand those results. First, through a Beverton-Holt two-species juvenile-
adult model, we confirm that these findings do not depend on the specifics of
density-dependence or life cycles, and obtain analytical expressions
explaining how this coexistence emerging from stage structure can occur.
Second, we show using a community-level sensitivity analysis that such
emergent coexistence is robust to perturbations of parameter values. Finally,
we ask whether these results extend from two to many species, using
simulations. We show that they do not, as coexistence emerging from stage
structure is only seen for very similar life-history parameters. Such emergent
coexistence is therefore not likely to be a key mechanism of coexistence in
very diverse ecosystems, although it may contribute to explaining coexistence
of certain pairs of intensely competing species.

Keywords: coexistence; stage structure; competition; matrix population models

Correspondence: *<EMAIL_ADDRESS>$\dagger$
<EMAIL_ADDRESS>

##  1 Introduction

A basic tenet of demography and population ecology is that species vital
rates, such as survival and fertility, can vary with age or stage (Leslie,
1945; Cushing, 1998; Caswell, 2001). Surprisingly, little of this rich age-
structured population-level theory influences community-level coexistence
theory, that by contrast mostly builds on unstructured Lotka-Volterra
competition or consumer-resource models (Chesson, 2000; Barabás et al., 2016;
Letten et al., 2017). There have been, however, some calls to include more
demography into community ecology (e.g., Miller and Rudolf, 2011) and
coexistence studies in particular, as well as a number of empirically-driven
modelling studies doing so (Péron and Koons, 2012; Chu and Adler, 2015). They
show in general that the vital rates of the various life-stages do not react
in the same way to changes in the densities of juvenile and adult life-stages.
These differences may well create new avenues for coexistence.

Although mainstream coexistence theory often neglects the complexity of life
histories, some stage-structured coexistence theory has been previously
developed: several theoretical models, in past decades (Haigh and Smith, 1972;
Loreau and Ebenhoh, 1994) as well as in more recent literature (Moll and
Brown, 2008; Fujiwara et al., 2011) have considered complex life cycles in
which ontogenetic changes occur, leading to models where vital rates can be
differentially affected by competition. Such complexity can promote equalizing
mechanisms (Fujiwara et al., 2011) as well as niche differentiation (Moll and
Brown, 2008). A main observation of Moll and Brown (2008) was the occurrence
of so-called ‘emergent coexistence’ in two-species two-stages models, which
was defined as a coexistence equilibrium where competition on a single vital
rate would lead to exclusion, but the combination of competition processes
affecting two different vital rates leads to coexistence (e.g., competition
affecting juvenile survival leads to species 1 winning, competition affecting
fertility leads to species 2 winning, and yet both species coexist in the
model with both types of competition). This prompted the exciting idea that
community models with stage structure may foster more widespread coexistence
than classic Lotka-Volterra theory, suggesting stage structure as a key
missing ingredient in our explanations of diversity maintenance.

However, the results of Moll and Brown (2008) had two limitations that
prevented the generalization of this coexistence-enhancing role of stage
structure. First, their investigations were restricted to a particular life
cycle where juveniles and adults occupy separated habitats or niches (i.e.,
adults can only compete with adults and juveniles with juveniles), and these
complex life cycles (sensu Wilbur 1980) were even suggested to be potentially
responsible for the coexistence observed. Such segregation of competition
faithfully represents some systems with metamorphosis, such as odonates
(dragonflies and damselflies) or anurans (frogs), where there is an
ontogenetic niche shift (Werner and Gilliam, 1984). But such a model does not
correspond to say, birds and mammals where adults can impose strong
competition on the juveniles, as well as numerous other taxa of animals or
plants for which the niche shift after development is only partial, which
implies that the much larger adults can exert the bulk of the competition felt
by juveniles. Thus whether competition acting on multiple life-history
parameters could foster coexistence is still an open question for general life
cycles, and whether complex life cycles have anything special that promotes
coexistence is also something to uncover. Second, the analyses of Moll and
Brown (2008) relied entirely on numerical simulations, which means that while
we can observe emergent coexistence in the Ricker-based density-dependent
models that they considered, it is still difficult to understand fully the
phenomenon from a theoretical standpoint. Moreover, it is yet unclear if the
emergent coexistence equilibria highlighted by Moll and Brown (2008) are
structurally stable (robust to small deviations of the model framework, such
as a slightly different shape of density-dependence) and realistic (occurring
for a large range of parameter values). This highlights the need to
investigate further structured population models with additional model
structures and analytical techniques.

A further key question is whether emergent coexistence can occur in a
$S$-species context (where $S\gg 2$). Indeed, the emergent coexistence
observed by Moll and Brown (2008) resulted from a trade-off across the life-
history stages which allowed each species to the best competitor for a given
vital rate. While it is obvious that such a trade-off may occur with two
species, there are multiple ways to generalize such trade-offs in competitive
rankings to the $S$-species case, and not all of them necessarily lead to more
coexistence. In fact, unstructured Lotka-Volterra models suggest great caution
in generalizing results from two-species competition to many-species
competition, as the criteria for coexistence with many species usually
correspond to much broader niche separation between species than required for
two species, with intraspecific competition dominating (Barabás et al., 2016).
Similar phenomena might be expected for stage-structured competition models.

Here, we extend the results of Moll and Brown (2008) to a new model framework
using the canonical Beverton-Holt density-dependence (which is closer to the
continuous-time logistic model, Cushing et al., 2004) and life cycles where
juveniles can experience competition from adults. We present an analytical
invasion criterion applicable in this new model as well as to the
metamorphosis model structure numerically explored by Moll and Brown (2008).
This criterion allows us to better understand how and why emergent coexistence
occurs in both models, and is valid whenever the single-species models have
stable fixed points, which occurs for iteroparous life cycles. We also use the
sensitivity analysis developed by Barabás et al. (2014 _a_) to explore the
robustness of emergent coexistence equilibria. We find that emergent
coexistence is a general feature of structured two-species competition models
with two density-dependent vital rates, and that this equilibrium is not more
sensitive to perturbation of parameters than ‘classic’ coexistence equilibria
allowed by niche separation. We finally explore $S$-species models’ properties
using simulations of diverse community structures (number of initial species,
similarity of species) and for both life cycles. These many-species models
show the limitation of the emergent coexistence mechanism to explain
coexistence in highly speciose communities.

##  2 Models and Methods

###  2.1 Competition models

Our main model is a two-species difference equation with two life stages
(juvenile and adult) for each species. It is a development of models presented
in Fujiwara et al. (2011), which were restricted to the case where only one
matrix parameter per species (either fertility or survival) is affected by
densities. Our model extends this framework to multiple density-dependent
parameters per species, which provides a useful springboard to combine two
types of competition, affecting both fertility and survival rates. Our model
could also be viewed as an extension to iteroparous life-histories and delayed
maturation of the model of Cushing et al. (2007), who considered two-species
competition with semelparous life-histories and fast maturation. Our model
resembles Moll and Brown (2008)’s in that two vital rates are density-
dependent, and in that the life cycle includes iteroparity and a delayed
maturation; however, we consider like Fujiwara et al. (2011) and Cushing et
al. (2007) Beverton-Holt functions to model competition. These competition
models are sometimes called Leslie-Gower after Leslie and Gower (1958), and
are the closest discrete-time equivalents of continuous-time Lotka-Volterra
competition models (Cushing et al., 2004). They promote fixed points
equilibria, instead of cycles or chaos that can be generated by Ricker
functions at high fecundities (Neubert and Caswell, 2000), though see Cushing
et al. (2007) for semelparous life-histories and associated two-cycles. The
Beverton-Holt functions help to find analytically the equilibria whenever
multiple vital rates are density-dependent, which then fosters analytical
insights with invasion criteria.

In our main model (model 1), competition is generated by adults because it
felt more realistic for many animal species in which adults and juveniles
share the same habitat (in mammals and birds, see e.g. snowshoe hares
populations, Boutin, 1984 or avian scavengers populations, Wallace and Temple,
1987).

Without any competition, population densities of the two stages of a species
after one time step can be expressed (we arbitrarily pick species 1 for
illustrative purposes):

$n_{1j}(t+1)=n_{1j}(t)(1-\gamma_{1})s_{1j}+n_{1a}(t)f_{1}$ (1)

and

$n_{1a}(t+1)=n_{1j}(t)\gamma_{1}s_{1j}+n_{1a}(t)s_{1a}$ (2)

where $\gamma_{1}$, $f_{1}$, $s_{1j}$ and $s_{1a}$ respectively denote
maturation rate, fertility, juvenile survival rate, and adult survival rate of
species 1, all in $[0,1]$ except $f_{1}\in\mathbb{R}^{+}$. Additionally
requiring $\gamma_{i}>0$ and $s_{i,a}<1$ ($i\in\\{1,2\\}$) is necessary to
ensure that individuals do not get infinite youth or infinite lifespans,
respectively.

Then, we add competition on juvenile survival rates and fertility by making
these parameters density-dependent. In our first and main model, we consider
that competition is generated by adults of both species. Using the Beverton-
Holt model of competition, we get:

$f_{1}=\frac{\pi_{1}}{1+\alpha_{11}n_{1a}+\alpha_{12}n_{2a}}$ (3)

and

$s_{1j}=\frac{\phi_{1}}{1+\beta_{11}n_{1a}+\beta_{12}n_{2a}}$ (4)

where competition on fertility is composed of intra-specific competition
($\alpha_{11}n_{1a}$) and inter-specific competition ($\alpha_{12}n_{2a}$),
and similarly for juvenile survival competition. Here we denoted $\phi_{1}$
and $\pi_{1}$ for maximal juvenile survival and fertility of species 1.

With the vector $\mathbf{n}_{i}=(n_{ij},n_{ia})$, we have in matrix form
$\mathbf{n}_{i}(t+1)=\mathbf{A}_{i}(\mathbf{n}(t))\mathbf{n}_{i}(t)$. Our
model (model 1) can be written for one species with a projection matrix
$\mathbf{A}_{i}$, for instance with $i=1$,

$\mathbf{A}_{1}(\mathbf{n})=\begin{pmatrix}\frac{(1-\gamma_{1})\phi_{1}}{1+\beta_{11}n_{1a}+\beta_{12}n_{2a}}&\frac{\pi_{1}}{1+\alpha_{11}n_{1a}+\alpha_{12}n_{2a}}\\\
\frac{\gamma_{1}\phi_{1}}{1+\beta_{11}n_{1a}+\beta_{12}n_{2a}}&s_{1a}\end{pmatrix}=\begin{pmatrix}(1-\gamma_{1})s_{1j}(\mathbf{n})&f_{1}(\mathbf{n})\\\
\gamma_{1}s_{1j}(\mathbf{n})&s_{1a}\end{pmatrix}$ (5)

where $\alpha_{ij}$ and $\beta_{ij}$ are the competition coefficients
associated to the effect of species $j$ on species $i$ on fertility and
juvenile survival, respectively.

Although the mathematical derivations of equilibria do not impose that
$s_{i,a}>0$ (iteroparity), we additionally impose that condition because the
invasion and sensitivity analyses performed later require fixed point
equilibria, while semelparous life-histories also produce more varied
attractors, in particular 2-cycles (Cushing et al., 2007; Cushing and Henson,
2012).

We also considered a metamorphosis structure model (model 2) with a life cycle
similar to the one described by Moll and Brown (2008) and Cushing et al.
(2007), where juveniles compete only with juveniles and adults only with
adults. Unlike Moll and Brown (2008) but similarly to Cushing et al. (2007),
we chose Beverton-Holt density-dependence to model competition so that it can
still be compared to model 1. However, our model 2 also differ from Cushing et
al. (2007), who focus instead on a different life history context, with
semelparity ($s_{a}=0$), fast maturation ($\gamma=1$), and solely
intraspecific competition on fertility ($\alpha_{21}=\alpha_{12}=0$). In our
model 2, juvenile survival does not depend on adult densities (as life-stages
are separated) but depends on juvenile densities. We have a density-dependent
juvenile survival rate of species $1$ given by:

$s_{1j}^{M}(t)=\frac{\phi_{1}}{1+\beta_{11}n_{1j}(t)+\beta_{12}n_{2j}(t)}.$
(6)

This model 2 can be written for one species with a projection matrix
$\mathbf{A}_{i}^{M}$, for instance with $i=1$ and the same notations as for
model 1:

$\mathbf{A}_{1}^{M}(\mathbf{n})=\begin{pmatrix}\frac{(1-\gamma_{1})\phi_{1}}{1+\beta_{11}n_{1j}+\beta_{12}n_{2j}}&\frac{\pi_{1}}{1+\alpha_{11}n_{1a}+\alpha_{12}n_{2a}}\\\
\frac{\gamma_{1}\phi_{1}}{1+\beta_{11}n_{1j}+\beta_{12}n_{2j}}&s_{1a}\end{pmatrix}.$
(7)

We extended our models to $S$-species contexts, allowing each species to
potentially interact with every other. For model 1, we obtained the following
projection matrix (for species $i$):

$\mathbf{A}_{i}(\mathbf{n})=\begin{pmatrix}\frac{(1-\gamma_{i})\phi_{i}}{1+\sum_{k=1}^{S}\beta_{ik}n_{k,a}}&\frac{\pi_{i}}{1+\sum_{k=1}^{S}\alpha_{ik}n_{k,a}}\\\
\frac{\gamma_{i}\phi_{i}}{1+\sum_{k=1}^{S}\beta_{ik}n_{k,a}}&s_{i,a}\end{pmatrix}.$
(8)

The extension for model 2 is similar:

$\mathbf{A}_{i}^{M}(\mathbf{n})=\begin{pmatrix}\frac{(1-\gamma_{i})\phi_{i}}{1+\sum_{k=1}^{S}\beta_{ik}n_{k,j}}&\frac{\pi_{i}}{1+\sum_{k=1}^{S}\alpha_{ik}n_{k,a}}\\\
\frac{\gamma_{i}\phi_{i}}{1+\sum_{k=1}^{S}\beta_{ik}n_{k,j}}&s_{i,a}\end{pmatrix}.$
(9)

###  2.2 Scenarios of coexistence and associated parameter sets

We chose three scenarios of coexistence for illustration and sensitivity
analysis, matching those described by Moll and Brown (2008): (1) a classic
coexistence state that is indicated by both $\alpha$ and $\beta$ coefficients
(equivalent to the usual intra $>$ interspecific competition coefficients in
the factorised formulation of the unstructured Lotka-Volterra model, e.g.
Letten et al. 2017), (2) an emergent coexistence state where the competition
coefficients associated with the two vital rates indicate exclusion by the
opposite species, and (3) a coexistence state where there is a priority effect
indicated by competition coefficients acting on one vital rate while
competition acting on the other vital rate indicates coexistence. For
simplicity, we refer to these two last scenarios as ‘emergent coexistence’
since coexistence is not guaranteed by the conditions verified separately by
each set of coefficients. We note that coexistence scenario (3) does not
strictly correspond to an emergent outcome as defined by Moll and Brown (2008)
because competition coefficients for one of the two vital rates in fact
already indicate coexistence.

Although we consider below arbitrary parameters in Invasion analysis, we also
selected for numerical applications some parameter sets. To obtain parameter
sets promoting emergent coexistence, we set a trade-off between $\alpha_{ij}$
and $\beta_{ij}$ for both species (large $\alpha_{ij}$ is associated to low
$\beta_{ij}$ and vice versa).

To be able to compare the three scenarios, we set intraspecific competition
coefficients to a constant and only changed interspecific competition
coefficients. Note that we also chose the values of interspecific competition
coefficients such that we obtain the same densities at equilibrium across the
three coexistence scenarios. The equilibrium densities have been obtained
through simulation of the dynamics over 3000 time steps. No closed-form
solutions for the equilibrium densities are known for all possible values of
parameters, but some analytical progress can be made in special cases
(Supplement A). Parameter values for the three corresponding parameter sets
are given in Table 1.

Set | $\pi_{1}$ | $\pi_{2}$ | $\gamma_{1}$ | $\gamma_{2}$ | $\phi_{1}$ | $\phi_{2}$ | $s_{1a}$ | $s_{2a}$ | $\alpha_{12}$ | $\alpha_{21}$ | $\beta_{12}$ | $\beta_{21}$
---|---|---|---|---|---|---|---|---|---|---|---|---
1 | 30 | 25 | 0.7 | 0.8 | 0.5 | 0.4 | 0.5 | 0.6 | 0.05 | 0.06 | 0.06 | 0.06
2 | 30 | 25 | 0.7 | 0.8 | 0.5 | 0.4 | 0.5 | 0.6 | 0.02 | 0.112 | 0.125 | 0.01
3 | 30 | 25 | 0.7 | 0.8 | 0.5 | 0.4 | 0.5 | 0.6 | 0.043 | 0.035 | 0.155 | 0.165
Table 1: Parameter sets chosen to obtain equal equilibrium densities but
corresponding to contrasted scenarios of coexistence. We use intraspecific
coefficients $\alpha_{ii}=\beta_{ii}=0.1$.

The three parameter sets led to the same densities at equilibrium with
$n_{1}=200\pm 1$ and $n_{2}=130\pm 1$ (these have been fine-tuned by changing
interspecific interactions while keeping other parameters constant).

To check what coexistence outcomes were indicated by the competition
coefficients associated to fertility and juvenile survival, we used the
invasion criteria given by Fujiwara et al. (2011), valid for models with a
single density-dependent vital rate affected by competition (either fertility
or juvenile survival). For example, the invasion criteria of the model with
only fertility affected by competition (i.e., with only $\alpha$ coefficients)
is given by:

$\mathcal{R}_{\alpha}(i)=\frac{R^{\prime}_{i}}{R^{\prime}_{j}}\frac{\alpha_{jj}}{\alpha_{ij}}$
(10)

where

$R^{\prime}_{i}=\frac{\pi_{i}}{f_{i}^{(c)}}-1\\\ $ (11)

and

$f^{(c)}_{i}=\frac{1}{\gamma_{i}s_{j}}(1-s_{a})\left(1-s_{j}+\gamma_{i}s_{j}\right).$
(12)

This criterion $\mathcal{R}_{\alpha}$ reflects invasion outcomes if we only
had density-dependence on fertility, $\mathcal{R}_{\beta}$ does the same if we
only had density-dependence on juvenile survival. These invasion criteria aim
to confirm rigorously if, based on $\alpha$ competition coefficients only, or
$\beta$ coefficients only, we should conclude to coexistence or exclusion (the
well-known rule intrasp. competition $>$ intersp. competition only works in
certain parameterisations of unstructured two-species Lotka-Volterra models,
not structured ones). As shown in Table 2, the three parameter sets do
correspond to the three mechanisms mentioned above. Parameter set 1
corresponds to mutual invasibility in models with a single vital rate that is
density-dependent. Parameter set 2 matches emergent coexistence since 1
excludes 2 when competition is on fertility ($\alpha$ coefficients) while 2
excludes 1 when competition is solely on juvenile survival ($\beta$
coefficients). With parameter set 3, we have classical coexistence through
mutual invasibility indicated by the model with competition on fertility only
($\alpha$ coefficients) and a priority effect in the model with competition on
juvenile survival only ($\beta$ coefficients).

Parameter set | $\mathcal{R}_{\alpha}$ | $\mathcal{R}_{\beta}$ | Outcome
---|---|---|---
Set 1 |  2.29 / 1.46 |  1.76 / 1.58 | Classical coexistence
Set 2 |  5.72 / 0.78 | 0.84 / 9.47 | Emergent coexistence
Set 3 |  2.60 / 2.43 | 0.68 / 0.57 | Emergent coexistence
Table 2: Invasion criteria for species 1/species 2, evaluated for models where
a single vital rate is affected by competition. $\mathcal{R}_{\alpha}$ and
$\mathcal{R}_{\beta}$ are invasion criteria respectively for models with
competition only on fertility and juvenile survival. Invasion criteria whose
values are larger than 1 are in bold.

###  2.3 Invasion analysis of exclusion equilibria

We devised analytical invasion criteria for our model combining competition on
fertility and juvenile survival, starting with model 1 where all competition
is generated by adults. These differ from the invasion criteria of the
previous section in that both $\alpha$ and $\beta$ are now potentially non-
zero. We were inspired by the work of Cushing (1998, 2008) who studied the
stability of the exclusion equilibria (with the focal invading species absent)
to describe the ability of the focal species to invade a community. We
therefore looked for conditions equivalent to an unstable exclusion
equilibrium, which corresponds to the fact that the absent species can
actually invade. All equations are presented in Appendix 1.

We used the aggregated parameters $C=(1-\gamma)\phi$ and
$D=\frac{\pi\gamma\phi}{1-s_{a}}$ in the single-species model. By considering
the model with a single species (since at the time of invasion, the invader is
almost absent), we found an expression for the inherent net reproductive
number, $C+D=(1-\gamma)\phi+\frac{\pi\gamma\phi}{1-s_{a}}$. This number must
exceed 1 to have a viable population of the resident species. We then obtained
the densities at a stable positive equilibrium for the single-species model:

$\left\\{\begin{array}[]{r c l}n_{j}^{*}&=&(1-s_{a})\left(\frac{1+\beta
n_{a}^{*}}{\gamma\phi}\right)n_{a}^{*}\\\ n_{a}^{*}&=&\frac{(\alpha
C-\alpha-\beta)+\sqrt{(-\alpha
C+\alpha+\beta)^{2}-4\alpha\beta(-C-D+1)}}{2\alpha\beta}\end{array}\right.$
(13)

provided that

$C+D=(1-\gamma)\phi+\frac{\pi\gamma\phi}{1-s_{a}}>1.$ (14)

We then evaluated the stability of the exclusion equilibrium in the two-
species model: if it is locally asymptotically stable, the invader converges
to zero density; if not, the invader invades. To do this, we calculated the
eigenvalues of the Jacobian of the 2 species system evaluated at the exclusion
equilibrium. We found the following condition for stability when species 2 is
absent and species 1 is resident:

$\frac{C_{2}}{1+\beta_{21}n_{1a}^{*}}+\frac{D_{2}}{(1+\beta_{21}n_{1a}^{*})(1+\alpha_{21}n_{1a}^{*})}<1$
(15)

with $C_{2}=(1-\gamma_{2})\phi_{2}$ and
$D_{2}=\frac{\pi_{2}\gamma_{2}\phi_{2}}{1-s_{2a}}$ as denoted previously. This
last expression (eq. 15) gives us an invasion criterion, in the sense that if
the expression is larger than 1, species 2 is able to invade the community
when rare.

We applied the exact same method to model 2 with metamorphosis structure,
where juveniles compete with juveniles and adults compete with adults
(Supplement B).

We found analytical expressions for the boundaries of coexistence domains over
$\alpha_{12}$ and $\alpha_{21}$. These expressions are respectively given by
equations 16 and 17 that are derived directly from the invasion criteria:

$\displaystyle\alpha_{21}^{*}$
$\displaystyle=\left(\frac{D_{2}}{(1+\beta_{21}n_{1a}^{*})-C_{2}}-1\right)\frac{1}{n_{1a}^{*}}$
(16) $\displaystyle\alpha_{12}^{*}$
$\displaystyle=\left(\frac{D_{1}}{(1+\beta_{12}n_{2a}^{*})-C_{1}}-1\right)\frac{1}{n_{2a}^{*}}$
(17)

If both conditions $\alpha_{12}$<$\alpha_{12}^{*}$ and
$\alpha_{21}$<$\alpha_{21}^{*}$ are met, then there is coexistence. If both
conditions are violated, there is a priority effect, and if one condition is
met and the other is violated, there is exclusion of one of the species.

To visualize a case of emergent coexistence outcomes as did Moll and Brown
(2008) for their Ricker competition model, we represent the coexistence
domains over $\alpha_{12}$ and $\alpha_{21}$ with fixed values of $\beta_{12}$
and $\beta_{21}$, for model 1. We stress that here, these coexistence domains
arise from analytical expressions of the boundaries rather than numerical
simulations of the dynamical systems. We computed $\alpha_{12}^{*}$ and
$\alpha_{21}^{*}$ for the three parameter sets mentioned in section 2.2 where
$\beta$ coefficients indicated respectively coexistence, exclusion of species
1 and a priority effect.

###  2.4 Sensitivity analysis

To examine the sensitivity of equilibrium densities to changes in parameters
under different scenarios of coexistence, we applied the community-level
sensitivity analysis developed by Barabás et al. (2014 _a_) to model 1. The
computation of fixed point densities, required to perform the analysis, was
done numerically. With iteroparity, only stable fixed points were observed,
both in single-species and two-species competition models, as shown in
Supplement C. The two-species competition models fixed points do not have
closed-form expressions in the general case—although this may sound surprising
(but see Meszéna et al. 1997), it was indeed easier to compute invasion
criteria than to find analytically the two-species fixed point, that may be
analytically tractable only in special cases (see Supplement A). We sum up
here the notations and expressions used to compute the sensitivity. Explaining
all the steps would require paraphrasing too much of Barabás et al. (2014
_a_), so we simply refer the reader to their paper for more details. We use
$n_{i}=\mathbf{n}_{i}\mathbf{q}_{i}$ where $\mathbf{n}_{i}$ is the vector of
population density of the species $i$ at equilibrium for the various stages
and $\mathbf{q}_{i}$ a vector of weights given to stages. We considered only
the case where $\mathbf{q}_{i}$ is a vector of ones. We also need to define
$\mathbf{p}_{i}=\frac{\mathbf{n}_{i}}{n_{i}}$ the proportion of each stage of
species $i$. We denote $\mathbf{A}_{i}$ the projection matrix of species $i$
with its eigenvalues $\lambda^{k}_{i}$ and its left and right eigenvectors
$\mathbf{v}_{i}^{k}$ and $\mathbf{w}_{i}^{k}$. We also define the dominant
eigenvalue associated with its eigenvector : $\lambda_{i}$ and
$\mathbf{w}_{i}$. We normalized the vectors to have $|\mathbf{w}_{i}|=1$ and
$\mathbf{v}_{i}^{j}\mathbf{w}_{i}^{k}=\delta_{jk}$ with $\delta_{jk}=1$ if
$j=k$ and $0$ otherwise. Barabás et al. (2014 _a_) used regulating variables
$\mathcal{R}_{\mu}$ which correspond for example for species 1 in our model 1
to $\mathcal{R}_{\alpha,1}=\alpha_{11}n_{1a}+\alpha_{12}n_{2a}$ and
$\mathcal{R}_{\beta,1}=\beta_{11}n_{1a}+\beta_{12}n_{2a}$. Note that the
notation for regulating variables of Barabás et al. (2014 _a_) is similar to
that used for invasion criteria in single-species models that we previously
introduced, even though ecological meanings differ greatly. As regulating
variables are not used later in the text, there is no risk of confusion and we
simply follow the conventions of Barabás et al. (2014 _a_) in this section.

We then computed the sensitivity of population densities of both species to
perturbations of each parameter using the following expressions, where $E$ is
the perturbed parameter:

$\frac{dn_{i}}{dE}=\sum^{S}_{j=1}(M_{ij})^{-1}g_{j}$ (18)

where the matrix $M_{ij}$ is given by

$M_{ij}=\sum_{\mu,\nu}\left(\mathbf{v}_{i}\frac{\partial\mathbf{A}_{i}}{\partial\mathcal{R}_{\mu}}\mathbf{w}_{i}\right)\left(\delta_{\mu\nu}-\frac{\partial\mathcal{G}_{\mu}}{\partial\mathcal{R}_{\nu}}\right)^{-1}\frac{\partial\mathcal{R}_{\nu}}{\partial\mathbf{n}_{j}}\mathbf{p}_{j}$
(19)

and the vector $g_{j}$ is given by

$g_{j}=\mathbf{v}_{j}\frac{\partial\mathbf{A}_{j}}{\partial
E}\mathbf{w}_{j}+\sum_{\sigma,\rho}\left(\mathbf{v}_{j}\frac{\partial\mathbf{A}_{j}}{\partial\mathcal{R}_{\sigma}}\mathbf{w}_{j}\right)\left(\delta_{\sigma\rho}-\frac{\partial\mathcal{G}_{\sigma}}{\partial\mathcal{R}_{\rho}}\right)^{-1}\frac{\partial\mathcal{G}_{\rho}}{\partial
E}.$ (20)

The function $\mathcal{G}_{\mu}(\mathcal{R}_{\nu},E)$, introduced by Barabás
et al. (2014 _a_) to simplify notations, is given by

$\mathcal{G}_{\mu}(\mathcal{R}_{\nu},E)=\sum_{j=1}^{S}\left(\frac{n_{j}}{\mathbf{q}_{j}\mathbf{w}_{j}}\frac{\partial\mathcal{R}_{\mu}}{\partial\mathbf{n}_{j}}\sum^{s_{j}}_{k=2}\frac{1}{\lambda_{j}-\lambda_{j}^{k}}\times\left(\mathbf{w}^{k}_{j}-\frac{\mathbf{q}_{j}\mathbf{w}^{k}_{j}}{\mathbf{q}_{j}\mathbf{w}_{j}}\mathbf{w_{j}}\right)\otimes\mathbf{v}_{j}^{k}\right)\mathbf{A}_{j}(\mathcal{R}_{\nu},E)\mathbf{w}_{j}.$
(21)

With these expressions, we were able to compute the sensitivity to
perturbations of parameter values for the three parameter sets / contrasted
scenarios of coexistence. These have near-equal densities at equilibrium (see
section 2.2), to allow for meaningful comparisons.

###  2.5 $S$-species simulations

We simulated the dynamics of larger communities (5, 10 and 40 species) with
parameters chosen to favour emergence of coexistence through life-history
complexity, using the extension of model 1 and of model 2 in $S$-species
contexts (eqs. 8 and 9). Our idea to extend the mechanism behind emergent
coexistence from two to an arbitrary number of species was to create opposite
competitive hierarchies in the competition coefficients related to fertility
($\alpha$) and juvenile survival ($\beta$). We therefore devised a mechanism
by which a species strongly competitive on fertility would be weakly
competitive on juvenile survival and reciprocally. This was done through a
negative correlation between the coefficients $\alpha_{ij}$ and $\beta_{ij}$.
The negative correlation was parameterized using a bivariate normal
distribution and a covariance between the two variables equal to
$-0.9\sigma^{2}$, where $\sigma^{2}$ is the variance of both $\alpha_{ij}$ and
$\beta_{ij}$.

To draw the full set of parameter values, we have used:

  * •

Log-normal fertilities
$\pi_{i}\sim\text{Log-}\mathcal{N}(\mu_{\pi},\sigma_{\pi}^{2})$

  * •

Beta-distributed survival and transition probabilities
$\phi_{i}\sim\mathcal{\text{Beta}}(a_{\phi},b_{\phi})$,
$s_{a,i}\sim\mathcal{\text{Beta}}(a_{s_{a}},b_{s_{a}})$ and
$\gamma_{i}\sim\mathcal{\text{Beta}}(a_{\gamma},b_{\gamma})$

  * •

Normally distributed intraspecific competition coefficients
$\alpha_{ii}\sim\mathcal{N}(\mu_{\alpha},\sigma_{\alpha}^{2})$ and
$\beta_{ii}\sim\mathcal{N}(\mu_{\beta},\sigma_{\beta}^{2})$

  * •

As mentioned above, when $i\neq j$ we use interspecific coefficients drawn as
$(\alpha_{ij},\beta_{ij})\sim\mathcal{N}_{2}(\mu_{\alpha\beta},\Sigma_{\alpha\beta})$

with $\Sigma_{\alpha\beta}=\begin{pmatrix}\sigma^{2}&\rho\sigma^{2}\\\
\rho\sigma^{2}&\sigma^{2}\end{pmatrix}$ and $\rho=-0.9$. We have
$\mu_{\alpha\beta}=0.8\begin{pmatrix}\mu_{\alpha}\\\ \mu_{\beta}\end{pmatrix}$
to ensure a strong interspecific competition, which seemed necessary to see
possible effects of trade-off between competitive coefficients. Meta-parameter
values determining the abovementioned parameter distributions are given in
full in Appendix 3. We chose 3 sets for the distribution-level parameters
($\mu_{\pi},\sigma_{\pi},a_{\phi},b_{\phi},...$), corresponding to 3 different
levels of variance in the vital rates that allowed us to have either species
with very close vital rates (i.e., low variance of the parameter distributions
across species) or more variable vital rates (medium and high variation, see
Appendix 3).

We simulated the dynamics and computed the number of extant species after 3000
time steps, which allowed to reach an equilibrium state (see Appendix 3).
Then, for each of the three main parameter sets, we created new, replicated
parameter sets where we permuted all the inter-specific coefficients in order
to remove the correlation between $\alpha_{ij}$ and $\beta_{ij}$, while
keeping the same inter-specific competition values. We performed these
permutations 200 times and represented the histogram of the number of species
that coexist after 3000 time steps for the 200 permutated datasets. This
allowed us to compare scenarios where $\alpha_{ij}$ and $\beta_{ij}$ are
negatively correlated to the null hypothesis of no correlation between them
(which renders the emergence of coexistence through stage structure
impossible).

##  3 Results

###  3.1 Invasion analysis

First of all, the analytical expressions for the boundaries of coexistence
domains obtained thanks to the invasion analysis (eq. 16 and 17) can be
interpreted to understand how the trade-off between competition coefficients
actually promotes coexistence. In equation 16 (and similarly in equation 17),
as $n_{1a},C_{2},D_{2}$ are independent of inter-specific competition
coefficients, a decrease of $\beta_{21}$ (the effect of species 1 on species
2’s juvenile survival) on the right hand side will mechanically increase the
value of $\alpha^{*}_{21}$ on the left-hand site (the boundary value for
$\alpha_{21}$, the effect of species 1 on species 2’s fertility). What does
this mean ecologically? When the effect of species 1 on species 2’s juvenile
survival $\beta_{21}$ decreases, $\alpha_{21}$ can increase without species 2
getting extinct. This means that as one type of competition is lowered, the
other type of competition can be increased while having both species
coexisting. Thus we can have coexistence with low $\beta_{21}$ and high
$\alpha_{21}$, and reciprocally, using the other equation, with high
$\beta_{12}$ and low $\alpha_{12}$. That is, we have coexistence through
opposite competitive hierarchies, which can be obtained for the full range of
competition coefficients.

To better visualize coexistence properties implied by our invasion analysis,
we produced coexistence outcome domains for model 1. These are similar to Moll
and Brown (2008)’s but are computed analytically. They are produced using the
newly derived invasion criterion (eq. 15). Figures 1A and 1B are used for
reference (with dotted lines corresponding to the domains boundaries of figure
1A, which materializes what coexistence outcomes are predicted based on
$\alpha$ coefficients). We found, as shown in Figs. 1C, 1D and 1E, that
emergent coexistence outcomes can be observed in our model. This is best seen
in Fig. 1D, where the coexistence domain in our model combining $\alpha$ and
$\beta$ coefficients (in light grey) is extended far to the right of the
reference dashed line, well beyond what would be indicated based on $\alpha$
coefficients only (Fig. 1A). Thus we have coexistence with competition on both
vital rates when competition on a single vital rate ($\alpha$) would instead
have indicated exclusion of species 2, while $\beta$ indicates exclusion of
species 1. The three parameter sets highlighted in the Models and Methods
section are marked by asterisks, as panels C, D, E actually include a larger
span of parameter space since $(\alpha_{12},\alpha_{21})$ do vary.

In Figs. 1C and 1E, we do not find strictly emergent coexistence sensu Moll
and Brown (2008), as they restricted this label to coexistence which is not
indicated by either $\alpha$ or $\beta$ but is obtained by combining both
types of competition. We do find other emergent outcomes though (e.g. in Fig.
1E: exclusion of one of the species when $\alpha$ or $\beta$ would indicate
priority effect and coexistence, respectively).

Figure 1: Outcomes of the model predicted by invasion analysis depending on
the coefficients of the inter-specific competition on fertility ($\alpha_{12}$
and $\alpha_{21})$ with $\beta$ coefficients indicating coexistence (C),
exclusion of species 1 (D) and a priority effect (E). The dotted lines
correspond to the outcome boundaries of panel A which is used as a reference.
The three asterisks stand for the locations of the three parameter sets
corresponding to contrasted coexistence scenarios.

We also reproduced Fig. 1 for model 2 where habitats of juveniles and adults
are separated so that they cannot compete with each other, as in the
metamorphosis structure model studied by Moll and Brown (2008) (see Supplement
B for details). The results are very similar.

###  3.2 Sensitivity analysis

Sensitivity and elasticity computation at the equilibrium for the three
scenarios provides information relative to how much a small perturbation can
change the equilibrium densities, and therefore how robust coexistence is to
change in parameter values: if a small change in a parameter can make the
density of a species greatly decrease, it is not robust. With some
approximation (a linear one), it also provides information on how a large
perturbation in parameter values could generate exclusion of one species by
moving the equilibrium densities of this species to 0 (see for this the
coexistence regions in Appendix 2). We represented elasticity
$\frac{E}{n_{i}}\frac{dn_{i}}{dE}$, i.e., a relative response of the density
$n_{i}$ of species $i$ to a relative perturbation on $E$. The elasticity
values obtained from the sensitivity analysis performed on our three parameter
sets are given in Fig. 2.

Figure 2: Elasticity of the equilibrium population densities to changes in
parameters for the 3 parameter sets of Table 2.

We found here that the emergent coexistence that occurs when there are
opposite competitive hierarchies on fertility and juvenile survival (parameter
set 2) is as sensitive as the ‘classical’ (niche-based) coexistence that
occurs when the competition coefficients on both vital rates indicate
coexistence (parameter set 1). However, when coexistence stems from
competition on one vital rate indicating a priority effect and the other vital
rate indicating niche-based coexistence (parameter set 3), the resulting
equilibrium is more sensitive to perturbation, according to our elasticity
analysis. Note that sensitivity/elasticity values are formally valid only
close to the fixed point equilibrium. The coexistence regions representation
(Barabás et al. 2014 _b_ , presented in Appendix 2), which stretches a bit
that validity, confirmed the results of Fig. 2 by showing similar regions for
the ‘classical’ and the ‘emergent’ coexistence, together with a smaller domain
for the third scenario (priority effect on $\beta$ \+ coexistence on $\alpha$
$\rightarrow$ coexistence). This is likely due to the destabilizing influence
of priority effects.

###  3.3 $S$-species simulations

The species richnesses at the end of the simulations for the extended model 1
with $S$ species (see eq 8) are given in Fig. 3 (see Supplement D for results
with extended model 2). Each panel displays a histogram of the number of
species remaining in the community after 3000 time steps, for the original
parameter set including a negative correlation between $\alpha_{ij}$ and
$\beta_{ij}$, as well as for the 200 permutated parameter sets where this
correlation is set to zero.

Figure 3: Histograms of the number of species persisting at equilibrium for
200 permutations of inter-specific competition coefficients. The bar in gray
indicates the number of species persisting at time $t=3000$ time steps for the
original $\alpha_{ij}$ and $\beta_{ij}$ with
$\text{Corr}(\alpha_{ij},\beta_{ij})<0$. From left to right, columns
correspond the communities with initially 5, 10 and 40 species. From top to
bottom, the lines correspond to situations with low (A-B-C), average (D-E-F)
and high (G-H-I) variance on the vital rates.

We checked that the negative correlation between $\alpha_{ij}$ and
$\beta_{ij}$ indeed induces opposite hierarchies of competition in sets of
coefficients associated to fertility and juvenile survival (i.e., a negative
correlation between $\text{rank}(\alpha_{ij})$ and $\text{rank}(\beta_{ij})$,
see Appendix 3). In spite of this, the parameter sets with the above-mentioned
opposite competitive hierarchies do not systematically lead to richer
communities than the corresponding permutated parameter sets (where
$\text{Corr}(\alpha_{ij},\beta_{ij})=0$). For very close parameters between
species (low variances in parameter distributions, Figs. 3A, 3B and 3C),
opposite competitive hierarchies seem to have a positive effect on the species
richness, when compared to the null hypothesis of zero correlation between
$\alpha_{ij}$ and $\beta_{ij}$. However, the assumption of a community
composed solely of very similar species for all life-history parameters (as
opposed to similar averaged birth and death rates, as in neutral theory) is
rather unrealistic. To sum up, we showed that the negative correlation between
$\alpha_{ij}$ and $\beta_{ij}$, even if it indeed induces opposite competitive
hierarchies on the two vital rates, does not seem to allow, for reasonably
variable distributions of parameter values, an increase of the number of
species coexisting at the equilibrium. In other words, two-species coexistence
through differing competitive hierarchies on life-history parameters does not
scale up to many-species.

##  4 Discussion

In this article, we studied competitive interactions between species with
structured life-cycles, with multiple vital rates that can depend on
population densities. Using structured models for two competing species, Moll
and Brown (2008) described a form of coexistence whereby opposite competitive
hierarchies in the competition coefficients associated to two vital rates
could yield coexistence in the model where these two density-dependent vital
rates combine—while models with a single density-dependent vital rate would
predict exclusion of the inferior competitor. To show the existence of such
emergent coexistence, they used Ricker functions for density-dependence and a
particular interaction setup where competition between adults and juveniles is
forbidden, as in some systems with metamorphosis, where adults and juveniles
occupy separated habitats. This brought into question whether these results
could be generalized to other functional forms and life cycles. We considered
as a baseline an alternative and perhaps more frequent life cycle setup, where
adults can affect both adults and juveniles, so that both fertilities and
juvenile survival rates can depend upon adult densities of both species. We
first found, using structured two-species models, that emergent coexistence
can be generalized to Beverton-Holt functional forms and our new life cycle,
in addition to the life cycle with metamorphosis considered by Moll and Brown
(2008) (see Supplement B). We provided an in-depth explanation of how emergent
coexistence can occur, using analytical invasion criteria (facilitated by the
use of Beverton-Holt functions). Then, we explored the robustness of such
emergent coexistence equilibria with a community-level sensitivity analysis
and we found that a perturbation on parameters had not more effect on
equilibrium densities in the case of emergent coexistence than in a
‘classical’ coexistence allowed by niche separation. Finally, we numerically
explored the possibility for such mechanisms of emergent coexistence to
manifest in $S$-species models (with $S>>2$). We observed that opposite
competitive hierarchies acting on fertility vs juvenile survival, allowed by a
negative correlation between the two sets of competition coefficients, did not
seem to generate substantial increases in species richness in our simulated
communities. Importantly, this was true both for our model 1 where all
competition is generated by adults, as well as model 2 where adults only
compete with adults and juveniles with juveniles–had Moll and Brown (2008)
considered more species they would have found similar results to us.

We now dive deeper into each of the abovementioned results. As we succeeded to
reproduce the emergent coexistence of Moll and Brown (2008) with a general
invasion framework and an additional, frequently observed life cycle, we think
that coexistence emerging through opposite competitive hierarchies on
different life history traits may be a general feature of two-species two-
life-stage models. We found nothing in complex life histories with
habitat/niche shifts that specifically promotes this effect of stage structure
on coexistence. This reinforces the importance of structured models to study
interactions between species and their role in shaping community dynamics
(Miller and Rudolf, 2011), as we find that the set of life cycles where stage
structure affects coexistence is considerably extended. Furthermore, by
finding an analytic invasion criterion, we ensured that emergent coexistence
is not dependent on specific sets of vital rate parameters but rather a
general feature of two-species two-stages models, depending on competition
hierarchies. The robustness of such coexistence equilibria to perturbation of
parameters has then been confirmed by the community-level sensitivity
analysis, as the emergent coexistence scenario was not more sensitive than the
classical coexistence scenario (where coexistence was indicated by $\alpha$
and $\beta$ coefficients separately). To our knowledge, there has been little
implementation in practice of the discrete-time sensitivity analysis for
coupled nonlinear systems developed by Barabás et al. (2014 _a_). Here, it has
allowed us to explore in more depth the properties of equilibria, as we could
not easily find an expression for equilibrium densities (they are defined by
the intersection of two conic sections as shown in Supplement A, and
parameters are so interwoven in those expressions that even if analytical
solutions were to be found, these might be too complex to interpret). The
robustness of the emergent coexistence equilibria to perturbations further
suggests that such coexistence may be likely even when environmental
stochasticity is added to two-species two-stages models. It is therefore an
explanation to keep in mind when seeing coexisting pairs of species in the
wild in spite of coexistence hierarchies on a single resource axis or vital
rate indicating exclusion.

However, two-species coexistence, even with many stages, can be quite
different from many-species coexistence. Our study of $S$-species communities
did not show a positive effect of opposite competitive hierarchies in the two
sets of competition coefficients, except in communities where species were
very similar in all their parameters—but these communities are probably not
very lifelike. We have done additional explorations of many-species models,
using slightly different correlation structures, such as negative correlations
between interspecific $\alpha$ and $\beta$ coefficients within each pair of
species (see Appendix 3), rather than simply a negative correlation between
pooled $\alpha_{ij}$ and $\beta_{ij}$ values for all pairs of species. This
stronger condition on negative correlations between competitiveness on the two
vital rates had no effect on the species richness of the final community. In
sum, we showed that the properties promoting emergence of coexistence through
stage structure in the 2-species models could hardly be extended to many-
species models where all species interact together (albeit sometimes weakly).
These results resonate with those of simpler, unstructured models (Barabás et
al., 2016): criteria needed to promote $S$-species coexistence with $S\gg 2$
are much more stringent than those for two species, and require a more
important separation of niches (manifesting in the ratio between intra and
interspecific competition coefficients).

Finally, we should keep in mind that even if our results are both robust and
supported by analytical formulas, some assumptions could be debated. While we
considered iteroparous life histories, Cushing et al. (2007) considered
semelparous life histories with more variable dynamics, which creates
additional possibilities for coexistence, with 2-cycle attractors where
juveniles of species 1 co-occur with adults of species 2 and vice versa. The
extension of these semelparous life histories to our new life cycle would be
quite interesting. More generally, Cushing et al. (2007), Moll and Brown
(2008) and ourselves considered a modified Lotka-Volterra framework, with
Lotka-Volterra competition coefficients modulating the effect of adult or
juvenile densities on several vital rates. However, we have no mechanistic
derivations for such models (a derivation through MacArthur-style consumer-
resource model with separation of time scales could perhaps be attempted).
Structured mechanistic models of resource-based competition or competition for
space might inform on the realism of the $\alpha_{ij}$ and $\beta_{ij}$
correlation structures considered here (Goldberg and Landa, 1991; Loreau and
Ebenhoh, 1994; Tilman, 1982; Qi et al., 2021). That said, we expect that even
with a mechanistic perspective, our main finding will likely be robust: while
coexistence emerging from life-history complexity is quite easily seen in two-
species pairs, it is much more difficult to obtain in many-species models. In
speciose communities, some degree of niche separation is required to prevent
species to interact too much with each other. This could be promoted by many
weak-few strong interactions, possibly combined with some interaction network
structures for each stage (e.g. Kinlock, 2021) or density-dependent vital
rate, as the network structure is likely to be nonrandom and could differ
among stages.

#### Acknowledgements

We thank Sam Boireaud for an exploratory numerical study of some two-species
models, as well as Olivier Gimenez for discussions and Coralie Picoche for
comments on the manuscript. We are grateful to three reviewers for their
constructive input.

## Statements and Declarations

#### Code accessibility

Computer codes written for analyses are available at GitHub:
https://github.com/g-bardon/StructuredModels_Coexistence

#### Contributions

GB wrote the computer code, produced the figures, and performed numerical
analyses. Analytic formulas were derived by FB and GB. The first draft was
written by GB with help from FB, both authors contributed equally to
subsequent versions.

#### Competing interests and funding

The authors have no relevant financial or non-financial interests to disclose.
Funding for Gaël Bardon’s internship was provided by the French National
Research Agency through ANR Democom (ANR-16-CE02-0007) to Olivier Gimenez.

## References

  * Barabás et al. (2014 _a_) Barabás, G., G. Meszéna, and A. Ostling. 2014 _a_.  Fixed point sensitivity analysis of interacting structured populations.  Theoretical population biology 92:97–106.
  * Barabás et al. (2014 _b_) Barabás, G., L. Pásztor, G. Meszéna, and A. Ostling. 2014 _b_.  Sensitivity analysis of coexistence in ecological communities: theory and application.  Ecology Letters 17:1479–1494.
  * Barabás et al. (2016) Barabás, G., M. J. Michalska-Smith, and S. Allesina. 2016.  The effect of intra- and interspecific competition on coexistence in multispecies communities.  The American Naturalist 188:E1–E12.
  * Boutin (1984) Boutin, S. A. 1984.  The effect of conspecifics on juvenile survival and recruitment of snowshoe hares.  The Journal of Animal Ecology pages 623–637.
  * Caswell (2001) Caswell, H. 2001.  Matrix Population Models: Construction, Analysis, and Interpretation. 2nd edn Sinauer Associates.  Inc., Sunderland, MA .
  * Chesson (2000) Chesson, P. 2000.  Mechanisms of maintenance of species diversity.  Annual review of Ecology and Systematics 31:343–366.
  * Chu and Adler (2015) Chu, C., and P. B. Adler. 2015.  Large niche differences emerge at the recruitment stage to stabilize grassland coexistence.  Ecological Monographs 85:373–392.
  * Cushing (2006) Cushing, J. 2006.  Nonlinear semelparous leslie models.  Mathematical Biosciences & Engineering 3:17–36.
  * Cushing (2008) ———. 2008.  Matrix models and population dynamics.  Mathematical biology 14:47–150.
  * Cushing and Li (1989) Cushing, J., and J. Li. 1989.  On ebenman’s model for the dynamics of a population with competing juveniles and adults.  Bulletin of Mathematical Biology 51:687–713.
  * Cushing (1998) Cushing, J. M. 1998.  An introduction to structured population dynamics, volume 71.  SIAM.
  * Cushing and Henson (2012) Cushing, J. M., and S. M. Henson. 2012.  Stable bifurcations in semelparous leslie models.  Journal of biological dynamics 6:80–102.
  * Cushing et al. (2004) Cushing, J. M., S. Levarge, N. Chitnis, and S. M. Henson. 2004.  Some discrete competition models and the competitive exclusion principle.  Journal of difference Equations and Applications 10:1139–1151.
  * Cushing et al. (2007) Cushing, J. M., S. M. Henson, and L.-I. Roeger. 2007.  Coexistence of competing juvenile–adult structured populations.  Journal of Biological Dynamics 1:201–231.
  * Fujiwara et al. (2011) Fujiwara, M., G. Pfeiffer, M. Boggess, S. Day, and J. Walton. 2011.  Coexistence of competing stage-structured populations.  Scientific Reports 1:107.
  * Goldberg and Landa (1991) Goldberg, D. E., and K. Landa. 1991.  Competitive effect and response: hierarchies and correlated traits in the early stages of competition.  Journal of Ecology 79:1013–1030.
  * Haigh and Smith (1972) Haigh, J., and J. M. Smith. 1972.  Can there be more predators than prey?  Theoretical Population Biology 3:290–299.
  * Kinlock (2021) Kinlock, N. L. 2021.  Uncovering structural features that underlie coexistence in an invaded woody plant community with interaction networks at multiple life stages.  Journal of Ecology 109:384–398.
  * Leslie and Gower (1958) Leslie, P., and J. Gower. 1958.  The properties of a stochastic model for two competing species.  Biometrika 45:316–330.
  * Leslie (1945) Leslie, P. H. 1945.  On the use of matrices in certain population mathematics.  Biometrika 33:183–212.
  * Letten et al. (2017) Letten, A. D., P.-J. Ke, and T. Fukami. 2017.  Linking modern coexistence theory and contemporary niche theory.  Ecological Monographs 87:161–177.
  * Loreau and Ebenhoh (1994) Loreau, M., and W. Ebenhoh. 1994.  Competitive exclusion and coexistence of species with complex life cycles.  Theoretical Population Biology 46:58–77.
  * Meszéna et al. (1997) Meszéna, G., I. Czibula, and S. Geritz. 1997.  Adaptive dynamics in a 2-patch environment: a toy model for allopatric and parapatric speciation.  Journal of Biological Systems 5:265–284.
  * Miller and Rudolf (2011) Miller, T. E., and V. H. Rudolf. 2011.  Thinking inside the box: community-level consequences of stage-structured populations.  Trends in ecology & evolution 26:457–466.
  * Moll and Brown (2008) Moll, J. D., and J. S. Brown. 2008.  Competition and coexistence with multiple life-history stages.  The American Naturalist 171:839–843.
  * Neubert and Caswell (2000) Neubert, M. G., and H. Caswell. 2000.  Density-dependent vital rates and their population dynamic consequences.  Journal of Mathematical Biology 41:103–121.
  * Péron and Koons (2012) Péron, G., and D. N. Koons. 2012.  Integrated modeling of communities: parasitism, competition, and demographic synchrony in sympatric ducks.  Ecology 93:2456–2464.
  * Qi et al. (2021) Qi, M., N. DeMalach, T. Sun, and H. Zhang. 2021.  Coexistence under hierarchical resource exploitation: the role of r*-preemption tradeoff.  arXiv preprint arXiv:1908.08464 .
  * Tilman (1982) Tilman, D. 1982.  Resource Competition and Community Structure. (MPB-17), Volume 17.  Princeton University Press.
  * Wallace and Temple (1987) Wallace, M. P., and S. A. Temple. 1987.  Competitive interactions within and between species in a guild of avian scavengers.  The Auk 104:290–295.
  * Werner and Gilliam (1984) Werner, E. E., and J. F. Gilliam. 1984.  The ontogenetic niche and species interactions in size-structured populations.  Annual review of ecology and systematics 15:393–425.
  * Wilbur (1980) Wilbur, H. M. 1980.  Complex life cycles.  Annual review of Ecology and Systematics 11:67–93.

##  Appendix Appendix 1 Invasion analysis

We develop here the method to obtain invasion criteria in two-species two-life
stages models including competition on two vital rates, such as model 1
combining competition on fertility and juvenile survival.

We first place ourselves at the exclusion equilibrium where species 1 lives
without species 2:

$\tilde{n}=\begin{pmatrix}n_{1j}^{*}\\\ n_{1a}^{*}\\\ 0\\\ 0\\\
\end{pmatrix}.$ (A1)

To compute the equilibrium densities $n_{1j}^{*}$ and $n_{1a}^{*}$, we have to
solve the following system corresponding to the model with a single species
(we dropped here species-specific indices for simplicity):

$\left\\{\begin{array}[]{r c l}n_{j}&=&(1-\gamma)s_{j}(n)n_{j}+f(n)n_{a}\\\
n_{a}&=&\gamma s_{j}(n)n_{j}+s_{a}n_{a}\end{array}\right.$ (A2)

with fertility $f(n)$ and juvenile survival $s_{j}(n)$ affected by intra-
specific competition only, according to the relations

$f(n)=\frac{\pi}{1+\alpha n_{a}}\quad\text{and}\quad
s_{j}(n)=\frac{\phi}{1+\beta n_{a}}.$ (A3)

We now express $n_{j}$ according to $n_{a}$:

$\displaystyle n_{j}$ $\displaystyle=(1-s_{a})n_{a}\frac{1}{\gamma s_{j}(n)}$
(A4) $\displaystyle=(1-s_{a})\left(\frac{1+\beta
n_{a}}{\gamma\phi}\right)n_{a}$ (A5)

and then, as we have by assumption $n_{a}\neq 0$ and $n_{j}\neq 0$, we obtain

$\displaystyle n_{j}$ $\displaystyle=(1-\gamma)s_{j}(n)n_{j}+f(n)n_{a}$
$\displaystyle\Leftrightarrow 1$ $\displaystyle=(1-\gamma)\frac{\phi}{1+\beta
n_{a}}+\frac{\pi}{1+\alpha n_{a}}\frac{n_{a}}{n_{j}}$
$\displaystyle\Leftrightarrow 1$ $\displaystyle=(1-\gamma)\frac{\phi}{1+\beta
n_{a}}+\frac{\pi\phi\gamma}{(1-s_{a})(1+\beta n_{a})(1+\alpha n_{a})}$
$\displaystyle\Leftrightarrow(1+\beta n_{a})(1+\alpha n_{a})$
$\displaystyle=(1-\gamma)\phi(1+\alpha n_{a})+\frac{\pi\gamma\phi}{1-s_{a}}$

that leads to the polynomial of degree 2:

$\alpha\beta
n_{a}^{2}+(\alpha+\beta-\alpha(1-\gamma)\phi)n_{a}+1-(1-\gamma)\phi-\frac{\pi\gamma\phi}{1-s_{a}}=0.$
(A6)

We denote $C=(1-\gamma)\phi$ and $D=\frac{\pi\gamma\phi}{1-s_{a}}$ and we find
the solutions of the polynomial equation:

$n_{a}^{*}=\frac{(\alpha C-\alpha-\beta)\pm\sqrt{(-\alpha
C+\alpha+\beta)^{2}-4\alpha\beta(-C-D+1)}}{2\alpha\beta}.$ (A7)

We search for a positive solution and we have $C=(1-\gamma)\phi<1$ which
implies $\alpha C-\alpha-\beta<0$ assuming positive interaction coefficients
($\alpha>0$ and $\beta>0$). Positivity is then equivalent to:

$\displaystyle(\alpha C-\alpha-\beta)+\sqrt{(-\alpha
C+\alpha+\beta)^{2}-4\alpha\beta(-C-D+1)}>0$ (A8)
$\displaystyle\iff-4\alpha\beta(-C-D+1)>0$ (A9) $\displaystyle\iff C+D>1.$
(A10)

The equivalence is provided by the fact that the left term has its square
inside the square root. We define the inherent net reproductive number as
$C+D=(1-\gamma)\phi+\frac{\pi\gamma\phi}{1-s_{a}}$ that must exceed 1 to have
a viable species. We finally obtain the densities at a stable positive
equilibrium for the model with a single species:

$\left\\{\begin{array}[]{r c l}n_{j}^{*}&=&(1-s_{a})(\frac{1+\beta
n_{a}^{*}}{\gamma\phi})n_{a}^{*}\\\
n_{a}^{*}&=&\frac{(\alpha(1-\gamma)\phi-\alpha-\beta)+\sqrt{(-\alpha(1-\gamma)\phi+\alpha+\beta)^{2}-4\alpha\beta(-(1-\gamma)\phi-\frac{\pi\gamma\phi}{1-s_{a}}+1)}}{2\alpha\beta}\end{array}\right.$
(A11)

provided that

$C+D=(1-\gamma)\phi+\frac{\pi\gamma\phi}{1-s_{a}}>1.$ (A12)

We have then obtained the densities at the exclusion equilibrium. We have now
to evaluate the stability of the equilibrium $(n_{j}^{*},n_{a}^{*},0,0)$: if
it is locally asymptotically stable, the excluded species cannot invade the
community starting from very small density. Conversely, if it is unstable,
this means that the excluded species can invade provided that equilibria of
the single-species model are restricted to stable fixed points. This last
restriction is important. Indeed, if the single-species model has a stable
two-cycle or an even more complex fluctuating attractor, the fixed point will
already be unstable, and a proper invasion analysis should instead be done on
this more complex attractor. Two-cycles tend to occur for $s_{a}=0$ (Cushing
et al., 2007), as shown in Supplement C. So long as $s_{a}>0$ this model does
not show 2-cycles—so far as we have seen numerically and following the
argumentation of Cushing et al. (2007) for model 2. The combination $s_{a}=0$
and $\gamma=1$ renders the projection matrix primitive, which in turns creates
possibilities for bifurcations towards two-cycles directly from the trivial
single-species (0,0) equilibrium (Cushing and Li, 1989; Cushing, 2006). We
therefore restrict ourselve to the case $s_{a}>0$, although the equations
derived above in this Appendix for the single-species fixed point are also
valid whenever $s_{a}=0$.

To evaluate the (un)stability of the exclusion equilibria in the two-species
system, we have to compute the eigenvalues of the Jacobian of the 2 species
system evaluated at the exclusion equilibrium. The full system iterated over a
time step reads:

$\left\\{\begin{array}[]{r c
l}n_{1j}(t+1)&=&(1-\gamma_{1})\frac{\phi_{1}}{1+\beta_{11}n_{1a}(t)+\beta_{12}n_{2a}(t)}n_{1j}(t)+\frac{\pi_{1}}{1+\alpha_{11}n_{1a}(t)+\alpha_{12}n_{2a}(t)}n_{1a}(t)\\\
n_{1a}(t+1)&=&\gamma_{1}\frac{\phi_{1}}{1+\beta_{11}n_{1a}(t)+\beta_{12}n_{2a}(t)}n_{1j}(t)+s_{1a}n_{1a}(t)\\\
n_{2j}(t+1)&=&(1-\gamma_{2})\frac{\phi_{2}}{1+\beta_{22}n_{2a}(t)+\beta_{21}n_{1a}(t)}n_{2j}(t)+\frac{\pi_{2}}{1+\alpha_{22}n_{2a}(t)+\alpha_{21}n_{1a}(t)}n_{2a}(t)\\\
n_{2a}(t+1)&=&\gamma_{2}\frac{\phi_{2}}{1+\beta_{22}n_{2a}(t)+\beta_{21}n_{1a}(t)}n_{2j}(t)+s_{2a}n_{2a}(t).\end{array}\right.$
(A13)

We place ourselves at the abovementioned exclusion equilibrium where species 1
dominates the community and species 2 is absent. The Jacobian evaluated at
this point is a $4\times 4$ matrix and it has the following triangular block
form (see also Cushing 2008):

$\begin{pmatrix}B_{1}&B_{2}\\\ 0&B_{3}\\\ \end{pmatrix}.$ (A14)

Therefore, we only need to know the eigenvalues of the $2\times 2$ matrices
$B_{1}$ and $B_{3}$. The $B_{1}$ matrix corresponds to the Jacobian of the
system in the absence of species 2, so provided that the single-species
equilibrium is stable, we can infer that it has a spectral radius lower than
1. Unfortunately, although all simulations that we have done, even for very
large fertilities or low competition coefficients, provide a stable
equilibrium for $C+D>1$, it is not straightforward to show local stability
analytically. Indeed matrix $B_{1}$ writes:

$\begin{pmatrix}(1-\gamma)\frac{\phi}{1+\beta
n_{a}}&\frac{1-\gamma}{\gamma}\left(\frac{-\beta n_{a}}{1+\beta
n_{a}}\right)(1-s_{a})+\frac{\pi}{1+\alpha n_{a}}\left(1-\frac{\alpha
n_{a}}{1+\alpha n_{a}}\right)\\\ \gamma\frac{\phi}{1+\beta n_{a}}&\frac{-\beta
n_{a}}{1+\beta n_{a}}+s_{a}\end{pmatrix}$ (A15)

with $n_{a}$ given by eq. (A11), making the use of Jury conditions (used below
on the simpler $B_{3}$ block) difficult here, although the determinant
simplifies somewhat. We therefore conjecture the stability of the single-
species equilibrium, with a spectral radius of $B_{1}$ below 1, which we deem
reasonable in this context. The $B_{3}$ matrix is given by:

$B_{3}=\begin{pmatrix}(1-\gamma_{2})\frac{\phi_{2}}{1+\beta_{21}n_{1a}^{*}}&\frac{\pi_{2}}{1+\alpha_{21}n_{1a}^{*}}\\\
\gamma_{2}\frac{\phi_{2}}{1+\beta_{21}n_{1a}^{*}}&s_{2a}\\\ \end{pmatrix}.$
(A16)

The $B_{3}$ matrix block, which is part of the Jacobian, also corresponds to
the projection matrix of species 2 when invading species 1, which implies that
the exclusion equilibrium is not stable when species 2 has a growth rate
larger than 1 (i.e., the modulus of $B_{3}$’s leading eigenvalue is larger
than 1). The stability of this two-dimensional system can be investigated
thanks to the Jury conditions, given by the equations (A17, A18, A19):

$\displaystyle 1-\operatorname{tr}J(\mathbf{n})+\det J(\mathbf{n})$
$\displaystyle>0$ (A17) $\displaystyle 1+\operatorname{tr}J(\mathbf{n})+\det
J(\mathbf{n})$ $\displaystyle>0$ (A18) $\displaystyle 1-\det J(\mathbf{n})$
$\displaystyle>0$ (A19)

with $J(\mathbf{n})$ the Jacobian matrix of the two-compartment system whose
stability is being evaluated. If these three conditions hold, the system is
stable. If the first condition (equation (A17)) is violated, one of the
eigenvalues of $J(\mathbf{n})$ is larger than 1, which means in our context
that the excluded species can invade the community. The violations of the
conditions (A18) and (A19) correspond to the creation of limit cycles or
invariant loops (Neubert and Caswell, 2000). We then check these conditions on
the $B_{3}$ matrix. We have:

$\displaystyle\operatorname{tr}B_{3}$
$\displaystyle=(1-\gamma_{2})\frac{\phi_{2}}{1+\beta_{21}n_{1a}^{*}}+s_{2a}$
$\displaystyle\det B_{3}$
$\displaystyle=(1-\gamma_{2})\frac{\phi_{2}}{1+\beta_{21}n_{1a}^{*}}s_{2a}-\gamma_{2}\frac{\phi_{2}}{1+\beta_{21}n_{1a}^{*}}\frac{\pi_{2}}{1+\alpha_{21}n_{1a}^{*}}.$

We start with the third condition (A19):

$1-(1-\gamma_{2})\frac{\phi_{2}}{1+\beta_{21}n_{1a}^{*}}s_{2a}+\gamma_{2}\frac{\phi_{2}}{1+\beta_{21}n_{1a}^{*}}\frac{\pi_{2}}{1+\alpha_{21}n_{1a}^{*}}>0.$
(A20)

This equation will always be satisfied with biologically meaningful parameters
and strictly competitive interaction between species:

  * •

The third term of (A20) is positive because each parameter is positive,

  * •

$(1-\gamma_{2})\frac{\phi_{2}}{1+\beta_{21}n_{1a}^{*}}s_{2a}<1$ because
$0\leq\gamma\leq 1$, $0\leq\phi\leq 1$ and all $\beta$’s are positive.

Since $\operatorname{tr}B_{3}>0$, the first Jury’s condition implies the
second ((S29) $\Rightarrow$ (S30)). Then, we only have to check the first Jury
condition to describe the stability of the system:

$\displaystyle
1-(1-\gamma_{2})\frac{\phi_{2}}{1+\beta_{21}n_{1a}^{*}}+s_{2a}+(1-\gamma_{2})\frac{\phi_{2}}{1+\beta_{21}n_{1a}^{*}}s_{2a}-\gamma_{2}\frac{\phi_{2}}{1+\beta_{21}n_{1a}^{*}}\frac{\pi_{2}}{1+\alpha_{21}n_{1a}^{*}}>0$
$\displaystyle(1-\gamma_{2})\frac{\phi_{2}}{1+\beta_{21}n_{1a}^{*}}-(1-\gamma_{2})\frac{\phi_{2}}{1+\beta_{21}n_{1a}^{*}}s_{2a}+\gamma_{2}\frac{\phi_{2}}{1+\beta_{21}n_{1a}^{*}}\frac{\pi_{2}}{1+\alpha_{21}n_{1a}^{*}}<1-s_{2a}$
$\displaystyle(1-\gamma_{2})\phi_{2}(1-s_{2a})<(1-s_{2a})(1+\beta_{21}n_{1a}^{*})-\frac{\gamma_{2}\phi_{2}\pi_{2}}{1+\alpha_{21}n_{1a}^{*}}$
$\displaystyle(1-\gamma_{2})\phi_{2}\frac{1}{1+\beta_{21}n_{1a}^{*}}+\frac{\gamma_{2}\phi_{2}\pi_{2}}{1-s_{2a}}\frac{1}{(1+\beta_{21}n_{1a}^{*})(1+\alpha_{21}n_{1a}^{*})}<1.$

We finally have the following condition for the stability of the exclusion
equilibrium:

$\frac{C_{2}}{1+\beta_{21}n_{1a}^{*}}+\frac{D_{2}}{(1+\beta_{21}n_{1a}^{*})(1+\alpha_{21}n_{1a}^{*})}<1$
(A21)

with $C_{2}=(1-\gamma_{2})\phi_{2}$ and
$D_{2}=\frac{\pi_{2}\gamma_{2}\phi_{2}}{1-s_{2a}}$.

This last expression (eq. 15) gives us an invasion criteria that is larger
than 1 if, conversely, species 2 is able to invade the community when rare.

Note that if we remove the competition on one vital rate by setting all the
competition coefficients associated to this vital rate to 0, we re-obtain the
invasion criteria given by Fujiwara et al. (2011) for their model where a
single vital rate was affected by competition. This ensures some internal
coherence to the analytical results.

##  Appendix Appendix 2 Coexistence regions obtained through sensitivity
analysis for the three scenarios of coexistence

To have another, complementary representation of the sensitivities, we
estimated coexistence regions for the three parameter sets using the
definition given by Barabás et al. (2014 _b_). Note that we arbitrarily use
the term region in place of domain to separate the regions found with this
sensitivity analysis from the domains of coexistence found through invasion
analysis. The method consists in using the sensitivities at the equilibrium
for each parameter in order to determine the region of parameter space where
species both persists despite changes in parameter values. It is calculated by
finding the smallest (positive and negative) perturbation that would lead to
the extinction of one of the species from the value of sensitivities and
densities at equilibrium, and makes the (crude) approximation that
sensitivities are not only valid close to the fixed point but also far away
from it. The coexistence regions for the three parameter sets are presented in
Fig. A1.

Figure A1: Parameter values and coexistence regions of the three parameter
sets of Table 2 (A - set 1, B - set 2, C - set 3) chosen to promote different
scenarios of coexistence.

We see here smaller coexistence regions in the third scenario, consistent with
this equilibrium being more sensitive (less robust) to perturbations.

##  Appendix Appendix 3 Additional simulations with $S$ species, $S>2$

The meta-parameter sets, used in the parameter distributions from which the
vital rate parameters for each species have been drawn, are provided in Table
A1.

Distribution variance | $\mu_{\pi_{i}}$ | $\sigma_{\pi_{i}}$ | a | b | $\mu_{\alpha_{ii}}$ | $\sigma_{\alpha_{ii}}$ | $\mu_{\beta_{ii}}$ | $\sigma_{\beta_{ii}}$ | $\mu_{ij}$ | $\sigma^{2}$
---|---|---|---|---|---|---|---|---|---|---
high | 3 | 0.5 | 5 | 5 | 0.05 | 0.01 | 0.05 | 0.01 | 0.04 | 0.02
average | 3 | 0.05 | 50 | 50 | 0.05 | 0.01 | 0.05 | 0.01 | 0.04 | 0.02
small | 3 | 0.005 | 500 | 500 | 0.05 | 0.001 | 0.05 | 0.001 | 0.04 | 0.02
Table A1: (Meta)-Parameter sets used to draw parameter values to simulate
community dynamics with $(a,b)$ parameters of the Beta distribution for
$\phi$, $s_{a}$ and $\gamma$ ($a=a_{\phi_{i}}=a_{s_{a,i}}=a_{\gamma_{i}}$ and
$b=b_{\phi_{i}}=b_{s_{a,i}}=b_{\gamma_{i}}$). The last two columns correspond
respectively to the mean of inter-specific coefficients (for both
$\alpha_{ij}$ and $\beta_{ij}$) and to the variance of the inter-specific
coefficients (with the covariance between $\alpha_{ij}$ and $\beta_{ij}$ equal
to $-0.9\times\sigma^{2}$).

For initially 40 species, a small variance of parameters across species, and
opposite competitive hierarchies, we obtained the dynamics illustrated in Fig.
A2.

Figure A2: Population dynamics of a community with initially 40 species and
opposite competitive hierarchies on fertility and juvenile survival. All
species have very similar vital rates due to a low variance in the
distribution used to draw the parameter. Each line corresponds to one species.

The population dynamics highlighted in Fig. A2 demonstrates that 3000 time
steps are enough to reach a stable number of coexisting species.

We verified our hypothesis that a negative correlation between $\alpha_{ij}$
and $\beta_{ij}$ lead to opposite hierarchies in sets of coefficients
associated to fertility and juvenile survival by plotting the rank of
$\alpha_{ij}$ against the rank of $\beta_{ij}$. For 40 species, these ranks
are plotted in Fig. A3.

Figure A3: Rank of $\beta_{ij}$ versus rank of $\alpha_{ij}$ under a negative
correlation of $\beta_{ij}$ and $\alpha_{ij}$ of $-0.9$.

The negative correlation between ranks in Fig. A3 confirmed that the negative
correlation between $\alpha_{ij}$ and $\beta_{ij}$ was creating the expected
opposite competitive hierarchies.

However, opposite competitive hierarchies may not be sufficient to promote
reciprocal exclusion (e.g., exclusion of 1 by 2 in the fertility competition
model and of 2 by 1 in the juvenile survival competition model) of species
within each or most pairs of species. To make the comparison between 2 and
S-species model, we have to check that we consider situations that are fully
comparable, i.e., where species would exclude each other out in a two-species
contest on either $\alpha$ or $\beta$ competition. We therefore checked if the
trade-off between being competitive on $\alpha$ and $\beta$ could promote such
reciprocal exclusion or priority effects (in models with a single parameter
that is density dependent) within pairs of species. We computed the invasion
criteria of such simple models (where only one vital rate is affected by
competition) for each pair of species composing the communities. For a
community of initially 40 species and a small variance in parameters across
species, we found the coexistence outcomes in simpler models highlighted in
Fig. A4.

We found that for a community of 40 species, the proportions of pairs of
species where a priority effect or a reciprocal exclusion are indicated (by
the simpler models with density-dependence on a single vital rate) are low
when we assume only a negative correlation between $\alpha_{ij}$ and
$\beta_{ij}$.

Figure A4: Frequencies of outcomes predicted by pairwise invasion criteria for
a parameter set with 40 species with very close parameters and a negative
correlation between $\alpha_{ij}$ and $\beta_{ij}$ with $\phi=0.9\sigma$

In fact, when we simply use a negative correlation to create opposite
competitive hierarchies, we draw our couple of values ($\alpha_{ij}$,
$\beta_{ij}$) in a Gaussian distribution, and then most of the values are
close to their means, which does not allow to obtain pairs of species where
each species is strongly competitive on a different vital rate, and exclude
the other if the competition was on this vital rate only. Therefore, even if
there are opposite competitive hierarchies within pairs of species modelled,
this setup is unlikely to generate situations where species exclude each other
when considering single-vital-rate-competition.

To consider S-species scenarios where species can exclude each other out if
competition was solely on $\alpha$ or solely on $\beta$ (hereafter referred to
as ‘reciprocal exclusion’), we changed our structure of correlation to
generate situations where within most pairs of species, each species is
strongly competitive on a different vital rate (e.g. species 1 excludes
species 2 in the model with only fertility competition and vice-versa in the
model with only juvenile survival competition).

Our method to achieve reciprocal exclusion of pairs (for single-vital-rate-
competition) consisted in drawing pairs $(a_{ij},b_{ij})$ from a bivariate
normal distribution with $\mu_{a_{ij}}$<$\mu_{b_{ij}}$ and
$\text{Corr}(a_{ij},b_{ij})<0$. We then always assign $(a_{ij},b_{ij})$ to
$(\alpha_{ij},\beta_{ij})$ and $(b_{ji},a_{ji})$ to
$(\alpha_{ji},\beta_{ji})$. In this way, for each pair of species $(i,j)$,
species $i$ is competitive on $\alpha$ and species $j$ is competitive on
$\beta$. We acknowledge that this may be hard to grasp, and suggest to the
reader to use the code to create an example.

We computed the invasion criteria for each pair of species and counted how
often each invasion scenario in single-vital-rate-competition models, which we
represented in the histogram of Fig. A5.

Figure A5: Frequencies of outcomes predicted by pairwise invasion criteria -
for singe-vital-rate-competition models - for a parameter set with 40 species
with very close parameters and bimodal distribution of interspecific
coefficients drawn to promote situations of opposite exclusion within each
pair of species.

We succeeded to generate reciprocal exclusion within most pairs of species.
Moreover, we showed that a permutation of competition coefficient was
sufficient to remove the structure of reciprocal exclusion (Fig. A6).

Figure A6: Frequencies of outcomes predicted by pairwise invasion criteria for
a parameter set with 40 species, with very close parameters and permuted
competition coefficient that promoted before permutation situations of
reciprocal exclusion within each pair of species.

However, as for our simple negative correlation between competitive ranks for
$\alpha$ and $\beta$, this new parameter set where reciprocal exclusion occurs
within most pairs of species does not allow to significantly increase the
number of extant species at $t=3000$ time steps (Fig. A7).

Figure A7: Histograms of the number of species persisting at equilibrium, with
initially 40 species, and for 200 permutations of inter-specific competition
coefficients. The bar in gray indicates the number of species persisting at
equilibrium for the original $\alpha_{ij}$ and $\beta_{ij}$ with
$\text{Corr}(\alpha_{ij},\beta_{ij})<0$ and reciprocal exclusion for each pair
of species.

Therefore, even if competition is highly structured in a way that makes
species exclude each other out when a single vital rate is density-dependent,
which should greatly promote situations of emergent coexistence, we did not
observe a positive effect of this structure on species richness in a community
larger than $S=2$. We therefore conclude that emergent coexistence is highly
unlikely in many-species contexts.

##  Appendix Supplement A Fixed points for the two-species models: equations
and analytical progress in a simple case

To compute the densities of both species at equilibrium for model 1, we have
to solve the following system:

$\left\\{\begin{array}[]{r c
l}n_{1j}&=&(1-\gamma_{1})\frac{\phi_{1}}{1+\beta_{11}n_{1a}+\beta_{12}n_{2a}}n_{1j}+\frac{\pi_{1}}{1+\alpha_{11}n_{1a}+\alpha_{12}n_{2a}}n_{1a}\\\
n_{1a}&=&\gamma_{1}\frac{\phi_{1}}{1+\beta_{11}n_{1a}+\beta_{12}n_{2a}}n_{1j}+s_{1a}n_{1a}\\\
n_{2j}&=&(1-\gamma_{2})\frac{\phi_{2}}{1+\beta_{22}n_{2a}+\beta_{21}n_{1a}}n_{2j}+\frac{\pi_{2}}{1+\alpha_{22}n_{2a}+\alpha_{21}n_{1a}}n_{2a}\\\
n_{2a}&=&\gamma_{2}\frac{\phi_{2}}{1+\beta_{22}n_{2a}+\beta_{21}n_{1a}}n_{2j}+s_{2a}n_{2a}.\end{array}\right.$
(S1)

We start by expressing $n_{1j}$ according to $n_{1a}$:

$\displaystyle
n_{1j}\left(1-\frac{(1-\gamma_{1})\phi_{1}}{1+\beta_{11}n_{1a}+\beta_{12}n_{2a}}\right)=\frac{\pi_{1}}{1+\alpha_{11}n_{1a}+\alpha_{12}n_{2a}}$
$\displaystyle
n_{1j}=\frac{(1+\beta_{11}n_{1a}+\beta_{12}n_{2a})\pi_{1}}{(1+\alpha_{11}n_{1a}+\alpha_{12}n_{2a})(1+\beta_{11}n_{1a}+\beta_{12}n_{2a}-(1-\gamma_{1})\phi_{1}}n_{1a}.$

Then, we insert the latter expression in the expression for $n_{1a}$:

$\displaystyle
n_{1a}=\frac{\gamma_{1}\phi_{1}}{1+\beta_{11}n_{1a}+\beta_{12}n_{2a}}\frac{(1+\beta_{11}n_{1a}+\beta_{12}n_{2a})\pi_{1}}{(1+\alpha_{11}n_{1a}+\alpha_{12}n_{2a})(1+\beta_{11}n_{1a}+\beta_{12}n_{2a}-(1-\gamma_{1})\phi_{1}}n_{1a}+s_{1a}n_{1a}.$

We place ourselves in the case where both species coexist. Therefore,
$n_{1a}\neq 0$ and we have:

$\displaystyle(1+\alpha_{11}n_{1a}+\alpha_{12}n_{2a})(1+\beta_{11}n_{1a}+\beta_{12}n_{2a}-(1-\gamma_{1})\phi_{1})=\frac{\gamma_{1}\phi_{1}\pi_{1}}{1-s_{1a}}$
$\displaystyle(1+\alpha_{11}n_{1a}+\alpha_{12}n_{2a})(1+\beta_{11}n_{1a}+\beta_{12}n_{2a}-C_{1})=D_{1}$
$\displaystyle\frac{C_{1}}{1+\beta_{11}n_{1a}+\beta_{12}n_{2a}}+\frac{D_{1}}{(1+\beta_{11}n_{1a}+\beta_{12}n_{2a})(1+\alpha_{11}n_{1a}+\alpha_{12}n_{2a})}=1.$

This last expression correspond to the growth rate of the species 1 that must
be equal to 1 at the equilibrium.
By symmetry, we have the same equation for the other species :

$\displaystyle(1+\alpha_{22}n_{2a}+\alpha_{21}n_{1a})(1+\beta_{22}n_{2a}+\beta_{21}n_{1a}-C_{2})=D_{2}$
$\displaystyle\frac{C_{2}}{1+\beta_{22}n_{2a}+\beta_{21}n_{1a}}+\frac{D_{2}}{(1+\beta_{22}n_{2a}+\beta_{21}n_{1a})(1+\alpha_{22}n_{2a}+\alpha_{21}n_{1a})}=1.$

Finally, to compute the equilibrium densities, we have to solve the following
system for $n_{1a}$ and $n_{2a}$ :

$\left\\{\begin{array}[]{r c
l}(1+\alpha_{11}n_{1a}+\alpha_{12}n_{2a})(1+\beta_{11}n_{1a}+\beta_{12}n_{2a}-C_{1})=D_{1}\\\
(1+\alpha_{22}n_{2a}+\alpha_{21}n_{1a})(1+\beta_{22}n_{2a}+\beta_{21}n_{1a}-C_{2})=D_{2}\end{array}\right.$
(S2)

which can also be written in the following form with $x=n_{1a}$ and $y=n_{2a}$
:

$\left\\{\begin{array}[]{r c
l}\alpha_{11}\beta_{11}x^{2}+\alpha_{12}\beta_{12}y^{2}+(\alpha_{11}\beta_{12}+\alpha_{12}\beta_{11})xy+(\beta_{11}+\alpha_{11}-C_{1}\alpha_{11})x+(\beta_{12}+\alpha_{12}-C_{1}\alpha_{12})y+(1-C_{1}-D_{1})=0\\\
\alpha_{22}\beta_{22}y^{2}+\alpha_{21}\beta_{21}x^{2}+(\alpha_{22}\beta_{21}+\alpha_{21}\beta_{22})xy+(\beta_{22}+\alpha_{22}-C_{2}\alpha_{22})y+(\beta_{21}+\alpha_{21}-C_{2}\alpha_{21})x+(1-C_{2}-D_{2})=0.\end{array}\right.$
(S3)

This system corresponds to two bivariate quadratics equations (quadrics).
Therefore, the solutions of the system are the intersection points of the two
conic sections described by each equation. These conic sections are themselves
the intersection of the surface of a cone with a plane. We can classify the
conic sections by looking at the discriminant of the matrix representing the
conic sections.

With the standard notation of a quadratic equation:

$Ax^{2}+Bxy+Cy^{2}+Dx+Ey+F=0.$ (S4)

We use the usual matrix representation of the conic sections that have the
following general form:

$\mathbf{x^{T}}A_{Q}\mathbf{x}$ (S5)

with

$\mathbf{x}=\begin{pmatrix}x\\\ y\\\ 1\end{pmatrix}$ (S6)

and

$A_{Q}=\begin{pmatrix}A&B/2&D/2\\\ B/2&C&E/2\\\ D/2&E/2&F\end{pmatrix}$ (S7)

the matrix of the quadratic equation. Then, we denote the matrix of the
quadratic form $A_{33}$:

$A_{33}=\begin{pmatrix}A&B/2\\\ B/2&C\end{pmatrix}$ (S8)

which is obtained by removing the 3rd column and the 3rd line to the $A_{Q}$
matrix.

With these standard notations, we know that the conic section is degenerated
if $\det A_{Q}=0$ and is non-degenerated otherwise. Then, if the conic section
is non-degenerated, we have the following classification :

  * •

the section is a hyperbola if and only if $\det A_{33}<0$

  * •

the section is a parabola if and only if $\det A_{33}=0$

  * •

the section is an ellipse if and only if $\det A_{33}>0$

If the conic section is degenerated ($\det A_{Q}=0$), the conic section
corresponds to:

  * •

2 intersecting lines if and only if $\det A_{33}<0$

  * •

2 parallel straight lines if and only if $\det A_{33}=0$

  * •

A single point if and only if $\det A_{33}>0$

Firstly, we look at the determinant of the $A_{Q1}$ matrix for the first
equation of the system, corresponding to the first species:

$A_{Q1}=\begin{pmatrix}\alpha_{11}\beta_{11}&\frac{\alpha_{11}\beta_{12}+\alpha_{12}\beta_{11}}{2}&\frac{\beta_{11}+\alpha_{11}(1-C_{1})}{2}\\\
\frac{\alpha_{11}\beta_{12}+\alpha_{12}\beta_{11}}{2}&\alpha_{12}\beta_{12}&\frac{\beta_{12}+\alpha_{12}(1-C_{1})}{2}\\\
\frac{\beta_{11}+\alpha_{11}(1-C_{1})}{2}&\frac{\beta_{12}+\alpha_{12}(1-C_{1})}{2}&1-C_{1}-D_{1}\end{pmatrix}.$
(S9)

We can easily find an expression of the determinant of the matrix:

$\det A_{Q1}=(\alpha_{11}\beta_{12}-\alpha_{12}\beta_{11})D_{1}.$ (S10)

We have $D_{1}\neq 0$, therefore we have $\det
A_{Q1}=0\Leftrightarrow\alpha_{11}\beta_{12}=\alpha_{12}\beta_{11}$.

We place ourselves in the case of
$\alpha_{11}\beta_{12}=\alpha_{12}\beta_{11}$. For the first equation of the
system, we have :

$A_{33,1}=\begin{pmatrix}\alpha_{11}\beta_{11}&\frac{\alpha_{11}\beta_{12}+\alpha_{12}\beta_{11}}{2}\\\
\frac{\alpha_{11}\beta_{12}+\alpha_{12}\beta_{11}}{2}&\alpha_{12}\beta_{12}\end{pmatrix}.$
(S11)

We can see that the determinant of the matrix $A_{33,1}$ is equal to 0
provided $\alpha_{11}\beta_{12}=\alpha_{12}\beta_{11}$. Then we find that the
conic section corresponds to 2 parallel straight lines (a degenerate parabola)
when $\alpha_{11}\beta_{12}=\alpha_{12}\beta_{11}$.

In the case where $\alpha_{11}\beta_{12}\neq\alpha_{12}\beta_{11}$, the conic
section is non-degenerated, we have to look at the sign of $\det A_{33}$. For
the first equation, we have:

$\displaystyle\det A_{33,1}$
$\displaystyle=\alpha_{11}\beta_{11}\alpha_{12}\beta_{12}-(\frac{\alpha_{11}\beta_{12}+\alpha_{12}\beta_{11}}{2})^{2}$
$\displaystyle=\frac{1}{4}(\alpha_{11}\beta_{12}-\alpha_{12}\beta_{11})^{2}.$

The discriminant is then always positive, which indicates that the conic
sections are hyperbola. As the system is symmetric, we have the same
conclusions for species 2.

The analytical computation of the intersection of the conic sections, in the
general case, requires a tedious analysis that consists in degenerating one
conic section to find more easily the expression of the intersection. Notably,
it requires solving of a third degree equation that would lead to analytical
results which are probably so complex that they might be impossible to
interpret. We have not found a way to prove the uniqueness of a strictly
positive intersection point but our graphical observations on multiple
parameter sets seemed to show that we can find only one or zero positive point
of intersection.

In order to go further on the analysis of the intersection point, we placed
ourselves in the case where $\alpha=\beta$. This case corresponds to a model
where adult and juvenile parameters react in the same way to increases in
adult densities. Even if it can be biologically inaccurate, it provides us a
simple framework to understand how the conic sections can intersect. With
$\alpha=\beta$, we fulfill the condition to have degenerated conic sections,
therefore the computation of the intersection points is easier. The equation
of the degenerated conic section for species 1 becomes in this case:

$\alpha_{11}^{2}x^{2}+2\alpha_{11}\alpha_{12}xy+\alpha_{12}^{2}y^{2}+\alpha_{11}(2-C_{1})x+\alpha_{12}(2-C_{1})y+(1-C_{1}-D_{1})=0.$
(S12)

By separating the terms in $x$ and $y$, we can find the following expressions
:

$y=-\frac{\alpha_{11}}{\alpha_{12}}x+\frac{\alpha_{11}}{\alpha_{12}}\left(\sqrt{\frac{1}{4}C_{1}^{2}+D_{1}}-1+\frac{C_{1}}{2}\right)$
(S13)

and

$y=-\frac{\alpha_{11}}{\alpha_{12}}x+\frac{\alpha_{11}}{\alpha_{12}}\left(-\sqrt{\frac{1}{4}C_{1}^{2}+D_{1}}-1+\frac{C_{1}}{2}\right)$
(S14)

which correspond to two parallel lines. We have symmetrical equations for
species 2, which allow us to compute the intersections between the two pairs
of lines. We find that for both species, only one of the two lines cross the
section where $x>0$ and $y>0$. Therefore, there is at maximum one intersection
point which has both coordinates positive. This point has for coordinates:

$\displaystyle x$
$\displaystyle=\frac{\alpha_{12}\alpha_{21}}{\alpha_{22}\alpha_{12}-\alpha_{11}\alpha_{21}}\left(\frac{\alpha_{22}}{\alpha_{21}}\left(\sqrt{\frac{1}{4}C_{2}^{2}+D_{2}}-1+\frac{C_{2}}{2}\right)-\frac{\alpha_{11}}{\alpha_{12}}\left(\sqrt{\frac{1}{4}C_{1}^{2}+D_{1}}-1+\frac{C_{1}}{2}\right)\right)$
(S15) $\displaystyle y$
$\displaystyle=\frac{\alpha_{11}\alpha_{22}}{\alpha_{22}\alpha_{12}-\alpha_{11}\alpha_{21}}\left(\sqrt{\frac{1}{4}C_{1}^{2}+D_{1}}+\frac{C_{1}}{2}-\sqrt{\frac{1}{4}C_{2}^{2}+D_{2}}-\frac{C_{2}}{2}\right).$
(S16)

We can search the conditions to have $x>0$ and $y>0$.
We assume $\alpha_{22}\alpha_{12}>\alpha_{11}\alpha_{21}$.
The condition to have $y>0$ is:

$\sqrt{\frac{1}{4}C_{1}^{2}+D_{1}}+\frac{C_{1}}{2}>\sqrt{\frac{1}{4}C_{2}^{2}+D_{2}}+\frac{C_{2}}{2}$
(S17)

and to have $x>0$, we need:

$\displaystyle\frac{\alpha_{22}}{\alpha_{21}}\left(\sqrt{\frac{1}{4}C_{2}^{2}+D_{2}}-1+\frac{C_{2}}{2}\right)>\frac{\alpha_{11}}{\alpha_{12}}\left(\sqrt{\frac{1}{4}C_{1}^{2}+D_{1}}-1+\frac{C_{1}}{2}\right)$
(S18)
$\displaystyle\alpha_{22}\alpha_{12}\left(\sqrt{\frac{1}{4}C_{2}^{2}+D_{2}}-1+\frac{C_{2}}{2}\right)>\alpha_{11}\alpha_{21}\left(\sqrt{\frac{1}{4}C_{1}^{2}+D_{1}}-1+\frac{C_{1}}{2}\right)$
(S19)

which is valid provided the condition to have $y>0$.

Therefore, we obtained the conditions to have a fixed point whose coordinates
are both positive, which correspond to the conditions of coexistence, for the
model where $\alpha=\beta$.

The search for fixed points in our model combining density-dependence on
fertility and juvenile survival have allowed us to better understand how a
strictly positive equilibrium (interior solution) can occur. With the simpler
model where $\alpha=\beta$, we found additional conditions for coexistence
than with the invasion analysis. However, the analytical results remain
difficult to interpret and/or restricted to the simplest models.

##  Appendix Supplement B Invasion analysis for a model with metamorphosis
life history structure

We performed our invasion analysis on the model with a life-history structure
similar to Moll and Brown (2008) and Cushing et al. (2007), where juveniles
and adults are separated, i.e. juveniles are in competition with juveniles and
adults with adults. The projection matrix of this model is given in equation
(S36). The exact same method as in section 2.3 was used to compute the
invasion criteria with this life-history structure.

In the metamorphosis model, with our notations and Beverton-Holt competition
functions, the model can be written, for species 1:

$A_{1}^{M}=\begin{pmatrix}(1-\gamma_{1})\frac{\phi_{1}}{1+\beta_{11}n_{1j}+\beta_{12}n_{2j}}&\frac{\pi_{1}}{1+\alpha_{11}n_{1a}+\alpha_{12}n_{2a}}\\\
\gamma_{1}\frac{\phi_{1}}{1+\beta_{11}n_{1j}+\beta_{12}n_{2j}}&s_{1a}\end{pmatrix}.$
(S20)

We applied the invasion analysis as for our combined model. Firstly, we search
the densities of exclusion equilibria. Then, we have the following system for
the model with a single species (we drop here the time indices for clarity and
because the model is evaluated at the fixed point):

$\left\\{\begin{array}[]{r c l}n_{j}&=&(1-\gamma)\frac{\phi}{1+\beta
n_{j}}n_{j}+\frac{\pi}{1+\alpha n_{a}}n_{a}\\\
n_{a}&=&\gamma\frac{\phi}{1+\beta n_{j}}n_{j}+s_{a}n_{a}.\end{array}\right.$
(S21)

We can easily express $n_{j}$ according to $n_{a}$:

$\displaystyle n_{a}=\frac{\gamma\phi}{1+\beta n_{j}}n_{j}+s_{a}n_{a}$
$\displaystyle\Leftrightarrow\gamma\phi
n_{j}=n_{a}(1-s_{a})+n_{a}(1-s_{a})\beta n_{j}$
$\displaystyle\Leftrightarrow(\gamma\phi-
n_{a}(1-s_{a})\beta)n_{j}=n_{a}(1-s_{a})$

so that

$n_{j}=\frac{(1-s_{a})n_{a}}{\gamma\phi-(1-s_{a})\beta n_{a}}.$ (S22)

We replace $n_{j}$ by its expression in the first equation of the system
(S21), that leads to the polynomial of degree 2:

$E\alpha\beta n_{a}^{2}+\left(\alpha(1-C)+\beta(\pi+E)\right)n_{a}+1-C-D$
(S23)

where $C=\phi(1-\gamma)$, $D=\frac{\pi\gamma\phi}{1-s_{a}}$ and
$E=\frac{(1-\gamma)(1-s_{a})}{\gamma}$.

We find the same condition for inherent net reproductive number than for our
model 1. It is given by $C+D>1$ that must be true to have a viable species.
This is due to the similarity of the equations to those of model 1, which
mainly differ in that factor $E$ now appears in front of the degree 2 term. We
then have positive densities for the exclusion equilibrium where one species
is absent:

$\left\\{\begin{array}[]{r c
l}n_{j}^{*}&=&\frac{(1-s_{a})n_{a}^{*}}{\gamma\phi-(1-s_{a})\beta
n_{a}^{*}}\\\
n_{a}^{*}&=&\frac{-\left(\alpha(1-C)+\beta(\pi+E)\right)+\sqrt{\left(\alpha(1-C)+\beta(\pi+E)\right)^{2}-4E\alpha\beta(1-C-D)}}{2E\alpha\beta}\end{array}\right.$
(S24)

provided

$C+D=(1-\gamma)\phi+\frac{\pi\gamma\phi}{1-s_{a}}>1.$ (S25)

We have then obtained the densities at the exclusion equilibrium. We have now
to evaluate the stability of this equilibrium: if it is locally asymptotically
stable, the excluded species cannot invade the community starting from very
small density. Conversely, if it is unstable, this means that the excluded
species can invade provided that equilibria of the single-species model are
restricted to stable fixed points. This last restriction is key. Indeed, if
the single-species model has a stable two-cycle or an even more complex
fluctuating attractor, the fixed point will already be unstable, and a proper
invasion analysis should instead be done on this more complex attractor. Two-
cycles tend to occur for $s_{a}=0$ (Cushing et al., 2007), as shown in
Supplement C. So long as $s_{a}>0$ this model does not show 2-cycles—so far as
we have seen numerically and following the argumentation of Cushing et al.
(2007) for model 2. We therefore restrict ourselves to the case $s_{a}>0$,
although the equations derived above in this Appendix for the single-species
fixed point are also valid whenever $s_{a}=0$.

To evaluate the (un)stability of the exclusion equilibria in the two-species
system, we have to compute the eigenvalues of the Jacobian of the 2 species
system evaluated at the exclusion equilibrium. The full system iterated over a
time step reads:

$\left\\{\begin{array}[]{r c
l}n_{1j}(t+1)&=&(1-\gamma_{1})\frac{\phi_{1}}{1+\beta_{11}n_{1j}(t)+\beta_{12}n_{2j}(t)}n_{1j}(t)+\frac{\pi_{1}}{1+\alpha_{11}n_{1a}(t)+\alpha_{12}n_{2a}(t)}n_{1a}(t)\\\
n_{1a}(t+1)&=&\gamma_{1}\frac{\phi_{1}}{1+\beta_{11}n_{1j}(t)+\beta_{12}n_{2j}(t)}n_{1j}(t)+s_{1a}n_{1a}(t)\\\
n_{2j}(t+1)&=&(1-\gamma_{2})\frac{\phi_{2}}{1+\beta_{22}n_{2j}(t)+\beta_{21}n_{1j}(t)}n_{2j}(t)+\frac{\pi_{2}}{1+\alpha_{22}n_{2a}(t)+\alpha_{21}n_{1a}(t)}n_{2a}(t)\\\
n_{2a}(t+1)&=&\gamma_{2}\frac{\phi_{2}}{1+\beta_{22}n_{2j}(t)+\beta_{21}n_{1j}(t)}n_{2j}(t)+s_{2a}n_{2a}(t).\end{array}\right.$
(S26)

We place ourselves at the abovementioned exclusion equilibrium where species 1
dominates the community and species 2 is absent. The Jacobian evaluated at
this point is a $4\times 4$ matrix and it has the following triangular block
form:

$\begin{pmatrix}B_{1}&B_{2}\\\ 0&B_{3}\\\ \end{pmatrix}.$ (S27)

Therefore, we only need to know the eigenvalues of the $2\times 2$ matrices
$B_{1}$ and $B_{3}$. The $B_{1}$ matrix corresponds to the Jacobian of the
single-species dynamical system, which has always been found to be stable
numerically, thus for reasons explained in more detail in Appendix 1 we assume
the eigenvalues have a modulus smaller than 1. The $B_{3}$ matrix is given by:

$B_{3}=\begin{pmatrix}(1-\gamma_{2})\frac{\phi_{2}}{1+\beta_{21}n_{1j}^{*}}&\frac{\pi_{2}}{1+\alpha_{21}n_{1a}^{*}}\\\
\gamma_{2}\frac{\phi_{2}}{1+\beta_{21}n_{1j}^{*}}&s_{2a}\\\ \end{pmatrix}.$
(S28)

The $B_{3}$ matrix corresponds to the projection matrix of species 2, which
implies that the exclusion equilibrium is not stable when species 2 has a
growth rate larger than 1 (i.e., the modulus of $B_{3}$’s leading eigenvalue
is larger than 1). The stability of this two-dimensional system can be
investigated thanks to the Jury conditions, given by the equations (S29, S30,
S31):

$\displaystyle 1-\operatorname{tr}J(\mathbf{n})+\det J(\mathbf{n})$
$\displaystyle>0$ (S29) $\displaystyle 1+\operatorname{tr}J(\mathbf{n})+\det
J(\mathbf{n})$ $\displaystyle>0$ (S30) $\displaystyle 1-\det J(\mathbf{n})$
$\displaystyle>0$ (S31)

with $J(\mathbf{n})$ the Jacobian matrix of the two-compartment system whose
stability is being evaluated. If these three conditions hold, the system is
stable. If the first condition (equation (S29)) is violated, one of the
eigenvalues of $J(\mathbf{n})$ is larger than 1, which means in our context
that the excluded species can invade the community. The violations of the
conditions (S30) and (S31) correspond to the creation of limit cycles (Neubert
and Caswell, 2000). We then check these conditions on the $B_{3}$ matrix. We
have:

$\displaystyle\operatorname{tr}B_{3}$
$\displaystyle=(1-\gamma_{2})\frac{\phi_{2}}{1+\beta_{21}n_{1j}^{*}}+s_{2a}$
$\displaystyle\det B_{3}$
$\displaystyle=(1-\gamma_{2})\frac{\phi_{2}}{1+\beta_{21}n_{1j}^{*}}s_{2a}-\gamma_{2}\frac{\phi_{2}}{1+\beta_{21}n_{1j}^{*}}\frac{\pi_{2}}{1+\alpha_{21}n_{1a}^{*}}.$

We start with the third condition (S31):

$1-(1-\gamma_{2})\frac{\phi_{2}}{1+\beta_{21}n_{1j}^{*}}s_{2a}+\gamma_{2}\frac{\phi_{2}}{1+\beta_{21}n_{1j}^{*}}\frac{\pi_{2}}{1+\alpha_{21}n_{1a}^{*}}>0.$
(S32)

This equation will always be satisfied with biologically meaningful parameters
and strictly competitive interaction between species:

  * •

The third term of (S32) is positive because each parameter is positive,

  * •

$(1-\gamma_{2})\frac{\phi_{2}}{1+\beta_{21}n_{1j}^{*}}s_{2a}<1$ because
$0\leq\gamma\leq 1$, $0\leq\phi\leq 1$ and all $\beta$’s are positive.

Since $\operatorname{tr}B_{3}>0$, the first Jury’s condition implies the
second ((S29) $\Rightarrow$ (S30)). Then, we only have to check the first Jury
condition to describe the stability of the system:

$\displaystyle
1-(1-\gamma_{2})\frac{\phi_{2}}{1+\beta_{21}n_{1j}^{*}}+s_{2a}+(1-\gamma_{2})\frac{\phi_{2}}{1+\beta_{21}n_{1j}^{*}}s_{2a}-\gamma_{2}\frac{\phi_{2}}{1+\beta_{21}n_{1j}^{*}}\frac{\pi_{2}}{1+\alpha_{21}n_{1a}^{*}}>0$
$\displaystyle(1-\gamma_{2})\frac{\phi_{2}}{1+\beta_{21}n_{1j}^{*}}-(1-\gamma_{2})\frac{\phi_{2}}{1+\beta_{21}n_{1j}^{*}}s_{2a}+\gamma_{2}\frac{\phi_{2}}{1+\beta_{21}n_{1j}^{*}}\frac{\pi_{2}}{1+\alpha_{21}n_{1a}^{*}}<1-s_{2a}$
$\displaystyle(1-\gamma_{2})\phi_{2}(1-s_{2a})<(1-s_{2a})(1+\beta_{21}n_{1j}^{*})-\frac{\gamma_{2}\phi_{2}\pi_{2}}{1+\alpha_{21}n_{1a}^{*}}$
$\displaystyle(1-\gamma_{2})\phi_{2}\frac{1}{1+\beta_{21}n_{1j}^{*}}+\frac{\gamma_{2}\phi_{2}\pi_{2}}{1-s_{2a}}\frac{1}{(1+\beta_{21}n_{1j}^{*})(1+\alpha_{21}n_{1a}^{*})}<1.$

We finally have the following condition for the stability of the equilibrium:

$\frac{C_{2}}{1+\beta_{21}n_{1j}^{*}}+\frac{D_{2}}{(1+\beta_{21}n_{1j}^{*})(1+\alpha_{21}n_{1a}^{*})}<1$
(S33)

with $C_{2}=(1-\gamma_{2})\phi_{2}$ and
$D_{2}=\frac{\pi_{2}\gamma_{2}\phi_{2}}{1-s_{2a}}$.

We obtain the coexistence outcome domains given in Fig. S1 with the $\beta$
coefficients given in Table S1. The results are qualitatively similar to those
presented in the main text with adults $\rightarrow$ adults + juveniles.

Figure S1: Outcomes of the model with a metamorphosis structure predicted by invasion analysis, depending on the coefficients of inter-specific competition on fertility ($\alpha_{12}$ and $\alpha_{21})$ with $\beta$ coefficients indicating coexistence, exclusion of species 1 and priority effect. The dotted lines correspond to the outcome boundaries of panel A which is used as a reference. Figure | $\beta_{12}$ | $\beta_{21}$
---|---|---
C | 0.06 | 0.06
D | 0.11 | 0.05
E | 0.12 | 0.12
Table S1: $\beta$ competition coefficients used to obtain contrasting
situations of coexistence.

##  Appendix Supplement C Dynamics of single-species and two-species models

Here we explore the dynamics of single- and two-species models in four
contrasted life-history configurations (identical to those previously
considered by Neubert and Caswell 2000 in single-species models where one
vital rate has Ricker density-dependence):

  * •

$s_{a}>0,\gamma<1$ which corresponds to the parameters used for illustration
in the main text

  * •

$s_{a}>0,\gamma=1$: iteroparity with fast maturation

  * •

$s_{a}=0,\gamma<1$: semelparity and delayed maturation

  * •

$s_{a}=0,\gamma=1$: semelparity and fast maturation.

For single-species models (which we use when performing invasion analysis at
the exclusion equilibria), only the last case (i.e. semelparity and fast
maturation) leads to 2-cycles with the remainder of parameters as in the main
text (see figures below). This stays true when multiplying the fecundity by
10—a choice motivated by the fact that high fecundities generally enhance
cycling. For the metamorphosis model, the semelparity and fast maturation case
leads to a fixed point when initial conditions are strictly positive for
juveniles and adults. Indeed, Cushing et al. (2007) covers the semelparous,
fast maturation version of metamorphosis model and showed that while the
nontrivial equilibrium for the single-species model is stable on the positive
cone, it is not so when starting on the boundary with either $n_{j}=0$ or
$n_{a}=0$, which can generate so-called synchronous two-cycles (Theorem 2.1,
Cushing et al. 2007).

For two-species models, the case $s_{a}=0,\gamma=1$ generates 2-cycles but
$s_{a}=0,\gamma<1$ can also generate 2-cycles for some parameter sets, and
some different kinds of cycles (see below).

Figure S2: Population dynamics for the single species model and competition
generated by adults (model 1) with $n_{j}(0)=50$ and $n_{a}(0)=50$ (in cyan)
and for $n_{j}(0)=0$ and $n_{a}(0)=50$ (in purple). Similar dynamics are
obtained with both initial conditions for model 1.  Figure S3: Population
dynamics for the single species model and metamorphosis structure with
$n_{j}(0)=50$ and $n_{a}(0)=50$ (in cyan) and for $n_{j}(0)=0$ and
$n_{a}(0)=50$ (in purple). A stable fixed point is always obtained with
$n(0)=(50,50)$ for model 2.  Figure S4: Population dynamics for the two-
species models 1, where competition is generated by adults, with the three
parameter sets built to generate three different scenarios of coexistence. We
consider 4 different life history configurations, always with the initial
conditions $n_{i}(0)=(0,50)$ to allow for all accessible 2-cycles.

For the reader interested in dynamical aspects, we recommend to consult
Cushing et al. (2007). To facilitate the reading, we make a correspondence
below between their parameter notations and ours. Their model writes

$\begin{split}J_{t+1}=b_{1}\frac{1}{1+d_{1}A_{t}}A_{t}\\\
A_{t+1}=s_{1}\frac{1}{1+J_{t}+c_{1}j_{t}}J_{t}\\\
j_{t+1}=b_{2}\frac{1}{1+d_{2}a_{t}}a_{t}\\\
a_{t+1}=s_{2}\frac{1}{1+c_{2}J_{t}+j_{t}}j_{t}\end{split}$ (S34)

with a state-space vector $(J_{t},A_{t},j_{t},a_{t})$ where one of the species
is written in capitals.

Cushing et al. (2007) is therefore a special case of

$\mathbf{A}_{1}^{M}(\mathbf{n})=\begin{pmatrix}\frac{(1-\gamma_{1})\phi_{1}}{1+\beta_{11}n_{1j}+\beta_{12}n_{2j}}&\frac{\pi_{1}}{1+\alpha_{11}n_{1a}+\alpha_{12}n_{2a}}\\\
\frac{\gamma_{1}\phi_{1}}{1+\beta_{11}n_{1j}+\beta_{12}n_{2j}}&s_{1a}\end{pmatrix}$
(S35)

that is,

$\mathbf{A}_{1}^{M}(\mathbf{n})=\begin{pmatrix}0&\frac{\pi_{1}}{1+\alpha_{11}n_{1a}+0\times
n_{2a}}\\\
\frac{\phi_{1}}{1+\beta_{11}n_{1j}+\beta_{12}n_{2j}}&0\end{pmatrix}.$ (S36)

Note that
$\alpha_{11}=d_{1},\alpha_{12}=\alpha_{21}=0,\alpha_{22}=d_{2},\pi_{1}=b_{1},\pi_{2}=b_{2},\beta_{11}=1,\beta_{12}=c_{1},\beta_{21}=c_{2},\beta_{22}=1,\phi_{1}=s_{1},\phi_{2}=s_{2}$.

### Why is this sub-case of the “metamorphosis life cycle” model more amenable
to analytical insights?

#### Going to $t+2$

$n_{1j}(t+2)=\frac{\pi_{1}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}n_{1a}(t+1)}}{1+\alpha_{11}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}n_{1a}(t+1)}+0\times
n_{2a}(t+1)}$ (S37)

with
$\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}n_{1a}(t+1)=\frac{\phi_{1}n_{1j}(t)}{1+\beta_{11}n_{1j}(t)+\beta_{12}n_{2j}(t)}$.
This leads to

$n_{1j}(t+2)=\frac{\pi_{1}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\frac{\phi_{1}n_{1j}(t)}{1+\beta_{11}n_{1j}(t)+\beta_{12}n_{2j}(t)}}}{1+\alpha_{11}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\frac{\phi_{1}n_{1j}(t)}{1+\beta_{11}n_{1j}(t)+\beta_{12}n_{2j}(t)}}}$
(S38)

which, multiplying by $1+\beta_{11}n_{1j}(t)+\beta_{12}n_{2j}(t)$ the
numerator and denominator, leads to

$\begin{split}n_{1j}(t+2)=\frac{\pi_{1}\phi_{1}n_{1j}(t)}{1+\beta_{11}n_{1j}(t)+\beta_{12}n_{2j}(t)+\alpha_{11}\phi_{1}n_{1j}(t)}\\\
=\frac{\pi_{1}\phi_{1}n_{1j}(t)}{1+(\beta_{11}+\alpha_{11}\phi_{1})n_{1j}(t)+\beta_{12}n_{2j}(t)}.\end{split}$
(S39)

This equation gives $n_{1j}(t+2)$ as a function of $n_{1j}(t)$ and $n_{2j}(t)$
and is equivalent to the classical 2-species Leslie-Gower model (except for
the doubled timestep), and therefore we know its behaviour (fixed point), as
shown by Cushing et al. (2007). Because $n_{1a}(t+1)$ only depends on
$n_{1j}(t)$ and $n_{2j}(t)$ we can deduce that it won’t change the dynamical
behaviour.

##  Appendix Supplement D Simulations of metamorphosis model structure (model
2) extended to a $S$-species context

The species richness at the end of simulations for the extended model 2 with
$S$ species are given in Fig. S5.

Figure S5: Histograms of the number of species persisting at equilibrium for
200 permutations of inter-specific competition coefficients. The bar in gray
indicates the number of species persisting at time $t=3000$ time steps for the
original $\alpha_{ij}$ and $\beta_{ij}$ with
$\text{Corr}(\alpha_{ij},\beta_{ij})<0$. From left to right, columns
correspond the communities with initially 5, 10 and 40 species. From top to
bottom, the lines correspond to situations with high, average and low variance
on the vital rates.

We find here that with a slightly different competition structure (adults only
compete with adults and juveniles with juveniles), results are only slightly
different from those of model 1. They go in the same general direction: it is
difficult to obtain a positive effect of opposite competitive hierarchies on
species richness at equilibrium. One could have believed that there was
something special with Moll and Brown (2008)’s life cycle structure with
separated adults from juveniles, which would have made many-species
coexistence despite competition a little bit easier. We found the opposite:
for this particular set of simulations, we have not been able to observe more
coexistence with opposite competitive hierarchies on $\alpha$ and $\beta$,
even for very similar species. This only strengthens the generality of our
finding that opposite competitive hierarchies are unlikely to be a general
explanation to many-species coexistence.