# Heterogeneity of Interaction Strengths and Its
Consequences on Ecological Systems

Zachary Jackson, BingKan Xue

(Department of Physics, University of Florida, Gainesville, FL, United States)

###### Abstract

Ecosystems are formed by networks of species and their interactions.
Traditional models of such interactions assume a constant interaction strength
between a given pair of species. However, there is often significant trait
variation among individual organisms even within the same species, causing
heterogeneity in their interaction strengths with other species. The
consequences of such heterogeneous interactions for the ecosystem have not
been studied systematically. As a theoretical exploration, we analyze a simple
ecosystem with trophic interactions between two predators and a shared prey,
which would exhibit competitive exclusion in models with homogeneous
interactions. We consider several scenarios where individuals of the prey
species differentiate into subpopulations with different interaction
strengths. We show that in all these cases, whether the heterogeneity is
inherent, reversible, or adaptive, the ecosystem can stabilize at a new
equilibrium where all three species coexist. Moreover, the prey population
that has heterogeneous interactions with its predators reaches a higher
density than it would without heterogeneity, and can even reach a higher
density in the presence of two predators than with just one. Our results
suggest that heterogeneity may be a naturally selected feature of ecological
interactions that have important consequences for the stability and diversity
of ecosystems.

##  1 Introduction

Traditional physical systems involve interactions between objects
characterized by universal coupling constants, such as the Newtonian constant
of gravitation and the Coulomb constant for electrostatic forces. This is true
whether the objects are subatomic particles, electric charges, or celestial
bodies. The objects can differ in the amount of charge or mass they have, but
otherwise interact in the same way with each other. This is in contrast to
biological systems, where individual organisms are described by a large number
of traits, be it morphological, metabolic, or behavioral, which are not all
identical even between individuals of the same species [1, 2, 3]. Any of these
traits can affect the way an individual interacts with other organisms and the
environment. Thus, for biological interactions, the coupling constants (or
“interaction strengths”) themselves can be heterogeneous among individuals.
Such interaction strengths can also vary in time due to behavioral changes,
seasonal variations, or different life stages of an organism.

However, when trying to apply dynamical models to biological systems, the
heterogeneity among individuals is often left out, so that a whole population
is treated as having the same interaction strengths. A classic example is the
Lotka-Volterra model of trophic interaction between two species, a predator
and a prey. The predation rate is assumed to be proportional to the population
densities of both species (in a quadratic form that loosely resembles some of
the physical interactions mentioned above). This may be true if each
population is homogeneous, so that only the population size or density is
relevant. However, if we take into account the individual variation in various
external or internal traits, the interaction strength between the predator and
prey need not be the same for every individual. Intraspecific trait variation
can have significant ecological effects [4, 1]. Thus, using an “average”
interaction strength and ignoring the heterogeneity among individuals may
cause models to miss important features. Here we address this problem by
analyzing simple ecological models that demonstrate nontrivial consequences of
heterogeneity in interaction strengths.

As a case study, we analyze an ecosystem with “exploitative competition”
between two predators feeding on the same prey species. Traditional models of
such ecosystems treat each species as being homogeneous, so that the
interaction strength between each pair is constant. We will incorporate
heterogeneity in how the prey species interact with the predators, which can
happen if there is trait variation among the prey population. Our model allows
us to study several kinds of heterogeneity, including what we call “inherent”,
“reversible”, and “adaptive” heterogeneities, depending on whether the prey
phenotype is determined at birth, can change reversibly and stochastically, or
adapts plastically to the density of predators. In all cases, we show that
heterogeneous interactions lead to new phenomena in the population dynamics of
species. Without heterogeneity, the expected outcome of the ecosystem is that
one of the predators drives the other to extinction, a phenomenon known as
“competitive exclusion” [5, 6, 7]. However, if the interaction strength varies
among the prey species, it turns out that the system can stabilize in a state
where both predators coexist. For some range of parameters, we observe
emergent facilitation between the predators [8], such that the presence of one
predator allows the other predator to persist. Moreover, the prey species can
reach a higher abundance than when it interacts with only one predator,
suggesting that heterogeneity in interactions may be an evolutionarily
favorable feature for the population.

##  2 System with homogeneous interactions

As a null model, we first consider an ecosystem with homogeneous interactions.
The two predators and one prey species are modeled by a Lotka-Volterra
dynamical system,

$\displaystyle\dot{A}$ $\displaystyle=A\,(\alpha_{AC}\,C-\beta_{A})$ (1a)
$\displaystyle\dot{B}$ $\displaystyle=B\,(\alpha_{BC}\,C-\beta_{B})$ (1b)
$\displaystyle\dot{C}$
$\displaystyle=C\,(\beta_{C}-\alpha_{CC}\,C-\varepsilon_{A}\,\alpha_{AC}\,A-\varepsilon_{B}\,\alpha_{BC}\,B)$
(1c)

Here $A$, $B$, and $C$ represent the density of each species in a spatially
well-mixed system, where $A$ and $B$ are the predators and $C$ is the prey
(Fig. 1a). Parameters $\beta_{A}$ and $\beta_{B}$ are the death rates of $A$
and $B$, respectively; $\alpha_{AC}$ and $\alpha_{BC}$ are the predation rates
of $A$ and $B$ on $C$. Species $C$ has a birth rate $\beta_{C}$ and an
intraspecific competition strength $\alpha_{CC}$. The $\varepsilon$’s are the
efficiency by which an amount of $C$ consumed is converted to the increase of
the predator population. We can rescale the variables by
$t\leftarrow\beta_{C}\,t$,
$A\leftarrow\frac{\varepsilon_{A}\alpha_{CC}}{\beta_{C}}A$,
$B\leftarrow\frac{\varepsilon_{B}\alpha_{CC}}{\beta_{C}}B$, and
$C\leftarrow\frac{\alpha_{CC}}{\beta_{C}}C$, so that the equations above
simplify to

$\displaystyle\dot{A}$ $\displaystyle=A\,(a\,C-a_{0})$ (2a)
$\displaystyle\dot{B}$ $\displaystyle=B\,(b\,C-b_{0})$ (2b)
$\displaystyle\dot{C}$ $\displaystyle=C\,(1-\,C-a\,A-b\,B)$ (2c)

where the rescaled parameters are $a=\frac{\alpha_{AC}}{\alpha_{CC}}$,
$a_{0}=\frac{\beta_{A}}{\beta_{C}}$, $b=\frac{\alpha_{AC}}{\alpha_{CC}}$,
$b_{0}=\frac{\beta_{B}}{\beta_{C}}$.

Figure 1: An ecosystem with two predators ($A$, $B$) and one prey ($C$). (a) A
homogeneous system where all individuals of a species have the same
interaction strength with another species. Arrows represent trophic
interaction pointing from prey to predator. (b) Heterogeneity is added by
splitting the prey into two types ($C_{1}$, $C_{2}$), each with their own
interaction strengths with the predators. Double-sided arrow represents the
exchange of individuals between the prey subpopulations. (c) An equivalent
description using the total prey population $C$ and the prey type composition
$\lambda$. Arrows represent “effective” interaction strengths that depend on
$\lambda$.

There are 4 equilibrium states of this system, which are labeled $P_{O}$,
$P_{C}$, $P_{A}$, and $P_{B}$, as shown in Figure 2. They all belong to the
surface $a_{0}A+b_{0}B=C(1-C)$ restricted to the non-negative octant of the
$A$-$B$-$C$ space. $P_{O}$ is the point where all three species have zero
population sizes, which is an unstable equilibrium. $P_{C}$ is where only $C$
persists; $P_{A}$ and $P_{B}$ are where $A$ or $B$ coexists with $C$,
respectively. In general, only one of $P_{A}$, $P_{B}$, and $P_{C}$ can be
stable for a given set of parameters. This demonstrates the competitive
exclusion principle, by which two consumers (predators) of the same resource
(prey) cannot coexist. The persistent predator is the one that has a lower
$R^{*}$ value (i.e., the minimal prey density required to sustain a predator,
$R^{*}_{A}=a_{0}/a$ and $R^{*}_{B}=b_{0}/b$). Only if the parameters are fine-
tuned, such that the two predators have equal $R^{*}$, will all three species
coexist. In such a fine-tuned system, there is a continuum of possible
equilibria that form a line attractor $\mathcal{L}$ shown in Figure 2.

Figure 2: $A$-$B$-$C$ space showing the locations of the equilibrium points.
$P_{O}$ (black) is an unstable equilibrium at the origin where all species go
extinct. $P_{C}$ (yellow) is where only $C$ survives. $P_{A}$ (red) and
$P_{B}$ (blue) represent the persistence of one predator ($A$ and $B$,
respectively) and the prey. These points each lie on a parabola in the $B=0$
and $A=0$ planes, respectively. $P_{N}$ (green) is a point where all three
species coexist, which lies on a line $\mathcal{L}$ (purple) that is a line
attractor for $\lambda=\lambda^{*}$. All five possible equilibria are on the
surface $a_{0}A+b_{0}B=C(1-C)$, shown in gray.

##  3 Heterogeneous Interactions

Now consider a generalized model in which the species interactions are
heterogeneous. A natural way of introducing heterogeneity in the system is by
having a species diversify into subpopulations with different interaction
strengths. This way of modeling heterogeneity is useful as it can describe
different kinds of heterogeneity. For example, the subpopulations could
represent polymorphic traits that are genetically determined or result from
plastic response to heterogeneous environments. A population could also be
divided into local subpopulations in different spatial patches, which can
migrate between patches and may face different local predators. We can also
model different behavioral modes as subpopulations that, for instance, spend
more time foraging for food or hiding from predators. We will study several
examples of these different kinds of heterogeneity after we introduce a common
mathematical framework.

We will focus on the simple case where only the prey species splits into two
types, $C_{1}$ and $C_{2}$, as illustrated in Figure 1b. The situation is
interesting when the predators $A$ and $B$ prefer different prey types such
that, without loss of generality, $A$ consumes $C_{1}$ more readily while $B$
prefers $C_{2}$ (i.e., $a_{1}/a_{0}>b_{1}/b_{0}$ and
$b_{2}/b_{0}>a_{2}/a_{0}$, which is equivalent to the condition that the
nullclines of $A$ and $B$ cross, see Section 4.1). The arrows between $C_{1}$
and $C_{2}$ in Figure 1b represent the exchange of individuals between the two
subpopulations, which can happen by various mechanisms considered below. Such
exchange as well as intraspecific competition between $C_{1}$ and $C_{2}$
result from the fact that the two prey types remain a single species.

The system is now described by an enlarged Lotka-Volterra system with four
variables, $A$, $B$, $C_{1}$, and $C_{2}$:

$\displaystyle\dot{A}$
$\displaystyle=\alpha_{A1}AC_{1}+\alpha_{A2}AC_{2}-\beta_{A}A$ (3a)
$\displaystyle\dot{B}$
$\displaystyle=\alpha_{B1}BC_{1}+\alpha_{B2}BC_{2}-\beta_{B}B$ (3b)
$\displaystyle\dot{C_{1}}$
$\displaystyle=C_{1}(\beta_{C}-\alpha_{CC}C)-\varepsilon_{A}\alpha_{A1}C_{1}A-\varepsilon_{B}\alpha_{B1}C_{1}B-\sigma_{1}C_{1}+\sigma_{2}C_{2}$
(3c) $\displaystyle\dot{C_{2}}$
$\displaystyle=C_{2}(\beta_{C}-\alpha_{CC}C)-\varepsilon_{A}\alpha_{A2}C_{2}A-\varepsilon_{B}\alpha_{B2}C_{2}B+\sigma_{1}C_{1}-\sigma_{2}C_{2}$
(3d)

Here we assume that the prey types $C_{1}$ and $C_{2}$ have the same
intraspecific competition strength, but different interaction strengths with
$A$ and $B$. The parameters in these equations and their meanings are listed
in Table 1. For the convenience of analysis, we will transform the variables
$C_{1}$ and $C_{2}$ to another pair of variables $C$ and $\lambda$, where
$C\equiv C_{1}+C_{2}$ is the total population of $C$ as before, and
$\lambda\equiv C_{2}/(C_{1}+C_{2})$ represents the composition of the
population. After this transformation and rescaling of some variables
(described in Appendix), the new dynamical system can be written as:

$\displaystyle\dot{A}$
$\displaystyle=A\,\big{(}C(a_{1}(1-\lambda)+a_{2}\lambda)-a_{0}\big{)}$ (4a)
$\displaystyle\dot{B}$
$\displaystyle=B\,\big{(}C(b_{1}(1-\lambda)+b_{2}\lambda)-b_{0}\big{)}$ (4b)
$\displaystyle\dot{C}$
$\displaystyle=C\,\big{(}1-C-A(a_{1}(1-\lambda)+a_{2}\lambda)-B(b_{1}(1-\lambda)+b_{2}\lambda)\big{)}$
(4c) $\displaystyle\dot{\lambda}$
$\displaystyle=\lambda(1-\lambda)\,\big{(}A(a_{1}-a_{2})+B(b_{1}-b_{2})\big{)}+\eta_{1}(1-\lambda)-\eta_{2}\lambda$
(4d)

Here, $a_{i}$ and $b_{i}$ are the (rescaled) feeding rates of the predators on
the prey type $C_{i}$; $a_{0}$ and $b_{0}$ are the death rates of the
predators as before; $\eta_{1}$ and $\eta_{2}$ are the exchange rates of the
prey types (Table 1). The latter can be functions of other variables,
representing different kinds of heterogeneous interactions that we will study
below.

Original | Rescaled | Meaning
---|---|---
$\alpha_{A1}$ | $a_{1}$ | consumption rate of $C_{1}$ by $A$
$\alpha_{A2}$ | $a_{2}$ | consumption rate of $C_{2}$ by $A$
$\alpha_{B1}$ | $b_{1}$ | consumption rate of $C_{1}$ by $B$
$\alpha_{B2}$ | $b_{2}$ | consumption rate of $C_{2}$ by $B$
$\alpha_{CC}$ | $1$ | intraspecific competition rate of $C$
$\beta_{A}$ | $a_{0}$ | death rate of $A$
$\beta_{B}$ | $b_{0}$ | death rate of $B$
$\beta_{C}$ | $1$ | birth rate of $C$
$\varepsilon_{A}$ | $1$ | biomass conversion of $C$ to $A$
$\varepsilon_{B}$ | $1$ | biomass conversion of $C$ to $B$
$\sigma_{1}$ | $\eta_{1}$ | exchange rate from $C_{1}$ to $C_{2}$
$\sigma_{2}$ | $\eta_{2}$ | exchange rate from $C_{2}$ to $C_{1}$
Table 1: Parameters of the model (before and after rescaling) and their
meanings.

The variable $\lambda$ can be considered an internal degree of freedom within
the $C$ population. In all models we study below, the dynamics of $\lambda$ is
such that it stabilizes to a special value $\lambda^{*}$ (a bifurcation
point), as shown in Figure 3. As a result, a new equilibrium point $P_{N}$
appears in the $A$-$B$-$C$ space (on the line $\mathcal{L}$, see Figure 2), at
which all three species coexist. If $\lambda$ were fixed at any other value,
there would be exclusion of one of the predators (Figure 3a), just as in the
case of homogeneous interactions. Thus, with heterogeneous interactions, the
prey composition $\lambda$ dynamically adjusts itself to the value
$\lambda^{*}$ at which a new coexistence phase emerges, without having to
fine-tune the interaction strengths. The exact conditions for this phase to be
stable are detailed in the Appendix.

Figure 3: (a) Equilibrium population of each species $X=A$, $B$, or $C$, with
$\lambda$ held fixed at different values. Solid curves represent stable
equilibria and dashed curves represent unstable equilibria. The vertical
dashed line is where $\lambda=\lambda^{*}$, which is also the bifurcation
point. Notice that the equilibrium population of $C$ is maximized at this
point. (b) Time series of $\lambda$ for systems with each kind of
heterogeneity. All three systems stabilize at the same $\lambda^{*}$ value,
which is the bifurcation point in panel (a). Parameters used here are
$(a_{0},a_{1},a_{2},b_{0},b_{1},b_{2},\rho,\eta_{1},\eta_{2},\kappa)=(0.25,0.5,0.2,0.4,0.2,0.6,0.5,0.05,0.05,50)$.

###  3.1 Inherent heterogeneity

We first consider a scenario where individuals of the prey species are born in
two types with a fixed ratio, such that a fraction $\rho$ of the newborns are
$C_{2}$ and $(1-\rho)$ are $C_{1}$. This could describe dimorphic traits, such
as in sexual dimorphism [9, 10] where males and females are produced in a
genetically determined ratio regardless of the composition of the population.
We call this “inherent” heterogeneity, because individuals are born with a
certain type and cannot change in later stages of life. The prey type given at
birth will determine the individual’s interaction strength with the predators.
This kind of heterogeneity can be described by Eq. (4d) with
$\eta_{1}=\rho(1-C)$ and $\eta_{2}=(1-\rho)(1-C)$ (see Appendix).

Figure 4: Phase diagrams showing regions of parameter space identified by the
stable equilibrium points. Yellow region represents $P_{C}$, red represents
$P_{A}$, blue represents $P_{B}$, and green represents $P_{N}$. (a–d)
Different values of the newborn composition $\rho$ for the model of inherent
heterogeneity. In the extreme cases of $\rho=0$ and $1$ the prey is
homogeneous and there is no coexistence of the predators. Parameters used here
are $(a_{0},a_{1},b_{0},b_{2})=(0.3,0.5,0.4,0.6)$.

The stable equilibrium of the system can be represented by phase diagrams that
show the identities of the species at equilibrium. We plot these phase
diagrams by varying the parameters $a_{2}$ and $b_{1}$ while keeping $a_{1}$
and $b_{2}$ constant. As shown in Figure 4, the equilibrium state depends on
the parameter $\rho$. In the limit $\rho=0$ or $1$, we recover the homogeneous
case because only one type of $C$ is produced. The corresponding phase
diagrams (Figure 4a,d) contain only two phases where one of the predators
persist, illustrating the competitive exclusion principle. For intermediate
values of $\rho$, however, there is a new phase of coexistence that separates
the two exclusion phases.

There are two such regions of coexistence, which touch at a middle point and
open toward the bottom left and upper right, respectively. The middle point is
at $(a_{2}/a_{0}=b_{2}/b_{0},b_{1}/b_{0}=a_{1}/a_{0})$, where the feeding
preferences of the two predators are identical (hence their niches fully
overlap). Towards the origin and the far upper right, the predators consume
one type of $C$ each (hence their niches separate). The coexistence region in
the bottom left is where the feeding rates of the predators are the lowest
overall. There can be a region (yellow) where both predators go extinct, if
their primary prey type alone is not enough to sustain each predator. As shown
in Figure 5, increasing the productivity of the system by increasing the birth
rate ($\beta_{C}$) of the prey eliminates this extinction region, whereas
lowering productivity allows the extinction region to subsume the lower
coexistence region.

Figure 5: Phase diagrams showing how the coexistence region changes with
varying productivity (prey birth rate $\beta_{C}$). (a–d) The re-scaled
predator death rates $a_{0}$ and $b_{0}$ are increased in proportion to each
other, which is equivalent to decreasing $\beta_{C}$ (see Appendix for
parameter transformation). The lower left region of coexistence (green) is
absorbed into the prey-only region (yellow) as productivity is decreased.
Parameters used here are $(a_{1},b_{2},\rho)=(0.5,0.6,0.5)$.

The new equilibrium is not only where the predators $A$ and $B$ can coexist,
but also where the prey species $C$ grows to a larger density than what is
possible for a homogeneous population. This is illustrated in Figure 3b, which
shows the equilibrium population of $C$ if we hold $\lambda$ fixed at
different values. The point $\lambda=\lambda^{*}$ is where the system with a
dynamic $\lambda$ is stable, and also where the population of $C$ is
maximized. That means the population automatically stabilizes at the optimal
composition of prey types. Moreover, the value of $C^{*}$ at this coexistence
point can even be larger than the equilibrium population of $C$ when there is
only one predator $A$ or $B$. This will be discussed further in Section 4.4.
These results suggest that heterogeneity in interaction strengths can
potentially be a strategy for the prey population to leverage the effects of
multiple predators against each other to improve survival.

###  3.2 Reversible heterogeneity

We next consider a scenario where individuals can switch their types. This
kind of heterogeneity can model reversible changes of phenotypes, i.e., trait
changes that affect the prey’s interaction with predators but are not
permanent. For example, changes in coat color or camouflage [11, 12, 13],
physiological changes such as defense [14], and biomass allocation among
tissues [15, 16]. One could also think of the prey types as subpopulations
within different spatial patches, if each predator hunts at a preferred patch
and the prey migrate between the patches [17, 18]. We can model this
“reversible” kind of heterogeneity by introducing switching rates from one
prey type to the other. In Eq. (4d), $\eta_{1}$ and $\eta_{2}$ now represent
the switching rates per capita from $C_{1}$ to $C_{2}$ and from $C_{2}$ to
$C_{1}$, respectively. Here we study the simplest case where both rates are
fixed.

In the absence of the predators, the natural composition of the prey species
given by the switching rates would be
$\rho\equiv\eta_{1}/(\eta_{1}+\eta_{2})$, and the rate at which $\lambda$
relaxes to this natural composition is $\gamma\equiv\eta_{1}+\eta_{2}$.
Compared to the previous scenario where we had only one parameter $\rho$, here
we have an additional parameter $\gamma$ that modifies the behavior of the
system. Figure 6 shows phase diagrams for the system as $\rho$ is fixed and
$\gamma$ varies. We again find the new equilibrium $P_{N}$ where all three
species coexist. When $\gamma$ is small, the system has a large region of
coexistence. As $\gamma$ is increased, this region is squeezed into a border
between the two regions of exclusion. However, this is different from the
exclusion we see in the limit $\rho\to 0$ or $1$ for the inherent
heterogeneity, where the borders are horizontal or vertical (Figure 4). Here
the predators exclude each other despite having a mixture of prey types in the
population.

Figure 6: Phase diagram showing how the coexistence region changes with the
switching rates $\eta_{1}$ and $\eta_{2}$. (a–d) The coexistence region
shrinks as the switching rates increase. When these rates are large, $\lambda$
is effectively held constant, so the system is equivalent to a homogeneous
system and coexistence disappears. Parameters used here are
$(a_{0},a_{1},b_{0},b_{2})=(0.3,0.5,0.4,0.6)$

This special limit can be understood as follows. For a large $\gamma$,
$\lambda$ is effectively set to a constant value equal to $\rho$, because it
has a very fast relaxation rate. In other words, the prey types exchange so
often that the population always maintains a constant composition. In this
limit, the system behaves as if it were a _homogeneous_ system with effective
interaction strengths $a_{\text{eff}}=(1-\rho)\,a_{1}+\rho\,a_{2}$ and
$b_{\text{eff}}=(1-\rho)\,b_{1}+\rho\,b_{2}$. As in a homogeneous system,
there will be competitive exclusion between the predators instead of
coexistence. This demonstrates that having a constant level of heterogeneity
is not sufficient to cause coexistence. The overall composition of the
population must be able to change dynamically as a result of population growth
and consumption by predators.

An interesting behavior is seen when we examine a point inside the shrinking
coexistence region as $\gamma$ is increased. Typical trajectories of the
system for such parameter values are shown in Fig. 7. As $\gamma$ increases,
the system relaxes to the line $\mathcal{L}$ quickly, then slowly crawls along
it towards the final equilibrium point $P_{N}$. This is because increasing
$\gamma$ increases the speed that $\lambda$ relaxes to $\lambda^{*}$, and when
$\lambda\to\lambda^{*}$, $\mathcal{L}$ becomes marginally stable. Therefore,
the attraction to $\mathcal{L}$ in the perpendicular direction is strong, but
the attraction towards the equilibrium point along $\mathcal{L}$ is weak. This
leads to a long transient behavior that makes the system appear to reach no
equilibrium in a limited time [19, 20]. It is especially true when there is
noise in the dynamics, which will cause the system to diffuse along
$\mathcal{L}$ with only a weak drift towards the final equilibrium (Fig. 7).
Thus, the introduction of a fast timescale (quick relaxation of $\lambda$ due
to a large $\gamma$) actually results in a long transient.

Figure 7: Trajectories of the system projected in the $A$-$B$ plane, with
parameters inside the coexistence region (by holding the position of $P_{N}$
fixed). As $\gamma$ increases, the system tends to approach the line
$\mathcal{L}$ quickly and then crawl along it. The grey trajectory is with
independent Gaussian white noise ($\sim\mathcal{N}(0,0.5)$) added to each
variable’s dynamics. Noise causes the system to diffuse along $\mathcal{L}$
for a long transient period before coming to the equilibrium point $P_{N}$.
Parameters used here are
$(a_{0},a_{1},a_{2},b_{0},b_{1},b_{2})=(0.2,0.8,0.5,0.2,0.6,0.9)$.

###  3.3 Adaptive heterogeneity

A third kind of heterogeneity we consider is the change of interactions in
time. By this we mean an individual can actively change its interaction
strength with others in response to certain conditions. This kind of response
is often invoked in models of adaptive foraging behavior, where individuals
choose appropriate actions to maximize some form of fitness [21, 22]. For
example, we may consider two behaviors, resting and foraging, as our prey
types. Different predators may prefer to strike when the prey is doing
different things. In response, the prey may choose to do one thing or the
other depending on the current abundances of different predators. Such
behavioral modulation is seen, for example, in systems of predatory spiders
and grasshoppers [23]. Phenotypic plasticity is also seen in plant tissues in
response to consumers [24, 25, 26].

This kind of “adaptive” heterogeneity can be modeled by having switching rates
$\eta_{1}$ and $\eta_{2}$ that are time-dependent. Let us assume that the prey
species tries to maximize its population growth rate by switching to the more
favorable type. From Eq. (4c), we see that the growth rate of $C$ depends
linearly on the composition $\lambda$ with a coefficient
$u(A,B)\equiv(a_{1}-a_{2})A+(b_{1}-b_{2})B$. Therefore, when this coefficient
is positive, it is favorable for $C$ to increase $\lambda$ by switching to
type $C_{2}$. This can be achieved by having a positive switching rate
$\eta_{2}$ whenever $u(A,B)>0$. Similarly, whenever $u(A,B)<0$, it is
favorable for $C$ to switch to type $C_{1}$ by having a positive $\eta_{1}$.
In this way, the heterogeneity of the prey population will be constantly
adapting to the predator densities. We model such adaptive switching by making
$\eta_{1}$ and $\eta_{2}$ functions of the coefficient $u(A,B)$, e.g.,
$\eta_{2}(u)=1/(1+\mathrm{e}^{-\kappa u})$ and $\eta_{1}(u)=1-\eta_{2}(u)$.
The sigmoidal form of the functions means that the switching rate in the
favorable direction for $C$ is turned on quickly, while the other direction is
turned off. The parameter $\kappa$ controls the sharpness of this transition.

Phase diagrams for the system with different values of $\kappa$ are shown in
Figure 8. A larger $\kappa$ means the prey adapts its composition faster and
more optimally, which causes the coexistence region to expand. In the extreme
limit, the system changes its dynamics instantaneously whenever it crosses the
boundary where $u(A,B)=0$, like in a hybrid system [27]. Such a system can
still reach a stable equilibrium that lies on the boundary, if the flow on
each side of the boundary points towards the other side [28]. This is what
happens in our system and, interestingly, the equilibrium is the same three-
species coexistence point $P_{N}$ as in the previous scenarios. The region of
coexistence turns out to be largest in this limit (Figure 8d).

Figure 8: Phase diagrams showing how the coexistence region changes when the
switching rates $\eta_{i}$ dynamically adapt to the predator densities, so as
to maximize the growth rate of $C$. (a–d) As the sharpness $\kappa$ of the
sigmoidal decision function is increased, the region of coexistence expands.
In the limit $\kappa\to\infty$, the condition for coexistence becomes quite
simple (see Appendix). Parameters used here are
$(a_{0},a_{1},b_{0},b_{2})=(0.3,0.5,0.4,0.6)$.

Our results again demonstrate that the coexistence of the predators can be
viewed as a by-product of the prey’s strategy to maximize its own benefit. The
time-dependent case studied here represents a strategy that involves the prey
evaluating the risk posed by different predators. This is in contrast to the
scenarios studied above, where the prey population passively creates
phenotypic heterogeneity regardless of the presence of the predators. These
two types of behavior are analogous to the two strategies studied for
adaptation in varying environments, i.e., sensing and bet-hedging [29, 30].
The former requires accessing information about the current environment to
make optimal decisions, whereas the latter relies on maintaining a diverse
population to reduce detrimental effects caused by environmental changes. Here
the varying abundances of the predators play a similar role as the varying
environment. From this point of view, the heterogeneous interactions studied
here can be a strategy of the prey species that is evolutionarily favorable.

##  4 Discussion

Exploitative competition is a basic motif in the modeling of trophic
interactions. However, the exclusion of all but one consumer that share a
common resource as implied by such models is at odds with the high diversity
of species seen in natural ecosystems. This is known as the “paradox of the
plankton” [31], for which many explanations have been considered [6, 32, 33,
34, 35]. Our results suggest a new mechanism for species coexistence in these
ecosystems through heterogeneous interactions between the predators and prey.
Our model is general enough to describe many types of trait differentiation
within a species, including phenotypic polymorphism and switching, spatial
localization and migration, as well as behavioral changes and foraging
strategies. Our results are related to the following ecological concepts that
have been studied previously.

###  4.1 Resources competition and nullcline analysis

If we think of the two subpopulations of the prey as two resources, then the
competition between the two predators (consumers) can be analyzed using
nullclines. A nullcline is a contour line in the space of both resource
densities ($C_{1}$ and $C_{2}$), along which the net growth rate of a consumer
($\dot{A}$ or $\dot{B}$) is zero. As illustrated in Figure 9, there is one
nullcline for each consumer ($A$ and $B$). The nullcline is a generalization
of the $R^{*}$ value for a single (homogeneous) resource, and represents the
minimal combination of resource levels required for a consumer to sustain its
population. For resources that vary independently, the persistence of the
consumers can be determined from this picture. For example, if the resource
levels are above the $A$ nullcline and below the $B$ nullcline, then the $A$
population can grow but the $B$ population will decline, leading to the
exclusion of $B$ by $A$. By such analysis, the two consumers can coexist only
if the resource levels are precisely at the intersection of the nullclines
[36].

Figure 9: A schematic for nullcline analysis showing the prey subpopulations
as two resources. Nullclines for $A$ and $B$ are shown in black. Total
population of $C$ is represented by a line with slope $-1$. Prey composition
$\lambda$ is represented by the angle of a line going through the origin. The
equilibrium values of $C^{*}$ and $\lambda^{*}$ at three-species coexistence
is determined by the intersection of the nullclines. However, for our models
of heterogeneous interactions, the dynamics of $\lambda$ and its stability at
$\lambda^{*}$ are undetermined by this picture.

The nullcline analysis helps determine the location of the coexistent
equilibrium in the parameter space. However, the stability of this equilibrium
point cannot be derived from the picture, because in our model the two
resource types interact with each other dynamically. Indeed, they are
subpopulations of the same prey species, exchanging fluxes of individuals and
sharing a common carrying capacity. The ratio of the prey types ($\lambda$) is
determined by the angle of lines going through the origin (Figure 9). The
dynamics of $\lambda$ and of the total prey population $C$ both depend on the
densities of $A$ and $B$ (Eqs. (4c,4d)) and are not captured in the picture
here. Nevertheless, the nullcline analysis does show that the equilibrium
values of $C^{*}$ and $\lambda^{*}$ in the coexistence phase do not depend on
the parameters such as $\rho$ and $\gamma$, and in fact do not depend on the
form of heterogeneity at all (see Appendix). Beyond that, the main point of
our model is to demonstrate that the coexistent equilibrium is in fact stable,
and is robust to the different kinds of heterogeneity that we studied.

###  4.2 Emergent fitness equalization vs. niche separation

It is important to note that the coexistence of the predators shown here does
not result from the separation of niches, but rather from the convergence of
the predators to a common fitness, similar to the situation studied in [37].
The separation of niches would mean that each predator consumes only one prey
type. Here the niches of the predators overlap because they consume both prey
types but with different preferences. In [37] and in the models that we study
here, an internal variable of the prey population (the phenotype composition
in our case, and the defense level in [37]) is adjusted to ensure equal
fitness of the predators. Here fitness is used to mean the $R^{*}$ values of
the predators for a given internal trait value. The convergence of fitness is
emergent in that it is not the result of fine-tuning the predators’
consumption preferences (i.e., interaction strengths), but rather a result of
the dynamics of the prey composition that automatically stabilizes at the
critical value where $R^{*}$’s are equal. Thus, the coexistence of predators
observed in our model is an emergent consequence of the prey having
heterogeneous interaction strengths with the predators.

###  4.3 Emergent facilitation and trait-mediated indirect effects

Our model also demonstrates an interesting effect that one predator can have
on the other, namely facilitation [8, 38]. It happens when one predator is
better off in the presence of another even though they compete for the same
resources. We can distinguish between two levels of facilitation: Strong
facilitation is such that one predator would go extinct if the other predator
is absent but can survive when both are present. Weak facilitation is such
that the abundance of one predator is increased when the other predator is
present rather than absent. Strong facilitation was studied in [8] in the case
of adult and juvenile prey with body size differentiation. In that system,
facilitation emerges as a result of two predators specializing in different
prey stages and the prey stages competing with each other. Our system shares
the feature that two predators prefer different subpopulations of the same
prey. We found that our models exhibit both strong and weak facilitation, as
show in Fig. 10(a,b) for the case of inherent heterogeneity. Intuitively, the
consumption of one prey type by one predator allows the other prey type to
flourish, thus helping the other predator.

Figure 10: Phase diagrams showing facilitation of predators by their
competitor, and promotion of prey by predators. (a,b) Parameter regions of
strong and weak facilitation, with (a) $B$ facilitating $A$ and (b) $A$
facilitating $B$. Color represents the difference between population levels of
each predator in the presence and absence of the other. (c,d) Parameter
regions showing promotion of the prey by the addition of a second predator,
with (c) $B$ promoting $C$ and (d) $A$ promoting $C$. Color represents the
difference in population levels of the prey $C$ in the presence of both or
only one predator. Parameter values used here are
$(a_{0},a_{1},b_{0},b_{2},\rho)=(0.25,0.6,0.5,0.6,0.55)$.

We can interpret such emergent facilitation as a trait-mediated indirect
effect [39, 40]. The interaction between the predator $A$ and the prey $C$
causes a shift in the composition of $C$, which can be considered a trait of
$C$ on the population level. This trait change then affects the overall
interaction between $C$ and the other predator $B$. Therefore, in a coarse-
grained picture where the subpopulations of $C$ are lumped together (Figure
1c), it is as if one species $A$ modifies the effective interaction strength
$b_{\textrm{eff}}(\lambda)$ between the other two species $B$ and $C$, and
similarly for $B$ that modifies the interaction between $A$ and $C$. Here the
interaction modification by one predator can be either detrimental to the
other predator, or beneficial as in the cases of facilitation.

###  4.4 Multiple-predator effects and emergent promotion of prey

In addition to the facilitation between predators shown above, our model also
exhibits a surprising “promotion” of the prey population by a predator. By
this we mean that the equilibrium population of the prey is higher in the
presence of both predators than in the presence of just one. This will be an
extreme case of subadditive effects on the prey from multiple predators [41,
42], to the extent that the total effect from two predators is lower than that
from one predator alone. Figure 10(c,d) show the regions in parameter space
where there is promotion of the prey $C$ by predator $A$ and $B$ respectively.
This surprising effect arises when both predators prefer to consume the same
prey type, and the promoting predator has a stronger preference than the
other. Heuristically, if we start with only the latter predator in equilibrium
with the prey, then adding the promoting predator would push the composition
of the prey away from the preferred type. This compositional shift then
inhibits both predators and allows the prey to reach a larger population at
the new equilibrium. Mathematical conditions for such emergent promotion to
happen are given in the Appendix.

##  5 Conclusion

Heterogeneity is natural in biological systems, as each individual organism
possesses a large number of traits that are influenced one way or the other by
developmental noise or environmental variation. Every individual interacts
with other individuals and the environment somewhat differently due to its
unique combination of traits. The heterogeneity in such interactions are
generally treated in one of two ways. It is ignored in some cases where the
system is treated as being well-mixed or homogeneous, such as in the classic
Lotka-Volterra model [43, 44]. In other cases the system is treated as
“disordered” such that the interaction strengths are randomly drawn from some
probability distribution [45, 46, 47]. Here, we have addressed the
heterogeneity of biological interactions in a more structured framework that
considers subpopulations of a species with a dynamic composition. We have
demonstrated nontrivial consequences of having heterogeneous interactions. We
expect such ecosystems to exhibit stronger persistence and richer diversity.
We find that one predator can facilitate another both qualitatively and
quantitatively through trait-mediated indirect effects. And we show that a
prey species can change its population composition to benefit its own growth
under multiple predators. These effects are apparent in our model with only
one species of a small ecosystem exhibiting a single differentiating trait. In
real natural ecosystems, there is a massively higher dimensionality in both
the number of traits within a species that can vary among individuals and the
number of species that can exhibit such trait differentiation. We expect such
high dimensionality to lead to more manifestations of the effects illustrated
here.

We have focused here on biological and ecological systems, but heterogeneous
interactions can have nontrivial effects on dynamical systems in general. Here
we have seen new attractors appear where they could not in a homogeneous
system. We have seen new timescales introduced to the system that can result
in slow convergence to the equilibrium. We have even seen the transformation
of the dynamical system from a simple differential system to one which
exhibits flow switching depending on the state of the system [48]. It is less
common for traditional systems in physics to have heterogeneous interaction
constants. However, non-uniform parameters can be important in systems
involving a large number of mesoscopic components. In such systems, the
parameters of the components are often non-identical, and the shape or width
of the parameter distribution can have strong effects on the behavior of the
system. This can be seen in systems of coupled oscillators [49, 50, 51] where
the spread of the natural frequencies of the oscillators can lead to different
degrees of synchronization. The effect of inhomogeneous parameters can also be
seen in systems of colloids [52], where differences in size and mass among
large numbers of colloids can lead to intermittent behavior with periods of
ordered lattice configuration and more random motion. These are examples of
systems where differences among interacting units of a system can have
qualitative effects on the overall system behavior. As we have shown,
incorporating such differences between the units can be important for
predicting the outcome of these complex dynamics.

##  Appendix A Methods

Our model of two predators and one prey that differentiates into two types is
described by Eqs. (4a–4d) in the main text. Changing variables from $C_{1}$
and $C_{2}$ to $C\equiv C_{1}+C_{2}$ and $\lambda\equiv C_{2}/(C_{1}+C_{2})$
leads to the equations:

$\displaystyle\dot{A}$
$\displaystyle=A\big{(}\big{(}\alpha_{A1}(1-\lambda)+\alpha_{A2}\lambda\big{)}C-\beta_{A}\big{)}$
(A1a) $\displaystyle\dot{B}$
$\displaystyle=B\big{(}\big{(}\alpha_{B1}(1-\lambda)+\alpha_{B2}\lambda\big{)}C-\beta_{B}\big{)}$
(A1b) $\displaystyle\dot{C}$
$\displaystyle=C\big{(}\beta_{C}-\alpha_{CC}C-\varepsilon_{A}\big{(}\alpha_{A1}(1-\lambda)+\alpha_{A2}\lambda\big{)}A-\varepsilon_{B}\big{(}\alpha_{B1}(1-\lambda)+\alpha_{B2}\lambda\big{)}B\big{)}$
(A1c) $\displaystyle\dot{\lambda}$
$\displaystyle=\lambda(1-\lambda)\big{(}\varepsilon_{A}(\alpha_{A1}-\alpha_{A2})A+\varepsilon_{B}(\alpha_{B1}-\alpha_{B2})B\big{)}+(1-\lambda)\sigma_{1}-\lambda\sigma_{2}$
(A1d)

To simplify, we then rescale the variables by:

$\displaystyle
t\leftarrow\beta_{C}t,\;A\leftarrow\frac{\varepsilon_{A}\alpha_{CC}}{\beta_{C}}A,\;B\leftarrow\frac{\varepsilon_{B}\alpha_{CC}}{\beta_{C}}B,\;C\leftarrow\frac{\alpha_{CC}}{\beta_{C}}C$

and rename the parameters as:

$\displaystyle
a_{0}=\frac{\beta_{A}}{\beta_{C}},\;a_{1}=\frac{\alpha_{AC_{1}}}{\alpha_{CC}},\;a_{2}=\frac{\alpha_{AC_{2}}}{\alpha_{CC}},\;b_{0}=\frac{\beta_{B}}{\beta_{C}},\;b_{1}=\frac{\alpha_{BC_{1}}}{\alpha_{CC}},\;b_{2}=\frac{\alpha_{AC_{2}}}{\alpha_{CC}},\;\eta_{1}=\frac{\sigma_{1}}{\beta_{C}},\;\eta_{2}=\frac{\sigma_{2}}{\beta_{C}}$

After these transformations, we arrive at Eqs. (4a-4d) in the main text.

In the case of inherent heterogeneity, the dynamical equations for $C_{1}$ and
$C_{2}$ are:

$\displaystyle\dot{C_{1}}$
$\displaystyle=(1-\rho)\beta_{C}C-\alpha_{CC}C_{1}C-\varepsilon_{A}\alpha_{A1}C_{1}A-\varepsilon_{B}\alpha_{B1}C_{1}B$
(A2a) $\displaystyle\dot{C_{2}}$
$\displaystyle=\rho\beta_{C}C-\alpha_{CC}C_{2}C-\varepsilon_{A}\alpha_{A2}C_{2}A-\varepsilon_{B}\alpha_{B2}C_{2}B$
(A2b)

Using the same transformations and rescaling as above, these equations become
the same as Eqs. (4c-4d) in the main text with $\eta_{1}=\rho(1-C)$ and
$\eta_{2}=(1-\rho)(1-C)$.

The locations of the equilibrium points in the $A$-$B$-$C$ space are:

$\displaystyle P_{O}$
$\displaystyle=\bigg{(}0,0,0,\frac{\eta_{1}}{\eta_{1}+\eta_{2}}\bigg{)}$
$\displaystyle P_{C}$
$\displaystyle=\bigg{(}0,0,1,\frac{\eta_{1}}{\eta_{1}+\eta_{2}}\bigg{)}$
$\displaystyle P_{A}$
$\displaystyle=\bigg{(}\frac{a_{1}(1-\lambda)+a_{2}\lambda-
a_{0}}{(a_{1}(1-\lambda)+a_{2}\lambda)^{2}},0,\frac{a_{0}}{a_{1}(1-\lambda)+a_{2}\lambda},\lambda\bigg{)}$

where $\lambda$ is the solution to the cubic equation

$\displaystyle\lambda(1-\lambda)(a_{1}-a_{2})(a_{1}(1-\lambda)+a_{2}\lambda-
a_{0})+(\eta_{1}(1-\lambda)-\eta_{2}\lambda)(a_{1}(1-\lambda)+a_{2}\lambda)^{2}=0$

And $P_{B}$ has the same expression except with parameters associated with
species $B$. The location of $P_{N}$ depends on the form of $\eta_{i}$’s, but
it always lies on the line $\mathcal{L}$ given by
$\\{a_{0}A+b_{0}B=C^{*}(1-C^{*}),C=C^{*}\\}$. For fixed $\eta_{i}$’s,

$\displaystyle P_{N}$
$\displaystyle=\Bigg{(}\frac{1}{a_{0}}\bigg{(}\frac{\tilde{b}_{2}-\tilde{b}_{1}}{\tilde{a}_{1}-\tilde{a}_{2}+\tilde{b}_{2}-\tilde{b}_{1}}C^{*}(1-C^{*})-\bigg{(}\frac{\eta_{1}}{\tilde{a}_{1}-\tilde{b}_{1}}-\frac{\eta_{2}}{\tilde{b}_{2}-\tilde{a}_{2}}\bigg{)}\bigg{)},$
$\displaystyle\frac{1}{b_{0}}\bigg{(}\frac{\tilde{a}_{1}-\tilde{a}_{2}}{\tilde{a}_{1}-\tilde{a}_{2}+\tilde{b}_{2}-\tilde{b}_{1}}C^{*}(1-C^{*})+\bigg{(}\frac{\eta_{1}}{\tilde{a}_{1}-\tilde{b}_{1}}-\frac{\eta_{2}}{\tilde{b}_{2}-\tilde{a}_{2}}\bigg{)}\bigg{)},C^{*},\lambda^{*}\Bigg{)}$

where

$\displaystyle C^{*}$
$\displaystyle=\frac{\tilde{a}_{1}-\tilde{a}_{2}+\tilde{b}_{2}-\tilde{b}_{1}}{\tilde{a}_{1}\tilde{b}_{2}-\tilde{a}_{2}\tilde{b}_{1}}$
$\displaystyle\lambda^{*}$
$\displaystyle=\frac{\tilde{a}_{1}-\tilde{b}_{1}}{\tilde{a}_{1}-\tilde{a}_{2}+\tilde{b}_{2}-\tilde{b}_{1}}$

with $\tilde{a}_{1}=a_{1}/a_{0}$, $\tilde{a}_{2}=a_{2}/a_{0}$,
$\tilde{b}_{1}=b_{1}/b_{0}$, and $\tilde{b}_{2}=b_{2}/a_{0}$. The solutions
for the inherent case are recovered by making the substitutions for $\eta_{1}$
and $\eta_{2}$ above, and for the adaptive case by removing the $\eta$ terms.

The stability of these equilibrium points is determined by the following
rules, with the physical region $\mathcal{W}\equiv\\{A\geq 0,B\geq 0,0\leq
C\leq 1,0\leq\lambda\leq 1\\}$:

  * •

If only the points $P_{O}$ and $P_{C}$ are in $\mathcal{W}$, then $P_{C}$ is
stable.

  * •

If only $P_{O}$, $P_{C}$ and $P_{A}$ (or $P_{B}$) are in $\mathcal{W}$, then
$P_{A}$ (or $P_{B}$) is stable.

  * •

If $P_{O}$, $P_{C}$, $P_{A}$, $P_{B}$ are all in $\mathcal{W}$ but not
$P_{N}$, then one of $P_{A}$ and $P_{B}$ with a smaller $C^{*}$ is stable.

  * •

If $P_{N}$ is in $\mathcal{W}$ then it is stable, or there is a small limit
cycle enclosing it that is stable.

For completeness we list here the exact conditions for the regions of
coexistence between the predators. The expression is simplest for the case of
instant ($\kappa\to\infty$) adaptive heterogeneity:

$\displaystyle\big{(}(\tilde{a}_{2}<\tilde{b}_{2})\land(\tilde{b}_{1}<\tilde{b}_{2})\land(C^{*}(P_{N})<1)\big{)}\lor\big{(}(\tilde{a}_{2}>\tilde{a}_{1})\land(\tilde{b}_{1}>\tilde{a}_{1})\big{)}$

For the other kinds of heterogeneity, the condition for coexistence is given
by

$\displaystyle(A^{*}(P_{N})>0)\land(B^{*}(P_{N})>0)\land(1>C^{*}(P_{N})>0)\land(1>\lambda^{*}(P_{N})>0)$

In the case of inherent heterogeneity, the requirement for strong facilitation
of predator $B$ by predator $A$ is that $P_{C}$ is stable in the system when
$A=0$ and $P_{N}$ is stable when $A$ is present. For weak facilitation it is
required that $P_{B}$ and $P_{N}$ are stable, respectively, and
$B^{*}(P_{N})>B^{*}(P_{B})$. In order to see emergent promotion, in addition
to $P_{N}$ being stable, we need
$a_{1}/a_{0}>b_{1}/b_{0}>b_{2}/b_{0}>a_{2}/a_{0}$ for $A$ to promote $C$, and
$b_{1}/b_{0}>a_{1}/a_{0}>a_{2}/a_{0}>b_{2}/b_{0}$ for $B$ to promote $C$.

## References

  * [1] Daniel I Bolnick, Priyanga Amarasekare, Márcio S Araújo, Reinhard Bürger, Jonathan M Levine, Mark Novak, Volker H W Rudolf, Sebastian J Schreiber, Mark C Urban, and David A Vasseur.  Why intraspecific trait variation matters in community ecology.  Trends Ecol Evol, 26:183–192, 2011.
  * [2] Andrew Siefert and et al.  A global meta-analysis of the relative extent of intraspecific trait variation in plant communities.  Ecology Letters, 18(12):1406–1419, 2015.
  * [3] Amy B Forsythe, Troy Day, and William A Nelson.  Demystifying individual heterogeneity.  Ecology Letters, 24(10):2282–2297, 2021.
  * [4] Daniel I. Bolnick, Richard Svanbäck, James A. Fordyce, Louie H. Yang, Jeremy M. Davis, C. Darrin Hulsey, and Matthew L. Forister.  The ecology of individuals: Incidence and implications of individual specialization.  The American Naturalist, 161(1):1–28, 2003.
  * [5] Garrett Hardin.  The competitive exclusion principle.  Science, 131(3409):1292–1297, 1960.
  * [6] Peter Chesson.  Mechanisms of maintenance of species diversity.  Annual Review of Ecology and Systematics, 31(1):343–366, 2000.
  * [7] Dang H. Nguyen and George Yin.  Coexistence and exclusion of stochastic competitive lotka-volterra models.  Journal OF Differential Equations, 262(3, 1):1192–1225, 2017.
  * [8] A. M. De Roos, T. Schellekens, T. Van Kooten, and L. Persson.  Stage-specific predator species help each other to persist while competing for a single prey.  Proceedings of the National Academy of Sciences, 105(37):13930–13935, 2008.
  * [9] David W Frayer and Milford H Wolpoff.  Sexual dimorphism.  Annual Review of Anthropology, 14(1):429–473, 1985.
  * [10] Douglas J Emlen, John Hunt, and Leigh W Simmons.  Evolution of sexual dimorphism and male dimorphism in the expression of beetle horns: phylogenetic evidence for modularity, evolutionary lability, and constraint.  the american naturalist, 166(S4):S42–S68, 2005.
  * [11] Martin Stevens and Sami Merilaita.  Animal camouflage: current issues and new perspectives.  Philosophical Transactions of the Royal Society B: Biological Sciences, 364(1516):423–427, 2009.
  * [12] Yang Niu, Hang Sun, and Martin Stevens.  Plant camouflage: Ecology, evolution, and implications.  Trends Ecol Evol, 33(8):608–618, 2018.
  * [13] Mohamed Noor, Robin Parnell, and Bruce Grant.  A reversible color polyphenism in american peppered moth (biston betularia cognataria) caterpillars.  PloS one, 3:e3142, 2008.
  * [14] Harry Mickalide and Seppe Kuehn.  Higher-order interaction between species inhibits bacterial invasion of a phototroph-predator microbial community.  Cell Systems, 9(6):521–533, 2019.
  * [15] P.E.H. Minichin, M.R. Thorpe, and J.F. Farrar.  Short-term control of root: shoot partitioning.  Journal of Experimental Botany, 45(5):615–622, 1994.
  * [16] J. J. Gedroc, K. D. M. McConnaughay, and J. S. Coleman.  Plasticity in root/shoot partitioning: Optimal, ontogenetic, or both?  Functional Ecology, 10:44–50, 1996.
  * [17] Arnold J. Nicholson and V. A. Bailey.  The balance of animal populations.—part i.  Journal of Zoology, 105:551–598, 1935.
  * [18] Alan Hastings and Susan Harrison.  Metapopulation dynamics and genetics.  Annual review of Ecology and Systematics, 25(1):167–188, 1994.
  * [19] Alan Hastings, Karen C. Abbott, Kim Cuddington, Tessa Francis, Gabriel Gellner, Ying-Cheng Lai, Andrew Morozov, Sergei Petrovskii, Katherine Scranton, and Mary Lou Zeeman.  Transient phenomena in ecology.  Science, 361(6406):eaat6412, 2018.
  * [20] Andrew Morozov, Karen Abbott, Kim Cuddington, Tessa Francis, Gabriel Gellner, Alan Hastings, Ying-Cheng Lai, Sergei Petrovskii, Katherine Scranton, and Mary Lou Zeeman.  Long transients in ecology: Theory and applications.  Physics of Life Reviews, 32:1–40, 2020.
  * [21] Peter A. Abrams.  Adaptive foraging by predators as a cause of predator-prey cycles.  Evolutionary Ecology, 6, 1992.
  * [22] Peter A. Abrams.  Implications of flexible foraging for interspecific interactions: lessons from simple models.  Functional Ecology, 24(1):7–17, 2010.
  * [23] Andrew P. Beckerman, Maria Uriarte, and Oswald J. Schmitz.  Experimental evidence for a behavior-mediated trophic cascade in a terrestrial food chain.  Proceedings of the National Academy of Sciences, 94(20):10735–10738, 1997.
  * [24] Anurag A. Agrawal.  Induced responses to herbivory and increased plant performance.  Science, 279(5354):1201–1202, 1998.
  * [25] Charlotte T. Lee, Tom E. X. Miller, and Brian D. Inouye.  Consumer effects on the vital rates of their resource can determine the outcome of competition between consumers.  The American Naturalist, 178(4):452–463, 2011.
  * [26] Charlotte T. Lee and Brian D. Inouye.  Mutualism between consumers and their shared resource can promote competitive coexistence.  The American Naturalist, 175(3):277–288, 2010.
  * [27] Rafal Goebel, Ricardo G. Sanfelice, and Andrew R. Teel.  Hybrid dynamical systems.  IEEE Control Systems Magazine, 29(2):28–93, 2009.
  * [28] Aleksej F. Filippov.  Differential equations with discontinuous righthand sides.  In Mathematics and Its Applications, 1988.
  * [29] Edo Kussell and Stanislas Leibler.  Phenotypic diversity, population growth, and information in fluctuating environments.  Science, 309(5743):2075–8, 2005.
  * [30] Matina C Donaldson-Matasci, Carl T Bergstrom, and Michael Lachmann.  When unreliable cues are good enough.  The American naturalist, 182(3):313–27, 2013.
  * [31] G. E. Hutchinson.  The paradox of the plankton.  The American Naturalist, 95(882):137–145, 1961.
  * [32] Gyorgy Barabas, Matthew J. Michalska-Smith, and Stefano Allesina.  The effect Of intra- and interspecific competition on coexistence in multispecies communities.  American Naturalist, 188(1):E1–E12, 2016.
  * [33] Peter Chesson.  Updates on mechanisms of maintenance of species diversity.  Journal of Ecology, 106(5):1773–1794, 2018.
  * [34] Jonathan M Levine, Jordi Bascompte, Peter B Adler, and Stefano Allesina.  Beyond pairwise mechanisms of species coexistence in complex communities.  Nature, 546(7656):56–64, 2017.
  * [35] Serguei Saavedra, Rudolf P Rohr, Jordi Bascompte, Oscar Godoy, Nathan JB Kraft, and Jonathan M Levine.  A structural approach for understanding multispecies coexistence.  Ecological Monographs, 87(3):470–486, 2017.
  * [36] David Tilman.  Resource Competition and Community Structure.  Princeton University Press, 1982.
  * [37] Ellen van Velzen.  Predator coexistence through emergent firness equalization.  Ecology, 101:e0295, 2020.
  * [38] André M. de Roos.  Effects of life history and individual development on community dynamics: A review of counterintuitive consequences.  Ecological Research, 35(6):930–946, 2020.
  * [39] J. Timothy Wootton.  The nature and consequences of indirect effects in ecological communities.  Annual Review of Ecology and Systematics, 25(1):443–466, 1994.
  * [40] Takayuki Ohgushi, Oswald Schmitz, and Robert D. Holt, editors.  Trait-Mediated Indirect Interactions: Ecological and Evolutionary Perspectives.  Cambridge University Press, 2012.
  * [41] Andrew Sih, Göran Englund, and David Wooster.  Sih a, englund g, wooster d.. emergent impacts of multiple predators on prey.  Trends Ecol Evol, 13:350–355, 1998.
  * [42] Michael W. McCoy, Adrian C. Stier, and Craig W. Osenberg.  Emergent effects of multiple predators on prey survival: the importance of depletion and the functional response.  Ecology Letters, 15(12):1449–1456, 2012.
  * [43] Alfred J. Lotka.  Contribution to the theory of periodic reactions.  The Journal of Physical Chemistry, 14(3):271–274, 1910.
  * [44] Alfred J. Lotka.  Elements of physical biology.  Nature, 116:461, 1925.
  * [45] Robert M. May.  Will a large complex system be stable?  Nature, 238:413, 1972.
  * [46] Dominique Gravel, François Massol, and Mathew A. Leibold.  Stability and complexity in model meta-ecosystems.  Nature Communications, 7:12457, 2016.
  * [47] Matthieu Barbier, Claire de Mazancourt, Michel Loreau, and Guy Bunin.  Fingerprints of high-dimensional coexistence in complex ecosystems.  Phys. Rev. X, 11:011009, 2021.
  * [48] Albert C.J. Luo.  A theory for flow switchability in discontinuous dynamical systems.  Nonlinear Analysis: Hybrid Systems, 2(4):1030–1061, 2008.
  * [49] Steven H Strogatz.  From kuramoto to crawford: exploring the onset of sychronization in populations of coupled oscillators.  Physica D: Nonlinear Phenomena, 142:1–20, 2000.
  * [50] Henning Ulrichs, Andreas Mann, and Ulrich Parlitz.  Synchronization and chaotic dynamics of coupled mechanical metronomes.  Chaos, 19:043120, 2009.
  * [51] Zachary Jackson and Kurt Wiesenfeld.  Emergent, linked traits of fluctuation feedback systems.  Phys. Rev. E, 104:064216, 2021.
  * [52] Guram Gogia and Justin C. Burton.  Emergent bistability and switching in a nonequilibrium crystal.  Phys. Rev. Lett., 119:178004, 2017.