Mass extinction is a phenomenon in the history of life on Earth when a considerable number of species go extinct over a relatively short period of time. The magnitude of extinction varies between the events, the most well known are the "Big Five" when more than one half of all species got extinct. There were many extinctions with a smaller magnitude too. It is widely believed that the common trigger leading to a mass extinction is a climate change such a global warming or global cooling. There are, however, many open questions with regard to the effect and potential importance of specific factors and processes. In this paper, we develop a novel mathematical model that takes into account two factors largely overlooked in the mass extinctions literature, namely, (i) the active feedback of phytoplankton to the climate through changing the albedo of the ocean surface and (ii) the species's adaptive evolutionary response to a climate change. We show that whether species goes extinct or not depends on a subtle interplay between the scale of the climate change and the rate of the evolutionary response. We also show that species's response to a fast climate change can exhibit long transient dynamics (false extinction) when the species population density remain at a low value for a long time before recovering to its safe steady state value. Finally, we show that the distribution of extinction frequencies predicted by our model is generally consistent with the fossil record. Mass extinction is a phenomenon that several times punctuated the functioning of global biota during 541 Myrs (the Phanerozoic Eon) of the recorded history of the life on Earth [1,3,46]. Species get extinct all the time at some background rate (eventually being replaced by new, emerging species); this is regarded as a normal course of macroevolution [22,50,57]. However, sometimes the species extinction rate exceeded the average rate by more than an order of magnitude resulting in disappearing, over a relatively short period of time, of 50-75% of all the existing species [3,4]: a mass extinction occurred. Mass extinctions have been attracting a considerable and ever-growing scientific attention over the last several decades [2,13,15,36,40,49]. It is now widely believed that the main trigger leading to a mass extinction is a sufficiently large (overcritical) perturbation of the CO 2 cycle [43,44,45]. This can occur for different reasons, e.g. a massive volcano eruption called a LIP (Large Igneous Province) [4] or a collision of Earth with a large bolide [2]. An overcritical perturbation of the CO 2 cycle then leads to species extinction through a variety of mechanisms or pathways, e.g. global warming or cooling, ocean acidification and ocean anoxia [15]. While a considerable progress in understanding mass extinctions has been made over the last two decades (e.g. see the references above), yet many questions and issues remain [53]. One potentially important issue is the effect of the active feedback of a vegetation -in particular, phytoplankton -to a climate change. Studies on mass extinctions tend to regard species as passive entities that merely accept an environmental change (and get extinct if the climate change makes the environment too harsh). However, this is not always true as, in fact, some taxa can fight back to transform the environment according to their needs. On small temporal and spatial scales, such species are called "ecosystem engineers" [17,23], on a global scale this effect is summarized by the Gaia concept [28,30,31]. For instance, phytoplankton, when present in high densities, is known to be have a feedback on climate [7], however it remains unclear whether such effect can be extended to a global scale, e.g. to accelerate or slow down a mass extinction. Another open question is the role of species phenotypic plasticity: species can adapt to an environmental change, at least to a certain extent [38,39]. It seems obvious that plasticity can slow down the extinction rates through a directed change of species traits, hence affecting the extinction magnitude; however, this issue has been largely unaddressed in the literature. The inherent deficiency of the fossil record, see [53] for details, severely impedes further progress in understanding mass extinctions, in particular in distinguishing between the effects of different factors and mechanisms. This is a situation where, arguably, mathematical modelling can be used to partially compensate for the shortcomings of the data [53]. Surprisingly, the challenges of understanding mass extinctions have been largely overlooked by applied mathematicians. The literature concerned with mass extinctions modelling mostly focuses on a statistical analysis of the data, e.g. revealing the power laws in the frequency distribution [35,36,49] or attempting to reveal a pattern in the extinction timing [33,41]. There are very few papers attempting at a dynamical modelling of mass extinctions [8,10,42,54]. We mention here that there is a considerable literature dealing with modelling of extinctions of individual species (e.g. [24,27,56] but there is a huge epistemological gap between extinction of a particular species at a particular location and a mass extinction that can wipe out more than 50% of species globally. Upscaling the local/specific processes to a global scale is a highly nontrivial problem. In this paper, we develop a novel model that attempts to describe mass extinctions by considering an interplay between an active feedback of vegetation (particularly, phytoplankton) to a climate change and a differential species evolutionary response to the corresponding environmental change. The paper is organised as follows. In Section 2, we develop the baseline single-species model that includes the phytoplankton feedback on the global energy balance and the dependence of the plankton growth rate on the temperature (but not an evolutionary response). In Section 3, we investigate the properties of the baseline model to reveal possible species extinctions resulting from bifurcations of the steady states. In Section 4, the baseline model is extended to include species's adaptive evolutionary response to a climate change. We show that extinction or persistence of the species depends on the interplay between the magnitude of the climate change and the evolutionary rate.
We also show that the model exhibits long transient behaviour (false extinctions), so that species recovery from small densities can take a very long time. In Section 5, we turn the single species model to a multi-species one. We investigate how the extinction frequencies depend on the values of the key parameters or their distributions and in Section 5.1 we endeavour to tentatively compare the simulation results with the fossil record on mass extinctions. Finally, in Section 6 we discuss our findings and outline possible directions of future work. Following a widely used conceptual approach [9,11,29,48], in this study we assume that the state of the Earth climate can be described by a single variable, i.e. the average Earth temperature. Correspondingly, for the dynamics of the global climate we use the zerodimensional Budyko-Sellers model described by the following equation [11]: where T is the average surface temperature, t is time and  is a coefficient known as thermal inertia. In the right-hand side of Eq. ( 1), the first and second terms describe, respectively, the outgoing light emission (described by the Stefan-Boltzmann law with the effective emissivity e) and the fraction of the incoming solar radiation that falls on the Earth surface where S is the surface albedo, I 0 quantifies the total amount of solar energy and  0 is a coefficient. Obviously, any change in the Earth surface albedo will change the balance between the incoming and outgoing energy fluxes, thus resulting in an increase or decrease of the average temperature and hence potentially leading to a global climate change. One factor that may change the albedo is known to be vegetation cover. While the effect of the terrestrial vegetation is still under-investigated, the effect of the aquatic vegetation -i.e. phytoplankton -is well established [7]: the presence of phytoplankton tends to decrease the ocean surface albedo; see Fig. 1a. Since about 70% of Earth is covered by oceans, any global tendency in phytoplankton density to increase or decrease is likely to have a significant effect on the Earth temperature. For our baseline model, we assume that, as well as the temperature, the ocean phytoplankton can also be described by a single variable, i.e. the average phytoplankton density in the upper ocean layer -say, u(t) -hence considering phytoplankton as a single 'meta-species'. Apparently, this disregards not only the geographical variations in the plankton density but also the presence of different plankton species; the effect of the latter will be discussed below. Arguably, in such schematic, conceptual model the dynamics of phytoplankton as a whole can be modelled similarly to the dynamics of individual populations, e.g. by considering the logistic growth [14,18,51]: where  is the mortality rate and K is the carrying capacity, i.e. the maximum amount of phytoplankton that can be sustainably contained by a unit volume of the sea water. Note that the phytoplankton reproduction rate B(T ) depends on the water temperature, for which there is considerable empirical evidence. It is well known that the growth rate increases for low and intermediate temperatures,in particular in the range 0 < T < 30  [26]. The temperature dependence for higher water temperatures (e.g. for T > 40  ) is less known but it seems intuitively obvious that it should be decreasing (as, ultimately, too high temperatures will kill phytoplankton altogether). Thus, if considered over a sufficiently broad temperature range, function B(T ) should be bell-shaped . In order to complete the model, we now need a specific parametrization of functions S(u) and B(T ). While it is clear that the dependence of the surface albedo on the phytoplankton density is a monotonously decreasing function (cf. Fig. 1a), its precise shape is unknown. Since our goal here is to build a model which is qualitative rather than quantitative (i.e. which predicts the tendencies but not necessarily specific numbers), instead of looking for a precise shape, we use the following hypothetical function: where s 0 , s 1 and  1 are coefficients. Note that, while the exponential decay in (3) is hypothetical, the marginal values S(0) and S() (and hence coefficients of s 0 and s 1 ) can be estimated based on available data for the value of the ocean surface albedo for a clean water and for the turbid water with high plankton density, For the dependence of the phytoplankton growth rate on the water temperature, we use the following parametrization: where b 1 , T 0 and  2 are coefficients. It is readily seen that B(T ) defined by function (4) has a unique maximum at T m = (T 0 / 2 ) 1/2 and decreases to zero for |T  T m | 1. Note that the first exponent in ( 4) is essentially the Arrhenius law describing the generic dependence of chemical reactions (i.e., in our case, metabolic reactions) on the temperature, cf. [26]. Thus, we arrive at the following temperature-phytoplankton system: where we have introduced coefficients a = e and b =  0 I 0 4 for the notations simplicity. We now investigate the properties of the model (5)(6), particularly the number and stability of its steady states and how that may depend on the model's parameters. The steady states of the dynamical system (5-6) are the solutions of the following algebraic system: Equations (7)(8) define the two (null) isoclines of the system, which we will call the temperature isocline and the phytoplankton isocline, respectively. The temperature isocline is therefore given by The phytoplankton isocline consists of two parts, i.e., of the following curve and the straight line u = 0. It does not seem possible to solve the system (7-8) analytically. However, given the fact the system's equilibria are the intersection points of the two isoclines, it is possible to reveal the tendencies leading to the system's bifurcations (e.g. the disappearance or emergence of the steady states) by revealing the change of the relative position of the isoclines (i.e. the relative position of the corresponding curves) in response to a change in parameter values.
For the reasons that will be explained in Section 4, we consider b 1 as the main controlling parameter but we also reveal the effect of some other parameters such as s 0 , s 1 and . We now consider s 1 as the controlling parameter. The change in the relative position of the isoclines and the corresponding bifurcation structure of the system are shown in In the case where either s 0 or  are considered as controlling parameters, the system shows a sequence of qualitative changes in its properties similar to that shown in Fig. 3; we therefore do not show those results here for the sake of brevity. Note that the bifurcation values of parameter b 1 are not universal but depend on the values of other parameters, in particular on . As an example, Figure (4) shows the critical values of b 1 corresponding to the saddle-node bifurcation in the system (5-6) as a function of  obtained for two different values of s 1 and having other parameters fixed at the same hypothetical values as above ( = 1, a = 1, b = 0.3, s 0 = 0.1,  1 = 3,  2 = 1, and T 0 = 1). In order to obtain the temperature and the the phytoplankton density as functions of time, system (5-6) is solved numerically. Figure 5 shows the results obtained for parameter values b = 0.3, a = 1, T 0 = 1,  2 = 1,  = 0.1,  1 = 3, s 0 = 0.1, s 1 = 0.95,  = 1 and a few different values of b 1 . In Fig. 5a, we choose b 1 = 1.98 which is a subcritical value, i.e. when the system possesses a stable positive equilibrium . In this case, the initial conditions fast approach their asymptotical positive steady state values. Figures 5b and 5c are obtained for supercritical values b 1 = 1.96 and b 1 = 1.93, respectively, i.e. after the saddle-node bifurcation occurs, so that there are no positive equilibria. Correspondingly, in the course of time the system eventually approaches the semi-trivial steady state ( T0 , 0), hence resulting in the phytoplankton extinction. Interestingly, for the values of b 1 below the bifurcation value but close to it, the extinction does not actually happen until after a rather long period of time when the state variables change very slowly, cf. Fig. 5b. This type of dynamics is known as a 'long transient' [16,34] and, in our case, is a result of the ghost attractor [21,52,55]. The stable equilibrium has disappeared in the saddle-node bifurcation but, in the part of the phase pane where it was before the bifurcation, the phase flow remain slow. Thus, the system spends a long time there, which, ultimately (for a parameter value sufficiently close to its bifurcation value), may give an impression of a steady state -a ghost attractor. For a smaller value of b 1 , i.e. further away from its bifurcation value, the long transient becomes shorter (cf. Fig. 5c) and eventually disappears. One important phenomenon that is entirely missed by the above analysis is possible species adaptation to the environmental changes. Indeed, in the model (5-6), all parameters are fixed and that means that the species traits remain unchanged in the course of the dynamics. Meanwhile, there is considerable evidence showing that, in many cases, populations can evolve sufficiently to survive in new or altered environments, e.g. resulting from the climate change or human-related activities such as biological invasions, pollution, etc. [5,19,32,37]). Models of adaptive evolution are numerohoffmann1991,us and any overview of them is beyond the scope of this study. Here we use the approach based on the conceptual model by Gomulkiewicz & Holt [12] who showed that a modification of a given species's trait in response to a sudden environmental change can be described by an exponential transition of the corresponding parameter from its 'old' (before change) value to a new one. In our model, a change in the temperature affects the phytoplankton growth rate. We therefore assume that the phytoplanckton evolves with time in order to reach an equilibrium with the changed environment. Specifically, to account for this adaptive evolution we assume that parameter b 1 is the following function of time: where and  is the parameter quantifying the rate of adaptive evolution. Thus, b old is the parameter value that ensures that, before the climate change, the phytoplankton is at the positive steady state. The effect of a change in b 1 can be understood in As a result, the plankton density is not at at equilibrium any more and hence undergoes the dynamics as prescribed by the Eqs. (5)(6). A big question is then whether this dynamics brings the phytoplankton to a new positive steady state (i.e. another point at the upper steady branch of the bifurcation curve) or to extinction. Apparently, in terms of our approach the answer to the above question depends on two factors. The first factor is how far the sudden environmental change pushed the bifurcation curve away, i.e. how far the before-change value of the phytoplankton density is away from the new position of the bifurcation curve. The second factor is the evolutionary rate: intuitively, a high enough evolutionary rate should help the species to avoid extinction. To illustrate the effect of the first factor, Fig. 6 shows the solution of the system (5-6, 11) obtained numerically for parameter values b new = 2.1,  = 0.01 and different values of b old , other parameters are the same as above. (Note that, for the convenience of simulations, we quantify the difference between the pre-change and the post-change equilibrium values of parameter b 1 by varying b old rather than b new .) We readily observe that there is a critical value of the difference  = b new  b old .
For  less than the critical value, the system overcomes the effect of the climate change so that, after an initial decay in the values of T and u, they approach their new positive steady state values; see Figs. 6a and 6b. However, for  larger than the critical value, after a certain period of time when T and u show almost no changes, they eventually approach the semi-trivial steady state where phytoplankton goes to extinction and the temperature stabilizes at a lower value; see Figs. 6c and 6d. Interestingly, as well as in our baseline model without species adaptation, for the values of  close to its critical value, the final transition of the system to its asymptotic state, whether it is to the positive state or to the semi-trivial (extinction) state, does not happen until after a rather long period of a quasi-steady state dynamics when T and u do not show any visible changes (cf. Figs. 6b and 6c). This is another example of the long transient dynamics [16,34] shown by our global climate-vegetation model, cf. In conclusion of this section, we look into the properties of the long transient dynamics shown in Fig. 6b. Specifically, we are interested to know how the duration of the transient -say, y -to which we refer as the recovery time (the time taken by the system dynamics to recover the phytoplankton density from the very small values to which it initially falls) depends on the closeness of parameter b old to its bifurcation value: the property known in the literature as the scaling law, see [34] and further references there. Figure 9a shows the recovery time obtained for different values of b old . The results indicate that when the value of b old approaches its critical value, the recovery time grows to infinity; however, the rate of the increase (the scaling law) is unclear. In order to reveal the scaling law, in Fig. 9b the results are shown on semi-logarithmic scale and are well approximated by a straight line. We therefore arrive at the following scaling law: where a and b are certain coefficients. Thus, the long transient's duration grows rather slowly compared to the exponential and power laws usually observed in nonlinear systems [34]. In this section, we attempt to abate one of the main limitations of our model, namely, that the phytoplankton density is described by a single variable while in reality phytoplankton consists of many different species. Correspondingly, we now consider an extension of the model where different phytoplankton species are included explicitly. Let phytoplankton community consist of m species, then the state of system is described by m variables, u = (u 1 , u 2 , . . . , u m ), where u k is the kth species density. As above, we assume that the growth of each species is described by the logistic equation but with different (species-specific) parameters. The species compete for common resources; we assume that the competition affects species per capita growth rates but not carrying capacities. Thus, we arrive at the following model: (cf. [20,25]) where B k is the kth species's per capita growth rate, and R = (R 1 , . . . , R n ) are resources, e.g. nutrients [20]. In this paper, for the sake of simplicity we consider the case where resources are aplenty and hence their availability is not a limiting factor for the species population dynamics. Therefore, the dependence of B k on R in Eqs. ( 13) can be dropped. Instead, we focus on the effect of the temperature, for which we consider the same parametrization as above, see Fig. 1b and Eq. ( 4), but with some of the coefficients being now species-specific. Thus, we arrive at the following multi-species temperature-phytoplankton system: where L(u) = m k=1  1,k u k . The species can evolve to adapt to the climate changes so that where we assume that parameters b old,k and  k are species-specific but b new is the same for all species. The choice of b old,k and  k requires a brief discussion. These parameters quantify certain biological traits of phytoplankton species. The traits of currently existing species have developed in the course of evolution. Evolution acts through a random selection aiming to maximize the species fitness. Because of this inherent randomness, a particular biological parameter can often be regarded as a random number drawn from a certain probability distribution [47]. This approach is particularly relevant in the case of a community consisting of many similar species [6]. Correspondingly, for the values of b old,k and  k , we consider them being randomly distributed, that is, the arrays {b old,k , k = 1, . . . , m} and { k , k = 1, . . . , m} both consist of m random numbers. One can expect that the properties of the multi-species system (14-15), in particular the fraction of species that goes extinct in response to a sudden climate change, depend on the properties of the probability distributions. As we will show it below, this is indeed the case. We begin with the example where b old,k are the same for all species but  k are described by a normal distribution. Figure 10 shows the results obtained for a community of m = 10 species. The coloured curves visualize the ten corresponding trajectories in the (b 1 , u) plane obtained for b old = 0.903 and ten different random values of  k . Simulations are repeated three times, as shown in panels (a), (b) and (c). Since in each simulation the array of ten random values for  k is different, the results are different too. In particular, the fraction of species that goes extinct changes between the simulation runs. In Fig. 10a, five out of ten species go extinct while in Figs. 10b and 10c it is, respectively, nine and two species go extinct.
Thus, by repeating the above simulations sufficiently many times, one can obtain a distribution of extinction frequencies. That is, by repeating the simulations N times, in N 1 simulation runs only one species goes extinct, in N 2 runs two species go extinct, etc., with In order to understand how the agreement between the simulation results and the data could be improved, we first consider a different distribution of  k . Figure 11b shows the distribution of extinction frequencies obtained for  k being exponentially distributed. We readily observe that the shape of the distribution does not change much, remaining bellshaped, with only a slightly different position of the maximum. A key to resolving this difficulty (i.e. the apparent disagreement with fossil record) lies in the understanding that, generally speaking, some other parameters of the model should also be regarded as random values. Specifically, we now proceed to a more general case where b old is also random. Figure 12 shows the distribution of extinction frequencies obtained in system ( 14 where b and w are parameters. Figures 12a, 12b and 12c shows the results obtained for w = 1.8 and b = 2, 2.1, and 2.2, respectively. In all three cases, m = 100 values of  are distributed exponentially with mean=0.03. Thus, changing constant parameter b old to an array of random numbers changes the results dramatically. Instead of a bell-shaped distribution of extinction frequencies (as shown in Fig. 11), we know obtain a monotonously decreasing distribution, which is in a qualitative agreement with the fossil record [58], see also below. In this section, we endeavour to compare our simulation results to historical data on extinctions in a more explicit way. Fossil record contains data on 163 extinction events (including the "Big Five") with the extinction magnitude, i.e. the estimated percentage of all species that went extinct, ranging approximately from 5% to 68%. The thick black curve in Fig. 13 shows these extinctions events in the ranked order, i.e. from the largest to the smallest. We now recall that different values of b old can be regarded as a climate change of different magnitude. Correspondingly, in order to provide a direct comparison between the simulation results and the data, we now repeat the simulations with 163 different values of b old and show the obtained fractions of extinct species in the ranked order, i.e. from largest to smallest. The results are shown in We readily observe that, although our results are in a general agreement with the data, e.g. predicting mass extinctions of the magnitude comparable with those shown by the fossil record, there are apparent disagreements too. In particular, our results overestimate the frequency of large-scale extinctions (with more than 60% of species going extinct) and slightly under-estimate the frequency of small extinctions. Possible reasons for this disagreement are discussed in the next section. Mass extinctions were an important part of the history of life on Earth [3,46,58]. It is widely believed that the ultimate reason for a mass extinction is a sufficiently large climate change [4,53]. This, in turn, can be caused by a sufficiently large (over-critical) perturbation of CO 2 cycle [43]. A large climate change can eventually lead to species extinctions and considerable biodiversity loss on the global scale through a variety of specific mechanisms or pathways, the most common one being a significant change in the average Earth temperature resulting in global warming or global cooling. While there has been a great progress over the last two decades in understanding the main causes and triggers leading to a mass extinctions, e.g. see [1,3,4,15,33,43,44,58], many questions remain open. In particular, it is well known that not every climate change in the Earth history resulted in mass extinction; thus, there must exist factors or feedbacks that can attenuate the effect of the change. Due to the intrinsic deficiencies of paleontological data (for details and a discussion of this issue see Section 3 in [53]), it does not seem possible to identify such factors or feedbacks based solely on the analysis of fossil record. Mathematical models are needed; mechanistic process-based models create a 'virtual laboratory' where specific hypotheses can be tested and various scenarios investigated. Mathematical modelling has long become a part of a standard toolbox in ecology as well as environmental sciences more generally; surprisingly, application of models in paleontology is still rare. In this paper, we have developed a new mathematical model with a potential to describe a mass extinction using the phytoplankton community as a proxy to broader taxa. In particular, our model reveals how the extinction magnitude can be affected by the climatevegetation coupling (through phytoplankton active feedback to the climate) and by the species evolutionary response to the climate change: the two factors that are apparently important yet have been largely overlooked in the existing literature. We have shown that our model (see Eqs. (14)(15)(16)) predicts the distribution of extinction frequencies (e.g. as in Figs. 12 and 13) that is consistent with the data. In spite of the general consistency, there are some apparent differences between the ranked order of extinctions predicted by our model and that shown by the fossil record. The simulated distribution is described by a steeper curve, hence overestimating the frequency of large extinctions and underestimating the frequency of small ones. However, we believe that this is not an inherent flaw of the model and the agreement with the data can, in principle, be improved. In particular, here we recall that the results depend significantly on the type of parameter b old . Considering it as an exponentially distributed random value (instead of a fixed single value) leads to a qualitative change in the extinctions frequency distribution changing it from a bell-shaped distribution to a monotonous one, cf. Figs. 11 and 12. Thus, one can expect that a different choice of the probability distribution for b old (and perhaps for ) may cause further changes in the extinctions frequency distribution making it closer to the data.
This should become a focus of future work. At the same time, we also want to mention that the fossil record gives only, at best, a partial view of the actual extinction magnitude. The matter is that the fossil record contains predominantly data on hard-bodied species (e.g. mollusks and vertebrates) while soft bodied species (to which most of phytoplankton species belong) usually disappear without leaving any trace. How the distribution of extinctions (cf. thick black curve in Fig. 13) may change if data on soft-bodied species are included is an entirely open issue. This may partially explain the disagreement between our model and the fossil record. A model allowing for a direct comparison with the fossil record has to include explicitly a hard-bodied taxa, e.g. some zooplankton species. That should become another focus of future work. Since most of the results in the main text were obtained numerically, the question arises as to how a change in parameter values may change the system's properties. We are mostly interested to the effect of variations in , as this parameter describes the rate of species's adaptive evolution. Correspondingly, below we repeat some of the simulations of Section 4, specifically those shown in Figs. 6, Since most of the results in the main text were obtained numerically, the question arises as to how a change in parameter values may change the system's properties. We are mostly interested to the effect of variations in , as this parameter describes the rate of species's adaptive evolution. Correspondingly, below we repeat some of the simulations of Section 4, specifically those shown in Figs. 6,