Biodiversity, extinctions and evolution of ecosystems with shared resources

TL;DR

This paper analyzes how stable ecological networks with shared resources form through extinctions, deriving relations for species coexistence, extinction impacts, and trait concentration, with a focus on a stress parameter influencing biodiversity.

Contribution

It introduces a novel analytical framework linking species traits, resource characteristics, and environmental stress to ecosystem stability and biodiversity outcomes.

Findings

Stable ecosystems emerge via extinctions driven by environmental stress.

Derived analytical relations for species coexistence and extinction impacts.

Identified a key parameter predicting trait concentration and biodiversity levels.

Abstract

We investigate the formation of stable ecological networks where many species share the same resource. We show that such stable ecosystem naturally occurs as a result of extinctions. We obtain an analytical relation for the number of coexisting species and find a relation describing how many species that may go extinct as a result of a sharp environmental change. We introduce a special parameter that is a combination of species traits and resource characteristics used in the model formulation. This parameter describes the pressure on system to converge, by extinctions. When that stress parameter is large we obtain that the species traits concentrate at some values. This stress parameter is thereby a parameter that determines the level of final biodiversity of the system. Moreover, we show that dynamics of this limit system can be described by simple differential equations.

Equations

$$\frac{d{x}_i}{dt}=x_i(-r_i+{\varphi }_i(v)-{\gamma }_i{x}_i),i=1,\dots ,M$$

$$\frac{dv}{dt}=D(S_0-v)-\sum _{i=1}^M{c}_i{x}_i{\varphi }_i(v),$$

$${\varphi }_i(v)=a_i\Phi (v,K_i),\Phi (v,K)=\frac{v}{K_i+v}.$$

$$x_i(0)={\bar{x}}_i,v(0)=v_0.$$

$$0\le x_i(t)\le \frac{{\bar{x}}_i\exp ({\bar{a}}_{i}t)}{1+{\bar{x}}_i{\gamma }_i{\bar{a}}_i^{-1}(\exp ({\bar{a}}_{i}t)-1)},$$

$$0\le v(t)\le S_0(1-\exp (-Dt))+v(0)\exp (-Dt).$$

$$\frac{d{y}_i}{dt}=y_i(-r_i+a_i-{\gamma }_i{y}_i),y_i(0)={\bar{x}}_i.$$

$$X_i(v)=(\Phi (v,K_i)-p_i{)}_{+},p_i=r_i/a_i,$$

$$F_M(v)=F_M(v,b,K,p)=\sum _{i=1}^M{R}_i(v,b,K,p),$$

$$R_i(v,b,K,p)=b_i\Phi (v,K_i)X_i(v),b_i=c_i{\gamma }_i^{-1}{a}_i^2.$$

$$D(S_0-v)=F_M(v,b,K,p),$$

$$x_i=a_i{\gamma }_i^{-1}{X}_i(v),i=1,\dots ,M$$

$$\underset{i}{\max }(\Phi (S_0,K_i)-p_i)>0.$$

$$v=S_0-\frac{1}{D}\sum _{i=1}^M{b}_i\Phi (v,K_i)X_i(w)=G(v,w),$$

$$v=V(v),v\in [0,S_0].$$

$$v^{(k+1)}=V(v^{(k)}),k=0,1,\dots ,\text{mbox}and{v}^{(0)}=0.$$

$$v^{(2)}=G(v^{(2)},S_0),$$

$$0=v^{(0)}<v^{(2)}\le v^{(4)}\le \cdots \cdots \le v^{(3)}\le v^{(1)}=S_0.$$

$$\hat{v}=\underset{k\to \infty }{\lim }{v}^{(2k)}\text{mbox}and\check{v}=\underset{k\to \infty }{\lim }{v}^{(2k+1)}.$$

$$0<v^{(2)}\le \hat{v}\le v_{eq}\le \check{v}\le S_0.$$

$$\hat{v}=S_0-\frac{1}{D}\sum _{i=1}^M{b}_i\Phi (\hat{v},K_i)X_i(\check{v})$$

$$\underset{t\to \infty }{\mathrm{liminf}}v(t)\ge \check{v},\underset{t\to \infty }{\mathrm{limsup}}v(t)\le \hat{v}$$

$$\underset{t\to \infty }{\mathrm{liminf}}{x}_i(t)\ge X_i(\check{v}),\underset{t\to \infty }{\mathrm{limsup}}{x}_i(t)\le X_i(\hat{v}),$$

$$\lambda +D+G(\lambda )=0,$$

$$G(\lambda)=\sum_{i\in{\bf I}_{eq}}c_{i}\big{(}x_{i}\phi_{i}^{\prime}(v_{eq})+\phi_{i}(v_{eq})\frac{x_{i}\phi_{i}^{\prime}(v_{eq})}{\lambda+P_{i}(v_{eq})}\big{)}$$

$$Re\lambda +D+ReG=0,$$

$$Re\ G=\sum_{i\in{\bf I}_{eq}}c_{i}\big{(}x_{i}\phi_{i}^{\prime}(v_{eq})+\phi_{i}(v_{eq})\frac{x_{i}\phi_{i}^{\prime}(v_{eq})(Re\lambda+P_{i}(v_{eq}))}{|\lambda+P_{i}(v_{eq})|^{2}}\big{)}.$$

$$Re\lambda \le -D-\sum _{i\in {\mathbf{I}}_{eq}}{c}_i{x}_i{\varphi }_i^{\prime }(v_{eq})\text{mbox}orRe\lambda \le -\underset{i\in {\mathbf{I}}_{eq}}{\min }{P}_i(v_{eq}).$$

$${\bar{x}}_i\in \mathcal{X}(\bar{X},{\sigma }_X)$$

$$I_M=S_e(t)\cup S_v(t)t\ge 0.$$

Full text

I Introduction

Thermodynamics and Statistical Physics allow us to describe equilibrium states of systems consisting of many particles by a few variables. This approach is very effective for many physical and chemical applications, and particularly, when we are dealing with closed systems. It would be very tempting to find such reduced macroscopic descriptions for large ecological and economic systems.

Recently, in the paper Barab a very reduced description is proposed for a class of ecological models. It is based on an analog of the mean field theory and aims to describe possible bifurcations. In the present paper, we suggest a variant of a reduced description for other large class of ecosystems, where many species compete for a few resources. Our second goal is to consider one of the most intriguing puzzles in ecosystem theory, namely, the biodiversity problem: why can such a large number of similar species share the same habitat and how to estimate this number. Aquatic ecosystems is a typical example where many different phytoplankton species do coexist. The principle of competitive exclusion Volterra ; Hardin asserts that different species sharing the same resource cannot coexist and it predicts that an assembly of competing species will converge to a single species. In general, competition models show that the number of species that can coexist in equilibrium cannot be greater than the number of limiting factors Til77 ; Zeeman . However, hundreds of species of phytoplankton coexist although only nitrate, phosphate, light and carbon are resources regulating phytoplankton growth Hu61 ; Roy . The reasons of this fact are poorly understood although the problem rests under a great attention of ecologists (see, for example, Moll ; Harb ; Record , among many others) and they have suggested numerous different approaches to that problem based on game theory, chaos, stochastics, space inhomogeneities, turbulence etc.

In this paper, we consider a model of an ecological system where many species share the very same resource. Our dynamical equations are close to eqs. considered in HuWe99 but we extend the model of HuWe99 to take into account species extinctions and self-limitation effects (which are important, in particular, for plankton populations, Roy ). There is a number of works devoted to evolution in a multi species context as food-webs, for example, Loreau ; Loe ; Drossel ; McKane as well as the effect of species extinctions Kondo ; Gall . These papers have used the Lotka Volterra model where species interactions slowly evolve in time (that may be connected with foraging Kondo or an adaptation of species behaviour) and species extinctions are possible when the abundances attain a critical threshold. In our model, all parameters are random but, opposite to Kondo ; Gall , fixed in time and the primary effect on evolution is species extinctions.

The main results are as follows. We introduce a specific numerical characteristics, which we will call the stress parameter $P_{stress}$. That characteristics is a natural dimensionless multiplicative combination of some main parameters involved in the model formulation, namely, the resource turnover rate, the maximum supply of the resource, self-limitation coefficient and averaged specific growth rates. The stress parameter appears, in a natural way, as a result of a model rescaling and it can be interpreted as a magnitude of selection pressure on ecosystem species induced by the interaction between the ecosystem and its environment. For example, $P_{stress}$ is large if the resource amount or the resource turnover rate is small. For large values of the stress parameter, the model exhibits an effect of convergence to similarity as in the competitive Lotka-Volterra systems studied earlier in Sim . In contrast to Sim , we obtain a complete analytical description of the system behaviour. We show that the model with or without extinction thresholds are sharply different. Namely, the simpler model without critical extinction threshold exhibits a global stability, when all positive trajectories converge to the same equilibrium independent of initial state. The model with extinctions dynamics is fundamentally non-predictable. The trajectories tend to different equilibria, and these final states depend on initial data. The biodiversity level, which we observe after a long evolution, can be expressed via the initial number of species, the stress parameter and other ecosystem characteristics in an explicit way. This convergence to similarity in a species trait can be called a ”concentration effect”. However, this concentration effect does not mean that all species are completely identical: other traits defined by species parameters can still differ.

The concentration effect leads to interesting phenomena. Let us suppose that initially an ecosystem contains a number of species with random parameters and consider the limit system, which is a result of a long evolution. We compare how a sharp change of environment can affect initial and final systems. By our analytical relations one can estimate the number of species that may go to extinct. The limit system is more stable than the initial one.

The limit ecosystem, arising as a result of a long evolution, has interesting properties: its dynamics is governed by a simple differential equation of the second order. This equation describes a nonlinear oscillator with a friction and a memory. If the friction is small and the memory is negligible the dynamics of this oscillator is defined by a Hamiltonian system. Formation of this universal limit system is a result of species extinctions and a selection pressure on some species parameters, namely those important for species survival.

II Population dynamics

We consider the following system of equations:

[TABLE]

[TABLE]

where

[TABLE]

Here $x_{i}$ are species abundances, $M$ is the number of species, $v$ is resource amount, $D$ is the resource turnover rate, $S_{0}$ is the maximum supply of resource $v$. In total there are five species specific parameters: $r_{i}$ are species mortality, $c_{i}>0$ is the fraction of resource consumption by individuals of the $i$ species, $\gamma_{i}>0$ defines species specific self-limitation, the coefficients $a_{i}>0$ are species specific growth rates, and $K_{i}>0$ are species specific resource constants indicating reduction of resource effect by half. For $\gamma_{i}=0$ this system have been used to study the plankton paradox HuWe99 . Following Roy we assume $\gamma_{i}>0$ since it is known that self-limitation is essential for large ecosystems Alles1 ; AllPR and that plankton and plant ecosystems can induce effects leading to self-limitation Roy . We complement the system (1), (2) with the initial conditions

[TABLE]

III Global stability for model without extinctions

Here we show that the Cauchy problem (1), (2) and (4) has a positive solution for all positive Cauchy data. Furthermore, we study both the stability and large time behavior of solutions.

Proposition I. Solution $(x(t),v(t))$ of (1), (2) with initial data $v(0)\geq 0,\bar{x}_{i}=x_{i}(0)\geq 0$ is defined for all positive $t$, and it satisfies the estimates

[TABLE]

where $\bar{a}_{i}=a_{i}-r_{i}$, and

[TABLE]

Proof. Since $\phi_{i}(v)<a_{i}$, we have $x_{i}(t)\leq y_{i}(t)$, where $y_{i}(t)$ is the solution of the Cauchy problem

[TABLE]

Solving this equation, we obtain (5). Estimate (6) follows from the non-negativity of the term $\sum_{i=1}^{M}c_{i}\;x_{i}\;\phi_{i}(v)$.

III.1 Global stability

Let

[TABLE]

where $z_{+}=\max\{z,0\}$, and let also

[TABLE]

where

[TABLE]

Here the quantity $R_{i}$ can be interpreted as a consuming rate of $i$th species and $F_{M}$ is the sum of all consuming rates.

Here $a=(a_{1},...,a_{M})$, $b=(b_{1},...,b_{M})$, $p=(p_{1},...,p_{M})$ and $K=(K_{1},...,K_{M})$. If $v$ is a non-negative root of the equation

[TABLE]

then

[TABLE]

and $v$ is an equilibrium point of the system (1), (2). We assume here and in what follows that

[TABLE]

Then the function $F_{M}$ is non-negative for $v\geq 0$, $F_{M}(v)=0$ and $F_{M}(S_{0})>0$ due to (11). This implies that equation (9) has a unique non-negative solution, which belongs to the interval $(0,S_{0})$. We denote this solution by $v_{eq}$.

Let us rewrite equation (9) in the following way. Consider first the relation

[TABLE]

with $w\in[0,S_{0}]$. Since $G(v,w)$ is decreasing in $w$ from $S_{0}$ to something which is smaller than $S_{0}$, for each $w$ the above equation has a unique solution $v=V(w)$. One can verify that the function $V$ is non-decreasing and continuous, $V(0)=S_{0}$ and $V(S_{0})>0$. We can consider (9) as the following fixed-point equation

[TABLE]

In order to describe large time behavior of the system (1), (2), we consider the following iterative procedure of solving (12):

[TABLE]

Then $v^{(1)}=S_{0}$ and $v_{2}$ is the solution to

[TABLE]

which is positive due to assumption (11). Since $V$ is non-decreasing, we have

[TABLE]

We put

[TABLE]

Clearly,

[TABLE]

Moreover,

[TABLE]

and the same relation holds if $\hat{v}$ and $\check{v}$ are exchanged. Now we can formulate our main result about the large time behavior of solutions to (1), (2).

Theorem I Let $(x(t),v(t))$ be a solution of (1), (2) with positive initial data. Then

[TABLE]

and

[TABLE]

$i=1,\ldots,M$.

For the proof of this theorem see Appendix.

Note that $\hat{v}=\check{v}$ if $dV/dw>-1$, which is true when, for example, $D$ or $\gamma_{0}=\min_{i}\gamma_{i}$, are sufficiently large. Indeed, if $dV/dw\in(-1,0]$, the operator $v\to V(v)$ defined on $[0,S_{0}]$ is a contraction and therefore the iterations $v^{(k)}$ converge to the same limit. This observation implies the following

Corollary I

*For sufficiently large $D>0$ or $\gamma_{0}>0$ all the solutions $(x(t),v(t))$ of (1), (2) with positive initial data converge, as $t\to\infty$, to the unique equilibrium point defined by eqs. (9) and (10). *

III.2 Local stability

Consider now the problem of stability of equilibrium states $(x_{1},\ldots,x_{M},v_{eq})$ defined by (9), (10). Denote by ${\bf I}_{eq}$ the set of indices $i$ for which $\phi_{i}(v_{eq})-r_{i}>0$ and by $N_{eq}$ the number of such indices. Then $x_{i}>0$ when $i\in{\bf I}_{eq}$.

One can show that the eigenvalues of the linear approximation of (1), (2) at the equilibrium point $(x_{1},\ldots,x_{M},v_{eq})$ satisfies the equation (see Appendix):

[TABLE]

where

[TABLE]

and $P_{i}(v)=\phi_{i}(v)-r_{i}$.

Let us show that $Re\lambda<0$. In fact, taking the complex conjugate to (17) and summing these equations we have

[TABLE]

where

[TABLE]

This implies that

[TABLE]

Thus the equilibrium point $(x_{1},\ldots,x_{M},v_{eq})$ is locally stable for all $D$.

IV Extinctions

System (1), (2) does not take into account species extinctions due to extinction thresholds. Here we present a model describing this effect. The system thereby handles the evolution to the final set of species. We follow KVV with essential simplifications since we do not take into account the emergence of new species. We start from random values of the model parameters.

Main parameters of our model in this section are the coefficients $r_{i},K_{i},a_{i}$ and $\gamma_{i}$. Let us introduce the vector parameter ${\bf P}_{i}=(r_{i},K_{i},a_{i},\gamma_{i})$. Note that $c_{i}$ is a species specific parameter not necessary to include in this analysis and assumed to be fixed.

Let ${\bf P}=(P_{1},P_{2},P_{3},P_{4})$ be a random vector with a probability density function $\xi({\bf P})$. This means that the values ${\bf P}_{i}$ are defined by random sampling, i.e., the parameters of the species are random independent vectors ${\bf P}_{i}$ that are drawn from the cone ${\bf R}^{4}_{+}=\{{\bf P}:\ P_{1}>0,P_{2}>0,P_{3}>0,P_{4}>0\}$ by the density $\xi$. Our assumption to $\xi$ can be formulated as follows:

Assumption I. The probability density function $\xi$ is a continuous function with a support, which has a compact closure in the positive cone ${\bf R}^{4}_{+}$.

The function $\xi$ is positive on $S_{\xi}$, where $S_{\xi}$ is an open and bounded set. The closure of $S_{\xi}$ we denote by $\bar{S}_{\xi}$. Assumption I implies that the mortality rates do not approach zero and resource consumption is restricted. It is supposed that initial data $\bar{x}_{i}=x_{i}(0)$ are random mutually independent numbers drawn according to a density distribution

[TABLE]

with the mean $\bar{X}$ and the deviation $\sigma_{X}$. The random assembly of the species defines an initial state of the ecosystem for $t=0$.

In order to describe species extinction we introduce a small positive parameter $X_{ext}$ being an extinction threshold. We represent the set of indices $I_{M}=\{1,2,...,M\}$ as a union of the two disjoint sets:

[TABLE]

Here $S_{e}(t)$ is the set of indices of species which exist at the time $t$, and $S_{v}(t)$ is the set of indices of species, which have disappeared by the moment $t$. Let $N(t)$ denote the number of species in $S_{v}(t)$ at the moment $t$, $N(0)=M$. We assume that $S_{e}(0)=\{1,2,...,M\}$, and $S_{v}(0)=\emptyset$. In our model the species with abundance $x_{k}$ vanishes at the moment $t_{*}$ if $x_{k}(t_{*})=X_{ext}$ and $x_{k}(t)>X_{ext}$ for $t<t_{*}$. The parameter $X_{ext}$ can be interpreted as a threshold for species abundances.

The time evolution of the sets $S_{e}(t)$ and $S_{v}(t)$ can be described as follows.

(A) if the $k$-th species vanishes at a certain moment $t_{*}$, i.e. $k\in S_{e}(t)$ for $t<t_{*}$ and $x_{k}(t_{*})=X_{ext}$, then the index $k$ moves from $S_{e}(t)$ to $S_{v}(t)$ at this moment $t=t_{*}$ and we put $x_{k}(t)=0$ for all $t>t_{*}$;

(B) we assume that the evolution stops at the moment $t_{end}$, if at this moment $S_{e}(t)=\emptyset$.

With modifications described above, eqs. (1), (2) define the dynamics as follows. Within each time interval $(t_{*},T_{*})$ between the subsequent species extinctions the dynamical evolution of $x_{i}(t)$ is defined by the system (1), (2).

The quantity $N(t)$ is a piecewise constant decreasing function, therefore, there exists a limit

[TABLE]

where $N_{f}$ is the number of species, which survived to the limit state (note that it is possible that $N_{f}=0$).

Let us introduce the parameters

[TABLE]

and assume that

[TABLE]

for some $i$. Condition (22) means that the resource supply is large enough for existence of a positive equilibrium.

V Dynamics of the model with extinctions

By (20) there exists a time moment $T_{f}$ such that all extinctions have occurred and thus we can use Theorem I and its corollary for the remaining species. According to section IV $S_{e}(T_{f})$ is the set of indices corresponding to the species, which exist for all $t>0$. That set contains $N_{f}=N(T_{f})$ indices. We modify equation (9) as follows:

[TABLE]

where

[TABLE]

Corollary II For sufficiently large $D>0$ or $\gamma_{0}>0$ all the solutions $(x(t),v(t))$ of (1), (2) with positive initial data converge, as $t\to\infty$, to an equilibrium point defined by eqs. (10), (23), and (24). That equilibrium depends on the set of remaining species $S_{e}(T_{f})$.

The assertion follows from the arguments at the beginning of this section and Corollary I.

Note that the set $S_{e}(T_{f})$ depends on initial data, therefore, in contrast to Theorem I, we have a number of possible final equilibria. To show this, let us consider the following situations. Let $M=3$ and for $X_{ext}$ all three species survive, thus, $N_{eq}=3$.

Let $X_{ext}>0$ and $x_{3}(0)=X_{ext}+\kappa$, where $\kappa>0$ is a small number. We assume that $x_{1}(0)-X_{ext}$ and $x_{2}(0)-X_{ext}$ are not small. Suppose moreover that $D(S_{0}-v(0))-F_{M}(v(0))<0$ and $|D(S_{0}-v(0))-F_{M}(v(0))|>>\kappa$. Then it is clear that $x_{3}$ will go extinct within a short time period and thus $3$-th species is not involved in the set $S_{e}(T_{f})$ of final equilibria. If $\kappa$ is not small, then the set contains the $3$-th species.

VI Concentration of species traits

Let us consider the case of arbitrary parameter values, supposing that the initial number of species $M>>1$. For each $\epsilon>0$ let us denote by $W_{\epsilon}(z)$ the set of the points in $\bar{S}_{\xi}$, which lie in the ball of radius $\epsilon$ centered at $z=(r,K,a,\gamma)$. The $\epsilon$-neighborhood $W_{\epsilon}(B)$ of a subset $B\subset\bar{S}_{\xi}$ is the union of $\epsilon$-neighborhoods $W_{\epsilon}(z)$ taken over all the points $z\in B$.

In the set $\bar{S}_{\xi}$ we introduce the partial order $\leq_{e}$: $(r_{i},K_{i},a_{i},\gamma_{i})\leq_{e}(r_{j},K_{j},a_{j},\gamma_{j})$ if $a_{i}\leq a_{j}$, $r_{i}\geq r_{j}$, $K_{i}\geq K_{j}$ and $\gamma_{i}\geq\gamma_{j}$.

Consider the points $z_{*}=(r_{*},K_{*},a_{*},\gamma_{*})$, which are maximal with respect to the order $>_{e}$ in the set $\bar{S}_{\xi}$. Since that set is closed, bounded from below with respect to $K,r,\gamma$, and bounded from upper with respect to $a$, the set $B_{*}$ of the points $z_{*}$ is not empty. It is clear that $B_{*}$ is a subset of the boundary $\partial S_{\xi}$.

Theorem III (Concentration of traits). *Let Assumption I and (22) hold and $\epsilon>0$ be a number. Then the parameters $a_{i},r_{i},\gamma_{i}$ and $K_{i}$ of species $x_{i}$ such that $x_{i}(t)>X_{ext}$ for all $t>0$ lie in the domain $W_{\epsilon}(B_{*})$ with the probability $Pr_{M}(\epsilon)$ such that $Pr_{M}(\epsilon)\to 1$ as $M\to+\infty$. *

Proof can be found in Appendix.

If the set $B_{*}$ is a singleton (i.e. consists of a single point), then we have the concentration trait effect, i.e., all essential parameters of ecosystem become almost identical as a result of extinctions. Note that the set $B_{*}$ is a singleton in the case when $S_{\xi}$ is a box, i.e.,

[TABLE]

The set $B_{*}$ can have a more complicated structure, it may be a union of isolated points or a curve. Also note that even in the singleton case species may differ in coefficients $c_{i}$.

VII Limits of biodiversity in stress environment

The following assertion gives us an information on limits of biodiversity for arbitrary parameter values and our results are valid for arbitrary system dynamics: we do not use here no assumptions on existence of globally attracting equilibria. Remind that $N_{f}$ is the number of species, which survive as $t\to+\infty$, i.e., the corresponding abundances $x_{i}(t)\geq X_{ext}$ for all $t\geq 0$.

Proposition II The number $N_{f}$ is bounded by a constant independent of $M$, namely

[TABLE]

where $[x]$ denotes the integer part of $x$, $c_{0}=\min c_{i}$, and

[TABLE]

Proof. First we use an idea from Hofbauer . Let $\langle F\rangle_{T}=T^{-1}\int_{0}^{T}F(s)ds$ be the average of a function $F$ on $[0,T]$. The average of $F$ on $[0,+\infty)$ we denote by $\langle F\rangle$. By averaging of (2) one obtains

[TABLE]

Since the left-hand side here tends to [math] as $T\to+\infty$ eq. (26) leads to

[TABLE]

that in turn entails the estimate

[TABLE]

where $K_{*}=\max_{i\in S_{e}(T_{f})}K_{i}$. Consider the equation in (1) with the index $i\in S_{e}(T_{f})$ for which $K_{i}=K_{*}$. Dividing both sides there by $x_{i}$, averaging and using that $x_{i}$ is bounded and separated from zero by $X_{ext}$, we get

[TABLE]

Hence

[TABLE]

This together with (28) leads to (25).

To find more precise estimates we assume that coefficients $c_{i},a_{i},\gamma_{i}$ and $r_{i}$ satisfies

[TABLE]

[TABLE]

where $a,c,\gamma,r$ are characteristic values of the corresponding coefficients, $C_{\pm}$ are positive constants independent of $M,a,c,\gamma,r$. Let us introduce the stress parameter by

[TABLE]

To simplify the statement, we also suppose that $K_{i}=K$. The general assertion on the trait concentration can be formulated as follows:

Proposition III Suppose Assumption I and condition (22) hold. Let $i,j$ be two indices such that the corresponding species abundances $x_{i}(t),x_{j}(t)$ satisfy $x_{i}(t)>X_{ext},x_{j}(t)>X_{ext}$ for all $t\geq 0$. Then

[TABLE]

where $C_{0}>0$ does not depend on $a,c,\gamma,r$ and $X_{ext}$.

Proof. Consider the species such that $x_{i}(t)>X_{ext}$ for all $t\geq 0$. The corresponding set of indices we denote by $S_{e}$. Averaging equations (1) for the species $x_{i}$ with $i\in S_{e}$ we obtain the following relation:

[TABLE]

Furthermore, we divide (1) on $x_{i}$ and average the obtained equation that gives

[TABLE]

The Cauchy inequality implies $\langle x_{i}^{2}\rangle\geq\langle x_{i}\rangle^{2}$. Therefore, (34) and (35) entail

[TABLE]

where, for brevity, we use notation $\Phi=\Phi(v,K)$. By (2) we have

[TABLE]

We observe that

[TABLE]

By the above identity and (34) one has

[TABLE]

That relation and (36), (37) lead to the inequality

[TABLE]

that, by (29), can be rewritten as follows:

[TABLE]

where $\beta_{i}=c_{i}a_{i}\gamma_{i}(ca\gamma)^{-1}$ are bounded coefficients independent of $a,\gamma,c$. Estimate (39) entails

[TABLE]

for some $C_{2}>0$, which is independent of $\gamma,a,c,r$. By (35) $\Phi-p_{i}$ is positive and hence the index ${+}$ in (40) can be removed. Combining (40) for $l=i$ and $l=j$ and taking into account that $\langle\Phi\rangle>p_{0}+\delta_{0}$ one has (33).

Let us derive an estimate of $N_{f}$ via $P_{stress}>>1$ and the average $\langle\Phi\rangle$. We suppose that $p_{i}=p_{0}+(i-1)\Delta p$, $\Delta p<<1$ and $\beta_{i}=\beta=O(1)$. Then estimate (39) implies

[TABLE]

This calculation is consistent with numerical simulations. For $\gamma=0.00001$, when all other parameters have the order $1$, we obtain a strong concentration effect (see Fig. 1). Computations were made for a population of $M=50$ species, where random parameters chosen as explained above.

As a measure of the trait concentration, we can use the quantity

[TABLE]

Then the initial $Var(p)\approx 0.6$ but for $4$ remaining species with large abundances we have $Var(p)=0.07$. We see that these four species are abundant whereas all other species are extremely rare. Note that these asymptotic results can be generalized in the case of different $K_{i}$.

VII.1 Mass extinctions: an analytical approach

Relation (25) allows us to describe, in an analytical way, mass extinctions. Mass extinctions may result as a consequence of a sharp change of some environmental parameter. It is natural to assume that climate variations or other abiotic ones can reduce the resource supply level $S_{0}$. Assume for example that this reduction is $\Delta S_{0}>0$ and the new resource supply $S_{new}=S_{0}-\Delta S_{0}$ satisfy

[TABLE]

The last equation implies that for sufficiently large $\Delta S_{0}$ even all species may go extinct. We therefore refer this level of the resource $S_{new}$ as a catastrophic level.

Note that these analytical results show that there are interesting phenomena. Firstly, let us compare two ecosystems. One is a random assembly of many species where the variation $Var(p)$ defined by (42) is large, and the other ecosystem is a result of long evolution leading to the concentration, i.e., $Var(p)$ is small. We find that the concentrated system is more stable with respect to variations in the resource. Namely, a sharp change of $S_{0}$ will kill many more species in the first ecosystem then in the second one. Secondly, assume that catastrophes do occur several times yet with a fairly long time in between. Each catastrophe will reduce the biodiversity yet with less and less probability since the concentration effect becomes stronger and stronger.

To investigate more realistic situations when $K_{i}$ are different and the parameters are random, we performed numerical simulations described in the following section.

VIII Numerical simulations

In numerical simulations, the parameters are chosen as follows. The coefficients $a_{i}$ are independent and identically distributed (i.i.d.) random quantities such that each $\ln a_{i}$ is a normally distributed number with the mean $E_{a}=1$ and the standard deviation $\sigma_{a}=0.2$. This means that each $a_{i}$ has the same log-normal distribution, $a_{i}\in\ln{\bf N}(E_{a},s_{a})$.

Similarly, the coefficients $r_{i}$ are i.i.d. random quantities, $\ln r_{i}$ is normally distributed on $[0,1]$ number with the mean $E_{r}=0.1$ and the standard deviation $\sigma_{r}=0.03-\ln(20)$. The parameters $D=0.1,K=4,S=30$ and $\gamma_{i}=\gamma=0.001$. The coefficients $c_{i}$ are random numbers uniformly distributed on $[0,1]$ and normalized in such a way that $\sum c_{i}=1$. The initial data $X_{i}$ are i.i.d. random numbers distributed log-normally, $X_{i}=\exp(y_{i})$, $y_{i}\in{\bf N}(E_{y},s_{y})$ with parameters $E_{y}=1,s_{y}=0.3$. We suppose that all $c_{k},r_{j}$ and $a_{l}$ are mutually independent.

For these random species communities and $N=50$ we observe oscillations and then a convergence to an equilibrium (see Fig. 2).

By simulations we have considered the dependence of biodiversity and concentration trait effect on the stress parameter for populations with random parameters $K_{i},p_{i}$. The results are consistent with analytic considerations of the previous section and can be illustrated by Fig. 3.

In the numerical simulations of biodiversity we can also observe the trait concentration. For the example illustrated by Fig. 3 the variation $Var(p)$ decreases very strongly as a result species extinctions and this reduction increases as the stress parameter increases.

IX Dynamics of limit ecosystem

According to Theorem III, if the set $B^{*}$ consists of a single point, for $M>>1$ the limit system (that appear as a result of many extinctions) has the property $p_{i}\approx p,K_{i}\approx K$, where $p,K$ are some parameter values. To understand dynamics of that system, we consider system (1), (2) in the case $p_{i}=p,K_{i}=K$. Let us introduce a new variable $Q=-pt+\int_{0}^{t}\Phi(v(s),K)ds.$ The variable $Q$ is an analogue of ”quality of life” introduced in Volterra for the linear case $\Phi(v,K)=v$. This case is studied in Koz1 . Results of Koz1 can be extended to our limit model. We seek solutions to eqs. (1) in the form

[TABLE]

where $C_{i}$ are new unknowns. From (1) one obtains

[TABLE]

that gives

[TABLE]

By solving these equations, we find

[TABLE]

Using the last relation, by (2) one obtains

[TABLE]

where $P=dQ/dt+p$.

After some straight forward computations eqs. (1), (2) reduce to the system

[TABLE]

[TABLE]

where

[TABLE]

Eq. (45) describes a nonlinear oscillator with a damping term and nonlinearities with a time delay. Note that $f$ depends on initial data $\bar{x}_{i}$. So, we see that the limit ecosystem can be considered as a nonlinear oscillator with a friction and a memory. The oscillator state is determined by two variables: $P$ and $Q$. The first variable is a difference between the normalized species consuming rate $\Phi(v,K)$ and the normalized species mortality rate, i.e., admits a biological interpretation. This variable can be called Malthusian parameter. The second variable $Q$ does not admit a simple explicit interpretation. It is a generalization of the quality of life introduced by Volterra’s. Note that $Q$ is the integral of $P$, i.e., it can be considered as an integral Malthusian parameter. Since this is a parameter expressing a trait over long time we can call $P$ the sustainable Malthusian parameter.

Eq. (45) can be simplified in two cases: for $\gamma=0$ and for bounded times, $t<<\ ln(\gamma^{-1})$ (an initial stage) and for $t>>1/\gamma$ (large times, the final stage). In the first case from (47) we have

[TABLE]

We obtain an oscillator, which is a perturbed Hamiltonian integrable system without memory. In this case for small $D$ the solutions of (45) tend to the equilibrium in an oscillating manner ( see Fig. 2).

X Conclusion

In this paper, we have investigated a model of ecosystems exploiting a single resource and interacting with the environment.

Until May’s seminal works May1 ; May2 , ecologists believed that large complex ecosystems, involving a larger number of species and interconnections, are stable. May May1 ; May2 considered a community of $S$ species with connectance $C$ that measures the number of realized links with respect to the number of all possible links. R. May’s analysis of local stability of an equilibrium gave quite revolutionary results that inspired a great discussion. It was shown that for large systems with random interaction parameters the instability can occur for large $C$. More connected communities are more unstable. This approach is developed in Alles1 ; Alles2 , where more complicated networks with interactions of different types (predator-prey, amensalism, mutualism, competition) were studied.

All these fundamental results hold under the assumption that, at an equilibrium, ecosystems have a random structure, namely, the entries of the matrix, that defines the linearization of system at the equilibrium, are distributed according to smooth densities, for example, Gaussian ones.

In this paper, we use a similar assumption but on the initial choice of species traits. The initial distribution of species traits is defined by continuous densities with non-empty supports, i.e., roughly speaking, the species traits are distributed homogeneously in a domain. We show that in the evolution process the distribution of species traits becomes more concentrated when ecosystem evolves under a stress or as a result of species extinctions. During the evolution process the domain of species trait localization shrinks. That small domain of localization means that species become more and more similar (as in Sim ). In contrast to Sim , we do not use any specific assumptions on the adaptation of system parameters. We have found a parameter, which defines the stress level. This parameter depends on the supply level, turnover rate, and resource consuming intensity.

The most interesting effect of species trait concentration is as follows. For large times a stable and simple limit ecosystem appears just because of most species goes extinct under stress. For large times ecosystem dynamics and extinctions of species under stress produces a self-organized community consisting of species with close consumer efficiencies (note that these species can be different in other traits). In some cases, dynamics of this limit community can be described by a simple equation, which describes a nonlinear oscillator with a friction and a memory, which is close to a Hamiltonian system. We have found an asymptotic approach to study this system. Note that the reduction mechanism differs from previously found one in Barab ; Koz1 . In Barab a mean field approach is applied to complex ecosystems and gene networks. This approach exploits the system topology, when species (genes) can interact with many others. Complicated systems of equations were reduced to a single differential equation of the first order. Such equations do not exhibit time oscillations whereas our equation simulates a perturbed nonlinear oscillator and it can describe slowly decreasing oscillations. In Koz1 a reduction to Hamiltonian systems is also based on topological properties of interactions in ecosystems. So, the reduced descriptions of complex systems proposed in Barab ; Koz1 can be called a topological one. In contrast to Barab ; Koz1 , in the present paper the reduced description is based not only on the system topology (i.e., the fact that species share the same resource) but also on others phenomena: extinctions and selection by a tough environment (which can be measured by the stress parameter).

These results can be useful for understanding why ecosystems where species feeds on few resources can have a large biodiversity, and how mass extinctions depend on environment and ecosystem parameters. The intriguing effect is that we observe a picture similar to statistical physics: the state of a ecosystem which arises as a product of a long evolution can be described by two quantities $P,Q$ having a biological interpretation. Namely, $P$ can be called Malthusian parameter, and this quantity determines a balance between mean mortal and growth rates. The second quantity $Q$ can be named sustainable Malthusian parameter, and it can be obtained by integrating $P$ over time.

We have computed analytically the number of finally coexisting survived species and how to this number depends on main ecosystem parameters (the resource supply, the mortality rates, the resource turnover etc). It is shown that the main quantity that determines final biodiversity is the stress parameter.

XI Acknowledgements

The authors are thankful to Referee for useful remarks. The second author was supported by Linkoping University and by the Government of the Russian Federation through mega-grant 074-U01, and RFBR grant 16-01-00648.

XII Appendix

XII.1 Proof of formula (17)

Assume that all $P_{i}(v_{eq})>0$. Then the eigenvalue problem for the linear part of the right-hand side of (1), (2) at the equilibrium point $(x_{1},\ldots,x_{M},v_{eq})$ has the form

[TABLE]

[TABLE]

where $(X_{1},\ldots,X_{M},V)$ is an eigenvector corresponding to the eigenvalue $\lambda$. Solving the first system with respect to $X_{i}$ and inserting the solution into the second equation we obtain

[TABLE]

and

[TABLE]

Since $\gamma_{i}x_{i}=P_{i}(v_{eq})$ at the equilibrium point we arrived at (17).

If some of $P_{i}(v_{eq})$ are non-positive the corresponding terms in (49)-(51) are zeros and we again arrive at (17).

XII.2 Proof of Theorem I on global stability of positive solutions

We apply a special method based on the theory of decreasing operators in Banach spaces (see Ban and references therein) that allows us to prove this assertion without any additional assumptions. This approach is applicable here due to special properties of monotonicity of our problem.

Let us rewrite (12) as follows:

[TABLE]

where the operator $V$ is described in subsection III.1. We remind that $V(\bar{v})$ is a decreasing function in $\bar{v}$.

Our next step is to rewrite system (1), (2) as an integral equation for unknown function $v(t)$. Let $w(t)$ be a given non-negative, continuous, bounded function on $[0,\infty)$ having a limit $\bar{w}$ at infinity. We can resolve eqs. (1) (with $v$ replaced by $w$) following section IX. As a result, we obtain

[TABLE]

where

[TABLE]

and

[TABLE]

One can verify that for $x_{i}(0)>0$

[TABLE]

where $X_{i}$ is defined in subsect. III.1

Next, we can solve eq. (2) with respect to $v$, where $x_{i}$ is given by (53) and $v(0)=v_{0}$. We denote this solution by ${\bf V}(t)={\bf V}(w(\cdot))(t)$. We cannot write this solution explicitly but we need in what follows only some of its properties. First, this solution is a decreasing function with respect to $x_{i}$ and consequently with respect to $w$. Second,

[TABLE]

where $\bar{v}=V(\bar{w})$. Thus the unique solution to the problem (1), (2) with the Cauchy data (4) can be obtained by solving the following fixed point problem

[TABLE]

and then

[TABLE]

To solve the equation $v={\bf V}(v(\cdot))$ in the class of bounded, continuous, non-negative functions (denoted by ${\bf B}$), we use the following iterations

[TABLE]

Then

[TABLE]

(here $\leq$ denotes the partial order on ${\bf B}$: $v\leq u$ if $u(t)\leq v(t)\ \forall t\in[0,T]$).

To show the convergence of the odd and even iterations, we observe that we can consider the fixed point equation (54) on a finite interval $(0,T)$. Now the operator ${\bf V}:C[0,T]\rightarrow C[0,T]$ is compact and hence the odd and even terms of sequences converge on $[0,T]$ for each $T$. We introduce their limits

[TABLE]

Then ${\bf V}(\check{V})=\hat{V}$ and ${\bf V}(\hat{V})=\check{V}$. Let $\hat{x}_{i}$ be given by (53) with $w=\hat{V}$ and $\check{x}_{i}$ be given by (53) with $w=\check{V}$. Then the vector function $(\check{x}_{1},\ldots,\check{x}_{M},\hat{V})$ satisfies the problem

[TABLE]

and the functions $(\hat{x}_{1},\ldots,\hat{x}_{M},\check{V})$ are solutions of

[TABLE]

Moreover, the last two systems have the same Cauchy data. Taking differences we obtain a homogeneous Cauchy problem for $(\hat{x}_{1}-\check{x}_{1},\ldots,\hat{x}_{M}-\check{x}_{M},\hat{V}-\check{V})$ and by uniqueness for the Cauchy problem we obtain that $\check{V}=\hat{V}$.

Let us turn to the asymptotic behaviour of the fixed-point solutions. Let $\bar{v}_{k}=\lim_{t\to\infty}v_{k}(t)$. Then

[TABLE]

This proves inequalities (15) and (16) and completes the proof of Theorem I.

XII.3 Proof of Theorem III

It proceeds in three steps.

Step 1: Monotonicity of species abundances.

Consider a point $\bar{z}=(a_{i},r_{i},\gamma_{i},K_{i})$, which are not contained in $W_{\epsilon}(B_{*})$, and the corresponding species population $x_{i}(t)$. Suppose that for all $t\geq 0$ we have

[TABLE]

Consider $j$-th species with parameters $(a_{j},r_{j},\gamma_{j},K_{j})$ and the species abundance $x_{j}(t)$. We assume that

[TABLE]

Then

[TABLE]

Indeed, let us consider equations for $x_{i},x_{j}$:

[TABLE]

[TABLE]

If (57) is violated then there is a time moment $t_{1}>0$ such that

[TABLE]

But

[TABLE]

due to the first inequality in (60) and (56). The last inequality contradicts the second inequality in in (60), thus, (57) is proved.

Inequality (57) shows that if the species $x_{i}$ survives for all times, then all the species with parameters satisfying (56) also survive for all $t>0$.

Step 2: a priori boundness of biodiversity. Here we use Proposition II. The number $N_{s}$ of species, which survive for all times, a priori bounded by the system parameters and does not depend on $M$ as $M\to\infty$. We refer the corresponding set of species parameters as ${\mathcal{P}}_{s}$. Due to Prop. II,

[TABLE]

where $C>0$ is independent of $M$.

Step 3. Let us consider the $\epsilon$-neighborhood $W_{\epsilon}(B_{*})$. Suppose there exists a point $\bar{z}\notin W_{\epsilon}(B_{*})$. The initial data $x_{i}(0)$ for the corresponding species we denote by $\bar{x}_{i}$. Then, according to Step 1, the set ${\mathcal{P}}_{s}$ contains all points $z$ from $W_{\epsilon}(B_{*})$ such that $z\geq_{e}\bar{z}$. We denote the set of such points by $W_{\epsilon,\bar{z}}(B_{*})$. Note that due to the conditions to the set $S_{\xi}$ (see Assumption I), the set $W_{\epsilon,\bar{z}}(B_{*})$ contains a small open ball. Therefore, since $\xi$ is positive on the interior of $S_{\xi}$ (see Assumption I), we have

[TABLE]

The number $\delta_{\epsilon,\bar{z}}$ is independent of $M$. Consider the event $E=AB$ where $A$ is the event that the species parameters lie in $W_{\epsilon,\bar{z}}(B_{*})$ and $B$ is the event that initial data $x_{i}(0)>\bar{x}_{i}\ \forall i$. The events $A$ and $B_{s}$ are independent and $Prob(A)>0$ due to the above estimate for $J$. According to hypothesis on the random choice of $x_{i}(0)$ we also have $Prob(B)>0$. Therefore, $Prob(E)=q>0$.

Consider the event $E_{M,N_{s}}$ that among $M$ species there are not more than $N_{s}$ species such that the corresponding species parameters lie in $W_{\epsilon,\bar{z}}(B_{*})$ and that initial data $x_{i}(0)>\bar{x}_{i}\ \forall i$. The probability of $E_{M,N_{s}}$ can be computed by the Bernoulli relation, and we have

[TABLE]

We see that $E_{M,N_{s}}\to 0$ as $M\to\infty$ and the Theorem III is proved.