Global stability of the coexistence equilibrium for a general class of models of facultative mutualism

ABSTRACT Many models of mutualism have been proposed and studied individually. In this paper, we develop a general class of models of facultative mutualism that covers many of such published models. Using mild assumptions on the growth and self-limiting functions, we establish necessary and sufficient conditions on the boundedness of model solutions and prove the global stability of a unique coexistence equilibrium whenever it exists. These results allow for a greater flexibility in the way each mutualist species can be modelled and avoid the need to analyse any single model of mutualism in isolation. Our generalization also allows each of the mutualists to be subject to a weak Allee effect. Moreover, we find that if one of the interacting species is subject to a strong Allee effect, then the mutualism can overcome it and cause a unique coexistence equilibrium to be globally stable.


Introduction
A mutualism is an interaction between a pair of species that is beneficial for both parties and that is commonly manifested as an increase in an individual's capacity to grow, survive or reproduce [16]. Mutualisms are many and varied [18], including associations between pollinators and plants [31], seed dispersers and plants [26], mycorrhizae and plants [23], and cleaner species and their hosts [9]. By degree of dependency, a mutualist can be classified as facultative when it benefits from the interaction yet can survive on its own, or obligate when it requires an interaction as it is not able to survive without the other species. Moreover, mutualisms can be protective when a species gets physical protection by the other species (such as when clown fish hide in sea anemones) or resource-based when a species receives a resource (such as when bees collect nectar). Most likely, it is because of this diversity and complexity of mutualistic interactions that work on predatory and competitive interactions has always dominated in ecology.
A number of mathematical models have been proposed to describe mutualistic interactions, starting with the familiar Lotka-Volterra model [31]. It is a legacy of this pioneering mutualism model that a majority of its successors focus on facultative mutualisms, are composed of two structurally identical equations and aim to prevent populations from increasing without bounds as a result of reinforcing effects of their mutual benefits [17,18,31]. Once a stable coexistence equilibrium is shown to exist, a mathematical question of biological interest that typically follows is establishing its uniqueness and/or global stability. As one may imagine, looking for uniqueness and/or global stability of a coexistence equilibrium in a mutualism model (but also in any other similarly complex model) is a difficult task.
The high diversity of mutualistic interactions calls for a similar diversity of mathematical models [17,18,31]. For instance, the Lotka-Volterra mutualism model represents the situation in which the per capita effects of a (facultative) mutualism are density independent On the contrary, Wolin and Lawlor [32] introduced the following model to accommodate the situation in which a (facultative) mutualism has a higher per capita impact at higher densities In both models (1) and (2), r i is the intrinsic growth rate and K i is the carrying capacity of species, i = 1,2 when alone. In addition, b 12 and b 21 are positive constants that scale the beneficial effect of one species on the other. Hence, if one species is missing, the other grows in a logistic fashion. However, the mutualistic association is manifested differently in models (1) and (2), which translates into different signs of the interaction terms. In particular, while model (1) is a special case of a general model for some positive functions a 1 , a 2 , f 1 , f 2 , model (2) can rather be generalized as ], It would certainly be a good move if the issues of uniqueness and/or global stability of a coexistence equilibrium were resolved for some general classes of models of mutualism such as these.
In this direction, Georgescu and Zhang [12] formulated sufficient conditions for the global stability of a unique coexistence equilibrium for a general model of mutualism, in the form under various trade-offs between monotonicity and sign conditions of the (not necessarily positive) functions a 1 , a 2 , f 1 , f 2 , g 1 , g 2 and their combinations. Note that model (5) also includes models (1) and (2) as its special cases. Lyapunov functionals extending those of [30] were used to explore global stability in this case. Georgescu et al. [13] further enlarged applicability of the abstract functionals introduced in [12] to a model of mutualism with restricted growth rates proposed in [14] and to the versions of models (1) and (2) with the logistic growth rates replaced by the Richards ones. Despite a variety of proposed models of mutualism, there is still potential for developing others. As we have already emphasized above, most models of mutualism assume a symmetry of the interaction effect. That is, the mutualistic effect is commonly modelled using structurally identical equations [17,18,31], as is also exemplified by models (1) and (2). Although this may suffice for certain systems, in other cases the mutualists have different life histories, such as when representing a plant and an animal [31]. Not only may interaction terms differ structurally for each mutualist, but also beneficial effects of mutualism may impact different fitness components in each species. For instance, in the sea anemone-anemone fish mutualism, waste ammonia from the fish feed the symbiotic algae that are found in the anemone's tentacles [22]. Similarly, by consuming parasitic bugs on the zebras' or rhinos' skin, the oxpeckers acquire food and the zebras or rhinos get rid of pests [4]. Last but not least, the spider crabs prefer shallow waters and seek algae that would grow on their back. This way the crabs are camouflaged and thus more likely to escape predation, while the algae obtain an extra resource which is another place to live at [4]. Of course, the same fitness component may be affected in both mutualists, such as in the case of brood-site pollination mutualism [25] or when some ants nest inside thorns of trees of the genus Acacia; in exchange for shelter, the ants protect acacias from attack by herbivores [15]. Therefore, any sufficiently general class of mutualism models should also allow for equations that are structurally different and cover cases in which different fitness components are affected in different species.
Apart from some exceptions such as [21], all published models of mutualism assume that individual species grow in a logistic manner when alone. In the recent decades, however, biologists have started to acknowledge that populations may be subject to an Allee effect, a concept originally studied by Allee [1,2]. Allee effects occur when mean individual fitness declines with decreasing population density when the population is rare and can occur due to a plethora of causes, the most prominent being enhanced difficulties of individuals to find mates, escape predation or avoid inbreeding when population densities become low [6,7,11,29].
Allee effects can be classified as weak or strong. In both cases, the per capita population growth rate declines with decreasing population density when the population is rare. Weak Allee effects occur when, despite this decline, the per capita growth rate never becomes negative in low-density populations [6,29]. As with the logistic growth, populations always attain the carrying capacity, yet at a slower rate when rare. On the other hand, strong Allee effects occur when the per capita growth rate becomes negative once the population density falls below a critical value termed an Allee threshold [6,29]. Populations above the Allee threshold attain the carrying capacity, while those below the Allee threshold go extinct. For instance, the population model exhibits a strong Allee effect with the Allee threshold A; populations with densities below A go extinct, while those above A attain the carrying capacity K. An example of a model with a weak Allee effect is provided by the following modification of the classical logistic model: For p > 1, the per capita growth rate f (x) = rx p−1 (1 − x/K) is indeed positive and increasing at small densities. Since there is no Allee threshold in this model, the Allee effect is weak. If p = 1, the model (7) represents the classical logistic equation. Other single-species models of Allee effects have been proposed in [6,8,24]. An important role of mutualistic models such as the Lotka-Volterra one and many of its successors is to propose generalizations that hold for many specific mutualistic interactions. Ability to establish uniqueness and/or global stability of a coexistence equilibrium for a general class of mutualism models is certainly a decent step in this direction. The aim of this article is thus twofold. First, we propose a general class of models of facultative mutualism that covers a broad set of existing mutualistic models, including Equations (1) and (2), and allows for asymmetry of the mutual species effect (i.e. mutualistic effects need not be modelled in the same way for each species) and/or populations subject to some Allee effects. Second, we provide mathematical proofs for this class of models concerning the boundedness of their solutions and the global stability of the unique coexistence equilibrium whenever it exists.
The paper is organized as follows. In Section 2, we propose our general class of facultative models of mutualism and prove our main theoretical results. In Section 3, we provide several biologically motivated examples to illustrate how our general results may apply to specific models of mutualism. A discussion on the biological implications of our results as well as on future research avenues is given in Section 4.

The general model
Before developing and analysing the models of mutualism, we first introduce a singlespecies model which describes the dynamics of a species in the absence of mutualism. We assume that the dynamics of each species in isolation is modelled as where a and f are positive real-valued continuous functions, which shall be elaborated upon later on. In what follows, a is to be interpreted as a growth function, and f as a selflimiting function. We now introduce the following growth assumption: This ensures that the population approaches the carrying capacity K from every positive initial density, and therefore covers both logistic-like population dynamics and dynamics of a population subject to a weak Allee effect. This assumption will be modified later on when treating strong Allee effects in mutualistic interactions.

Mutualisms that affect the self-limiting function
Let us now consider the following abstract model of a mutualistic interaction: in which r 1 , r 2 > 0, a 1 , a 2 are positive functions continuously differentiable on R and f 1 , f 2 are positive functions continuously differentiable on R 2 . We shall also assume that the expressions a 1 (x 1 ) − f 1 (x 1 , 0) and a 2 (x 2 ) − f 2 (x 2 , 0) satisfy the growth assumption (G) with K 1 and K 2 as the respective single-species carrying capacities. We now introduce the following consistency assumptions: is strictly decreasing in x 1 and a 2 (x 2 )/f 2 (x 2 , x 1 ) is strictly decreasing in x 2 . (C3) a 1 (x 1 )/f 1 (x 1 , αx 1 ) and a 2 (x 2 )/f 2 (x 2 , αx 2 ) are eventually decreasing for all α > 0.
That is, where M 1 and M 2 are positive real numbers that may depend on α.
The common meaning of the assumptions (C1) and (C2) is that if one species is kept constant, then the other should still follow a logistic-like or weak Allee effect growth in which the ratio between the growth and self-limiting functions decreases as the population density increases. The meaning of (C3) is somewhat similar, supposing this time that the ratio between the population densities is kept constant. In particular, the assumption (C3) has intentionally a weaker monotonicity property which allows our boundedness theorems formulated and proved below to hold also for mutualism models including strong Allee effects.
We now define the mutualism-specific assumption: Noting that f 1 and f 2 appear with negative signs in the model (9), the assumption (M) characterizes the fact that the model (9) indeed represents a mutualistic interaction, since an increase in the density of a population benefits the other species, for all population densities.
Due to the assumption (C3), for all α > 0 there exist finite non-negative numbers R 1 (α), R 2 (α) given by Due to the assumption (M), R 1 , R 2 are increasing as functions of α, although they may not necessarily be strictly increasing (it shall be observed later on that one or both of them may even be identically 0). Also, by the assumption (G), R 1 (0) < 1 and R 2 (0) < 1.
The existence of a coexistence equilibrium E * of the model (9) can be characterized in terms of a joint quantitative property of R 1 and R 2 , as follows.

Lemma 2.1:
If the assumption (C3) holds for all positive x 1 and x 2 and there exists a coexistence equilibrium E * of the model (9), then there is α > 0 such that R 1 (α) < 1 and be a coexistence equilibrium of the model (9). Because of that = 1, and the conclusion follows using the assumption (C3).
This result suggests a characterization of whether solutions of the model (9) are bounded through conditions imposed on R 1 and R 2 . We first establish the boundedness of solutions of the model (9). Proof: First, we shall prove that any solution that starts with positive initial conditions will not have a limit point on any of the semi-axes. Due to (M) and (G), it follows that

Theorem 2.2: Let the assumptions
provided that x 1 , x 2 stay positive and under their corresponding carrying capacities. Consequently, solutions starting with positive initial conditions cannot get too close to any of the positive semi-axes and the system (9) is uniformly persistent.
To show that the solutions are bounded from above, we consider the function For a given t > 0, we have two possibilities.
. Then x 2 (t) ≤ αx 1 (t) and, due to (M) Since the expression inside the square brackets has limit R 1 (α) − 1 < 0, it follows that Again, since the expression inside the square brackets has limit Even though αx 1 and x 2 are differentiable, the function h may not necessarily be so for all t > 0 (it is, however, continuous for all t ≥ 0). Specifically, h may not be differentiable for values of t such that hx 1 (t) = x 2 (t). However, the following result [20, Lemma 1.2.1] provides a condition for the monotonicity of continuous functions using only Dini derivatives. From the above discussion, it is seen that Consequently, h is uniformly bounded from above. Since h(t) := max{αx 1 (t), x 2 (t)}, it follows that x 1 and x 2 are also uniformly bounded from above. Moreover, all solutions starting in the positive quadrant eventually enter the compact region of this quadrant defined by max{αx 1 , x 2 } ≤ max{αM 1 , M 2 } regardless of initial population densities. Finally, by Theorem 2.8.6 of Bhatia and Szegö [3] (see also [28], Theorem D.3), there is a coexistence equilibrium E * .
A well-known sufficient condition to rule out the existence of periodic orbits of a dynamical system in a simply connected region of the plane is given in the following theorem.

Theorem 2.4 (Dulac criterion):
Let D be an open and simply connected region of the plane and consider the dynamical system given by Let the functions f and g be continuously differentiable on D and suppose that there is a continuously differentiable function ϕ on D (called a Dulac function) such that (∂/∂x 1 )(ϕf ) + (∂/∂x 2 )(ϕg) is either strictly positive or strictly negative almost everywhere on D. Then D does not contain any periodic solution of system (10).
We are now able to state and prove our stability result concerning a coexistence equilibrium.
Proof: Under the assumptions of this theorem, the global stability of E * follows from ruling out periodic solutions using the Dulac criterion with the following Dulac function .
Using assumption (C2), the differential expression appearing in the Dulac criterion evaluates as Therefore, no periodic solution exists in the first quadrant for the model (9) and E * is thus globally asymptotically stable. Note also that if the assumption (C3) holds for all positive x 1 and x 2 , the existence of a unique coexistence equilibrium E * implies that all solutions of the model (9) are uniformly bounded, by Lemma 2.1 and Theorem 2.2.

Remark 2.1:
In some sense, α is a critical ratio, although in this context one may find a continuum of critical ratios rather than a single one. If the ratio between densities of the second and first species equals α, both species experience a decrease at large population densities, since negative rates surpass positive ones. If this ratio falls below α, to some β < α, then, due to the monotonicity of R 1 , the first species, which is now 'excedentary', will have to decrease its density, since R 1 (β) ≤ R 1 (α) < 1. Similarly, if this ratio goes above α, to some γ > α, then due to the monotonicity of R 2 , it will be the turn of the second species to be 'excedentary' and to decrease its density, since R 2 (1/γ ) ≤ R 2 (1/α) < 1. That is, the inequalities R 1 (α) < 1 and R 2 (1/α) < 1 help create a balance mechanism, which keeps unboundedness in check.
The following result shows that the boundedness condition from the previous theorem is necessary to guarantee that all solutions are bounded (and thus to avoid the unrealistic reinforcing of the mutualistic effects). Theorem 2.6: Let the assumptions (C1) and (C3) hold. If there is α > 0 such that R 1 (α) > 1 and R 2 (1/α) > 1, then solutions of the model (9) are unbounded if the initial population density is large enough.

Proof:
We have already seen that solutions of the model (9) are bounded from below by a strictly positive constant. Hence we can assume that there exists ξ > 0 such that x 1 (t) > ξ and x 2 (t) > ξ. Seeking for a contradiction, we assume that the function is bounded from above, i.e. h(t) < M for a certain M > 0 and for any t. With these assumptions we have the following two cases for each t > 0.
Case 1. h(t) = αx 1 (t). One then has, by (M) and the first part of (C1) Assuming that x 1 (0) > M 1 > 0 the last term is positive and decreasing in x 1 . Passing to the limit as x 1 → ∞ and using (C3), it follows that Case 2. h(t) = x 2 (t). One then has, again by (M) and the first part of (C1), that Similarly, if x 2 (0) > M 2 > 0, passing to the limit as x 2 → ∞ in the last term and using (C3), it follows that This shows that h /h is bounded from below by a strictly positive constant C > 0 with as long as h(0) > max{αM 1 , M 2 }. While h is not differentiable everywhere in the classical sense, as an increasing function it is differentiable almost everywhere and the following weak version of the Fundamental Theorem of Calculus holds We thus have Letting t → ∞ in this inequality, it follows that h(t) → ∞. This contradicts the assumption that h is bounded from above. Hence h and, therefore, both x 1 and x 2 are unbounded.

Remark 2.2:
It can be seen that the conditions for boundedness and unboundedness are mutually exclusive. To this purpose, note that, since R 1 and R 2 are increasing Consequently, if there is any α such that R 1 (α) < 1 and R 2 (1/α) < 1, there cannot be any β such that R 1 (β) > 1 and R 2 (1/β) > 1, and vice versa.

Remark 2.3:
The boundedness conditions derived in our theorems can be expressed in the form of a single threshold for a broad class of mutualism models, including most of the models already analysed in the literature. For example, this is always the case when the limits R 1 (α) and R 2 (α) come from rational expressions, having the form Then the boundedness condition is equivalent to finding α > 0 such that If one of the constants C 1 , C 2 is 0, then the boundedness condition is satisfied regardless of the value of the other constant. For example, if C 1 = 0 then R 1 (α) = 0 for any α and R 2 (1/α) can be made less than 1 by choosing α large enough. Similarly, if C 2 = 0 then R 2 (1/α) is always zero and R 1 (α) can be made less than 1 by choosing α small enough. If both R 1 and R 2 are identically zero, then the boundedness follows without any condition on the existence of a suitable α.
Finally, if C 1 , C 2 > 0, then condition which is satisfied if C Similarly, the unboundedness condition is equivalent to finding α > 0 such that Obviously, the unboundedness condition cannot be satisfied if any of C 1 , C 2 is null. If C 1 , C 2 > 0, then condition which is satisfied if C then the boundedness condition can be expressed as R < 1, while the unboundedness condition can be expressed as R > 1.

Mutualisms that affect the growth function
We now show that the approach used to discuss stability of Equation (9) can be easily adapted for the situation in which the mutualism affects the growth function a rather than the self-limiting function f. The main theorems presented above still hold and, since their proofs are very similar, we only state here the specific assumptions, thresholds and the Dulac function to be used in the argument for non-existence of periodic solutions.
The consistency assumptions are now as follows: are eventually decreasing for all α > 0, and the mutualism-specific assumption is Hence, the limits R 1 (α), R 2 (α) are now given by and the Dulac function corresponding to the model (13) is ϕ(x 1 , x 2 ) := 1/x 1 x 2 a 1 (x 1 ) a 2 (x 2 ), which leads to for any positive x 1 and x 2 .

Models with asymmetric mutualism effects
We can also have a combination of the situations shown in the models (9) and (13), in which x 2 reduces the death rate of x 1 and x 1 increases the birth rate of x 2 , which allows for the modelling of interacting species with very different vital parameters, as follows: Here, the assumptions on a 1 and f 1 are those used in the analysis of model (9) while for f 2 and a 2 we consider the assumptions from Equation (13). Again, the bounding theorems presented above hold with very similar proofs and the Dulac function is ϕ(x 1 ,

Mutualism model with strong Allee effects
We now assume that a strong Allee effect affects one of the species involved in a mutualistic interaction. Specifically, we shall replace the growth assumption (G) with the following strong Allee effect assumptions: This ensures that a species with a strong Allee effect goes extinct if its density drops below the Allee threshold A and that A is now an unstable equilibrium of the model (8). If both species are subject to a strong Allee effect, then it is easy to see that the extinction state will always be a locally stable attractor. Hence, global stability of any coexistence equilibrium can only be addressed if just one of the species exhibits a strong Allee effect. In this case, the benefit this species gets from mutualism can be sufficiently strong to prevent its extinction when below A. We will establish a general condition for this to occur.
Let us consider the model (9) and assume that x 1 undergoes a strong Allee effect with the Allee threshold A 1 and the carrying capacity K 1 , while x 2 obeys the growth assumption (G) with the carrying capacity K 2 . In the absence of mutualism this system will have the following five equilibria (0, 0), (A 1 , 0), (K 1 , 0), (0, K 2 ), (A 1 , K 2 ), and (K 1 , K 2 ). Of these, (0, K 2 ) and (K 1 , K 2 ) are stable. In the presence of mutualism, assuming all other consistency conditions (C1), (C2) and (C3), if there exist a unique coexistence equilibrium of the model (9), then this equilibrium will be globally stable provided that all solutions are uniformly bounded in the positive quadrant. The boundedness arguments shown in the main two theorems of the previous section hold also here, with the caveat that we need to establish a condition which ensures that there is no limit point for any solution on a semi-axis.
While x 2 cannot go extinct due to the growth assumption (G), x 1 may approach zero because of the strong Allee effect assumptions (S1)-(S4). This may happen at the boundary equilibrium (0, K 2 ). Evaluating the Jacobian of model (9) at this equilibrium we obtain ⎛ We see that the above matrix has two eigenvalues The second eigenvalue is negative by the growth assumption (G). Hence, (0, K 2 ) is unstable provided that This represents the necessary and sufficient condition that species x 2 prevents the extinction of species x 1 and annihilates its strong Allee effect. If all other necessary conditions for the theorems in the previous section are met, then the coexistence equilibrium is globally stable.

Application to specific mutualism models
Here we apply our general results of the previous section to several specific, biologically motivated mutualism models, thus demonstrating broad utility of our approach.

Wolin-Lawlor model with weak Allee effect
Weak Allee effects are thought to be relatively common in nature, so it is natural to consider mutualism models with one or both species possessing a weak Allee effect. The following modification of model (2) allows for coverage of weak Allee effects In the absence of mutualism (i.e. for b 21 = b 12 = 0), Equation (16) decouple and each of them is of type (7). Therefore, if any of the exponents p 1 , p 2 is greater than 1, then the respective species is subject to a weak Allee effect, while if both exponents are equal to 1, then there is no Allee effect and the model (16) reduces to the model (2). It is easy to see that existence and uniqueness of a coexistence equilibrium E * of model (16) is assured if and only if b 12 b 21 < 1. If this is the case, then and consequently This fits the framework of the general model (4) with R 1 (α) = b 12 α, R 2 (α) = b 21 α and R = b 12 b 21 . As a consequence, if R < 1 then E * is globally stable.

q-Logistic Wolin-Lawlor model with weak Allee effect
The classical logistic model of population growth with a linear per capita population growth rate is an idealization that is only rarely met in nature. Rather, relationships between the per capita population growth rate and population density are commonly non-linear, with convex relationships being more numerous [27]. A common model used to describe such non-linear relationships is x = rx(1 − (x/K) q ) for a positive parameter q [27]. We thus consider the following mutualism model with the generalized logistic term and possibly with a weak Allee effect assuming that K i > 0, p i ≥ 1, q i > 0 and n i > 0 for i = 1,2. For this model, we have and consequently The consistency condition (C3) here requires that n 1 ≤ q 1 and n 2 ≤ q 2 . If n 1 < q 1 then R 1 (α) = 0 for all α > 0, while if n 1 = q 1 we get R 1 (α) = b 12 α q 1 . Similarly, if n 2 < q 2 then R 2 (α) = 0 for all α > 0, while if n 2 = q 2 then R 2 (α) = b 21 α q 2 . It then follows that if n 1 < q 1 or n 2 < q 2 then the boundedness condition is satisfied for all other parameter values, while if n 1 = q 1 and n 2 = q 2 the boundedness condition is satisfied provided that We provide an illustrative numerical example. Let p 1 = p 2 = 1 + 10 −4 , q 1 = n 1 = 3, and q 2 = n 2 = 1. For K 1 = 15, K 2 = 35, b 12 = 0.75, and b 21 = 0.7 we see that R = 0.64 and there is a unique equilibrium, which is globally stable. On the contrary, for K 1 = 40, K 2 = 55, b 12 = 2, and b 21 = 1 we have R = 1.38 and the solutions are unbounded. In Figure 1 we show the phase portraits corresponding to these two situations, with r 1 = 3 × 10 −5 and r 2 = 1.05 × 10 −4 in panel (A) and r 1 = 8 × 10 −5 and r 2 = 1.65 × 10 −4 in panel (B). Note that the parameters r 1 and r 2 do not affect equilibrium stability but we set their values in order to plot the respective phase portraits.
In Figure 2 we show the region of a selected two-parameter space in which the solutions of the model (17) are bounded and we have a unique globally stable equilibrium. The shaded region corresponds to R < 1. In panel (A), we compare the relative effect of the terms describing the mutualism effect, b 12 and b 21 . In panel (B), we compare the effect of a limiting resource term of one species (q 1 ) with the mutualistic effect of the other species (b 12 ). All other parameters are as considered for Figure 1.
We note that in the case p 1 = p 2 = 1 we obtain a q-logistic Wolin-Lawlor model with no Allee effect. If we further set q 1 = q 2 = n 1 = n 2 = q > 0, there is a unique positive equilibrium E * = (x * 1 , x * 2 ), with Figure 1. Phase portrait of model (17).

Model with non-polynomial thresholds R 1 and R 2
For most published mutualism models, the limits R 1 and R 2 are powers of α. Nevertheless, our generalization does not require this particular form and for illustrative purposes we consider one example for which these limits have a different form, as follows: with K 1 < 1 and K 2 < 1. Then The existence of a real number α such that R 1 (α) < 1 and R 2 (1/α) < 1 is now equivalent to R := ln K 1 ln K 2 < 1 and the solutions of model (18) are unbounded if the opposite inequality holds. For instance, if K 1 = 0.3 and K 2 = 0.55 then R = 0.72, the solutions are bounded and we have a unique globally stable equilibrium (Figure 3(a)). On the other hand, if K 1 = 0.23 and K 2 = 0.43 and all the other parameters are the same then R = 1.24 and all solutions are unbounded (Figure 3(b)). We have r 1 = 0.002 and r 2 = 0.0033 in panel (A) and r 1 = 0.0026 and r 2 = 0.0042 in panel (B). In Figure 4 we plot the region of global stability when K 1 and K 2 vary.  Global stability region for model (18), in which R(K 1 , K 2 ) < 1.

Model with mutualistic effects at low densities
In the previously considered mutualistic models (except Equation (18)), the mutualistic effect can be factored such that f i (x 1 , x 2 ) is separated as a product of a function depending only on x 1 and a function depending only on x 2 . Our general results allow for models in which one cannot perform such a factorization. The following model exemplifies this situation and has a biological appeal in describing the case when mutualism takes effect if the beneficiary species is at low densities and wanes as it grows. More specifically, the larger x 1 is the lower is its benefit gained from the presence of x 2 and vice versa. It is seen that and consequently The consistency condition (C3) requires in this case that q 1 ≥ n 1 and q 2 ≥ n 2 . However, R 1 (α) = R 2 (α) = 0 for all α > 0, and the boundedness condition is thus satisfied. This is consistent with our previous interpretation: if mutualism wanes when the beneficiary species is at larger densities, there is less need to impose a further condition to prevent unboundedness.
Since the growth assumption (G) implies that R 1 (0) < 1 and R 2 (0) < 1, using the ideas of Theorems 2.2 and 2.6 and Remark 2.3 one may obtain the following result:

Models with saturating functional response for mutualistic interaction
TheLotka-Volterra mutualism model (1) is a model of mutualism that fits the general formulation (3). Not surprisingly, our general theory gives for it the boundedness (and global stability) condition b 12 b 21 < 1, which is also the existence condition for its unique positive equilibrium. Here we consider another model that fits the general formulation (3). It is a modification of the Lotka-Volterra model that considers a saturating functional response for the mutualistic interaction. This model is commonly used in theoretical studies on mutualistic interactions and assumes that the effects of mutualists saturate as their density becomes high [19,31]. Here we consider a version proposed in [14] which has the form Consider a species subject to a strong Allee effect that follows model (6) when alone. A possible mutualism model is then In this setting and it is clear that condition (15) cannot hold. On the other hand, our boundedness theorem still applies. Since it follows that model (23) has bounded solutions whenever On the other hand, if R > 1 then all solutions are unbounded except those starting in the basin of attraction of the boundary equilibrium (0, K 2 ), which describes presence of a strong Allee effect in spite of the mutualism. On the contrary, the following example describes a situation in which mutualism can overcome an existing strong Allee effect, condition (15) holds and there is a unique coexistence equilibrium that is globally stable. The mutualism model is as follows: In this case The per capita growth rate of the first species in the absence of mutualism, a 1 (x 1 ) − f 1 (x 1 , 0), has a critical point at x 1 = A 1 − 1. Hence, in order to satisfy the strong Allee effect assumptions (S1)-(S4), we need to set the parameters A 1 and B 1 so that This ensures that, in the absence of mutualism, the first species is subject to a strong Allee effects and has a hump-shaped per capita growth rate that equals zero at two different positive values of x 1 . Since  In addition, the strong Allee effect is eliminated in the first species if and only if condition (15) holds, which in this particular case requires that

model (24) has bounded solutions if and only if
We provide a numerical example to show that conditions (25) and (26) can indeed be satisfied simultaneously. Choosing A 1 = 3, B 1 = 7, b 12 = 0.1, and b 21 = 0.7 we have that R = 0.21. If K 2 = 3, then the strong Allee effect in the first species is maintained and the boundary equilibrium (0, K 2 ) remains locally stable ( Figure 5(a)). On the other hand, if K 2 = 13 then (0, K 2 ) becomes a saddle and the first species will never go extinct since mutualism prevents the strong Allee effect (Figure 5(b)). In Figure 6, we show the region given by the parameter that controls a vital rate of the first species (A 1 ) and the parameter describing the beneficial effect of the other species (b 12 ) in which the strong Allee effect in the first species is eliminated.

Conclusions
In this article, we extend the work of Georgescu et al. [12,13] and provide a general class of models of facultative mutualism that covers many of the commonly used mutualism models. Most importantly, we study boundedness of solutions of models in this general class which, from the biological perspective, is an important necessary condition for any mutualism model to be considered plausible. The main theorems we provide in this respect allow for a much greater flexibility in designing these models. Formulations of the mutualistic interaction can now be more complex (i.e. do not require a factorization in the state variables) and each species can be modelled under different assumptions (e.g. one can benefit from a reduction in mortality and the other from an increase in fertility). Our generalization is therefore suitable for modelling mutualistic systems in which the interacting species are very different with respect to their vital parameters but also to their demographic assumptions. In addition, our generalization covers Allee effects as a part of mutualism models. This is certainly relevant, since at small population sizes or low population densities Allee effects (and especially weak Allee effects) are known to be relatively common (see, for instance, [5]). Mutualism models with Allee effects form just a specific case of our general model class, so they can be treated without any additional overhead provided that some general assumptions are met. Therefore, the most important biological implication of the results established herein is that many more systems than before are now known to be stable and that the presented models may become more adequate descriptions of mutualisms in models of more complex food webs. Virtually all published mutualism models, including all we present in this article, follow a tradition set by Lotka and Volterra and have a logistic-like model or its Allee effect variant as their core element. Because of that, there is no clear distinction between birth and death rates in these models and both positive (or growth) and negative (or self-limiting) model terms comprise a piece of each. An alternative, more mechanistic approach would be to start with a model in which birth and death rates are explicit right from the beginning, such as where the positive functions b and d represent the per capita birth and death rates, respectively. Mutualistic effects could then be classified as those that affect births, deaths, or both terms. Note that the single-species model (8) we start from has the same structure as the model (27). Therefore, if any mutualistic extension of the model (27) satisfies our assumptions, the general results we derive would apply also to this class of models. In addition, with birth and death rates clearly defined, the ratios R 1 and R 2 which enter all our general results would technically represent births to deaths ratios at given population densities and would thus be somewhat akin to reproductive ratios that characterize infection dynamics in epidemiology or population dynamics in life history theory. See also [10] for further discussion on related matters. Some limitations of our work still exist, which may serve as potential avenues for further research. First of all, the consistency condition (C2) is far more restrictive than the similar condition (C1). This is due to our choice of the Dulac function used in our global stability argument. Presumably, a more elaborate Dulac function may allow a weaker condition in place of (C2). On the other hand, this may require more information about the structure of functions a i and f i and so there is likely to be a trade-off between how general these functions are and how less restrictive the Dulac function can be. In this paper, we place more importance on the solution boundedness results (required for biological plausibility) and on generality of mutualism models. While the derived conditions on bounded and unbounded solutions are mutually exclusive and in many cases threshold-like, the issue of uniqueness of a coexistence equilibrium is a separate question and does not follow directly from these conditions. It would certainly be of interest if a generalization were made so that uniqueness of a coexistence equilibrium follows from further suitable conditions on the right-hand sides of our generalized mutualism models. We expect to pursue these questions in a follow-up paper.