Open Research Online Biome stability and fragmentation under critical environmental temperature change

Critical phenomena in the climate system can cause drastic changes in the state of planetary ecosystems as well the entire biosphere. There also are mechanisms through which the biosphere can make an eﬀect on climate. In this manuscript, we study the nonlinear dynamics of the interaction of the climate system with the biosphere by linking an energy balance climate model to diﬀerent species competition models. We develop an asymptotic approach to these models and investigate how migration strengthens biome stability and biodiversity. Moreover, we derive relations describing biome boundary shifts under global warming (or cooling) and check those relations against paleo data on plant biome location. Finally, the models demonstrate that critical rates of changes in the environmental temperature dynamics may shift biome stability. We show that there


Introduction
We aim to propose a model describing nonlinear interaction between planetary climate and a biosphere consisting of many ecosystems integrated into biomes where many species compete for a few of the same resources. The existence and stability of such food webs on a global scale is a key problem in ecology [1,2]. It becomes even more pressing to understand in the context of extinctions and mass extinctions in such systems under climate changes [3]. Resource competition models link the population dynamics of competing species with the dynamics of the resources. As it was mentioned in [4] an attractive feature of resource competition models is that they use the biological traits of species to predict the time evolution of competition. In fact, many rigorous results [5][6][7] show that, in a general situation, a single species survives, while to obtain the coexistence of many species one needs very special assumptions made to the species parameters (mortalities and resource consumption rates).
Note that ecosystems present a class of dynamical systems exhibiting complex dynamics, bifurcations, and strange attractors (see, for example, [8] and [9] for an analytical proof). On the other hand, climate also presents a complicated dynamical system [10]. Recent observations have shown that climate change may be a leading factor influencing ecosystem behavior [11]. To simplify climate modeling, Energy Balance Models (EBMs) were proposed [12][13][14][15] that allow us to describe current climate change and paleoclimate. The stochastic effects in such models are investigated in [16]. Other interesting phenomena in these models arise as a result of multistability, see [17,18].
Here, we develop a mathematical approach to investigate effects in nonlinear climate-biosphere interaction by a model that takes into account spatial heterogeneity, species migration, and self-limitation, and as well simple climate dynamics. We describe the biosphere with the help of models like the paradigmatical MacArthur model, and climate by the classical EBM with thermal conductivity. The coupling of two systems goes through planetary albedo [19]. We find that migration increases the growth rate function by a term that is proportional to τ −1 m , where τ m is a time needed for a species to migrate through the whole species habitat. Further, ecosystem multistationarity may lead to a fast biome formation and heterogeneity under climate variations when migration effects are weak.
We obtain equations describing a shift of the ecosystem boundaries under environmental temperature changes. There exist two regimes, slow and fast climate temperature variations. For heterogeneous ecosystems, slow climate variations (relatively uniform over habitats) can weakly change species habitats. On the other hand, too great a heterogeneity in the ecosystem can cause instability under random and periodic climate oscillations so one can envision that there is an optimal level of biosphere heterogeneity. Finally, we confirm the GAIA hypothesis that climate and biosphere can evolve together in a convergent manner: the biosphere has become increasingly diversified and spatially heterogeneous due to climate change that, in turn, has increased the stability of both climate and biosphere.
The paper is organized as follows. In the next section, we describe different classical models of ecosystems where competing species share resources and extended models that take into account climatic factors, migration, and selflimitation effects. In Section 3 it is shown that for large turnover rates the long term behavior of all these models can be described by asymptotic solutions. These solutions show that all these models can be reduced to a special multispecies Lotka-Volterra (LV) model that is well studied [9,[20][21][22]. In addition, we restate some results on the LV system, where we show how to realize prescribed polynomial dynamics by an LV system using a nonlinear oscillator concept. Section 4 states an averaging method that allows us to estimate biodiversity for systems without migration. To obtain these estimates we need no information on the attractor structure, they hold for any oscillating large time dynamics. Further, in Section 5 we extend these estimates on spatially extended systems with migration effects. It allows us to obtain the following important result: the effects of migration can be described by adding a term to the growth function that is proportional to the time it takes for the species to migrate throughout the entire ecosystem. At last, in Sections 6, 7, and 7.1, we describe how environmental temperature warming affects species habitats and apply analytical results to interpret data on evolution of the vegetation biomes of Europe from 16 4 million years ago (Neogene period). We show that there arises two cases: F and S. The case F corresponds to a fast climate variation with respect to species adaptation rate and S a slow climate variation. In the case of F the species have no time to adapt, and the habitat boundaries move at the same speed as the isotherms. In the second case, the speed of habitat boundary movement may be much less than the speed of isotherm change. The analysis of paleodata by our approach (both cases F and S) shows that biome evolution under cooling was rather fast than slow.
The main result is that a complicated spatially heterogeneous biosphere split into a number of small different ecosystems shows high stability with respect to slow climate variations but that the biosphere may be volatile with respect to rapid climate changes. So, a rate-induced (R) tipping point [23] may appear.

Models
In this section, we describe some conceptual climate and biosphere models and connections between them.

Energy balance climate models
The energy balance models are basic for such climate applications as temperature warming and paleoclimate [15]. In the simplest case when we consider the average temperature of the Earth's surface, the main equation has the form [13] : where λ is thermal inertia, T is the averaged surface temperature, t is time, and A is the average albedo of the earth's surface. On the right-hand side, the first term is the outgoing emission and the second term represents the incoming Sun's radiation. Generally, incoming radiation to the Earths surface is modified by a parameter, µ 0 , to allow for variations in the irradiance per unit area, I 0 (the solar constant), or for long-term variations of the planetary orbit. The outgoing emission depends on the fourth power of temperature, the effective emissivity e and a Stefan−Boltzmann constant σ.
This model can be coupled with biosphere's dynamics as follows. The complete averaged albedo A can depend on the biosphere state. For simplicity, we consider ecosystems in which species are competing for several resources. We will consider a more complicated model than in [19], which is spatially extended and takes into account species migration.
To take into account the dependence of T on latitude, we use the classical EBM model [12][13][14] and following [15] we insert a term describing a heat conduction into (1) and extend (1) equation to a PDE adding a spatial dimension.
Let y be a latitude, y ∈ (−π/2, π/2). Then the corresponding extended equation where D is the thermal conductivity coefficient, and ∆ B is the Laplace-Beltrami operator, which in the one-dimensional case, where we take into account the latitude only, has the form: where y is the latitude, which ranges from −π/2 (the South Pole) via y = 0 (the equator) to π/2 (the North Pole) (if y ∈ (0, π), where 0 is the North pole, then sin y should stand instead of cos φ). The boundary conditions at the poles have the form ∂T (y) ∂y = 0, y = ±π/2.

Standard ecosystem model
We consider the following system of equations that extends the known model [8,24], taking into account diffusion in space (species migration) along latitude (representing species response to global temperature change): Here u = (u 1 , u 2 , . . . , u N ) is the vector of species abundances, u i = u i (x, t), is a vector of resource amounts. The domain Ω ⊂ R d , where d = 1, 2 is bounded and has a smooth boundary ∂Ω. In eqs.
(4), (5) unknown v k is the resource of k-th type consumed by all ecosystem species, µ i is the species mortalities, D k > 0 are resource turnover rates, S k is the supply of the resource v k , D k v k (S k − v k ) is the supply rate for k-th resource and c ik > 0 is the content of k-th resource in the i-th species. We suppose that resource supplies S k are slow smooth functions of x, t: where j > 0 are small parameters. Such a slow dependence may describe variations of resource supplies in time and space that can appear as a result of climate and environmental changes. For example, one can suppose that for plant ecosystems S k depend on the temperature T and precipitation level P , mean annual values of these parameters can slowly evolve in time: T = T (τ ) and P = P (τ ). So, such ecological factors can be included in our model.
Conditions (7) and (8) can be interpreted as a generalization of the von Liebig law, where where a i and b ij are positive coefficients, i = 1, . . . , N . The coefficient a i is the maximal level of the resource consumption rate by i-th species and coefficients b ij , where j = 1, . . . , n, defines the sharpness of the consumption curve φ i (v).
Note that for a biome of relatively small area we can neglect curvature effects and replace the operator ∆ B by the usual Laplace operator ∆. Our results can be extended to the general case in a quite straightforward way.

MacArthur consumer competition model
The MacArthur model is another paradigmatical model of ecology for systems with sharing resources, discussed recently in detail in [25]. It is similar to the standard model but here the growth rate functions and resource turnover coefficients are linear. In our notation the spatially extended MacArthur model reads: where c ik is the successful encounter rate of species i searching for resource k, r i is the maintenance cost or threshold consumption level for growth, w k is the per capita weight or nutritional value of each resource, e i is the quantity of nutritional value required for the reproduction of a given species.

Contemporary niche model and models with self-limitation
The other variant of the MacArthur model is proposed in [25], where the turnover rates do not include the factors v k : We will show below that these models without migration encounter some difficulties when we try to explain high biodiversity in some ecosystems. We can improve the situation by incorporating self-limitation effects. For example, the previous model can be modified as follows:

Reduction to special Lotka-Volterra system
Here, in the case of a large turnover D k , we reduce all the models to special Lotka-Volterra (LV) systems that are well-studied in [9,20]. First, we consider the resource competition model. Assume that D k ∝ D, where D is a large parameter that allows us to proceed a time-scale analysis. We setD k = D k /D and, to balance the first term on the right-hand side of eq. (5), we put in that whereṽ k are new unknowns. We substitute the relations (20) for v k into our system (4) and (5). Then, by matching principal terms and neglecting terms which have the order O(D −2 ), one obtains where and 2 , τ are defined by relation (6).
A given u, system of eqs. (22) forṽ k has an asymptotic solutioñ Removing terms << D −1 and << 2 and substituting (20) into eq. (21) we see that system (22) takes the form where and where A, B are matrices of size N × n and n × N , respectively. Note that by the change u = Dû we can replace the large factors D k >> 1 by normalized turnoversD k = D k /D that gives a LV system without small parameters, where Analogous results can be obtained for the MacArthur and contemporary models.
Let us consider, for example, the MacArthur model. In fact, we note that for D >> 1 the resource v k is close to S k , thus we can replace the coefficients In a similar way, we then obtain that this system (14)- (15) reduces to the LV system (24), where the factorization (26) takes the place of the parameters So, we conclude that, in a sense, all models are asymptotically equivalent to LV systems for large supply turnovers. If we take into account self-limitation effects then

Biodiversity in models without migration
In this section, we present estimates of ecosystem biomass and biodiversity.
These estimates are based on the concept of permanency and the averaging method. We need no precise information on the kind of dynamics, i.e., estimates hold for chaotic, oscillating, or convergent dynamics. In fact, in Appendix 1 it is shown that the dynamics of ecosystems can be complicated.
Recall the fundamental concept of permanency [20]. A system is permanent (i.e., ecologically stable) if there exist positive constants δ and C such that This implies δ < u i (t) < C for all sufficiently large t. Let f (t) be a function such that f (t) > 0 for all t and for some C one has f (t) < C for all t > 0. Let us denote by Suppose that our ecosystem is permanent and the resources do not vanish: for some δ R > 0. Note moreover that the equations (15) for We consider the MacArthur model with self-limitation coefficients γ i > 0.
We divide the i-th equation (18) with d i = 0 by u i and proceed averaging for the resulting equations. Then the permanency condition implies that Therefore, using (29) and (30) one obtains from (18) Repeating the same trick for v equations one has These averaging equations for averages of u i and v k lead to interesting consequences. Let first γ i = 0, i.e., there are no self-limitations. Then relations (31) produce a linear algebraic system for means v k : For generic w k , c ik , this system has a solution only if N ≤ n, i.e., the species number is not more than the resource number. This means that the competitive exclusion principle works. However, many real ecosystems consist of many species sharing a few resources. To avoid this paradox, we can take into account migration (see below) or self-limitation. Let us consider the second case. Then for γ i > 0 eqs. (31) and (32) lead to the following system of nonlinear equations for averaged resources v k ∈ (0, S k ): where z + = max{0, z}.
Similar equations can be obtained for the resource competition model. For Systems (34) and (35) and in the case (35) the quantity N i is the number of indices i such that

Migration and biodiversity
The averaging method and results of the previous section allow us to estimate the effects of migration (diffusion) species on biodiversity. Consider i-th equation (4) under conditions that the system (4), (5) is permanently uniform in space, i.e. there exist positive constants δ and C such that Let us divide the both sides of this equation on u i and integrate over all x ∈ Ω and t ∈ (0, T ). The averages over x and t we denote by and let us set Then we obtain Let us consider the first term. Using the boundary conditions (12) we observe Therefore, (37) is equivalent to For large D k we can replace φ(v) with φ(S). The first term admits the following rough estimate by a scale analysis. Let us consider, say, an ecosystem where u i (x, y) depends on latitude y, which changes from L 0 to L 1 . Then the average of (∇u i /u i ) 2 can be roughly estimated as (L 1 − L 0 ) −2 and thus can be evaluated as a time needed for the species migration across the entire ecosystem (along the latitude direction). So, we obtain that when we take into account the migration effects, the main eq. (35) should be modified as follows: where the correction τ i is equal to The quantity τ i can be interpreted as a migration rate of the species across the entire ecosystem. We observe that for systems with self-limitation the migration always increases biodiversity and this effect can be described by an additional term in the growth function.

How climate affects ecosystem boundaries
In this section, we use the results of previous sections and the niche theory approach. Using eqs. (39) we can describe how climate affects biodiversity and rate tipping point effects. The main idea can be outlined as follows. We consider two modes of interaction. If climate changes at the same rate as the adaptation of species, or more slowly, then the movement of biome boundaries is rather slow. If the climate (for example, the average annual temperature) changes too quickly, then even a large enough biome can disappear for several million years.
Thus, one can formulate a hypothesis that biosphere-climate interaction can lead to an R-tipping point.
To simplify the statement, let us consider a simple model, where the main resource S depends on is the temperature T a averaged over a year and over the ecosystem, but in principle we can take into account precipitation, light, nutrients, and others. Note that, temperature is usually viewed as an environmental factor rather than a resource. We can easily extend our basic model to include such a parameter, which is important for defining an ecological niche. In addition, in some works, temperature is considered as a resource (see, for example, the work [26]).
We set, to simplify, S = S 1 = T a . We model the growth function by a Gaussian approximation where a i , σ i > 0 are parameters.
This growth function corresponds to a species for which the optimal temperature equals T opt,i . Suppose that parameters T opt,i are distributed according to a probabilistic density function (pdf), for example, a Gaussian N(E T , σ T ) one with parameters E T , σ T or a log-normal density. Then Eq. (36) reduces to where r i = c −1 i µ i a i . By adjusting parameters T opt,i , r i , E T , σ T we can obtain a complex distribution of biodiversity in space. In the next section, we use the obtained results to analyze paleodata of the plant biomes.
6.1. Slow movement of species habitat boundaries under climate variations.

General relations
Let us consider the ecosystem boundary motion under a climate impact.
The main idea to define the ecosystem boundary is as follows. Consider, say, a plant ecosystem (for example, forest or tundra). In these ecosystems we have different sets of species. Let us consider a species growing in tundra, its local growth function at the point (x, y), where x is the longitude and y is the latitude, is defined by the expression One can expect what are the most favorable conditions for species survival available in the species habitat center whereas at the habitat edge the growth function vanishes.
Let us define optimal T opt,i that gives the maximal species abundance. Using results of Sections 4 and 5 we obtain that these optimal values T opt,i are defined by the relation Then one has We insert (42) into the last equations that gives finally where we use the following notation for averages: The relation (45) defines the center of the Hutchinson box (see [27]). If an interval of evolution goes on for sufficiently long enough and the dependency T a (x, y, t) on time t is slow, one could expect that almost all the species in the ecosystem have growth functions with optimal parameters defined by (45). We will see below that this conclusion is consistent with certain paleodata.
One can think that the boundary is located at the points (x, y) where Φ i (x, y) = r i . Then we obtain the following equation for the habitat boundary: where b i are constants which can be found by empiric data. If we suppose that only the temperature plays a role in the system, and there is no species adaptation, we obtain the simplest equation i.e. the habitat boundaries coincide with isotherms. This result permits us to check a hypothesis that cooling or warming are the main factors determining the evolution of the species habitats.
Note that there are two different cases: (F), when we have a fast climate variation with respect to species adaptation rate and the case (S) corresponds to a slow climate change. In the case F the constantsb i do not change because the species have no time to adapt, and then the motions of habitat boundaries go with the same speed that the isotherm moves. In the second case, the speed of habitat boundary motion may be less than the isotherm speed.

Fast and slow climate changes: a comparison
First, we consider an approximation of the temperature profile on the latitude y.

Temperature profile
We follow the manuscript [15] (see Chapter 5 in there). By linearization of the term with T 4 , we can reduce eq. (2) to a linear equation. This can be resolved by the Legendre polynomials that give an approximation where T 0 > 0 and T 1 are constants. When climate changes these constants very slowly evolve. We assume that they depend on the slow time τ . Further, we consider a more general approximation where α > 0 corresponds to a globally uniform warming and α < 0 describes a globally uniform cooling. In the coming subsection we apply the relation (49) to compare the cases F and S.

Rapid climate variations
Suppose the species area occupies the domain, where the latitude y ranges in the interval (y − (τ ), y + (τ )) and T a depends on y, t only: T a = T a (y, t). Then for the right habitat boundary y + (τ ) one has ∂T a (y + , t) ∂t where V b = dy + /dτ is the velocity of the boundary. Using (49) we see that eq.
(50) reduces to where ψ (y) = dψ dy . This quantity may not be small for a small habitat size L = y + − y − . From the last equation we obtain a rough estimate for the biome boundary shift ∆y under the temperature change: where ∆T clim is a change of temperature under the climate variation and the constant const ≈ 2, ∆T biome is a change of the local temperature along the latitude within the entire biome (from the biome center to boundaries).

Slow climate variations
Let us consider the case S, where the analysis is more sophisticated. Then instead of eq. (47) we have where according to (45) In the general case for V b one obtains a complicated equation, but if we consider an approximation (49) and suppose that σ i is large enough, then the boundary shift velocity V b does not depend on the climate change rate α: So, one can expect that in the case S the biome is stable under climate change impact and biodiversity increases the stability.

Handling real data
In this section, we use real data to apply the approach stated in previous subsections. As an example, we use the vegetation biomes of Europe from 16 4 million years ago (Neogene period). This time interval is characterised by warmer-than-present climate conditions, but with a general cooling trend towards present day [28,29]. Against this background global cooling trend, the type of vegetation biome across the majority of Europe remains relatively unchanged with a dominance of the warm-temperate mixed forest [29][30][31][32][33]. Progressive loss of species over this time interval is known, especially amongst plants requiring permanently humid and warm conditions [34].

Paleo-biome data was reconstructed from paleobotanical sites from across
Europe covering the Langhian Zanclean stages (16 3.6 million years ago).
These come from [29][30][31][32]  Co-existence Approach [35], although some use probability density functions [36]. The nearest living relative techniques use the climate tolerance of a fossil plants extant relatives to reconstruct past climates. The overlapping envelope of all fossil taxas climate tolerance is then taken as the reconstruction for a paleobotanical assemblage [35,36].
In Section 6, we obtained a method to compute the complex distribution of biodiversity in space. Here, we use the obtained data to find the biodiversity distribution of warm-temperate mixed forests (wtf), which is shown in Fig. 1.

Likelihood approach
Let us consider a set of species with random parameters T i distributed according to a probability density function (pdf) ρ(T ). For example, we can take the normal ρ, where T i are distributed according to Norm(ET, s T ), where Norm(a, s) denotes the normal law with the mean a and the standard deviation s (another variant is log-normal distribution, where log T i ∈ Norm(a, s)).
We define niche optimal parameters by the maximal likelihood principle [37] (this standard procedure is outlined in Appendix 3).

Data analysis and comparison with the theory
The first result is that the theoretical relation (45), that defines the Hutchinson niche center, is confirmed by the likelihood estimate described above. It can be illustrated by Fig. 2 Our analysis proceeds as follows: we divide the entire time interval into two equal intervals (i) and (ii). We select several points corresponding to time points from the first interval with the highest latitudes (say, 10 points) and find the average latitude. We consider the average of these latitudes as the boundary of the biome in the first time interval. We do the same for the second time interval. Calculating the difference between these averages, and dividing that difference by the time lag between the points, we get a rough estimate of the speed of movement of the biome boundaries.
As mentioned at beginning of Section 7, paleo-biome data are taken from [29][30][31][32] see also references cited in these works. We split that time period 16  The main conclusions of these computations are as follows: i. slow climate variations are not so dangerous for the biosphere; ii. the biodiversity and heterogeneity help survival because small areas are warmed up more evenly.
So, biosphere heterogeneity can help mitigate uniform temperature warming. However, too strong heterogeneity makes ecosystems unstable against fast periodical and random fluctuations [39]. In the next section, we show describe a possible mechanism for heterogeneity formation, which can be emerged via a fragmentation process. This mechanism is a result of the existence of competing We see that the most probable values obtained by experimental data via the maximum likelihood estimate are close to zero that is consistent with estimate (45). ecodynamical attractors.

A dynamical mechanism of biosphere fragmentation.
In the previous section, a mechanism of ecosystem extension (shrinking) under a climate effect was investigated. This mechanism describes slow noncatastrophic changes that reduce the evolution of biome boundaries. There is another mechanism connected with bifurcations producing fast changes. This second mechanism leads to a relatively fast biome formation.
We use an asymptotical theory developed in [40] assuming that migration effects are weak: d i << 1. This approach works for all reaction-diffusion systems with small diffusion coefficients if the global attractor of the corresponding shorted system (without diffusion terms) consists of a few local attractors, i.e., there is a multistationarity. We have seen above that multistationarity is possible even for a single resource (see remark (iii) to Theorem 9.1).
Let us consider the system (24). Let first d i = 0, i.e., diffusion is absent.
Suppose that this shortened LV system has local attractors A l , l = 1, ..., m A , m A > 1. Then one expects that there exists separatrice sets S lk that separate attraction basins B l of attractors A l . We suppose that these sets (which are attractor basin boundaries) are unions of N −1 dimensional submanifolds in the u-phase space R N (but here it is necessary to note that in the case of competing chaotic attractors there is a possibility that our assumption is violated, see [41,42]). Now we consider the mapŪ : x →ū(x) defined on the domain Ω by the initial data for system (24). Let us define the set Ω lk as preimages of S lk under the mapŪ : Basic theorems of differential topology assert (see [43]) that if this map is, in a sense, generic, and S lk is a manifold having codimension 1 in R N then the set Ω lk also is a manifold of codimension 1 in the domain Ω. This means that for dimΩ = 1 that set is a point and for dimΩ = 2 it is a curve separating two subdomains Ω k and Ω l , which are preimages of attraction basins for competing local attractors A k and A l , respectively. As a simple example, we can consider the well studied scalar Ginzburg-Landau equation [44] u t = d∆u + u − u 3 , Then we obtain that for large t u(x, t) approaches to A l if the initial datā u(x) ∈ B l and u(x, t) approaches to A k ifū(x) ∈ B k . Therefore, the sets Ω lk can be interpreted as boundaries between different ecosystems, which can emerge within a short time period τ b , which, when species migration is weak and coefficients d i are small, can be estimated as follows: This logarithmic estimate is obtained in [45], see also Appendix 2. The dynamics within ecosystems can be quite different. This mechanism of biosphere fragmentation into biomes is illustrated by Fig 3.

Conclusions
In this paper, we find a possibility of chaos and oscillations as a result of global climate warming. We study biome formation and fragmentation. We conclude that small biomes are disadvantageous because they are too unstable under fluctuations. For the ecosystems, we find that the rate of reversal occurs due to climate change.
In summary, the main results are: 1. It is shown that three main models of ecosystems where species compete Finally, we can formulate a hypothesis that climate and the biosphere evolved together in a convergent manner: the biosphere has become increasingly diversified and spatially heterogeneous due to climate change, which increased the stability of both climate and the biosphere.
Understanding biome fragmentation is important in the research of paleodata. To further complicate our understanding, the majority of terrestrial fossil sites have a high degree of uncertainty in their dating resolution (e.g. [46]).
Attempts have been made to get around this via likely biomes in warm-wet and cold-arid orbital phases [31] and through high-resolution studies [47]. These are not feasible for all paleobotanical sites and therefore restrict our geographical understanding of the response of biomes to climate change. We expect that the proposed approach can help to reconstruct paleobiomes.

Acknowledgments
The reported study was funded by the the Russian We are grateful to anonymous Referees for useful and interesting comments.

Appendix 1. Dynamics of reduced LV systems without migration
In Section 3, it is shown that all the models can be reduced to a special Lotka-Volterra (LV) system. This system has interesting properties.

Lotka-Volterra system with many resources
We remove diffusion and self-limitation terms, and ignore the dependence on slow variables X, τ . Then LV system (24) takes the form in which we have N species with populations u i for i = 1 to N . We consider this system in the positive cone R N > = {u = (u 1 , ..., u N ) : u i > 0}. Notice that this cone is invariant under dynamics (57). Below we assume that initial data for (57) always lie in this cone: We also assumer for certain µ k (S), k = 1, . . . , n. This condition is necessary to provide coexistence of many species. For a biological interpretation of this assumption see [9].
Note that this condition holds if φ i are linear functions of S k .
Under the above assumptions, eqs. (57) can be represented as a model with n resources: where This means that all the growth coefficients S i depend on resources R j (u) that are linear functions of u. The coupling with climate energy balance models can go via dependence of A and B on temperature.

Dynamical complexity
In [9] the following Theorem is proved: Theorem 9.1. Let F (l) , l = 1, ..., p, be C 1 -vector fields on B n directed inward on ∂B n and having compact invariant hyperbolic sets I (l) . Then there exist a positive integer N , µ ∈ R n , matrices A, B of sizes N × n and n × N respectively and C (l) ∈ R N > such that system (60) has compact invariant sets K(C (l) ) ⊂ Q n (C (l) ), which are homeomorphic to I (l) . These sets are hyperbolic for the flow S t LV | K(C (l) ) and, moreover, the flows S t LV | K(C (l) ) and S t F (l) | I (l) are orbitally topologically equivalent.
In other words, if a finite family of hyperbolic dynamics is given, a sufficiently large Lotka-Volterra model with appropriate parameters can generate this family by a variation of initial data. These hyperbolic dynamics may be chaotic.

Remarks.
(i) For C sufficiently close to C (l) the manifold Q n (C) also contains a compact hyperbolic invariant set K(C) homeomorphic to I (l) . (iv) The main idea of the proof of Theorem 9.1 is to use the Volterra variables The sums of exponents in the right hand side can approximate any smooth functions in compact domain and moreover by (62) we can reproduce all polynomial systems.
Let us consider an example of how to realize a prescribed polynomial dynamics by an LV system. The case of the Lorenz system is considered in [9] where N = 11 species are used to represent the Lorenz equations. We consider another fundamental system We consider this system in the domains D C = {x : |x 1 | < C 1 , |x 2 | < C 2 }, where C 1 , C 2 > 0 are constants. This system describes a nonlinear oscillator, so, we observe here periodical oscillations and also homoclinic solutions x 1 = ± tanh α(t − t 0 ) , where α = 1/ √ 2, which correspond to oscillations with infinite period. Existence of such homoclinic solutions shows that small perturbations of system (63) can exhibit a chaotic behaviour (heteroclinic chaos) [48]. The most natural way is to make a variable change in eqs. (63) by x i = exp(b i q i ) − C i that immediately leads to a vector field of the form (62). In fact, then we obtain that eqs. (63) take the form It is clear that we can adjust such N , A il , µ l and B lj that system (63) with n = 2 coincides with (64) and (65). The minimal possible species number is N = 5, and entries A jl take values −b 2 , −b 1 , b 2 , 2b 1 , 3b 1 .

Appendix 2.
Following [45], let us outline the derivation of the estimate (56). For simplicity let us consider a single reaction diffusion equation where x ∈ Ω, d << 1,ū is a smooth function and we set the zero Neumann condition at the boundary of Ω. we assume that the pure reaction dynamics, defined by the Cauchy problem for v = v(x, t) (where space variable x is a parameter) v generates different competing local attractors, say, rest points v eq = u ± . We note that if we have at least two local attractors for (67) and initial dataū(x) lie in different attraction basins for different value x, then |v x (x, t)| increases as exp(a 0 t), where a 0 > 0 is a constant corresponding to the maximal positive Lyapunov exponent for a saddle set S corresponding to two competing attractors (an escape rate). Therefore, the term d∆v grows as d exp(2a 0 t). This term becomes essential for times t of the order τ b ∝ ln d. In fact, for t > τ b we must take into account in (66) not only the reaction part but also the diffusion part. The quantity τ b is a characteristic time needed to create a spatially heterogeneous structure or a traveling wave connecting states v = u + and v = u − .
To justify these estimates, we put u = v +ũ in (66) and by (66), (67) one obtains the following estimate for a unknown correctionũ: where C, a are positive constants depending on f [45]. From this last inequality we see that the norm of perturbation U (t) remains small on the interval O(ln d).
Note that if there are no competing attractors for (67), then we can take a = 0 and under some conditions we obtain uniform in t estimate for U (t).
Moreover, these arguments can be extended on the case of reaction-diffusion systems [45].

Appendix 3. MLE
Given observed data the maximum likelihood estimation (MLE) finds the parameters of an assumed probability distribution. To obtain those parameters we maximize the likelihood function L so that, under our statistical model, the observed data is most probable. Our list of parameters is P = {φ c , ET, s T }.
Here ψ(x, y) is the probability that a species located at (x, y) survives in the Hutchinson box (niche) with the parameters P (see about the Hutchinson box concept in [27]). For the points (x, y) belonging to warm temperate forest biome this probability is ψ(x, y, P) = H(x,y,P) For the points (x, y), which are not belong to warm temperate forest biome, this probability is ψ(x, y, P) = 1 − H(x,y,P) ρ(T )dT.
Practically we can compute the integrals (69) and (71) by the Monte-Carlo method sampling T i by the pdf ρ.