Evolution and extinction can occur rapidly: a modeling approach

The model for the self-development of a large ecosystem

In the model, the state of the ecosystem is characterized by the vector X(t) = (x₁(t), …, x_n (t)), where x_i(t) is the density of individuals, i—species with a self-reproduction potential at the moment t, n—number of species.

Further on, we will designate the totality of individuals with the potential for self-reproduction as biota. Let W_i denote the law of self-reproduction of biota of species i, and Y_i denote the law of the loss of self-reproduction potential, which we call the law of mortality. In general, the functional form of the laws of self-reproduction and mortality, W_i and Y_i, is determined by the structure of the external environment and connections that exist between individuals of the ecosystem and their environment. That is, W_i and Y_i are functions that depend on time t, the density of species X(t), and a number of parameters U = {u }. Then, self-development of the ecosystem can be described by the following system of differential equations: (1)${\displaystyle \frac{d{x}_i}{dt}\mathit{=}{k}_{in,i}+x_i\cdot \left(W_i-Y_i\right),i=1,\dots ,n.}$

Here, k_in,i denotes the rate of the species i biota inflow from the outside into the modeling area.

The model for the evolution of ecosystem parameters

In system Eq. (1), self-development is controlled by the laws of self-reproduction and mortality, W_i and Y_i. We assume that evolution of an ecosystem is the process of changing the values of a subset of parameters U = {u }, which increases the adaptability of species.

Further on, parameters that change their values over time will be designated as evolving parameters. Without loss of generality, we can denote them U = {u }. Thus, parameters from U are not constants, but functions of time. As a result, in order to calculate system Eq. (1), it is necessary to know the form of functions u(t)∈U. To find their form, we derive differential equations that satisfy u(t) ∈U. First, we note that if value of the parameter equals u(t) at the current time, then its rate of change over time satisfies the general equation ${\displaystyle \frac{du}{dt}=E_u\cdot u,u\in U,}$ where, E_u is the rate at which the parameter u changes.

Without limiting generality, we assume that mutational variability of a particular parameter is a result of mutations that occur during reproduction of individuals of a single specific species. Let σ(u) denote this species. Further, we take into account that E_u value is formed as a result of two successive processes: (1) change in the value of the evolving parameter due to mutations that arise during reproduction of individuals of a particular species; (2) evolutionary selection of mutations that increase the adaptability of a given species. The probability of a change in the u value at the time t is denoted by P_u, and by Ω_u: the rate at which the parameter u changes under the influence of evolutionary selection. Accordingly, E_u = P_u⋅Ω_u and the rate of evolutionary change u(t) satisfies the differential equation: ${\displaystyle \frac{du}{dt}=P_u\cdot {\Omega }_u\cdot u.}$

But, in general case, P_u and Ω_u are not constant values, but depend on the state of the ecosystem. Let us determined the type of P_u from the fundamental properties of living systems, that is the presence of inaccuracies, mutations, during genome replication in the process of self-reproduction. As a result, we obtain that P_u value is proportional to the specific intensity of reproduction of an individual: ${\displaystyle P_u=p_u\cdot W_{\sigma \left(u\right)},}$ where, p_u is the constant for mutation-associated change in the value of the evolving parameter u.

In general, p_u may not be a constant value, since replication accuracy may vary in the course of evolution, but for the purposes of our work, no further specification is required. Let us assume that p_u is a constant value.

To determine Ω_u, we assume that at each time moment t the adaptability of the σ(u) species to the environmental conditions can be estimated by the value of the “fitness functional” F_σ(u). Let them be known for all species. Now we consider that not all mutations are useful. Only those mutations that make a positive contribution to increasing the value of the “fitness functional” F_σ(u) are being fixed in the genome. Hence we express Ω_u through F_σ(u). To do this, we note that in the absence of a mutation in the parameter u in a short time h the current value F_σ(u) changes by the value $h\frac{d{F}_{\sigma \left(u\right)}}{dt}\left(t,U^{\prime }\left(t\right),u\left(t\right)\right)$.

Here, U′(t) denotes the vector of evolving parameters with excluded parameter u. Imagine now that a mutation occurred at the time t that introduced a variation in the u(t) value and it became equal to u(t) + Δu. Then at the moment t + h we have $h\frac{d{F}_{\sigma \left(u\right)}}{dt}\left(t,U^{\prime }\left(t\right),u\left(t\right)+\Delta u\right)$.

From here we obtain an increase in the “fitness functional” F_σ(u), which is achieved due to the mutational change in the evolving parameter u_j: ${\displaystyle \begin{array}{c}{\Omega }_u={\lim }_{\Delta u\to 0}li{m}_{h\to 0}\left[\frac{h\left(\frac{d{F}_{\sigma \left(u\right)}}{dt}\left(t,U^{\prime }\left(t\right),u+\Delta u\right)-\frac{d{F}_{\sigma \left(u\right)}}{dt}\left(t,U^{\prime }\left(t\right),u_{i,j}\right)\right)}{h\Delta u}\right]\hfill \\ ={\lim }_{\Delta u\to 0}\left[\frac{\frac{d{F}_{\sigma \left(u\right)}}{dt}\left(t,U^{\prime }\left(t\right),u+\Delta u\right)-\frac{d{F}_{\sigma \left(u\right)}}{dt}\left(t,U^{\prime }\left(t\right),u_{i,j}\right)}{\Delta u}\right]\hfill \\ =\frac{\partial }{\partial u}\left(\frac{d{F}_{{}_{\sigma \left(u\right)}}}{dt}\right).\hfill \end{array}}$

As a result, for each parameter u, we obtain an evolution equation of the general form (2)${\displaystyle \frac{du}{dt}\mathit{=}{p}_u\cdot u{W}_{\sigma \left(u\right)}\cdot \frac{\partial }{\partial u}\left(\frac{d}{dt}{F}_{\sigma \left(u\right)}\right).}$

Note that Eq. (2) does not include the “fitness functional” F_σ(u), but its total time derivative $\frac{d{F}_{\sigma \left(u\right)}}{dt}$.

That is the difference between Eq. (2) and analogous equation previously used in the models of the evolution of ecosystems (Sznajd-Weron & Weron, 2001). The set of equations Eq. (2) taken for all evolving parameters u ∈ U represents a model for the evolution of ecosystem parameters. Together with system (1), they represent a model for the evolutionary self-development of a large ecosystem. For clarity, we shall write it down (3)${\displaystyle \left\{\begin{array}{c}\frac{d{x}_i}{dt}=k_{in,i}+x_i\cdot \left(W_i-Y_i\right),i=1,\dots ,n,\hfill \\ \frac{du}{dt}\mathit{=}{p}_{u}u\cdot W_{\sigma \left(u\right)}\cdot \frac{\partial }{\partial u}\left(\frac{d}{dt}{F}_{\sigma \left(u\right)}\right),u\in U.\hfill \end{array}\right.}$

The law of biota self-reproduction

Let us assume that individuals of a “transient” population have a mixed (sexual and asexual) type of reproduction and sexual reproduction requires the meeting of two individuals. Then we describe the specific self-reproduction rate of a transit type biota by the classical logistic model, which assumes a simple positive feedback between diversity and its growth rate (parents—more descendants), and negative feedback based on the limited ecological space available for reproduction: (4)${\displaystyle W=\left(k_{aSx}+k_{Sx}x\right)\left(1-\frac{x}{C}\right)}$

In Eq. (4) k_aSx and k_Sx are rate constants for asexual and sexual reproduction, factor x corresponds to the sexual reproduction, factor (1 − x/C) describes the negative effect of biota density on the self-reproduction rate, C—maximum possible of biota density.

We used a logistic model because density dependence is a general tendency and fundamental principle of population ecology with potentially regulatory effects on population growth whereby the population has a propensity to increase when small and decrease when large (Hixon, Pacala & Sandin, 2002; Knape & De Valpine, 2012).

Also, according to AV Markov, hypothesis that the dynamics of the Phanerozoic marine biota calculated by traditional methods (without amendments) adequately reflects real changes in biodiversity has not been unproved and remains the most convenient and reliable basis for meaningful biological interpretations (Markov & Korotaev, 2007, p. 4).

The specific rate of mortality of the biota of a transit type we describe by the following function: (5)${\displaystyle Y=k_d+D\cdot H.}$

In Eq. (5), k_d denotes the specific rate of mortality, which is determined by environmental conditions. The second addendum (D⋅H) sets the specific rate of mortality, which is determined by the internal laws of the ecosystem’s functioning. We write it as multiplication of the rate constant for degradation D by the function H, which determines the mechanism of the influence of biota density on the mortality of individuals of the “transient” population. Assume that H is a bounded monotonically increasing function of the biota density. This assumption is implemented through the function: (6)${\displaystyle H=\frac{k_d\cdot x\cdot e^{\frac{x}{K_E}-1}}{K_D+x\cdot e^{\frac{x}{K_E}-1}}.}$

In Eq. (6) $r=e^{\frac{x}{K_E}-1}$denotes the coefficient of “auto-aggression”, K_D and K_E –parameters. The coefficient of specific “auto-aggression” in H(x) describes density-induced stress, which decreases the rate of reproduction (Ostfeld, Canham & Pugh, 1993; Wauters et al., 2004; Zhao et al., 2019) that may occur irrespective of taxon, life history and body size (Odden et al., 2014). In the model its value is monotonic, nonlinear, increases with increasing biota density x from 0 to 1: at zero density there is no “aggression”—H(0) = 0; at high density the coefficient of specific “aggression” tends to 1 − H(∞) = 1.

General form of the model of evolution of a “transient” population

Thus, the model (3) describes evolution as a process of ecosystem self-development (population of a “transit” type), during which there is a local increase in the adaptability of its individuals to the environmental conditions due to mutational variability and natural selection. The model has the following characteristics:

• mixed (asexual/sexual) type of reproduction,

• logistic law of reproduction W,

• exponentially rational law H of the influence of population density on the average life span of an individual,

• two evolving parameters C and D,

• function of adaptability F_fitness = W,

• influx or migration of live biota from the outside, k_in.

Why adaptability function F_fitness = W?

A priori, there are two natural measures for assessing the adaptation degree of a “transient” population to habitat conditions: W reflects the relative, specific efficiency of its self-reproduction, and xW—the absolute rate of its reproduction. We called them “fitness functionals”, F₁ and F₂, respectively. We have previously shown (Likhoshvai et al., 2017) that F₁ “fitness functional” allows to describe both the stable functioning of populations with asexual reproduction, for example, cyanobacteria that persist and develop for at least two billion years (Butterfield, 2007), and cyclic dynamics of the Phanerozoic marine biota (Rohde & Muller, 2005; Melott, 2008; Lieberman & Melott, 2007; Melott & Bambach, 2014; Melott & Bambach, 2017; Roberts & Mannion, 2019), which is characterized by sexual reproduction. Based on this, we assume that in model (3) F ≡W.

Further, if we write Eqs. (2) of model (3) for each evolving parameter, reproduction parameter C and mortality parameter D, and assume n = 1, we have the following evolutionary model of a “transient” population: (7)${\displaystyle \left\{\begin{array}{c}\frac{d}{dt}x=k_{in}+x\left(W-Y\right),\hfill \\ \frac{d}{dt}C=p_{C}WC\frac{\partial }{\partial C}\left(\frac{d}{dt}F\right),\hfill \\ \frac{d}{dt}D=p_{D}WD\frac{\partial }{\partial D}\left(\frac{d}{dt}F\right),\hfill \end{array}\right.}$ where x(t) is the density of transitive individuals and W and Y are the laws of self-production and mortality, respectively, which are determined based on the Eqs. (4) and (5). The parameters p_C and p_D stand for the rate constants for mutation-associated change in the values of the “reproduction” parameter C and “mortality” parameter D. There is no doubt that parameters p_C and p_D also changed values during evolution. But in the present study, we consider them to be constant. Since $\frac{d}{dt}F=W_x{x}_t+W_C{C}_t$, partial derivatives $\frac{d}{dt}F$ with respect to the parameters C and D have the following form: ${\displaystyle \begin{array}{c}\frac{\partial }{\partial C}\left(\frac{d}{dt}F\right)=\frac{\partial }{\partial C}\left[W_x{x}_t+W_C{C}_t\right],\hfill \\ \frac{\partial }{\partial D}\left(\frac{d}{dt}F\right)=\frac{\partial }{\partial D}\left[W_x{x}_t+W_C{C}_t\right],\hfill \end{array}}$ where ${\displaystyle \begin{array}{c}\frac{d}{dx}W=W_x=k_{Sx}-\frac{k_{aSx}+2k_{Sx}x}{C},\hfill \\ \frac{d}{dC}W=W_C=\left(k_{aSx}+k_{Sx}x\right)\frac{x}{C^2},\frac{\partial }{\partial C}\left(\frac{d}{dC}W\right)=W_{CC}=-2\left(k_{aSx}+k_{Sx}x\right)\frac{x}{C^3},\hfill \\ \frac{\partial }{\partial C}\left(\frac{d}{dx}W\right)=W_{xC}=\left(k_{aSx}+2k_{Sx}x\right)\frac{1}{C^2},\hfill \\ H=\frac{k_d\cdot x\cdot e^{\frac{x}{K_E}-1}}{K_D+x\cdot e^{\frac{x}{K_E}-1}},\frac{d}{dx}H=H_x=k_d\cdot K_D\frac{\left[\left(1+\frac{x}{K_E}\right)\cdot e^{\frac{x}{K_E}-1}\right]}{{\left(K_D+x\cdot e^{\frac{x}{K_E}-1}\right)}^2}.\hfill \end{array}}$

Let us find the solution of the functions C_t and D_t (the derivation is given in the Supplementary Information 1). After performing the calculations, we arrive at the following model of evolution: (8)${\displaystyle \begin{array}{c}\left\{\begin{array}{c}{x}_t=k_{in}+x\left(W-Y\right)\hfill \\ C_t=-p_{C}W\sum _{i=1}^{\infty }i\frac{u_i+Y{w}_i}{C^i}\hfill \\ D_t=p_{D}WDH\sum _{i=0}^{\infty }\frac{w_i}{C^i}\hfill \end{array}\right.\hfill \\ F=W,Y=DH+k_d,\hfill \\ W=\left(k_{aSx}+k_{Sx}x\right)\left(1-\frac{x}{C}\right),\frac{d}{dx}W=W_x=k_{Sx}-\frac{k_{aSx}+2k_{Sx}x}{C},\hfill \\ \frac{d}{dC}W=W_C=\left(k_{aSx}+k_{Sx}x\right)\frac{x}{C^2},\frac{d}{dC}\left(\frac{d}{dC}W\right)=W_{CC}=-2\left(k_{aSx}+k_{Sx}x\right)\frac{x}{C^3},\hfill \\ \frac{d}{dC}\left(\frac{d}{dx}W\right)=W_{xC}=\left(k_{aSx}+2k_{Sx}x\right)\frac{1}{C^2},\hfill \\ H=\frac{k_d\cdot x\cdot e^{\frac{x}{K_E}-1}}{K_D+x\cdot e^{\frac{x}{K_E}-1}},\frac{d}{dx}H=H_x=k_d\cdot K_D\frac{\left[\left(1+\frac{x}{K_E}\right)\cdot e^{\frac{x}{K_E}-1}\right]}{{\left(K_D+x\cdot e^{\frac{x}{K_E}-1}\right)}^2},\hfill \\ u_0=\left(k_{Sx}\right)\left(k_{in}+x\left(\left(k_{aSx}+k_{Sx}x\right)\right)\right),\hfill \\ u_1=-\left[k_{in}\left(k_{aSx}+2k_{Sx}x\right)+x\left(k_{aSx}+k_{Sx}x\right)\left(k_{aSx}+3k_{Sx}x\right)\right],\hfill \\ u_2=x^2\left(k_{aSx}+2k_{Sx}x\right)\left(k_{aSx}+k_{Sx}x\right),\hfill \\ u_3=-p_{C}x{\left(k_{aSx}+k_{Sx}x\right)}^2{u}_1,\hfill \\ u_i=p_{C}x{\left(k_{aSx}+k_{Sx}x\right)}^2\left[\left(i-3\right)x{u}_{i-3}-\left(i-2\right)u_{i-2}\right],i=4,\dots \hfill \\ w_0=-k_{Sx}x,w_1=\left(k_{aSx}+2k_{Sx}x\right)x,w_2=0,w_3=-p_{C}x{\left(k_{aSx}+k_{Sx}x\right)}^2{w}_1,\hfill \\ w_i=p_{C}x{\left(k_{aSx}+k_{Sx}x\right)}^2\left[\left(i-3\right)x{w}_{i-3}-\left(i-2\right)w_{i-2}\right],i=4,\dots \hfill \end{array}}$

Numerical calculations of the model Eq. (8), depending on the type of reproduction of the transit population (k_aSx, k_Sx) and the values of the evolutionary parameters (p_C, p_D) is shown below (see also Supplementary Information 2).

Cyclic and intermittent dynamics in the model of evolution of a “transient” population

The key information provided by numerical calculations was observation that only stationary dynamics x(t) is observed in model Eq. (8) at k_Sx=0 (no sexual type of reproduction) (Fig. 1A), that is, constancy (stasis) of diversity is established in the course of evolution of a “transient” population. At the same time, at k_aSx=0 (no asexual type of reproduction) in model Eq. (8) for certain values of the parameters p_C and p_D, which are the rate constants for mutation-associated change in the values of the reproduction parameter C and mortality parameter D, undamped oscillatory dynamics of x(t) variation was observed (Figs. 1B and 1C). Moreover, this dynamics became chaotic at certain values of the parameters p_C and p_D (Fig. 1D).

It is significant that oscillatory dynamics is generated due to only four factors: (1) presence of sexual reproduction; (2) negative impact of population density on the rates of reproduction and mortality of biota individuals; (3) genetic variation caused by mutations that occur during replication; (4) evolutionary selection directed towards increasing the adaptability of the individuals to habitat conditions.

Therefore, the results presented in Figs. 1B–1D already allow us to explain the oscillatory (cyclic) dynamics of changes in species diversity observed in marine biota in the last 500 million years without resorting to external influences.

However, the calculations obtained allow us to move forward. Namely, the same factors can explain such aspects of punctuated evolution observed in the fossil record as extinction catastrophes, rapid growth phases and stasis phases of species diversity. Indeed, analysis of the model revealed that in the oscillatory kinetics presented in Figs. 1C and 1D, qualitatively different phases of the x(t) value evolution are clearly distinguished. Let us consider them in more detail. Any arbitrary interval including a full cycle, similar to the one marked with a red bracket in Fig. 1C, is reasonable for this. Such cycle presents a typical scenario of the evolutionary turn of ecosystem development, in the sense that qualitative properties observed over a given time interval are reproduced in the course of evolution.

Four phases of the parameter x(t) evolution derived from the model Eq. (8) calculation at p_C = 0.0015 and p_D = 0.0205 are shown in Fig. 2. It can be clearly seen that f1 phase is characterized by a fast, on the scale of calculations, almost vertical decrease in x(t) to almost zero. This phase qualitatively corresponds to the degree of mass extinctions observed on Earth at least three times: 65, 250, 370 million years ago (see Fig. 1C in: Rohde & Muller, 2005). f2 phase is characterized by the low x(t) value. This phase qualitatively corresponds to the stages of life development that follow mass extinctions. This phase is characterized by a very low biodiversity level and very slow growth over a certain time interval. Then follows the phase of increasing species diversity, which can be divided into two stages. x(t) value rapidly increases during the first stage (Fig. 2, f3 phase). f3 phase of oscillatory dynamics qualitatively corresponds to the stages of explosive growth in the terrestrial biodiversity. The second stage is characterized by a high x(t) value, which continues to increase slowly (Fig. 2, f4 phase). f4 phase qualitatively corresponds to those stages of life development on Earth, in which a great variety of life forms and a relatively low growth rate are observed. Usually such stages of evolution are described as stasis.

In the model, losses and gains in species diversity inevitably repeat an unlimited number of times with approximately equal time intervals.

It should be noted that, depending on the calculation, f4 phase may consist of one or more segments of a gradual change in the x(t) value. In the case of several segments (see Fig. 1C), there are areas with a rapid decrease in x(t) value between them. However, the minimum achieved in these areas is much higher than the minima achieved in f1 phase; therefore, the observed decrease in x(t) can be interpreted as local catastrophe. It corresponds to the mild evolution that was shown in the calculation of model Eq. (8) with parameter values p_C = 0.0025, p_D = 0.045 (see Fig. 1B). Here, all evolutionary development consists of only f4 phase realization, in which oscillations of the parameter x(t) occur.

Thus, the modeling results have shown that if the efficiency of reproduction and mortality in a population depends on its density and the most adapted individuals, the genetic diversity of which is a result of genome replication errors during self-reproduction, are being selected, then these conditions are sufficient for the formation of cyclical intermittent dynamics of biodiversity in a living system with sexual type of reproduction.

Therefore, both cyclic and intermittent dynamics of changes in the diversity of organisms with sexual reproduction that make up the Earth’s global ecosystem observed in the last 500 million years may well be explained solely by internal laws of self-development with no impact of external factors.

The question arises—what properties of the system Eq. (8) determine the nature of the evolution of the x(t) curve?

Global extinction in the history of life as a bistability phenomenon

In Eq. (8), the “transient” biota density determined by the x(t) variable. The value of x(t) parameter in model Eq. (8) at each moment in time depends on the ratio between the rates of reproduction and mortality described by the functions W and Y. If we plot these functions at different times of evolution x(t), then we can note that curves W and Y can intersect up to three times (see Figs. 3B–3I). Each intersection corresponds to a stationary state, which can have the first equation of the dynamical system Eq. (8) at the fixed values of the parameters C and D.

Let us analyze the specific calculations of model Eq. (8) that are implemented at p_C = 0.0015, p_D = 0.02111 in the interval t∈[32000,33400] conv. time units in more detail. Transition from phase f2 to phase f3 and then to phase f4 occurs in the evolutionary cycle of biota development within the specified interval (see Fig. 3A). That is, it qualitatively corresponds to the stages of life development that follow global extinction.

Calculation of the model Eq. (8) at t = 32,000 is shown in Fig. 3B. At this time moment x(t)∼0.3096 and the first equation of system Eq. (8) has two stable stationary states: x_min ∼ 0.3 and x_max ∼ 26 with an unstable stationary state x_mdl between them. Since x_min∼x(t)<x_mdl at the given time moment, current x(t) position lies in the attraction region of the stationary state x_min. This means that if we fix values of C and D at the given time moment, x (t) will tend to the stationary state x_min. Such scenario is interpreted as a state of minimum species richness (corresponds to f2 phase). The fact that at the same time, system Eq. (8) has another stable stationary state x_max does not affect the state of the system, since x(t) value lies outside the region of its attraction and the system cannot get into it without external influences. Therefore, we can assume that at time t = 33000, the x_min stationary state is manifested, and the x_max stationary state is not manifested.

At every current moment, system Eq. (8) is in a position that is in the attraction area of the manifested stationary state, whose parameters determine the state of the system, and it is the increase in the integral adaptiveness of this system that drives the evolutionary progress.

In the model, this is achieved via mutational change in the C and D values and by fixing only those values that lead to an increase in adaptability. Indeed, the following calculation at t = 33100 (Fig. 3C) shows that, in the course of time, system Eq. (8) evolves toward an increase in the number and diversity of biota—x(t) value increases, while the attraction area of the stationary state x_min decreases; x_min and x_mdl stationary states converge, merge and disappear at some point within the interval t∈[33100,33140], so that we observe only one stationary state x_max at t = 33140 (Fig. 3D). Current x(t) value falls into the attraction area of this stationary state.

Since value of the stationary state x_max >  > x(t), an explosive increase in x(t) value is observed and x(t) value increases by an order of magnitude (Figs. 3D and 3E) and becomes larger than x_max (Figs. 3E and 3F) over a very short period of time (∼ per one conventional unit). As a result, f2 phase changes to f3 phase, during which the x(t) value reaches a new quantitative level. This event is accompanied by oscillations of x(t) with significant amplitude and small period (Fig. 3A), which arise as a result of the instability of the stationary state x_max (Fig. 3G). However, after a certain number of oscillations (Figs. 3H and 3I), the stability of the stationary state x_max is restored, the system goes into phase 4 and further evolution proceeds in the attraction area of the stationary state x_max.

Subsequently, evolution of the system Eq. (8) proceeds within the f4 phase, which can be characterized as a stasis phase (Fig. 4). During phase f4, x is quite high and, in the framework of the analysis, increases approximately twofold (Figs. 4B and 4D), and evolution proceeds in the attraction area of the stationary state x_max (Figs. 4B–4D). f4 phase is the longest and is not qualitatively homogeneous. It is divided into subphases called f4s—the monotonic growth of x(t) value, and f4n, which can be interpreted as local extinction. However, the exact number of subphases cannot be predicted due to the chaotic dynamics of evolution (Fig. 1D).

The exact duration of each phase cannot be predicted with a set of values p_C = 0.0015, p_D = 0.02111, since the dynamics observed in model Eq. (8) is chaotic (Fig. 1D).

Thus, we have demonstrated that such features of punctuated evolution observed in the fossil record as global extinctions and phases of rapid growth and biodiversity stasis are a reflection of the emergence of bistability (two distinct stable states) in a self-organizing system, which is the Earth’s biota. Transition from one stable state to another is accompanied by a rapid change in the parameters of the system, which may reflect the uneven pace of evolution observed in phylogenetic studies (Nichol, Rowe & Fitch, 1993; Pagel, Venditti & Meade, 2006; Wolf et al., 2006; Palmer et al., 2012). The reason for the transition from one state to another is the selection of the most adapted individuals, the genetic diversity of which is a result of genome replication errors during self-reproduction.