Large ecosystems in transition: bifurcations and mass extinction

# Large ecosystems in transition: bifurcations and mass extinction

Ivan Sudakov Sergey A. Vakulenko Dubrava Kirievskaya Kenneth M. Golden University of Dayton, Department of Physics, 300 College Park, SC 111, Dayton, OH 45469-2314 USA Institute of Problems in Mechanical Engineering, Russian Academy of Sciences, Bolshoy pr., 61, V.O., St. Petersburg 199178, Russia University ITMO, Kronverkskiy pr., 49, St. Petersburg 197101, Russia University of Utah, Department of Mathematics, 155 S 1400 E RM 233, Salt Lake City, UT 84112-0090, USA
###### Abstract

We propose a model of multispecies populations surviving on distributed resources. System dynamics are investigated under changes in abiotic factors such as the climate, as parameterized through environmental temperature. In particular, we introduce a feedback between species abundances and resources via abiotic factors. This model is apparently the first of its kind to include a feedback mechanism coupling climate and population dynamics. Moreover, we take into account self-limitation effects. The model explains the coexistence of many species, yet also displays the possibility of catastrophic bifurcations, where all species become extinct under the influence of abiotic factors. We show that as these factors change there are different regimes of ecosystem behavior, including a possibly chaotic regime when abiotic influences are sufficiently strong.

###### keywords:
multispecies ecosystems, dynamical systems, feedback mechanisms, Lotka-Volterra model, bifurcations, mass extinction.
\newdefinition

rmkRemark \newproofpfProof \newproofpotProof of Theorem LABEL:thm2

## 1 Introduction

Models of ecosystems form an important class of dynamical systems generating complex dynamics, bifurcations and strange attractors (Ulanowicz and Kemp, 1979). However, modeling these large systems is made difficult by rapid, large scale biological evolution and gaps in observations to use for comparison. Also, there is uncertainty in how to set up reliable experiments on such ecosystems.

Recent observations have shown that climate change may be a leading factor influencing ecosystem behavior (Walther, 2010). Large multispecies marine ecosystems are sensitive indicators of climate change (Doney et al., 2012; Kedra et al., 2015). As a key part of the global ecosystem, they influence climate feedback processes and possible tipping points (Selkoe et al., 2015). A well studied example is the ocean ecosystem, where phytoplankton are the main resource for many species. Phytoplankton populations play an important role in the dynamics of the climate system through the oceanic carbon cycle – by removing about half of all carbon dioxide from the atmosphere during photosynthesis (Field et al., 1998). Previous studies (Arhonditsis and Brett, 2004; Travers et al., 2007) have shown that phytoplankton communities respond to climate warming through changes in diversity and productivity. However, it was recently determined (Toselandet al., 2013) that changing the climate temperature directly impacts the chemical cycles in plankton, affecting the system as much as nutrients and light.

We consider here a model of a large ecosystem where many species share few resources. It extends the model of phytoplankton species competition in (Huisman and Weissing, 1999), by taking into account that the resources depend on environmental factors, in particular, climate, as well as self-limitation and competition effects. Our aim is to explore the connections among complexity of the temporal behavior, biodiversity, and the structure of the climate–ecosystem interaction.

Note that competition may occur as a result of the following mechanism (Roy and Chattopadhyay, 2007). There are a number of species of phytoplankton which have the ability to produce some toxic or inhibitory compounds. These toxic materials compensate for the competitive disadvantages among phytoplankton species which leads to self-limitation effects. Moreover, resource levels may depend on the environment via temperature or greenhouse gas concentration.

Many mathematical models (Hofbauer and Sigmund, 1988; Takeuchi, 1996; Zeeman, 1995) show that only a single species can survive in an ecosystem for certain fixed parameters. Biologically, this is the competitive exclusion principle. In the framework of the phytoplankton model, it is known as the so-called plankton paradox studied in many interesting works (Hutchinson, 1961; Tilman, 1977; Huisman and Weissing, 1999; Irigoien et al., 2004). In particular, it is sometimes observed in nature that numerous species can coexist while depending on the same resource, even though competition tends to exclude species. In fact, in contrast to the exclusion principle, we observe here the coexistence of many plankton species sharing the same niche and resources. Numerical simulations (Hutchinson, 1961; Tilman, 1977; Huisman and Weissing, 1999) have shown that in such systems chaos and unpredictable behavior occur. In (Hsu et al., 1977; Smith, 1981) it was shown that temporal variability of the nutrient supply can lead to coexistence of species.

The environment may alter the distribution and abundance of the species in a population. Such effects have been studied in terms of internal processes within the population, like competition for resources and conditions for chemical reactions. However, current models have not been linked to feedback with the environment. Feedback between a population and the environment can occur as a result of changes in abiotic factors such as temperature, nutrient concentrations, and light intensity.

The main results of this paper show that the population dynamics depends sharply on feedback with the environment. For simplicity, hereinafter we refer to this as climate – ecosystem feedback. If the abiotic factor is temperature , for example, then it is natural to talk about the feedback between an ecosystem and the climate system, which can be parameterized as a function of a rate of change of the resource supply with respect to temperature. If that feedback is negative – where species abundance decreases resources – then an ecosystem can support a number of species and the dynamics is relatively simple (non-chaotic and non-periodic). If the feedback is positive – where species abundance increases resources – then for a sufficiently large feedback level there are possible mass extinctions which occur suddenly, and moreover, there are possible chaotic or periodic dynamics.

The paper is organized as follows. In the next section we formulate the standard model of species coexistence and the extended model, which takes into account climatic factors. Further, in section 3 we prove a general assertion on the existence of an attractor for this model. In section 4 it is shown that for large turnover rates the system admits an asymptotic solution and, under additional assumptions, can be reduced to the Lotka-Volterra model (Vakulenko, 2013; Kozlov and Vakulenko, 2013). This model is well studied (Hofbauer and Sigmund, 1988; Takeuchi, 1996; Zeeman, 1995, 1998) and known results allow us to describe the influence of climate and climate warming in large ecosystems (see section 5). In section 6, for the case of a single resource, we show that the global attractor consists of equilibria and derive an equation for the species abundances. This investigation is aimed at describing the influence of climate on biodiversity.

## 2 Models of large ecosystems

### 2.1 Standard model

Consider the following model of an ecosystem with species, which extends the model of resource competition in (Huisman and Weissing, 1999):

 dxidt=xi(−ri+ϕi(v)−N∑j=1γijxj),1≤i≤N, (1)
 dvdt=D(S−v)−N∑j=1cjxjϕj(v), (2)

where

 ϕj(v)=ajvKj+v,aj, Kj>0, (3)

is the specific growth rate of species as a function of the availability of the resource (also known as Michaelis-Menten’s function), are species abundances, are the species moralities, is the resource turnover rate, is the supply concentration of the resource, and can be interpreted as the supply rate. The dynamics of the species depend on the availability of the resource, which in turn depends on the rate of resource supply and the amount of resource used by the species.

The coefficient is the content of the resource in the -th species. The constants define how different species share resources. Note that if all then the equation for becomes trivial and for large times , i.e., the resource equals the resource supply. We consider this system in the non-negative cone: , , where . The coefficients are specific growth rates and the are self-saturation constants.

We assume that the . The terms define self-regulation of species populations that restricts their abundances. In the case with these terms describe competition between species. These effects can appear as a result of an ability to produce some toxic or inhibitory compounds (Roy and Chattopadhyay, 2007). However, we admit the possibility of mutualistic interactions, in which case . Assumptions on are formulated below, at the beginning of Section 3.

For the case of resources, we have the more complicated equations

 dxidt=xi(−ri+ϕi(v)−N∑k=1γikxk),1≤i≤N, (4)
 dvjdt=Dj(Sj−vj)−N∑k=1cjkxkϕk(v),1≤j≤M, (5)

where , and the are smooth functions. We consider general satisfying the conditions

 ϕj(v)∈C1,0≤ϕj(v)≤C+, (6)

where is a positive constant, and

 ϕk(v)=0,∀k,v∈∂RM+, (7)

where denotes the boundary of the cone . Condition (6), in particular, means that forms a uniform upper bound for the . We assume that This model is widely used for primary producers like phytoplankton, and can also be applied to describe competition for terrestrial plants (Tilman, 1977).

When for all this system is equivalent to those in works where the plankton paradox is studied (Huisman and Weissing, 1999). The choice and for allows us to take into account self-limitation effects, which is important in these systems, as shown by Roy and Chattopadhyay (2007).

Below we use the notation . We define the scalar product in together with the corresponding norm by

 ⟨f,g⟩C=N∑j=1Cjfjgj,||f||2C=⟨f,f⟩C. (8)

This scalar product is defined for -component vectors and depends on non-negative coefficients , .

### 2.2 Extended standard model with climate influence

We extend the system (4) and (5) to describe potential effects connected with an influence of the climate. In fact, temperature has a significant effect on the maximum growth rate of phytoplankton (Richardson et al., 2000), and can be considered as a crucial factor in population dynamics.

For one and two species (), a model of climate influence was proposed in (Sekerci and Petrovskii, 2015). We consider the case of arbitrary . In certain aspects, however, our model is simpler than in (Sekerci and Petrovskii, 2015). In particular, we do not account for zooplankton and, therefore, do not take into account possible predator-prey interactions in an explicit form.

Let us assume that the resource supplies can depend on the environmental parameters, for example, temperature : . In turn, may depend on species abundances, for example, via albedo (Chapin et al., 2002). We assume, for simplicity, that this effect is linear:

 T=¯T+ΔT,ΔT=N∑k=1μkjxj, (9)

where are coefficients and is a reference temperature corresponding to the albedo of the ecosystem environment, such as the upper ocean, without the ecosystem influence. If the temperature variations induced by the species are small, we have

 Sk=¯Sk(¯T)+ΔSk+O(ΔT2),ΔSk=N∑k=1bkj(¯T)xj,k=1,...,M, (10)

where . If all we are dealing with purely positive feedback (then species abundance increases resources), and if all one has purely negative feedback.

There is, however, an interesting case where some of the coefficients are positive numbers and others are negative (mixed feedback). For mixed feedback a cumulative effect of the climate-ecosystem feedback on the resource supplies may be small since the different terms in may cancel each other. On other hand, when the signs of the alternate, but these coefficients are sufficiently large, there may be complicated large time behavior. We discuss this problem in more detail in Section 5.

There are also possible alternative physical mechanisms leading to relations like (10). An important resource for phytoplankton is oxygen (Sekerci and Petrovskii, 2015). The production of oxygen is proportional to the phytoplankton concentration and depends on temperature .

Finally, the extended model takes the form

 dxidt=xi(−ri+ϕi(v)−N∑j=1γijxj),1≤i≤N, (11)
 dvjdt=Dj(Sj(x)−vj)−N∑k=1cjkxkϕk(v),1≤j≤M, (12)

where

 Sk(x)=¯Sk+N∑k=1bkjxj,k=1,...,M. (13)

This model is an approximation of the model with temperature dependent only up to the terms of order .

In the next section we show that under some assumptions this model is well posed.

## 3 General properties of the model

Let us first describe some sufficient conditions which guarantee that systems (1), (2), (4), (5) and (11), (12) are dissipative and have an attractor, and recall some basic notions. Since there are variations in the definition of attractor, for correctness, we follow (Hofbauer and Schreiber, 2010).

Let us consider the Cauchy problem defined by eqs. (11), (12) and positive initial data in (16) below. The solution with initial data (where the superscript denotes transpose) is unique and is defined for all (see Lemma 1). We then obtain the map defining a global semiflow , in a cone , which serves as a phase space.

Given an interval and a set , let . We denote . A set is invariant if for all , and forward invariant if for all . The omega limit set is the intersection of all over all . Given a forward invariant set a subset of is an attractor for the semiflow restricted to provided there exists an open neighborhood of such that .

The stable set of a compact invariant set is defined by

 Ws(K)={z∈C:ω(z)≠∅ and ω(z)⊂K}.

In other words, the stable set of consists of points where trajectories enter inside the set , and stay in for large times

The semiflow is dissipative if there exists an attractor such that . In other words, for dissipative semiflows the attractor is a minimal invariant set, which attracts all points. If the attractor consists of a single isolated point, then this point is stable in the standard Liapunov sense.

Define the matrix with the entries to satisfy one of the following conditions:

Assumption 1A. The matrix with the entries has a positive dominant diagonal:

 γii−∑j=1,...,N,j≠i|γij|=κi>01≤i≤N. (14)

Assumption 1B. The matrix has non-negative entries

 γij≥0,γii>0,1≤i,j≤N. (15)

Assumption 1A means that species self-regulation is stronger than species interaction, while assumption 1B implies that all species in our ecosystem compete. Let us show that the solutions to (11), (12) exist, and that they are non-negative and bounded.

Lemma 1. Assume the functions satisfy (6). Let us consider for eqs. (11), (12) the Cauchy problem with positive initial data for and positive initial resources

 xi(0)>0,vj(0)>0,∀i∈{1,…,N},∀j∈{1,…,M}. (16)

Then, if either assumption 1A or 1B holds, solutions of this Cauchy problem exist for all , are positive and bounded for large times , that is,

 00, (17)

where is a positive constant, , and

 0

where

 Vj(t)=¯Sj+¯bjX(t),¯bj=N∑i=1(bji)+.

Proof. For a proof, see the Appendix.

Due to boundness of solutions for large we then obtain the following corollary.

Theorem. Under the conditions of the previous lemma, system (11), (12) defines a global semiflow in the cone . This semiflow is dissipative and has a compact attractor.

## 4 Asymptotic approach

Our next step is to find asymptotic solutions of the system in (11) and (12), where the are defined by (10). We consider the case of large . Note that a reduction to a Lotka-Volterra system described below also holds for bounded and large resource supplies . To simplify the statement, we assume that for all . Let us make the change of variables

 vk=Sk(x)−~vk,τ=Dt. (19)

System (11) and (12) then takes the form

 dxidτ=ϵxi(−ri+ϕi(S(x)−~v)−N∑j=1γijxj), (20)
 d~vjdτ=−~vj−ϵUj(x,~v), (21)

where , and

 Uj(x,v)=N∑k=1cjkϕk(S(x)−~v)+N∑k=1bjk(ϕk(S(x)−~v)−rk−∑klγklxl). (22)

For small equations (20) and (21) form a typical system with slow variables and fast variables . We can find an asymptotic solution of (21), which has the form

 ~vj=ϵUj(x,0)+O(ϵ2). (23)

Finally, for the species abundances we obtain

 dxidt=xi(ϕi(S(x))−ri−N∑j=1γijxj)+O(ϵ). (24)

## 5 Qualitative analysis of large time behavior

If the coefficients are small, i.e., the feedback between the resource supply and the climate is weak, then the system (24) can be simplified by the Taylor expansion

 ϕi(S(x))=ϕi(¯S)+∑l=1,...,M∑j=1,...,N∂ϕi∂Sl(¯S)bljxj+....

Removing terms quadratic in , equation (24) reduces to the Lotka -Volterra system

 dxidt=xi(Ri−N∑j=1Aijxj). (25)

where

 Ri=ϕi(¯S)−ri,Aij=γij−M∑l=1ailblj, (26)

and

 ail=∂ϕi∂Sl(¯S). (27)

The Lotka-Volterra systems are very well studied (see, for example, (Hofbauer and Sigmund, 1988; Takeuchi, 1996)) and we can use these results to help understand how climate warming can affect ecosystems. We assume that 1B holds and consider the two limiting cases, the “weak climate” (WC) regime and the “strong climate” (SC) regime. The WC case corresponds to weak climate influence, where the ecosystem-climate interaction via the coefficients is much weaker than the competition effects associated with the coefficients . This means that all the , where is a characteristic magnitude of the entries .

In the SC case (regime of strong climate influence; coefficients determining climate feedback are stronger than the coefficients that define species interaction), we assume that .

In the WC case, system (25) is close to so-called competitive systems, which are well studied (Hirsch, 1985; Smith, 1981; Smith and Thieme, 1991; Zeeman, 1995, 1998). Under some conditions (Hirsch, 1985; Hofbauer and Sigmund, 1988; Zeeman, 1998) these systems exhibit no stable periodic or chaotic regimes: almost all trajectories converge to equilibria, which will be investigated in section 6 for the case of a single resource.

Consider the SC case. We set for all . Then equations (25) represent a Lotka-Volterra system of a special structure. An analysis (Hofbauer and Sigmund, 1988) shows that, for general , no more than species can coexist – an expression of the competitive exclusion principle. Mathematically this means that if then for some either the corresponding or as , i.e., the system is not permanent (Hofbauer and Sigmund, 1988). However, if the condition

 Ri=M∑k=1aikθk,∀i=1,...,N (28)

for some is fulfilled, then it is possible that all species can coexist. In this case system (25) can be studied by an idea proposed by Volterra (Volterra, 1931). We introduce new variables , named the qualities of life in (Volterra, 1931), where . Then eq. (25) reduces to a system involving only the variables (Kozlov and Vakulenko, 2013):

 dqjdt=Gj(q), (29)
 Gj(q)=−θi+N∑i=1bjiCiexp(−M∑j=1aijqj), (30)

where the are arbitrary positive constants. The species abundances can be expressed via by

 xi=Ciexp(−M∑j=1aijqj),i=1,…,N.

Note that and therefore the vector field depends on initial data and the species number . So, system (29) completely determines dynamics of .

The main results on system (29) can be outlined as follows (see (Kozlov and Vakulenko, 2013) for more details). Let be a compact connected domain in with a smooth boundary, be a compact smooth field on , and be a real number. Then there exist a number and coefficients and such that the corresponding field approximates in the domain in -norm with accuracy . This approximation result implies that system (29) with variables can generate all structurally stable dynamics in dimension . In particular, due to the Theorem on Persistence of hyperbolic sets (Ruelle, 1989), system (29) can exhibit all (up to topological orbital equivalences) hyperbolic dynamics, including periodic and chaotic, including for example, the Smale horseshoe, Anosov flows, etc.

Under condition (28) we find that the time behavior of solutions of system (25) depends sharply on . Assume that . Note that this assumption looks natural since it means that increases as a resource supply increases.

If it is possible that all species survive in an equilibrium state, and may be large. Although periodic and chaotic trajectories are impossible, we can observe multistability (coexistence of many equilibria).

For and of different signs, system (25) can have time periodic solutions and for this system can produce time chaotic solutions (we can then obtain all possible hyperbolic invariant sets of dimension ). If all or all we have no complex behavior for the trajectories and they are convergent. Therefore, the most interesting situation arises in the biodiversity case when have different signs. Finally, we conclude that in the SC regime there are possible chaotic phenomena and periodic oscillations if there exist at least three resources .

In the next subsection we will study the case and we will see that in this case Andronov-Hopf bifurcations are possible.

### 5.1 Bifurcations, complexity and biodiversity

If there exists a positive climate-ecosystem feedback, and , then time periodic (for ) or even chaotic (for ) behavior, as well as complicated bifurcations, can occur.

We consider two cases: (a single resource) and , and investigate the existence of different bifurcations, in particular, the Andronov-Hopf bifurcations. If there are possible saddle-node, pitchfork, and transcritical bifurcations, but the Andronov-Hopf does not occur. The main climate effect in the case is a destruction of the ecosystem under climate forcing that can be described as follows. Let us consider a population consisting of species with random parameters, and denote . We can assume, for example, that the parameters and in (3) and in (10) are normally distributed random variables. The equilibria are defined by roots of equation .

Let us consider system (29) for . Let be a steady state for this system, and we define a matrix with entries

 Mlj=∂Gl∂qj(Q1,Q2).

We introduce vectors and

 Ea(Q)(j)=col(a1jexp(−a11Q1−a12Q2),...,aNjexp(−aN1Q1−aN2Q2)).

Then we obtain

 Mkl=⟨b(k),Ea(Q)(l)⟩C,l,k∈{1,2}.

An Andronov-Hopf bifurcation occurs if the trace of the matrix changes its sign as the bifurcation parameter goes through a critical value and if the determinant of is positive at this critical value. Using the notation in (8) we obtain

 DetM=M11M22−M12M21, (31)
 TrM=⟨b(1),Ea(Q)(1)⟩C+⟨b(2),Ea(Q)(2)⟩C. (32)

These relations allow us to see connections between bifurcations, feedback, and diversity. First let us observe that components of the vectors are always positive. Note that if the climate influence is absent, then all the components of are negative, and it is clear that does not change its sign. Thus in this case the Andronov-Hopf bifurcations are absent. The same fact holds if all the climate-ecosystem feedbacks are negative. For purely positive or mixed feedbacks these bifurcations are possible under additional conditions. In order to find a biological meaning of these conditions, we define as the angles between the vectors and . We then have

 ϕlj(C)=⟨b(l),Ea(Q)(j)⟩C||b(l)||−1C||Ea(Q)(j)||−1C.

Then the condition reduces to

 ϕ11(C)ϕ22(C)>ϕ12(C)ϕ21(C). (33)

The condition implies that and have opposite signs. Then (33) means that and also have opposite signs. If all the species affect the climate in a similar manner (the coefficients have the same signs) then all the quantities have the same sign. Therefore, Andronov-Hopf bifurcations are impossible in this case.

We conclude that not only feedback positivity but also biodiversity and a complex ecosystem structure support complicated time periodic behavior. Moreover, all bifurcation conditions depend on the initial data . From a biological point of view, this means that bifurcation effects have a ”memory”, i.e., they depend on the choice of initial data.

## 6 Equilibria

The aim of this section is to show that the cases of negative (NF) and positive (PF) feedback between climate and ecosystem are markedly different. In the NF case, positive equilibria with many species can exist. In the PF case, such equilibria vanish for some critical feedback level; this can be interpreted as a mass extinction. We compute this critical level.

On the attractor structure, one can say more for the particular case of system (11) and (12), where we have a single resource, . We use eqs. (1), (2), where and thus , where . Let us set . These quantities are important characteristics of species. The species with smaller have a greater chance to survive. Moreover, using an analogue of (13) for the case a single resource , we assume that depends on as follows:

 S(x)=¯S+N∑k=1bkxk,

where are the coefficients.

### 6.1 Equation for equilibrium resource value

Moreover, for simplicity, let us set

 γij=γiδij,γi>0. (34)

In this case numerical simulations show that all trajectories tend to equilibria. As was pointed out by V. Kozlov, using the theory of decreasing operators and an assumption that increases in the resource and , one can prove this fact by analytic methods (a detailed analysis of this question will appear in future work, since the proof is quite involved). The resting points of systems (1) and (2) can be found as follows.

Setting in (1), we obtain . This gives the following nonlinear equation for :

 D(¯S−¯v)=G(v), (35)

where

 G(v)=N∑j=1ajγ−1j(cjajϕ(v)−bj)(ϕ(¯v)−ρj)+. (36)

We have obtained a complicated equation with non-smooth nonlinearities. An important characteristic of the solutions is , the number of positive involved in the sum in the right hand side of (39). The number can be interpreted as biodiversity.

Note that, for any , in the NF case a solution with always exists under the following condition:

 ϕ(S)>ρ0=minjρj. (37)

Indeed, observe that is a decreasing function of , while is increasing. The solution is given by an intersection of the curve and the right line , which exists if (37) holds.

Moreover, the same geometrical argument shows that the resource is an increasing function of . Therefore, in the case of negative feedback the biodiversity is larger (if a solution exists). However, for negative that are too large, the positive solution does not exist.

Consider a large ecosystem with random parameters . We suppose that and are selected randomly according to a distribution with probability density function , which is positive on some open interval .

Assertion. Consider the case of negative feedback ( ). If

 ϕ(¯S)>R0, (38)

then for any there exists a positive solution of eq. (35) with biodiversity such that as .

Proof. The existence of solutions is obvious from geometrical arguments (see remarks on the monotonicity of and above). To show that is large for , we observe that for any fixed and such that , the interval contains points , with as . For large we seek a solution of (35) in the form , where . Since (38) holds, such a solution exists. The number approximately equals the number for , and the assertion is proved.

In the PF case this assertion, in general, does not hold. Using the arguments from the proof, we note that all species die if the following relation holds:

 R0N∑j=1a2jγ−1jcj

This relation shows that mass extinction inevitably arises if the are sufficiently large.

Results on a numerical solution of equation (35) are discussed below. They confirm that mass extinctions are possible as the feedback magnitude increases.

### 6.2 Numerical results

In the general case equation (35) for equilibria can be resolved numerically for . We choose the coefficients in equation (35) as follows. The positive coefficients are random numbers subject to log-normal distributions. This means that are distributed normally, , where is the mean and is the deviation. The same distribution is taken for , with the parameters and .

We assume that the and are distributed normally, namely, and , where , and is the magnitude of the feedback level. The other parameters were taken as follows: and .

The results are shown by Fig.1. Comparison of the two plots shows that when the number of species increases, so does the likelihood of a sharp drop in species number as the climate changes and feedback processes grow stronger. These findings are consistent with analytical results. Biodiversity grows with the feedback parameter , until at some critical level we observe a mass extinction.

## 7 Conclusions

In this paper, a resource model for a system of many coexisting species is proposed. It is a generalization of the well known model in (Huisman and Weissing, 1999), takes into account species self-regulation and a dependence on the environment, and is the first model of an ecosystem with many species and feedback which couples climate and population dynamics. Such conceptual models describe a simple and easily understandable mechanism for resource competition. For the case of fixed parameters, a general assertion on attractor existence for this model is proved. One of the sufficient conditions for the existence of an attractor is that species self-regulation is stronger than species competition.

Climate-ecosystem feedbacks are an important problem in terms of uncertainty in predictions and modeling future climate change. The proposed model allows us not only to investigate climate-ecosystem feedbacks for large ecosystems, but also to show that coexistence of many species feeding on a few resources is possible. In the case of positive feedback in the ecosystem-climate interaction, the numerical results show a possibility of catastrophic bifurcations, when all (or almost all) species become extinct under the impact of climate warming. The ecosystem biodiversity increases with the magnitude of positive feedback , but at some critical level of feedback, a mass extinction occurs. For negative climate-ecosystem feedback we observe smaller biomass and biodiversity values, but we do not observe catastrophes. Note that in the contemporary world, human impact on the climate system can possibly lead to positive feedback in the above context.

To investigate more complicated situations, where complex dynamics may be possible, we have considered the case of just a few resources. We find asymptotic solutions for the case of a large resource turnover. This allows us to reduce this system to the Lotka-Volterra model, which is well studied. The existence of two sharply different regimes of ecosystem behavior is proven: the weak climate regime (WC), and the strong climate regime (SC). This behavior depends on a parameter that determines the intensity of ecosystem-climate interactions. Note that this analytical result is consistent with experimental data (Crampton et al., 2016), where it is shown that two distinct regimes of extinction dynamic are present in the major marine plankton group. Results in (Crampton et al., 2016) suggest that the dominant, primary controls on extinction were abiotic (environmental), which corresponds to the SC case.

In the SC case we do not observe complicated dynamical effects when the ecosystem – climate interactions involve only negative or only positive feedback loops. However, if the ecosystem – climate interaction involves terms of different signs, then there are possible Andronov-Hopf bifurcations, time periodic behavior for the case of two resources, and chaotic behavior for more than three resources. We conclude that not only feedback positivity, but also biodiversity and a complex ecosystem structure (when different species affect climate differently creating positive and negative feedback ecosystem – climate loops) support complicated temporal dynamics of the ecosystem.

For the case of a single resource the ecosystem equilibria can be described implicitly. We find these equilibria by a nonlinear equation for the equilibrium resource level. We show that, due to self-limitation effects, the system can support equilibria with a number of species sharing the same single resource.

## Acknowledgments

The authors are grateful to the anonymous reviewers and the guest editor for interesting remarks and comments, which significantly improved the paper. We are thankful to Prof. M. L. Zeeman (Bowdoin College) for very helpful comments on the systems we study here. We also express our gratitude to Prof. V. Kozlov (Linkoping University) for his help.

This study was funded by RFBR, through research projects No.16-34-00733 mol_a and No.16-31-60070 mol_a_dk. We gratefully acknowledge support from the Government of the Russian Federation through mega-grant 074-U01, as well as from the Division of Mathematical Sciences and the Division of Polar Programs at the U.S. National Science Foundation (NSF) through Grants DMS-0940249 and DMS-1413454. We are also grateful for support from the Office of Naval Research (ONR) through Grant N00014-13-10291. Finally, we would like to thank the NSF Math Climate Research Network (MCRN) as well for their support of this work.

## Appendix

We state here the proof of Lemma 1. The proof proceeds in the following steps.

Step 1. Positivity of the follows from the fact that the -th right hand side of system (4) is proportional to , thus, , where is a function.

Step 2. Let us prove that . Assume that this fact is violated. Then there exists an index and a time such that

 vj0(t0)=0,dvj0dt≤0,vj(t0)≥0,∀ j. (40)

Condition (7) entails the term equaling zero. Then we substitute these inequalities into the -th equation (12) and obtain a contradiction.

Step 3. Let us prove estimate (17). First let us suppose that assumption 1B is satisfied. Let . Let us estimate for large . Let be an index such that . According to (6) the are uniformly bounded by . Therefore within any open interval , where is fixed, one has

 dEdt≤ERi0,Ri0≤C+−γE(t), (41)

where due to assumption (14) on .

In the case 1A we note that

 N∑j=1γi0jxj≥γiiE−∑j≠i0|γi0j|xj≥κE,

and we have an inequality analogous to (41):

 dEdt≤ERi0,Ri0≤C+−κE(t). (42)

Note that the sequence of intervals is not bounded and these intervals cover all since, according to the Lemma, the solutions exist for all .

Inequality (41) implies that where is the solution to the Cauchy problem

 dXdt=X(C+−γ0X),X(0)=maxixi(0), (43)

where equals in the case 1B and in the case 1A. Let . If , then equation (43) shows that for all and (17) follows. If , then equation (43) shows that for all . By the change of variables we obtain that and thus

 d~Xdt=−γ0(X0+~X)~X≤−γ0X0~X,

which implies , and we obtain (17).

Step 4. Having (17), we can prove (18). Indeed, using the non-negativity of the and , one obtains

 dvjdt≤Dj(Sj(x(t))−vj).

Therefore,

 vj(t)=exp(−Djt)(vj(0)+∫t0Sj(x(s))exp(Djs)ds)

which yields

 vj(t)≤exp(−Djt)vj(0)+maxs∈[0,t]Sj(x(s)).

Here . These two last inequalities imply , which completes the proof.

## References

• Arhonditsis and Brett (2004) Arhonditsis, G.B. and Brett, M.T., 2004. Evaluation of the current state of mechanistic aquatic biogeochemical modeling. Mar. Ecol. Prog. Ser. 271, 13-26.
• Chapin et al. (2002) Chapin, F.S., Matson, P.A., Mooney H.A., 2002. Principles of Terrestrial Ecosystem Ecology. Springer Science and Business Media.
• Crampton et al. (2016) Crampton J.S, Roger A., Coopera R.A, Peter M., Sadler P.M, Foote M., 2016. Greenhouse-icehouse transition in the Late Ordovician marks a step change in extinction regime in the marine plankton, PNAS, 113, 1498-1503.
• Doney et al. (2012) Doney, S.C., Ruckelshaus, M., Duffy, J.E., Barry, J.P., Chan, F., English, C.A., Galindo, H.M., Grebmeier, J.M., Hollowed, A.B., Knowlton, N., Polovina, J., Rabalais, N.N., Sydeman, W.J., Talley, L.D., 2012. Climate change impacts on marine ecosystems. Ann. Rev. Mar. Science 4, 11-37.
• Field et al. (1998) Field, C.B., Behrenfeld, M.J., Randerson, J.T., Falkowski, P., 1998. Primary production of the biosphere: integrating terrestrial and oceanic components. Science 281, 237-240.
• Hirsch (1985) Hirsch, M., 1985. Systems of differential equations that are competitive or cooperative. II: Convergence almost everywhere, SIAM J. Math. Anal., 16. 423-439.
• Hofbauer and Sigmund (1988) Hofbauer, J., Sigmund, K., 1988. Evolutionary Games and Population Dynamics. Cambridge University Press, Cambridge.
• Hofbauer and Schreiber (2010) Hofbauer, J and Schreiber S. J., 2010. Robust permanence for interacting structured populations, Journal of Diff. Equations, 248. 1955-1971.
• Hsu et al. (1977) Hsu, S.B., Hubbell, S., Waltman, P., 1977. A mathematical theory for single-nutrient competition in continuous cultures of microorganisms. SIAM J. Appl. Math., 32, 366-383.
• Huisman and Weissing (1999) Huisman, J., Weissing, F.J., 1999. Biodiversity of plankton by species oscillations and chaos, Nature 402, 407-410.
• Hutchinson (1961) Hutchinson. G.E., 1961. The paradox of the plankton. Am. Nat., 95, 137-145.
• Irigoien et al. (2004) Irigoien, X., Huisman, J., Harris, R.P., 2004. Global biodiversity patterns of marine phytoplankton and zooplankton. Nature 429: 863-867.
• Kedra et al. (2015) Kedra, M., Moritz, C., Choy, E., David, C., Degen, R., Duerksen, S., Ellingsen, I., Gorska, B., Grebmeier, J., Kirievskaya, D., van Oevelen, D., Piwosz, K., Samuelsen, A., Weslawski, J., 2015. Status and trends in the structure of Arctic benthic food webs, Polar Res., 34, 23775.
• Kozlov and Vakulenko (2013) Kozlov, V., Vakulenko, S., 2013. On chaos in Lotka-Volterra systems: an analytical approach. Nonlinearity 26 , 2299-2314.
• Legrand et al. (2003) Legrand, C., Rengefors, K., Fistarol, G.O., Graneli, E., 2003. Allelopathy in phytoplankton : Biochemical, ecological and evolutionary aspects. Phycologia 42, 406-419.
• Richardson et al. (2000) Richardson, T.L., Gibson, C.E., Heaney, S.I., 2000. Temperature, growth and seasonal succession of phytoplankton in Lake Baikal, Siberia. Freshw. Biol., 44, 431-440.
• Roy and Chattopadhyay (2007) Roy, S., Chattopadhyay, J., 2007. Towards a resolution of ’the paradox of the plankton’: A brief overview of the proposed mechanisms. Ecol. Complex. 4, 26-33
• Ruelle (1989) Ruelle, D., 1989. Elements of differential dynamics and bifurcation theory. Academic Press, Boston.
• Sekerci and Petrovskii (2015) Sekerci, Y., Petrovskii, S., 2015. Mathematical modelling of plankton-oxygen dynamics under the climate change. Bull. Math. Biol. 77, 2325-2353.
• Selkoe et al. (2015) Selkoe, K.A., Blenckner, T., Caldwell, M.R., Crowder, L., Erickson, A., Essington,T.E. et al. 2015. Principles for managing marine ecosystems prone to tipping points. Ecosystem Health and Sustainability 1, 1–18.
• Smith (1981) Smith, H.L., 1981. Competitive Coexistence in an Oscillating Chemostat. SIAM J. Appl. Math. 40, 498-522.
• Smith and Thieme (1991) Smith, H.L., Thieme, H.R., 1991. Convergence for strongly order-preserving semiflows, SIAM J. Math. Anal. 22, 1081-1101.
• Takeuchi (1996) Takeuchi, Y., 1996. Global Dynamical Properties of Lotka-Volterra Systems. World Scientific, Singapore.
• Tilman (1977) Tilman, D., 1977. Resource competition between platonic algae: an experimental and theoretical approach. Ecology 58, 338-348.
• Toselandet al. (2013) Toseland, A., S.J., Daines, Clark, J.R., Kirkham, A., Strauss, J., Uhlig, C., Lenton, T.M., Valentin, K., Pearson, G.A., Moulton, V., Mock, T., 2013. The impact of temperature on marine phytoplankton resource allocation and metabolism. Nat. Clim. Chang. 3, 979-984.
• Travers et al. (2007) Travers, M., Shin, Y.J., Jennings, S., Cury, P., 2007. Towards end-to-end models for investigating the effects of climate and fishing in marine ecosystems. Prog. Oceanogr. 75, 751-770.
• Ulanowicz and Kemp (1979) Ulanowicz, R.E. Kemp, W.M.,1979. Prediction, Chaos and ecological perspective. In: Efraim Halfon, (ed.), Theoretical Systems Ecology. Academic Press, NY.
• Vakulenko (2013) Vakulenko, S., 2013. Complexity and Evolution of Dissipative Systems. An Analytical Approach. Berlin, Boston: De Gruyter.
• Zeeman (1998) Van den Driessche P. and Zeeman, M. L., 1998, Three dimensional competetive Lotka-Volterra systems with no periodic orbits, Siam J. Appl. Math., 58, 227-234.
• Volterra (1931) Volterra V., 1931. Lecons sur la théorie mathématique de la lutte pour la vie. Paris: Gauthier-Villars. Reissued 1990, Gabay, J., ed.
• Walther (2010) Walther, G.R., 2010. Community and ecosystem responses to recent climate change. Phil. Trans. R. Soc. B 365, 2019-2024.
• Zeeman (1995) Zeeman, M. L., 1995. Extinction in Competitive Lotka-Volterra Systems, Proceedings of the AMS, 123, 87-96.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters