text:
# Study of color connections in 𝑒⁺𝑒⁻ annihilation

## I Introduction

Strong interaction processes in high-energy reactions, e.g. $`e^+e^{}`$hadrons, are generally described by two distinct sectors: the perturbative and hadronization phases. The perturbative phase is well described by perturbative quantum chromodynamics (PQCD), while the hadronization one is non-perturbative and currently can only be described by hadronization models Lund Webber . A natural problem arises as how to link these two phases, which is beyond the capability of PQCD. In present hadronic event generators Jetset Herwig , the color connections are assumed to be the color string or the color cluster chain. This is true only in the case of $`N_C\mathrm{}`$, because with infinitely many colors the probability that two or more partons have accidentally the same color is zero, and then the way to connect the partons by a string or a cluster chain becomes unique C.Friberg1997 Q.Wang2000 . The present hadronization models work quite well, which shows that the large $`N_C`$ limit does reflect some features of the real world. However, with only three colors in nature, the color structure of a multiparton state at the end of the parton cascade is copious and complex Q.Wang2000  \- maltoni . The one that is used in the current string and cluster models is called the color singlet chain (CC) state, where in each string piece or cluster the color charge from one parton is connected to its anti-color from the other, so that the string piece or the cluster is color neutral. This results in a color neutral flow in the string or the cluster system Q.Wang2001 maltoni . PQCD calculation shows that, when projected onto the color space of the final partonic state, the CC state occurs with a probability less than 1 and the probability decreases as the growing number of partons Q.Wang2001 Y.Jin . Part of the rest probability goes to such states as the so called color separate (CS) singlet, where a multiparton system separates itself into several color singlet subsystems and each subsystem can hadronize independently C.Friberg1997 Q.Wang2000  T.Sjostrand1994 L.Lonnblad1996 .

The presence of CS states can be easily understood in $`e^+e^{}q\overline{q}+ngh^{}s`$ process. If two gluons have opposite colors, they can form a color singlet subsystem or a closed string. When the parton number is large, there are many possibilities for two or more gluons to form CS states Q.Wang2000 Q.Wang1996 . As is shown in Ref. Q.Wang2000 , CC and CS states are not orthogonal to each other in perturbative sense. They belong to two different complete sets of color singlets in color space of the final partonic state and are equivalent in the context of PQCD. However, they lead to different color connections on the hard-soft interface and may give rise to different hadronic states through the subsequent hadronization process Q.Wang2000 Q.Wang2001 . Whether the hadronization chooses CC or CS states as its starting point is determined by non-perturbative QCD dynamics. This question become critical in ultra-relativistic heavy ion collisions, since it is almost impossible to identify a unique color singlet chain in a bulk of partons.

The most efficient and practical way to study effects of CS connections on hadronic events is to modify the available event generators by substituting the CC connection with the CS one. Then we can use the new program to generate hadronic events and to investigate if there are any deviations from those generated with the default (CC) connection. An explicit example is given in Ref. Q.Wang2000 , where a phenomenological CS model is proposed by incorporating CS connection into JETSET. The results show that there are no significant differences between CS and CC connections in terms of global properties of unbiased events Q.Wang2000 . This indicates that the available data for these properties does not rule out CS possibility. In other words, color connections beyond the traditional CC one, such as the CS connection discussed above, are also possible states on the interface between the hard and soft sectors in hadronic processes. To further confirm such an observation, it is natural to seek new observables in special events which may be sensitive to color connections. In this paper, we first give other properties of unbiased events except those in Q.Wang2000  and show these two different color connections can both describe data. Then we study some properties for three-jet events and compare them with the recent OPAL OPAL2002  and DELPHI data DELPHI2000 . We find no evidence that CC is the unique color connection on the hard-soft interface. However, our investigation shows that the CC and CS states can lead to larger differences for some selected three-jet events where jets are well separated and the jet resolution scale, $`y_{cut}`$, is smaller. This could be put to further experimental test, e.g., by re-analysis of LEP data or in future Giga-Z experiment at linear collider.

The outline of this paper is as follows. In section 2 we give a brief description of the CS model Q.Wang2000 . We present in section 3 some results about unbiased events. In section 4 we study three-jet events and propose observables sensitive to color connections. A summary of results is given in section 5.

## II The color separate state model

Our PQCD analysis shows that up to $`O(1/N_C^2)`$ a partonic system with $`q\overline{q}`$ and $`n`$ gluons can be decomposed into two sub-singlets, one is the CC state made of $`q\overline{q}`$ and $`m`$ gluons, the other is the CS state made of $`(nm)`$ gluons Q.Wang2000 . The former forms an open string stretched between the quark and the anti-quark spanning $`m`$ gluons, while the latter forms a closed string by the rest gluons. These CS states formed in the partonic system to $`O(1/N_C^2)`$ are called leading CS states. Based on this analysis a phenomenological model is constructed in the following way: (1) Only leading CS states are produced; (2) A CS state is produced with the relative weight given by the T-measure B.Andersson1986 . We implement the model into a Monte-Carlo program based on JETSET7.4. In the program we just replace the default way of color connection at the end of parton cascade in JETSET7.4, i.e. the CC connection, with the CS-allowed connection without touching all other parts. The CS cluster made of two or more gluons in a CS state hadronizes as closed strings Lund .

Using this program, we have generated CS unbiased events and studied such global features as thrust, sphericity, oblateness and aplanarity. They are compared with those from CC unbiased events. The comparison shows that there are no significant differences between CS and CC unbiased events in terms of these global observables. In fact such a result is not surprising, because for unbiased events the global properties are mainly determined by PQCD parton cascade process and are not sensitive to hadronization details. This is nothing but the well-known property of local parton-hadron duality B.Andersson1989 . We therefore need to make further investigation in other observables more sensitive to color connections.

## III Investigation of unbiased events

### III.1 Identified hadron multiplicities and their momentum spectra

Two basic observables for hadronic events are multiplicities of identified hadrons and their momentum distributions. Using 6 million Monte Carlo events we calculate the average multiplicities and the inclusive momentum spectra for $`\pi ^\pm `$,$`K^\pm `$,$`p(\overline{p})`$ in unbiased events. The statistical uncertainties for these observables are much less than those in data (about 1 - 2 order of magnitude lower than data). Our results for multiplicities together with data Particle2000  at $`\sqrt{s}=10GeV,29GeV`$ and $`91GeV`$ are listed in Table 1-3. The inclusive momentum spectra for $`\pi ^\pm `$, $`K^\pm `$, $`p(\overline{p})`$ at same energies are shown in Fig. 1a-c respectively. The data are taken from Ref. Particle2000 . The results show that CS unbiased events agree with data at the same level as CC ones.

### III.2 Ratio of baryon to meson and baryon-antibaryon correlation

It is well known that in $`e^+e^{}`$ annihilation most of the measured baryons are directly produced. Even if they are decay products of primary baryons, they keep more of their rapidities and momenta from their parents than most mesons do. Therefore, the spectra of the baryons in particular the correlation between baryon (B) and anti-baryon ($`\overline{B}`$) could be a criteria to test different hadronization mechanism ZG1997 . Here we investigate if this correlation can distinguish different color connections before hadronization.

In Lund string fragmentation model, baryon production is described by diquark or popcorn mechanism. The diquark production is controlled by the parameter $`qq/q`$, the ratio of the production rate of $`qq`$ or $`\overline{q}\overline{q}`$ to that of $`q\overline{q}`$ in color field Lund . The parameter $`qq/q`$ controls the ratio of baryon to meson. Six million Monte Carlo events are produced and the statistical uncertainty is much reduced as compared to data. The calculated ratios of baryon to meson at $`\sqrt{s}=10GeV`$, 29GeV and 91GeV are listed in Table 4. We see that there is no significant difference between CS and CC connections. For the CC connection, i.e. default setting in JETSET, one has to take different values of $`qq/q`$ at different energies to fit data ($`qq/q=0.06`$ at $`\sqrt{s}=10GeV`$, and $`qq/q=0.1`$ at other two energies). On the other hand, one can use an energy-independent value $`qq/q=0.1`$ for CS events, and the results agree with data quite well.

The popcorn mechanism is realized by the popcorn parameter:

$$\rho =\frac{P(BM\overline{B})}{P(B\overline{B})+P(BM\overline{B})},$$
(1)
where $`\rho =1`$ means that a meson is always produced between B and $`\overline{B}`$, while in the other limit $`\rho =0`$ it is not at all. The popcorn parameter $`\rho `$ also describes the $`B\overline{B}`$ correlation. Here we take $`\rho =0.5`$ and calculated the $`p\overline{p}`$, $`\mathrm{\Lambda }\overline{\mathrm{\Lambda }}`$ and $`p\overline{\mathrm{\Lambda }}/\mathrm{\Lambda }\overline{p}`$ rapidity correlations for unbiased events at $`Z^0`$ pole. The results are shown in Fig. 2a-c. Both CS and CC connections agree with data DELPHI1997  at the same level.

Besides the global properties studied in Q.Wang2000 , we have investigated other properties for unbiased events with CS or CC color connections. The results obtained in both cases are consistent with data. For some results the agreement is even better in CS events.

## IV Investigation of three-jet events at $`Z^0`$ pole

Since different color connections lead to different color strings and may cover different phase space, one expects that the CS state should have some effects on the energy and momentum flow of final hadrons. So the impact of the CS connection could be observable in some specific events. Currently a large number of three-jet events, where a gluon-jet can be identified, are available at $`Z^0`$ resonance energy at LEP I. It is highly possible that the gluon-jet could be the fragmentation product of gluonic color singlet clusters which are CS singlets. Therefore, we expect that some sensitive observables should be out there in a specified window of phase space around the gluon jet. In the following, we will study the properties of three-jet events at $`Z^0`$ resonance energy to probe the possible effect of CS states. We note that a similar work has been done for two-jet events SY2002 .

### IV.1 Selection of three-jet events

We use Durham jet algorithm Catani1991  to select three-jet events. In the algorithm, particles are grouped into jets based on a cutoff variable $`y_{cut}`$ defined as follows. For each pair of particles $`i`$ and $`j`$ with energies $`E_i`$ and $`E_j`$ respectively and their opening angle $`\theta _{ij}`$, one defines $`y_{ij}`$ as

$$y_{ij}\frac{2\mathrm{m}\mathrm{i}\mathrm{n}(E_i^2,E_j^2)(1\mathrm{cos}\theta _{ij})}{E_{vis}^2},$$
(2)
where $`E_{vis}`$ is the total visible energy in the event. Given a cutoff value $`y_{cut}`$, one can put all particles into one jet where all $`y_{ij}`$’s for any two particles inside the jet are smaller than $`y_{cut}`$, while $`y_{ij}>y_{cut}`$ for the case with one inside the jet and the other outside it. The number of particles in a jet is then the multiplicity of the jet. Increasing $`y_{cut}`$ makes the fraction of multi-jet events lower because $`y_{cut}`$ is actually a jet resolution parameter.

We have calculated the fraction of three-jet events as function of $`y_{cut}`$. Both color connections agree with dataP.D.Acton1992 . We have also calculated the charged particle multiplicity for three-jet events at some values of $`y_{cut}`$ and compared our results with data DELPHI1992 . There is no sizable difference between CS and CC events.

### IV.2 Particle multiplicities and momentum spectra of identified hadrons in Y events

Recently OPAL collaboration has measured the charged particle multiplicity $`N_{ch}`$ for three-jet light-quark events of the so-called Y shape from $`Z^0`$ decay as function of the angle $`\theta `$ between the two jets with lowest energies OPAL2002 . The jets are defined by Durham jet-finder. The resolution $`y_{cut}`$ is adjusted separately for each tagged $`uds`$ event so that three jets can be exactly reconstructed. They applied this procedure instead of using a fixed resolution value just to avoid introducing a bias into the gluon jet. The jets are labelled 1,2 and 3 in such an order that jet 1 carries the highest energy. Only those events, the so-called Y events, are kept where the angles between the most energetic jet and the rest two are the same up to $`3^{}`$. We are particularly interested in these data which could be used to test our CS model. We choose such Y events from 10 million Monte Carlo samples at $`Z^0`$ pole accommodating the same experimental condition. The calculated particle multiplicities varied with angle $`\theta `$ are shown in Fig. 3. One can see that there is some difference between two types of color connections and that results from the CC connection seem to agree with data better than those from the CS one. However, such a difference is not significant enough to discriminate the two in terms of data. The multiplicity distribution calculated for these Y events with $`0^{}<\theta <120^{}`$ is shown in Fig. 4. Besides these results, we also calculated other observables measured in Ref. OPAL2002 . All these results fail to show any significant difference between the two types of connections.

DELPHI collaboration has measured momentum spectra of identified hadrons in quark jets of all species ($`duscb`$) together with gluon jets for Y eventsDELPHI2000 . The three-jet events are identified by Durham algorithm with $`y_{cut}=0.015`$. The jets are also labelled in descending order of jet energies. The Y events are obtained with $`\theta _2,\theta _3[150^{}15^{},150^{}+15^{}]`$, where $`\theta _2`$ and $`\theta _3`$ are angles between the highest energetic jet and other two jets with lower energies. In order to compare with data of this experiment, we calculate the momentum spectra of identified hadrons in quark and gluon jets for Y events with the same constraints as in the experiment. The samples are 10 million Monte Carlo events at $`Z^0`$ pole. The results are shown in Fig. 5. One sees that it is hard to distinguish the CS from the CC connection by comparison with data. So up to now, no evidence has ruled out the CS connection as a possible candidate state on the hard-soft interface.

However, the above results imply that some differences between the CC and CS connections are present. In the next section we will investigate whether they could be amplified by choosing suitable $`y_{cut}`$ and the angles between jets.

### IV.3 Observables in special three-jet events

One expects that more sensitive observables come from the phase space around the gluon jet. A normal method to identify the gluon-jet in experiments is to label three jets by their energies in descending order $`E_1E_2E_3`$. The most energetic jet (jet 1) is usually recognised as the quark-jet (or antiquark-jet), and the one with the smallest energy (jet 3) as the gluon jet. We denote the angle between the first and the second jet by $`\theta _{12}`$ and that between the first and the third one by $`\theta _{13}`$.

We look at three-jet events with restricted range of $`\theta _{12}`$ and $`\theta _{13}`$. We calculate charged particle multiplicities whose sensitivity to different color connections varies with $`\theta _{12}`$ and $`\theta _{13}`$. We find that only those events where jet 3 is well separated from jet 1 and 2 exhibit considerable difference in multiplicity from color connections. Considering both the statistics and the sensitivity to color connection, we choose $`100^{}<\theta _{12},\theta _{13}<160^{}`$. We also find that the difference increases with decreasing $`y_{cut}`$. Here we use a fixed value of the resolution parameter $`y_{cut}`$ instead of one for each three-jet event.

Now we focus on the gluon jet in three-jet events. We calculate the multiplicity, the invariant mass and the longitudinal momentum of the gluon jet as functions of $`y_{cut}`$. Compared with those from CC connections, the results from CS ones are found to be smaller for the multiplicity and the invariant mass, and larger for the longitudinal momentum. We also find that the difference between CC and CS is more significant for smaller $`y_{cut}`$, as shown in Fig.6. We select three-jet events from 10 million Monte Carlo samples at $`Z^0`$ pole. Our numerical results are given with constraint $`100^{}<\theta _{12},\theta _{13}<160^{}`$ at $`y_{cut}=0.0005`$. The fraction of selected three-jet events in the total ones is about $`0.502\%`$. We give in Fig. 7 the multiplicity distribution. The result shows that two different ways of color connections do bring larger differences in the final hadronic state. The shapes of two multiplicity distributions look similar, but the peak of the distribution for the CS case locates at a smaller value of $`n_{ch}`$ than that for the CC case, implying a lower average multiplicity: $`n_{ch}_{CS}=18.22`$, $`n_{ch}_{CC}=19.41`$, and the difference is $`\mathrm{}n_{ch}=1.19`$. The reason is obvious: in CS events several groups of gluons form color singlet clusters and each hadronizes independently, which makes the effective energy to produce hadrons smaller than that in CC events.

Considering that the average longitudinal momentum of the gluon-jet for CS events is larger than that for CC ones, we define the following observable:

$$T=\frac{\underset{i}{}p_{zi}}{\underset{i}{}|𝐩_i|},$$
(3)
where $`𝐩_i`$ and $`p_{zi}`$ are the 3-momentum and its longitudinal component of a final state hadron in the gluon-jet respectively. Here longitudinal means parallel to the gluon-jet axis. The observable $`T`$ is introduced to amplify the difference of CC and CS events. The distribution of observable $`T`$ with both CS and CC connections is given in Fig. (8a) and shows that T is a sensitive observable to color connections. We also calculated the spectra of the charged particle multiplicity, the rapidity and the invariant mass in the gluon-jet, as shown in Fig.(8b-d). One sees that now the predicted distinction between CC and CS connections is obvious. The reason is that the formation of closed strings in a gluon jet with the CS connection takes more momentum along the the direction of the mother parton, i.e. the gluon, than the CC one, which makes the gluon-jet thinner. This point can be further strengthened by the spectra of the polar angle and the longitudinal momentum for hadrons in the gluon jet. We note that all results in Fig.8 are consistent with each other.

## V Summary

In this paper, the comparison of our results with the available data up to now shows that the CS state is a possible choice for the hard-soft interface in hadronic processes. This is consistent with PQCD analysis. Such a picture of colour connection is natural for a deconfined quark-gluon system, where a huge number of partons interact with each other and it is impossible to identify a unique color singlet chain. This paper is also an attempt for searching observables in special events where the differences between the CS and CC states are significant. The identification of these differences in experiments will deepen our understanding of the hard-soft interface in hadronic processes.

By a careful look at certain type of three-jet events with smaller $`y_{cut}`$, where angles between jets are restricted in a specific range, we have found more sensitive observables to color connections. The difference between CC and CS connections is amplified in the charged particle multiplicity and properties of the gluon-jet in these events. The distinction between CS and CC events can be put to test by re-analysing the LET-I data and/or in future experiments such as Giga-Z at linear collider.

### Acknowledgements

The authors thank Z.-T. Liang and Z.-G. Si for insightful discussions. This work is supported in part by the National Natural Science Foundation of China under grant No. 10075031 and 10205009.


format: text

text:
# Punctuated equilibria and 1/𝑓 noise in a biological coevolution model with individual-based dynamics

## I Introduction

Biological evolution offers a number of important, unsolved problems that are well suited for investigation by methods from statistical physics. Many of these can be studied using complex, interacting model systems far from equilibrium DROS01 . Areas that have generated exceptional interest among physicists are those of coevolution and speciation. A large class of coevolution models have been inspired by one introduced by Bak and Sneppen BAK93 , in which species with different levels of “fitness” compete, and the least fit species and those that interact with it are regularly “mutated” and replaced by new species with different, randomly chosen fitness. Models in this class exhibit avalanches of extinctions and appear to evolve towards a self-organized critical state BAK87 ; BAK88 . Although such models may be said to incorporate Darwin’s principle of “survival of the fittest,” they are artificial in the sense that mutation and selection are assumed to act collectively on entire species, rather than on individual members of their populations.

In reality, mutations are changes in the genotypes of individual organisms that are introduced or passed on during reproduction. These changes in the genotype affect the phenotype (physical and behavioral characteristics of the organism and its interactions with other organisms), and it is on the level of the phenotypes of individuals that competition and selection act. A number of evolution models (see, e.g., Ref. DROS01  for a review from a physicist’s point of view) therefore take as their basis a genome in the form of a string of “letters” as in Eigen’s quasi-species model of molecular evolution EIGE71 ; EIGE88 . Depending on the level of the modeling, the letters of the genomic “alphabet” can be a (possibly large) number of alleles at a gene locus, or the four nucleotides of a DNA or RNA sequence BAAK99 . However, the size of the alphabet is not of great importance in principle, and it is common in models to use a binary alphabet with the two letters 0 and 1 (or $`\pm 1`$) DROS01 ; EIGE71 ; EIGE88 ; BAAK99 .

We believe a fruitful approach to the study of coevolution is one in which selection is provided by interspecies interactions along lines commonly considered in community ecology, while genetic variation is provided by random mutations in the genomes of individual organisms. An early attempt in this direction is the coupled $`NK`$ model with population dynamics introduced by Kauffman and Johnsen KAUF91 ; KAUF93 , but thus far not many similar models have been studied. Recently, Hall, Christensen, and collaborators HALL02 ; CHRI02 ; COLL03  introduced a model they called the “tangled-nature” model, in which each individual lives in a dynamically evolving “fitness landscape” created by the populations of all the other species. Here we consider a simplified version of this model, in which no individual is allowed to live through more than one reproduction cycle. This restriction to nonoverlapping generations enables us to both study the model in detail by linear stability analysis and to perform very long Monte Carlo simulations of its evolutionary behavior. In short, the model consists of populations of different species, on which selection acts through asexual reproduction rates that depend on the populations of all the other species via a constant, random interaction matrix that allows mutualistic, competitive, and predator/prey relations. In addition to these direct interactions, all individuals interact indirectly through competition for a shared resource. The competition keeps the total population from diverging. Genetic variety is provided by a low mutation rate that acts on the genomes of individual organisms during reproduction, inducing the populations to move through genotype space. The resulting motion is intermittent: long periods of stasis with only minor fluctuations are interrupted by bursts of significant change that is rapid on a macroscopic time scale, reminiscent of what is known in evolutionary biology as punctuated equilibria GOUL77 ; GOUL93 ; NEWM85 . A preliminary report on some aspects of the work reported here is given in Ref. RIKV03A .

The main foci of the present paper are the structure and stability of communities and the statistical properties of the dynamical behavior on very long (“geological”) time scales. The rest of the paper is organized as follows. In Sec. II we describe our model in detail, including the detailed Monte Carlo algorithm used for our simulations. In Sec. III we discuss the properties of the fixed points of the population for the mutation-free version of the model. Many of these properties can be understood analytically within a simple mean-field approach. A full, probabilistic description of our model is also possible, although the mathematics is somewhat involved. In Appendix A we provide, for the sake of completeness, the key equations of this approach. In Sec. IV we give a detailed report on our large-scale Monte Carlo simulations and the numerical results, together with a discussion of their relations to the fossil record and to other theoretical models. Finally, in Sec. V, we summarize our results and give our conclusions and suggestions for future studies.

## II Model and algorithm

As mentioned in Sec. I, we use a simplified version of the tangled-nature model introduced by Hall, Christensen, and collaborators HALL02 ; CHRI02 ; COLL03 . It consists of a population of individuals with a genome of $`L`$ genes, each of which can take one of the two values 0 or 1. Thus, the total number of different genotypes is $`𝒩_{\mathrm{max}}=2^L`$. We consider each different genotype as a separate species, and we shall in this paper use the two terms interchangeably. This is justified by the idea that in the relatively short genomes we can consider computationally, each binary “gene” actually represents a group of genes in a coarse-grained sense.

In our version of the model, the population evolves in discrete, nonoverlapping generations (as in, e.g., many insects), and the number of individuals of genotype $`I`$ in generation $`t`$ is $`n_I(t)`$. The total population is $`N_{\mathrm{tot}}(t)=_In_I(t)`$. In each generation, the probability that an individual of genotype $`I`$ produces a litter of $`F`$ offspring before it dies is $`P_I(\{n_J(t)\})`$, while the probability that it dies without offspring is $`1P_I`$. Although the fecundity $`F`$ could be quite complex in reality (e.g., a function of $`I`$ on the average, but random both in time and for each individual), we here take a simplistic approach and assume it is a constant, independent of $`I`$ and $`t`$. The main difference from the model of Refs. HALL02 ; CHRI02 ; COLL03  is that the generations are nonoverlapping in our model, while individuals in the original model may live through several successive reproduction cycles. This simplification facilitates theoretical analysis and enables us to make significantly longer simulations than those reported in Refs. HALL02 ; CHRI02 ; COLL03 . In the “opposite” direction, the earlier model could be generalized to include non-trivial “age structures,” since individuals living through several cycles can be assigned an age, with age-dependent survival and reproductive properties LESL45 ; PIEL69 ; CHAR94 ; PENN95 ; HOWA01 ; ZIA02 ; TUZE01A ; CHOW03A ; CHOW03B .

As in the original model HALL02 ; CHRI02 ; COLL03 , the reproduction probability $`P_I`$ is taken as

$$P_I(\{n_J(t)\})=\frac{1}{1+\mathrm{exp}\left[_JM_{IJ}n_J(t)/N_{\mathrm{tot}}(t)+N_{\mathrm{tot}}(t)/N_0\right]}.$$
(1)
Here, the Verhulst factor $`N_0`$ VERH1838  represents an environmental “carrying capacity” that might be due to limitations on shared resources, such as space, light, or water. It prevents the total population from indefinite growth, stabilizing it at $`O\left(N_0\right)`$. The interaction matrix $`𝐌`$ expresses pair interactions between different species such that the element $`M_{IJ}`$ gives the effect of the population density of species $`J`$ on species $`I`$. Thus, a mutualistic relationship is represented by both $`M_{IJ}`$ and $`M_{JI}`$ being positive, while both being negative models a competitive relationship. If they are of opposite signs, we have a predator-prey situation. To concentrate attention on the effects of interspecies interactions, we follow Refs. HALL02 ; CHRI02 ; COLL03  in setting the self-interactions $`M_{II}=0`$. The offdiagonal elements of $`M_{IJ}`$ are randomly and uniformly distributed on $`[1,1]`$ as in Ref. CHRI02 . The interaction matrix is set up at the beginning of the simulation and is not changed later, a situation that corresponds to quenched disorder in spin-glass models FISC91 . We note that we have not attempted to give $`𝐌`$ a particularly biologically realistic form. Some possible modifications are discussed in Sec. V.

In this model there is a one-to-one correspondence between genotype (the $`I`$th specific bit string) and phenotype (the $`I`$th row and column of $`𝐌`$). Thus, the phenotype specifies both how the $`I`$th species influences the other species that are either actually or only potentially present in the community (the $`I`$th column) and how it is influenced by others (the $`I`$th row). The reproduction probability $`P_I`$ provides selection of the “most fit” phenotypes according to these “traits” (matrix elements) and the populations of the other species present in the community. The effect (or lack thereof) of a particular trait depends on the community in which the species exists, just as a cheetah’s superior speed is only relevant to survival if fast-moving prey is available.

In each generation, the genomes of the individual offspring organisms are subjected to mutation with probability $`\mu /L`$ per gene and individual. Mutated offspring are re-assigned to their new genotypes before the start of the next generation. This provides the genetic variability necessary for evolution to proceed.

The analytic form of $`P_I(\{n_J(t)\})`$, Eq. (1), is by no means unique. For instance, one could use a “soft dynamic” RIKV02  in which the effects of the interactions and the Verhulst factor factorize in $`P_I`$. Alternatively, one could dispense with the exponential and use linear or bilinear relationships instead, as is most common in population biology MURR89 . Investigations of the effects of such modifications are left for the future.

Our evolution algorithm proceeds in three layers of nested loops.

* Loop over generations $`t`$.
* Loop over the $`𝒩(t)`$ populated genotypes $`I`$ in generation $`t`$.
* Loop of length $`n_I(t)`$ over the individuals of genotype $`I`$. Each individual produces $`F`$ offspring for generation $`t+1`$ with probability $`P_I(\{n_J(t)\})`$, or dies without offspring with probability $`1P_I`$. In either case, no individual survives from generation $`t`$ to generation $`t+1`$.
* Loop of length equal to the total number of offspring of genotype $`I`$ generated in loop 3a, attempting to mutate each gene of each individual offspring with probability $`\mu /L`$ and moving mutated offspring to their respective new genotypes for generation $`t+1`$.

## III Linear stability analysis

Though neither our model nor the simulations are deterministic, a number of their gross properties can be understood in terms of a mean-field approximation that ignores statistical fluctuations and correlations. The time evolution of the populations is then given by the set of difference equations,

$$n_I(t+1)=n_I(t)FP_I(\{n_J(t)\})[1\mu ]+(\mu /L)F\underset{K(I)}{}n_{K(I)}(t)P_{K(I)}(\{n_J(t)\})+O(\mu ^2),$$
(2)
where $`_{K(I)}`$ runs over the $`L`$ genotypes $`K(I)`$ that differ from $`I`$ by one single mutation (i.e., the Hamming distance HAMM86  $`H_{KI}=1`$). The corrections of $`O(\mu ^2)`$ correspond to multiple mutations in single individuals. Naturally, a full investigation of this equation is highly non-trivial, even in the absence of random noise. In particular, Eq. (2) can be regarded as a logistic map in a $`2^L`$-dimensional space. Logistic maps are known to admit, in general, both fixed points and cycles of non-trivial periods MURR89 . To keep both analysis of simulations and theoretical investigations manageable, we choose parameters in such a way that we can focus our attention on fixed points.

Of course, what we simulate are stochastic processes. However, as long as the mutation rate is well below the error threshold for mutational melt-down DROS01 ; EIGE71 ; EIGE88 ; BAAK99 ; HALL02 ; COLL03 , new successful mutations can become fixed in the population before another successful mutant arises. This results in a separation of time scales between the ecological scale of fluctuations within fixed-point communities and the evolutional scale of the durations of such communities DOEB00 , which are known as quasi-steady states (QSS) or quasi evolutionarily stable strategies (q-ESS) HALL02 ; CHRI02 . By ignoring mutations, we can make some analytical progress and gain some insight into the nature and stability properties of these QSS.

### III.1 $`𝒩`$-species fixed points

A fixed point of Eq. (2), defined by

$$n_I^{}(t+1)=n_I^{}(t)n_I^{},$$
(3)
is characterized by having only $`𝒩`$ ($`2^L`$) non-vanishing $`n_I^{}`$. (Henceforth, an asterisk will signify a quantity as a fixed-point value.) Following Ref. ROBE74 , we denote a fixed point as feasible if $`n_I^{}>0`$ for all $`𝒩`$ values of $`I`$. Corresponding to the coexistence of $`𝒩`$ genotypes, such a point will be simply referred to as an “$`𝒩`$-species fixed point.” When mutations are ignored, Eq. (2) reduces to

$$n_I(t+1)=n_I(t)FP_I(\{n_J(t)\}).$$
(4)
Specializing to the form of $`P_I`$ given by Eq. (1), we see that all single-species fixed points are trivially “identical,” with $`n_I^{}=N_{\mathrm{tot}}^{}=N_0\mathrm{ln}(F1)`$. This somewhat unrealistic result is just a consequence of our choice of $`M_{II}=0`$ and can be avoided by lifting this restriction. However, the absence of self-interactions places no restriction on the main purpose of our work – the exploration of the effects of random but time independent interspecies interactions.

Proceeding to the $`𝒩2`$ cases, the existence of an $`𝒩`$-species fixed point depends critically on the submatrix $`\stackrel{~}{𝐌}`$, with matrix elements $`M_{IJ}`$ in which both $`I`$ and $`J`$ are among the $`𝒩`$ species in question. (We shall use the tilde to emphasize that a quantity corresponds to an $`𝒩`$-species subspace, rather than to the full, $`2^L`$-dimensional genotype space.) In particular, if $`\stackrel{~}{𝐌}`$ is non-singular, then a unique fixed point exists, as we show below. On the other hand, singular $`\stackrel{~}{𝐌}`$’s may result in a variety of “degenerate” cases. We provide just two examples to illustrate the mathematical richness of singular $`\stackrel{~}{𝐌}`$’s. If all elements vanish, then the behavior in this subspace is highly degenerate, with the *total* population given again by $`N_{\mathrm{tot}}^{}=N_0\mathrm{ln}(F1)`$, regardless of $`𝒩`$. However, the fractions of each species,

$$\rho _I^{}n_I^{}/N_{\mathrm{tot}}^{},$$
(5)
are completely undetermined, an understandable consequence of having dynamically indistinguishable species. Another example of a singular interaction submatrix is $`\stackrel{~}{𝐌}=\left(\begin{array}{cc}0& m>0\\ 0& 0\end{array}\right)`$ in an $`𝒩=2`$ subspace. Then, Eq. (4) will drive one of the species to extinction, so that no fixed point with both $`n_I>0`$ (stable or unstable) can exist in this two-dimensional subspace, and the system collapses to $`𝒩=1`$. For the remainder of this article, we shall study analytically only the dynamical behavior of interacting species with nonsingular $`\stackrel{~}{𝐌}`$’s.

To proceed, we insert Eq. (3) into Eq. (4) and arrive at the $`𝒩`$ equations

$$FP_I(\{n_J^{}\})=1.$$
(6)
With our choice of $`P_I`$, these lead to

$$\underset{J}{}M_{IJ}\rho _J^{}=\frac{N_{\mathrm{tot}}^{}}{N_0}\mathrm{ln}(F1).$$
(7)
(The analysis in this section remains valid even if $`M_{II}0`$, although $`M_{II}=0`$ is the case explicitly considered elsewhere in this paper.) Note that the right-hand side of Eq. (7) is *independent* of $`I`$. Since $`\stackrel{~}{𝐌}`$ is non-singular, we define its inverse by $`\stackrel{~}{𝐖}`$, with elements $`\stackrel{~}{W}_{IJ}`$:

$$\stackrel{~}{W}_{IJ}(\stackrel{~}{𝐌}^1)_{IJ}.$$
(8)
Then Eq. (7) can be inverted to give

$$\rho _I^{}=\left[\frac{N_{\mathrm{tot}}^{}}{N_0}\mathrm{ln}(F1)\right]\underset{J}{}\stackrel{~}{W}_{IJ}.$$
(9)
The only unknown in this equation, $`N_{\mathrm{tot}}^{}`$, can now be found via the normalization condition, $`\rho _I^{}=1`$, which leads to

$$N_{\mathrm{tot}}^{}/N_0=\mathrm{ln}(F1)+1/\stackrel{~}{\mathrm{\Sigma }},$$
(10)
where

$$\stackrel{~}{\mathrm{\Sigma }}\underset{IJ}{}\stackrel{~}{W}_{IJ}.$$
(11)
Putting these into Eq. (9), we have the explicit form of the fixed-point populations:

$$n_I^{}=N_0\left[\mathrm{ln}(F1)+\frac{1}{\stackrel{~}{\mathrm{\Sigma }}}\right]\frac{_J\stackrel{~}{W}_{IJ}}{\stackrel{~}{\mathrm{\Sigma }}}.$$
(12)
Although Eq. (12) appears to provide fixed-point values for *any* choice of control parameters, we emphasize that it is applicable only for a limited range of $`F`$ and $`\stackrel{~}{𝐌}`$. The subtlety lies in its stability properties within the $`𝒩`$-species subspace. First, we remind the reader that, even in the case of $`𝒩=1`$, there are a variety of behaviors. Solutions may have a stable fixed point with either monotonic or oscillatory decay of small deviations, or they may show bifurcations, period doubling, and chaos MURR89 . In the present study, we are interested in the effects of the interspecies interactions, and we therefore choose to focus on systems with monotonically decaying fluctuations in the noninteracting limit. This means that small fluctuations about the fixed point must decay as $`\delta (t+1)/\delta (t)0,1`$. For the single-species and $`𝐌=\mathrm{𝟎}`$ cases, this restriction is easily translated to the condition $`2<F4.5`$ for the fecundity. We use $`F=4`$ in all of our simulations.

Next, through a straightforward linear stability analysis around the $`𝒩`$-species fixed point, we obtain a condition that represents decay of all small perturbations, i.e., all the eigenvalues of the matrix

$$\frac{n_I(t+1)}{n_J(t)}|_{\{n_I^{}\}}$$
(13)
must have real parts, lying between $`1`$ and 1. Carrying out the differentiation and using Eq. (4) to simplify the result, we find this matrix to be of the form $`\mathrm{𝟏}+\stackrel{~}{𝚲}`$, where $`\mathrm{𝟏}`$ is the $`𝒩`$-dimensional unit matrix and the elements of $`\stackrel{~}{𝚲}`$ are given by

$$\stackrel{~}{\mathrm{\Lambda }}_{IJ}=\left(1\frac{1}{F}\right)\rho _I^{}\left(M_{IJ}\mathrm{ln}(F1)2/\stackrel{~}{\mathrm{\Sigma }}\right).$$
(14)
Our criterion translates to the requirement that all the eigenvalues of $`\stackrel{~}{𝚲}`$ must have real parts that lie in $`2,0`$. Fixed points with this property will be called “internally stable.” Unfortunately, this criterion cannot be made more explicit. Given $`F`$ and a set of $`M_{IJ}`$, $`\stackrel{~}{𝚲}`$ must be constructed using Eqs. (8)–(11) and diagonalized. The matrix $`\stackrel{~}{𝚲}`$ is recognized as what is known in ecology as the community matrix MURR89  of the $`𝒩`$-species fixed point for the discrete-time dynamic defined by Eq. (4).

It has been shown, both numerically GARD70  and analytically MAY72 , that the proportion of large, random matrices for which all eigenvalues have a negative real part vanishes as the proportion of nonzero matrix elements increases. This has been used as an argument that highly connected ecosystems are intrinsically unstable MAY72 , contrary to ecological intuition. However, the matrix that must be studied to determine the internal stability of a fixed point is the community matrix $`\stackrel{~}{𝚲}`$, which has a complicated relationship to the interaction matrix $`𝐌`$ and should not be expected to have a simple element distribution. In fact, for some bilinear population dynamics models there is numerical evidence that most feasible fixed points are also internally stable ROBE74 ; GILP75 . However, the relations between connectivity and stability have not yet been fully clarified and are still being discussed ROZD01 ; WILM02 .

Of course, issues of internal stability are somewhat academic for simulations. In practice an internally unstable fixed point could not be observed for more than a brief time, especially since the populations would be driven away from such fixed points by the noise due to both the birth/death process and the mutations.

### III.2 Stability against other species

Since mutations are essential for the long-time evolutionary behavior, the “external stability” properties of the fixed-point communities are important. Even if the population corresponds to an internally stable $`𝒩`$-species fixed point, mutations will generate small populations of “invader” species, i.e., genotypes outside the resident $`𝒩`$-species community. Denoting such invaders by the subscript $`i`$ and linearizing Eq. (4) about the $`𝒩`$-species fixed point, we see that the important quantity is the multiplication rate of the invader species in the limit of vanishing $`n_i/N_{\mathrm{tot}}`$,

$$\frac{n_i(t+1)}{n_i(t)}=\frac{F}{1+(F1)\mathrm{exp}\left[\left(1_{JK}M_{iJ}\stackrel{~}{W}_{JK}\right)/\stackrel{~}{\mathrm{\Sigma }}\right]}.$$
(15)
To obtain this result we exploited Eqs. (10) and (12). Explicitly, the condition for stability against the invader, $`n_i(t+1)/n_i(t)<1`$, reduces to the requirement that the argument of the exponential function in Eq. (15) must be positive. We also note that, if $`M_{iJ}=0`$ for all $`J`$ in the resident community, then the multiplication rate of the invader equals unity. In population biology the Lyapunov exponent $`\mathrm{ln}\left[n_i(t+1)/n_i(t)\right]`$ is known as the invasion fitness of the mutant with respect to the resident community DOEB00 ; METZ92 . From Eq. (15) it becomes clear that the success of an invader depends not only on its direct interactions $`M_{iJ}`$ with each of the resident species, but also on the interactions between the resident species through the inverse interaction submatrix $`\stackrel{~}{W}_{JK}`$.

## IV Simulation results

The model described in Sec. II was studied by Monte Carlo simulations with the following parameters: $`L=13`$, $`F=4`$, $`N_0=2000`$, and $`\mu =10^3`$. The random matrix $`𝐌`$ (with zero diagonal and other elements randomly distributed on $`[1,1]`$) was chosen at the beginning of each simulation run and then kept constant for the duration of the run. In this regime both the number of populated species $`𝒩(t)`$ and the total population $`N_{\mathrm{tot}}(t)N_0\mathrm{ln}(F1)2200`$ are substantially smaller than the number of possible species, $`𝒩_{\mathrm{max}}=2^L=8192`$. This appears biologically reasonable in view of the enormous number of different possible genotypes in nature. (But see further discussion in Sec.IV.5 and Appendix B.) In a QSS, the average number of mutant offspring of species $`I`$ in generation $`t+1`$ is approximately $`\mu n_I(t)`$. Thus, with the mutation rate used here, each of the dominant species will produce of the order of one mutant organism per generation. As shown in Sec. IV.2, during a QSS most of these mutants become extinct after one generation. However, due to the small genome size, the same mutant of a species with a large population will be re-generated repeatedly by mutation from the parent species.

A high-population genotype with its “cloud” of closely related low-population mutants could be considered a quasi-species in the sense of Eigen, with the high-population genotype as the “wildtype” DROS01 ; COLL03 ; EIGE71 ; EIGE88 ; BAAK99 . An alternative interpretation of the model is therefore as one of the coevolution of quasi-species COLL03 , in which the successful invasion of a resident community by a mutant represents a speciation event in the lineage of the parent genotype.

### IV.1 General features

Most of our simulations were started with a small population of 100 individuals of a randomly selected single genotype, corresponding to the entry into an empty ecological niche by a small group of identical individuals. However, runs starting from a random number of small populations of different species give essentially the same results. Generally, the initial species are not likely to be stable against mutants, and they usually become extinct within 100 generations. To eliminate any short initial transients, most of the quantitative analyses presented here are based on time series from which the first 4096 generations were removed. The simulated quantities were recorded every 16 generations. This time resolution was chosen to be just larger than the average time it would take for the descendants of a single individual of a mutant genotype $`i`$ to completely replace a resident genotype $`J`$ of $`N_0=2000`$ individuals in the case of a maximally aggressive mutant, $`M_{iJ}=+1`$ and $`M_{Ji}=1`$. This time, which is obtained by numerical solution of Eq. (15), is about 15 generations and represents a “minimum growth time” for the model. Sampling at shorter intervals would mostly add random noise to the results, while sampling on a much coarser scale could miss important evolutionary events. It is easy to show, by solving Eq. (15) analytically for short times, that the growth time and thus the optimal sampling interval increases logarithmically with $`N_0`$. Most quantitative results in this paper are based on averages over 16 independent simulation runs of $`2^{25}=\mathrm{33\hspace{0.17em}554\hspace{0.17em}432}`$ generations each COMP .

We define the diversity of the population as the number of species with significant populations, thus excluding small populations of mainly unsuccessful mutants of the major genotypes. Operationally we define the diversity as $`D(t)=\mathrm{exp}\left[S\left(\{n_I(t)\}\right)\right]`$, where $`S`$ is the information-theoretical entropy SHAN48 ; SHAN49 ,

$$S\left(\{n_I(t)\}\right)=\underset{\{I|\rho _I(t)>0\}}{}\rho _I(t)\mathrm{ln}\rho _I(t)$$
(16)
with

$$\rho _I(t)=n_I(t)/N_{\mathrm{tot}}(t).$$
(17)
This measure of diversity is known in ecology as the Shannon-Wiener index KREB89 . It is different from the definition of diversity as the number of populated species (known in ecology as the species richness KREB89 ) that was used in Refs. HALL02 ; CHRI02 . The entropy-based measure significantly reduces the noise during QSS, caused by unsuccessful mutations.

The Shannon-Wiener diversity index is shown in Fig. 1(a) for a simulation of $`10^6`$ generations, together with the normalized total population, $`N_{\mathrm{tot}}(t)/[N_0\mathrm{ln}(F1)]`$. We see quiet periods during which $`D(t)`$ is constant except for small fluctuations, separated by periods during which it fluctuates wildly. The total population $`N_{\mathrm{tot}}(t)`$ is enhanced relative to its noninteracting fixed-point value during the quiet periods, while it decreases toward the vicinity of this value during the active periods.

To verify that different quiet periods indeed correspond to different resident communities $`\{n_I\}`$, we show in Fig. 2 the genotype labels $`I`$ (integers between 0 and $`2^L1`$, corresponding to the decimal representation of the genotype bit string) versus time for the populated species. The populations of the different species are indicated by the grayscale (by color online). We see that, in general, the community is completely rearranged during the active periods, so that the quiet periods can be identified with the QSS of the ecology. An exception is afforded by the event near $`t=1.5\times 10^5`$ generations in Figs. 1(a) and 2(a). In this instance the populations of the dominant species decrease, and the total population is for a brief time spread over a larger number of species. However, the dominant species “regain their footing,” and the original resident community continues for approximately another 50 000 generations. This situation is reminiscent of rare events in nucleation theory RIKV94 , in which a fluctuation of the order of a critical or even a supercritical droplet nevertheless may decay back to the metastable state.

### IV.2 Stability of communities against invaders

To investigate the stability properties of individual communities against invaders, we chose from a particular simulation run of $`10^6`$ generations ten different QSS of durations longer than 20 000 generations. The genotype labels and fixed-point populations that characterize these QSS are given in the first four columns of Table 1. The population-weighted average of the Hamming distances $`H_{IJ}`$ between pairs of genotypes in the community,

$$H=\underset{I}{}\underset{J>I}{}\rho _I\rho _J\left(\frac{1}{1\rho _I}+\frac{1}{1\rho _J}\right)H_{IJ},$$
(18)
and the corresponding standard deviation,

$$\sigma _H=\sqrt{\frac{𝒩}{𝒩1}\underset{I}{}\underset{J>I}{}\rho _I\rho _J\left(\frac{1}{1\rho _I}+\frac{1}{1\rho _J}\right)\left(H_{IJ}H\right)^2},$$
(19)
are shown in the fifth and sixth column, respectively. Even though the community as a whole moves far through the $`2^L`$-dimensional population space, the different QSS communities are seen to retain the property that they consist of relatively close relatives (which of course are all descendants of the single, initial genotype). In Sec. IV.5 we demonstrate that $`H`$ and $`\sigma _H`$ remain in the range shown in Table 1, even during very long simulations.

The averages and standard deviations of the offdiagonal interaction-matrix elements $`M_{IJ}`$ between the community members are shown in the seventh and eighth column of Table 1, respectively. They show that the QSS communities are strongly mutualistic, as was also observed for the stable states in Ref. KAUF91 . In contrast, the matrix elements of 22 feasible, but otherwise randomly chosen, communities FEAS ; TREG79  were found to be approximately uniformly distributed on $`[1,1]`$.

In Fig. 3 we show histograms of the multiplication rates Eq. (15) (i.e., the exponential function of the invasion fitness) of mutants that differ from the resident species by one single mutation (nearest-neighbor species, $`\mathrm{min}_JH_{iJ}=1`$), and by two mutations (next-nearest-neighbor species, $`\mathrm{min}_JH_{iJ}=2`$) ENDNOTE . Only a very small proportion of the nearest-neighbor mutants have a multiplication rate above unity \[Fig. 3(a)\]. (In our sample, this very small proportion came from just one of the ten QSS considered.) Among the next-nearest neighbors, on the other hand, a not insignificant proportion may be successful invaders \[Fig. 3(b)\]. The picture for third-nearest neighbors (not shown) is essentially the same as for the next-nearest neighbors. We also considered the multiplication rates for neighbors of the two most long-lived QSS observed in our simulations, which lasted about $`1.0\times 10^7`$ and $`1.4\times 10^7`$ generations, respectively. They both had no nearest neighbors with multiplication rates above unity, and the proportions for second-and third-neighbors were significantly smaller than those included in Fig. 3. In both parts of Fig. 3 are also shown the multiplication rates for neighbors of the 22 randomly chosen, feasible communities introduced in the previous paragraph. In that case, approximately half of the neighbors at any Hamming distance have multiplication rates above unity.

Summarizing the results from Table 1 and Fig. 3, we see that a long-lived QSS community is characterized by strongly mutualistic interactions and is surrounded by a “protection zone” of closely related genotypes that are very unlikely to successfully invade the resident community. The evidence from the two most long-lived QSS indicates that the average lifetime of a QSS is positively correlated with the extent of the protection zone. By contrast, a randomly chosen, feasible community has an approximately uniform distribution of $`M_{IJ}`$ over $`[1,1]`$. It also typically has no protection zone and so is more vulnerable to invasion than a QSS. It is clear from the fact that none of our randomly chosen feasible communities turned out to be a QSS that QSS are relatively rare in this model, even among feasible communities. We believe this may be a favorable condition for continuing evolution, as the ecology can move rapidly from QSS to QSS through series of unstable communities.

### IV.3 QSS lifetime statistics

The protection zone surrounding a QSS community acts much like the free-energy barrier that separates a metastable state in a physical system from the stable state or other metastable states. In both cases a sequence of improbable mutations or fluctuations is required to reach a critical state that will lead to major rearrangement of the system. It is thus natural to investigate in detail the lifetime statistics of individual QSS in the way common in the study of metastable decay RIKV94 ; RIKV94A . To this effect we started simulations from each of the ten QSS communities discussed in Sec. IV.2 and ran until the overlap function with the initial community,

$$𝒪(0,t)=\frac{_In_I(0)n_I(t)}{\sqrt{_In_I(0)^2_Jn_J(t)^2}},$$
(20)
became less than 0.5, at which time the system was considered to have escaped from the initial state. (We note that the overlap function defined here is simply the cosine of the angle between two unit vectors in the space of population vectors $`\{n_I\}`$.) The precise value of the cutoff used is not important, as long as it is low enough to exclude fluctuations during the QSS RIKV94A . The composition of the population at this time (the “exit community”) was then recorded and the run terminated. Each individual initial QSS was simulated for 300 independent escapes. The means and standard deviations of the individual lifetime distributions are shown in the last two columns of Table 1. The mean lifetimes were found to range over about one order of magnitude, from approximately 10 000 to 80 000 generations. In most cases the individual lifetime distributions were exponential (for which the standard deviation equals the mean), or they had an exponential component that described the behavior in the short-time end of the range of observed lifetimes. In about half of the QSS studied, a tail of very long lifetimes beyond the exponential part (indicated by a standard deviation significantly larger than the mean) was also observed. We further found that an initial QSS does not always escape via the same exit community. Typically, the genotypes in the exit community that are not present in the initial QSS differ from those in the original community by a Hamming distance of two or three. This picture is quite consistent with the stability properties of the QSS discussed in Sec. IV.2.

The range of the lifetimes of the ten QSS discussed above was limited: for a period to be identified by eye as a QSS, its lifetime could not be too short, while the lifetimes were bounded above by the length of the simulation run. The variation of a factor of ten within these limits indicates that the dynamical behavior of the system may display a very wide range of time scales. Another indication to this effect is provided by the details of the activity within periods that appear as high-activity in Figs. 1(a) and 2(a). Such detail is shown in Figs. 1(b) and 2(b). Statistically, the picture in these expanded figures is similar to the one seen on the larger scale, with shorter quiet periods punctuated by even shorter bursts of activity. These observations, which are similar to Refs. HALL02 ; CHRI02 , suggest statistical self-similarity at least over some range of time scales. Indications of self similarity have also been seen in fossil diversity records SOLE97 .

The suggestion of self-similarity that emerges from Figs. 1 and 2 makes it natural to investigate the statistics of the durations of active periods and QSS over a wider range of time scales than that represented by the ten QSS included in Table 1. From Fig. 1 we see that a reasonable method to make this distinction is to observe the magnitude of the entropy changes, $`|dS(t)/dt|`$, and to consider the system to be in the active state if this quantity is above some suitably chosen cutoff. This cutoff can best be determined from a histogram of $`dS/dt`$, such as shown in Fig. 4(a). It shows that the probability density of the entropy derivative consists of two additive parts: a near-Gaussian one corresponding to population fluctuations caused mostly by the birth/death process in the QSS, and a second one corresponding to large changes during the active periods. From this figure we chose the cutoff as 0.015, which was used to classify every sampled time point as either active or quiet. Normalized histograms for the durations of the active and quiet periods are shown in Fig. 4(b). About 97.4% of the time is spent in QSS. The active periods are seen to be relatively short, and their probability density is fitted well by an exponential distribution NOTEa . On the other hand, the lifetimes of the QSS show a very broad distribution CHRI02  – possibly a power law with an exponent near $``$2 for durations longer than about 200 generations. We note that there is some evidence of power laws with exponents near $``$2 in the distributions of several quantities extracted from the fossil record, such as the life spans of genera or families, and the sizes of extinctions SOLE96 ; NEWM96 ; SHIM02 ; SHIM03 . However, the fossil data are sparse and extend over no more than one or two decades in time, and they can be fitted almost as well by exponential distributions SOLE96 . Although our data extend over a much wider range of time scales than the paleontological data, indisputable evidence of a power law remains to be established. Though a few QSS of the order of one million generations will certainly appear in every run of $`2^{25}`$ generations, more definite conclusions about the statistics of such long QSS must await simulations an order of magnitude longer. Another feature in Fig. 4(b) is that the distribution of QSS durations changes to a smaller slope below about 200 generations. A similar effect has been observed in the fossil record of lifetimes of families SHIM02 ; SHIM03 .

### IV.4 Power spectral densities

From the discussion in Sec. IV.3 it is clear that by using a cutoff for the intensity of some variable which is large in the QSS and small in the active periods (but in both cases with fluctuations of unknown distribution), it is impossible to classify the periods unambiguously. For example, by increasing the cutoff one can make appear as a single QSS period what previously appeared as two successive, shorter QSS periods separated by a short active period. Since the probability of encountering an extremely large fluctuation in a QSS increases proportionally to the length of the QSS, this effect will affect the longest periods most severely. The problem with determining a suitable cutoff can be avoided by, instead of concentrating on period statistics, calculating the power spectral density (PSD) PSD ; PRES92  of variables such as the Shannon-Wiener diversity index or the total population. PSDs for these two quantities are shown in Fig. 5. The PSD for the diversity $`D(t)`$, shown in Fig. 5(a), appears inversely proportional to the frequency $`f`$ ($`1/f`$ noise or flicker noise MILO02 ) for $`f>10^4`$, goes through a crossover regime of about one decade in $`f`$ where it is $`f^\alpha `$ with $`\alpha >1`$, and then appears to return to $`\alpha 1`$ for $`10^6<f<10^5`$. For $`f<10^6`$ our data are insufficient to determine the PSD unambiguously, and much longer simulations would be necessary NOTEb . In the PSD for $`N_{\mathrm{tot}}(t)`$, shown in Fig. 5(b), the substantial population fluctuations due to the birth/death process during the QSS periods produce a large noise background which interferes with the interpretation at high frequencies. However, for $`f<10^5`$ the behavior is also consistent with $`1/f`$ noise.

On time scales much longer than the mean duration of an active period, the time series for the diversity can be approximated by constants during the QSS, interrupted by delta-function like spikes corresponding to the active periods (see Fig. 1). In this limit, the relation between a very wide distribution of the QSS durations $`\tau `$, described by a long-time power-law dependence of the probability density $`p(\tau )\tau ^\beta `$, and the low-frequency behavior of the PSD is known analytically PROC83 :

$$P(f)\{\begin{array}{ccc}f^{(\beta 1)}\hfill & \text{for}\hfill & \beta <2\hfill \\ 1/\left(f|\mathrm{ln}f|^2\right)\hfill & \text{for}\hfill & \beta =2\hfill \\ f^{(\beta 1)(3\beta )}\hfill & \text{for}\hfill & 2<\beta <3\hfill \\ |\mathrm{ln}f|\hfill & \text{for}\hfill & \beta =3\hfill \\ \mathrm{const}.\hfill & \text{for}\hfill & \beta >3\hfill \end{array}$$

Thus, the approximate $`1/f`$ behavior of the PSDs in Fig. 5 is consistent with the approximate $`1/\tau ^2`$ behavior of the probability density for the QSS durations, shown in Fig. 4(b).

As well as power-law distributions, PSDs that go as $`1/f^\alpha `$ with $`\alpha 1`$ have been extracted from the fossil record SOLE97 . However, such observed PSDs extend only over one to two decades in frequency (corresponding to power-law probability densities extending only over one or two decades in time), and more recent work indicates that the $`1/f`$ spectra obtained in Ref. SOLE97  may be artifacts of the analysis method NEWM99 . Although $`1/f`$ noise, at least in some frequency interval, is a property of our model and is possibly also seen in other models of macroevolution BAK93  (but see Refs. DAER96 ; DAVI01  for conflicting opinions on its presence in the Bak-Sneppen model), whether or not it is really present in the fossil record remains an open question.

### IV.5 Stationarity and effects of finite genome size

An important issue in evolutionary biology is whether or not the evolving ecosystem is stationary in a statistical sense. In Fig. 6 we show several diversity-related measures, averaged over a moving time window and independent simulation runs. These quantities are: the species richness $`𝒩(t)`$, the number $`𝒩_2(t)`$ of genotypes with $`n_I2`$, the Shannon-Wiener index $`D(t)`$, the total population $`N_{\mathrm{tot}}(t)`$, and the average Hamming distance $`H`$ between genotypes in the community and its standard deviation $`\sigma _H`$. As seen in the figure, none of these quantities show any signs of a long-time trend over $`2^{25}`$ generations. This is consistent with fossil evidence for constant diversity MAYN89 . However, it is in disagreement with simulations of the original version of the tangled-nature model HALL02 , in which a slow growth in species richness was observed. This discrepancy might possibly result from the different forms of the interaction matrix used in the two studies, but it may also be due to the relatively short time series used in Ref. HALL02 .

For our model, a simple, combinatorial phase-space argument (containing neither the individual population sizes nor the mutation rate) indicates that the most probable number of major genotypes should indeed be limited. The argument goes as follows. A community of $`𝒩`$ genotypes can be chosen in $`\left(\genfrac{}{}{0pt}{}{2^L}{𝒩}\right)`$ different ways and is influenced by $`𝒩(𝒩1)/2`$ different pairs of $`M_{IJ}`$ and $`M_{JI}`$. We let the probability that the pair of interactions is suitable to forming a stable community be $`q`$ (if the requirement is simply that both $`M_{IJ}`$ and $`M_{JI}`$ should be positive, then $`q=1/4`$ with our choice of $`𝐌`$). Thus, a stable $`𝒩`$-species community can be formed in approximately

$$\mathrm{\Omega }(𝒩)=\left(\genfrac{}{}{0pt}{}{2^L}{𝒩}\right)q^{𝒩(𝒩1)/2}$$
(21)
ways, which obeys the recursion relation

$$\mathrm{\Omega }(𝒩+1)=\frac{2^L𝒩}{𝒩+1}q^𝒩\mathrm{\Omega }(𝒩)$$
(22)
with $`\mathrm{\Omega }(1)=2^L`$. This recursion relation can either be used to find the most probable value $`𝒩^{}`$ numerically, or one can easily obtain an estimate valid for $`2^L𝒩`$: $`𝒩^{}L\mathrm{ln}2/\mathrm{ln}(1/q)`$. The value numerically obtained with $`q=1/4`$ is $`𝒩^{}=6`$, but the close correspondence to the results shown in Fig. 6 may be fortuitous since the dependence of $`𝒩^{}`$ on $`q`$ is only logarithmic.

Another question of interest is whether the limited size of the genome leads to “revisiting” of genotypes and communities. The answer is affirmative and indicates the importance of further studies of the effects of the genome size on the long-time dynamical behavior. For genotypes the question of revisits is easy to answer, and the results for a single run of $`2^{25}`$ generations are shown in Fig. 7. As we see, in less than $`3\times 10^6`$ generations, at least two individuals of every genotype have appeared simultaneously (curve labeled $`n_I>1`$). By the end of the run, almost every genotype has enjoyed being a major species with $`n_I1001`$ at least once. Plots on a log versus linear scale (not shown) indicate that the curves are reasonably well approximated by exponential approaches to $`𝒩_{\mathrm{max}}`$, indicating that genotypes appear to be nearly randomly visited and revisited at a constant rate dependent on $`n_I`$.

For communities, the question is more difficult. If a community is simply defined as an internally stable $`𝒩`$-species community, then an estimate for their total number can be obtained from $`\mathrm{\Omega }(𝒩)`$ of Eq. (21) as described in Appendix B. The result with the parameter values used in this paper is of the order of $`10^{12}`$. However, unstable communities are also briefly visited, especially during active periods, and a more inclusive way of counting communities would be a coarse-graining procedure based on the overlap function defined in Eq. (20). The nonzero populations are recorded at $`t=0`$, and the overlap with this initial community is monitored until at some $`t^{}>0`$ it falls below a suitably chosen cutoff $`𝒪_{\mathrm{cut}}`$. The set $`\{n_I(t^{})\}`$ is then recorded as a new community with $`t^{}`$ as its starting time. The process is repeated, now comparing with the newly recorded community. This procedure creates a list of communities and their starting times. To address the issue of revisits, we next scan this list to extract those communities that represent revisits to a previously visited community, or “prototype community.” We compare each community $`\{n_I(t)\}`$ in the original list with the previously found communities $`\{n_I(\stackrel{~}{t})\}`$, sequentially in order of increasing $`\stackrel{~}{t}<t`$. If none of the overlaps is greater than $`𝒪_{\mathrm{cut}}`$, $`\{n_I(t)\}`$ is added to the list of prototypes $`_\mathrm{p}(t_\mathrm{p})`$ with $`t_\mathrm{p}=t`$. If an overlap greater than $`𝒪_{\mathrm{cut}}`$ is found at some $`\stackrel{~}{t}`$, we stop the comparisons and add this community to the list of revisits $`_\mathrm{r}(t_\mathrm{r})`$. The starting time of this revisited community, $`t_\mathrm{r}=t`$, is associated in the list $`_\mathrm{r}(t_\mathrm{r})`$ with the unique $`t_\mathrm{p}=\stackrel{~}{t}`$, the starting time of the associated prototype. As the system evolves, the number of items in each list increases monotonically. Clearly, the cutoffs should be sufficiently small that communities that differ only by fluctuations inside a QSS are not considered different, but sufficiently large to avoid counting significantly different communities as identical. From inspection of the overlap fluctuations in some of the QSS included in Table 1, we found that $`𝒪_{\mathrm{cut}}`$ in the range of 0.90 to 0.95 is optimal. Estimates of the total numbers of communities (stable and unstable) that can in principle be distinguished by this method are obtained in Appendix B. Depending on details of the assumptions, the estimates vary between $`10^{21}`$ and $`10^{28}`$ for these cutoffs.

The results of the procedure described above with $`𝒪_{\mathrm{cut}}=0.90`$ and $`0.95`$ are shown in Fig. 8 for a single run of $`2^{25}`$ generations. The upper pair of curves in Fig. 8(a) corresponds to $`𝒪_{\mathrm{cut}}=0.95`$, and the lower pair to $`𝒪_{\mathrm{cut}}=0.90`$. In each pair, the upper curve shows the total number of communities in the lists $`_\mathrm{p}\left(t_\mathrm{p}\right)`$ and $`_\mathrm{r}\left(t_\mathrm{r}\right)`$, while the lower curve shows just the number of prototypes in $`_\mathrm{p}\left(t_\mathrm{p}\right)`$. For $`𝒪_{\mathrm{cut}}=0.95`$, about 40% of the communities are seen to be revisits, while the proportion is about 30% for $`𝒪_{\mathrm{cut}}=0.90`$. In both cases, the curve for the number of prototypes remains approximately linear, indicating that the supply of previously unvisited communities is nowhere near to being exhausted, even for such a long run. In view of the enormous numbers of available communities estimated above, this result is reasonable. QSS appear as plateaus in the curves showing the numbers of prototypes. Figure 8(b) provides a different perspective (for $`𝒪_{\mathrm{cut}}=0.95`$ only). Each revisit is represented by a point with $`t_\mathrm{r}`$ as abscissa and $`t_\mathrm{p}`$ as ordinate. Thus, every point lies strictly below the diagonal. Long-lived QSS appear as large gaps and horizontal segments (e.g., the one near $`10^7`$ generations). The inset shows a detail of $`10^6`$ generations near the diagonal around $`1.65\times 10^7`$. As we see, there are many points just below the diagonal, representing the fluctuations around the (many short-lived) QSS. By contrast, points far below the diagonal represent “throwbacks” to the vicinity of earlier prototype communities. Note that the density of points is much higher just below the diagonal, implying that a large portion of the revisited communities are “fluctuation related.” To highlight these differences, we show the cumulative probability of the ratio $`t_\mathrm{p}/t_\mathrm{r}`$ in Fig. 8(c), in which the lower curve corresponds to the data in Fig. 8(b). For the case of $`𝒪\mathrm{cut}=0.95`$ (lower curve), we see that about 60% of the revisits can be regarded as “fluctuation related” (with $`t_\mathrm{p}t_\mathrm{r}`$), while the rest are “throwbacks.” Roughly, the latter component appears to be distributed uniformly over all earlier times. The upper curve shows the corresponding result for $`𝒪_{\mathrm{cut}}=0.90`$. Corresponding to using a more coarse-grained covering of state space, it naturally displays a larger proportion of throwbacks. Not surprisingly, there is no sharp distinction between these two components of the revisited communities since even a rough partition depends on the details of coarse graining. Nevertheless, we can conclude that the dynamics produces a steady stream of essentially new communities drawn from the vast supply of possibilities.

## V Summary and Conclusions

In this work we have studied, by linear stability analysis and large-scale Monte Carlo simulations, a simplified version of the tangled-nature model of biological coevolution, recently introduced by Hall, Christensen, and collaborators HALL02 ; CHRI02 ; COLL03 . Selection is provided by interspecies interactions through the reproduction probability $`P_I`$ \[Eq. (1)\], which corresponds to a nonlinear population-dynamics model of the community ecology, while the genetic variability necessary for evolution is provided by a low rate of mutations that act on individual organisms during reproduction.

At the low mutation rate studied here, the model provides an intermittent, statistically self-similar behavior, characterized by periods of relative calm, interrupted by bursts of rapid turnover in genotype space. During the quiet periods, or quasi-steady states (QSS), the population consists of a community of a relatively small number of mutualistically interacting genotypes. The populations of the individual genotypes, $`n_I(t)`$, fluctuate near a stable fixed point of a deterministic mean-field, mutation-free version of the model \[Eq. (4)\]. During the active periods, the system moves through genotype space at a rate that is rapid on a macroscopic (“geological”) time scale, although of course finite on a microscopic (“ecological”) scale. These periods of rapid change are characterized by large fluctuations in the diversity and an overall reduction of the total population. In long simulation runs of $`2^{25}`$ generations, the ecosystem spends about 97.4% of the time in QSS. The time series produced by the model are statistically stationary, and there is no evidence that any particular quantity is being optimized as the system moves through genotype space. In that sense, the dynamical behavior can probably best be described as a neutral drift DROS01 . Overall, the dynamical behavior of this model resembles closely the punctuated-equilibria mode of evolution, proposed by Gould and Eldredge GOUL77 ; GOUL93 ; NEWM85 .

Investigation of duration statistics for the quiet periods shows a very wide distribution with a power-law like long-time behavior characterized by an exponent near $`2`$. Consistent with this result, the power spectral density of the diversity shows $`1/f`$ noise. While there are claims that similar statistics characterize the fossil record SOLE97 ; SOLE96 ; NEWM96 ; SHIM02 ; SHIM03 , this is still a contested issue SOLE96 ; NEWM99 . At best, observations of power-law distributions and $`1/f`$ noise in the fossil record extend over no more than one or two decades in time or frequency, and it must remain an open question whether this is the optimal interpretation of the scarce data available.

Due to the absence of sexual reproduction, our model can at best be applied to the evolution of asexual, haploid organisms such as bacteria. It should also be noted that no specific, biologically relevant information has been included in the interaction matrix $`𝐌`$. In particular, this fact may be responsible for our QSS being strongly dominated by mutualistic relationships. The absence of biologically motivated detail in $`𝐌`$ is both a strength and a weakness of the model. Its strength lies in reinforcing the notion of universality in macroevolution models, e.g., power law behaviors and $`1/f`$ noise. By the same token, its weakness lies in the lack of biological detail, to the point that comparison with specific observational or experimental data is difficult. Clearly, the detailed effects of interspecies interactions on the macroevolutionary behavior in models similar to the one studied here represents an important field of future research. Examples include the importance of the connectivity of the interaction matrix, correlated interspecies interactions SEVI04 , and interaction structures corresponding to food webs with distinct trophic levels.

Despite all these caveats, we find it encouraging that such a simple model of coevolution with individual-based dynamics can produce punctuated equilibria, power-law distributions, and $`1/f`$ noise consistent with current theories of biological macroevolution. We believe future research should proceed in the direction pointed out by this and similar models. This entails combining stochastic models from community ecology with models of mutations and sexual reproduction at the level of individual organisms, and investigating the consequences of more biologically realistic interspecies interactions.

## Acknowledgments

We are grateful for many helpful comments on the manuscript by J. A. Travis, and we acknowledge useful conversations and correspondence with I. Abou Hamad, G. Brown, T. F. Hansen, N. Ito, H. J. Jensen, A. Kolakowska, S. J. Mitchell, M. A. Novotny, H. L. Richards, B. Schmittmann, V. Sevim, U. Täuber, and J. C. Wilgenbusch.

P. A. R. appreciates the hospitality of the Department of Physics, Virginia Polytechnic Institute and State University, and the Department of Physics and Astronomy and ERC Center for Computational Sciences, Mississippi State University.

The research was supported in part by National Science Foundation Grant Nos. DMR-9981815, DMR-0088451, DMR-0120310, and DMR-0240078, and by Florida State University through the Center for Materials Research and Technology and the School of Computational Science and Information Technology.

## Appendix A Master equation

A complete description of the stochastic process can be given in terms of a master equation, which specifies the evolution of $`𝒫(\stackrel{}{n},t)`$, the probability that the system is found with composition $`\stackrel{}{n}`$ at time $`t`$. Here, $`\stackrel{}{n}\{n_1,n_2,\mathrm{},n_{𝒩_{\mathrm{max}}}\}`$, and $`𝒩_{\mathrm{max}}=2^L`$ is the number of individuals of species $`I`$. In our case, $`L=13`$ and $`𝒩_{\mathrm{max}}=8192`$. Similar to the main text, we define $`N_{\mathrm{tot}}_In_I`$ and let $`N_0`$ be the carrying capacity.

We write the probability for an individual of species $`I`$ to survive to reproduce as

$$P_I\left(\stackrel{}{n}\right)=\left\{1+\mathrm{exp}\left[\frac{N_{\mathrm{tot}}}{N_0}\underset{J}{}M_{IJ}\frac{n_J}{N_{\mathrm{tot}}}\right]\right\}^1.$$
(23)
The main difference between this expression and Eq. (1) lies in the interpretation. Here, $`\stackrel{}{n}`$ is a “coordinate variable” in the $`𝒩_{\mathrm{max}}`$-dimensional space, in contrast to $`n_I\left(t\right)`$ being just a point in this space.

To proceed, we define the symbol $`\left[\genfrac{}{}{0pt}{}{n_I}{m_I}\right]`$ as the rate for survival (from $`n_I`$ individuals to $`m_I`$ “mothers”). Since each individual is given a chance to survive according to $`P_I`$, we have

$$\left[\genfrac{}{}{0pt}{}{n_I}{m_I}\right]=\frac{n_I!}{m_I!\left(n_Im_I\right)!}\left(P_I\right)^{m_I}\left(1P_I\right)^{n_Im_I},$$
(24)
which has the form of a binomial distribution. Next, each mother gives rise to $`F`$ offspring. However, due to mutations, not every offspring is of the same genotype as the mother. Although it is possible to have mutants with a genotype differing from the mother by more than one bit, we restrict ourselves here to a simpler version, namely mutant genotypes that can differ only by one bit. Since our simulations typically involve $`\mu 10^3`$, this restriction should not lead to serious difficulties. Given that only one bit may be flipped, there are $`L+1`$ possible varieties of offspring for each maternal genotype. We introduce the notation

$`b_{J,0}\text{for the number of offspring from mother }J\text{ with no mutations}`$
$`b_{J,k}\text{for the number of offspring from mother }J\text{ with the }k\text{th bit flipped}`$
We now define the multinomial-like symbol

$$\left[\genfrac{}{}{0pt}{}{Fm_J}{b_{J,0},b_{J,1},\mathrm{},b_{J,L}}\right]=\frac{\left(Fm_J\right)!}{\left(b_{J,0}\right)!}\left(1\mu \right)^{b_{J,0}}\underset{k=1}{\overset{L}{}}\frac{1}{\left(b_{J,k}\right)!}\left(\frac{\mu }{L}\right)^{b_{J,k}},$$
(25)
which is the probability that the $`Fm_J`$ offspring are distributed into the specific set $`\{b_{J,0},b_{J,1},\mathrm{},b_{J,L}\}`$. The last ingredient needed is the connection matrix

$$\mathrm{\Delta }_G^{J,k}=\{\begin{array}{cc}1& \text{if genotype }G\text{ is }J\text{ with the }k\text{th bit flipped}\\ 0& \text{otherwise}\end{array}$$
(26)
so that the number of offspring born into species $`G`$ due to mutations is

$$\stackrel{~}{b}_G\underset{J,k}{}\mathrm{\Delta }_G^{J,k}b_{J,k}.$$
(27)
The final result is

$$𝒫(\stackrel{}{n}^{},t+1)=\underset{\stackrel{}{n},\stackrel{}{m},\left\{b\right\}}{}\underset{G}{}\delta \left(n_G^{}b_{G,0}\stackrel{~}{b}_G\right)\underset{J}{}\left[\genfrac{}{}{0pt}{}{Fm_J}{b_{J,0},b_{J,1},\mathrm{},b_{J,L}}\right]\underset{I}{}\left[\genfrac{}{}{0pt}{}{n_I}{m_I}\right]𝒫(\stackrel{}{n},t).$$
(28)
In principle, once the dynamics of $`𝒫`$ is found, the time dependence of various quantities (e.g., averages and correlations) can be computed. For instance,

$$n_I_t\underset{\stackrel{}{n}}{}n_I𝒫(\stackrel{}{n},t)$$
(29)
is the average number of individuals of species $`I`$ at time $`t`$. From this very detailed description, we see that the evolution of various quantities can be derived. However, these evolution equations are extremely complex in general and progress is usually possible only by making a “mean-field” approximation, in which all correlations are neglected. Thus, averages of products of $`\stackrel{}{n}`$ are replaced by the products of the averages. As a result, a non-linear evolution equation for $`\stackrel{}{n}_t`$ results. In a subsequent paper ZIA04 , we shall show that Eq. (2) is precisely the result of such an approach (with $`\stackrel{}{n}_t`$ denoted by $`\{n_I\left(t\right)\}`$). We shall also demonstrate that, for a QSS, $`𝒫\left(\stackrel{}{n},t\mathrm{}\right)`$ is well approximated by a Gaussian (for $`\mu ,1/N_01`$), the center of which is just the fixed point given by Eq. (12). The width of this Gaussian is governed in part by the community matrix $`\stackrel{~}{𝚲}`$ of Eq. (14). In addition to this (systematic, “drift”) part, there is a part governing the noise correlations. Between the two, we can compute distributions of step sizes ($`\stackrel{}{n}^{}\stackrel{}{n}`$). For a long-lived QSS, it is easy to compile data so that quantitative comparisons with these predictions can be made.

## Appendix B Counting communities

In this Appendix we obtain estimates for the numbers of different communities that could in principle be observed in an infinitely long simulation. Even our most conservative estimate is vastly larger than the numbers observed in a simulation run of $`2^{25}`$ generations in Sec. IV.5.

If a community is simply defined as an internally stable $`𝒩`$-species community, then an estimate for their total number can be obtained from $`\mathrm{\Omega }(𝒩)`$ of Eq. (21). Since $`\mathrm{\Omega }(𝒩)`$ is quite sharply peaked around $`𝒩^{}`$, we can estimate that $`_𝒩\mathrm{\Omega }(𝒩)\mathrm{\Omega }(𝒩^{})2^{L𝒩^{}}q^{𝒩^{}(𝒩^{}1)/2}\left(2^{a/2}\right)^{aL^2}`$ with $`a=\mathrm{ln}2/\mathrm{ln}(1/q)`$. For $`q=1/4`$ and $`L=13`$, this gives $`\mathrm{\Omega }(𝒩^{})5.2\times 10^{12}`$, while direct summation of $`\mathrm{\Omega }(𝒩)`$ with the same parameters gives $`_𝒩\mathrm{\Omega }(𝒩)=8.5\times 10^{11}`$ with 46% due to the contribution from $`𝒩=𝒩^{}=6`$ and convergence by $`𝒩=11`$. As it does not include unstable communities, this would serve as our most conservative estimate of the total number of communities that could be visited in an infinitely long simulation.

To estimate the total number of different communities that can be distinguished by the overlap cutoff method described in Sec. IV.5, we note that the overlap threshold $`𝒪_{\mathrm{cut}}`$ is simply the cosine of the angle between two vectors in population space, $`𝒪_{\mathrm{cut}}=\mathrm{cos}\theta `$. For convenience, we define $`ϵ1𝒪_{\mathrm{cut}}`$, such that $`ϵ`$ is small ($`0.1`$ and $`0.05`$). Then, $`\theta =\mathrm{arccos}(1ϵ)\sqrt{2ϵ}`$. The fact that the number of populated species at any one time is near $`𝒩^{}`$, means that the direction of the population vector $`\{n_I\}`$ to a reasonable approximation lies within the first hyper-octant of one of $`\left(\genfrac{}{}{0pt}{}{𝒩_{\mathrm{tot}}}{𝒩^{}}\right)`$ $`𝒩^{}`$-dimensional hyperspheres. The total number of communities that can be distinguished with a cutoff $`𝒪_{\mathrm{cut}}`$ is the ratio of the total $`𝒩^{}`$-dimensional solid angle covered by the hyper-octants to that of a hyper-cone subtended by $`\theta `$. The solid angle of the hyper-octants is

$$V_{\mathrm{tot}}\left(\genfrac{}{}{0pt}{}{𝒩_{\mathrm{tot}}}{𝒩^{}}\right)\frac{S_𝒩^{}}{2^𝒩^{}},$$
(30)
where $`S_d=2\pi ^{d/2}/\mathrm{\Gamma }(d/2)`$ is the surface area of a $`d`$-dimensional unit sphere. The solid angle of the hyper-cone is

$$C(ϵ)=S_{𝒩^{}1}_0^\theta (\mathrm{sin}\beta )^{𝒩^{}2}S_{𝒩^{}1}\frac{(2ϵ)^{(𝒩^{}1)/2}}{𝒩^{}1}$$
(31)
for small $`ϵ`$. Using Stirling’s approximation for the gamma functions in $`S_d`$, $`\mathrm{\Gamma }(z)(z/e)^z\sqrt{2\pi /z}`$ and the identity $`lim_n\mathrm{}\left(1\frac{1}{n}\right)^n=e^1`$, we estimate the total number of different communities to be

$$M(ϵ)\left(\genfrac{}{}{0pt}{}{𝒩_{\mathrm{max}}}{𝒩^{}}\right)\frac{\sqrt{2\pi 𝒩^{}ϵ}}{(2ϵ)^{𝒩^{}/2}}.$$
(32)
This estimate is reasonable only as long as $`ϵ`$ is large enough to exclude minor fluctuations within QSS. Observation of overlap fluctuations in some of the QSS included in Table 1 indicates that $`𝒪_{\mathrm{cut}}`$ in the range of 0.90 to 0.95 is optimal. With the parameters used in this work, Eq. (32) yields $`M10^{21}`$ for $`𝒪_{\mathrm{cut}}=0.90`$ and $`M10^{22}`$ for 0.95. If, instead of using $`𝒩^{}=6`$ we assume that the number of genotypes at any time is near 8, as indicated by the top curve in Fig. 6, these estimates are increased by six orders of magnitude.