Biodiversity and robustness of large ecosystems
Vladimir Kozlov, Sergey Vakulenko, Uno Wennergren and Vladimir Tkachev
The self-archived postprint version of this journal article is available at Linköping
University Institutional Repository (DiVA):
http://urn.kb.se/resolve?urn=urn:nbn:se:liu:diva-150011
N.B.: When citing this work, cite the original publication.
Kozlov, V., Vakulenko, S., Wennergren, U., Tkachev, V., (2018), Biodiversity and robustness of large ecosystems, Ecological Complexity, 36, 101-109. https://doi.org/10.1016/j.ecocom.2018.07.003
Original publication available at:
https://doi.org/10.1016/j.ecocom.2018.07.003
Copyright: Elsevier
Biodiversity and robustness of large ecosystems
Vladimir Kozlova, Serge Vakulenkob, Uno Wennergrenc, Vladimir Tkacheva,∗
aDepartment of Mathematics, Link¨oping University, Sweden
bInst. for Mech. Eng. Problems, Russian Acad. Sci., St. Petersburg, ITMO University, St. Petersburg and Saint Petersburg
State University of Industrial Technology and Design
cDepartment of Physics, Chemistry, and Biology, Link¨oping University, Sweden
Abstract
We study the biodiversity problem for resource competition systems with extinctions and self-limitation effects. Our main result establishes estimates of biodiversity in terms of the fundamental parameters of the model. We also prove the global stability of solutions for systems with extinctions and large turnover rate. We show that when the extinction threshold is distinct from zero, the large time dynamics of system is fundamentally non-predictable. In the last part of the paper we obtain explicit analytical estimates of ecosystem robustness with respect to variations of resource supply which support the R∗ rule for a system
with random parameters.
Keywords: Foodweb, Biodiversity, Global stability, Extinction threshold, Ecological networks, the R∗rule
1. Introduction
1
Existence and stability of large foodwebs, where many species share a few of resources, is one of key
2
problems in ecology [4, 21, 9] as well as extinctions and mass extinctions in such systems under climate
3
changes [14]. In this paper, we consider the model initiated in [11], [10] describing an ecological system,
4
where several (many) species compete or fight for few limited resources.
5
The most typical examples are plant or plankton ecosystems. Sunlight, water, nitrogen, phosphorus
6
and iron are all abiotic essential resources for phytoplankton and plant species. Resource competition
7
models link the population dynamics of competing species with the dynamics of the resources. As it was
8
mentioned in [7] an attractive feature of resource competition models is that they use the biological traits
9
of species to predict the time evolution of competition. In fact, many rigorous results [5, 19, 16] show
10
that, in general situation, a single species survives and to obtain coexistence of many species one needs
11
very special assumptions to species parameters (mortalities and resource consumption rates). This paradox
12
(the so-called paradox of plankton [9, 4]) has received a great attention in past decades [17, 15]. Several
13
ways to explain the extreme diversity of phytoplankton communities have been proposed. In particular, the
14
proposed mechanisms include spatial and temporal heterogeneity in physical and biological environments,
15
horizontal turbulence of ocean, oscillation and chaos generated by several internal and external causes, stable
16
coexistence and compensatory dynamics under fluctuating temperature in resource competition, and
toxin-17
producing phytoplankton [17,15]. Although the mathematical problem has been studied for more than two
18
decades it is still far from to be well-understood. The most of available results do not give explicit estimates
19
of biodiversity in terms of the fundamental observable ecosystem parameters (such as species mortality rates,
20
rates of resource consumptions, resource supply and resource turnover rate).
21
The main goal of this paper is to present such estimates. To this end we consider dynamical equations
22
are close to the model equations, which considered in the seminal paper [7] but extend that model in the two
23
∗Corresponding author
Email addresses: vladimir.kozlov@liu.se (Vladimir Kozlov), vakulenfr@gmail.com (Serge Vakulenko), Uno.Wennergren@liu.se (Uno Wennergren), vladimir.tkatjev@liu.se (Vladimir Tkachev)
aspects. First, we take into account self-limitation effects (which are important for plankton populations
24
[15] and to explain stability of large foodwebs [1, 2]). Roughly speaking when we introduce a weak
self-25
limitation we replace equations of Maltus type on Verhulst type equations. Second, following [11] we take
26
into account species extinction thresholds, however, in contrast to [11, 18] we consider here the case of a
27
few resources. Mathematically, our approach with extinction thresholds and self-limitation terms can be
28
considered as a regularization of resource competition models.
29
Our main results can now be formulated as follows. A summary of the mathematical framework and the
30
global stability results established earlier in [12], [10] for the model with zero extinction threshold is collected
31
in Sections2 and3. In that case, a complete description of the system large time behaviour is obtained for
32
systems with sufficiently large turnover rates and without extinctions. More precisely, the model exhibits
33
the global stability: all positive trajectories converge to the same equilibrium state. This result holds due to
34
two principal properties of our system. First, the system has a typical fast/slow structure for large turnover
35
rates. Second, the system obeys a monotonicity property: if resources increase then species abundances also
36
increase. We recall the principal ideas of the proof at the end of section 3.2.
37
Next, if one allows even small positive extinction threshold, the ecosystem behaviour exhibits new
inter-38
esting effects. We study this in sections 4 and 5 below. We establish a weaker stability result: the limit
39
equilibrium state still exists but it depends on the initial ecosystem state. This in particular implies that
40
there can a priori exist several distinct equilibrium states.
41
In section 5, we establish explicit upper and below estimates of biodiversity expressed in terms of the
42
fundamental ecosystem parameters (such as species mortalities, resource consuming rate etc.). Remarkably,
43
the obtained estimates are universal for small extinction thresholds and self-limitation parameters. We point
44
out that these results use no assumptions on the system dynamics and do not use our theorem on global
45
stability.
46
In the part of the paper, we study large ecosystems with random fundamental parameters. Here the
47
main assumption is that the system dynamics has no oscillating or chaotic regimes. Note that it follows from
48
Theorem3.2, that the assumption is automatically holds if, for example, turnovers rates are large enough.
49
Recall that the R∗ rule (also called the resource-ratio hypothesis) is a hypothesis in community ecology
50
that attempts to predict which species will become dominant as the result of competition for resources. It
51
predicts that if multiple species are competing for a single limiting resource, then species, which survive
52
at the lowest equilibrium resource level, outcompete all other species [20], [3]. In section 6 we obtain a
53
complete description of parameters for survived species and establish the validity of the R∗ rule for systems
54
with random parameters. We show that if the resources are limited and initially the number of species is
55
sufficiently large then only species with the fitness which is close to maximal one can survive. In our model,
56
the fitness is determined as the resource amounts available for an organism.
57
Finally, in section 7, we study sensitivity of those states with respect to a change of environmental
58
parameters. This allows us to essentially extend recent results of [14]. Namely, not only the magnitude
59
of environmental changes and their rates determine how much species will extinct but also the achieved
60
biodiversity level, and some species parameters. For example, ecosystems where the species parameters are
61
localized at some values are less stable than ecosystems with a large species parameter variation.
62
The basic notation
63
x(t) = (x1(t), . . . , xM(t)) the vector of species abundances
v(t) = (v1(t), . . . , vm(t)) the vector of resource abundances
µi, γi the mortality and the self-limitation constant of species i
Dj, Sj the turnover rate and the supply of resource vj
cij the content of resource j in species i
φi the specific growth rates of species i
Kij the half-saturation constant for resource j of species i, page3
Xext(i) the extinction threshold of species i, page5
Ne(t) the number of species which exist at the time t, page6
2. Preliminaries
64
Given x, y ∈ Rn we use the standard vector order relation: x ≤ y if x
i ≤ yi for all 1 ≤ i ≤ n, x < y if
x ≤ y and x 6= y, and x y if xi< yi for all i; Rn+ denotes the nonnegative cone {x ∈ Rn : x ≥ 0} and for
a ≤ b, a, b ∈ Rn
[a, b] = {x ∈ Rn: a ≤ x ≤ b} is the closed box with vertices at a and b.
65
We consider the following system of equations: dxi dt = xi(φi(v) − µi− γixi), i = 1, . . . , M, (1) dvk dt = Dk(Sk− vk) − M X i=1 ckixiφi(v), k = 1, . . . , m. (2)
Here x = (x1, x2, . . . , xM) is the vector of species abundances and v = (v1, . . . , vm) is a vector of resource
66
amounts, where vk is the resource of k-th type consumed by all ecosystem species, µi are the species
67
mortalities, Dk > 0 are resource turnover rates, Sk is the supply of the resource vk, and cik > 0 is the
68
content of k-th resource in the i-th species. The coefficients γi> 0 describe self-limitation effects [15], [11],
69
[10].
70
We consider general φj which are bounded, non-negative and Lipshitz continuous
|φj(v) − φj(˜v)| ≤ Ljkv − ˜vk (3)
and
φk(v) = 0, for all k and v ∈ ∂Rm+. (4)
We use the norm notation kxk = max1≤i≤m|xi|.
71
Furthermore, we shall assume that each φk(v) is a non-decreasing function of each variable vj in RM+.
72
This assumption means that as the amount of j-th resource increases all the functions φl also increase.
73
Conditions (4) and (3) can be interpreted as a generalization of the well known von Liebig law, where
φi(v) = rimin n v1 Ki1+ v1 , . . . , vm Kim+ vm o (5)
where ri and Kij are positive coefficients, i = 1, . . . , M . Here, ri is the maximal level of the resource
74
consumption rate by i-th species and Kij is the half-saturation constant for resource j of species i.
75
The Liebig law can be considered as a generalization of Holling functional response (Michaelis-Menten kinetics) for the case of many resources. It assumes that the species growth is determined by the scarcest resource (limiting factor). In particular, the Liebig law can be applied to ecosystem models for resources such as sunlight or mineral nutrients, for example, for plant ecosystems. For the case of a single resource m = 1 and v = v1∈ R it reduces to the Holling response. In this case, a typical example of φi satisfying all
above conditions is
φi(v) =
riv
Ki+ v
, i = 1, . . . , M. (6)
For γi = 0 system (1), (2) was considered in the studies of the plankton paradox, see, for example, [7].
Following [15] and [10] we assume γi> 0 since it is known that self-limitation is essential for large ecosystem
[2,1] and plankton or plant ecosystems can induce effects leading to self-limitation. We complement system (1), (2) by non-negative initial conditions
x(0) = ¯x, v(0) = ¯v, (7)
where
¯
3. Estimates and equilibria
76
In this section, we study the stability and large time behavior of solutions to the Cauchy problem (1),
77
(2) and (7). Using the standard partial order relations in Rm, we write v ≤ w if v
i ≤ wi for each i, and
78
v w if vi< wi for each i. We also denote z+= max{z, 0}.
79
3.1. Boundedness of solutions
80
Let S = (S1, . . . , Sm). The proof of the following estimates can be found in [10].
81
Proposition 3.1. Solution (x, v) of (1), (2) with initial data (8) is well defined for all positive t, and it satisfies the estimates
0 ≤ xi(t) ≤ ¯ xiexp(¯ait) 1 + ¯xiγi¯a−1i (exp(¯ait) − 1) , i = 1, . . . , M, (9) where ¯ai= φi(S) − µi, and 0 ≤ vk(t) ≤ Sk(1 − exp(−Dkt)) + ¯vkexp(−Dkt), k = 1, . . . , m. (10)
Furthermore, if φk(S) ≤ µk for some k then limt→∞xk(t) = 0.
82
In what follows, we make the following natural assumption:
φi(S) − µi> 0, for all i = 1, . . . , M . (11)
Indeed, if φi(S) − µi ≤ 0 for a certain index i then according to Proposition3.1xi(t) → 0 as t → ∞, thus,
83
species i will surely go extinct and therefore can be excluded from the analysis.
84
3.2. Equilibrium resource values and convergence to equilibria
85
Let E denote the set of nonnegative equilibrium points (stationary solutions) (x, v) of (1)-(2). It is straightforward to see that (0, S) ∈ E. This point expresses the equilibrium resource availabilities in the absence of any species. Furthermore, it was shown in [10] that under assumption (11), for an arbitrary (x, v) ∈ E such that (x, v) 6= (0, S) there holds
x > 0 and 0 v S. (12)
Among all equilibrium points in E, we shall distinguish the special ones defined as follows: an equilibrium point (xeq, veq) is called special if veq is a solution of the fixed point problem
Dk(Sk− veqk ) = Fk(veq), 1 ≤ k ≤ m, (13) where Fk(v) = M X i=1 ckiφi(v) γi (φi(v) − µi)+,
and xeq is uniquely determined by
xeqi = 1 γi
(φi(veq) − µi)+. (14)
Note that Fkcan be interpreted as total consuming ratesγ1ickiφi(v)(φi(v)−µi)+over all the species ≤ i ≤ M . 86
As a corollary of the monotonicity of φi, it can be shown that the set of special equilibrium points (xeq, veq)
87
is nonempty for any choice of the fundamental parameters of the model, see [10] for the proof. If m = 1 then
88
one can easily prove that there always exists a unique special point, see [11]. But, if m ≥ 2 and Di or γi
89
are small enough, there can exist several special points. On the other hand, under some natural conditions,
90
there exists a unique special equilibrium point. More precisely, we have the following result.
Theorem 3.2 (Corollary 7.2 in [10]). Assume that ρ := max 1≤i≤m M X j=1 cijLj(2φj(S) − µi) Diγj ≤ 1 (15)
Then all the solutions (x(t), v(t)) of (1), (2) with positive initial data v(0) > 0, x(0) > 0 converge, as t → ∞,
92
to the unique special equilibrium point defined by (13) and (14).
93
In Section 3.2we will see that the unique point (xeq, veq) is the right equilibrium point for the system
94
(1)-(2), which attracts all trajectories. The proof of Theorem 3.2 is given in our recent paper [10]. Our
95
approach is based on an iteration technique and two side estimates for period-two-points of nonincreasing
96
maps and relies on the monotone properties of our system. More precisely, eliminating in an appropriate
97
way the x-variables, we reduce the analysis to certain integral equations for v = (v1, . . . , vm). It turns
98
out that these integral equations give rise to a nonincreasing operator in an appropriate Banach space of
99
functions 0 ≤ v ≤ S. By the monotonicity, the solution v is estimated by the iterated sequences ¯v(2n) and
100
¯
v(2n+1), which are shown to converge, say to limits v
even and vodd, respectively. The second part of the
101
proof establishes the estimates on the difference veven− vodd in terms of ρ in (15) above. This implies that
102
if the turnover rates Di are large enough then the limits veven and vodd coincide, which yields the required
103
stability result.
104
Alternatively, the main idea of the proof can be explained as follows. If the turnover rate D = mini{Di}
105
are large enough then our system has the standard structure typical for systems with slow/fast variables.
106
More precisely, in that case the resources vi evolve fast in time whereas the species abundances are slow
107
variables. Using this property we make substitution vi= Si− ˜vi, where ˜vi are new unknowns. This yields
108
an a priori estimate |˜vi| = O(D−1). Therefore, given xi(t) we can solve the differential equations for vi
109
by iterations. The a priori estimates also imply that the first approximation for xi can be determined by
110
certain simplified equations which do not involve v. Proceeding further, we get an approximation for vi of
111
order O(D−2), etc. This yields the convergence of the procedure by standard methods of theory of invariant
112
manifolds for slow/fast systems.
113
4. Model with extinction threshold
114
4.1. Formulation of the model
115
System (1)-(2) does not take into account species extinctions due to extinction thresholds. To describe
116
this effect, we follow [11] with certain simplifications. In order to describe species extinction we introduce
117
small positive parameters Xext(i) being an extinction threshold and we will consider only solutions such that
118
xi(0) > X (i)
ext. Estimates of X (i)
ext can be obtained by stochastic models, for example, [13], however, these
119
models are rather complicated. The approach suggested below simplifies the problem and still allows us to
120
obtain nontrivial effects.
121
Let Se(t) denote the set of indices of those species which exist at the time t and let Ne(t) denote the
122
cardinality of Se(t). Then Se(0) = {1, 2, . . . , M } and Ne(0) = M .
123
We say that xk disappears at the moment t∗ if
xk(t∗) = X (k)
ext and xk(t) > X (k)
ext for t < t∗.
In this case we remove the index k from the set Se(t), t ≥ t∗. The parameters Xext(k) can be interpreted as a
124
thresholds for species abundances.
125
With modifications described above, equations (1) and (2) define the dynamics as follows. Within each time interval (t∗, T∗) between the subsequent species extinctions, the dynamical evolution of xi(t) is
determined by the system (1), (2). We obviously have
therefore Ne(t) is a piecewise constant decreasing function, hence there exist limits
Se(t) → Se, N (t) → Ne, as t → +∞, (16)
where Ne is the cardinality of the set Se. Note that the set Se and its cardinality Nedepend on the initial
data as it is shown in [11]. By (16) there exists a time moment Te such that all extinctions have occurred
and thus we can use Theorem3.2for the remained species. For t > Te system (1), (2) can be rewritten as
follows: dxi dt = xi(φi(v) − µi− γixi), i ∈ Se, (17) dvk dt = Dk(Sk− vk) − X i∈Se ckixiφi(v), k = 1, . . . , m, (18) where xi(t) > X (i) ext, ∀t > Te. (19)
This system describes large time dynamics of ecosystem when all extinctions are finished. Note that the set
126
of remaining species Se may be different for different initial data of original system (1), (2).
127
First let us consider the particular case of (17)-(18) is when all species disappear, i.e. Se = ∅. Then
(17)-(18) amounts to
dvk
dt = Dk(Sk− vk) k = 1, . . . , m,
hence the solution converges to the supply equilibrium state (x, v) = (0, S). The dynamic in this case is trivial, therefore, we assume in what follows that
Se6= ∅.
An analysis of the dependence of Seand Neon the initial data seems to be rather complicated, therefore,
we consider instead the maximal possible number Nmax of species which may survive. More precisely, we
define
Nmax= max Ne (20)
where the maximum is taken over all possible sets Se⊂ {1, . . . , M } for which the problem (17)-(19) has a
128
solution. We present some estimates of Nmax in section5 below.
129
4.2. Dynamics of the model with extinction thresholds
130
To obtain equilibria for system (17), (18), we should take into account that xi(t) > X (i)
ext for all t > Tf.
Therefore, equations for equilbria take the form
Xi(veq) = γi−1(φi(veq) − µi)+,i, i= γiX (i) ext. (21) and D(S − veq) = Fext(veq, b, K, p), (22) where Fext(v, b, K, p) = X i∈Se Ri(v, b, K, p). (23)
Here and in what follows, we denote by z+,δ the cut-off function
z+,δ=
0 if z < δ; z if z ≥ δ;
Since the dynamics after the moment Tf is completely determined by differential equations (17), (18),
131
we obtain by Theorem3.2:
132
Corollary 4.1. Suppose that Seis not empty and condition (15) holds. Let (x, v) be solution of (17), (18)
133
and (19). Then this solution converges, as t → ∞, to an equilibrium point satisfied (21), (22), and (23).
5. Estimates of biodiversity
135
In the first three subsections of this section we derive general estimates of biodiversity which make no
136
assumptions on large time behaviour of the system. In section 5.4, we obtain the biodiversity estimates
137
which are asymptotically sharp. In the remaining part, we also establish estimates under assumption that
138
there are not oscillating or chaotic regimes, and each trajectory converge to an equilibria.
139
5.1. The general case
140
In the model with extinctions we are able to estimate the maximal biodiversity, i.e. the maximal possible
141
value Neexpressed by means of the ecosystem parameters. For this purpose we apply the averaging procedure
142
similar to that considered in [11] for m = 1.
143
Let f (t) be a continuous function which is uniformly bounded on [0, +∞). Then its t-average value hf i is defined by hf it= 1 t Z t 0 f (s)ds.
Let x = (x1, . . . , xNe), v = (v1, . . . , m) be a solution of (17)-(19). Applying the t-average to (17) and (18)
and using the boundedness of x and v, we obtain
hxiφi(v)it− µihxiit− γihx2iit= ξi(t) (24) and Dk(Sk− hvkit) = X i∈Se ckihxiφi(v)it+ ˜ξk(t), (25)
where ξi(t) = 1t(xi(t) − xi(0)), ˜ξk(t) = 1t(vk(t) − vk(0)). Since xi and vk are nonnegative and bounded from
144
above, we obtain limt→+∞ξi(t) = limt→+∞ξ˜k(t) = 0.
145
Next, since xi(t) is bounded from above and xi(t) > X (i) ext> 0, we have ηi(t) := x0 i xi t =1 t Z t 0 dxi(s) xi(s) =1 t ln xi(t) xi(0) → 0 as t → +∞,
therefore dividing the left-hand and right-hand sides of (17) by xi(t) followed by the t-average we get
hφi(v)it= µi+ γihxiit+ ηi(t), (26)
where ηi(t) → 0 as t → +∞. Furthermore, by the Cauchy inequality we have hxii2t ≤ hx2iit, hence we derive
from (24) that
hxiφiit≥ (µi+ γiX (i)
ext)hxiit+ ξi(t). (27)
Relations (26) and (27) allow us to obtain a general estimate of consumed resources vk expressed only
146
in the fundamental parameters of the main system (not involving γi and X (i) ext).
147
Lemma 5.1. If the von Liebig law (5) holds, we have for sufficiently large t and all ≤ k ≤ m
hvkit> Vk(Se) := max i∈Se
µiKik
ri− µi
. (28)
Proof. Let T > 0 be chosen such that |ηi(t)| < γiX (i)
ext for all t > T . Next, note that the function φi(v)
defined by (5) is concave as the minimum of concave functions. Therefore, combining Jensen’s inequality with (26) we obtain
φi(hvit) ≥ hφi(v)it> µi, ∀t > T
which implies the desired estimate (28).
Remark 5.2. Note that in the case m = 1 one has hvit> ¯V = min i∈{1,...,M } λi, (29) where λi= µiKi ri− µi . (30)
The parameters λi represent break-even concentrations and appear in analysis of resource competition
149
ecosystems as important species characteristics (see [5] and the references therein).
150
Now, combining (28) with (27) it follows from (25) that for large t
Dk(Sk− Vk) ≥ X i∈Se cki(µi+ γiX (i) ext)hxiit+ ξk(t), (31)
where ξk(t) → 0 as t → +∞. Let us introduce the following averages:
θk(Se) := 1 Ne lim sup t→∞ X i∈Se cki(µi+ γiX (i) ext)hxiit,
and also define
¯ θk = min Se6=∅ θk(Se), V¯k = min Se6=∅ Vk(Se),
where we take the minima over all possible sets Se6= ∅. Then (31) implies
151
Proposition 5.3. The number of survived species satisfies
Ne≤ min k Dk(Sk− ¯Vk) ¯ θk .
Note that the latter estimate involves only the main observable fundamental ecosystem parameters.
152
Furthermore, the obtained estimate is universal in the sense that it is valid without any assumptions on Dk,
153
the number of resources and the ecosystem dynamics.
154
5.2. Modifications of the main estimate
155
Estimate (25) can be simplified for small γiX (i)
ext µi. In that case, one has
Ne≤ min k Dk(Sk− ¯Vk) ˆ θk , where ˆ θk = min Se6=∅ 1 Ne lim sup t→∞ X i∈Se ckiµihxiit.
Taking into account that xi(t) > X (i)
ext, we also have a rough estimate
Ne≤ min k Dk(Sk− ¯Vk) Zk , (32) where Zk = min Se6=∅ Ne−1 X i∈Se cikµiX (i) ext.
5.3. Estimates of biodiversity from below
156
Now we want to estimate the biodiversity from below. We consider the case of a single resource m = 1. To this end, let us introduce the set
BM = n i ∈ {1, . . . , M } : riV¯ Ki+ S > µi for all i o ,
where ¯V is defined by (29). The cardinality of a set A will be denoted by |A|.
157
Proposition 5.4.
158
a) Consider the dynamic without extinctions defined by (1) and (2), where φi are defined by von Liebig’s
159
law (5). Then the number of species such that lim inft→+∞xi(t) > 0 is not less than |BM|.
160
b) In the model with extinctions, there holds Ne= 0 for some initial data (even if |BM| > 0).
161
So, this claim shows in particular that models with and without extinctions have completely different
162
behaviour. Remarkably, the obtained biodiversity estimate does not involve γi, while the proof makes use
163
the fact that γi> 0.
164
Proof. Consider a). We use (26) that gives
γihxiit= hφi(v)it− µi+ ξi(t) (33)
where ξi(t) → 0 as t → ∞. Since Kriv
i+v >
ri
Ki+Sv, we find from (33) that
γihxiit≥
ri
Ki+ S
hvit− µi+ ξi(t).
Now we use (29) and then a) is obtained from the last inequality as t → ∞.
165
Let us consider b). Let v(0) be small enough such that
φi(v(0)) − µi= −2κi < 0
holds for all i = 1, . . . , M . Note that vl(t) < DlSlt + vl(0). Then there exists τ > 0 such that
φi(v(t)) − µi< −κi< 0 t ∈ (0, τ ).
Therefore,
xi(t) < xi(0) exp(−κit), t ∈ (0, τ )
Therefore, if all xi(0) are sufficiently close to X (i)
ext all the species extinct.
166
To illustrate the dependence of the biodiversity the fundamental parameters, let us consider an example (see also Example 2 below). Let Ki = K for all i = 1, . . . , M and ai = µi/ri. Then λi = K(1 − ai)−1.
Therefore, according to Proposition5.4all species xj with
aj 1 +
S
K <i=1,...,Mmin (1 − ai)
−1 (34)
survive. The set of such species can have the maximal cardinality M . Indeed, let ai = 12 + ˜ai, where
0 < ˜ai< 12 and let us assume that S/K < 1. Then the right hand side of (34) will be less than 2. Therefore,
the species with the property
˜ ai< 2 1 + S K −1 −1 2
do certainly survive. Thus, if all ˜ai are small enough, all the species survive.
5.4. Three examples
168
We illustrate the obtained results by simple examples relating them to the (unified) neutral theory.
169
Recall that neutrality means that at a given trophic level in a food web, species are equivalent in birth rates,
170
death rates, dispersal rates and speciation rates, when measured on a per-capita basis. Mathematically this
171
means that all basic parameters are almost equal, see [6]. In the all example we suppose that Dk are large.
172
Example 1. For large d = minkDk system (13) becomes a fixed point problem with a contraction operator,
therefore its solution veqis uniquely determined and vkeq= Sk− D−1k
M
X
i=1
ckiφi(S)γi−1(φi(S) − µi)++ O(d−2). (35)
By Theorem3.2all solutions (x, v) of the system (1), (2) with the Cauchy data (8) have a limit vk(t) → Sk+ O(d−1), k = 1, . . . , m,
xi(t) →
φi(S) − µi
γi
+ O(d−1).
So, if d is sufficiently large the system is persistent. Moreover, the number Nmax is equal approximately to
the cardinality {i : φi(S) − µi γi > Xext(i)} .
Example 2 (The case of a single resource). In the case m = 1 we set φi =Kriv
i+v and suppose that
µi= µ, Ki = K, ri= r, γi= γ, C = N¯ e−1 Ne
X
i=1
ci
and Xext(i) = Xext. Let us introduce parameters
p = µ r, ˜ S = S K, = γXext r , R = KDγ r2C¯1. (36)
We consider veq= Ku, where u is to be determined. Then (13) becomes
R( ˜S − u) = NeF (u), (37) where F (u) = u 1 + u u 1 + u− p +, . Let u= p + 1 − p − . (38)
We seek a nontrivial solution of (37) such that u > u. Equation (37) implies
N∗(veq) ≤ Ne≤ N∗(veq) + 1, (39) where N∗(v) = " R(S − v)(K + v) Kv(K+vv − p)+, # (40)
and [x] is the floor of x. Note that the function N∗(veq) is decreasing in veq. Then since u > u , relations
(38) and (39) allows us to conclude that the maximal possible biodiversity Nmax satisfies
N∗(Ku) < Nmax< N∗(Ku) + 1. (41)
In the Section6 we shall see that this maximum of biodiversity is realizable for some initial data.
Example 3 (The multi-resource case). We assume that m ≥ 2 and that the ecosystem is in an equilibrium state defined by (21), (13) and (23). We also assume that we are in the neutral position, i.e.
Kki= K, ri= r, µi= µ, γi= γ, Xext= X (i)
ext. (42)
Then setting φ(z) := K+zrz we obtain
φi(v) = min
k φ(vk) = φ(w), (43)
where w := minkvk = vk∗. Then the equilibrium abundances and resources are determined respectively by
xi(w) = γ−1(φ(w) − µ)+,γXext, 1 ≤ i ≤ Ne, (44) vkeq(w) = γ−1Dk Ne X i=1 ckiφ(w)(φ(w) − µ)+,γXext, 1 ≤ k ≤ m. (45) Since w = vkeq∗(w) = vk∗, we obtain w = γ−1Dk∗ Ne X i=1 ck∗iφ(w)(φ(w) − µ)+,γXext. (46)
The system (44), (45) and (46) determines an equilibrium state depending on the index k∗∈ {1, 2, . . . , m}.
174
For any fixed k∗ = 1, 2, . . . , m, we solve (44), (45) and (46). If vkeq < w is valid for all k 6= k∗ then k∗ is
175
found.
176
Using change of the variables w = Ku, where u is a new variable, we obtain from (46) a similar relation (37) as in Example 2, where all parameters but R are defined by (36), and R = KDk∗γr
−2C¯−1 m , where ¯ Cm= M−1 M X i=1 ck∗,i. Then N∗(w) ≤ Ne≤ N∗(w) + 1, (47) where N (v) is defined by (40). 177
6. Asymptotically sharp estimate of biodiversity for systems with random parameters
178
In the remaining part of the paper we consider ecosystems with random parameters for large values of
179
M . Moreover, we suppose that the system is in an equilibiria state, i.e. oscillating large time regimes are
180
absent.
181
First let us suppose that at the initial moment there are M 1 species with different (possibly random)
182
parameters and that the resource supplies Sk are limited. We are going to address the following problem:
183
How many species Newill survive?
184
It is clear that Ne is a priori bounded by M . Below we obtain an asymptotically sharp estimate of Ne.
We shall assume that m = 1 and the ecosystem has a globally convergent equilibrium state, v = veq. We
also assume that φi are defined by the von Liebig law (5). Then, if the i-th species is survived for all times
t > 0, we have riveq Ki+ veq − µi> γiX (i) ext. (48) Consequently, veq> βi, βi= (µi+ γiX (i) ext)Ki ri− µi− γiX (i) ext , (49)
and
ri> µi+ γiX (i)
ext. (50)
Note also that for γi = 0 or X (i)
ext = 0 the numbers βi coincide with the parameters λi defined by (30).
185
They determine survival of species in classical resource competition models without extinction thresholds
186
and without self-limitation [5].
187 We introduce parameters ηi= D−1ci(γiX (i) ext+ µi)X (i) ext, i = 1, . . . , M (51)
and describe the choice of random values of the model parameters. Under assumptions (5), the main
188
parameters of our model are the coefficients µi, γiand the vectors K(i)= (Ki1, . . . , Kim), where i = 1, . . . , M ,
189
and r = (r1, . . . , rM). As before, we consider UM = µ, K, r, c, Γ ∈ R2M (m+2). Note that ci are species
190
specific parameters which not necessary to be included in analysis, so we suppose that ci are fixed.
191
Let UM be a random vector with a probability density function ξ(UM). This means that the values UM
192
are defined by random sampling, i.e. the parameters of the species are random independent vectors that are
193
drawn from the cone R2M (m+2)+ by the density ξ.
194
Our basic assumption to ξ then can be formulated as follows:
195
Assumption 6.1. The probability density function ξ is a continuous function with a compact support in Rn+. Moreover, as N → ∞ the parameters βi are distributed on (0, +∞) with the smooth probability density
ρ0(β) such that
supp ρ0(β) ⊂ (βmin, βmax). (52)
The parameters ηi are distributed on (0, +∞) with a continuous probability density ρ1(η). The densities ρ0
196
and ρ1 are mutually independent.
197
From Assumption 6.1 it follows that the function ξ is positive on Sξ, where Sξ is an open bounded
set. Assumption6.1also yields that the mortality rates do not approach zero and resource consumption is restricted. It is supposed that initial data ¯xi = xi(0) are random mutually independent with the density
distribution
¯
xi∈ X ( ¯X, σX)
with the mean ¯X and the deviation σX. The random assembly of the species defines an initial state of the
198
ecosystem for t = 0.
199
Let us order βi so that
βmin≤ β1≤ β2≤ . . . ≤ βM ≤ βmax. Then βNe < v eq< β Ne+1. Therefore, we obtain βNe+1> S − D −1 Ne X i=1 ci riveq Ki+ veq ( riv eq Ki+ veq − µi)γi−1, (53) βNe < S − D −1 Ne X i=1 ci riveq Ki+ veq ( riv eq Ki+ veq − µi)γi−1, (54)
where by Assumption6.1one has that δβ= βNe+1− βNe= O(M
−1) as M → +∞. Therefore, the last two
200
inequalities are asymptotically exact.
201
Lemma 6.2. Let us define the probability Prβ by
Prβ= Prob{βn− βmin< M−1/3}
Then for sufficiently large M
Prβ > 1 − Cnexp(−cnM−1/6),
where Cn, cn are positive constants.
Proof. Let us consider the probability Prk,δ that a random sample consisting of M numbers βi containing
exactly k numbers within the interval J = [βmin, βmin + δ], where δ is a small positive number. The
probability that a random βi distrubuted according to the density ρ0 lies within J is
pδ =
Z βmin+δ
βmin
ρ0(s)ds.
Since ρ0is smooth, the probability can be estimated as pδ < cδ2. Let us define
Prk,δ= M k pkδ(1 − pδ)M −k. Then Prk,δ< Ck(M pδ)kexp(−0.5M pδ),
where the constants Ck> 0 do not depend on M . We set δ = M−1/3. Then
Prk,δ< CkM−k/6exp(−c(k)M−1/6). (55) We observe that Prβ= 1 − n X k=1 Prk,δ> 1 − n max k∈{1,....,n}Prk,δ. (56)
For fixed n and M large enough one has
max
k∈{1,....,n}Prk,δ= Prn,δ.
Combining this with (55) and (56) yields that for sufficiently large M
Prβ> 1 − nPrn,δ > 1 − Cn(M )−n/6exp(−cnM−1/6).
The latter estimate proves (6.2).
203
Proposition 5.3shows that Ne M and ρ0 is a smooth density. Therefore, by Lemma6.2we conclude
that βNe− βmin< M −1/3 and ( riveq Ki+ veq − µi) = X (i) ext+ O(M −1/3),
with a probability exponentially close to 1 for large M . This yields the following relation which is asymp-totically sharp: βmin+ Ne X i=1 θi= S + O(M−1/3).
Since the densities ρ0 and ρ1 are mutually independent, the sum in the left hand side of the last equation
can be replaced by Nehθi, where hθi = Eθ is the averaged value of θ. Finally one has
Ne= (S − βmin)hθi−1+ O(M−1/3). (57)
The latter estimate considerably refines our previous result [11]. Moreover, we conclude that the final state of ecosystem originally consisting of many species with random parameters can be described by the so-called R∗rule. The essential parameters βi are well localized in a narrow domain that follows from (53) and
(54) (see Fig. 1) that confirms the R∗rule. Recall that the R∗rule (also called the resource-ratio hypothesis) is a hypothesis in community ecology that attempts to predict which species will become dominant as the result of competition for resources. It predicts that if multiple species are competing for a single limiting
-2 -1 0 1 2 3 4 5
logarithm of consuming rate r
-4 -3 -2 -1 0 1 2 3 4 logarithm of mortality
Figure 1: This plot shows the case when an ecosystem in a neutral state where all parameters are localized. We assume that the population parameters are subject to log-normal laws. The corresponding averages E ln r = 1 and E ln µ = −0.5. These parameters are shown by the star point in the center of a circle of the small radius 0.5. All the species parameters are localized inside this circle. The half-plane ΩS consists of the points that lie below the upper right straight line LS defined by eq.
y = x − S. As it follows from our analytical arguments, if species parameters lie in ΩSthen that species survives; otherwise it
would be extinct. For S = 1 all species survive. If resource supply S slightly diminishes, the line LS shifts to the bottom and
only the half species rest inside the small circle the second half of the species disappears.
resource, then species, which survive at the lowest equilibrium resource level, outcompete all other species [20]. Note that a generalized resource ratio value βi for small γiX
(i) ext is
βi≈ λi,
where λi are defined by (30), and it involves the mortality, sharpness of consuming rates, and consuming
204
rates. Note, however, that the proof of (57) cannot be extended for m > 1.
205
To conclude this section, let us remark that (57) improves (32). Indeed, Ne from (57) satisfies (32).
206
Moreover, a difference between those estimate is that in (57) the nominator S is replaced by S − βmin,
207
therefore, for large resource supply these estimates coincide.
208
7. Extinctions and mass extinctions
209
7.1. Mass extinctions
210
Let us revisit our Examples1 and 2 above, which describe the localized case (Assumption 6.1), where parameters of all species are localized in a narrow domain that consistent with the R∗-rule. Note that the nontrivial solution u 6= 0 of eq. (37) exists under condition u< ˜S thus
µ + γXext
r − µ − γXext
< S/K.
The violation of this condition leads to a mass extinction, when all the species disappear. In [11], the
211
following localization effect is described: in a long time evolving species population, which was randomly
212
assembled at the initial moment (see Assumption 6.1 above), all species parameters converge to a limit
213
value. Then one can show that this localized population state has an interesting property:
214
All species die together: As a result of a long evolution, initially randomly assembled population
215
approaches to a neutral fitness invariant state, where, if an extinction occurs then the extinction is a mass
216
extinction, i.e. all species die simultaneously.
217
So, we conclude that
1. extinctions are stronger if the variations of species parameters are smaller, i.e. localization of
param-219
eters is higher;
220
2. for well localized parameters (σµ and σrare small) there is a critical domain of parameter values when
221
a small variation in the resource supply Sk can lead to a mass extinction.
222
These properties can be illustrated by Fig. 1.
223
We assume m = 1 and that coefficients ci, ri, γi and µi satisfy
C−r <ri< C+r, C−c < ci < C+c, 1 ≤ i ≤ M, (58)
C−γ <γi< C+γ, C−µ < µi< C+µ, 1 ≤ i ≤ M (59)
where a, c, γ, r are characteristic values of the corresponding coefficients, C±are positive constants
indepen-224
dent of M, r, c, γ, µ. Let us introduce the stress parameter by
225
Pstress=
cr2
γDS.
Let us find a relation between that parameter and the biodiversity robustness, which we define as a coefficient Rb in the relation
∆Ne
Ne
= R−1b ∆S
S , (60)
where ∆S < 0 is a small variation of the resource supply S caused for instance by a climate change, and ∆Ne ≤ 0 is the corresponding variation of the number of coexisting species Ne. We suppose that relation
(57) is fulfilled. Then comparing (57) and (60) we see that
hθi = M−1
M
X
i=1
θi,
where θi are defined by (51). By (51) we note that in the case of large biodiversity Ne 1 the coefficient
226
hθi of same order that SPstress. In fact, if all r,γi and µi of the same order (as it was supposed above, in
227
(58) and (59)), then Xext has order r/γ. Substituting this relation in (51), we obtain that Rbhas the same
228
order SPstress. Note that Pstress is a dimensionless parameter.
229
8. Conclusions
230
The human activity affects ecosystems restricting their resources and producing climatic changes. The
231
key question is about consequences of that activity, it will be an abrupt change of biodiversity comparable
232
with famous mass extinctions such as the end-Permian extinction (252 million years ago (Ma)), and even
233
more severe the Phanerozoic (the past 542 Ma), or we will observe a relatively smooth decline of biodiversity.
234
The history of the Earth system show that some changes were gradual, but others were sharp and
235
produced catastrophic mass extinctions. To understand that wait us in Future we should investigate
ecosys-236
tem biodiversity, find key factors affecting this biodiversity and estimate the sensitivity of ecosystems with
237
respect to environmental parameters.
238
A number of works are devoted to biodiversity problem beginning with [21] (see, for example, [4,1,9,20,
239
7,8] this list does not pretend to be complete). However, mathematical methods developed up to now, does
240
not allow us to obtain explicit estimate of ecosystem biodiversity via the fundamental and experimentally
241
measurable ecosystem parameters. The main difficulty is competitive exclusion principle. According to
242
[5] for the case of one fixed resource supply we obtained two situations. The species survival depends on
243
”break-even” parameters λi (see also (30) below). If λ1is less than all the rest λi only the species number
244
1 survives. To obtain coexistence of m species observed in many natural resource competition ecosystems
245
([9]) we must set λ1= λ2 = . . . = λm, where we, however, can take m in an arbitrary way (only to be less
246
than general species number). To get around this difficulty a number of different approach were suggested,
247
for example, dynamical chaos in systems with more than 2 resources [7, 8], to take into account temporal
variations of parameters etc. (see for example [17] for an overview) but these ideas does not permit us to
249
obtain an explicit analytic estimates for the number of coexisting species.
250
In this paper, we develop further two relatively new ideas: namely we take into account weak
self-251
limitation effects [15] and extinction thresholds [11]. Although the modified models seem to be more
sophis-252
ticated, these effects allow us to regularize the problem mathematically even in the case of a single resource
253
and also to obtain explicit biodiversity estimates for diverse scenarios.
254
We also extend results [14] and show that for plant or plankton systems, the magnitude of the future
255
catastrophe depends not only on the level of climate variations. The structure of ecosystems also influences
256
the number of species that go to extinction. We establish estimates of biodiversity via the system parameters
257
such as the mortality rate, the consuming rate, teh sharpness of consuming rate, self-limitation and the
258
critical species abundance. Furthermore, we show that the ecosystems where these parameters are well
259
localized are less stable that the ecosystems with large variations in these parameters. We also investigate
260
how biodiversity depends on structure climate-ecosystem interactions.
261
In summary, or main results are:
262
1. The dynamics of model without extinctions is determined for large resource turnovers: we show that
263
all trajectories converge to a unique equilibrium, system has no memory and forgets initial data
264
completely.
265
2. For large resource turnovers, the dynamics of the model with extinctions is stable in a weaker sense:
266
any trajectory goes to an equilibria which now may depend on the initial data. Thus, the ecosystem
267
has memory and the final equilibria depend on initial data.
268
3. For models with extinctions, the ecosystem biodiversity can be explicitly estimated by a relation
269
involving only measurable parameters (turnover and mortality rates, resource supplies and the species
270
content coefficients). This estimate is asymptotically sharp under certain natural assumptions and in
271
the case when the neutral theory correctly describes the ecosystem state. We show also that such a
272
situation naturally occurs when the ecosystem is under stress, i.e. originally there exists species with
273
randomly distributed parameters and the most of them disappear as a result of resource limitations.
274
4. In the stress situation, we also show that the final state of the system (after extinctions) can be
275
determined according to the R∗ rule.
276
5. The ecosystem robustness is investigated and conditions for mass extinctions are obtained. The
robust-277
ness depends, in particular, on the dimensionless stress parameter introduced in [11]. Systems with a
278
large stress parameter is less sensitive with respect to resource decrease, however, if an extinction is
279
happened. Such an event is catastrophic and destroys an essential part of that ecosystem.
280
Acknowledgments
281
The authors would like to thank the reviewers for their detailed comments and suggestions for the
282
manuscript. We are also grateful to L. Mander for a question on extinction problem that have stimulated
283
the present research. The second author was supported by Link¨oping University, by Government of Russian
284
Federation, Grant 08-08 and via the RFBR grant 16-01-00648.
285
References
286
[1] S. Allesina. Ecology: the more the merrier. Nature, 487(7406):175–176, 2012. 287
[2] S. Allesina and Si Tang. Stability criteria for complex ecosystems. Nature, 483(7388):205–208, 2012. 288
[3] R.A. Fisher, J.H. Bennett, and H. Bennett. The Genetical Theory of Natural Selection: A Complete Variorum Edition. 289
OUP Oxford, 1999. 290
[4] Garrett Hardin. The competitive exclusion principle. Science, 131(3409):1292–1297, 1960. 291
[5] S.B. Hsu. A survey of constructing Lyapunov functions for mathematical models in population biology. Taiwanese J. 292
Math., 9(2):151–173, 2005. 293
[6] S.P. Hubbell. The Unified Neutral Theory of Biodiversity and Biogeography (MPB-32). Princeton University Press, 2001. 294
[7] J. Huisman and F.J. Weissing. Biodiversity of plankton by species oscillations and chaos. Nature, 402(6760):407–410, 295
1999. 296
[8] J. Huisman and F.J. Weissing. Biological conditions for oscillations and chaos generated by multispecies competition. 297
Ecology, 82(10), 2001. 298
[9] G. E. Hutchinson. The paradox of the plankton. Amer. Natur., 95(882):137–145, 1961. 299
[10] V. Kozlov, V. Tkachev, S. Vakulenko, and U. Wennergren. Global stability and persistence of complex foodwebs. submitted 300
to Bull. of Math. Biol., 2017. 301
[11] V. Kozlov, S. Vakulenko, and U. Wennergren. Stability of ecosystems under invasions. Bull. Math. Biol., 78(11):2186–2211, 302
2016. 303
[12] V. Kozlov, S. Vakulenko, and U. Wennergren. Biodiversity, extinctions, and evolution of ecosystems with shared resources. 304
Phys. Rev. E, 95:032413, Mar 2017. 305
[13] Ingemar N˚asell. Extinction and quasi-stationarity in the verhulst logistic model. J. Theor. Biology, 211(1):11 – 27, 2001. 306
[14] Daniel H. Rothman. Thresholds of catastrophe in the earth system. Science Advances, 3(9), 2017. 307
[15] S. Roy and J. Chattopadhyay. Towards a resolution of the paradox of the plankton: A brief overview of the proposed 308
mechanisms. Ecological complexity, 4(1):26–33, 2007. 309
[16] H. L. Smith and P. Waltman. The theory of the chemostat, volume 13 of Cambridge Studies in Mathematical Biology. 310
Cambridge University Press, Cambridge, 1995. Dynamics of microbial competition. 311
[17] U. Sommer and B. Worm, editors. Competition and Coexistence, volume 161 of Ecological Studies. Springer-Verlag Berlin 312
Heidelberg, 2002. 313
[18] Ivan Sudakov, Sergey A. Vakulenko, Dubrava Kirievskaya, and Kenneth M. Golden. Large ecosystems in transition: 314
Bifurcations and mass extinction. Ecological Complexity, 32:209 – 216, 2017. Uncertainty in Ecology. 315
[19] D. Tilman. Resources: A graphical-mechanistic approach to competition and predation. Amer. Natur., 116(3):362–393, 316
1980. 317
[20] David Tilman. Resource competition and community structure. Princeton University Press, 1982. 318
[21] Vito Volterra. Le¸cons sur la th´eorie math´ematique de la lutte pour la vie. Les Grands Classiques Gauthier-Villars. 319
[Gauthier-Villars Great Classics]. ´Editions Jacques Gabay, Sceaux, 1990. Reprint of the 1931 original. 320
