Criticality, adaptability and early-warning signals in time series in a discrete quasispecies model

Abstract

Complex systems from different fields of knowledge often do not allow a mathematical description or modeling, because of their intricate structure composed of numerous interacting components. As an alternative approach, it is possible to study the way in which observables associated with the system fluctuate in time. These time series may provide valuable information about the underlying dynamics. It has been suggested that complex dynamic systems, ranging from ecosystems to financial markets and the climate, produce generic early-warning signals at the “tipping points,” where they announce a sudden shift toward a different dynamical regime, such as a population extinction, a systemic market crash, or abrupt shifts in the weather. On the other hand, the framework of Self-Organized Criticality (SOC), suggests that some complex systems, such as life itself, may spontaneously converge toward a critical point. As a particular example, the quasispecies model suggests that RNA viruses self-organize their mutation rate near the error-catastrophe threshold, where robustness and evolvability are balanced in such a way that survival is optimized. In this paper, we study the time series associated to a classical discrete quasispecies model for different mutation rates, and identify early-warning signals for critical mutation rates near the error-catastrophe threshold, such as irregularities in the kurtosis and a significant increase in the autocorrelation range, reminiscent of 1/f noise. In the present context, we find that the early-warning signals, rather than broadcasting the collapse of the system, are the fingerprint of survival optimization.

Keywords

time series / complexity / early-warning signals / quasispecies / 1/f noise / optimization

Introduction

Complexity theory offers a new paradigm to study intricate compound systems in many different fields of knowledge, such as physics, biology, physiology and ecology. These complex systems are composed of many mutually interacting components that induce new emerging properties in the whole system that are not present in the individual components. One particular interesting example from ecology is socially organized animals, such as ant colonies (Solé et al., 1993a, 1993b; Miramontes, 1995) or starling flocks (Ballerini et al., 2008; Cavagna et al., 2010), where each individual animal interacts directly only with its nearest neighbors. These local interactions appear to be ajusted in such a way that somehow long-range correlations emerge, which permit the group to act as an organized social “super-organism.” Interestingly, the emergence of such long-range correlations and the coexistence of different scales appear to be universal features for complex systems in many fields of knowledge, suggesting that a unifying mathematical framework might exist that describes this phenomenology. One example of such a unifying concept is auto-organized criticality, which states that some complex systems converge toward a critical state where long-range correlations emerge spontaneously (Kauffman, 1995; Bak, 1996).

On the other hand, the apparent similarities between different complex systems could be an illusion and future understanding might demand a detailed analysis and modeling for each specific case. Often, however, these complex systems are too complicated to allow for a precise mathematical description and it is desirable to develop alternative approaches that overcome this limitation. One way out is to study the signals or the time series that the system produces, i.e. the way in which a particular observable associated with the system fluctuates in time. The main assumption is that the statistical properties of these series provide information about the underlying dynamics of the system and- more particularly- about its possible critical states (Fossion et al., 2010; Landa et al., 2011). Recently, it has been proposed that complex systems from many different fields of knowledge produce in their time series early-warning signals that are generic, and that are able to predict an impending collapse or phase transition of the system (Scheffer et al., 2001, 2009, 2012; Carpenter et al., 2011). Where at first sight it might appear strange that ecological systems such as populations and physical systems such as the global climate produce similar signals, it results that near a critical point the dynamics of a complex system has universal properties, regardless of the differences in the details of each system. One example of such a universal property is that after a perturbation complex systems in a critical state take longer to return to their equilibrium state. This phenomenon is called critical slowing down, which can be recognized as an increase in autocorrelation and variance of the associated time series. This leads to a “reddening” (a favoring of the lower frequencies) of the corresponding Fourier power spectrum. Other characteristics are flickering and an increased skewness of the distribution function.

Evolving populations are examples of a complex system, where robustness and evolvability–that is, a balance between resisting and allowing change in the internal states to adapt to different internal and external conditions- are characterized by multiple spatial and time scales (Lenski et al., 2006). In the context of populations of cells, these properties are intrinsically related with an optimal critical state at which cells regulate and control their phenotype (Aldana et al., 2007). This dynamics can be studied theoretically within the framework of the quasispecies model, where one finds an error-catastrophe threshold that is associated with an optimization of the genetic variability of the species (Bull et al., 2005), and which is observed experimentally in RNA viruses (Crotty et al., 2001). Briefly, the error-catastrophe model states that the transfer of genetic information between subsequent generations in an individual population depends on the mutation rate μ. More specifically, there is a critical mutation rate μ_crit below which genetic information is passed on between generations but for which the population does not evolve in time. Conversely, in populations with mutation rates higher than μ_crit, evolution is promoted, but there is a loss of hereditability between the different generations. Given the fact that viruses proliferate when they evolve more rapidly than the immune system of their host can recognize them, RNA viruses might thus self-organize and converge toward this critical mutation threshold μ_crit where both functions of robustness and evolvability are optimized. The quasispecies model and the notion of criticality inspire antiviral treatments in two ways. On the one hand, a treatment that enhances the mutation rate pushes the virus over the error threshold after which it automutates into a harmless entity, an effect wich is also called “lethal mutagenesis” or “mutational meltdown” (Eigen, 2002; Anderson et al., 2004; Jonsson et al., 2005). On the other hand, an attenuation of the mutation rate converts the virus in a less adaptable species, which makes it unable to run ahead of the innate immune response or tu circumvent vaccination protection (Lauring and Andino, 2010).

The quasispecies model has been extensively studied in the literature and additional interesting mechanistic and biological properties have emerged (Wilke et al., 2001). Phase transitions have been reported in the context of the quasispecies model (Solé et al., 1999b; Solé, 2003). However, traditionally, the quasispecies model is defined in terms of differential equations, which only permit to study the smooth time evolution of the population mean of the quasispecies. The time series of internal fluctuations of the population can be taken into account by discretizing the system and by using difference equations (Solé et al, 1999a, 1999b; Solé, 2003). However, to date, the full potential of the time-series approach has not been exploited, and in particular, the model has not been studied from the perspective of early-warning signals. Our main hypothesis is that the statistical properties of time-series measurements can be an indicator for identifying critical regions in a dynamical sense. Two advantages of this approach are the following: 1) we overcome the need of having detailed mechanistic information of the systems; and 2) this view constitutes an immediate link with experimental data of complex systems. Early-warning signals identifying critical points such as tipping points are usually interpreted as “negative” signals, which indicate that the system is becoming increasingly unsustainable. In the present contribution however, we will identify early-warning signals for time series of the quasispecies model at the error threshold. Thus, in our case, early-warning signals rather serve as “positive” signals and indicate that the system has optimized its possibilities of survival. The early-warning signals we find are discontinuities in the behavior of distribution moments such as variance, kurtosis and skewness. Additionally, whereas the framework of early-warning signals usually focuses on a reddening of the power spectrum, which corresponds to an increase in strength of the autocorrelations, we argue that near the error threshold the power spectrum approximates a P(f)~1/f behavior, which implies a maximization of the range of the autocorrelations (Keshner, 1982; Fossion et al., 2010). Overall, these findings are in agreement with our suggestion in a previous article that a 1/f power spectrum can provide a fingerprint for the critical state of dynamical systems undergoing a phase transition (Landa et al., 2011).

Materials and methods

A discrete quasispecies model

The quasispecies model is a description of the process of the Darwinian evolution of self-replicating entities, and is used to model the competition between variety-increasing processes such as mutation as well as variety-reducing effects such as genetic drift and selection (Bull et al., 2005). Put simply, a quasispecies is a heterogeneous group or cloud of related genotypes that exist in an environment of nonzero mutation rate μ, where a large fraction of offspring is expected to contain one or more mutations relative to the parent. This is in contrast to the biologic concept of a species, which from an evolutionary perspective corresponds to a homogeneous group based on a single and stable genotype, of which most of the offspring will be genetically accurate copies. The quasispecies model is useful mainly in providing a qualitative understanding of the evolutionary processes of self-replicating macromolecules such as RNA or DNA or simple asexual organisms such as bacteria or viruses (viral quasispecies). In (Solé et al., 1999a, 1999b; Solé, 2003), a computational simplified quasispecies model is proposed, based on the Moran evolutionary model. Here, a set of rules dictates how a population of individuals, corresponding to different genotypes, evolves discretely in time. The total number of individuals in the population remains constant, but the number of individuals corresponding to specific genotypes can fluctuate. Also, subsequent generations are allowed to overlap. The genotypes are represented as bit-strings, that is, binary sequences of 1's and 0s with a length ν = 15, where each genotype is associated with a specific fitness. Based on this first approach, we propose an alternative recipe for the dynamical population model, see Fig. 1, based on the even more common Wright-Fisher evolutionary model, which also considers a constant population size, but consisting of non-overlapping generations. Thus, the populations at different time steps correspond to successive generations, so that there is no overlap between parent-sequences and offspring-sequences. In general, both models behave in similar fashion (Wakeley, 2006), however, our model computationally requires fewer time steps to obtain similar results in comparison with the Solé model. In the Supporting Information online, the algorithm of the model of Solé et al. is compared with the algorithm used in this contribution. Briefly, the rules used in this contribution are summarized as follows. For the initial population, i = 1, a number of S = 100 bit-strings is generated consisting of random 1's and 0s. The next population at time i + 1 is initially empty, and is iteratively filled with selected sequences from the previous population at time i. Selection is based on fitness, and can be simulated in a variety of ways. The simplest approach is to use a “needle in a haystack” fitness landscape, with an equal fitness for all sequences, except for a “master sequence,” which is the most fit. The master sequence is arbitrarily taken as the binary sequence consisting of only digits 1, see Fig. 2.

A candidate sequence is chosen at random from the population i, and is selected for copying to the population i + 1 with a probability of P = 100% if it is the master sequence, and with a probability of P = 5% if it is not the master sequence. The copy process involves a per-digit mutation rate μ, which is the probability for each digit to switch from 0 to 1, or from 1 to 0. These steps are repeated until the population at time i + 1 contains S = 100 bit strings. The time evolution of the population is followed over a total of N = 10⁴ iterations. The count of the number of master sequences in the population as a function of time, n(i), serves as a time series, and is found to be wildly fluctuating, because the number of master sequences from one population to the next can change drastically. For each mutation rate μ, the algorithm is repeated 100 times, and the percentage of trials P_m is computed where at least one master series has survived in the population at the end of the time evolution. Although intuition would expect a smooth, continuous decay of of this probability with mutation rate μ, there is actually a sharp decay close to a given critical rate, μ_crit, very much like a phase transition, where μ serves as a control parameter and P_m acts as the order parameter, see Fig. 3. This phenomenon is known as the error catastrophe (Bull et al., 2005). It describes the breakdown of hereditability and the transition from a quasispecies, that is, a heterogeneous set of mutant sequences, centered around the master sequence in sequence space at low mutation rates, toward a random set of strings, at high mutation rates, see Fig. 4.

Time-series analysis and early-warning signals

Complex dynamical systems, ranging from ecosystems to financial markets and the climate, can have critical thresholds–so-called tipping points–at which a sudden shift to a contrasting dynamical regime may occur. Predicting such critical points before they are reached, proves to be extremely difficult because models of complex systems are usually not accurate enough to reliably indicate where critical thresholds may occur. A more realistic approach may be “listening” to the system and to study its time series in a statistical way. It appears that certain generic symptoms may occur in time series of a wide class of systems as they approach a critical point, regardless of the differences in the details of each system. Early-warning signals can be defined in terms of statistical properties of the fluctuations of these time series. These properties can be static, i.e. related to the frequency and size distribution of the fluctuations, or dynamic, i.e. how the fluctuations follow each other up and how they are correlated in time, see e.g. (Scheffer et al., 2001, 2009, 2012; Drake and Griffen, 2010) for an overview. Most studies have concentrated on theoretical models, but the theoretical early-warning signals have been confirmed also in whole-ecosystem experiments (Carpenter et al., 2011), although their application in real systems might not always be obvious (Boettiger and Hastings, 2013). In this work, we will study the time series of the evolution of the master population n(i), and, more specically, whether its statistical properties differ for subcritical, critical and supercritical values of the mutation rate μ_crit. We are interested in the following early-warning signals:

Flickering and bimodality

Flickering is a phenomenon that can be seen when a critical point is approached, where different competing states coexist. If fluctuations are strong enough, then the system moves back and forth between the basins of attraction of the two alternative attractors (Berglund and Gentz, 2002). Such behavior is considered to be an early-warning signal, because the system may shift permanently to the alternative state if the underlying slow change in conditions persists, moving it eventually to a situation with only one stable state. Statistically, flickering can be observed in the frequency distribution of states as bimodality (reflecting the two alternative regimes), see e.g. the bimodal population distribution P_μ(k) for μ_crit in Fig. 4. Flickering also leads to increased variance and skewness in the frequency distribution.

Standard deviation and coefficient of variation

When a dynamical system is driven closer to its critical point, there is an increasing delay in recuperation rate of the system toward its steady-state, a phenomenon which is called critical slowing down. As a consequence, local fluctuations decay more slowly, can accumulate and grow larger in size. Thus the standard deviation of the state variable increases when the tipping point is approached (Carpenter and Brock, 2006), and we expect a similar behavior for the standard deviation of time series n(i) of the discrete quasispecies model near the error threshold. The standard deviation is defined as,

(1)\sigma =\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}{\displaystyle \sum _{i=\mathrm{1}}^N{\left(n(i)-\overline{n}\right)}^{\mathrm{2}}}},

where \mathrm{\Delta }n=n(i)-\overline{n} is a measure of the size of a fluctuation around the mean count \overline{n}. In the following, we will find that \overline{n} varies in an important way with the mutation rate μ. However, the standard deviation σ should always be understood in the context of the mean of the data, especially when data sets with different units or widely different means are compared. One solution is to consider the coefficient of variation (CV) (Bedeian and Mossholder, 2000),

(2)CV={\displaystyle \frac{\sigma }{\overline{n}}}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}{\displaystyle \sum _{i=\mathrm{1}}^N{\left({\displaystyle \frac{n(i)-\overline{n}}{\overline{n}}}\right)}^{\mathrm{2}}},}

of which the value is independent of the unit in which the measurement has been taken, so it is a dimensionless number. A disadvantage of CV is that when the mean value is close to zero, CV approaches infinity and is hence sensitive to small changes in the mean. We will indeed find that in our model the mean number of master-sequence counts \overline{n} approaches zero for large (supercritical) mutation rates μ. We will use as an alternative to CV a rescaling of all time series considered, so that they all fluctuate in the same dynamic range, between a predetermined minimum n_min and maximum n_max, and then we will compare the standard deviation σ of these rescaled time series. The standard deviation now quantifies the variability of the time series, independently from the mean value around which the time series fluctuates, which we will find to serve as a better early-warning signal of the critical error threshold.

Skewness and kurtosis

The asymmetry of the fluctuations is expected to increase near the tipping point, and hence the frequency distribution of the corresponding time series will be skewed. This asymmetry effect can be quantified by the skewness (Guttal and Jayaprakash, 2008), which is the third standardized moment of the frequency distribution,

(3)Skew(n)={\displaystyle \frac{{\displaystyle \sum _{i=\mathrm{1}}^N{\left(n(i)-\overline{n}\right)}^{\mathrm{3}}}}{N{\sigma }^{\mathrm{3}}}}.

In this work, we will find that also the kurtosis serves as an early-warning signal. The kurtosis or peakedness is the fourth standardized moment of the frequency distribution,

(4)Kurt(n)={\displaystyle \frac{{\displaystyle \sum _{i=\mathrm{1}}^N\left(n(i)-\overline{n}\right)}}{N{\sigma }^{\mathrm{4}}}}.

In this definition, the kurtosis and the skewness of the Gaussian distribution are Skew = 0 and Kurt = 3. A high kurtosis distribution has a sharper peak and longer, fatter tails, such that more of the variance is the result of infrequent extreme deviations; the opposite is true for a low kurtosis distribution, that has a more rounded peak and shorter, thinner tails, and where most of the variance comes from frequent modestly sized deviations (De Carlo, 1997). A zero skewness value indicates that the values are evenly distributed on both sides of the mean, typically but not necessarily implying a symmetric distribution. The skewness value can also be positive or negative, or even undefined. Qualitatively, a negative skewness indicates that the tail on the left side of the probability density function is longer than the right side and the bulk of the values (possibly including the median) lie to the right of the mean. A positive skew indicates that the he bulk of the values lie to the left of the mean (Doane and Seward, 2011).

Correlation strength

Because slowing down causes the intrinsic rates of change in the system to decrease, the state of the system at any given moment becomes more and more like its past state, with a consequent increase of the autocorrelation in the fluctuation pattern. The increase can be in the strength of the correlations and/or in their range, and both can be measured in the time domain by means of the autocorrelation function (Williams, 1997),

(5)C(\tau )={\displaystyle \frac{{\displaystyle \frac{\mathrm{1}}{N}}{\displaystyle \sum _{i=\mathrm{1}}^N\left(n(i)-\overline{n}\right)\left(n(i+\tau )-\overline{n}\right)}}{{\displaystyle \frac{\mathrm{1}}{N}}{\displaystyle \sum _{i=\mathrm{1}}^N{\left(n(i)-\overline{n}\right)}^{\mathrm{2}}}}},

where n(i)-\overline{n} measures the size of the fluctuation of the time series around its mean \overline{n} at time i, n(i+\tau )-\overline{n} corresponds with the size of the fluctuation after a time lag τ, and {\sigma }^{\mathrm{2}}=(\mathrm{1}/N){\displaystyle \sum {\left(n(i)-\overline{n}\right)}^{\mathrm{2}}}is the variance of the time series. The autocorrelation function C(τ) is normalized in such a way that C(0) = 1 and C(\infty )\to \mathrm{1}. The literature on early-warning signals focuses on the correlation strength. This strength can be quantified with the lag-1 autocorrelation coefficient C(1), which measures how strongly nearest neighbors are correlated on the average (i.e. neighbors that are a time lag τ = 1 apart). The correlation strength is found to rise drastically when the system approaches its critical point (Livina and Lenton, 2007; Biggs et al., 2009; Scheffer et al., 2009). The amount of autocorrelation strength in a time series can be estimated as well in the frequency domain by means of the power spectrum, which is the square of the Fourier transform of the time series (or the inverse Fourier transform of the autocorrelation function, according to the Wiener-Khinchin theorem). The power spectrum shows the frequencies that contribute to the time series, together with their relative importance. In the framework of early-warning signals, a “reddening” of the power spectrum is observed near the critical point, i.e. an increase in the contribution of low frequencies (Kleinen et al., 2003; Livina and Lenton, 2007; Biggs et al., 2009; Carpenter and Brock, 2010). Near the critical point, long-range correlations are supposed to emerge, giving rise to a fractal time series with a power spectrum that behaves as a power law (Bak et al., 1987; Fossion et al., 2010; Landa et al., 2011),

(6)P(f)\propto f^{\beta },

where β = 0 for a random time series (white noise), β<0 for a correlated time series (reddened spectrum) (Halley, 1996; Halley and Inchausti, 2004), and β>0 for an anti-correlated time series (bluish spectrum) (Peng et al., 1993). The spectrum is in general broadband, displaying its power-law behavior over several frequency scales. The spectral density exponent,

(7)\beta ={\displaystyle \frac{\mathrm{\Delta }\log P}{\mathrm{\Delta }\log f}},

can be seen as a measure of the correlation strength present in the time series, as a larger \left|\beta \right| implies that more power \mathrm{\Delta }\log P is concentrated in a smaller frequency range \mathrm{\Delta }\log f (Fossion et al., 2010). The autocorrelation function C(τ) and the power spectrum P(f) are strictly only defined for stationary time series (\beta \ge -\mathrm{1}) and theoretically can only be used to study linear correlations. However, C(τ) and P(f) are found to approximate well the behavior of the ensemble mean, if different trials of the same phenomenon are taken into account (Flandrin, 1989). On the other hand, in contrast to the autocorrelation function, the mutual information function is applicable also to nonlinear and non-stationary time series (Addison, 1997). In the Supporting Information online we present results for the mutual information function that confirm the results obtained with the autocorrelation function and the power spectrum.

Correlation range (memory)

Correlation strength as a parameter for criticality has a drawback. This parameter works well if a dynamical system is undergoing a transition in between two different non-correlated regimes: then correlations show up only at the critical point, and the correlation strength rises abruptly at this point and drops immediately after. If, instead, a dynamical system undergoes a transition from a non-correlated regime to a strongly correlated one, then the correlation strength would be rising monotonously for the whole transition without any distinctive behavior for the critical point. On the other hand, only at the critical point, long-range correlations are expected to emerge, and correlation range or memory could thus offer a good fingerprint for criticality. It can be shown analytically (Keshner, 1982) and numerically (Fossion et al., 2010) that random time series (\beta =\mathrm{0}) have zero memory, strongly correlated time series such as brownian noise (\beta =-\mathrm{2}) have litle memory, whereas the memory effect maximizes for moderately correlated times series such as 1/f noise (\beta =-\mathrm{1}). This memory or correlation range can be quantified with the autocorrelation function C(τ) of Eq. (5) by the time interval τ₀ for which the autocorrelations decay to zero, i.e. C(τ₀) = 0. Remark that some authors choose as memory the time intervals τ_1/e or τ_1-1/e for which the the autocorrelation function decays below the thresholds 1/e or 1-1/e (Rosenstein et al., 1993; Hausdorff et al., 1999). In this article, we will find that another way to quantify the memory effect of a time series is the lag-<t>autocorrelation coefficient, i.e. the average correlation strength of neighbors that are a time lag τ =<t>apart, where<t>is the mean duration of fluctuations in the time series. The average fluctuation duration can be estimated in the following way,

(8)\langle t\rangle ={\displaystyle \frac{N}{\langle f\rangle }},

where<f>is the average fluctuation rate over the whole broadband frequency range of the Fourier power spectrum, from f_min to f_max (Rosenstein et al., 1993; Fossion et al., 2010),

(9)\langle f\rangle ={\displaystyle \frac{{\displaystyle \sum _{f=f_{\min }}^{f_{\max }}f\cdot P(f)}}{{\displaystyle \sum _{f=f_{\min }}^{f_{\max }}P(f)}}}.

In the following, we will find that the memory measures {\tau }_{\mathrm{0}} and C(<t>) have a very similar behavior.

Results

Time series: time and frequency analysis

In what follows, we will apply the method of early-warning signals to the time series n(i) generated with the discrete quasispecies model described above for different values of the mutation rate parameter μ. As mentioned in the introduction, in the present case of quasispecies time-series, the early-warning signals identified appear to indicate the optimization of the systems’s survival possibilities. In Fig. 5, we present, for each of the three regions of the phase diagram, for μ<0.12 (subcritical regime), for \mathrm{0.12}\le \mu \le \mathrm{0.15} (critical regime) and for 0.15<μ (supercritical regime), a representative time series n(i), and the corresponding power spectrum P(f), autocorrelation function C(τ), and histogram P(n).

For small (subcritical) values of the mutation-rate parameter (shown is the case for μ = 0.05, left panels of Fig. 5), the time series n(i) is continuous and corresponds to random fluctuations, also called white noise, which is confirmed by the the power spectrum (\beta \approx \mathrm{0}). Shown are the power spectrum for a single trial of the time series, and the mean power spectrum averaged over the 100 trials of the ensemble for μ = 0.05. The mean power spectrum shows the average behavior of the ensemble more clearly. The correlation functions for all realizations of the time series decay very steeply and approximately correspond to a Dirac delta function, C(τ)\approxδ(τ-0). The correlation range or memory of the time series is very small, {\tau }_{\mathrm{0}}\approx {\mathrm{10}}^{\mathrm{0.5}}\approx \mathrm{3} iterations. The histogram P(n) of the time series corresponds to a symmetric and Gaussian distribution. All the statistical measures agree that in this subcritical regime, the time series n(i) corresponds to random Gaussian fluctuations around a mean \overline{n}.

For intermediate (critical) values of the mutation rate (shown is the case for μ = 0.13, middle panels of Fig. 5), the time series is intermittent, and bunches with nonzero master-sequence counts (n(i)\ne \mathrm{0}) alternate with time intervals with zero master-sequence counts (n(i) = 0). The intermittent time series describes the situation where a lineage of master sequences survives for several hundreds or thousands of generations in the population, and then dies out, to resurrect many iterations later. This intermittency can be a manifestation of the flickering behavior that is proposed as one of the early-warning signals. Intermittent systems are known sometimes to lead to 1/f behavior (Manneville, 1980; Procaccia and Schuster, 1983; Schuster and Just, 2005). In the present case, the power spectrum consists of two flat regions and two inclined regions, but overall approximates a 1/f behavior (Azhar and Gopala, 1992; Hausdorff and Peng, 1996). The autocorrelation function maximizes its memory in this critical regime, {\mathrm{10}}^{\mathrm{2.5}}\le {\tau }_{\mathrm{0}}\le {\mathrm{10}}^{\mathrm{3.0}}(300 to 1000 iterations). The histogram of the time series is bimodal, corresponding to a normal distribution and a superposed delta function. As in the subcritical regime, the normal distribution corresponds to the fluctuating behavior of the agglomerations or bunches. We will find the delta-function distribution to be characteristic of the supercritical regime; here the delta function corresponds to the contribution of the zero counts (n = 0) in the time series. Bimodal distributions have been suggested as early-warning signals, and a similar histogram to the present one has been described in (Carpenter and Brock, 2006).

For high (supercritical) values of the mutation rate (shown is the case for μ = 0.16, right panels of Fig. 5), the time series is a peak train, with rather strong correlations at high frequencies or short time intervals (\beta \approx -\mathrm{1.5}), and no correlations at low frequencies or large time intervals (\beta \approx \mathrm{0}), as can be seen in the power spectrum. The autocorrelation function decays rather quickly, with a memory effect of {\mathrm{10}}^{\mathrm{1.5}}\le {\tau }_{\mathrm{0}}\le {\mathrm{10}}^{\mathrm{2.0}} (30 to 100 iterations), which is an order of magnitude less than in the critical regime, but an order of magnitude more than in the subcritical regime. The histogram of the time series is dominated by a peak at n = 0, whereas higher counts n\ge \mathrm{1} are almost insignificant.

Correlations: strength vs. range

One of the proposed early-warning signals is the increased autocorrelation near the critical point, and literature focuses primarly on the increased autocorrelation strength. In Fig. 5, a reddening of the power spectra can be appreciated with increasing μ, with \beta \approx \mathrm{0} for subcritical mutations rates, \beta \approx -\mathrm{1} for critical mutation rates, and \beta \approx -\mathrm{1.5} for supercritical mutation rates and high frequencies. Here, correlation strength as quantified by the reddening of the power spectrum (more negative β) thus fails to indicate the critical region. This is confirmed by the behavior of the other measure of the autocorrelation strength, the lag-1 autocorrelation coefficient C(1), see Fig. 6(A), where the evolution of its minimum, mean and maximum is given over an ensemble of 100 trials, as a function of the mutation rate μ. It can be seen that C(1) is very small in the subcritical region, and that it rises slowly with μ. In the critical region, the correlation strength rises rapidly, and reaches a maximum C(\mathrm{1})\approx \mathrm{0.9} near μ = 0.13. In the supercritical region, the correlation strength drops slowly with μ, but remains rather large for the supercritical mutation rates considered. For zero mutation rate, μ = 0, the correlation strength is rather high,C(\mathrm{1})\approx \mathrm{0.6}, because without mutations, the model looses a significant part of its randomness.

In the case of a transition between two random regimes, correlation strength might serve as a good fingerprint of criticality, but this is not the case here, and we will try correlation range (memory) instead. In Fig. 6(B), we show the lag-<t>autocorrelation coefficient C(<t>) in function of the mutation rate μ. It can be seen that the behavior is much more abrupt, with a steep rise when entering the critical región \mathrm{0.12}\le \mu \le \mathrm{0.15} and a steep downfall on leaving. Also, C(<t>) corrects the anomalous behavior of C(1) at μ = 0. An argument for the statement that C(<t>) quantifies correlation range rather than strength can be found in comparison with Fig. 6(C), where the memory parameter τ₀ is shown as a function of mutation rate μ. The behavior of the parameters C(<t>) and τ₀ is very similar, in the sense that they behave very abruptly when entering and leaving the critical region. In the Supplementary Information online, one can appreciate that the correlation range or memory of the mutual information function also maximizes for critical values of the mutation rate.

Moments of the time series

In Fig. 7(A), we present the standard deviation σ as a function of the mutation rate μ, averaged over 100 trials of the time series. In the subcritical region, σ rises with μ until it reaches a maximum for\mu \approx \mathrm{0.04}, after which σ drops again. In the critical region, σ peaks abruptly for μ = 0.13. In the supercritical region, σ falls off quickly to zero. In the present case, the standard deviation used as an early-warning signal could lead to misleading conclusions, as the local maximum of σ near μ = 0.04 could be interpreted as a sign of criticality. However, the standard deviation from time series n(i) at different mutation rates μ cannot be compared directly, because the mean \overline{n}, around which the time series fluctuates, varies with the mutation rate μ, as can be appreciated in the time series of Fig. 5. In the inset of Fig. 7(A), the time series n(i) were first rescaled, so that they all fluctuate in the same dynamic range between n_min = 0 and n_max = 100, and then the standard deviation was calculated, so that σ is independent of the mean of the time series. Now, in the subcritical region, σ rises slowly with μ, rises very steeply in the critical region, maximizing at μ = 0.13, to drop again abruptly for larger μ. Finally, in the supercritical region, σ approaches 0. We find that the rescaled standard deviation indeed seems to function as an early-warning signal, and that its sudden rise in the critical region is due to the impending phase transition.

Figures 7(B) and 7(C) show, as a function of the mutation rate μ, the evolution of skewness Skew and kurtosis Kurt, averaged over 100 trials of the time series. It can be seen that for unrealistically small or large values of the mutation rate, Skew and Kurt have extreme values. This is understandable, because of the frequency distributions behaving as delta functions P(n) = δ(n–100) and P(n) = δ(n–0), for μ = 0 and μ = 0.20, respectively. For more moderate and realistic mutation rates, 0<μ<0.20, the probability distribution relaxes smoothly to a symmetric Gaussian distribution with symmetric fluctuations around the mean, so that also the kurtosis and the skewness behave smoothly and monotonously and Kurt\to \mathrm{3} and Skew\to \mathrm{0}. However, in the critical region, and specifically for μ = 0.130, Skew and Kurt behave irregularly (a sudden change of the monotonous behavior and/or change of sign), exactly as would be expected for an early-warning signal, see the insets of Figs. 7(B) and 7(C), which show the irregular behavior in more detail.

Discussion

To study natural phenomena, one approach is to use mechanistic models where the main problem is to define the internal interactions of the system and then to fine-tune the parameters in such way that the dynamics of the system is reproduced as well as possible. In complex systems, composed of many interacting components, this might be a hazardous task. A more “humble” approach might be to “listen” to the system, i.e. to carry out a thorough statistical study of the associated time series, which in principle permits to understand the underlying dynamics. This is the philosophy followed in the framework of early-warning signals where the focus is on complex systems, ranging from ecosystems to financial markets and the climate, that can be in multiple steady-states. When an external control parameter is varied, the system moves from an original state to an alternative state. It is possible that the behavior of the system is nonlinear as a function of the control parameter, in which case the transition can be sudden and drastic, even if the variation of the control parameter is smooth. In this case, one speaks about the “collapse” of the system. Time series, associated to the system, exhibit characteristic and universal fingerprints in their statistics when the system is near the transitional point, and hence can be used to predict the impending change or collapse. Early-warning signals are usually interpreted as negative signals, that predict that the system is at the brink of unsustainability (such as catastrophic shifts in wildlife populations, a systemic market crash, an abrupt shift in ocean circulation, etc.) and that measures need to be taken urgently to revert the system back to its original stable state.

Self-organizing systems are a particular type of complex systems that are supposed to converge spontaneously toward a critical state in between alternative phases, and where also long-range correlations emerge. Flocks of starling birds can serve as an example, since the birds can tune their individual behavior in between two extremes: imitating perfectly their neighbors, or following exactly their own free will. In the first case, the flock would behave as a rigid flying object, whereas in the second case the flock would resemble a random, ideal gas of bird-particles. Only in between these extremes the flock behaves as a complex organism. We have illustrated these concepts by studying time series in the context of a discrete quasispecies model, where there is competition between the random force of mutation and the deterministic force of selection. When the selection force is dominant (i.e. small mutation rates), the population is a heterogeneous collection of genotypes centered around the most fit genotype, which is called the “master sequence.” Genetic information is passed on between subsequent generations, but the population composition is rather static and not very adaptable to the changing circumstances. When the mutation force is dominant (i.e. high mutation rates), the population is a random mixture of genotypes and its composition is constantly changing, but transfer of genetic information between the generations ceases to exist. There is a critical range of moderate mutation rates where robustness and evolvability are simultaneously optimized, and which corresponds to the quasispecies regime. It has been suggested that RNA viruses self-organize their mutation rate into this regime, where they change quickly enough to escape recognition by the immune system without however losing hereditability. One of the main results of this work is that we identified for time series associated to the critical quasispecies regime the early-warning signals previously proposed in literature. Whereas usually early-warning signals are interpreted in a “negative” way, i.e. a collapse or unsustainability of the system, in the present context of a self-organized complex system the early-warning signals serve as “positive” signals, i.e. they are the fingerprint of a complex system that has optimized its functionality within its own boundaries.

One of the early-warning signals proposed in the literature is the increase in autocorrelation, where the focus is on an increase of the autocorrelation strength. In this contribution, we have found that it might be more useful to focus on an increase in the correlation range, and we used two different quantifications for the correlation range. More in particular, the correlation range maximizes for critical values of the control parameter and drops abruptly for parameter values outside of the critical region. The standard deviation σ is another important early-warning signal. However, one cannot compare the standard deviation of different time series if their mean values differ in a significant way, which is exactly the case in this work for the time series n(i) as a function of the mutation rate μ. This problem can be solved if the time series are first rescaled, so that all time series fluctuate in the same dynamical range, before the standard deviation is calculated. Alternatively, one can use the coefficient of variation, which is the standard deviation divided by the mean, if at least the mean is much larger than zero. In literature, asymmetry of the frequency distribution and skewness has been suggested as an early-warning signal. Here we find that also the kurtosis or peakedness can contribute as an early-warning signal.

References

[1] Addison P (1997). Fractals and chaos: An illustrated course. Bristol and Philadelphia: Institute of Publishing

[2] Aldana M, Balleza E, Kauffman S, Resendiz O (2007). Robustness and evolvability in genetic regulatory networks. J Theor Biol, 245(3): 433-448

[3] Anderson J P, Daifuku R, Loeb L A (2004). Viral error catastrophe by mutagenic nucleosides. Annu Rev Microbiol, 58(1): 183-205

[4] Azhar M A, Gopala K (1992). Clustering Poisson process and burst noise. Jpn J Appl Phys, 31(Part 1, No. 2A): 391-394

[5] Bak P (1996). How Nature works: the science of self-organized criticality. New York: Springer-Verlag.

[6] Bak P, Tang C, Wiesenfeld K (1987). Self-organized criticality: An explanation of the 1/f noise. Phys Rev Lett, 59(4): 381-384

[7] Ballerini M, Cabibbo N, Candelier R, Cavagna A, Cisbani E, Giardina I, Lecomte V, Orlandi A, Parisi G, Procaccini A, Viale M, Zdravkovic V (2008). Interaction ruling animal collective behavior depends on topological rather than metric distance: evidence from a field study. Proc Natl Acad Sci USA, 105(4): 1232-1237

[8] Bedeian A G, Mossholder K W (2000). On the use of the coe.cient of variation as a measure of diversity. Organ Res Methods, 3(3): 285-297

[9] Berglund N, Gentz B (2002). Metastability in simple climate models: pathwise analysis of slowly driven langevin equations. Stoch Dyn, 2(03): 327-356

[10] Biebricher C K, Eigen M (2005). The error threshold. Virus Res, 107(2): 117-127

[11] Biggs R, Carpenter S R, Brock W A (2009). Turning back from the brink: detecting an impending regime shift in time to avert it. Proc Natl Acad Sci USA, 106(3): 826-831

[12] Boettiger C, Hastings A (2013). Tipping points: From patterns to predictions. Nature, 493(7431): 157-158

[13] Bull J J, Meyers L A, Lachmann M (2005), Quasispecies made simple. Plos Comput Biol, 1: e61 (0450-0460).

[14] Carpenter S R, Brock W A (2006). Rising variance: a leading indicator of ecological transition. Ecol Lett, 9(3): 311-318

[15] Carpenter S R, Brock W A (2010). Early warnings of regime shifts in spatial dynamics using the discrete fourier transform. Ecosphere, 1(5): 10

[16] Carpenter S R, Cole J J, Pace M L, Batt R, Brock W A, Cline T, Coloso J, Hodgson J R, Kitchell J F, Seekell D A, Smith L, Weidel B (2011). Early warnings of regime shifts: a whole-ecosystem experiment. Science, 332(6033): 1079-1082

[17] Cavagna A, Cimarelli A, Giardina I, Parisi G, Santagati R, Stefanini F, Viale M (2010). Scale-free correlations in starling flocks. Proc Natl Acad Sci USA, 107(26): 11865-11870

[18] Crotty S, Cameron C E, Andino R (2001). RNA virus error catastrophe: direct molecular test by using ribavirin. Proc Natl Acad Sci USA, 98(12): 6895-6900

[19] DeCarlo L T (1997). On the meaning and use of kurtosis. Psychol Methods, 2(3): 292-307

[20] Doane D P, Seward L E (2011). Measuring skewness: A forgotten statistic? J Stat Educ, 19: 1-18

[21] Drake J M, Griffen B D (2010). Early warning signals of extinction in deteriorating environments. Nature, 467(7314): 456-459

[22] Eigen M (2002). Error catastrophe and antiviral strategy. Proc Natl Acad Sci USA, 99(21): 13374-13376

[23] Flandrin P (1989). On the spectrum of fractional brownian motion. IEEE Trans Inf Theory, 35(1): 197-199

[24] Fossion R, Landa E, Stránský P, Velázquez V, López Vieyra J C, Garduño, García D, Frank A (2010). Scale invariance as a symmetry in physical and biological systems: Listening to photons, bubbles and heartbeats. In: Benet L, Hess P O, Torres J M, Wolf K B, editors, Symmetries in Nature: Symposium in Memoriam Marcos Moshinsky (Cuernavaca, Mexico, 7 - 14 August), New York: AIP Conf. Proc., volume 1323: 74-90

[25] Guttal V, Jayaprakash C (2008). Changing skewness: an early warning signal of regime shifts in ecosystems. Ecol Lett, 11(5): 450-460

[26] Halley J M (1996). Ecology, evolution and 1/f -noise. Trends Ecol Evol, 11(1): 33-37

[27] Halley J M, Inchausti P (2004). The increasing importance of 1/f-noises as models of ecological variability. Fluct Noise Lett, 4(02): R1-R6

[28] Hausdorff J M, Peng C K (1996). Multiscaled randomness: A possible source of 1/f noise in biology. Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics, 54(2): 2154-2157

[29] Hausdorff J M, Zemany L, Peng C K, Goldberger A L (1999). Maturation of gait dynamics: stride-to-stride variability and its temporal organization in children. J Appl Physiol, 86(3): 1040-1047

[30] Jonsson C B, Milligan B G, Arterburn J B (2005). Potential importance of error catastrophe to the development of antiviral strategies for hantaviruses. Virus Res, 107(2): 195-205

[31] Kauffman S (1995). At home in the universo: The search for the laws of self-organization and complexity. New York, Oxford: Oxford University Press

[32] Keshner M S (1982). 1/f noise. Proc IEEE, 70(3): 212-218

[33] Kleinen T, Held H, Petschel-Held G (2003). The potential role of spectral properties in detecting thresholds in the earth system: application to the thermohaline circulation. Ocean Dyn, 53(2): 53-63

[34] Landa E, Morales I O, Fossion R, Stránský P, Velázquez V, Vieyra J C, Frank A (2011). Criticality and long-range correlations in time series in classical and quantum systems. Phys Rev E Stat Nonlin Soft Matter Phys, 84(1 Pt 2): 016224

[35] Lauring A S, Andino R (2010). Quasispecies theory and the behavior of RNA viruses. PLoS Pathog, 6(7): e1001005

[36] Lenski R E, Barrick J E, Ofria C (2006). Balancing robustness and evolvability. PLoS Biol, 4(12): e428

[37] Livina V N, Lenton T M (2007). A modified method for detecting incipient bifurcations in a dynamical system. Geophys Res Lett, 34(3): L03712

[38] Manneville P (1980). Intermittency self-similarity and 1/f spectrum in dissipative dynamical systems. J Phys (Paris), 41(11): 1235-1243

[39] Miramontes O (1995). Order-disorder transitions in the behavior of ant societies. Complexity, 1: 50-60

[40] Peng C K, Mietus J, Hausdorff J M, Havlin S, Stanley H E, Goldberger A L, (1993). Long-range anticorrelations and non-Gaussian behavior of the heartbeat. Phys Rev Lett, 70(9): 1343-1346

[41] Procaccia I, Schuster H (1983). Functional renormalization-group theory of universal 1/f noise in dynamical systems. Phys Rev A, 28(2): 1210-1212

[42] Rosenstein M T, Collins J J, Luca C J D (1993). A practical method for calculating largest lyapunov exponents from small data sets. Physica D, 65(1-2): 117-134

[43] Scheffer M, Bascompte J, Brock W A, Brovkin V, Carpenter S R, Dakos V, Held H, van Nes E H, Rietkerk M, Sugihara G (2009). Early-warning signals for critical transitions. Nature, 461(7260): 53-59

[44] Scheffer M, Carpenter S, Foley J A, Folke C, Walker B (2001). Catastrophic shifts in ecosystems. Nature, 413(6856): 591-596

[45] Scheffer M, Carpenter S R, Lenton T M, Bascompte J, Brock W, Dakos V, van de Koppel J, van de Leemput I A, Levin S A, van Nes E H, Pascual M, Vandermeer J (2012). Anticipating critical transitions. Science, 338(6105): 344-348

[46] Schuster H, Just W (2005) Deterministic chaos: an introduction. Weinheim: Wiley-VCH Verlag GmbH & Co. KGaA, 4th edition.

[47] Solé R V (2003). Phase transitions in unstable cancer cell populations. Eur Phys J B, 35: 117-123

[48] Solé R V, Ferrer R, González-García I, Quer J, Domingo E (1999b). Red queen dynamics, competition and critical points in a model of RNA virus quasispecies. J Theor Biol, 198(1): 47-59

[49] Solé R V, Manrubia S C, Benton M, Kauffman S, Bak P (1999a). Criticality and scaling in evolutionary ecology. Trends Ecol Evol, 14(4): 156-160

[50] Solé R V, Miramontes O, Goodwin B C (1993a). Oscillations and chaos in ant societies. J Theor Biol, 161(3): 343-357

[51] Solé R V, Miramontes O, Goodwin B C (1993b) Emergent behaviour in insect societies: Global oscillations, chaos and computation. In: Haken H, Mikhailov A, editors, Interdisciplinary Approaches To Nonlinear Complex Systems, Springer Series in Synergetics. Berlin Heidelberg: Springer-Verlag, volume 62:, 77-88

[52] Wakeley J (2006). Coalescent theory: An introduction. Greenwood Village, Colorado: Roberts & Co

[53] Wilke C O, Wang J L, Ofria C, Lenski R E, Adami C (2001). Evolution of digital organisms at high mutation rates leads to survival of the flattest. Nature, 412(6844): 331-333

[54] Williams G P (1997). Chaos theory tamed. Washington D.C.: Joseph Henry Press

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg