Next Article in Journal
Testing Equality of Several Distributions at High Dimensions: A Maximum-Mean-Discrepancy-Based Approach
Previous Article in Journal
Adaptive Sliding Mode Attitude Tracking Control for Rigid Spacecraft Considering the Unwinding Problem
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Facultative Mutualisms and θ-Logistic Growth: How Larger Exponents Promote Global Stability of Co-Existence Equilibria

1
Department of Mathematics, Technical University of Iaşi, Bd. Copou 11, 700506 Iaşi, Romania
2
School of Economics and Management, Changzhou Institute of Technology, Changzhou 213032, China
3
School of Innovation and Entrepreneurship, Changzhou Institute of Technology, Changzhou 213032, China
*
Author to whom correspondence should be addressed.
Mathematics 2023, 11(20), 4373; https://doi.org/10.3390/math11204373
Submission received: 21 September 2023 / Revised: 15 October 2023 / Accepted: 18 October 2023 / Published: 21 October 2023
(This article belongs to the Section Mathematical Biology)

Abstract

:
We investigate the stability of co-existence equilibria for two-species models of facultative mutualism for which birth and death are modeled as separate processes, with possibly distinct types of density dependence, and the mutualistic contributions are either linear or saturating. To provide a unifying perspective, we first introduce and discuss a generic stability framework, finding sufficient stability conditions expressed in terms of reproductive numbers computed at high population densities. To this purpose, an approach based on the theory of monotone dynamical systems is employed. The outcomes of the generic stability framework are then used to characterize the dynamics of the two-species models of concern, delineating between decelerating (lower-powered) and accelerating (higher-powered) density dependences. It is subsequently seen that accelerating density dependences promote the stability of co-existence equilibria, while decelerating density dependences either completely destabilize the system via promoting the unboundedness of solutions or create multiple co-existence equilibria.

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
1 x 1 d x 1 d t = b 1 μ 1 x 1 η 1 d 1 + ν 1 x 1 θ 1 + β 1 x 2 1 x 2 d x 2 d t = b 2 μ 2 x 2 η 2 d 2 + ν 2 x 2 θ 2 + β 2 x 1 .
In the above model, b i μ i x i η i , i = 1 , 2 , are the birth functions and d i + ν i x i η i , i = 1 , 2 , are the death functions, which are written as a combination of a density-independent term b i (or d i ) and a density-dependent term μ i x i η i (or ν i x i η i ). 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 β 1 and β 2 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
1 x 1 d x 1 d t = b 1 μ 1 x 1 η 1 d 1 + ν 1 x 1 θ 1 + γ 1 x 2 δ 1 + x 2 1 x 2 d x 2 d t = b 2 μ 2 x 2 η 2 d 2 + ν 2 x 2 θ 2 + γ 2 x 1 δ 2 + x 1 .
In (2), γ i represents the maximal benefit received by species i from its mutualistic partner, while δ i 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 b i d i > 0 , i = 1 , 2 .
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
d x 1 d t = r 1 x 1 A 1 x 1 K 1 p + r 1 b 12 K 1 x 1 x 2 d x 2 d t = r 2 x 2 A 2 x 2 K 2 p + r 2 b 21 K 2 x 1 x 2 ,
being shown that if p = 1 , 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 b 12 b 21 < 1 , while if p > 1 , then the system (3) has a unique co-existing equilibrium regardless of the value of b 12 b 21 . Also, this co-existing equilibrium is globally asymptotically stable in ( 0 , ) × ( 0 , ) 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 p < 1 .
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
1 x 1 d x 1 d t = g 1 ( x 1 , x 2 ) d 1 ( x 1 ) , 1 x 2 d x 2 d t = g 2 ( x 1 , x 2 ) d 2 ( x 2 ) ,
in which g is a positive, continuously differentiable function on R 2 and d is a positive, continuously differentiable function on R under the following consistency assumptions:
(C1)
g 1 ( x 1 , x 2 ) d 1 ( x 1 ) is decreasing as a function of x 1 for each x 2 and g 2 ( x 1 , x 2 ) d 2 ( x 2 ) is decreasing as a function of x 2 for each x 1 ;
(C2)
g 1 ( s 1 x , s 2 x ) d 1 ( s 1 x ) and g 2 ( s 1 x , s 2 x ) d 2 ( s 2 x ) are ultimately decreasing for all s 1 , s 2 > 0 ;
as well as the following self-limiting growth assumption:
(L)
There are K 1 , K 2 > 0 such that
g 1 ( x 1 , 0 ) d 1 ( x 1 ) ( x 1 K 1 ) < 0 , x 1 [ 0 , ) , x 1 K 1 , g 2 ( 0 , x 2 ) d 2 ( x 2 ) ( x 2 K 2 ) < 0 , x 1 [ 0 , ) , x 2 K 2 ;
and the following mutualistic assumption:
(M)
g 1 ( x 1 , x 2 ) is increasing in x 2 for each x 1 and g 2 ( x 2 , x 1 ) is increasing in x 1 for each x 2 .
Note that for concrete models, the g i and d i values are not necessarily uniquely defined, but they are only defined up to functions of x i . For this very reason, g i values may not represent the exact forms of the per capita mutualistic contributions (although they are certainly related to it) and d i 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
r ¯ 1 = g 1 ( 0 , 0 ) d 1 ( 0 ) , r ¯ 2 = g 2 ( 0 , 0 ) d 2 ( 0 )
and observe that r ¯ 1 , r ¯ 2 > 0 . Subsequently, using also (M), it is seen that the semi-trivial equilibria ( K 1 , 0 ) and ( 0 , K 2 ) 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
R 1 ( s 1 , s 2 ) = lim x g 1 ( s 1 x , s 2 x ) d 1 ( s 1 x ) , R 2 ( s 1 , s 2 ) = lim x g 2 ( s 1 x , s 2 x ) d 2 ( s 2 x ) .
The parameters R 1 and R 2 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 R 1 and R 2 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 α 1 , α 2 > 0 such that R 1 ( α 1 , α 2 ) < 1 and R 2 ( α 1 , α 2 ) < 1 .
(D)
lim inf x d 1 ( x ) > 0 , lim inf x d 2 ( x ) > 0 .
While (B) is a boundedness assumption, (D) essentially tells us that d i values act as limiting terms. In fact, d 1 , d 2 are usually increasing and unbounded in concrete situations, which easily implies (D). Also, although d 1 and d 2 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. 
Let us define
h ( t ) = max x 1 ( t ) α 1 , x 2 ( t ) α 2 .
Suppose that for a given t, one has that h ( t ) = x 1 ( t ) α 1 . Then,
x 2 ( t ) α 2 x 1 ( t ) α 1 .
Consequently,
1 x 1 d x 1 d t d 1 ( x 1 ) g 1 ( x 1 , α 2 x 1 α 1 ) d 1 ( x 1 ) 1 = d 1 ( x 1 ) g 1 ( α 1 x 1 α 1 , α 2 x 1 α 1 ) d 1 ( α 1 x 1 α 1 ) 1
Since
lim x 1 α 1 g 1 ( α 1 x 1 α 1 , α 2 x 1 α 1 ) d 1 ( α 1 x 1 α 1 ) = R 1 ( α 1 , α 2 ) < 1 ,
it follows that 1 x 1 d x 1 d t < δ 1 for some given δ 1 > 0 if x 1 is large enough. If h ( t ) = x 2 ( t ) α 2 instead, then one may infer via a parallel argument that 1 x 2 d x 2 d t < δ 2 if for some given δ 2 > 0 if x 2 is large enough.
Trouble is, h may not necessarily be differentiable for all t > 0 (specifically, it may not be differentiable for values of t such that x 1 ( t ) α 1 = x 2 ( t ) α 2 ). However, one may infer from the above that D + h ( t ) < δ h ( t ) for some δ > 0 when h ( t ) is large enough, which implies that h ( t ) is ultimately bounded from above. Consequently, so are both x 1 and x 2 . □
Notice that (C2), in fact, is used in our abstract framework only to ensure that R 1 ( α 1 , α 2 ) , R 2 ( α 1 , α 2 ) are well-defined for all α 1 , α 2 > 0 . If the values of R 1 ( α 1 , α 2 ) , R 2 ( α 1 , α 2 ) 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 ( 0 , ) 2 are ultimately bounded from below by positive constants (i.e., (4) is uniformly persistent).
Proof. 
From (M) and (L), it follows that
d x 1 d t x 1 g 1 ( x 1 , 0 ) d 1 ( x 1 ) r ¯ 1 2 x 1 > 0
as long as x 1 stays in some interval ( 0 , δ ] , δ < K 1 , from which lim inf t x 1 ( t ) δ and the species x 1 is uniformly persistent. The uniform persistence of x 2 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 E * , by Theorem 2.8.6 of Bhatia and Szegö [17]. The uniqueness and stability of E * , 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 E * is unique. Then, E * is globally asymptotically stable.
Proof. 
The global stability of the co-existence equilibrium E * follows from Lemma 1 by ruling out periodic solutions with the help of Dulac criterion. The Dulac function corresponding to (4) is φ ( x 1 , x 2 ) : = 1 x 1 x 2 d 1 ( x 1 ) d 2 ( x 2 ) , which leads to
x 1 [ φ ( x 1 , x 2 ) d x 1 d t ] + x 2 [ φ ( x 1 , x 2 ) d x 2 d t ] = = 1 x 2 d 2 ( x 2 ) x 1 g 1 ( x 1 , x 2 ) d 1 ( x 1 ) + 1 x 1 d 1 ( x 1 ) x 2 g 2 ( x 2 , x 1 ) d 2 ( x 2 ) < 0 .
Consequently, E * 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 x , y R n , x = ( x 1 , x 2 , x n ) , y = ( y 1 , y 2 , y n ) , we shall write x y if x i y i for 1 i n , x < y if x y , x y and x y if x i < y i for 1 i n . Also, given h : R + n R n , we shall denote by D h ( x ) its associated Jacobian matrix computed at x and say that h is cooperative if D h ( x ) is a Metzler matrix (all its off-diagonal components are non-negative) for all x R + n .
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 system
x i = x i f i ( x 1 , x 2 , , x n ) , 1 i n ,
under the following assumptions
(H1) 
F is cooperative;
(H2) 
F ( 0 ) 0 ;
(H3) 
D F ( y ) D F ( z ) when z y 0 .
F being the vector-valued function defined by F ( x ) = ( f i ( x ) ) 1 i n . If (5) possesses an equilibrium point in int ( R + n ) , then this equilibrium point is unique in int ( R + n ) and also globally asymptotically stable in int ( R + n ) .
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 g 1 ( 0 , 0 ) d 1 ( 0 ) > 0 , g 2 ( 0 , 0 ) d 2 ( 0 ) > 0 , inequalities which, as previously mentioned, are an outcome of (L). Also,
D F ( x 1 , x 2 ) = g 1 x 1 ( x 1 , x 2 ) d 1 ( x 1 ) g 1 x 2 ( x 1 , x 2 ) g 2 x 1 ( x 1 , x 2 ) g 2 x 1 ( x 1 , x 2 ) d 2 ( x 2 ) .
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.

3. The Model with Linear Mutualistic Contribution

3.1. Stability Results

We start with an analysis of the model (1). To fix the ideas, let us suppose that max { η 1 , θ 1 } = θ 1 > 1 , max { η 2 , θ 2 } = θ 2 > 1 (that is, for each species, at least one of the birth/death terms has an accelerating density dependence), min { η i , θ i } > 0 , i = 1 , 2 . Rearranging the right-hand sides of (1) and denoting r i = b i d i , i = 1 , 2 , we obtain the following system
1 x 1 d x 1 d t = r 1 + β 1 x 2 μ 1 x 1 η 1 + ν 1 x 1 θ 1 1 x 2 d x 2 d t = r 2 + β 2 x 1 μ 2 x 2 η 2 + ν 2 x 2 θ 2 ,
which fits the abstract framework (4) with
g 1 ( x 1 , x 2 ) = r 1 + β 1 x 2 , g 2 ( x 1 , x 2 ) = r 2 + β 2 x 1 d 1 ( x 1 ) = μ 1 x 1 η 1 + ν 1 x 1 θ 1 d 2 ( x 2 ) = μ 2 x 2 η 2 + ν 2 x 2 θ 2 .
This leads to
g 1 ( x 1 , x 2 ) d 1 ( x 1 ) = r 1 + β 1 x 2 μ 1 x 1 η 1 + ν 1 x 1 θ 1 , g 2 ( x 1 , x 2 ) d 2 ( x 2 ) = r 2 + β 2 x 1 μ 2 x 2 η 2 + ν 2 x 2 θ 2 ,
conditions (M), (D) and (C1) being then obviously satisfied. Let us then proceed with verifying that the other hypotheses of Theorem 1 are satisfied. It is seen that
g 1 ( x 1 , 0 ) d 1 ( x 1 ) = r 1 μ 1 x 1 η 1 ν 1 x 1 θ 1 ,
r 1 > 0 and the right-hand side is strictly decreasing, so the first half of (L) holds, K 1 being the unique positive root of g 1 ( x 1 , 0 ) d 1 ( x 1 ) = 0 . The second half holds by a similar argument. Since
g 1 ( s 1 x , s 2 x ) d 1 ( s 1 x ) = r 1 + β 1 s 2 x μ 1 ( s 1 x ) η 1 + ν 1 ( s 1 x ) θ 1 ,
it follows that
d d x g 1 ( s 1 x , s 2 x ) d 1 ( s 1 x ) = r 1 μ 1 η 1 ( s 1 x ) η 1 1 + ν 1 θ 1 ( s 1 x ) θ 1 1 β 1 s 2 μ 1 ( η 1 1 ) ( s 1 x ) η 1 ( μ 1 ( s 1 x ) η 1 + ν 1 ( s 1 x ) θ 1 ) 2 β 1 s 2 ( θ 1 1 ) ( s 2 x ) θ 1 ( μ 1 ( s 1 x ) η 1 + ν 1 ( s 1 x ) θ 1 ) 2 .
If η 1 1 , then
d d x g 1 ( s 1 x , s 2 x ) d 1 ( s 1 x ) 0 for all x > 0 .
If η 1 ( 0 , 1 ) , since in this case η 1 < θ 1 , it follows that
lim x r 1 μ 1 η 1 ( s 1 x ) η 1 1 + ν 1 θ 1 ( s 1 x ) θ 1 1 β 1 s 2 μ 1 ( η 1 1 ) ( s 1 x ) η 1 β 1 s 2 ( θ 1 1 ) ( s 2 x ) θ 1 =
and consequently, g 1 ( s 1 x , s 2 x ) d 1 ( s 1 x ) is ultimately decreasing. By a similar argument, g 2 ( s 1 x , s 2 x ) d 2 ( s 2 x ) is also ultimately decreasing, and consequently, (C2) is satisfied. In these settings, since θ 1 , θ 2 > 1 ,
R 1 ( α 1 , α 2 ) = lim x g 1 ( α 1 x , α 2 x ) d 1 ( α 1 x ) = lim x r 1 + β 1 α 2 x μ 1 ( α 1 x ) η 1 + ν 1 ( α 1 x ) θ 1 = 0 , R 2 ( α 1 , α 2 ) = lim x g 2 ( α 1 x , α 2 x ) d 2 ( α 2 x ) = lim x r 2 + β 2 α 1 x μ 2 ( α 2 x ) η 2 + ν 2 ( α 2 x ) θ 2 = 0
for all α 1 , α 2 > 0 . Also,
g 1 ( x 1 , 0 ) d 1 ( x 1 ) = r 1 μ 1 x 1 η 1 + ν 1 x 1 θ 1 g 2 ( 0 , x 2 ) d 2 ( x 2 ) = r 2 μ 2 x 2 η 2 + ν 2 x 2 θ 2 ,
so (L) holds (note that the right-hand sides are strictly decreasing and positive at 0), K 1 and K 2 being the unique roots of
r 1 μ 1 x 1 η 1 + ν 1 x 1 θ 1 = 0 , r 2 μ 2 x 2 η 2 + ν 2 x 2 θ 2 = 0 ,
respectively. From the above considerations, we see that the hypotheses of Theorem 1 are satisfied and obtain the following result.
Theorem 3. 
If max { η i , θ i } > 1 , i 1 , 2 , then (1) has no periodic solution and all its solutions are ultimately uniformly bounded from both above and below. Also, there exists a co-existence equilibrium E * which is globally asymptotically stable provided that it is unique.
Note that Theorem 3 does not require both birth/death terms for each species to have an accelerating density dependence; only one of them suffices. Naturally, in the given context, one does not have enough information on the shape of the nullclines to deduce the uniqueness of the co-existence equilibrium from a condition involving the maximal values of η i and θ i only. For further information regarding this aspect, let us observe that in our context, with the notations of Theorem 2,
F ( x 1 , x 2 ) = r 1 + β 1 x 2 μ 1 x 1 η 1 + ν 1 x 1 θ 1 r 2 + β 2 x 1 μ 2 x 2 η 2 + ν 2 x 2 θ 2 , D F ( x 1 , x 2 ) = μ 1 η 1 x 1 η 1 1 + ν 1 θ 1 x 1 θ 1 1 β 2 β 1 μ 2 η 2 x 2 η 2 1 + ν 2 θ 2 x 2 θ 2 1 ,
which imply that (H1) and (H2) are satisfied. For (H3) to hold, one needs η i 1 0 , θ i 1 0 , i 1 , 2 , which leads to the following result thanks to Theorem 2 and Lemmas 1 and 2.
Theorem 4. 
If max { η i , θ i } > 1 , min { η i , θ i } 1 , i 1 , 2 , then there exists a unique co-existence equilibrium E * of (1), which is globally asymptotically stable.
Note that Theorem 4 has stricter requirements on the birth/death terms: for each species, one of them should be accelerating and the other one should be at least linear. Consequently, the closer the system (1) is to manifesting “total” accelerating density dependence as far as the intrinsic dynamics of each species is concerned, the nearer to having an unique, globally asymptotically stable co-existence equilibrium it is.
We now provide an illustrative numerical example for the outcomes of Theorem 4. For Figure 1, r 1 = r 2 = 2 , β 1 = β 2 = 4 , μ 1 = μ 2 = ν 1 = ν 2 = 0.5 , η 1 = η 2 = 1.5 , θ 1 = θ 2 = 2 . In this situation, R 1 ( α 1 , α 2 ) = R 2 ( α 1 , α 2 ) = 0 for all α 1 , α 2 > 0 and both birth and death functions have accelerating density dependence. Consequently, (1) is covered by Theorem 4 and there is a unique co-existence equilibrium of (1), which is globally asymptotically stable.
Both Theorems 3 and 4 rely on the assumption that max { η i , θ i } > 1 , i 1 , 2 , so none of them apply to Gause–Witt mutualisms. A quick inspection of our argument shows, however, that η i = θ i = 1 , i 1 , 2 (i.e., (1) represents a classical Gause–Witt mutualism), then (M), (D), (C1), (C2), (L) still hold true,
R 1 ( α 1 , α 2 ) = β 1 α 2 ( μ 1 + ν 1 ) α 1 , R 2 ( α 1 , α 2 ) = β 2 α 1 ( μ 2 + ν 2 ) α 2
and one can find suitable α 1 , α 2 so that (B) holds provided that β 1 β 2 ( μ 1 + ν 1 ) ( μ 2 + ν 2 ) < 1 . Consequently, the hypotheses of Theorem 1 are satisfied and one obtains the following classical result concerning the stability of Gause–Witt mutualisms.
Theorem 5. 
If η i = θ i = 1 , i 1 , 2 , then there exists a unique co-existence equilibrium E * of (1), which is globally asymptotically stable, provided that
β 1 β 2 < ( μ 1 + ν 1 ) ( μ 2 + ν 2 ) .
For Figure 2, r 1 = r 2 = 2 , β 1 = β 2 = 2 , μ 1 = μ 2 = ν 1 = ν 2 = 1.5 , η 1 = η 2 = θ 1 = θ 2 = 1 . In this situation, (1) becomes a Gause–Witt mutualism covered by Theorem 5 since β 1 β 2 = 4 , ( μ 1 + ν 1 ) ( μ 2 + ν 2 ) = 9 , so β 1 β 2 < ( μ 1 + ν 1 ) ( μ 2 + ν 2 ) holds, and consequently, there is a unique co-existence equilibrium of (1), which is globally asymptotically stable.
Note also that by an argument similar to the one displayed above, one observes that Theorem 4 still holds if only one of max { η i , θ i } is > 1 , on the condition that the other η i and θ i are both equal to 1. Also, Theorem 4 slightly improves Theorem 3.1 of [15], to which it reduces when η 1 = θ 1 and η 2 = θ 2 (equal exponents for birth and death processes).

3.2. Alternate Proofs for the Existence and Uniqueness of the Co-Existence Equilibrium

In the above Theorems 3 and 4, the considerations relating to the existence and uniqueness of the co-existence equilibrium are somewhat obscured, as they are essentially an outcome of qualitative properties of monotone dynamical systems. For the sake of completeness, let us give elementary proofs for the existence and uniqueness of the co-existence equilibrium, respectively, under the hypothesis that max { η i , θ i } > 1 , i 1 , 2 .

3.2.1. Existence

The equilibrium conditions can be written as follows
x 2 * = μ 1 x 1 * η 1 + ν 1 x 1 * θ 1 r 1 β 1 , x 1 * = μ 2 x 2 * η 2 + ν 2 x 2 * θ 2 r 2 β 2
Let us consider
a 1 : [ 0 , ) R , a 1 ( x ) = μ 1 x η 1 + ν 1 x θ 1 r 1 β 1 , a 2 : [ 0 , ) R , a 2 ( x ) = μ 2 x η 2 + ν 2 x θ 2 r 2 β 2 .
Since
a 1 ( 0 ) = r 1 β 1 < 0 , lim x a 1 ( x ) = +
and a 1 is continuous, there is a unique x 1 0 > 0 such that a 1 ( x 1 0 ) = 0 . Now, the component x 1 * is a solution of the equation b 1 ( x ) = 0 , where
b 1 : [ x 1 0 , ) R , b 1 ( x ) = μ 2 μ 1 x η 1 + ν 1 x θ 1 r 1 β 1 η 2 + ν 2 μ 1 x η 1 + ν 1 x θ 1 r 1 β 1 θ 2 β 2 x r 2 .
Since
b 1 ( x 1 0 ) = β 2 x 1 0 r 2 < 0 , lim x b 1 ( x ) = +
and b 1 is continuous, there is at least an x 1 * such that b 1 ( x 1 * ) = 0 . Consequently, there is a corresponding x 2 * satisfying (8) and (1) has at least a positive equilibrium.

3.2.2. Uniqueness

We try to solve (8) as a system in the x 1 O x 2 plane, of the form
x 2 = a 1 ( x 1 ) , x 1 = a 2 ( x 2 ) .
One notes first that a 1 , a 2 are strictly increasing. If both a 1 , a 2 are convex, then the co-existence equilibrium is unique. Since
a i ( x ) = μ i η i ( η i 1 ) x η 1 2 + ν i θ i ( θ i 1 ) x θ i 2 β i , i 1 , 2 ,
it follows that min { η i , θ i } 1 , i 1 , 2 , which ensures the uniqueness of the co-existence equilibrium.

3.3. Further Local Stability Considerations

Up to now, we have been concerned mostly with the global stability of the co-existence equilibria but, of course, this is not the only possible outcome. Let us now be concerned with local stability properties instead. As seen above, if one of the exponents η i , θ i , i 1 , 2 , is < 1 , the co-existence equilibrium might not be unique, so it makes sense to further investigate the local stability of equilibria. Let us observe that the Jacobian associated to (1) is
J = r 1 ( η 1 + 1 ) μ 1 x 1 η 1 ( θ 1 + 1 ) ν 1 x 1 θ 1 + β 1 x 2 β 1 x 1 β 2 x 2 r 2 ( η 2 + 1 ) μ 2 x 2 η 2 ( θ 2 + 1 ) ν 2 x 2 θ 2 + β 2 x 1
and let E * = ( x 1 * , x 2 * ) be a co-existence equilibrium of (1). As the equilibrium relations lead to
r 1 + β 1 x 2 * ( μ 1 ( x 1 * ) η 1 + ν 1 ( x 1 * ) θ 1 ) = 0 , r 2 + β 2 x 1 * ( μ 2 ( x 2 * ) η 2 + ν 2 ( x 2 * ) θ 2 ) = 0 ,
it follows that
J E * = η 1 μ 1 ( x 1 * ) η 1 + θ 1 ν 1 ( x 1 * ) θ 1 β 1 x 1 * β 2 x 2 * η 2 μ 2 ( x 2 * ) η 2 + θ 2 ν 2 ( x 2 * ) θ 2 .
Noting that the trace of J E * is negative, it is seen that for E * to be stable, one needs det ( J E * ) to be positive. It then follows that a co-existence equilibrium E * = ( x 1 * , x 2 * ) of (1) is stable provided that
η 1 μ 1 ( x 1 * ) η 1 + θ 1 ν 1 ( x 1 * ) θ 1 η 2 μ 2 ( x 2 * ) η 2 + θ 2 ν 2 ( x 2 * ) θ 2 > β 1 β 2 x 1 * x 2 *
If min { η i , θ i } 1 , i 1 , 2 , then
η 1 μ 1 ( x 1 * ) η 1 + θ 1 ν 1 ( x 1 * ) θ 1 η 2 μ 2 ( x 2 * ) η 2 + θ 2 ν 2 ( x 2 * ) θ 2 μ 1 ( x 1 * ) η 1 + ν 1 ( x 1 * ) θ 1 μ 2 ( x 2 * ) η 2 + ν 2 ( x 2 * ) θ 2 = ( r 1 + β 1 x 2 * ) ( r 2 + β 2 x 1 * ) > β 1 β 2 x 1 * x 2 * ,
so (9) is satisfied, and in this case, any co-existence equilibrium is locally stable whenever it exists, giving further credence to the already observed fact that accelerating density dependence promotes the stability of equilibria. Note that condition min { η i , θ i } 1 , i 1 , 2 (just one of the two stability conditions employed in Theorem 4) alone does not ensure the existence of the positive equilibrium (for the Gause–Witt mutualism, for instance, one still needs β 1 β 2 < ( μ 1 + ν 1 ) ( μ 2 + ν 2 ) for that to hold).

4. The Model with Saturating Mutualistic Contribution

In this case, due to the saturation of mutualistic benefits at higher population densities, one would expect less demanding conditions for the exponents η i and θ i , i 1 , 2 , in order to obtain the existence and stability of the co-existence equilibrium. The reason is that those exponents appear in the removal rates d 1 and d 2 only and there are now weaker unboundedness tendencies coming from the mutualistic contributions to be balanced.

4.1. Stability Results

Let us suppose again that min { η i , θ i } > 0 , i = 1 , 2 . Rearranging the right-hand sides of (2) and denoting r i = b i d i , i = 1 , 2 , we obtain the following system
1 x 1 d x 1 d t = r 1 + γ 1 x 2 δ 1 + x 2 μ 1 x 1 η 1 + ν 1 x 1 θ 1 1 x 2 d x 2 d t = r 2 + γ 2 x 1 δ 2 + x 1 μ 2 x 2 η 2 + ν 2 x 2 θ 2 ,
which fits the abstract framework (4) with
g 1 ( x 1 , x 2 ) = r 1 + γ 1 x 2 δ 1 + x 2 , g 2 ( x 1 , x 2 ) = r 2 + γ 2 x 1 δ 2 + x 1 d 1 ( x 1 ) = μ 1 x 1 η 1 + ν 1 x 1 θ 1 d 2 ( x 2 ) = μ 2 x 2 η 2 + ν 2 x 2 θ 2 .
This leads to
g 1 ( x 1 , x 2 ) d 1 ( x 1 ) = r 1 + γ 1 x 2 δ 1 + x 2 μ 1 x 1 η 1 + ν 1 x 1 θ 1 , g 2 ( x 1 , x 2 ) d 2 ( x 2 ) = r 2 + γ 2 x 1 δ 2 + x 1 μ 2 x 2 η 2 + ν 2 x 2 θ 2 ,
conditions (M), (D) and (C1) being again obviously satisfied. Assumption (L) holds by exactly the same argument displayed in Section 3. Since
g 1 ( s 1 x , s 2 x ) d 1 ( s 1 x ) = r 1 + γ 1 s 2 x δ 1 + s 2 x μ 1 ( s 1 x ) η 1 + ν 1 ( s 1 x ) θ 1 ,
it follows that
d d x g 1 ( s 1 x , s 2 x ) d 1 ( s 1 x ) = r 1 μ 1 η 1 ( s 1 x ) η 1 1 + ν 1 θ 1 ( s 1 x ) θ 1 1 γ 1 s 2 2 x ( δ 1 + s 2 x ) 2 ( μ 1 ( s 1 x ) η 1 + ν 1 ( s 1 x ) θ 1 ) ( μ 1 ( s 1 x ) η 1 + ν 1 ( s 1 x ) θ 1 ) 2 .
Consequently,
d d x g 1 ( s 1 x , s 2 x ) d 1 ( s 1 x ) 0 for all x > 0
and (C2) is satisfied. In these settings,
R 1 ( α 1 , α 2 ) = lim x g 1 ( α 1 x , α 2 x ) d 1 ( α 1 x ) = lim x r 1 + γ 1 α 2 x δ 1 + α 2 x μ 1 ( α 1 x ) η 1 + ν 1 ( α 1 x ) θ 1 = 0 , R 2 ( α 1 , α 2 ) = lim x g 2 ( α 1 x , α 2 x ) d 2 ( α 2 x ) = lim x r 2 + γ 2 α 1 x δ 2 + α 1 x μ 2 ( α 2 x ) η 2 + ν 2 ( α 2 x ) θ 2 = 0
for all α 1 , α 2 > 0 . Also, (L) holds by exactly the same argument as above, which implies that the hypotheses of Theorem 1 are satisfied, and consequently, the following result holds.
Theorem 6. 
If min { η i , θ i } > 0 , i 1 , 2 , then (2) has no periodic solution and all its solutions are ultimately uniformly bounded from both above and below. Also, there exists a co-existence equilibrium E * which is globally asymptotically stable provided that it is unique.
Notice that for (2), the boundedness of solutions is not tied to the birth/death terms being accelerating for any species; just min { η i , θ i } > 0 , i 1 , 2 suffices. Not unexpectedly, saturating mutualistic benefits lead to a completely different picture: one in which no species grows unbounded under any circumstances not even when the per capita growth rates are decelerating.
Furthermore, with the notations of Theorem 2,
F ( x 1 , x 2 ) = r 1 + γ 1 x 2 δ 1 + x 2 μ 1 x 1 η 1 + ν 1 x 1 θ 1 r 2 + γ 2 x 1 δ 2 + x 1 μ 2 x 2 η 2 + ν 2 x 2 θ 2 , D F ( x 1 , x 2 ) = μ 1 η 1 x 1 η 1 1 + ν 1 θ 1 x 1 θ 1 1 γ 1 δ 1 ( δ 1 + x 2 ) 2 γ 2 δ 2 ( δ 2 + x 1 ) 2 μ 2 η 2 x 2 η 2 1 + ν 2 θ 2 x 2 θ 2 1 .
Again, it is seen that (H1) and (H2) are satisfied and for (H3) to hold, one needs η i 1 0 , θ i 1 0 , i 1 , 2 , which leads to the following stability result thanks to Theorem 2.
Theorem 7. 
If min { η i , θ i } 1 , i 1 , 2 , then there exists a unique co-existence equilibrium E * of (2), which is globally asymptotically stable.
From the above result, one sees that in spite of what happened with the boundedness of the solutions, which necessitated significantly less demanding hypotheses than in the case of linear mutualistic contributions, the stability picture still did not change, as only accelerating per capita growth rates ensures the global stability of the co-existence equilibrium.
For Figure 3, r 1 = r 2 = 2 , γ 1 = γ 2 = 4 , δ 1 = δ 2 = 10 , μ 1 = μ 2 = 0.8 , ν 1 = ν 2 = 0.5 , η 1 = η 2 = 2 , and θ 1 = θ 2 = 1.5 . In this situation, R 1 ( α 1 , α 2 ) = R 2 ( α 1 , α 2 ) = 0 for all α 1 , α 2 > 0 and both birth and death functions have accelerating density dependence. Consequently, (2) is covered by Theorem 7 and there is a unique co-existence equilibrium of (2), which is globally asymptotically stable.

4.2. Further Local Stability Considerations

As in Section 3.3 for the model with linear mutualistic contributions, we shall now be concerned with the situation in which a co-existence equilibrium of (2) is not necessarily globally asymptotically stable. Let us observe that the Jacobian associated to (2) is
J = r 1 ( η 1 + 1 ) μ 1 x 1 η 1 ( θ 1 + 1 ) ν 1 x 1 θ 1 + γ 1 x 2 δ 1 + x 2 γ 1 δ 1 x 1 ( δ 1 + x 2 ) 2 γ 2 δ 2 x 2 ( δ 2 + x 1 ) 2 r 2 ( η 2 + 1 ) μ 2 x 2 η 2 ( θ 2 + 1 ) ν 2 x 2 θ 2 + γ 2 x 1 δ 2 + x 1
and let E * = ( x 1 * , x 2 * ) be a co-existence equilibrium of (1). As the equilibrium relations lead to
r 1 ( μ 1 ( x 1 * ) η 1 + ν 1 ( x 1 * ) θ 1 ) + γ 1 x 2 * δ 1 + x 2 * = 0 , r 2 ( μ 2 ( x 2 * ) η 2 + ν 2 ( x 2 * ) θ 2 ) + γ 2 x 1 * δ 2 + x 1 * = 0 ,
it follows that
J E * = η 1 μ 1 ( x 1 * ) η 1 + θ 1 ν 1 ( x 1 * ) θ 1 γ 1 δ 1 x 1 * ( δ 1 + x 2 * ) 2 γ 2 δ 2 x 2 * ( δ 2 + x 1 * ) 2 η 2 μ 2 ( x 2 * ) η 2 + θ 2 ν 2 ( x 2 * ) θ 2 .
Noting that the trace of J E * is negative, it is seen that for E * to be stable, one needs det ( J E * ) to be positive. It then follows that a co-existence equilibrium E * = ( x 1 * , x 2 * ) of (1) is stable provided that
η 1 μ 1 ( x 1 * ) η 1 + θ 1 ν 1 ( x 1 * ) θ 1 η 2 μ 2 ( x 2 * ) η 2 + θ 2 ν 2 ( x 2 * ) θ 2 > γ 1 γ 2 δ 1 δ 2 x 1 * x 2 * ( δ 1 + x 2 * ) 2 ( δ 2 + x 1 * ) 2 .
One might be puzzled by the seemingly complicated expression of (11), since Theorem 7 proposes a much simpler condition for a much better outcome (global stability, rather than local). However, note that if min { η i , θ i } 1 , i 1 , 2 , then
η 1 μ 1 ( x 1 * ) η 1 + θ 1 ν 1 ( x 1 * ) θ 1 η 2 μ 2 ( x 2 * ) η 2 + θ 2 ν 2 ( x 2 * ) θ 2 μ 1 ( x 1 * ) η 1 + ν 1 ( x 1 * ) θ 1 μ 2 ( x 2 * ) η 2 + ν 2 ( x 2 * ) θ 2 = r 1 + γ 1 x 2 * δ 1 + x 2 * r 2 + γ 2 x 1 * δ 2 + x 1 * > γ 1 γ 2 x 1 * x 2 * ( δ 1 + x 2 * ) ( δ 2 + x 1 * ) > γ 1 γ 2 x 1 * x 2 * ( δ 1 + x 2 * ) ( δ 2 + x 1 * ) δ 1 δ 1 + x 2 * δ 2 δ 2 + x 1 * > γ 1 γ 2 δ 1 δ 2 x 1 * x 2 * ( δ 1 + x 2 * ) 2 ( δ 2 + x 1 * ) 2 ,
so (11) does hold in this case.

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 n × n square matrix A, the directed graph G ( A ) associated to it is a directed graph with vertices 1 , 2 , , n such that there exists an arc ( j , k ) leading from j to k, 1 j , k n , if and only if a j k 0 . A will then be called irreducible if G ( A ) 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 x = h ( x ) , where h : R + n R n is a continuously differentiable map. Assume that
(D1) 
h is cooperative on R + n and D h ( x ) is irreducible for any x R + n ;
(D2) 
h ( 0 ) = 0 and h i ( x ) 0 for all x R + n with x i = 0 , i = 1 , 2 , , n ;
(D3) 
s ( D h ( 0 ) ) > 0 .
Let φ ( t , y 0 ) be the maximally extended solution of y = h ( y ) initiating at y 0 . Then, either
1. 
For any y R + n 0 , lim t | φ ( t , y ) | = + , or alternatively
2. 
There exists y * 0 with h ( y * ) = 0 such that for any y with 0 < y y * , lim t φ ( t , y ) = y * . Moreover, for any y > 0 , lim inf t φ ( t , y ) y * .
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 ( 0 , ) 2 , the solutions will not reach the boundaries of the first quadrant, so irreducibility is needed only on ( 0 , ) 2 . Let us then assume the stronger version of (M) below.
(SM)
g 1 x 2 ( x 1 , x 2 ) > 0 , g 2 x 1 ( x 1 , x 2 ) > 0 on ( 0 , ) 2 .
Then, the required irreducibility property does hold, since the off-diagonal elements of the Jacobian are strictly positive. Also, (D2) is obvious and
D h ( 0 ) = r ¯ 1 0 0 r ¯ 2
so (D3) also holds. One then obtains the following result.
Theorem 9. 
Assume that (SM) and (L) hold. Let ψ ( t , y 0 ) be the maximally extended solution of (4) initiating at y 0 . Then either
1. 
For any y R + n 0 , lim t | ψ ( t , y ) | = + , or alternatively
2. 
There exists a co-existence equilibrium y * of (4) such that for any y with 0 < y y * , lim t ψ ( t , y ) = y * . Moreover, for any y > 0 , lim inf t ψ ( t , y ) y * .
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, r 1 = r 2 = 0.02 , β 1 = β 2 = 0.02 , μ 1 = μ 2 = 0.008 , ν 1 = ν 2 = 0.005 .
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, r 1 = r 2 = 2 , β 1 = 2.4 , β 2 = 4 , η 1 = η 2 = 0.5 , θ 1 = θ 2 = 1.1 , μ 1 = 25 , μ 2 = 0.5 , ν 1 = ν 2 = 0.5 . Notice that among the two co-existence equilibria y 1 * = ( 0.13 , 2.99 ) and y 2 * = ( 3.04 , 18.05 ) , the “smaller” one ( y 1 * ) is the stable one, and the “larger” one ( y 2 * ) 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 η i , θ i > 0 . Also, as previously mentioned, (C2) becomes redundant since one can compute R 1 ( α 1 , α 2 ) and R 2 ( α 1 , α 2 ) 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), R 1 ( α 1 , α 2 ) = R 2 ( α 1 , α 2 ) = 0 for all α 1 , α 2 > 0 whenever min { η i , θ i } > 0 , i = 1 , 2 , 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 min { η i , θ i } > 0 , i = 1 , 2 .
For the model with linear mutualistic contribution (1), the picture becomes more complicated. Note that R 1 ( α 1 , α 2 ) = R 2 ( α 1 , α 2 ) = 0 for all α 1 , α 2 > 0 if max { η i , θ i } > 1 , i = 1 , 2 , so only one of η 1 , θ 1 and η 2 , θ 2 , respectively, needs to be > 1 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 R 1 ( α 1 , α 2 ) and R 2 ( α 1 , α 2 ) to be both less than 1, not 0. Note also that R i ( α 1 , α 2 ) = whenever min { η i , θ i } < 1 , i = 1 , 2 , 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.

Author Contributions

Conceptualization, P.G.; Software, P.G. and H.Z.; Validation, H.Z.; Formal analysis, P.G.; Writing—original draft, P.G. and H.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Data is contained within the article.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Holland, J.N.; Bronstein, J.L. Mutualism. In Population Dynamics, Volume 3 of Encyclopedia of Ecology; Jorgensen, S.E., Fath, B.D., Eds.; Elsevier: Amsterdam, The Netherlands, 2008; pp. 2485–2491. [Google Scholar]
  2. Gause, G.; Witt, A. Behavior of mixed populations and the problem of natural selection. Am. Nat. 1935, 65, 596–609. [Google Scholar] [CrossRef]
  3. Moore, C.M.; Catella, S.A.; Abbott, K.C. Population dynamics of mutualism and intraspecific density dependence: How θ-logistic-like density dependence affects mutualistic positive feedback. Ecol. Model. 2017, 368, 191–197. [Google Scholar] [CrossRef]
  4. Sibly, R.M.; Barker, D.; Denham, M.C.; Hone, J.; Pagel, M. On the regulation of populations of mammals, birds, fish, and insects. Science 2005, 309, 607–610. [Google Scholar] [CrossRef] [PubMed]
  5. Gilpin, M.E.; Ayala, F.J. Global models of growth and competition. Proc. Natl. Acad. Sci. USA 1973, 70, 3590–3593. [Google Scholar] [PubMed]
  6. Smith, F.E. Population dynamics in Daphnia magna and a new model for population growth. Ecology 1963, 44, 651–663. [Google Scholar]
  7. Smith, D.W.; Cooper, S.D. Competition among cladocera. Ecology 1982, 63, 1004–1015. [Google Scholar] [CrossRef]
  8. Coulson, T.; Ezard, T.H.G.; Pelletier, F.; Tavecchia, G.; Stenseth, N.C.; Childs, D.Z.; Pilkington, J.G.; Pemberton, J.M.; Kruuk, L.E.B.; Clutton-Brock, T.H.; et al. Estimating the functional form for the density dependence from life history data. Ecology 2008, 89, 1661–1674. [Google Scholar] [CrossRef] [PubMed]
  9. Ribeiro, F.; Cabella, B.C.T.; Martinez, A.S. Richards-like two species population dynamics model. Theory Biosci. 2014, 133, 135–143. [Google Scholar] [CrossRef] [PubMed]
  10. Wang, D. Dynamic behaviors of an obligate Gilpin–Ayala system. Adv. Diff. Equ. 2016, 2016, 270. [Google Scholar] [CrossRef]
  11. García-Algarra, J.; Galeano, J.; Pastor, J.M.; Iriondo, J.M.; Ramasco, J. Rethinking the logistic approach for population dynamics of mutualistic interactions. J. Theor. Biol. 2014, 363, 332–343. [Google Scholar] [CrossRef] [PubMed]
  12. Georgescu, P.; Maxin, D.; Zhang, H. Threshold boundedness conditions for n-species mutualisms. Nonlinearity 2017, 30, 4410–4427. [Google Scholar] [CrossRef]
  13. Maxin, D.; Georgescu, P.; Sega, L.; Berec, L. Global stability of the coexistence equilibrium for a general class of models of facultative mutualism. J. Biol. Dyn. 2017, 11, 339–364. [Google Scholar] [CrossRef] [PubMed]
  14. Georgescu, P.; Zhang, H. Lyapunov functionals for two-species mutualisms. Appl. Math. Comput. 2014, 226, 754–764. [Google Scholar] [CrossRef]
  15. Georgescu, P.; Zhang, H.; Maxin, D. The global stability of coexisting equilibria for three models of mutualism. Math. Biosci. Eng. 2016, 13, 101–118. [Google Scholar] [CrossRef] [PubMed]
  16. May, R.M. Models of two interacting populations. In Theoretical Ecology: Principles and Applications; May, R.M., Ed.; Blackwell Scientific Publications: London, UK, 1976; pp. 78–104. [Google Scholar]
  17. Bhatia, N.P.; Szegö, G.P. Dynamical Systems: Stability Theory and Applications. In Lecture Notes in Mathematics; Springer: Berlin/Heidelberg, Germany, 1967; Volume 35. [Google Scholar]
  18. Smith, H. On the asymptotic behavior of a class of deterministic models of cooperative species. SIAM J. Appl. Math. 1986, 46, 368–375. [Google Scholar] [CrossRef]
  19. Zhao, X.-Q.; Jing, Z.-J. Global asymptotic behavior in some cooperative systems of functional differential equations. Canad. Appl. Math. Quart. 1996, 4, 421–444. [Google Scholar]
  20. Gómez, J.M.; Verdú, M. Mutualism with plants drivesprimate diversification. Syst. Biol. 2012, 61, 567–577. [Google Scholar] [CrossRef] [PubMed]
  21. Chomicki, G.; Weber, M.; Antonelli, A.; Bascompte, J.; Kiers, T. The impact of mutualisms on species richness. Trends Ecol. Evol. 2019, 34, 698–711. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Phase portrait of model (1) with accelerating density dependences. Several solutions corresponding to distinct initial data are set in black, the isoclines are set in blue and the red dot is the subsequent co-existence equilibrium.
Figure 1. Phase portrait of model (1) with accelerating density dependences. Several solutions corresponding to distinct initial data are set in black, the isoclines are set in blue and the red dot is the subsequent co-existence equilibrium.
Mathematics 11 04373 g001
Figure 2. Phase portrait of model (1) with linear density dependences, the stability condition (7) being satisfied. Several solutions corresponding to distinct initial data are set in black, the isoclines are set in blue and the red dot is the subsequent co-existence equilibrium.
Figure 2. Phase portrait of model (1) with linear density dependences, the stability condition (7) being satisfied. Several solutions corresponding to distinct initial data are set in black, the isoclines are set in blue and the red dot is the subsequent co-existence equilibrium.
Mathematics 11 04373 g002
Figure 3. Phase portrait of model (2) with accelerating density dependences. Several solutions corresponding to distinct initial data are set in black, the isoclines are set in blue and the red dot is the subsequent co-existence equilibrium.
Figure 3. Phase portrait of model (2) with accelerating density dependences. Several solutions corresponding to distinct initial data are set in black, the isoclines are set in blue and the red dot is the subsequent co-existence equilibrium.
Mathematics 11 04373 g003
Figure 4. Phase portrait of model (1) with linear density dependences, the reverse of the stability condition (7) being satisfied. Several solutions corresponding to distinct initial data are set in black.
Figure 4. Phase portrait of model (1) with linear density dependences, the reverse of the stability condition (7) being satisfied. Several solutions corresponding to distinct initial data are set in black.
Mathematics 11 04373 g004
Figure 5. Phase portrait of model (1) with mixing density dependences. Several solutions corresponding to distinct initial data are represented by differently colored curves starting from small empty circles. The carmin red dots are the co-existence equilibria, obtained as intersections of the two isoclines.
Figure 5. Phase portrait of model (1) with mixing density dependences. Several solutions corresponding to distinct initial data are represented by differently colored curves starting from small empty circles. The carmin red dots are the co-existence equilibria, obtained as intersections of the two isoclines.
Mathematics 11 04373 g005
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Georgescu, P.; Zhang, H. Facultative Mutualisms and θ-Logistic Growth: How Larger Exponents Promote Global Stability of Co-Existence Equilibria. Mathematics 2023, 11, 4373. https://doi.org/10.3390/math11204373

AMA Style

Georgescu P, Zhang H. Facultative Mutualisms and θ-Logistic Growth: How Larger Exponents Promote Global Stability of Co-Existence Equilibria. Mathematics. 2023; 11(20):4373. https://doi.org/10.3390/math11204373

Chicago/Turabian Style

Georgescu, Paul, and Hong Zhang. 2023. "Facultative Mutualisms and θ-Logistic Growth: How Larger Exponents Promote Global Stability of Co-Existence Equilibria" Mathematics 11, no. 20: 4373. https://doi.org/10.3390/math11204373

APA Style

Georgescu, P., & Zhang, H. (2023). Facultative Mutualisms and θ-Logistic Growth: How Larger Exponents Promote Global Stability of Co-Existence Equilibria. Mathematics, 11(20), 4373. https://doi.org/10.3390/math11204373

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop