Global stability and persistence of complex foodwebs

We develop a novel approach to study the global behaviour of large foodwebs for ecosystems where several species share multiple resources. The model extends and generalize some previous works and takes into account self-limitation. Under certain conditions, we establish the global convergence and persistence of solutions.


Introduction
To mathematically show the existence and stability of large foodwebs, large and complex as foodwebs in nature, is still one of the key problems in theoretical ecology. A specific part of this theoretical issue is that many species can share just a few resources (for example ocean ecosystems including thousands of phytoplankton species) yet the competitive exclusion principle [8,22] asserts that such foodwebs should not exist. To partly explain that paradox [13] showed that a system consisting of a single resource and three species can support chaotic dynamics where all species coexist. Another explanation of the paradox was proposed in [19] where self-limitation effects has been taken into account.
In this paper we add complexity to the work of [13] and [19] by extending the dynamical equations considered in [19] with self-limitation effects [19,2,3,1,14,15]) (see also a turbidostat model in [16]). We obtain a complete description of the large time behaviour of the system. In particular, we explore the range in the parameter space that leads the system to the global stable equilibrium point. Furthermore, we show that if the self-limitations exceed some critical values then the system exhibits either global stability or persistence, see Corollaries 7.2 and 7.3.
Traditionally, the Lyapunov function approach is used to establish global stability, see a recent review in [10]. In our case, however, an explicit information about equilibrium points is not available. Instead, we transform our problem to a finite dimensional nonlinear fixed point problem for an appropriate non-increasing operator. We show that the asymptotic behaviour of a generic solution to the initial problem is well controlled by iterations of the introduced operator. This allows us to derive explicit a priori estimates (see Theorem 7.1 below) and the global stability.
The paper is organized as follows. In Section 2, we present the model with self-limitations and in Section 3 we obtain some preliminary results. We review some elementary facts on period-two-points of non-increasing maps in Section 4 and discuss the structure and stratification of equilibrium points in Section 5. In particular, in Section 6 we consider the so-called special equilibrium points which are significant for the large time behaviour of the original dynamical system. Here we also define the corresponding finite dimensional fixed point problem. To study its dynamics and convergence we need to suitably polarize the fixed point problem. This allows us to establish bilateral estimates for the corresponding ω-limit set. The main result of this section is contained in Corollary 6.3, which gives a sufficient condition for the existence of a unique fixed point. In Section 7 we return to the main dynamical system formulate and prove the main results on the large behaviour of the original dynamic system. In particular, we obtain some explicit conditions when the system obeys strong persistence. Finally, in Section 8 we briefly discuss our results and relate them to some previous research.

The model
Given x, y ∈ R n we use the standard vector order relation: is the closed box with vertices at a and b.
We consider the model where the population dynamics of M species competing for m complementary resources is governed by chemostat-like equations Here x j (t) is species j abundance and v i (t) is the concentration of resource i at time t, µ j are the species mortalities, S i is the supply of resource i, c ij > 0 is the content of resource i in species j (growth yield constants), D i is the rate of exchange of resource i, resource turnover (or dilution) rate), γ j > 0 is a self-limitation constant of species j. We shall assume that the specific growth rates φ j are bounded Lipschitz functions subject to the following standard conditions: The most relevant for biological applications example of the specific growth functions φ j is given by the Monod equation and Liebig's 'law of the minimum' where r j is the maximum specific growth rate of species j and K ij is the halfsaturation constant for resource i of species j. Obviously, the functions (6) meet the above conditions. In absence of self-limitation (γ j = 0), the present model naturally appears in the bioengineering context [4] and was extensively studied for m ≤ 2 resources for both equal resource turnover rates µ j = D i = D in [11], [13], [12] and for different the removal rates µ j in [9], [10], [18], [21], see also a recent review in [20]. For a single resource m = 1, the dynamics of the standard model in absence of selflimitation is completely determined by the break-even concentrations R j defined as φ j (R j ) = µ j , see [4], [10]. For example, if the lowest break-even concentration achieves on a single species k then while lim t→∞ x j (t) = 0 for all j = k. However, if m ≥ 3, the behaviour becomes much more involved. Recent numerical simulations [13], [12] strongly support the possible chaos scenario for (m, M ) = (3,6) or (5,6). An important step was done by Li [17] who established the existence of the limit cycle for m = M = 3.

Preliminaries
In what follows, we shall assume that γ j > 0.
Proof. Note that by (1), x j (t) never vanishes unless x j (t) ≡ 0. In particular, by (1) x(t) ≫ 0 as long as x(t) is defined. Furthermore, if h i (x, v) denote the right hand side of (2) then by (4) > 0 for all admissible t > 0, see Proposition 2.1 in [6]. Similarly, h i (x, S) < 0 (unless x = 0) and v(0) ≤ S yields v(t) ≤ S, and thus (3) and (5) imply 0 ≤ φ j (v) ≤ φ j (S). This proves that R M + × [0, S] is an invariant subset for (1), (2), (3) . Furthermore, x j (t) ≤ y i (t), where y i (t) is the solution of the Cauchy problem This readily yields (8) and Proposition 3.1 shows that the extinction dynamics of (1), (2), (3) depends on the sign of φ j (S) − µ j : for species i to survive, its specific growth rate φ j (S) at the supply point S must exceed its specific mortality rate µ j . To eliminate the trivial extinctions, we shall assume in what follows that the survivability condition holds: (10) φ j (S) > µ j for all j.
Below we summarize some elementary observations which will be used throughout the paper.
Proof. An idea of the proof is clear from the figure below.
Then for any solution of Proof. By (b) F (z, t) is strictly decreasing in z for each t ≥ 0. It follows from the conditions (a)-(b) and the classical Clarke result [5] that z(t) in (c) is well-defined and local Lipschitz on [0, ∞). It follows from (a) that 0 < u(t) < S for all t ≥ 0. Now, two alternatives are possible: (i) either there exists T > 0 such that u(t) = z(t) for t ≥ T , or (ii) there exists t k ր ∞: u(t k ) = z(t k ). First let (i) hold and assume without loss of generality that u(t) < z(t) for t ≥ T . Then Combining (12) with the monotonicity of u(t) and (11) imples which implies the equality in (12). Next, if (ii) holds then lim k→∞ u(t k ) =z. Assume by contradiction that, for example,ū := lim sup t→∞ u(t) >z and let ξ k ր ∞ be a corresponding sequence where the lim sup is attained. Since one can redefine the sequence ξ k such that each ξ k becomes a local maximum of u.
Passing to limit as k → ∞ yields a contradiction.

Period-two-points of non-increasing maps
Let 0 ∈ D ⊂ R n + and G : D → D be an arbitrary map. Recall that a pair (a, b), a, b ∈ D, is called a period-two-point [7, p. 387 Any fixed point c ∈ Fix(G) gives rise to a trivial period-two-point (c, c).
Hereinafter, we assume that G is continuous and non-increasing in D, i.e. G(x) ≥ G(y) for any x ≤ y in D. Note that G is then automatically bounded: Since 0 ∈ D, the iterations u 0 = 0, u k := G k (0) ∈ D, k ≥ 1, are well-defined, u 1 ≥ u 0 = 0 (an a priori estimate) and u 2 = G(u 1 ) ≤ G(u 0 ) = u 1 (by virtue of the monotonicity of G). Hence, it follows by induction that This implies that the limits Thus obtained period-two-point is extremal as the following property shows.
In particular, and the box Proof. Since a ≥ u 0 = 0 and G is a non-increasing one has This readily yields (18). Then (19) follows from the fact that (c, c) is a period-twopoint for any c ∈ Fix(G). The last claim of the proposition follows immediately from the monotonicity of G and (17).
Proof. Let y 0 = y and y k = G k (y), k ≥ 1, hence (21) becomes Applying G we yields y 2 ≥ x 1 ≥ y 0 and x 0 ≥ y 1 ≥ x 2 . Proceeding by induction on k, we obtain by virtue of (14) This implies the existence of limits in (22). It also follows that G(a) = b and G(b) = a, hence (a, b) ∈ Fix 2 (G) and a ≤ x, y ≤ b. Combining with the extremal property (18) yields (22).

Stratification of equilibrium points
Let us denote by E the set of nonnegative equilibrium points (stationary solutions) of (1)- (2). It is natural to consider the standard stratification By the above, there exists x k = 0, therefore φ k (v) = 0 implying by (4) that v ∈ ∂R m + , thus φ j (v) = 0 for all j. Applying the stationary condition to (2) The converse is not necessarily true: if v is a solution of (25) and x is defined by (26) then (x, v) is an equilibrium point in E J ′ for some J ′ ⊂ J. Indeed, it might happen that φ j (v) ≤ µ j , i.e. x j = 0 for some j ∈ J. On the other hand, if J = ∅ then necessarily J ′ = ∅ because if x j = 0 for all j then (x, v) = (0, S), but F J (S) ≪ S in view of (4), a contradiction with (25).
To distinguish this situation, we denote by E J = the set of solutions (x, v) of (25) and (26).
Then E ∅ = E ∅ , and the above argument yields that for any J = ∅ Thus refined stratification J → E J still contains information about all equilibrium points but it has better properties than J → E J .
By Proposition 5.2, the set of special equilibrium points is nonempty. Note also that if (x, v) is an arbitrary equilibrium point of (1)-(2) with x ≫ 0 then it is necessarily a special one because by (1) φ j (v) > µ j for all j, hence x is determined by (29) and therefore v satisfies (28).
The set of special equilibrium points E M = Fix(F) reflects the complexity of large-time dynamics of the original system in the following sense. Theorem 7.1 below shows that if there exists a unique global stable equilibrium point of (1), (2), (3) then it is necessarily a special point (in this case, obviously, unique). Therefore the structure and the number of special equilibrium points plays a crucial role in the large-time dynamics of (1), (2), (3) .
Thus, it is naturally to expect that the global stability will be lost if the cardinality | Fix(F)| ≥ 2. Note that if m = 1 then Lemma 3.2 easily implies that Ω consists of exactly one point: | Fix(F)| = 1. However, if m ≥ 2, the situation is more subtle as the example below shows.
A careful analysis shows that for m = 2 there always holds | Fix(F)| ≤ 3 Liebig-Monod model (6). Furthermore, for any m ≥ 2, an argument similar to Example 6.1 yields | Fix(F)| ≥ m + 1 for certain sets of parameters. Now, let us turn to the fixed point problem (28). It is naturally to study solutions of (28) by virtue of iterations F k (0). But (28) is non-regular in the sense that already the second iteration F 2 (0) can be outside of [0, S]. Indeed, To refine iterations, we suitably polarize (28) to get a system with the same set of fixed points. Namely, given w ∈ [0, S] let us define V(w) ∈ [0, S] as the unique solution v of the system . . , w m ), i = 1, . . . , m.
Note that each equation of system (31) contains a single unknown variable v i , thus Lemma 3.2 implies that for all i a unique solution v i of (31) exists and 0 < v i ≤ S. Therefore 0 ≪ V(w) ≤ S. Also, by the survivability condition (10) Furthermore, the second part of Lemma 3.2 implies that V(w) is non-increasing: Now, if v solves (28) then by the uniqueness of solution of (31) one has Conversely, if v is a solution of (33) then it also solves (28). Thus, in the present setting, the fixed point problem (28) is completely equivalent to (33): The main advantage of V with respect to F is that by its definition, Now, with V in hands we apply the technique of Section 4. Namely, using the definition (16), we see that starting with u 0 = 0, the even and odd iterations converge respectively to In particular, and, furthermore, (0 V ,0 V ) possesses the extremal property in Proposition 4.1. In particular, it follows from (19) that This imemdiately yields Conversely, (35) implies that the cardinality of fixed points | Fix(F)| is an obstacle for the coincidence relation (36). Furthermore, Example 6.1 above shows that for certain values of parameters of our system one has | Fix(F)| > 1, thus, one cannot expect in general the coincidence in (36). Therefore it is important to know when (36) holds. One such sufficient condition is presented below.
In general one has from (32) that V(S) ≪ S and V 2 (S) ≤ S, hence the following simple bilateral estimates hold: The latter estimate (42) is optimal in general. Indeed, if cij Diγj are large enough, V(S) can be made arbitrarily small, for instance such that φ j (V(S)) ≤ µ j , which yields V 2 (S) = S and therefore (0 V ,0 V ) = (V(S), S).
In particular, Here x ∨ y (resp. x ∧ y) denote the vector whose ith coordinate is max(x i , y i ) (resp. min(x i , y i )).

Bilateral estimates
As it was pointed out before, Example 6.1 shows that a priori the asymptotic behaviour of solutions to (1)-(2) can be rather complicated for m ≥ 2. On the other hand, the result below shows that the global dynamics is completely controlled by the finite-dimensional fixed point problem (28) and the characteristic parameters in (17).
Proof. As the first step we reformulate the original system as an appropriate integral equation for unknown function v(t). Let w(t) : [0, ∞) → [0, S] be a a continuous vector-function with w(0) = v(0) and having a limit lim t→∞ w(t) =w. Then (47) solves (1) with v(t) replaced by w(t). Clearly, X (w)(t) is a nondecreasing functional of w, X (w)(0) = x(0) and one can readily verify that Next, let V (w)(t) denote the solution u(t) of the system below (obtained from (2) with x(t) replaced by (47)): Let , ∀T > 0. Next, note that V (w) is a non-increasing functional of w. Indeed, let 0 ≤ w(t) ≤ w(t) for all t ≥ 0, and let u i (t) and u i (t) be the corresponding solutions of (49). Denote by F i (u i (t), t) and F i ( u i (t), t) the right hand side of (49) corresponding to w(t) and w(t) respectively. Then the F i (z, t) and F i (z, t) satisfies the conditions of Lemma 3.3 and u i (0) = u i (0), therefore u i (t) ≥ u i (t) for all t, as desired.
Combining the obtained estimates with Corollary 6.3 implies the following global stability result. Numerical simulations in [13] show that certain solutions of the standard model with γ j = 0 and m ≥ 3 have periodic (chaotic) dynamics. Corollary 7.2 above shows that if the self-limitation constants γ j or dilution rates D i are large enough, the global behaviour of the modified model becomes stable for any choice of m and M .
In fact, one can chose the parameters of the system such that the strong persistence holds, see the corollary below. To present our result we need to define an analogue of the lowest break-even concentration R * in (7)  Note that by (10), R = ∅.  Proof. By (10), S ∈ R, therefore the number δ := sup{t ≥ 0 : (S 1 − t, . . . , S m − t) ∈ R} is well defined and positive. Since δ ≤ S ∞ , we have ρ 0 := δ/ S ∞ ≤ 1. If ρ ≤ ρ 0 then by Corollary 7.2 any solution with Cauchy data (3) converges to a unique equilibrium point 0 ≪ ξ ≪ S satisfying (28). We have for all 1 ≤ i ≤ m Therefore ξ ∈ R, implying by (58) and (29) that lim t→∞ x j (t) = X j (ξ) > 0 for all j, as desired.

Discussion
In this paper we established sufficient conditions for the global stability and persistence of a chemostat like model with self-limitations. For the Liebig-Mondoc model (6) one has L j = r j / min i {K ij } and φ j (S) ≤ r j . It is interesting to compare our result with simulations in [13] rigourously approved in [17], see especially Section 5 there. In that example, Huisman and Weissing assume in the present notation that m = M = 3, S j = 10, r j = 1, D j = 0.25 for all three species and matrices K ij and c ij be chosen as in [17, p. 38]. Then if γ j = 0 then Theorem 3.1 in [17] implies the existence of a nontrivial periodic oscillation. On the other hand, it follows from (38) that if γ j ≥ 1.64, i = 1, 2, 3, then any solution is global stable, for arbitrary positive initial data. In general, given arbitrary data, (38) explicitly defines the critical values γ * j such that the system is globally stable for γ j > γ * j .