# Evolution dynamics in terraced NK landscapes

## Abstract

We consider populations of agents evolving in the fitness landscape of an extended $`\mathrm{𝐍𝐊}`$ model with a tunable amount of neutrality. We study the statistics of the jumps in mean population fitness which occur in the ‘punctuated equilibrium’ regime and show that, for a wide range of landscapes parameters, the number of events in time $`t`$ is Poisson distributed, with the time parameter replaced by the logarithm of time. This simple log-Poisson statistics likewise describes the number of records in any sequence of $`t`$ independently generated random numbers. The implications of such behavior for evolution dynamics are discussed.

PACS numbers: 87.10.+e, 87.15.Aa, 05.40.-a   

Introduction. Evolutionary dynamics can be described as a stochastic process unfolding in a ‘fitness landscape’, a process which can be simulated by means of genetic algorithms . The dynamical behavior of such algorithms is controlled by a number of parameters as e.g. the population size, the rate of mutation and the strenght of selection.

A pervasive dynamical regime is the so-called ‘punctuated equilibrium’ or ‘epochal behavior’, where relevant measures of evolution, as e.g. the mean fitness, remain constants for long periods during which the fitness distributions of the individuals is strongly peaked. Occasionally, a fitter mutant appears and quickly ‘takes over’ the population (see e.g. Fig. 1). Real experiments performed on bacterial colonies evolving in a controlled environment have shown that the fitness and cell size increase at a decelerating rate . A similar slowing down is discussed by Kauffman in the ‘long jump’ dynamics of the $`\mathrm{𝐍𝐊}`$ model, while Aranson et al. find a logarithmic growth of the average fitness for quasi-species evolving in a rugged fitness landscape.

In a macroevolutionary context, Raup and Sepkoski suggested that the noticeable decay of the extinction rate might stem from the properties of an underlying optimization process. This idea was taken up in the ‘reset’ model, where jumps in the average fitness of populations are linked to fitness record achieved during evolution. Such a link between small and grand scale evolution is rather controversial: If macroevolutionary events are mainly driven extrinsically, e.g. by meteorite impacts, any patterns in the fossil record must ultimately stem from the mechanics of celestial bodies in the solar system. Conversely, if, as recently proposed by several authors, the fossil record is mainly shaped by complex interactions within the biotic system, macro-evolutionary patterns must emerge from population dynamics.

Our main interest lies in the statistical properties of the evolutionary jumps underlying fitness changes. Since the role of neutral mutations for evolution is well established, and since punctuated equilibria exist even in the absence of local fitness maxima, we chose to study evolution on ‘terraced’ landscapes with a tunable degree of neutrality.

We find that, on average, the number of jumps taking place in time $`t`$ grows proportionate to $`\mathrm{log}t`$, and that the rate of events consequently decays as $`1/t`$. Secondly and most importantly, we show that the record dynamics  provides a reasonable description of microevolution, with some limitations which are also outlined. Thirdly, we emphasize that power-law decays generically characterize the correlation functions of time series generated by a record driven dynamics.

Method. Each ‘genome’ in a population constitutes a point in an abstract configuration space usually called a fitness landscape . To construct such a landscape we use an elegant prescription due to Kauffman, the widely known $`\mathrm{𝐍𝐊}`$ model: We represent genomes as strings of $`𝐍`$ bits $`𝐱=(x_1,x_2,\mathrm{}x_N)`$, each being either $`0`$ or $`1`$. The fitness $`F(𝐱)`$ of configuration $`𝐱`$ is defined as

$$F(𝐱)=\frac{1}{𝐍}\underset{i=1}{\overset{𝐍}{}}f_i,$$
(1)   
where the contribution from site $`i`$, $`f_i`$, is a random function depending on $`x_i`$ and $`𝐊`$ other $`x_l`$’s. More precisely, $`f_i`$ is a random function of $`2^{𝐊+1}`$ arguments with values uniformly distributed in $`(0,1]`$. We let $`\mu _f`$ and $`\sigma _f`$ be the average and spread of the distribution from which the $`f_i`$ values are generated.

If $`𝐊`$ is zero, the change in fitness due to the change of one $`x_i`$ (a point mutation) is of order $`1/N`$, and the landscape may be regarded as smooth. By way of contrast, when $`𝐊=𝐍\mathrm{𝟏}`$ a single point mutation changes all the $`f_i`$’s, and the landscape becomes ‘rugged’. Intermediate cases correspond, of course, to intermediate $`𝐊`$ values.

Our version of the $`\mathrm{𝐍𝐊}`$ model is modified in two respects. Firstly, the sum in Eq.1 is shifted by $`\mu _f`$ and scaled by $`\sqrt{𝐍}`$ rather than $`N`$. This ensures that the distribution of fitness values keeps the same variance for any value of $`N`$. More importantly, we introduce tunable neutrality by discretizing the fitness values into ‘terraces’, according to the formula:

$$F^{}(𝐱)=\frac{\mathrm{nint}(j\sqrt{𝐍}[F(𝐱)\mu _f]+j\mu _f)}{j}.$$
(2)   
Here $`\mathrm{nint}`$ stands for the nearest integer function, and $`j`$ denotes the number of terraces. For small $`j`$ there are few broad terraces, while for $`j\mathrm{}`$ the fitness values approach a continuum and the original $`\mathrm{𝐍𝐊}`$ model is regained. Accordingly, one expects that, as seen e.g. in Fig. 2, the effect of varying $`j`$ should be stronger for small values of $`j`$.

In the simulations, a configuration $`𝐱`$ is cloned with a probability $`p\mathrm{exp}(\beta F(𝐱))`$. With probability $`u<<1`$, the cloned string undergoes a one-point mutation at a randomly chosen locus. Finally, a random and fitness independent deletion mechanisms is applied, which keeps the population size fluctuating around a fixed average $`n`$. Both the generation $`/`$ mutation part and the deletion part of the algorithm are performed sequentially. Subsequently, the information about the new fitness distribution is incorporated in the cloning probabilities. This entire process counts as one update and defines the unit of time.

Punctuated equilibrium dynamics requires that a fitter mutant be able to survive and spread in the population on a time scale short compared to the inverse mutation frequency. Therefore, $`u`$ should not be too high and $`\beta `$ should not be too low. Within these constraints $`u`$ should be as high as possible, in order to have a good statistics within the time window of the computation. Also, too high a $`\beta `$ value quenches the dynamics completely. These design consideration lead to the values of $`u`$ and $`\beta `$ used throughout the calculations.

A concise description of the dynamics is provided by the average $`\overline{f}(t)`$ of the distribution of fitness values through the population. Such a trajectory is shown in the main panel of Fig. 1 to consist of a number of flat plateaus separated by rather well defined jumps. The number $`m(t)`$ of jumps occurring in the interval $`[0,t)`$ is a stochastic process whose distribution can be sampled by repeating the simulations or, equivalently, by considering an ensemble of landscapes, where different trajectories are generated by independently updating each system. The ensemble average and variance of $`m`$ are denoted by $`E(m)`$ and $`\sigma ^2(m)`$ respectively, while $`a_m(t)`$ and $`v_m(t)`$ are the corresponding estimators. The notation describing the inputs and results of our simulations is summarized in Table 1.

Results. The basic qualitative features of the data are expressed by Fig. 1: Its left panel shows the mean fitness $`\overline{f}`$ of a single population on a semilogarithmic scale. The right panel details the behavior close to the evolutionary jump at $`t3\times 10^4`$ by depicting the distribution of fitness values right before, during and right after the jump. Nearly all strings have the same fitness in the initial and final situations, while the transition stage features two different fitness values. Punctuated equilibrium behavior and a very peaked fitness distribution are widely found in previous studies as well as in our simulations.

For a more quantitative data analysis, the statistical properties of the jumps and their associated waiting times must be studied. We let $`m_i(t)`$ be the number of jumps occurred at time $`t`$ in trajectory $`i`$, and consider the sample average:

$$a_m(t)=\frac{_{i=1}^sm_i(t)}{s},$$
(3)   
and the sample variance:

$$v_m(t)=\frac{_{i=1}^s(m_i(t)a_m(t))^2}{s1}$$
(4)   
as functions of time. To calculate error bars on $`a_m`$ and $`v_m`$ we need the corresponding standard deviations, which are $`\sigma (m)/\sqrt{s}`$ and $`\sqrt{(\sigma ^2(m^2)+8E^2(m)E(m^2)4E(m)(E^3(m)+E(m^3)))}/\sqrt{s}+𝒪(1/s)`$. The latter relation results from straightforward but rather tedious algebra. To lowest order in $`s^{1/2}`$ one may now replace the moments of $`m`$ with their sample estimators. In general, the relative errors on various quantities of interest are of order $`1/\sqrt{s}10\%`$.

Fig. 2 shows the average number of jumps, for a number of different parameters values, as a function of time. $`E(m)`$ is seen to grow logarithmically, with a strong $`𝐊`$ dependence slope, for ‘short’ log-times. The leveling off noticed at large times stems from the fact that, as the fitness increases, fitness improvements become progressively smaller. Eventually, they get lost in the noise, rather than triggering a jump. At this point record statistics and fitness evolution must part company.

For large $`t`$, the probability of $`n`$ records in a sequence of $`t`$ independently drawn random numbers is given by:

$$P_t(n)=\frac{(\mathrm{log}\lambda t)^n}{n!}t^\lambda .$$
(5)   
This is a Poisson distribution with $`\mathrm{log}t`$ replacing the time argument. The strength parameter $`\lambda `$ describes the possibility that many searches for records take place independently and in parallel and/or the situation where records remain undetected. In our systems we observed that several records were indeed lost in the noise, and did not trigger any evolutionary event.

A mathematically equivalent description of record dynamics is provided by the distribution of the variables

$$\mathrm{\Delta }_k=\mathrm{log}t_k\mathrm{log}t_{k1}=\mathrm{log}(t_k/t_{k1}).$$
(6)   
It follows from standard arguments that in a (log) Poisson process these $`\mathrm{\Delta }_k`$ are independent and have the common distribution:

$$P(\mathrm{\Delta }>x)=\mathrm{exp}(\lambda x).$$
(7)   
The empirical distribution of the $`\mathrm{\Delta }_k`$’s is shown on a semilogarithmic scale in Fig. 3 for $`𝐊=31`$ and for differing degrees of terracing. The decay of $`P(\mathrm{\Delta }>x)`$ seems quite well described by an exponential for $`x5`$, with the tail of the distribution falling off more rapidly, likely due to the inevitably poor sampling of large $`x`$ values in a finite time simulation. Banning the effect of the deviations from pure record statistics, $`\lambda `$ is, by Eqs. 5 and 7, the slope of $`a_m(t)`$ vs. $`\mathrm{log}t`$ as well as minus the slope of $`\mathrm{log}P(\mathrm{\Delta }>x)`$ vs $`x`$. We calculated these slopes from the data in both ways (cutting off the tails of the data), and obtained the following results. $`j=10^2,\lambda =0.54(0.34)`$; $`j=10^3,\lambda =0.49(0.39)`$; $`j=10^4,\lambda =0.58(0.44)`$; $`j=10^5,\lambda =0.54(0.35)`$; $`j=10^6,\lambda =0.62(0.47)`$; and $`j=10^9,\lambda =0.63(0.39)`$. The figures in parentheses stem from Fig. 3. There is a rough agreement, but certainly also considerable scatter in these data, with the log-wait time type of analysis yielding systematically lower figures. At this stage it is unclear whether the non-monotonic dependence of $`\lambda `$ on $`j`$ seen in Fig. 3 is a real effect - or just due to a combination of statistical fluctuations and systematic deviations from the ideal log-Poisson behavior.

Summarizing the results from Figs. 2 and 3, it appears that the value of $`𝐊`$ very strongly affects the average slope of the curves (i.e. the value of $`\lambda `$). The degree of terracing might also have an effect on the slope, albeit a minor one.

The independence of the different $`\mathrm{\Delta }_k`$’s implied by the record statistics was tested by calculating the correlation coefficients $`C_{ij}`$ between $`\mathrm{\Delta }_i`$ and $`\mathrm{\Delta }_j`$. In practice, we checked for $`C_{12}`$, $`C_{23}`$ and $`C_{13}`$. As expected, the highest degree of correlation was found for the relatively smooth landscape with $`𝐊=7`$. In this case, the $`C`$ values were close to $`0.4`$. For $`𝐊=31`$ the correlation coefficients were of the order of $`0.1`$, i.e. of the same order as the statistical sampling error.

To conclude the description of our data we plot in Fig. 4 the estimated variance $`v_m(t)`$ versus the estimated average $`a_m`$, for $`𝐊=31`$ and a number of terrace values. For a perfect agreement with the log-Poisson distribution the points should lie on a straight line of slope $`1`$. This is close to the observed behavior, except for the highest values (where, on the other hand, the statistics is poorest). A similar plot for $`𝐊=15`$ and $`𝐊=7`$ shows a systematic deviation from a straight line, with considerably less variance than in a purely random case. Additional details on the simulations and on the genetic algorithm utilized to produce them can be found in Ref. .

In summary, the dynamics of our evolutionary model is time inhomogeneous stochastic process, with a rate of events falling off as $`1/t`$. The log-Poisson statistics describes the data best for landscapes with large $`𝐊`$, with or without terraces. In this respect terraces have a minor effect on the dynamics.

Discussion. While the shape of the fossil record certainly reflects many different factors, including e.g. biogeography and external perturbations of the abiotic environment, the event statistics demonstrated here should have rather general implications for all models which do not completely dismiss the influence of population dynamics on macroevolution. If the evolutionary ‘jumps’ are the elementary events in any dynamics, possibly involving interactions between evolving species, using $`\mathrm{log}t`$ as the independent variable makes this dynamics appear as a stationary, markovian process, described e.g. by a master or Fokker-Planck equation. The eigenvalues (and eigenvectors) of the evolution equation describe the relaxation of any average of interest. Since the observational time window is usually narrow in $`\mathrm{log}`$ time, one is restricted to observing the decay of a single (or a few) relaxational mode(s). As an exponential function of $`\mathrm{log}`$ is a power-law, the above mechanism offers a generic explanation for the power-law like behavior found in several evolutionary patterns.

Aknowledgments. Both authors would like to thank Mark Newman and Richard Palmer for inspiring conversations and exchanges of ideas at the Santa Fe Institute of Complex Studies and at the Telluride Summer Research Institute. Parts of this work were commenced during a visit by A. P. at Duke University. A. P. owes a special thank to Richard Palmer for the kind hospitality extended to him and for the guidance he received during his stay. This project was partly supported by a block grant from Statens Naturvidenskabelige Forskningsråd.