{'fraction_non_alphanumeric': 0.030266427388183324, 'fraction_numerical': 0.00924903368304804, 'mean_word_length': 4.694045981528787, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Effects of correlated interactions in a biological coevolution model with individual-based dynamics

Models of coevolution have recently been studied, in which a species is defined by a genome in the form of a bitstring, and the interactions between species i and j are given by a fixed matrix with independent, randomly distributed elements M_{ij}. A consequence of the stochastic independence is that species whose genotypes differ even by a single bit may have completely different phenotypes. This is clearly unrealistic, as closely related species should be similar in their interactions with the rest of the ecosystem. Here we therefore study a model, in which the M_{ij} are correlated to a controllable degree. We calculate, both analytically and numerically, the correlation function for matrix elements M_{ij} and M_{kl} versus the Hamming distance between the bitstrings representing the species. We compare Monte Carlo simulations of coevolution models with uncorrelated and correlated interactions. In particular, we consider the lifetimes of individual species. The species-lifetime distribution is close to a power law with an exponent near -2 in both uncorrelated and correlated cases. The durations of quasi-steady states and power spectral densities for the diversity indices display noticeable differences. However, some qualitative features, like 1/f behaviour in power spectral densities for the diversity indices, are not affected by the correlations in the interaction matrix.


Introduction
During the last decade, ecological dynamics and biological evolution have become new areas of interest for statistical physicists [1]. This trend has merged with recent developments in the theory of complex networks [2], which has drawn increased attention to problems in ecological dynamics, such as the evolution and stability of food webs [1,3,4,5,6].
While many models of macroevolution are formulated explicitly at evolutionary (macroscopic) time scales, evolutionary dynamics in nature are driven by reproduction, mutations, and selection at the ecological (microscopic) scale of individual organisms. Although this means that individual-based models of macroevolution must span a dauntingly large range of timescales, this is now becoming possible with the aid of modern computers and algorithms. We therefore recently studied an individual-based biological coevolution model [7,8,9] with random interspecies interactions. This model is a simplified version of the tangled-nature model introduced by Jensen, et al. [10,11,12]. Individuals are represented by a haploid genome consisting of a string of bits. The individuals reproduce asexually and interact with individuals of other species through a random interaction matrix with stochastically independent elements, which determine how the reproduction probability of individuals of each species is affected by the presence and abundance of other species in the "ecosystem". During reproduction, an offspring individual can randomly mutate with a small but constant probability by having a single bit in its genome flipped. These mutations create perturbations in the ecosystem by introducing new species. A particular mutant may or may not succeed, depending on its interactions with the species already present in the system. Usually a mutant does not fit in well with the resident species, and so it quickly goes extinct. However, occasionally a successful mutant can cause avalanches of mass extinctions, destroying some or all of the other species through predation or competition. As a consequence, the individuals in this model ecosystem live in a dynamic "fitness landscape", which evolves due to changes in the populations of the resident species and the introduction of new species through mutations. The model displays quiet periods interrupted by bursts of high activity caused by mass extinctions that are triggered by the introduction of new mutants, reminiscent of the punctuated-equilibrium mode of evolution suggested by Eldredge and Gould [13,14,15]. The lifetime distribution for individual species is well described by a power law with exponent close to −2.
A rather unrealistic aspect of this model is the fact that the elements of the interaction matrix are completely uncorrelated. This means that a mutant in general will have completely different interactions with other species, and thus a completely different phenotype, than does its parent "wildtype". In the present paper, we therefore modify our previous model by introducing correlations into the interaction matrix. (Similar ideas have been pursued by Kauffman [16,17].) In the modified model, the correlations between matrix elements ensure that mutants are not completely different from their parents. Rather, their interactions with the other species are similar to those of the wildtype. The price to be paid is a significant reduction in the number of effectively independent phenotypes for a genome of a given length. In this paper we compare numerical results for this modified model with results for the corresponding uncorrelated model in order to assess the effects of including correlations. The results show that, although long-range correlations affect the dynamics of the system noticeably, some qualitative features remain unchanged.
The organization of the rest of this paper is as follows. In section 2 we describe the model and the algorithm used for the simulations. In section 3 we describe the construction of the correlated interaction matrix and briefly discuss some complications of this process. In section 4 we report our results, which are summarized in section 5. Derivations of our analytical results are provided in the Appendices.

Model
A species in our Monte Carlo (MC) simulations is represented by a bitstring genome of length L. This L-bit genome supplies a pool of 2 L possible species. These organisms are considered haploids. Thus, the terms "species", "species index", and "genotype" are synonymous in this paper.
Individuals in the simulation are allowed to give birth to F offspring per generation with a probability P i , where i denotes the species index. The reproduction is asexual (cloning), and it occurs in discrete, nonoverlapping generations. Thus, only offspring can survive to the next generation, and all individuals die at the end of the generation whether they reproduce or not. We use a mathematically convenient, nonlinear form for the probability of reproduction [7,8,10,11,18], given by where n i (t) is the population of species i at time t. The Verhulst factor N 0 represents the carrying capacity of the "ecosystem" and keeps the total population, N tot (t) = i n i (t), finite. M is the interaction matrix, in which a matrix element M ij represents the effect of the population density of species j on species i. Species i benefits from species j if M ij is positive, and it is inhibited by species j if M ij is negative. There is no interaction between species i and j when M ij is zero. The structure of M is discussed in the next section.
Although we begin the simulations with a population consisting of only one species, the system diversifies quickly because, in each generation, all individuals are exposed to mutations that create new species. The mutation process is modelled by flipping a randomly chosen bit in the individual's genome with a probability µ per generation per individual.
We note that the population dynamics used in the model is somewhat unrealistic, due to the absence of external energy resources and conservation of energy (or biomass) [7,8,9,10,11,12].

Form of the Matrix
The structure of the interaction matrix characterizes the dynamics of the system. Rikvold and Zia [7,8] studied the same system (except for a negligible probability of O(µ 2 ) for multiple mutations) using an interaction matrix with off-diagonal elements distributed randomly and uniformly over [−1, 1] and with the diagonal elements set to zero. We shall call this the uniform-uncorrelated model. For reasons that will become clear later, we instead use a Gaussian distribution to facilitate some comparisons. We call this model the Gaussian-uncorrelated (or just the uncorrelated) model.
The main focus of this paper is the effects of correlations in the interaction matrix. (Some preliminary results of this work were discussed in reference [18].) The motivation behind the introduction of correlations is the following. In reality, the interaction between two species x and z should be positively correlated with the interaction between species y and z if x and y are closely related. Therefore, one would expect a positive correlation between the matrix elements M xz and M yz .
A correlated matrix is generated using the following method. First, a random matrix with elements drawn from a Gaussian distribution with standard deviation σ 0 = 1/3 and mean zero is generated. (We chose this standard deviation for backward compatibility with earlier studies in which the matrix elements were uniformly distributed on [−1, 1].) Then, each matrix element in the random matrix is averaged over itself and its neighbours in genome space up to a Hamming distance n, where n is the averaging radius. Thus, the total number of terms in the average (or the number of points inside the averaging "hypersphere" with radius n in genome space) is Z n = n m=0 2L m . The average is multiplied by √ Z n in order to keep the standard deviation of the correlated matrix elements approximately the same as for the random ones. This process eliminates any possible effect by a change in the shape of the distribution of interaction strengths. Thus, a correlated matrix element M ij is given by where M 0 kl are the elements of the uncorrelated matrix, which are Gaussian distributed with standard deviation σ 0 , and H(ij; kl) is the Hamming distance between two matrix elements M 0 ij and M 0 kl . We use the city-block metric to define the Hamming distance between two matrix elements. Therefore, is the Hamming distance between bitstrings i and j. In fact, H(ij; kl) denotes the Hamming distance between the concatenated bitstrings ij and kl. (This concatenation notation is useful for Hamming distance calculations and will be used throughout the paper.) The diagonals are set to zero after the averaging process is finished. The model in which such a correlated interaction matrix is used is called the correlated model. Theoretical calculations of correlation functions are discussed in Appendix A.
There are two modelling issues that need to be discussed here. The first one is the interpretation of the mutation process. As we mentioned above, in the uncorrelated model, matrix elements M 0 xz and M 0 yz are not correlated, even if y is a mutant of x. (We call y a mutant of x if H(x, y) = 1.) Therefore, the set of interactions between the mutant y and other species can be completely different than those between the wildtype x and the other species. This lack of relationship between the mutant and the wildtype is not very realistic. Therefore, the mutation process in the uncorrelated model may, in some sense, be interpreted as the introduction of a completely new species to the system from a pool of species, reminiscent of migration. In the correlated model, the interaction constants of the mutant and the wildtype are correlated. However, whether these correlations are strong enough to describe those between a real mutant and its wildtype is arguable. The other issue is the non-antisymmetrical form of M, meaning M ij = −M ji . This feature of the model, along with the absence of an external energy source, gives rise to unrealistic numbers of mutualistic interactions in the system [7].
We have already made a preliminary comparison between the correlated model and the Gaussian-and uniform-uncorrelated models in reference [18] for short-range correlations, n = 1 and n = 2. The results suggested that the correlated and Gaussianuncorrelated models behave quite similarly, so that short-range correlations in M may not have a significant effect on the dynamics. In this paper, we extend our investigation to the effects of longer-range correlations with n > 2.

Distributions of Matrix Elements, Correlations Between Matrix Elements, and Complications in the Averaging Process
The averaging process mentioned above causes a few unexpected complications. The first one is a shift in the mean. Although the matrix elements are drawn from a distribution with mean zero, for a finite matrix, the mean, M 0 ij , while quite small, is always non-zero. The averaging alters the mean of the matrix elements, so that the means of the initial random and the averaged matrices can be considerably different. The reason for the shift is the following. Since each averaged element is multiplied by √ Z n to keep the standard deviation constant, the mean of the elements of the averaged interaction matrix becomes ij is a negligibly small number, multiplication by a large √ Z n could cause a considerable shift, M ij − M 0 ij , in the mean. In order to prevent this, we adjust the mean, M 0 ij , by subtracting M 0 ij /2 2L from each element of the initial random matrix. This modification minimizes the shift in the mean.
The other complication of the averaging process is a change in the standard deviation of M ij , even though we intend to keep it constant. As in the problem of the shift in the mean mentioned above, the standard deviation of the averaged matrix can be considerably different than that of the initial random matrix. As seen in figure 1, the matrices averaged for n = 1 have approximately the same Gaussian form with a standard deviation close to σ 0 , which practically overlaps with the M 0 ij distribution of the initial Gaussian-uncorrelated matrix. However, we found that long-range averaging (n L/2) could generate distributions with significantly different standard deviations, like the n = 7 curve included in figure 1. The reason for the change in the standard deviation is very similar to that for the shift in the mean mentioned above. This anomaly is discussed further in Appendix B.
We also checked the correlations between matrix elements to see the effect of averaging. The correlations between two matrix elements, M ij and M i ′ j ′ , depend on the overlap of the M 0 kl terms that occur in the average in equation (2). Theoretical calculations predict decaying correlation functions with steps due to the city-block metric, as shown in figure A1. These calculations are explained in Appendix A.
The correlation functions of the matrices with short-range correlations are in good agreement with the theory. However, as n is increased, the numerical correlation functions begin to deviate from their theoretical counterparts, as shown in figure 2. The numerical correlation functions of the highly correlated matrices are distorted and translated along the y-axis, compared to the theoretical ones. We believe that the cause of this anomaly is related to the cause of the anomaly in the standard deviations. This complication is discussed in Appendix B.

Simulation Results
Unless otherwise indicated, we performed sixteen independent runs for 2 25 =33 554 432 generations for each model using parameters similar to those in references [7] and [18]: genome length L = 8 bits, mutation rate µ = 10 −3 per individual per generation, carrying capacity N 0 = 2000, fecundity F = 4; and for the correlated model, averaging radii n = 1, 3, 5, and 7. Each run used a different M and a different sequence of random numbers. Each simulation started with 200 individuals of genotype 0. The system parameters were recorded every sixteen steps. In order to see the effect of the shape of the M 0 ij distribution, we first ran four sets of simulations of the uncorrelated model with a uniform M 0 ij distribution with standard deviation σ 0 , and with a Gaussian M 0 ij distribution with standard deviations σ 0 /2, σ 0 and 2σ 0 . For these systems with different M 0 ij distributions, we compared histograms of species lifetimes. The species lifetime is defined as the number of generations between creation and extinction of a species, during which its population is continuously positive.
As seen in figure 3, the lifetime distributions of the uncorrelated systems with different matrix-element distributions all exhibit approximately power-law like decay over seven decades in time, with an exponent near −2. Figure 3 also shows that a narrower M 0 ij distribution flattens the concavity between 10 and 10 4 generations in the lifetime distribution, while a wider one makes it deeper.
The lifetime distributions of the correlated models shown in figure 4 also exhibit a   figure 3, suggesting that the change in the standard deviation (or the shape) of the M ij distribution could perhaps be responsible for this anomaly.
In order to make sure that it is actually the correlations that cause the changes, we suppressed this anomaly: we ran another set of simulations for n = 7 in which the standard deviation of the M ij distribution was adjusted to σ 0 after averaging. The result is shown in the inset in figure 4. The two lifetime distributions in the inset almost completely coincide. They are both the averages of eight-run sets with n = 7 (using the same random number sequences), except that the M ij distribution of the second set is adjusted after averaging. The very small difference between these two curves shows that the effect of a change in the standard deviation of the M ij distribution after averaging is negligible. Thus, the correlations in M are largely responsible for the differences between the lifetime distributions of the correlated and uncorrelated models, shown in the main part of figure 4.  The same sets of random numbers were used for the "adjusted" eight-run set as for their "non-adjusted" counterparts. The distributions practically overlap, suggesting that the change in the standard deviation has much less effect on the lifetime distribution than the correlations in M.
The deviation trend in the lifetime distributions is not monotonic in the averaging radius (the n = 5 curve appears below the n = 3 curve). However, the error bars for the n = 3 and n = 5 curves overlap. The large error margin due to the anomaly in the correlation functions leaves some uncertainty. Nevertheless, even fairly longrange correlations in M (like n = 5 and 7) do not change the gross features like the approximate 1/t 2 behaviour of the lifetime distributions.
Another set of statistical measurements that we made for the durations of the quasisteady states (QSS) gave a somewhat conflicting result. The QSSs are the metastable periods, during which the community structure appears to be stationary [7,8]. They are interrupted by active periods, during which the community structure changes after the emergence of a successful mutant and resulting extinction of other species. In order to identify the QSSs, we recorded the entropy of the system every sixteen generations.
The information-theoretical entropy is given by [19,20], where ρ I (t) = n I (t)/N tot (t). The system is considered in a QSS, as long as the magnitude of the entropy difference, |S(t) − S(t − 16)|/16, is below a certain threshold, which we took as 0.015. (For an extensive discussion on how to obtain this value, see reference [7].) The QSS-duration distribution for the uncorrelated model is not significantly affected by the change in the standard deviation of the M ij distribution, except the density of the longest lived communities around 10 7 generations ( figure 5(a)). Considering the large error bars at this range, the QSS distributions of the uncorrelated systems nearly overlap. However, the correlations in M seem to increase the lifetimes of the QSS communities ( figure 5(b)).
In order to get more information on how the fluctuations in the system are affected by the correlations in M, we also calculated the power spectral density (PSD) of the Shannon-Wiener diversity index. The Shannon-Wiener diversity index [21] is defined as , where S is the information-theoretical entropy defined in equation (3). In contrast to the QSS statistics (figure 5(a)) (but consistent with the case for the species-lifetime, figure 4), the PSDs of the Shannon-Wiener indices of the uncorrelated systems with different M ij distributions do not overlap (figure 6(a)) [22,23]. Although the PSDs are not monotonic in the standard deviation of M ij , they have a similar 1/f like shape, indicating that the general characteristics of the system do not vary significantly with the M ij distribution. However, the correlations in M lead to a more pronounced 1/f like behaviour in the PSDs by increasing the relative weight of high-frequency fluctuations as the correlations grow stronger ( figure 6(b)). Although the deviations from the uncorrelated case are neither gradual nor monotonic in correlation strength, the result suggests that the correlations in M may in fact change some characteristics of the system, supporting the result obtained from the QSS duration statistics. Nevertheless, the qualitative 1/f behaviour of the PSDs is not affected by these changes.
The results discussed above may seem rather complicated, but we nevertheless believe some clear effects can be discerned.
There are clear differences between the results for the correlated and uncorrelated Gaussian matrix elements. The correlations further increase the proportion of very long QSS (figure 5b), but curiously they also lead to an increase of intermediate-lifetime species (figure 4). These effects do not seem very easy to reconcile, but we speculate that they may be due to a decoupling of the species lifetimes and the QSS durations that has been observed in an uncorrelated predator-prey version of this model [9]. The increased activity in the intermediate-time regime results in increased intensity in the higher-frequency half of the PSD (figure 6b). It should be noted that strong correlation effects are only seen when n becomes a significant proportion of the genome length L (here, typically n ≥ 3). Such long-range correlations are of very limited interest for the biological systems, and they only become significant in the present study because we,  for computational reasons, are working with the unrealistically short genome length, L = 8. With n ≪ L, the effects of correlations appear to be quite minor.

Summary and Conclusions
In this paper we have considered the effects of introducing correlations between the elements of the interaction matrix M that determines the interspecies interactions in an individual-based coevolution model [7,8,9], in which individuals are represented by a genome in the form of a bitstring of length L. The correlations were introduced by replacing each element M 0 ij in an uncorrelated interaction matrix with the average over all the elements M 0 kl such that the Hamming distance between the concatenated bitstrings ij and kl is less than or equal to n, and then reweighting the resulting matrix element M ij such as to maintain the standard deviation of its probability density unchanged.
In section 3 we calculated numerically for different n the correlation functions between modified matrix elements M ij and M kl as functions of the Hamming distance between ij and kl, and we compared these with the theoretical results derived in Appendix A. For short-range correlations (n < L/2) we found good agreement between the numerical and theoretical results, but for larger n there were significant discrepancies. These discrepancies were explained in Appendix B as a result of the reweighting of the correlated matrix elements that was performed in order to keep their standard deviation approximately equal to that of the uncorrelated ones.
Next, in section 4, we performed long kinetic Monte Carlo simulations of our coevolution model, both with correlated and uncorrelated interaction matrices. As an indicator of the similarity of the evolution processes we calculated the probability densities of the lifetimes of individual species in the two models. While we found statistically significant differences between the lifetime distributions in the two models (see figure 4), in both cases the distributions stayed near a 1/t 2 power law over near seven decades in time. However, the distributions for the durations of QSSs (figure 5) and the PSDs for the diversity index ( figure 6) showed that the uncorrelated and correlated systems have some different characteristics. Nevertheless, the overall behaviour of the system does not seem to change drastically, since the 1/f behaviour in the PSDs remain the same. This indicates that the correlations in M affect the long-term behaviour of the system in minor, but rather complicated ways. The correlation effects are sufficiently mild that it is permissible to draw conclusions from simulations of uncorrelated models. In this average, all cross terms, M 0 ab M 0 cd where a = c and/or b = d, are equal to zero. However, if a = c and b = d, then the term (M 0 ab ) 2 is equal to σ 2 0 , which is the variance of the elements of the random interaction matrix. Therefore, the correlation function becomes where N L,n (r) is the number of common terms in the two sums in equation (A.2). The calculation of N L,n (r) is explained below.
sums in equation (A.2) as the following. What is the total number of bitstrings that are n bits or less away from both bitstrings A and B, when A and B differ from each other by r bits? Equivalently, we want to find the number of identical bitstrings we can create by making n or less modifications on A and B each.
To make it easier to visualize the solution, we arrange the two bitstrings A and B so that the matching bits constitute the left-hand part of each bitstring. We shall call these 2L − r matching bits the clean part, and the r non-matching bits the dirty part ( figure A2(b)). We emphasize that this arrangement is only for visualization purposes, and not a part of the solution.
The contributions to N L,n (r) can be grouped in three cases: a) n ≥ r ≥ 0: First, we make the dirty parts identical by making i modifications on either A or B and r − i modifications on the other. There are r i=0 r i = r i=0 r r−i = 2 r ways of making the dirty parts identical. Since n > r, we can also make k = n − max(r − i, i) = min(n − r + i, n − i) modifications on the clean part, unless k > 2L − r, in which case the number of modifications exceeds 2L − r, the size of the clean part. Therefore, the correct upper limit for k should be min(2L − r, n − r + i, n − i). Thus, the total number of different configurations that can be obtained on the clean part is . Combining these, we obtain the contribution to N L,n (r) for the 0 ≤ r ≤ n case: b) 2n ≥ r > n: When r > n, we have to make at least r − n modifications on one of the dirty parts. Therefore, the contribution to N L,n (r) for the 2n ≥ r > n case is Note that the upper limit of the first sum is not r anymore, since we are allowed to make only n < r modifications in total. c) r > 2n: Since the two hyperspheres in figure A2(a) cannot overlap when r > 2n, the contribution to N L,n (r) in this case is zero. A plot of the theoretical correlation functions for L = 13 with n =1, 2, 3, 5, and 7 is shown in figure A1.

Appendix B. Deviations from the Theory in Correlation Functions of Highly Correlated Matrices
The theoretical calculations for the correlation function agree with the numerical results for matrices with short-range correlations. However, the numerically obtained correlation functions for highly correlated matrices differ significantly from the theory, as seen in figure 2. We here demonstrate why these deviations occur by calculating the variance of the correlated matrix elements (which is equal to the value of the correlation function at r = 0): Although we intended to keep the variance of the correlated matrix elements constant through multiplication by Z n (see equation (2)), equation (B.6) indicates that the deviations from the intended theoretical value σ 2 0 = (M 0 ij ) 2 scale by Z n . Our numerical studies showed that for Z n L/2, the term (Z n −1) M 0 ij M 0 kl 1≤H(ij;kl)≤2n −Z n M 0 ij 2 can no longer be neglected, and deviations from the theoretical value become significant.
We Although this calculation shows the deviation only for C(0) = σ 2 corr , deviations for r > 0 should be similar due to the relatively smooth change in the overlap of the averaging hyperspheres. Since the overlaps of the averaging hyperspheres are similar for adjacent r's, the deviations also should be similar for the adjacent C(r)'s. For example, the deviations at C(0) and C(1) should be similar in magnitude and direction. This argument explains why the numerical correlation functions look like their theoretical counterparts translated along the y axis, rather than randomly scattered around the expected data points.