Enabling stable coexistence by modifying the environment

In this work coexistence is studied using a model based on two classical population models: the quasispecies of Eigen [1] and the daisyworld presented by Watson and Lovelock [2]. It is assumed that species are able to modify the environment. We show that this ability enables the coexistence between, at most, two species in equilibrium. Given an initial population, the problem arises as to determine which of the many equilibrium populations, i.e. extinction, only one species or coexistence of two species, will be reached as a function of the species characteristics, specifically their capacity to modify the environment and the optimal growth rate. These results are obtained under general assumptions, which broadens its applicability to other fields as evolutionary biology and social sciences.


1
The question of whether an initial population formed by different species will persist on 2 time is of great relevance in Ecology and Evolutionary Biology [3][4][5]. However, to get 3 conclusive results seems to be elusive because of the so many factors that are involved 4 in the dynamic behavior of the population [6]. Besides, the individual properties of the 5 species have to be confronted with initial and boundary conditions. 6 From a dynamical point of view, this question is translated to which asymptotic 7 equilibrium the population will attain when multiple stable equilibria exist. Classical 8 qualitative analysis, using local properties of equilibria, is not enough to solve this 9 question [7,8]. Finding its solution requires a global approach that considers the 10 properties of all the species as well as the possible interactions among them and the 11 influence of the environment. Unfortunately, despite the large number of papers devoted 12 to study the relationship between biodiversity and stability in ecosystems [3,5], there There are several examples where the species ability to modify their environment are 23 specially relevant. Tumor cells are able to modify the surrounding tissues by segregating 24 chemical substances that increase the fitness of malign cells [13]. Soil modification by 25 microbial communities is being reported as one of the main driving factors of these 26 ecosystems [14]. Finally, climate change caused anthropogenically, i.e. induced by 27 human actuations, constitutes a dramatic example of the effect of modifying the 28 temperature by the species that inhabit the planet [15]. 29 In this paper we study a population of replicators that are able to modify their 30 environment [16]. The environment is here described by an unique scalar variable, that 31 we will call temperature. We assume that species fitness depend on this temperature 32 that changes over time as a function of the population distribution. Maximum fitness 33 for each species is reached at a particular optimal temperature value, decreasing when 34 temperature departs from this value. Species influence on temperature can be either 35 positive or negative and it is not coupled with their optimal temperature. Besides, total 36 population is bounded by a carrying capacity of the system. As a consequence, the 37 initial population undergoes a selective proccess that ends in an equilibrium population 38 formed by the survival species, if any. 39 It is worth noting that, as a consecuence of the feedback between the population and 40 its environment, the fitness of each species varies with time. There is no proper way to 41 rank the species out of the particular context they are placed. 42 For the sake of clarifying, the main assumptions of the model are listed below: 43 (i) A carrying capacity exists that bounds the total population, 44 (ii) The environment is described by a unique scalar function, its temperature, that is 45 a function of time,

46
(iii) The influence of each species on the temperature is proportional to its population 47 size,

48
(iv) The fitness of each species is a symmetrical function of temperature with a single 49 peak at its optimal temperature,

50
(v) Species only interact with each other indirectly through this temperature and the 51 resource (space) constraint.

52
The next section presents in detail the mathematical model. Results concerning with 53 the stability analysis of population of low diversity are obtained in the third section.

54
Fourth section presents the results obtained from the simulations of populations with 55 large initial biodiversity. We conclude and discuss these results, as well as their 56 implications, in the final section. 57 2 The mathematical model 58 We consider a population of error-free self-replicative species (replicators) I i for 59 i = 1, . . . , S. The size of the total population at time t is N (t) and, can be determined 60 as the sum of the population of each phenotype N i (t): Species I i has associated a fitness function f i that can be described by two real 62 numbers: T i and α i . The first one, T i , stands for its optimal growth temperature 63 whereas, α i denotes its influence on the environment. Both parameters, T i and α i , give 64 a measure of the survival probability of I i in each generation. The relative fitness f i 65 depends on the rest of species of the population through the global temperature T . In 66 particular, we assume that the relative fitness of copy I i is given by where r i is a non-negative parameter and f (x) is a function that exhibits an unique 68 maximum at x = 0 and decreases monotonically to 0 as x → ∞ (see Figure 1).

69
The intensive parameter T characterizes the environment. We suppose that each 70 species I i has a linear influence on T , weighted by the real parameter α i . Assuming 71 that external perturbations that could modify the value of T are negligible, the time 72 evolution of the global temperature is given by: In order to induce a selective process, we assume that the system has a maximum 74 carrying capacity K, the upper bound of the total population. Let the function of time, 75 s(t), be the available space at time t: then, the global fitness of species I i is given by: At each time step, the reproduction rate of every species is a function of both the 78 population size and the external temperature and becomes null as the population size 79 approaches its maximum capacity K.

80
The time evolution of each species population N i can be described by a system of 81 Ordinary Differential Equation (ODE): for i = 1, 2, . . . , S. The parameter δ i is the death rate of species I i . For sake of 83 simplicity, if not explicitely indicated, we will assume in what follows that r i = r j ≡ r 84 and δ i = δ j ≡ δ for all i, j. Under these assumptions, species differ each other by their 85 optimal temperature T i and their capacity to modify the external temperature α i .

86
Notice that, as all species have the same r value, we can redefine the fitness function f 87 (see 2) including the common factor r in it, i.e. from now on: It is convenient to normalize the equations with respect to the carrying capacity, K. 89 Let x i (t) = Ni(t) K for i = 1, . . . , S. In so doing, we have to redefine the system 90 parameters adequately: r * = r K and α * i = α i K, although, we keep the same notation 91 in the equations: The non-linear character of this ODE system prevents its analytical solution.

93
Nevertheless, a complete qualitative analysis for different cases has been carried out and 94 it is presented in the following.

95
First of all let us find the equilibrium points. There are two posibilities to cancel the 96 differential equation for species x i in system 8. Either x i = 0, or the parenthesis on the 97 right vanishes, i.e., for every k such thatx k = 0 The local stability character of each of these equilibrium points is given by the 99 Jacobian matrix associated to system 8, evaluated on them.

100
For everyx l = 0, the l-row in the Jacobian matrix has only one non-null entry at the 101 diagonal position J l, l . Consequently, the associated corresponding eigenvalues are: In order to compute the rest of the eigenvalues, we only consider the ODEs that 103 correspond to the valuesx j = 0. If we arrange these variables in the m last equations, 104 we form a non-trivial box of dimension (m + 1) × (m + 1) in the lower right part of the 105 Jacobian. The equilibrium points of this subsytem verify: Note that this system of algebraic equations is inconsistent for large values of m. In 107 practice, coexistence of more than two species is impossible, except in degenerate cases 108 of species with the same optimal temperature, which we will consider indistinguishable. 109 It is straightforward to prove that if T r = T s for a pair of species we can redefine the 110 system using linear combinations of x r and x s in such a way that one of the 111 combinations is stationary and the other behaves as a new species with an α value that 112 is a combination of α r and α s and the system behaves exactly as the one we are 113 studying with one variable less.

114
When a solution exists, the associated Jacobian is straightforwardly computed. The 115 entry J k, j para k, j = 1 . . . m is given by: The entries of the last column are (∀k = 1, . . . , m): T =T = C k and the last row: With this notation, the structure of the Jacobian is: In the following sections we will study separately each of the equilibria and 120 determine their stability conditions over the parameters, particularly, on T i and α i . Obviously, neutral situations in which more than one species have the same fitness are 133 possible but these species are considered as indistinguishable in this paper.

134
The largest fitness corresponds to the species with optimal temperature T i closer to 135 zero, f (min|T i |) and, therefore this species (lets assume index k) is the only one that 136 can take over the whole population. Its equilibrium size is given by: From the qualitative analysis explained in previous section one can inferred that, as 138 α i = 0 for all i, the last row of the Jacobian has only one non-null entry, J S+1,S+1 = −1. 139 Besides, forx i = 0, (i = k) the corresponding rows in the jacobian matrix have again 140 only a non-null entry, J i,i . So the eigenvalues directly emerge: As f is a monotonically decreasing function, all these eigenvalues are negative.

142
The row of the Jacobian corresponding to each surviving species reads as follows: and the corresponding eigenvalue is then: If f (min|T i |) < δ all eigenvalues are negative ant this is the asymptotically stable 145 state, whereas if δ < f (min|T i |), the final system tends towards extintion. Let us consider the case when the species are able to modify the environment (α i = 0) 149 and analyze the different possible stationary states. As we stated before, contrary to the 150 quasispecies model, now coexistence equilibria of at most two species can occur. Let us assume that only one species has a population different from 0, for instance I k . 158 That is,x k > 0 andx i = 0 for all i = k. The equilibrium point is obtained from the 159 system: The eigenvalues associated to i = k are The condition for this equilibrium point to be asymptotically stable is: This last condition implies that: which means that the value ofT = α kxk must be closer to T k than any other T i .

165
At this point, as we will see later, it is convenient to reorganize the species according 166 to their optimal temperatures, satisfying 167 T 1 < T 2 < ... < T k < ...
in a increasing succession as depicted in Figure 1.

168
The former condition reads now: The rest of eigenvalues are those of the submatrix in the Jacobian: The associated characterisitic polynomial is: The equilibrium point is asymptotically stable if the two roots have negative real parts. 170 This occurs when:

172
Let us assume that the unique two survival species are I k and I m . The values of the 173 equilibrium population are the solution of the algebraic system: The first obvious condition that must be satisfied is δ < f k (T ).

175
The right part of the second equation is satisfied ifT = 1 2 (T k + T m ) or T k = T m .

176
The second case can be reduced to the case studied in the previous section, again by 177 using linear combinations of these two species. One of the new variables is stationary 178 and the other behaves as a new species with an α value that is a linear combination of 179 α k and α m .

180
Let us focus on the first one:T = 1 2 (T k + T m ). Eigenvalues associate to i = k, m are: 181 Again, the condition for this equilibrium point to be asymptotically stable is: Then, which means that the value ofT = 1 2 (T k + T m ) must be closer to T k and T m than any 183 other T i and, therefore, T m y T k must be consecutives, i.e. m = k + 1 (Hence the 184 convenience of arranging the species according to their optimal temperatures, as we 185 stated before). Now, the resulting linear system: can be straightforwardly solved yielding: For these two concentrations to be simultaneously positiveT must have a value in and α k 1 − δ f k (T ) .

190
As it will be proven next, asymptotical stability is discarded in the case α k+1 > α k . 191 For α k > α k+1 , as both values in equation 19 must be positive 19, equation 18 implies: 192 The remaining three eigenvalues associated to this coexistence equilibrium are 193 obtained from the reduced Jacobian matrix: Notice that A > 0 and B < 0 as f k is a decreasing function atT (see Figure 1).

197
The characteristic polynomial associated to J 2 is: The Routh-Hurwitz criterion requires the following conditions to assure three roots with 199 negative real parts: This implies that: α k > α k+1 as stated before. In addition,

202
A(x k +x k+1 ) + B(α k+1xk+1 − α kxk ) > 0 2.  2. In order to describe the basins of attraction of the two 2-coexistence equilibrium solutions in the four species model with the parameter setup given in Table 1, we show a set of 1500 initial conditions and mark with red circles those that converge to the 2-coexistence between species I 3 and I 4 and with blue circles those initial conditions that converge to the 2-coexistence of species I 1 and I 2 . The rest of the parameter values are r = 1 and δ = 0.1.

Four species model 207
In order to illustrate the results obtained in this section, we explore the numerical 208 solution of several four species cases where multi-stability appears. This model of four 209 initial species exhibits, among others, two kind of multi-stabilities that do not occur in 210 the classical models: (i) bi-stability between two different coexistence equilibrium points 211 and (ii) tri-stability among two single species equilibria and one coexistence point.

212
The first situation occurs for the species described in Table 1 The final equilibrium these four 217 species will achieve depends on the initial population. In Figure 2 we show the 218 projection on the plane Π (whose coordinate axis are x 1 + x 2 and x 3 + x 4 ) of the 219 solution obtained for 1500 initial conditions, all with T (0) = 0. As it can be seen, the 220 basin of attraction of each 2-coexistence is diffusively separated by the hyperplane 221 x 1 + x 2 = x 3 + x 4 that is projected as a line in the plane Π.

222
Another example of a case with four species that exhibits tri-stability is depicted in 223 Table 2. The rest of the parameters is fixed as before.   In both cases, depending on the initial conditions the population tends to one of 232 these equilibrium points. A tridimensional projection of 1000 initial conditions is 233 depicted in Figure 3 where the coordinates axis are: x 1 , x 2 and x 3 + x 4 . As it can be 234 observed, the coexistence equilibrum has the largest basin of attraction (blue points).

235
On the other extreme, the single species equilibrium withx 1 > 0 has the smallest basin 236 of attraction (only seven red circles corresponding to 1-existencex 1 ocurrences out of 237 1000 initial conditions). In the previous section, we have studied the stability properties of a system formed by 240 few species, concretely, four. We have proven the existence of conditions that yield 241 multistability between single species and coexistence equilibria. Under these conditions, 242 the asymptotic equilibrium population depends on the initial populations of each of the 243 species that form the population. In this section, we explore further this dependence in 244 populations formed by a larger number of species. An important result to be pointed 245 out is that coexistence equilibria with more than two species does not exist. We want to 246 show in this section how the asymptotic behavior of this kind of populations depends on 247 three factors: (i) the number of species, (ii) the initial conditions and (iii) the confluence 248 of the properties of the species that initially form the population, i.e. their optimal 249 temperature T k and the rate of influence, α k , over the environment temperature T . In order to study the influence of the number of species that initially form the 252 population, we carry out a numerical integration of the system of differential equations 253 that describes the time evolution of each of the species population, x k (t) (see 8). As 254 before, we assume that each species I k is characterized by T k and α k . The values of 255 these parameters are taken randomly from intervals whose length is changed in each where the proportion of extinction is significant (see Figure 4). For the other T -intervals 267 this proportion is null. As it can be seen, the occurrence of 2-coexistence equilibria is  The capability of some species to modify the environment has been postulated as a 282 factor of stabilization that can promote biodiversity [3]. The coupling between biota 283 and the environement could endowe the ecosystem with homeostatic properties that 284 enables a better adaptation [5,17]. The Gaia hypothesis is mainly based on this 285 interaction [18,19]. Indeed, Evolution favours those species that are able to keep the 286 environment under specific conditions that are assumed good enough for Life.  [17,[20][21][22]. Even more 299 challenging is to explain the existence of stable populations when mutation enables the 300 appearance of species that can drive the environment towards destabilization.

301
This paper shows that when species are able to modify the environment, coexistence 302 can be enhanced. As a matter of fact, for populations with a large number of species, 303 coexistence is the most probable equilibrium population. However, if this modification 304 is not coupled with the own species characteristics then, this capability does not assure 305 its own persistence. The question about which species will survive under specific 306 conditions can not be answered exclusively by carrying out a stability analysis.

307
Qualitative stability analysis, for instance, based on the linearization around equilibria 308 only provides local information about the dynamics of small perturbations: if the 309 perturbation decays the equilibrium points are said to be asymptotically stable. On the 310 contrary, if the perturbation increases over time the equilibrium point is said unstable. 311 When the system exhibits multiple equilibrium points this analysis cannot solve the 312 aymptotic behaviour when the population starts from a particular initial condition.

313
Linear stability analysis says nothing about the attraction basin. Nevertheless, we have 314 got a useful result, proven in section 4, that states a necessary condition that must 315 satisfy two species that coexist in equilibrium: their temperatures must be consecutive. 316 Unfortunately, this condition is not conclusive and further investigation is required to 317 predict the species selected from a given initial population.

318
In order to solve this problem, we have applied a different approach that seeks to 319 determine the equilibrium point (asymptotically stable) which is achieved from an 320 initial population. In particular, we assume that the species optimal temperatures and 321 their rate of influence on the environment are randomly taken. The population starts 322 with a given number of species homogeneously populated and whose sum is well below 323 the carrying capacity of the system K. After a transient period large enough to assure 324 the system relaxation, we take note of the survival species and analyse their properties. 325 We check that all these equilibria are asymptotically stable by a qualitative analysis. As 326 expected, depending on the values of the external temperature the equilibrium 327 population may differ. We note that most of the equilibrium populations are 328 coexistence of two species that present a definite relationship between their optimal 329 temperature T k and the rate of influence on the environment α k , concretely T k α k < 0. 330 Few simulations yield with the survival of only one species. None, as expected, with 331 more than three species. approach is correct we have also simulated the dynamic of the population by using an 337 individual-based algorithm (data not shown). We simulate a discrete population in 338 which each inidividual can be chosen either to replicate or to die according to its 339 (global) fitness. We note that for this value of K the results of these simulations agree 340 in all cases with the numerical integration of the ODE system. On the contrary, we  The stability of populations is a key issue in Ecology. Grimm and Wissel (1997) stressed 346 the large number of definitions that have appeared in Ecology [23]. These definitions 347 are classified into six classes that resume the concepts frequently used in this field: 348 constancy, resilience, persistence, resistance, elasticity, and domain of attraction [24]. In 349 some point or another all of them appear in this work. Specifically, we handle with 350 persistence when studying the species that remain in the population after a period of 351 time. Mathematical stability, as shown in section 2, provides information about the 352 local resilience and, in addition, enables to estimate the domain of attraction of 353 equilibria [7]. In contrast, this stability analysis does not inform about the persistence 354 of the species composition of an initial population that change over time due to the 355 existence of multiple stable equilibria [25]. In other words, this approach can not predict 356 the equilibrium point to which the system will tend. The consideration that species can 357 modify the environment and, as a consequence, the composition of the population, 358 includes an extra difficulty to perform a qualitative analysis and prevent an analytical 359 solution of the problem. Numerical integration and simulations are instead appropiate 360 tools for solving this problem as it is shown in this paper.

361
An important subject to be handled in a forthcoming paper concerns the 362 evolutionary properties of this kind of ecosystems. We have assumed that the mutation 363 rate of the species is null and consequently we have avoided the appearance of new 364 species in the population. In this kind of models that consider the ability of species to 365 modify the enviroment, the appearance of a new species with influence on the 366 environment can have important consequences on the equilibrium properties [16]. The 367 induced new conditions change the relative fitness of the species and, therefore, changes 368 the final outcome of the selective process. This modification can stabilized the 369 population by endowing the species with a larger survival probability or, on the 370 contrary, it can destabilize the population by pushing the external temperature out of 371 the limits of survival of the existing species. In the latter case, this induced-variation 372 can drive the population to extinction. In any case, absolute stability does not exist in 373 these models. The possibility of modifying the environment is an additional factor that 374 contributes to keep open forever the fate of evolution.