Conservation of population size is required for self-organized criticality in evolution models

In the DG model, a directed weighted network is considered, where nodes and links represent species and interspecies interactions, respectively. The fitness of each species (node) is defined as the sum of its incoming links:

where $a_{ji}$ is the weight of the link from node $j$ to node $i$. The species goes extinct if its fitness becomes less than the threshold $f_{\rm th}=0$ and survives otherwise. At each time step, the network is updated by migration of a new species and the subsequent extinctions. When a new species migrates, links between the new species and extant species are made with probability $c$ for each direction, and its weight is randomly drawn from a Gaussian distribution with mean $0$ and standard deviation $1$. (See Fig. 2 in Murase et al. (2010c) or Fig. 1 in Shimada (2014).) Because of the introduction of the new species, some of the existing species, including the new one, may have a negative fitness. The species with the lowest negative fitness in the system is removed. Because of the extinction of this species, the fitnesses of the neighboring species change accordingly. The fitnesses are calculated again, and the species having the lowest fitness is removed if its fitness is negative. The removal of species continues until all the species in the system have a positive or zero fitness, causing an avalanche of extinctions. Since extinctions may or may not happen, the number of species $N$ changes with time. In this setting, the system reaches a statistically stationary state after a sufficiently long initialization period starting from an empty network.

We show the statistical properties of the DG model in Fig. 1(a). The distribution of the fitness $P(f)$ is shown in Fig. 1(a-1). It is positive only for $f\geq 0$ by model definition. The extinction size for this model is defined as the number of the species that went extinct in one time step. The distribution of extinction size follows a simple exponential function as shown in Fig. 1(a-2). The interval between extinctions, which we define as the number of steps between consecutive extinctions (or the number of consecutive steps having no extinction), shows an exponential distribution as well (Fig. 1(a-3)). The time series of the number of species for this model shows a $1/f^{2}$ power spectral density, which indicates the dynamics is characterized by a simple Ornstein-Uhlenbeck process. All these statistics indicate that the extinction events occur randomly, characterized by a Poisson process. A non-trivial aspect of this model is found in the lifetime distribution of species. The lifetime of a species is defined as the duration between its immigration and extinction. The species lifetime distribution $P(L)$ follows a stretched exponential function, $P(L)\propto\exp{\left(-(L/L_{0})^{\alpha}\right)}$, with the exponent $\alpha$ close to $1/2$. This functional form is explained by the modified Red-Queen hypothesis, where an age-independent and $N$-dependent mortality is assumed Murase et al. (2010c, 2015).

Contrary to the DG model, the BS model assumes that an ecosystem has a fixed number of species $N$ and each species $i$ has a single fitness value $f_{i}$. For each step, the species having the minimum fitness is removed from the system (extinction) and immediately replaced with a new one (migration). The fitness of the new species is uniformly randomly drawn from $[0,1]$, and the fitness of its neighbors are also updated to a randomly assigned value, modeling the interactions with the new species. Following Ref. Bak and Sneppen (1993), we assume that the fitness barrier that must be overcome before a species can go extinct is proportional to its fitness. This leads to the time before extinction of the least fit species, $f_{\rm min}$, and its replacement with a new species being of the form $\tau_{\rm ext}(f_{\rm min})=t_{0}\exp(f_{\rm min}/f_{0})$, where $t_{0}$ and $f_{0}$ are constants. Here, $f_{0}$ must be sufficiently small to produce a broad region of critical scaling in the temporal dynamics. In the following, we use $t_{0}=1$ without loss of generality. So the actual time scale proceeds by $\tau_{\rm ext}(f_{\rm min})$ in each step. Although species were assumed to be aligned along a one-dimensional lattice in the first paper where the BS model is originally proposed Bak and Sneppen (1993), we mainly focus on the random-neighbor (or mean-field) version of the modelFlyvbjerg et al. (1993), in which a new species interacts with randomly chosen species. This simplification is instrumental not only for analytical treatment Flyvbjerg et al. (1993); Fraiman (2018) but for comparison with the DG model.

A statistically stationary state is obtained after a sufficient number of iterations. For comparison, we show the results for the BS model in Fig. 1(b). The distribution of $f$, $P(f)$, is shown in Fig. 1(b-1). It shows a profile similar to a step function: a positive plateau is found for $f_{\rm th}<f<1$, while it drops to an infinitesimal value for $0<f<f_{\rm th}$. When we see the temporal dynamics of the extinction events, the dynamics shows an intermittent behavior characterized by a power-law inter-event time distribution (Fig. 1(b-3)), reminiscent of punctuated equlibria. Here, the inter-event time is defined as the duration between two consecutive extinctions, i.e., $\tau_{\rm ext}$. The extinction size is defined as the number of consecutive extinctions whose $f_{\rm min}$ is smaller than $f_{\rm th}$. From the viewpoint of the time series, the extinction size is defined as the size of the “correlated bursts” Karsai et al. (2012, 2018). We obtain the clusters of events, or the bursty train, by splitting the event sequence by interval $\Delta=\tau_{\rm ext}(f_{\rm th})$, and the number of events in each group is regarded as the extinction size. In Fig. 1(b-2), the extinction size is shown. It is well fitted by a power law of exponent $-3/2$.

The species lifetime can be defined as the time between a species’ appearance and its extinction although it has not been studied in the literature to our knowledge. From our simulations, as shown in Figs. 1(b-4), the species lifetime distributions shows a power law followed by an exponential distribution. The exponent for the power law part is $-1$, which is same as that for the inter-event time distribution. The typical time scale for the latter part is approximately $\tau_{\rm ext}(f_{\rm th})$. Since the neighbor species and its renewed fitness is randomly selected, each species always has a positive finite probability of having a fitness below the threshold, which leads to its extinction. Thus, the species lifetime distribution has an exponential part although the overall profile is $L^{-1}$.

We summarize the key differences in the model definitions between the BS and the DG models as follows.

1. 1.

   While the number of species in the DG model fluctuates around its equilibrium value, the BS model has a fixed number of species which is given as a model parameter.

2. 2.

   The extinction threshold $f_{\rm th}$ is predefined in the DG model. On the other hand, in the BS model, the threshold of the fitness is self-organized as a result of immigrations and extinctions.

3. 3.

   The time required for immigrations and extinctions are different. In the DG model, immigrations occur regularly at each time step. The extinction immediately happens when a species has a negative fitness. By contrast, in the BS model, the time required for an extinction is dependent on the fitness value, which is followed by an instantaneous immigration of a new species.

4. 4.

   The fitness is defined as the sum of the incoming interspecies interactions in the DG model while the fitness for the BS model is defined as an attribute of a species.

The first three factors are not mutually exclusive but dependent on each other. These factors are reminiscent of the difference between canonical and grand canonical ensembles. The canonical ensemble, where the number of particles is fixed and the chemical potential is obtained as an outcome, may be regarded as similar to the BS model, where the number of species is fixed and the threshold is obtained as an outcome. Similarly, the grand canonical ensembles, where the number of particles is obtained given the chemical potential, is analogous to the DG model, where the number of species is obtained given the extinction threshold. Thus, it is not easy to incorporate only one of these factors.

The fourth factor, however, is independently testable by modifying the BS model so that the fitness of the species is defined as the sum of the interactions. In the following section, we will first show that the results for this “link-based” BS model remain qualitatively the same as those for the original BS model, implying that Factor $4$ is not the key ingredient of the SOC. We will then generalize the model to interpolate between the link-based BS model and the DG model to see the effects of Factors $1$, $2$ and $3$.

We first investigate a variant of the BS model, where the fitness of a species is defined as the sum of the inter-species interactions. More specifically, the ecosystem is represented by a weighted directed network whose nodes and links represent species and interactions between them, respectively. At each time step, the species with the minimum fitness as well as its links are removed. Similar to the BS model, the extinction is followed by an immediate immigration of a new species, whose incoming and outgoing links are made with probability $c$ and their weights are randomly independently drawn from a Gaussian distribution with mean $0$ and variance $1$. Thus, the fitness of the other species undergoes changes both by the eliminations of the links and by the introduction of the new links. The time required for the extinction to occur is the same as that for the BS model, i.e., $\tau_{\rm ext}(f)=\exp(f/f_{0})$.

This modification does not alter the results significantly from the BS model. See Figs. 1(c) for the results. As shown in Fig. 1(c-1), the fitness distribution shows a sudden drop at a certain threshold $f_{\rm th}\approx 0$ similarly to the BS model although the shape of the distribution and the threshold value are different. The extinction size is defined similarly to the BS model. The distribution of the extinction size $P(s)$ shows a power law with exponent $-3/2$ as shown in Fig. 1(c-2), which is the same as that for the BS model. In Fig. 1(c-3), the inter-event time distribution $P(\tau)$ is shown. It shows an approximate power law with exponent around $-1.2$, although the profile is slightly concave on the log-log scale. The lifetime distribution, shown in Fig. 1(c-4), is also approximated by a power law with exponent $-1$. Thus, we conclude that the model shows SOC behavior, which is essentially the same as the BS model. The link-based definition of species fitness does not alter the picture fundamentally.

We propose the following model in order to study the effect of the key differences 1, 2, and 3 discussed in Section II.3. We generalize the DG model, incorporating some ingredients of the BS model. The system is represented as a weighted directed network and a species’ fitness is defined as the sum of its incoming links. Following the assumption made in the BS model, we assume that the time required for an extinction of a species depends on the fitness as $\tau_{\rm ext}(f)=\exp(f/f_{0})$, where $f_{0}$ is a constant determining typical scale. In addition to this, we also consider the time between migration events in order to control the number of species. We assume that the time required for a migration depends on the current number of species and a parameter $\mu$ which controls the fluctuations in $N$: $\tau_{\rm mig}(N)=\exp\left(\mu(N-N_{0})/f_{0}\right)$. Here, $N_{0}$ is an input parameter of the model, denoting the expected number of species. When $\mu=0$, the migration time is always $\tau_{0}=1$, irrespective of the current number of species. There is no external force controlling the number of species, which corresponds to the situation for the DG model. When $\mu$ is large enough, on the other hand, the number of species is controlled by an external force. Migrations occur frequently when $N<N_{0}$, putting species into the system more often, while migration occurs less frequently when $N>N_{0}$. Thus, the number of species is under stronger control yielding smaller fluctuations in the number of species around $N_{0}$. In the limit of $\mu\to\infty$, immigration occurs immediately when $N<N_{0}$, while it never occurs when $N>N_{0}$. This results in a situation that immigrations and extinctions occur one after the other as in the BS model where $N$ is fixed.

The other rules are kept the same as the DG model. When a migration of species occurs, a species whose interactions are randomly determined is added to the system. For each possible link between existing species, we make a link with probability $c$ and the weight of the links are drawn from a Gaussian distribution with mean $0$ and variance $1$. The algorithm to update the system is summarized as follows.

1. 1.

   Find the species with the minimum fitness in the system, $f_{\rm min}$.

2. 2.

   Calculate $\tau_{\rm ext}(f_{\rm min})$ and $\tau_{\rm mig}(N)$.

   • •

     If $\tau_{\rm ext}<\tau_{\rm mig}$, an extinction of the least fit species occurs, i.e., the species is removed from the system.
     The current time $t$ is increased by $\tau_{\rm ext}$.

   • •

     If $\tau_{\rm ext}\geq\tau_{\rm mig}$, a migration of a new species occurs.
     The current time $t$ is increased by $\tau_{\rm mig}$.

The model has a correspondence to the DG model when $\mu=0$. This is because species with negative fitness are removed during one migration as $\tau_{ext}<\tau_{mig}$ for $f_{min}<0$. As long as there is a species with negative fitness, extinctions continue. When all the species have a positive or zero fitness value, the system accepts an immigrant. Although the actual time scale may be different from the DG model, the long-time behavior is approximately the same for a sufficiently small $f_{0}$ since $\tau_{\rm ext}$ is negligibly small. When $\mu\to\infty$, the model has a correspondence to the BS model. If $N>N_{0}$, extinction of the least fit species always occurs since $\tau_{mig}\to\infty$. If $N<N_{0}$, on the other hand, migration of a new species instantly occurs. Thus, migration alternates with extinction, keeping the number of species nearly constant around $N_{0}$. In other words, the least fit species are removed and immediately replaced by a new species. This dynamics is essentially similar to the BS model.

Typical time series for this model with $\mu=0$ and $0.1$ are shown in Fig. 2. As expected, the number of species is controlled when $\mu=0.1$ while it fluctuates widely for $\mu=0$. Hereafter, the parameter $N_{0}$ is set to $950$, which is comparable to the average number of species for $\mu=0$, in order to keep the number of species when changing $\mu$ and eliminate any side-effect caused by a change in $N$.

In Fig. 3(a), the probability distribution of $f$ for various values of $\mu$ are shown. As shown in the figure, the distribution has positive values for $f>0$ while it drops quickly for $f<0$, similar to the link-based model. This profile does not change significantly by changing $\mu$ and the threshold value seems to be around zero.

We define the extinction size as the number of consecutive extinctions whose interval is less than $\Delta=\exp\left(f_{\rm th}/f_{0}\right)$. This definition is consistent with those in the DG and the BS models when $\mu=0$ and $\mu\to\infty$, respectively. Using $f_{\rm th}=0$, the extinction size is calculated for various values of $\mu$ and is shown in Fig. 3(b). When $\mu=0$, the extinction size distribution decays exponentially. As we increase $\mu$, the tail becomes heavier and an overall power-law distribution is found for a sufficiently large $\mu$. When $\mu=0.1$, the distribution is approximately fitted by a power law with exponent $-1.5$. Although there is a slight deviation from a power law for large $s$ ($\sim 10^{4}$) and a non-monotonic dependence on $\mu$ is found, it may be due to a finite-size effect since the extinction size is the same order as $N_{0}$ when it deviates from the power law.

The distribution of inter-event time, which is defined as the period between two consecutive extinctions, is shown in Fig. 3(c). When $\mu=0$, the distribution consists of an initial power-law decay for $\tau<1$ and an exponential decay for $\tau\geq 1$. The former part corresponds to the extinctions within a single migration event, which is represented as the data point for $\tau=0$ in Fig. 1(b-3). For $\tau\geq 1$, the curve decays exponentially. It is consistent with the results for the DG model. The tail of the distributions gets much broader and closer to a power law as we increase $\mu$. In the figure, a line indicating a power law $\tau^{-1.2}$ is shown as a gray dashed line. The distribution is fitted fairly well by a power law although the curve is slightly concave in a log-log plot, indicating the time series are bursty. (One may also find the burstiness from Fig. 2.) ¹¹1When $\mu=10^{-1}$, there are small spikes for every two decades. These come from the discreteness of the $\tau_{\rm mig}$.

Similar transient behavior is also observed in the species lifetime distribution $P(L)$. As shown in Fig. 3(d), the profile changes significantly with $\mu$. When $\mu$ is small enough, the curve shows a quick decay while it has heavy tails for larger $\mu$. The curves for small $\mu$ are well fitted by stretched exponential functions with exponent close to $1/2$, which is similar to the DG model. As $\mu$ increases, the curve becomes closer to a power law. For $\mu=10^{-1}$, the distribution is approximated by a power law $L^{-1}$.

Therefore, all these statistics indicate a crossover from a non-SOC behavior to SOC with increasing $\mu$, showing that the constraint on the number of species is an essential factor for the emergence of SOC.