1. Introduction
Mutualisms are interactions between two or more species which are beneficial for all, the benefits being materialized in the form of an enhanced capacity to survive, grow or reproduce (Holland and Bronstein [
1]). If all species involved can survive in isolation, then the mutualism is called facultative, while if the species can only survive in association, the mutualism is called obligate. Historically speaking, antagonistic interactions such as predation, parasitism and competition have received far greater attention than mutualisms, this lack of attention stemming, perhaps, from the adherence of the mainstream ecological thinking to the Darwinian paradigm of the survival of the fittest. The fact that even the comparatively simple model of mutualism of Gause and Witt [
2], obtained from a Lotka–Volterra predator–prey model with logistic growth by inverting the sign of the prey loss term, exhibits the unboundedness of solutions provided that the strength of the mutualistic interaction exceeds a certain threshold helped cement the idea that models of mutualism are prone to having built-in flaws.
The model of Gause and Witt relies on two linearity assumptions. First, the effect of each mutualist on the per capita growth rate of its partner is assumed to be directly proportional to its density. Second, the intraspecific competition term is assumed to have a linear density dependence, which leads to a logistic within-species dynamics. However, even though both assumptions are gross simplifications of reality, the blame on the unreasonable behavior of solutions is usually placed solely on the linearity of the functional response of each species toward its mutualistic partner, which is perhaps due to to the ubiquity of the logistic growth in theoretical ecology that leads to its widespread acceptance.
A more balanced view, namely that relaxing either of the linearity assumptions of [
2] is enough to prevent unrealistic behavior of the solutions, is taken by Moore et al. [
3]. In [
3], a linear mutualistic contribution is employed, coupled with a modified intraspecies competition term, being observed that a higher-powered intraspecific competition term is enough to balance this linear mutualistic contribution, regardless of the strength of the mutualism, while lower-powered order intraspecific competition terms tend to destabilize populations. A comprehensive analysis of the relative benefits brought on by mutualisms, the number of equilibria and their stability, primarily from a numerical viewpoint, has subsequently been performed in [
3], a set of predictions regarding the outcomes of mutualistic interactions being made even when births and deaths are modeled as separate processes with different intraspecific density dependences.
The first of the two-species models of mutualism with
-logistic birth and death functions and the linear per capita mutualistic contribution that were proposed and discussed in [
3] is as follows
In the above model,
,
, are the birth functions and
,
, are the death functions, which are written as a combination of a density-independent term
(or
) and a density-dependent term
(or
). Such a function is termed in [
3] as being an
accelerating function of density if the corresponding exponent is >1 and a
decelerating function if the corresponding exponent is <1. The constants
and
quantify the mutualistic support received by each species from the other.
Note that there is significant evidence for nonlinear per capita growth rates over a wide range of species and taxa, as presented in Sibly et al. [
4], suggesting that in each taxa, there are far more instances of decelerating functions (on average, 78%, higher in fish than in mammals, birds or insects). Nonlinear density dependence is often considered, from an evolutionary perspective, in connection with life history strategies. It is then theoretized that populations with high intrinsic growth rates (referred to as
r-selected populations) normally exhibit decelerating density dependence, since their survival probability diminishes at relatively low densities, while populations with low intrinsic growth rates (referred to as
K-selected populations) normally exhibit accelerating density dependence, since their survival probability diminishes at relatively high densities.
Regardless of that and of other specific evidence presented for both controlled laboratory populations with simple life histories (
Drosophila, Gilpin and Ayala [
5],
Daphnia, Smith [
6], Smith and Cooper [
7]) and species with more complex histories (Coulson et al. [
8]), models that relax the linear density dependence assumption have rarely been used to describe mutualistic interactions (Ribeiro et al. [
9], Wong [
10]).
A model that is formally related to (
1) although, from a modelling perspective, it does not explicitly distinguish between births and deaths as separate phenomena, has been introduced and analyzed in Garcia-Algarra et al. [
11] together with its multispecies generalization. Of primary concern in [
11] was assessing the resilience of the system when confronted with external perturbations, rather than studying the stability of equilibria, being determined via extensive numerical simulations following a stochastic approach that the system resilience is a function of the network structure, specifically of its nestedness.
Another model discussed in the Appendix A of [
3] allows for the saturation of mutualistic support at higher densities, having the form
In (
2),
represents the maximal benefit received by species
i from its mutualistic partner, while
is the half-saturation constant, measuring how fast the saturation process occurs, the birth and death functions retaining their significance. In what follows, we shall assume that both (
1) and (
2) describe facultative mutualisms and consequently
,
.
A consistent first step toward discussing well-posedness and stability considerations for mutualistic systems is establishing the boundedness of solutions, which has the added benefit of ruling out May’s “orgy of mutual benefaction”. In this regard, an approach toward establishing boundedness results (and consequently ruling out unrealistic model behavior) for mutualistic systems with an arbitrary number of species in terms of reproductive parameters computed at high population densities has been presented and illustrated in Georgescu et al. [
12]. In a certain significant particular case, characterizing the boundedness (or lack thereof) of a mutualistic interaction was reduced to computing the spectral abscissa of an associated limiting interaction matrix, yielding a hands-on tool to detect unrealistic model behavior.
The stability of co-existence equilibria is, however, not discussed in [
12] (in fact, it is outside the scope of that paper), and consequently, the impact of mutualistic benefits toward stabilizing ecological interactions cannot be completely quantified. A discussion in this regard, using a conceptually similar but notationally different approach, has been laid out in Maxin et al. [
13], and the interplay between mutualistic benefits and weak or strong Allee effects is also discussed via suitably modifying the per capita growth rates. It has been observed in [
13] that mutualisms can overcome even strong Allee effects and make a unique co-existence equilibrium globally asymptotically stable. Apart from that, Ref. [
13] allows for higher flexibility when modelling mutualistic interactions (for instance, one of the species can benefit from an increase in the reproductive rates, while the other can benefit from a reduction in mortality), accounting for the plethora of mechanisms on which real-life mutualisms are based.
A different approach toward discussing the stability of equilibria for certain models of mutualism in common use, based on Lyapunov’s second method, has been followed in [
14,
15]. Particularly, one of the models of concern in Georgescu et al. [
15] is the model with
-logistic growth (although the corresponding exponent is denoted by
p rather than by
) and linear functional response below
being shown that if
, then the system (
3) (which reduces in this case to Gause and Witt’s model) has a unique co-existing equilibrium if and only if
, while if
, then the system (
3) has a unique co-existing equilibrium regardless of the value of
. Also, this co-existing equilibrium is globally asymptotically stable in
whenever it exists. Note that no attempt at discussing the multiplicity or stability of the co-existing equilibrium is made in [
15] for the case
.
This relates to the findings of Moore et al. [
3] (accelerating density dependence for the intraspecific competition term tends to stabilize the system, although with comparatively less mutualistic benefit, defined as the difference in density in the presence and absence of the mutualist partner, respectively) and leads to the idea that, for a two-species interaction and under suitable conditions, the co-existence equilibrium is prone to be globally stable whenever it is unique. In fact, apart from characterizing the stability of the co-existence equilibria, [
3] is also concerned with their multiplicity, as (
1) has at most two of those. At this point, one should perhaps note that [
3] counts all equilibria together, including one trivial and two semi-trivial, so the reader interested only in co-existence equilibria should always subtract three from the equilibria count presented in [
3].
In what follows, we shall attempt to further establish an analytic framework to confirm the numerical findings of [
3] while also addressing the additional complications brought on by modeling birth and death as separate phenomena.
In
Section 2, we introduce a generic framework for studying the stability of two-species models of mutualism, sufficient conditions for the existence and global stability of the co-existence equilibria being obtained in terms of reproductive numbers computed at high population densities. In
Section 3, we particularize these findings for the model with linear per capita mutualistic benefits (
1) while also providing alternative proofs for the existence and uniqueness of the co-existence equilibrium and further local stability considerations. In
Section 4, we perform a parallel discussion for the model with saturating mutualistic benefits (
2), observing that the requirements on the birth and death functions are less demanding, since there is potentially less unboundedness potential from the mutualistic contributions to be mitigated.
While
Section 3 and
Section 4 are mostly concerned with global stability results,
Section 5 aims at dealing with the loss of global stability due to either unboundedness or multiple co-existence equilibria. It is then observed that a certain dichotomy property holds, which leads to one of the co-existence equilibria retaining an attractivity property once the boundedness of solutions is established. Finally,
Section 6 presents a discussion on the biological implications of our results along with concluding remarks.
2. Generic Boundedness and Stability Results via Reproductive Numbers
To set up an unified treatment of (
1) and (
2) from a stability viewpoint, let us consider the following abstract model
in which
g is a positive, continuously differentiable function on
and
d is a positive, continuously differentiable function on
under the following consistency assumptions:
- (C1)
is decreasing as a function of for each and is decreasing as a function of for each ;
- (C2)
and are ultimately decreasing for all ;
as well as the following self-limiting growth assumption:
- (L)
There are
,
such that
and the following mutualistic assumption:
- (M)
is increasing in for each and is increasing in for each .
Note that for concrete models, the
and
values are not necessarily uniquely defined, but they are only defined up to functions of
. For this very reason,
values may not represent the exact forms of the per capita mutualistic contributions (although they are certainly related to it) and
values may not represent the exact forms of the per capita removal rates. Assumptions (C1) and (C2) (which, although similar looking, have different scopes) essentially describe the “relative saturation” of the mutualistic benefits (as measured via the benefit-to-removal quotient) at high population densities. The logistic-like assumption (L) ensures that each species in isolation has self-limiting growth and that the mutualistic interaction is facultative, since the per capita growth rate of each species is strictly positive at small densities. Consequently, there is no danger of species extinction and we shall not concern ourselves with the stability of the trivial equilibrium, as this particular equilibrium is unstable by design. Let us then denote
and observe that
. Subsequently, using also (M), it is seen that the semi-trivial equilibria
and
are also a priori unstable, so we shall not concern ourselves with the stability of those, either. Assumption (C2) ensures that the following limits are well-defined
The parameters and can be thought of as reproductive numbers, much like the basic reproductive number of Mathematical Epidemiology (notice that its definition follows the same paradigm of new generation times the inverse of removal, with g values in place of F and d values in place of V). They are, however, computed at high population densities rather than in near-extinction conditions, since now the well-posedness is of primary concern (represented by the boundedness of solutions, which is essentially a non-local property), rather than the persistence or extinction of a certain compartment, as it happens in Mathematical Epidemiology. Note also that in the multiple-species context of our paper, one needs the reproductive quantities and computed on a per-species basis, rather than having a single reproductive number to rule them all, as it is the case in Mathematical Epidemiology.
We start by establishing the boundedness of solutions both from above and below. To this purpose, let us state additional assumptions.
- (B)
There are such that and .
- (D)
, .
While (B) is a boundedness assumption, (D) essentially tells us that values act as limiting terms. In fact, , are usually increasing and unbounded in concrete situations, which easily implies (D). Also, although and may not represent the exact forms of the removal rates, as outlined above, assumption (D) essentially ensures that no individual of any species exists in perpetuity due to the removal rates becoming negligible at higher population densities.
The next result establishes that under the assumptions listed herein, there is no “orgy of mutual benefaction”, as termed in May [
16], and the solutions of (
4) remain bounded.
Lemma 1. Assume that (M), (B), (D) and (C2) hold. Then, the solutions of (4) are ultimately uniformly bounded from above. Proof. Suppose that for a given
t, one has that
. Then,
Since
it follows that
for some given
if
is large enough. If
instead, then one may infer via a parallel argument that
if for some given
if
is large enough.
Trouble is, h may not necessarily be differentiable for all (specifically, it may not be differentiable for values of t such that ). However, one may infer from the above that for some when is large enough, which implies that is ultimately bounded from above. Consequently, so are both and . □
Notice that (C2), in fact, is used in our abstract framework only to ensure that , are well-defined for all . If the values of , can be obtained via a different argument, such as direct limiting, then (C2) becomes redundant.
To determine co-existence equilibria, simply proving boundedness is definitely not enough. One should also establish the long-term persistence of all species, ensuring that none of them become extinct, which is the purpose of the next result.
Lemma 2. Assume that (M) and (L) hold. Then, the solutions of (4) starting in are ultimately bounded from below by positive constants (i.e., (4) is uniformly persistent). Proof. From (M) and (L), it follows that
as long as
stays in some interval
,
, from which
and the species
is uniformly persistent. The uniform persistence of
holds by a similar argument. □
Note that if the hypotheses of both Lemmas 1 and 2 hold, then the solutions of (
4) ultimately reach a compact set of the positive quadrant and there exists (at least) a co-existence equilibrium
, by Theorem 2.8.6 of Bhatia and Szegö [
17]. The uniqueness and stability of
, however, remain to be discussed. While the former is not immediate, the latter can be obtained ruling out periodic orbits via Poincaré–Bendixson theory.
Theorem 1. Assume that (M), (C1), (C2), (L), (B), (D) hold and the co-existence equilibrium is unique. Then, is globally asymptotically stable.
Proof. The global stability of the co-existence equilibrium
follows from Lemma 1 by ruling out periodic solutions with the help of Dulac criterion. The Dulac function corresponding to (
4) is
, which leads to
Consequently, is globally asymptotically stable. □
From the above, we see that indeed, as numerically observed in Moore et al. [
3], the co-existence equilibrium is prone to be globally stable once it exists and is unique; of course, we are still left with the burden of proving its uniqueness. Note that the above Theorem 1 relies on a strictly two-dimensional argument, so it applies to mutualistic interactions between two species only, while Lemmas 1 and 2 remain valid,
mutatis mutandis, for mutualistic interactions between an arbitrary number of species.
We need now to find conditions for the uniqueness of the co-existence equilibrium. Although we shall later provide direct proofs, let us start by introducing notions and notations that will help us discuss this matter in a larger, more general context, of monotone dynamical systems. Given , , , we shall write if for , if and if for . Also, given , we shall denote by its associated Jacobian matrix computed at x and say that h is cooperative if is a Metzler matrix (all its off-diagonal components are non-negative) for all .
Let us now state an additional result that, apart from establishing the uniqueness of the co-existence equilibrium, is applicable to mutualistic interactions between an arbitrary number of species, the main point being that one may use the theory of monotone dynamical systems to establish the uniqueness and global stability of the co-existence equilibrium for n-dimensional systems under an additional assumption, namely a monotonicity condition involving the Jacobian.
Theorem 2 ([
18]).
Consider the systemunder the following assumptions- (H1)
F is cooperative;
- (H2)
;
- (H3)
when .
F being the vector-valued function defined by . If (5) possesses an equilibrium point in , then this equilibrium point is unique in and also globally asymptotically stable in . Regarding the significance of assumptions mentioned in the above theorem, if understood as being applied to an
n-species model of mutualism, while (H1) attests that (
5) describes indeed a mutualism, (H2) specifies that this mutualism is a facultative one, the main difficulty of using this result in concrete situations being to verify (H3).
Note that for our model (
4), (H1) corresponds to (M) and (H2) corresponds to
,
, inequalities which, as previously mentioned, are an outcome of (L). Also,
We shall now proceed with using Theorem 1 (together with its associate framework) and Theorem 2 to discuss the dynamics of the model with linear functional response (
1) and of the model with saturating functional response (
2), respectively, delineating the corresponding results in two distinct sections.
5. Dealing with Unboundedness and Multiple Co-Existence Equilibria
So far, we have been more concerned with the boundedness of solutions and global stability results (which assume the uniqueness of the co-existence equilibrium), which are associated with accelerating density dependences. Let us now focus more on unboundedness and decelerating density responses. In this regard, let us first state certain notions and notations for future reference.
Given a square matrix A, the directed graph associated to it is a directed graph with vertices such that there exists an arc leading from j to k, , if and only if . A will then be called irreducible if is strongly connected, i.e., any two distinct vertices are joined by an oriented path.
Let us first state the following dichotomy result.
Theorem 8 ([
19]).
Consider the system , where is a continuously differentiable map. Assume that- (D1)
h is cooperative on and is irreducible for any ;
- (D2)
and for all with , ;
- (D3)
.
Let be the maximally extended solution of initiating at . Then, either
- 1.
For any , , or alternatively
- 2.
There exists with such that for any y with , . Moreover, for any , .
Let us turn back our attention to (
4). Note that (M) implies that
h is cooperative. If both (M) and (L) hold then, as previously seen, (
4) is uniformly persistent; as far as we restrict our initial data to
, the solutions will not reach the boundaries of the first quadrant, so irreducibility is needed only on
. Let us then assume the stronger version of (M) below.
- (SM)
, on .
Then, the required irreducibility property does hold, since the off-diagonal elements of the Jacobian are strictly positive. Also, (D2) is obvious and
so (D3) also holds. One then obtains the following result.
Theorem 9. Assume that (SM) and (L) hold. Let be the maximally extended solution of (4) initiating at . Then either - 1.
For any , , or alternatively
- 2.
There exists a co-existence equilibrium of (4) such that for any y with , . Moreover, for any , .
Theorem 9 confirms the numerical findings of Moore et al. [
3] (pp. 193–195), namely the fact that adding mutualisms to populations with decelerating density dependence changes the dynamics in one of the following two ways.
- 1.
It completely destabilizes the dynamics, creating (all-round) unbounded population growth.
- 2.
It creates multiple interior equilibria, the “smallest” one enjoying a certain attractivity property.
Figure 4 presents the first type of outcome (unboundedness) for the model (
1) with linear density dependences when the reverse of the stability condition (
7) is satisfied (in some sense,
Figure 4 is a counterpart of
Figure 2, showcasing the role of the stability condition). For
Figure 4,
,
,
,
.
Figure 5 presents the second type of outcome (multiple co-existence equilibria) for the model (
1) with mixing density dependences (accelerating and decelerating). For
Figure 5,
,
,
,
,
,
,
,
. Notice that among the two co-existence equilibria
and
, the “smaller” one (
) is the stable one, and the “larger” one (
) is an unstable saddle.
Strictly speaking, the findings of Theorem 9 are valid for the fully accelerating density dependence case too, but they definitely are not optimal for this case. Specifically, the first situation is voided (since all solutions are bounded) and the attractivity property provided by the second situation is (much) weaker than the global stability property that is actually valid. Theorem 9 is then supposed to provide insights for the case in which at least one of the density dependences is decelerating and/or there are multiple co-existence equilibria. Regarding the latter situation, our simulations (admittedly not extensive) seem to confirm the findings of [
3], namely that there are at most two co-existence equilibria.
Let us now think again about the concrete models (
1) and (
2). Note that if (B), (D) and (C2) hold in addition to (SM) and (L), then the solutions of (
4) are ultimately uniformly bounded, as seen in Lemma 1, and consequently, the first possible outcome of Theorem 9 is ruled out. It is easy to see (and in part, it has already been seen, in fact) that (D), (SM) and (L) hold for both (
1) and (
2) whenever
. Also, as previously mentioned, (C2) becomes redundant since one can compute
and
via direct limiting. Consequently, one can think of (B) as a proxy to ultimate uniform boundeness and to the second outcome (the attractivity property) of Theorem 9.
As seen before, for the model with saturated mutualistic contribution (
2),
for all
whenever
,
, so (B) holds, and consequently, the first outcome (unboundedness) of Theorem 9 is ruled out. That is, (
2) can exhibit a loss of global stability for a co-existence equilibrium via the existence of other co-existence equilibria but not via unboundedness, and the attractivity property of Theorem 9 holds for (
2) whenever
,
.
For the model with linear mutualistic contribution (
1), the picture becomes more complicated. Note that
for all
if
,
, so only one of
,
and
, respectively, needs to be
in order to ensure that (B) holds, together with the attractivity property of Theorem 9. However, this condition is sufficient and not necessary, as one needs
and
to be both less than 1, not 0. Note also that
whenever
,
, which is a situation that is prone to unboundedness. Consequently, unlike (
2), the model with linear mutualistic contribution (
1) can exhibit a loss of global stability for a co-existence equilibrium both via the existence of other co-existence equilibria and via unboundedness.
6. Conclusions
In this paper, we have proposed a generic stability framework for studying the dynamics of two-species mutualisms, establishing sufficient conditions for the boundedness of solutions and stability of positive equilibria in terms of reproductive numbers computed at high population densities. This framework was used in order to provide a unifying perspective for the treatment of two models with linear and saturating mutualistic contributions, respectively, introduced in Moore et al. [
3], which represented the motivation for our research.
From a biological perspective, proving the boundedness of solutions is an important first step in establishing the plausibility of a model of mutualism. Furthermore, gaining insight about the stability of mutualisms can help toward the conservation of endangered species via the judicious release of mutualistic partners into the environment. In this regard, the models proposed and studied primarily from a numerical viewpoint in [
3] afford for greater flexibility, as the birth and death are modeled as separate processes with possibly distinct density dependences. Furthermore, the models of [
3] afford for greater realism, too, as the evidence for nonlinear per capita growth rates over a wide range of species and taxa presented in Sibly et al. [
4] calls for parsimony over the use of the classic logistic growth rate in spite of its ubiquity.
Mutualisms are also known to modulate biodiversity via a plethora of mechanisms, being argued in Gómez and Verdú [
20] that mutualistic partnerships are a major, potent driver of morphological and ecological diversification. In this regard, Chomicki et al. [
21] catalogued several ways through which mutualisms are able to enhance species diversification. While partner shifts can drive divergent selection, especially in plant–pollinator mutualisms, a species being able to access previously inaccessible resources via the mediation of the mutualistic partner brings an increase in ecological opportunity and a possible niche expansion. Mutualisms can increase range sizes, decreasing extinction rates and may accelerate specificity and drive lineage-specific incompatibility via deleterious mutations of the symbiont that are not compensated by host mutations. It then becomes of paramount importance to study the mechanisms of mutualism-facilitated stabilization and their robustness and resilience in order to promote species richness and coevolution. Our paper is just a step in this direction.
As numerically predicted in [
3], our results drastically delineate between accelerating and decelerating density dependences. While the former ensure the global stability of co-existence equilibria for both linear and saturating mutualistic contributions, the latter are prone to cause a loss of global stability for both types of contributions and even unboundedness for the linear ones.
There is a striking difference between the behavior of models with linear and saturating mutualistic contributions, too, as the requirements for boundedness of solutions are much less demanding for the model with saturating contributions. Simply put, this happens since in that case, there is less potential for unboundedness from the mutualistic contribution to mitigate.
When the sufficient conditions for the global stability of the co-existence equilibrium are not met, it is seen that a dichotomy property holds. Essentially, this property affirms that either unboundedness is a global phenomenon for the model of concern (i.e., any solution not starting in the origin grows unbounded) or there is a “smaller" co-existence equilibrium that still retains a certain attractivity property.
While one of our global stability results, Theorem 1, relies on a strictly two-dimensional argument, the rest of this paper does not, the other results being amenable to be extended to mutualistic models with an arbitrary number of species. We plan to do so in a subsequent paper.