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 generalizes some previous works and takes into account self-limitation. Under certain explicit conditions, we establish the global convergence and persistence of solutions.


Introduction
To mathematically show the existence and stability of large foodwebs, large and complex foodwebs in nature are still one of the key problems in theoretical ecology. A specific part of this theoretical issue is that many species can share more than just a few resources (for example ocean ecosystems including thousands of phytoplankton species) yet the competitive exclusion principle [8,24] 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 [21] where self-limitation effects have been taken into account. B Vladimir G. Tkachev vladimir.tkatjev@liu.se Vladimir Kozlov vladimir.kozlov@liu.se In this paper, we add complexity to the work of [13,21] by extending the dynamical equations considered in [21] with self-limitation effects [1][2][3]16,17,21]) (see also a turbidostat model in [18]). 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, see the explicit estimates in (37) and (38). Furthermore, we show that if the self-limitations exceed some critical values then the system exhibits either global stability or persistence, see Propositions 9 and 10.
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 1) and the global stability.
The paper is organized as follows. In Sect. 2, we present the model with self-limitations, and in Sect. 3, we obtain some preliminary results. We review some elementary facts on period-two-points of non-increasing maps in Sect. 4 and discuss the structure and stratification of equilibrium points in Sect. 5. In particular, in Sect. 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 Proposition 7, which gives a sufficient condition for the existence of a unique fixed point. In Sect. 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. In Sect. 8, we briefly discuss our results and relate them to some previous research. We finally mention that in our recent paper [14], we apply the results of the present paper to obtain explicit estimates of biodiversity for competition systems with extinctions.

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 i j > 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 i j is the half-saturation constant for resource i of species j. Obviously, the functions (6) meet the above conditions. In the 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][12][13] and for different removal rates μ j in [9,10,20,23], see also a recent review in [22]. For a single resource m = 1, the dynamics of the standard model in the absence of self-limitation is completely determined by the breakeven 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 [12,13] strongly support the possible chaos scenario for (m, M) = (3, 6) or (5,6). An important step was done by Li [19] who established the existence of the limit cycle for m = M = 3.

Preliminaries
In what follows, we shall assume that γ j > 0.
, v(t)) of (1), (2), (3) is well-defined and bounded for all t ≥ 0 and If φ j (S) = μ j , (8) should be replaced by 0 ≤ x j (t) ≤ 1+γ j x j (0)t . In particular, [0, S] is an invariant subset. Furthermore, if φ j (S) ≤ μ j for some i then lim t→∞ x j (t) = 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) This readily yields (8) and M i=1 c ki x j φ j (v) ≥ 0 yields the upper estimate in (9). Since Proposition 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: For the Monod-Liebig model (6), the survivability condition (10) is equivalent to 0 . . , R m j ) and R i j := μ j r j −μ j K i j are the resource requirement of a species j for a resource i [12].
Below we summarize some elementary observations which will be used throughout the paper.

Lemma 1 Let f (x), g(x) ≡ 0 be continuous nonnegative and non-decreasing maps
Proof An idea of the proof is clear from the figure below.
Then for any solution of 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: . First let (i) hold and assume without loss of generality that u(t) < z(t) for t ≥ T . Then hence u(t) is non-decreasing, therefore there exists Combining (12) with the monotonicity of u(t) and (11) implies 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. This yields 0 = u (ξ k ) = F(u(ξ k ), ξ k ), thus u(ξ k ) = z(ξ k ). 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 limits0 Thus obtained period-two-point is extremal as the following property shows.
In particular,0 and the box is invariant under the mapping G.
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-two-point for any c ∈ Fix(G). The last claim of the proposition follows immediately from the monotonicity of G and (17).
Then there exists (a, b) ∈ Fix 2 (G) such that 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 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 and J runs over all subsets of {1, 2, . . . , M}. The supply point S is the equilibrium resource availabilities in the absence of any species and obviously (0, S) is the only point in E ∅ :

Proposition 4 For an arbitrary
Proof If x = 0 then v = S, thus x > 0. If some v i = 0 then (4) yields φ j (v) = 0 for all j, hence by (2) Applying the stationary condition to (2) we see that v = S, a contradiction with v ∈ ∂R m + . Therefore, v S. and 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).
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 .
and x is given by By Proposition 5, 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 1 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 play 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 1 easily implies that Fix(F) consists of exactly one point: | Fix(F)| = 1. However, if m ≥ 2, the situation is more subtle as the example below shows (see also [15]).
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 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, F 2 i (0) = F i (S) becomes negative if γ j or D i are small enough (alternatively, c i j large enough).
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 Note that each equation of system (31) contains a single unknown variable v i , thus Lemma 1 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 1 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 Sect. 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 2. In particular, it follows from (19) that0 This immediately yields holds then there exists a unique special equilibrium point, i.e.
Conversely, (35) implies that the cardinality of fixed points | Fix(F)| is an obstacle for the coincidence relation (36). Furthermore, Example 1 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 c i j D i γ 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).

Proposition 8 For any
In particular, Here x ∨ y (resp. x ∧ y) denote the vector whose ith coordinate is min(x i , y i ) (resp. max(x i , y i )).
. . , w m ), it follows from the monotonicity of F and (31) that which yields the right inequality in (43). The left one follows by a similar argument from w ∧ v ≥ (w 1 , . . . , w i−1 , v i , . . . , w m ) and the fact that V(w) ≥ 0. Then (44) follows from 0 V(S) ≤ V 2 (S) ≤ S and (43).

Bilateral estimates
As it was pointed out before, Example 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). (x(t), v(t)) be the solution of (1), (2), (3). Then in notation of Sects. 4 and 5:

Theorem 1 Let
In particular, if0 V =0 V then all solutions of (1)-(2) converge to a unique special equilibrium point.
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 continuous vector function with w(0) = v(0) and having a limit lim t→∞ w(t) =w. Then solves (1) with v(t) replaced by w(t). Clearly, X (w)(t) is a non-decreasing function 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 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) satisfy the conditions of Lemma 2 and u i (0) = u i (0), therefore u i (t) ≥ u i (t) for all t, as desired. Furthermore, we claim that Indeed, rewrite (49) as Then Comparing the latter expression with (31), we conclude that z = V i (w) is the unique root ofF i (z) = 0 in [0, S]. Now, let 0 ≤ z i (t) < S be the unique solution of F i (z i (t), t) = 0, t ≥ 0. Suppose that t k ∞ realizesz := lim sup t→∞ z i (t). Then Similarly one shows that V i (w) = lim inf t→∞ z i (t). Thus, lim t→∞ z i (t) = V i (w) exists, as desired. Applying Lemma 3 yields (51).
In the present setting, if (x(t), v(t)) is the solution of (1), (2), Now we show that v(t) can be obtained as the limit of iterations As V is non-increasing and V (v) = v, one has Since V 2 is non-decreasing and by (50) T ] is compact, hence both the odd v 2k−1 (t) and even v 2k (t) terms converge in C 1 [0, T ], therefore the following limits are well-defined for any t ≥ 0: and Since v 0 = 0 ≤ v we also have by (52) and (x,v) and (x,v) solve, respectively, Taking the difference yields that (ξ, η) := (x −x,v −v) satisfies a homogeneous system of ODEs with bounded coefficients (recall that φ j are Lipschitz). Since (ξ, η) has the zero Cauchy data, we conclude by uniqueness for the Cauchy problem and (56) thatv(t) =v(t) = v(t) andx(t) =x(t) = x(t). In summary, for any fixed t > 0 one has Numerical simulations in [13] show that certain solutions of the standard model with γ j = 0 and m ≥ 3 have periodic (chaotic) dynamics. Proposition 9 shows that if the selflimitation 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 choose 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) for general response functions φ j . Let us consider the set R := {v ∈ [0, S] : φ j (v) > μ j for all j}. (58) Note that by (10), R = ∅.
In general, one has from (45), (42) and (44) the following explicit a priori estimate.

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 i j } and φ j (S) ≤ r j . It is interesting to compare our result with simulations in [13] rigorously proved in [19], 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 i j and c i j be chosen as in [19, p. 38]. Then if γ j = 0 then Theorem 3.1 in [19] implies the existence of a nontrivial periodic oscillation. On the other hand, it follows from (38) that if γ j ≥ 1.64, j = 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 .