text:
# Statistical properties of neutral evolution

## I Introduction

The molecular clock hypothesis, proposed by Zuckerkandl and Pauling about 40 years ago Zuckerkandl & Pauling (1962), has been a cornerstone in the foundations of molecular evolution. In comparing the sequences of homologous proteins, Zuckerkandl & Pauling (1962) observed that the rate of amino acid substitutions per site per year in a given protein family is approximately constant, independently of the pair of species compared. In particular, the substitution rate appeared to depend only weakly on the number of individuals in the population and on ecological variables, a property which is extremely useful for reconstructing phylogenetic trees from the comparison of protein sequences Saitou & Nei (1987); Thompson et al. (1994); Grauer & Li (2000). The constancy of the substitution rate was shown some years later to be an immediate consequence of the neutral theory of molecular evolution, proposed independently by Kimura (1968) and by King & Jukes (1969).

The proposal of the neutral theory raised heated controversies, in part because it challenged the expectation that most variations in the genetic material are driven by positive natural selection. At the present time, it is however generally recognized that both neutral and adaptive substitutions play an important role in the evolution of protein sequences, and methods are being developed to distinguish between them. Moreover, the neutral paradigm has been extended to include slightly deleterious mutations, as proposed by Ohta (1973, 1976). Models of slightly deleterious mutations are however rather complicate from the point of view of population genetics, and will not be considered here.

Despite this debate, less attention has been drawn to the question whether Kimura’s neutral model adequately describes the statistical properties of neutral evolution of proteins. This question is an important one, since an accurate knowledge of such expected properties would allow to better distinguish between neutral and adaptive substitutions and to improve methods to reconstruct evolutionary trees from sequence alignments.

Kimura’s model assumes that the overwhelming majority of mutations are either effectively neutral (their effect on the fitness is much smaller than the inverse of the effective population size) or lethal, in the latter case purged by negative selection. We will call this assumption the neutral hypothesis. Moreover, the rate of occurrence of neutral mutations is regarded as constant throughout evolution, independent of the current sequence. This is the homogeneity hypothesis. As a consequence, the rate of neutral substitutions is predicted to be constant, and the number of neutral substitutions being fixed after $`k`$ generations of evolution is predicted to be a Poissonian variable of mean value $`\mu xk`$, where $`\mu x`$ is the neutral mutation rate. At variance with this, the first tests performed by Ohta & Kimura (1971) on homologous proteins showed that the variance of the substitution distribution appears to be larger than that expected for a Poissonian distribution. However, since deviations seemed at first to be small, the Poissonian statistics was still regarded as a valid first approximation. More accurate later studies showed that such deviations are in fact not small Langley & Fitch (1973); Gillespie (1998). This discrepancy between the model and empirical observations was taken by Gillespie as evidence against the neutral hypothesis, and he favored the hypothesis that most substitutions in protein sequences are fixed by positive selection Gillespie (1991). On the other hand, Takahata (1987) proposed a modification of the neutral model, the fluctuating neutral space model, that can account for the non-Poissonian statistics of substitutions, still preserving the simplicity of the neutral model as the simplest population genetics model. This is one of the few instances in which the debate, instead of concentrating on the neutral hypothesis, reconsidered the way in which Kimura had modeled neutral mutations.

Recently, we have introduced a model in which the neutrality of mutations in protein sequences is explicitly tested by means of a computational model of protein folding Bastolla et al. (1999, 2000a, 2002a, 2002b). The statistical properties of the model are rather robust: We obtained the same qualitative behavior by using two folding models as different as a well designed lattice model Abkevich et al. (1994) and a model which allows recognition of native protein folds Bastolla et al. (2000b, 2001). In particular, we found that the hypothesis of homogeneity should be rejected and that the neutral mutation rates fluctuate broadly during evolution. As a consequence, the variance of the substitution process is much larger than that expected for a Poissonian distribution, in reasonable agreement with the statistics observed for most protein families.

Our model belongs to the “structural approach” to molecular evolution, which has been made possible by the recent advances in the understanding of the dynamics of protein folding and the thermodynamics of biomolecules. In particular, it is now possible to assess approximately the thermodynamic stability of biomolecules by computational methods. We have undertaken the challenge of exploring the applications of these new methods to the field of molecular evolution, where they can complement more traditional methods based on population genetics and on sequence comparison. The structural approach has been pioneered by Schuster and co-workers with a series of studies of neutral networks of RNA secondary structures Schuster et al. (1994); Huynen et al. (1996); Fontana & Schuster (1998) and it has been applied to proteins by several other groups Shakhnovich et al. (1996); Bornberg-Bauer (1997); Babajide et al. (1997); Govindarajan & Goldstein (1997, 1998); Bussemaker et al. (1997); Tiana et al. (1998); Mirny & Shakhnovich (1999); Bornberg-Bauer & Chan (1999); Dokholyan & Shakhnovich (2001); Xia & Levitt (2002).

In this paper, we study the statistical properties of neutral evolution of proteins as they are obtained from the Structurally Constrained Neutral (SCN) model that we have recently introduced Bastolla et al. (2002a). Our goal is to explore the consequences of these statistical properties on molecular evolution. We start the paper by reviewing Kimura’s model of neutral evolution, making explicit the homogeneity hypothesis and the arguments that seem to support its validity. In the next section, we review the SCN model, on which our results are based. Conservation of the protein fold during evolution is discussed in the following section, where we argue that fold conservation can be regarded as a consequence of the requirement of thermodynamic stability. We then investigate the origin of the strong fluctuations in the fraction of neutral neighbors that are found within the SCN model and that produce the non-Poissonian statistics of the substitution process. We define the neutral connectivities at each position in the protein and study their pairwise correlations. We find that even connectivities at positions far apart in the structure are strongly correlated. Remarkably, the only exceptions to this behavior are represented by structurally conserved positions, which are much less correlated with other positions and evolve in a more regular fashion. These spatial correlations are the counterpart of the temporal correlations that we described in an earlier work Bastolla et al. (2002b). In the following sections, we present some consequences of the broad distribution of neutral rates. We show that: (i) The neutral mutation rate is not constant as a function of the time interval considered, but decreases monotonously as larger time intervals are considered. (ii) The normalized variance of the substitution process increases with the time interval, and is in reasonable agreement with the data on the dispersion index for most proteins. (iii) Fluctuations in the neutral connectivities can influence the generation time effect. Our simulations are then used to test methods applied to estimate the number of substitutions in the divergence of two sequences. Finally, we discuss the impact of our findings on methods used to detect positive selection (as for instance the one proposed by McDonald & Kreitman (1991)), which are usually based on neutral models with constant neutral mutation rate. After summarizing our results in the Conclusions, the section ‘Materials and Methods’ reports the details of some technical points.

## II Kimura’s neutral model revisited

As discussed in the Introduction, Kimura’s neutral model is based on two assumptions. The neutral hypothesis is equivalent to the assumption that the most common mutations in protein sequences are either disruptive and eliminated by negative selection, or neutral, i.e. they leave the protein active and their effect on fitness is much smaller than the inverse of the effective population size. This mutational spectrum implies that protein sequences evolve on a neutral network, i.e. a set of sequences where the protein is active and which can be connected through point mutations. Fixation of slightly deleterious mutations, as well as advantageous mutations, are not included in the model. This is of course an important limitation of neutral models.

The second assumption is that the rate of appearance of neutral mutations is constant throughout evolution (homogeneity hypothesis). In a 1977 paper, Kimura commented that rate constancy may not hold exactly Kimura (1977), but he did not develop the consequences of the violation of this hypothesis any further.

The neutral model predicts that the rate of fixation of neutral mutations in an evolving population is equal to the rate of their appearance, independently of the population size, because the number of appearing mutations is proportional to the population size and the probability of their fixation is inversely proportional to it. Let us consider a given lineage evolving for a time interval $`t`$. The number of mutations appearing in this time interval is expected to obey a Poissonian statistics with mean value $`\mu t`$, where $`\mu `$ is the mutation rate. We call neutral connectivity the probability that one of such mutations is neutral, which is proportional to the connectivity of the neutral network, and denote it by $`x`$. The value of $`x`$ depends on the particular protein family chosen but it is assumed to be constant throughout evolution. As a consequence, the probability $`\mathrm{P}(S_t=n)`$ of the occurrence of $`n`$ neutral mutations within a time interval $`t`$ is

$$\mathrm{P}(S_t=n)=\underset{m=n}{\overset{\mathrm{}}{}}\mathrm{e}^{\mu t}\frac{(\mu t)^m}{m!}\left(\genfrac{}{}{0pt}{}{m}{n}\right)x^n(1x)^{mn},$$
(1)
where the factor $`\left(\genfrac{}{}{0pt}{}{m}{n}\right)x^m(1x)^{mn}`$ gives the probability that $`n`$ out of $`m`$ mutations are neutral. The summation can be performed exactly, yielding

$$\mathrm{P}(S_t=n)=\mathrm{e}^{\mu xt}\frac{(\mu xt)^n}{n!},$$
(2)
so that the number of neutral mutations also obeys a Poissonian statistics, with mean value $`\mu xt`$ and rate $`\mu x`$.

In principle the neutral connectivity $`x(𝐀)`$ depends on the amino acid sequence $`𝐀`$ considered, but Kimura’s model assumes that $`x(𝐀)`$ is the same for all protein sequences belonging to the same structural and functional family. This assumption can be justified with the following argument: Let $`x_i(𝐀)`$ be the fraction of possible mutations of sequence $`𝐀`$ at position $`i`$ which are neutral. Then the overall fraction of neutral mutations is just the average of this quantity over the $`N`$ positions in the amino acid chain,

$$x(𝐀)=\frac{1}{N}\underset{i}{}x_i(𝐀).$$
(3)
If the neutral connectivities at different positions are uncorrelated, or if only positions which are close in the native structure are correlated, then the variance of $`x(𝐀)`$ is proportional to $`1/N`$, which is very small for long protein chains, so that the approximation of constant $`x`$ is expected to be reasonably good. As we shall see in the following, this is not the case.

## III The Structurally Constrained Neutral Model

In the SCN model we assume the validity of the neutral hypothesis, but we do not make any assumption regarding the neutral mutation rate. Instead, we estimate explicitly the effect of a mutation on protein stability using an effective model of protein folding Bastolla et al. (2000b) which provides us with a genotype to phenotype mapping. In this respect, the rate of occurrence of neutral mutations is an outcome of the model. We show that this rate displays very broad fluctuations throughout evolution.

The SCN model addresses protein evolution at the level of a single sequence, and it does not take into account population dynamics. This simplification is justified by the fact that, within Kimura’s model, the substitution rate does not depend on the population size. Nevertheless, population size might influence the evolutionary process if the rate of neutral mutations shows broad fluctuations, as observed here. The explicit inclusion of population genetics into the model will be needed to investigate this interesting possibility.

### III.1 Estimating protein stability

In our model of protein folding, we evaluate effective conformational energies (temperature and pH dependent) using the effective energy function described by Bastolla et al. (2000b, 2001), which provides good thermodynamic properties for protein structures in the Protein Data Bank (PDB). For each sequence typically several millions alternative conformations of the same length are generated by gapless threading. The native conformation is identified as the lowest energy conformation. The energy function has several non-trivial properties: (i) It assigns lowest energy to the experimentally known native structures of basically all single chain protein sequences; (ii) It provides the native structure with a well correlated energy landscape (see below). (iii) For chains belonging to oligomeric proteins, for which interchain interactions are neglected, the experimentally known native state shows a deficit of stabilizing energy, well correlated with the amount of neglected free energy. (iv) The effective energy function that we use is able to estimate the folding free energy of proteins with two states kinetics reasonably well (U. Bastolla, unpublished).

There are two necessary conditions for thermodynamic stability (a) a low folding free energy, i.e. the unfolded state is less stable than the native state; (b) stability with respect to alternative misfolded states, i.e. the protein has a well correlated energy landscape, a property which turns out to be crucial in all simple models of protein folding. To enforce the first kind of stability, we estimate the folding free energy, which, in the simplest approximation, is just a linear function of the effective energy of the native state. The second condition is estimated through two computational parameters: The $`Z`$-score Bowie et al (1991); Goldstein et al. (1992), which measures the difference between the native energy and the average energy in units of standard deviation of the energy, and the normalized energy gap $`\alpha `$ Gutin et al. (1995); Bastolla et al. (1999), measuring the minimal value of the energy gap between an alternative conformation and the target one divided by their structural distance. A positive and large value of the $`\alpha `$ parameter ensures both that the target conformation has lowest energy and that the energy landscape is well correlated, in the sense that conformations very different from the native one have very high energy.

We impose stability conservation by requiring that the lowest energy state of the model coincides with the experimentally known native state and allowing its stability parameters, previously described, to be off not more than 1.5% with respect to the values of the corresponding PDB sequence.

### III.2 Exploring the neutral network

Simulations of protein evolution are performed starting from a protein sequence in the PDB. Evolution is constrained to viable sequences, which are sequences where the native structure is thermodynamically stable. A neutral network is a set of viable sequences which can be connected to the starting sequence through point mutations passing on other viable sequences. Thus sequences on a neutral network share the same protein fold and are evolutionarily connected.

For every amino acid sequence $`𝐀`$ in the neutral network we measure the fraction of neutral neighbors $`x(𝐀)(0,1]`$, which is the fraction of possible point mutations that are viable. Since this number defines the connectivity of the neutral network at point $`𝐀`$, we shall also call $`x(𝐀)`$ the neutral connectivity of sequence $`𝐀`$. In the framework of our model, this is the only property of a sequence which influences its evolution.

Subsequently visited sequences belonging to the neutral network constitute an evolutionary trajectory $`\{𝐀_1,𝐀_2,\mathrm{}\}`$. To each trajectory is associated the list of the corresponding neutral connectivities $`\{x(𝐀_1),x(𝐀_2),\mathrm{}\}`$.

### III.3 Substitution process

An amino acid substitution is controlled by two independent events: A random mutation, described as a Poissonian process as in the previous section, and an acceptance process which consists in testing whether the sequence is viable. The acceptance probability for a mutation taking place when the protein is in sequence $`𝐀`$ is given by $`x(𝐀)`$. As a result of the fluctuations in the neutral connectivities, the resulting substitution process is no longer Poissonian. For a given sequence of neutral connectivities $`\{x(𝐀_1),x(𝐀_2),\mathrm{}\}`$ we can compute the probability that the number $`S_t`$ of accepted mutations in a time interval $`t`$ is equal to $`n`$. This is just the product of the probability (Poissonian) that $`k`$ mutations take place in the time interval $`t`$ times the probability that $`n`$ of these are accepted,

$$\mathrm{P}(S_t=n)=\underset{k=n}{\overset{\mathrm{}}{}}\mathrm{e}^{\mu t}\frac{(\mu t)^k}{k!}P_{\mathrm{acc}}(n|k),$$
(4)
where the acceptance probability of $`n`$ mutations out of $`k`$ is given by

$$P_{\mathrm{acc}}(n|k)=\left(\underset{i=1}{\overset{n}{}}x(𝐀_i)\right)\underset{\{m_j\}}{}\underset{j=1}{\overset{n+1}{}}\left[1x(𝐀_j)\right]^{m_j}.$$
(5)
Here, the $`\{m_j\}`$ are all integer numbers between zero and $`kn`$ satisfying $`_{j=1}^{n+1}m_j=kn`$. In other words, the probability that a mutation is accepted is $`x(𝐀_1)`$ as long as the protein sequence is $`𝐀_1`$, $`x(𝐀_2)`$ as long as the sequence is $`𝐀_2`$, and so on.

Two kinds of random variables have to be distinguished. The first kind of variables represent the mutation and acceptance process. The average over these variables is indicated by angular brackets $`<>`$, and it is still a random variable dependent on the realization of the evolutionary trajectory. The average over evolutionary trajectories is indicated by an overline $`\overline{}`$. In the biological metaphor, the average over mutation and acceptance correspond to population averages and the average over different evolutionary trajectories correspond to averages over independent populations.

If all sequences have the same fraction of neutral neighbors $`x(𝐀)=x`$, the number of substitutions in a branch of length $`t`$ is Poissonian with mean $`\mu tx`$ and the substitution rate is equal to $`\mu x`$, as in Kimura’s model. If the variance of the neutral connectivity is not zero, the substitution distribution is more complicated and has to be computed numerically using the simulated evolutionary trajectories (see the section ‘Materials and Methods’).

## IV Results

An important property that we have shown to hold for the SCN model is that its statistical properties are robust: Their qualitative behavior does neither depend on the thresholds used to select viable sequences (and does even not change if lattice models are used instead of effective models of protein folding with real protein structures), nor it depends on the protein structures considered, although the studied proteins cover a broad spectrum of different biological activities. The protein structure only determines which positions play a structural role and are therefore more conserved than average positions, but the properties of structurally conserved residues, like the fact that their evolutionary rate is less dependent on the context of the sequence (see below), are general features for all structures. Although some properties that we determined may depend on the length of the protein, our data is insufficient to quantify such an effect.

We studied in this work eight protein folds, considering for the rubredoxin fold two different structures, one of a mesophilic and one of a thermophilic species. The nine proteins studied are: myoglobin (PDB code 1a6g), cytochrome c (PDB code 451c), lysozyme (PDB code 3lzt), ribonuclease (PDB code 7rsa), rubredoxin (mesophilic species: PDB code 1iro; thermophilic species: PDB code 1brf\_A), ubiquitin (PDB code 1u9a\_A), the TIM barrel (PDB code 7tim\_A), and kinesin (PDB code 1bg2). In what follows we describe general properties of the model which are qualitatively the same for all proteins studied. Unless otherwise stated, the results that we present refer to the myoglobin fold.

Although the model of protein stability that we use is only an approximate one, it is noteworthy that our most important qualitative results, the broad distribution of neutral connectivities and the overdispersion which is caused by correlations along neutral trajectories (see below), reproduce those of a former work in which we used a model of protein stability in some sense complementary to the present one Bastolla et al. (1999, 2000b). In that work the statistical mechanics of protein sequences was simulated using a lattice model to represent the ensemble of conformations. Such a treatment provides only a coarse-grained representation of protein structure (e.g., secondary structures are not well described). It does, however, allow for the enforcement of a rigorous criterion of stability, that has been then tested through extensive Monte Carlo simulations. The fact that qualitatively similar results are recovered using two complementary approximations to the protein folding problem makes us confident that our model captures some “universal” properties of protein evolution.

### IV.1 SCN and structural conservation

It is commonly observed that protein structures are much more conserved than amino acid sequences Holm & Sander (1996); Rost (1997), but there is no real evolutionary justification to impose such a rule as we do in the SCN model. In fact, the target of natural selection is protein function rather than protein structure, and the relationship between these two properties is not a simple one. It is well known that proteins with the same fold can perform quite different functions, and in some cases there is evidence that those proteins are evolutionarily related. In fact, proteins often acquire novel functions by gene duplication though maintaining a very similar structure. Moreover, the case of structurally unrelated proteins which perform the same function is not rare, a possible result of convergent evolution.

A more realistic version of the SCN model should impose conservation of protein thermodynamic stability, irrespective of the native structure. This, in our opinion, should be a requirement of any model of molecular evolution. We expect that such an improved version of the SCN model would show that the conservation of the native structure is a consequence of the stability requirement. The reason for this expectation is based on our simulations of the present SCN model and of a previous lattice version of it. In the SCN simulations we record the smallest value of the normalized energy gap $`\alpha `$ with the target conformation, over all neighbors of a sequence belonging to the neutral network. For seven out of the nine proteins that we studied, the smallest value of $`\alpha `$ remained positive for all of the order of $`10^5`$ examined sequences. Since a positive $`\alpha `$ means that the target conformation has the lowest energy, imposing that the target conformation is stable is effectively equivalent to imposing that the native (lowest energy) conformation is stable. Thus stability requirements alone suffice to guarantee the conservation of these protein structures. The two exceptions that we found are cytochrome c and the mesophilic version of rubredoxin, the two smallest proteins that we studied. For these proteins some sequences in the neutral network have neighbors where a structure different from the target state has lowest energy. However, this is not enough to ensure thermodynamic stability of this structure, since we still have to impose that all structures unrelated to it have much higher energy, while the energy of the target structure is very low. We shall investigate this subject in more detail in the future.

From the above considerations it turns out that imposing a well correlated energy landscape through a condition on the normalized energy gap makes it difficult to change the native structure maintaining the thermodynamic stability of the new structure. This result is very different from the one obtained in the study of neutral networks of RNA. In this case, Schuster and co-workers have shown through a computational study that the neutral networks of two different RNA secondary structures can be separated by just one point mutation Schuster et al. (1994). We think that this difference is not only due to the difference between protein folding and RNA folding, but also to the different model of thermodynamic stability applied. Schuster and co-workers require that the target structure has the lowest energy, while we impose additionally that it has a large folding free energy and a large normalized energy gap.

### IV.2 Spatial correlations

As we have seen, if the neutral connectivities $`x_i`$ of different positions are only weakly correlated, the overall neutral connectivity $`x`$ has very small fluctuations and the neutral model of Kimura is valid. However, we observed in a previous work that the connectivity distribution is quite broad, see Fig. 1, which implies that the statistics of the substitution process is very different from the Poissonian statistics. To investigate the origin of this important property, we examined the correlation matrix $`𝒞_{ij}`$, whose elements are the correlation coefficients of the neutral connectivities at positions $`i`$ and $`j`$ (see the section ‘Materials and Methods’).

Two interesting observations emerge. First, the correlations coefficients $`𝒞_{ij}`$ are positive and large (larger than 0.2) for most positions, irrespective of whether they are in contact or not in the native structure. Fig. 2 shows a comparison between the correlation matrix $`𝒞_{ij}`$ (upper right part of the central plot in the figure, see the section ‘Materials and Methods’) and the contact matrix of the native structure $`C_{ij}`$ (lower left part of the central plot of the figure). Second, there is a strong relation between correlation and conservation. Positions which are more conserved, as indicated by the neutral network average ($`\overline{x_i}`$) of their fraction of neutral neighbors, show only a weak correlation with other positions. We plot as black dots in the upper right part of Fig. 2 pairs of positions which are only weakly correlated, $`C_{ij}0.3`$. These dots arrange in horizontal and vertical lines at the conserved positions. The pattern is more clearly shown in Fig. 3, representing a scatter plot of the variability index $`\overline{x_i}`$ versus the correlation coefficient $`𝒞_i`$ between the neutral connectivity at position $`i`$ and the overall neutral connectivity (see the section ‘Materials and Methods’). It can be seen that the two quantities are strongly correlated.

Ptitsyn & Ting (1999) showed that there are two groups of conserved residues in the globin family. In the myoglobin sequence that we considered the first group is the heme-binding site, formed by the residues Leu29 (B10), Leu32 (B13), Phe33 (B14), Pro37 (C2), Phe43 (CD1), Phe46 (CD4), Leu61 (E4), His64 (E7), Val68 (E11), Leu89 (F4), His93 (F8), Ile99 (FG5), Leu104 (G5), and His142 (H19), where the number in parenthesis indicates the standard notation of Perutz et al. (1965) that specifies the position within an helix. The second group, whose conservation is structural and not functional Ptitsyn & Ting (1999), is formed by residues Val10 (A8), Trp14 (A12), Ile111 (G12), Leu115 (G16), Met131 (H8), and Leu13 (H12). In agreement with our argument, the residues in the second group are among those with the lowest values of $`𝒞_i`$, they are ranked 9th, 6th, 22th, 24th, 1st, and 14th, respectively. On the contrary, residues in the functionally conserved group are characterized by average values of $`𝒞_i`$, ranging in rank from from 12 to 132, with an average rank of 61. Even more interestingly, residue Met131, the one with the lowest $`𝒞_i`$, is anti-correlated (although weakly) with some residues in the functionally conserved group. This is the only case of significant anti-correlation. Leu69 and Leu76 are ranked 5th and 2nd, respectively. They have nearly zero correlation among themselves and with residues in the two groups. Interestingly, they are located between the two groups in the native structure. Residues His119 and Phe123 are ranked 3rd and 4th. They are in the loop between helices G and H and they are probably important to stabilize the second group of conserved residues. A similar observation holds for residues Val13 and Ala134, ranked 7th and 8th, respectively, since they are neighbors of residues Trp14 and Leu115. Taken together this data suggests that the analysis of cross-correlations complements Ptitsyn’s conservation analysis of structurally important residues Shakhnovich et al. (1996); Ptitsyn & Ting (1999).

### IV.3 Temporal correlations

To investigate the temporal correlations, we measured the auto-correlation function

$$C(x(𝐀_i),x(𝐀_{i+n}))=\frac{1}{m}\underset{k=1}{\overset{m}{}}\frac{x(𝐀_k)x(𝐀_{k+n})\overline{x}^2}{\sigma _x^2},$$
(6)
where an average is taken over all pairs of sequences connected through an evolutionary trajectory of exactly $`n`$ substitutions. The corresponding plot can be seen in Fig. 4. The auto-correlation function starts with the value one at $`t=0`$ and then decreases as subsequent substitutions make the sequences less and less correlated.

The characteristic time after which the correlation function reduces of a factor $`1/\mathrm{e}`$ can be estimated through an exponential fit of the first part of the correlation function. For all proteins examined, the characteristic length is two or three substitutions. Despite this fast decay, temporal correlations are responsible of many interesting properties of the model, from the large deviation of the distribution of the number of substitutions with respect to a Poissonian one to large variations of the neutral substitution rates in different evolutionary trajectories (see also Bastolla et al. (2002b)).

### IV.4 Origin of the correlations

The phenomenon that originates the spatial and temporal correlations that we have observed in the SCN model can be described as follows. We impose in our model that some stability indices are above some predefined thresholds in order to decide that a sequence is viable. Thus, a sequence where all of the stability parameters are very high will be much more tolerant to mutations than an average sequence: In this case, every positions will have a high neutral connectivity, since most amino acid changes at the position will maintain the protein viable. The only exceptions will be positions which perform a key structural role, which, as we have seen, are characterized by a much smaller neutral connectivity and are much less correlated with the neutral connectivity of the whole protein. The temporal correlations can be explained exactly in the same way: Since the thermodynamic properties change smoothly with changes in the sequences, they will remain very correlated after few changes in the sequences, so that also the neutral connectivities which derive from these stability parameters will show strong correlations.

To support our argument, we show results from a simulation where just one thermodynamic parameter has been considered to select sequences belonging to the neutral network. In the additional simulations presented in Fig. 5(a), only the native energy per residue has been used as a selection parameter. Accordingly, there is a very strong relationship between this thermodynamic parameter and the neutral connectivity $`x(𝐀)`$. In Fig. 5(b), both the normalized energy gap and the native energy have been used as selection parameters. Thus the relationship between thermodynamic parameters and neutral connectivity is less strong, but still it is clearly present.

### IV.5 Multiple substitutions

According to the molecular clock hypothesis, the number of substitutions in a time interval $`t`$ is proportional to the duration of the time interval, with the addition of random fluctuations, whose standard deviation is small with respect to the mean when $`t`$ is large. However, in the comparison of two protein sequences, only the number of sites occupied by different amino acids, $`N_\mathrm{d}`$, is accessible to measurement. From this quantity, the number of substitutions has to be calculated assuming some model of protein evolution to correct for possible occurrence of multiple substitutions at the same site.

The simplest such method was used already by Zuckerkandl & Pauling (1962). Assuming that all sites in a protein evolve at the same rate, and that the number of substitutions for each site follows a Poissonian distribution of average value $`K`$, the probability that one site did not change is $`\mathrm{exp}(K)`$. This probability can be estimated by the fraction of unchanged sites, $`1N_\mathrm{d}/N`$, so that the expected number of substitutions per site $`K`$ can be estimated as

$$K\mathrm{log}\left(1\frac{N_\mathrm{d}}{N}\right).$$
(7)
The assumption that the substitution rate is the same at all sites is stronger than the homogeneity assumption stating that the substitution rate is the same for all sequences. It is well known that functionally important residues change very rarely during evolution, often as the result of a change in the environment or a switch in the protein function, while other residues are much more free to mutate. Even within our model, which does not take into account functional constraints explicitly, some residues are much more difficult to mutate than others, since they are structurally important and belong for instance to the so-called folding nucleus of the protein. In Fig. 6 we show the distribution of the average fraction of substitutions at position $`i`$ which are neutral, for all 151 positions in the myoglobin fold.

As noticed by Nei & Kumar (2000) and by Kimura (1983), the Poissonian correction Eq. (7) remains valid also in case of varying evolutionary rates for small values of the fraction of changed positions, $`N_\mathrm{d}/N`$. Otherwise, one has to consider more sophisticated corrections. The most common method which takes into account variation of rate across sites consists in assuming that the rates are distributed according to a gamma distribution, which is the product of an exponential times a power law distribution Uzzell & Corbin (1971). It is then possible to compute analytically the probability that a site has not changed and, equating it to the observed frequency, to obtain the average number of substitutions per site as

$$K=a\left[\left(1\frac{N_\mathrm{d}}{N}\right)^{1/a}1\right],$$
(8)
where $`a`$ is the shape parameter of the gamma distribution Nei & Kumar (2000). The above results, Eq. (7), is recovered in the limit $`a\mathrm{}`$, as in this limit the gamma distribution becomes a delta distribution with vanishing variance. However, in realistic situations, $`a`$ tends to be small (typically smaller than four), so that the distribution of rates across different positions is broad and the estimated number of substitutions is much larger than the one obtained using Eq. (7).

In typical studies, the parameter $`a`$ is first estimated, for instance by maximum parsimony, and then used to obtain the number of substitutions through Eq. (8). Our simulations yield the number of substitutions for a given fraction of mutated positions, $`N_\mathrm{d}/N`$, so that, using our data, we can directly test the validity of Eq. (8).

Results of such an analysis are shown in Fig. 7. From Fig. 7(a) it is clear that two regimes are encountered. The first regime (large sequence similarity characterized by small $`N_\mathrm{d}/N`$) will be called the transient regime. In this case the average sequence similarity decreases on the average with time, measured as the number of substitutions, and the standard deviation is small. We can thus use the measured similarity to estimate the number of substitutions. This regime is illustrated in Fig. 7(b), where it is compared with various kinds of corrections for multiple substitutions. As it can been seen, the gamma corrections with small $`a`$ fits the data better than the Poissonian correction ($`a=\mathrm{}`$), although systematic deviations are appreciable. Only for lysozyme we find an optimal value of $`a`$ close to eight, whereas for all other proteins the optimal value, although not very precisely determined, lies in the range from one to five. For low sequence similarities, it is not allowed to neglect the probability $`p1/20`$ that a further mutation of a mutated site restores the amino acid initially present in the sequence. Taking this into account, Eq. (7) has to be modified to

$$K=\mathrm{ln}\left(1ppK\right)\mathrm{ln}\left(1\frac{N_\mathrm{d}}{N}p\right),$$
(9)
which, for small $`p`$, yields

$$K=(1p)\mathrm{ln}\left(1\frac{N_\mathrm{d}}{N(1p)}\right).$$
(10)
Similarly, Eq. (8) has to be modified to

$$\begin{array}{cc}\hfill \left(1p\frac{N_\mathrm{d}}{N}\right)=& \left(1+\frac{K}{a}\right)^a\hfill \\ & \times \left(1pp\frac{aK}{a+K}\right),\hfill \end{array}$$
(11)
which for small $`p`$ can be approximatly solved to yield

$$Ka(1p)\left[\left(1\frac{N_\mathrm{d}}{N(1p)}\right)^{1/a}1\right].$$
(12)
The previous formulas are recovered for $`p=0`$, while Eq. (10) is recovered as the infinite $`a`$ limit of the above equations. Notice that, for non-zero $`p`$, the expected number of substitutions now diverges when the sequence similarity $`1N_d/N`$ tends to $`p`$, and can not be computed anymore for smaller sequence similarity.

The second regime is called the stationary regime (right part of Fig. 7(a)). In this case, sequence similarity has reached a stationary distribution with respect to the number of substitutions, and it does not give us any information about the number of substitutions occurred in the evolutionary history. Clearly, in this case we cannot anymore identify the probability that a site has not changed, $`\mathrm{exp}(K)`$, with the frequency of identical sites, $`1N_\mathrm{d}/N`$, neglecting the fluctuations of the latter quantity, and consequently the above formulas loose their validity. The value of sequence similarity below which the transition to the stationary regime takes place depends on the sequence length $`N`$ and, more weakly, on the number of substitutions, $`KN`$, and is such that the stationary probability to observe a sequence similarity equal to $`N_\mathrm{d}/N`$ is roughly equal to $`1/KN`$. In Fig. 8 we plot the threshold of sequence similarity below which the standard deviation in the number of substitutions becomes of the same order of the average number of substitutions, and hence estimating the number of substitutions looses all meaning. Notice that, if some positions are conserved because of functional or other constraints, this threshold of similarity will be even higher. The data has been obtained from our simulations with $`KN=20000`$ substitutions for all nine proteins. The plot clearly shows the expected dependence of the threshold on the inverse square root of sequence length.

### IV.6 Time dependence of the substitution rate

In our analysis, two kinds of averages must be distinguished. We indicate by angular brackets $`<>`$ the average over the mutation and acceptance process for a given realization of the evolutionary trajectory, and by an overline $`\overline{}`$ the average over evolutionary trajectories. We determine the mean number of substitutions $`\overline{<S_t>}`$ within a time interval $`t`$, and show this quantity in Fig. 9. We also show in the same figure results based on annealed trajectories, obtained by extracting at random the values of $`x`$ at each substitution event according to the observed distribution $`P(x)`$ (cf. Fig. 1). In this case, the different $`x`$ along the trajectories are independent variables. Note that the annealed trajectories ‘interpolate’ between the Poissonian case (all $`x`$ are equal) and the correlated trajectories obtained through the SCN model: (i) The comparison between the annealed trajectories and the simple Poissonian case allows us to identify the effect of the broad distribution of the number of neutral neighbors, whereas (ii) the comparison between the actual and the annealed trajectories allows us to identify the effect of correlations.

As a result of the broad distribution of connectivities, the substitution rate $`\overline{<S_t>}/t`$ is not constant, as in the standard model of neutral evolution, but decreases as the time interval $`t`$ gets longer, as shown in Fig. 9 for both correlated and annealed trajectories. This can be qualitatively understood by the following consideration: In the annealed case, the time spans $`\tau `$ between subsequent substitutions are independent variables distributed with the density $`\overline{D(\tau )}=_0^1P(x)(\mu x)^1\mathrm{exp}(\mu x\tau )dx`$, whose average value is $`\overline{\tau }=_0^1P(x)(\mu x)^1dx`$. Thus, the average substitution rate $`\overline{<S_t>}/t`$ is not constant in time as in the Poissonian case. Initially, $`S_t`$ is a Poissonian variable with average rate $`\mu \overline{x}`$. At large time, however, the rate converges to the smaller value $`\overline{<S_t>}/t1/\overline{\tau }`$ (i.e., $`\overline{<S_t>}/(\mu t)(\overline{1/x})^1`$), since the process spends more and more time in sequences with small $`x(𝐀)`$.

### IV.7 Variances

The variance of the substitution process has been intensively studied in the first tests of the neutral theory. In those tests it was expected that the mean and the variance of the number of substitutions should be equal under neutral evolution, as the latter was expected to follow Poissonian statistics. The ratio between the variance and the mean, $`R=\mathrm{V}(S_t)/\mathrm{E}(S_t)`$, where $`\mathrm{V}`$ indicates variance and $`\mathrm{E}`$ indicates expectation value, is called dispersion index. Deviations of this quantity from unity indicate deviations from Poissonian statistics. In the framework of the SCN model, the broad distribution of neutral connectivities causes the substitution process to be overdispersed (its variance is larger than the mean), as we have shown previously Bastolla et al. (2002b). Here we investigate the relative contributions of the broad distribution and of the time correlations to the variance of the substitution process.

The variance of the substitution process can be decomposed in two components, one calculated for a fixed evolutionary trajectory (fixed population), the other taking into account the variance of different evolutionary trajectories,

$$\begin{array}{cc}\hfill \mathrm{V}(S_t)& =\mathrm{V}_\mu (S_t)+\mathrm{V}_x(S_t)\hfill \\ & =\left(\overline{<S_t^2>}\overline{<S_t>^2}\right)+\left(\overline{<S_t>^2}\overline{<S_t>}^2\right).\hfill \end{array}$$
(13)
The first term, $`\mathrm{V}_\mu `$, is the variance of the mutation and acceptance process for a fixed trajectory, averaged over evolutionary trajectories. The second term, $`\mathrm{V}_x`$, is the variance of the substitution rate with respect to different evolutionary trajectories. This term, which is not present in the standard neutral model, is responsible of the fact that the variance of the number of substitutions is typically larger than its mean value, in contrast with a Poissonian process.

We shall denote by $`R_\mu `$ the normalized mutational variance, $`R_\mu (S_t)=\mathrm{V}_\mu (S_t)/\mathrm{E}(S_t)`$, and by $`R_x`$ the corresponding normalized trajectory variance, $`R_x(S_t)=\mathrm{V}_x(S_t)/\mathrm{E}(S_t)`$. Thus, in the standard neutral model, one expects that $`R_\mu =1`$ and $`R_x=0`$, independent of the time interval $`t`$.

The results of the substitution process based on simulated trajectories are shown in Fig. 10 as full symbols. We also show as open symbols results based on annealed trajectories, which have the same distribution of neutral connectivities but lack the correlations. We start describing the latter results. As discussed in the previous section, for very short time interval the substitution distribution is practically Poissonian and one has $`R_\mu 1`$ and $`R_x0`$. As the time interval gets longer, the two variances grow only moderately: The mutation variance $`R_\mu `$ apparently reaches very soon a stationary value which is just slightly larger than one for all of the proteins that we studied, while the trajectory variance $`R_x`$ has not yet attained a stationary value when the number of mutation events is of the same order as the protein length, but still remains very small, of the order of the ratio between the variance and the average value of the neutral connectivities $`x`$.

Comparing the two ensembles of trajectories, we note that the presence of correlations has only a weak effect on the normalized mutation variance $`R_\mu `$, which remains very close to the value met for annealed trajectories. However, the normalized trajectory variance $`R_x`$ increases considerably in response to the correlations, as $`R_x1`$ for $`\mu t=200`$. Surprisingly, $`R_x`$ even grows with time, although more and more sequences are used to compute the mutational averages and one could expect that such averages approach typical values. Thus, surprisingly the large fluctuations between different trajectories caused by the strong correlations result in a peculiar phenomenon: Even averaging over an arbitrary long trajectory does not give values representative of typical trajectories.

The total dispersion index $`R=R_x+R_\mu `$, calculated for the trajectories obtained through the SCN model, becomes larger than $`2`$ already for $`\mu t100`$, in qualitative agreement with Gillespie’s results for a large set of proteins, most of which yielded dispersion indices in the range $`1<R5`$ Gillespie (1991). Hence, these large dispersion indices observed in protein evolution may be to a sizeable extent due to the correlations present in neutral trajectories in sequence space.

Notice however that, although the dispersion index is larger than one, still the main property of a Poissonian molecular clock is preserved, since, when the number of substitutions is large, the average number of substitutions $`\overline{<S_t>}`$ is much larger than the standard deviation $`\sqrt{\mathrm{V}(S_t)}`$. Thus fluctuations of the number of neutral substitutions do not hide completely the pattern arising from the evolutionary relationships.

### IV.8 Lineage effects and the generation time effect

The results presented in the previous sections imply the values of $`x`$ along one given evolutionary trajectory are very similar for several steps, because of the strong temporal correlations we identify. Hence, we will observe trajectories characterized by large substitution rates and trajectories characterized by small substitution rates, as a result of the broad and correlated fluctuations in the neutral connectivities. As an example, we show in Fig. 11 the substitution rates obtained from three different evolutionary trajectories. Note that the variance of the neutral rates $`<S_t>/t`$ across different evolutionary trajectories decreases with the time interval $`t`$, but rather slowly, so that even after a sizeable time interval the substitution rates of the three trajectories are still significantly different.

Since different trajectories can be interpreted as different species, different substitution rates in different trajectories can simulate lineage effects, i.e. variations of the substitution rate among different taxonomic groups. One such effect is the generation time effect: Since the natural time unit for measuring substitution events is the generation, species with longer generation time are expected to show a slower substitution rate. This prediction has been verified comparing for instance substitution rates in rodents and primates Britten (1986); Li et al. (1987); Grauer & Li (2000). Interestingly, the effect is significantly larger for synonymous substitutions (DNA changes which still code for the same amino acid) than for non-synonymous ones. Non-synonymous substitutions are superimposed with the large and correlated fluctuations in the substitution rate that we just described, and which could be strong enough to obscure the generation time effect. Thus, the statistics of neutral substitutions observed here could explain the quantitative difference between synonymous and non-synonymous substitutions.

### IV.9 Tests of neutrality

A better understanding of the statistical properties of neutral evolution can help to single out from a background of neutral substitutions the more interesting cases of positive selection as, for instance, changes in the protein function and responses to variations of the environment.

The best current bioinformatics method to identify such cases of positive selection is the one proposed by McDonald & Kreitman (1991), recently used with some modifications to study the evolution of Drosophila genes Smith & Eyre-Walker (2002); Fay et al. (2002). This method assumes implicitly that the neutral substitution process has a constant rate $`\mu x`$. Thus, the broad and correlated fluctuations in the neutral connectivities should be taken into account to provide a better null hypothesis to compare with.

McDonald and Kreitman’s method examines two samples of sequences from two related populations and distinguishes two classes of nucleotide differences in protein coding genes: intra-population differences (polymorphisms) and positions which are fixed in any two populations but differ between the two (fixed differences). Each class is further divided into synonymous (S) and non-synonymous (A) differences. Let us denote by $`D_{\mathrm{PS}}`$ and $`D_{\mathrm{PA}}`$ the number of synonymous and non-synonymous polymorphisms, respectively, and by $`D_{\mathrm{FS}}`$ and $`D_{\mathrm{FA}}`$ the number of synonymous and non-synonymous fixed differences, respectively. Under the neutral hypothesis, polymorphisms and substitutions are just the two faces of the same coin. Let us suppose that all synonymous mutations are neutral, and that all deleterious mutations are quickly removed from the population and do not show as polymorphisms. Under this hypothesis, the ratio $`D_{\mathrm{PA}}/D_{\mathrm{PS}}`$ of non-synonymous to synonymous polymorphisms can be interpreted as the fraction of neutral mutations, $`x`$, in the two populations. Analogously, $`D_{\mathrm{FA}}/D_{\mathrm{FS}}`$ can be interpreted to represent the fraction of neutral mutations, $`x`$, over an evolutionary trajectory. McDonald and Kreitman assume that the two quantities are equal under neutral evolution. Thus, in their method, the neutral expectation is that $`D_{\mathrm{PA}}/D_{\mathrm{PS}}`$ and $`D_{\mathrm{FA}}/D_{\mathrm{FS}}`$ are equal. Conversely, if the latter ratio is significantly larger than the former, this is interpreted as an evidence that some of the non-synonymous substitutions have been fixed by positive selection, since in this case fixation is much faster than for neutral mutations, and a larger number of substitutions can be fixed if they bring any selective advantage.

Let us now see how the neutral expectation must be modified taking into account the statistics of neutral substitutions. The ratio of non-synonymous to synonymous polymorphisms in two populations can be interpreted as the fraction of neutral mutations in the two populations, $`x(𝐀_1)`$ and $`x(𝐀_2)`$ respectively. On the other hand, the ratio $`D_{\mathrm{FA}}/D_{\mathrm{FS}}`$ should be interpreted as the average number of substitutions taking place in an evolutionary trajectory whose length $`t`$ is double of the time spent since the last common ancestor of the two populations,

$$\frac{D_{\mathrm{FA}}}{D_{\mathrm{FS}}}\frac{<S_t>}{\mu t},$$
(14)
where the average is the mutational average along the evolutionary trajectory which connects the two wild-type sequences $`𝐀_1`$ and $`𝐀_2`$. Now, if the number of substitutions between these two sequences is small with respect to the “correlation time” after which values of $`x(𝐀)`$ become uncorrelated, then we expect that the effective fraction of neutral mutations represented in Eq. (14) is close to both $`x(𝐀_1)`$ and $`x(𝐀_2)`$, and the neutral expectation of McDonald and Kreitman remains valid. Otherwise, we expect that $`D_{\mathrm{FA}}/D_{\mathrm{FS}}`$ is smaller than the larger between the two neutral fractions: As shown in Fig. 9, the effective substitution rate decreases as longer time intervals are considered, because the evolutionary process gets trapped in sequences with small neutral connectivity. As a rule of thumb, if $`x(𝐀_1)`$ and $`x(𝐀_2)`$ are very similar and the number of fixed substitutions $`D_{\mathrm{FA}}`$ is of the order of 10, then the neutral expectation of the McDonald and Kreitman method is likely to be valid, otherwise corrections should be considered.

## V Conclusions and perspectives

In this work we have proven that the statistical properties of neutral substitutions simulated through the SCN model, which imposes conservation of the thermodynamic stability of the protein, differ considerably from those predicted by the standard Kimura’s model, which assumes that the fraction of neutral mutations is constant throughout evolution. These differences have important consequences for several aspects of neutral evolution. Only taking into account accurate statistical properties of a neutral model allows to test the neutral hypothesis as the null hypothesis of molecular evolution, against which adaptive substitutions have to be distinguished.

We have shown that the differences between the standard neutral model and the SCN model arise because of correlations: (i) There are strong spatial correlations between the neutral connectivities of different positions along the protein chain, even if the corresponding amino acids are not close in space. As a consequence, the neutral connectivities of the entire sequences are broadly distributed, in contrast with the standard model. (ii) There are strong temporal correlations between the neutral connectivities of sequences nearby along a neutral trajectory. As a consequence, the substitution rates of different evolutionary trajectories, even averaged over a sizeable number of substitutions, can be significantly different.

We have illustrated a number of important consequences of these properties: (i) As a result of the broad distribution, the substitution rate is a decreasing function of the time interval considered, instead of being constant, and attains the stationary value $`\mu (\overline{1/x})^1`$ starting from the larger value $`\mu \overline{x}`$ at very small time intervals. (ii) As a result of the auto-correlations, even long evolutionary trajectories, corresponding to different species, can be characterized by rather different neutral substitutions rates. This fact can produce new lineage effects which may obscure the generation time effect, and it should be taken into account when testing the predictions of the neutral hypothesis. (iii) As a result of the broad variance of different evolutionary trajectories, the variance of the substitution process is not constant in time, and it is larger than that expected for a Poissonian process. Nevertheless, for large evolutionary separations the standard deviation of the number of substitutions is much smaller than the mean value, as for a Poissonian process, so that neutral substitutions can still be useful for reconstructing phylogenetic trees. (iv) Moreover, our simulations provide useful data to test methods for estimating the number of substitutions between two diverging sequences from their sequence similarity. Our results allow to estimate the threshold below which sequence similarity does not provide any information on the number of substitutions, and also show that the probability that the same amino acid arises as the result of two independent substitutions should be taken into account at low similarity. The neutral model does not make detailed predictions about this point, since one still needs to know the neutral mutation rates $`x_i`$ at different positions. It is important to note that Kimura’s model does not assume that such rates are all equal. Nevertheless, our simulations produce realistic $`K`$ vs $`N_\mathrm{d}/N`$ curves which can be used to test correction methods.

The structural approach to protein evolution has only just been initiated. Natural extensions of our model that will be considered in forthcoming work include consideration of the genetic code as well as modeling of the population genetics level and of possible fitness effects of changes in the thermodynamic parameters of the protein. We anticipate that the structural approach will provide key insights into different areas of molecular sciences, including directed evolution of enzymes, rational sequence design, phylogenetic reconstruction, protein fold prediction, and production of new materials for bio-nanotechnology.

## VI Materials and Methods

### VI.1 Protein model

We represent a protein structure by its contact matrix $`C_{ij}`$, where $`C_{ij}=1`$ if residues $`i`$ and $`j`$ are in contact and $`C_{ij}=0`$ otherwise. Two residues are considered in contact if any two of their heavy atoms are closer than $`4.5\text{Å}`$. The effective free energy associated to a sequence of amino acids $`𝐀`$ in the configuration $`𝐂`$ is approximated as a sum of pairwise contact interactions,

$$E(𝐀,𝐂)=\underset{i<j}{}C_{ij}U(A_i,A_j),$$
(15)
where $`A_i`$ labels one of the twenty amino acid types and $`U(,)`$ is a $`20\times 20`$ symmetric interaction matrix. Here we use the matrix derived by Bastolla et al. (2000b), which describes accurately the thermodynamic stability of a large set of monomeric proteins Bastolla et al. (2001).

Three remarks are needed: First, the effective energy parameters implicitly take into account the effect of the solvent and depend on temperature. They express free energies rather than energies. Second, the effective energy of a structure is defined with respect to a completely extended reference structure where no contacts are formed and which sets the zero of the energy scale. Third, one can derive from the database not the parameters $`U`$ themselves but the dimensionless quantities $`U/(k_\mathrm{B}T)`$. It is thus important to use dimensionless parameters to evaluate the stability of the protein model.

### VI.2 Candidate structures

We generate candidate structures for a protein sequence of $`N`$ residues by generating all possible gapless alignments of the sequence with structures in the PDB. This procedure is called threading. In this way, we generate many, typically of the order of $`10^6`$, protein-like structures per sequence. In the present context, threading is directly used to produce the contact maps of the candidate structures. In order to speed up the computations, we use a non-redundant subset of the PDB excluding proteins with homologous sequences, selected by Hobohm & Sander (1994).

### VI.3 The folding parameter $`\alpha `$

For a given sequence $`𝐀`$, the energy landscape is well correlated if all configurations of low energy are very similar to the configuration of minimal effective energy, $`C^{}(𝐀)`$. Structure similarity is measured by the overlap $`q(𝐂,𝐂^{})`$, counting the number of contacts that two structures have in common and normalizing it through the maximal number of contacts, so that $`q`$ is comprised between zero and one. In a well correlated energy landscape, the inequality

$$\frac{E(𝐀,𝐂)E(𝐀,𝐂^{})}{|E(𝐀,𝐂^{})|}\alpha (𝐀)\left(1q(𝐂,𝐂^{})\right)$$
(16)
holds, stating that the energy gap between each alternative structure $`𝐂`$ and the ground state $`𝐂^{}`$, measured in units of the ground state energy, is larger than a quantity $`\alpha (𝐀)`$ times the structural distance $`1q(𝐂,𝐂^{})`$. The dimensionless quantity $`\alpha (𝐀)`$, which is the largest quantity for which the above inequality holds, can be used to evaluate the folding properties of sequence $`𝐀`$. For random sequences, many different configurations have quite similar energy and hence $`\alpha (𝐀)0`$. In this case the energy landscape is rugged, the folding kinetics is very slow, and the thermodynamic stability with respect to variations in the solvent is very low. In contrast, computer simulations of well designed sequences have shown that, when $`\alpha (𝐀)`$ is finite, the folding kinetics is fast and the stability with respect to changes in the energy parameters as well as mutations in the sequence is very high.

Our algorithm computes the parameter $`\alpha (𝐀)`$ for a fixed target configuration $`𝐂^{}`$ and a large number of sequences $`𝐀`$. We thus indicate this parameter as $`\alpha (𝐀,𝐂^{})`$, since we do not know a priori that $`𝐂^{}`$ has lowest energy. Notice however that, if $`\alpha (𝐀,𝐂^{})`$ is positive, all alternative structures have higher energy than $`𝐂^{}`$. We impose that $`\alpha (𝐀,𝐂^{})`$ is larger than a positive threshold $`\alpha _{\mathrm{thr}}`$ for sequences $`𝐀`$ belonging to the neutral network.

### VI.4 The $`Z`$-score

The $`Z`$-score $`Z(𝐀,𝐂^{})`$ Bowie et al (1991); Goldstein et al. (1992) is a measure of the compatibility between a sequence $`𝐀`$ and a structure $`𝐂^{}`$, widely used in structure prediction. It depends on an effective energy function, and measures the difference between the energy of sequence $`𝐀`$ in configuration $`𝐂^{}`$ and its average energy in a set of alternative configurations, $`\{𝐂\}`$, in units of the standard deviation of the energy,

$$Z(𝐀,𝐂^{})=\frac{E(𝐀,𝐂^{})<E(𝐀,𝐂)>_𝐂}{\sqrt{<E(𝐀,𝐂)^2>_𝐂<E(𝐀,𝐂)>_𝐂^2}}.$$
(17)
When sequence $`𝐀`$ folds in structure $`𝐂^{}`$, its corresponding $`Z`$-score is very negative.

Given the above definition, one has still to specify how to select alternative structures. A possibility, often used for lattice models Mirny & Shakhnovich (1996) is to assume that alternative structures are maximally compact, randomly chosen structures, whose average energy can be estimated as $`<E(𝐀,𝐂)>_𝐂=Nc_{\mathrm{max}}<e(𝐀)>`$. Here, $`Nc_{\mathrm{max}}`$ is the maximal number of contacts of candidate structures and $`<e(𝐀)>`$ is the average energy of a contact, averaged over all possible contacts formed by sequence $`𝐀`$. This leads to introduce the parameter

$$Z^{}(𝐀,𝐂^{})=\frac{E(𝐀,𝐂^{})/Nc_{\mathrm{max}}<e(𝐀)>}{\sqrt{<e^2(𝐀)><e(𝐀)>^2}}.$$
(18)
The use of $`Z^{}`$ has two main advantages: (i) It makes the value of the $`Z`$-score much less sensitive to chain length $`N`$ and to the particular set of alternative structures used. (ii) The evaluation of $`Z^{}`$ is much faster than that of the $`Z`$-score. This is necessary in order to explore efficiently sequence space. We impose that $`Z^{}(𝐀,𝐂^{})`$ is larger than a positive threshold $`Z_{\mathrm{thr}}`$ for sequences $`𝐀`$ belonging to the neutral network.

### VI.5 Sampling the neutral network

Our algorithm explores the neutral network of a given protein starting from its PDB sequence $`𝐀_1`$ and iterating the following procedure: At iteration $`i`$, (i) the number $`X(𝐀_i)`$ of viable neighbors of sequence $`𝐀_i`$ is computed, and (ii) the sequence $`𝐀_{i+1}`$ is extracted at random among all the viable neighbors of $`𝐀_i`$. In this way we generate a stochastic process along the neutral network which simulates neutral evolution and looses memory of the initial sequence very fast.

Sequence $`𝐀`$ is regarded as viable if both parameters $`\alpha (𝐀,𝐂^{})`$ and $`Z^{}(𝐀,𝐂^{})`$ are above predetermined thresholds, chosen as 98.5% of the values of those parameters for the sequence in the PDB. This enforces conservation of the thermodynamic stability and folding capability of the native structure $`𝐂^{}`$. We verified that the observed behavior does not change qualitatively for thresholds between 95% and about 100% of the PDB values.

We impose strict conservation of the cysteine residues in the PDB sequence, and do not allow any residue to mutate to cysteine, since a mutation changing the number of cysteine residues by one would leave the protein with a very reactive unpaired cysteine that would probably affect its functionality. Accordingly, the maximal possible number of neighbors tested is $`X_{\mathrm{max}}=18(NN_{\mathrm{cys}})`$, where $`N`$ is the number of residues and $`N_{\mathrm{cys}}`$ is the number of cysteine residues in the starting sequence. The total number of viable point mutations, $`X(𝐀)`$, expresses the local connectivity of the neutral network. We normalize it by the total number of neighbors, $`X_{\mathrm{max}}`$, getting the fraction of neutral neighbors, $`x(𝐀)=X(𝐀)/X_{\mathrm{max}}(0,1]`$.

To compute $`x(𝐀)`$, we have to evaluate the $`\alpha `$ parameter for all sequences $`𝐀^{}`$ obtained through a point mutation of sequence $`𝐀`$. From Eq. (16) we note that the $`\alpha `$ parameter can be obtained from the configuration with the highest destabilizing power, i.e. the highest value of the energy gap divided by the structural distance from the native configuration. These change from sequence to sequence, but it is expected not to change very much for neighboring sequences. Thus, in order to speed up the computation of $`\alpha (𝐀^{},𝐂^{})`$, instead of considering all candidate configurations we consider only the 50 configurations with the highest destabilizing power (i.e. the energy gap divided by the structural distance from the native configuration) for sequence $`𝐀`$ and compute their mutated destabilizing power using sequence $`𝐀^{}`$. The $`\alpha `$ parameter is then obtained from the configuration with the highest destabilizing power. This procedure could slightly overestimate $`\alpha (𝐀^{},𝐂^{})`$ since not all configurations are used, but we have checked that the error introduced in the $`x`$ value is in all cases below 0.1%.

### VI.6 Cross-correlations

To obtain a deeper insight into the mechanism of neutral evolution, it is helpful to study the cross-correlations between the fraction of neutral neighbors for two given residues of a protein of $`N`$ residues. The analysis starts from $`N`$ individual trajectories $`\{x_i(𝐀_1),x_i(𝐀_2),\mathrm{}\}`$ of $`m`$ evolutionary steps for each residue $`i`$. We define the average fraction of neutral neighbors $`\overline{x_i}=(1/m)_{k=1}^mx_i(𝐀_k)`$ and the corresponding variance $`\sigma _i^2=\overline{x_i^2}\overline{x_i}^2`$ for each residue $`i`$. Then, we calculate the correlation matrix $`C_{ij}`$, where $`𝒞_{ij}=𝒞_{ji}`$ determines the cross-correlations between residue $`i`$ and $`j`$ and is defined as

$$𝒞_{ij}=\frac{1}{m}\underset{k=1}{\overset{m}{}}\frac{\left(x_i(𝐀_k)\overline{x_i}\right)\left(x_j(𝐀_k)\overline{x_j}\right)}{\sigma _i\sigma _j}.$$
(19)
Due to the enforced conservation of cysteine, these residues require a special treatment: If residue $`i`$ is cysteine, then $`𝒞_{ij}=𝒞_{ji}=0`$ for all $`j`$.

One needs to distinguish three different cases, (i) $`𝒞_{ij}=0`$ means that the fraction of neutral neighbors for residues $`i`$ and $`j`$ are uncorrelated, (ii) $`𝒞_{ij}>0`$ indicates that the fraction of neutral neighbors for residues $`i`$ and $`j`$ are correlated, and (iii) $`𝒞_{ij}<0`$ shows that the fraction of neutral neighbors for residues $`i`$ and $`j`$ are anti-correlated.

Alternatively, one may also study the cross-correlations between the fraction of neutral neighbors for one given residue and the fraction of neutral neighbors for the whole protein. Defining the average fraction of neutral neighbors for the whole protein $`\overline{x}=(1/m)_{k=1}^mx(𝐀_k)`$ and the corresponding variance $`\sigma ^2=\overline{x^2}\overline{x}^2`$, we calculate the correlation vector $`𝒞_i`$ which determines the cross-correlations between the fraction of neutral neighbors for residue $`i`$ and the fraction of neutral neighbors for the whole protein,

$$𝒞_i=\frac{1}{m}\underset{k=1}{\overset{m}{}}\frac{\left(x_i(𝐀_k)\overline{x_i}\right)\left(x(𝐀_k)\overline{x}\right)}{\sigma _i\sigma }.$$
(20)
Here, again three different cases need to be distinguished, (i) $`𝒞_i=0`$ indicates that the fractions of neutral neighbors for residues $`i`$ and for the whole protein are uncorrelated, (ii) $`𝒞_i>0`$ means that the fractions of neutral neighbors for residues $`i`$ and for the whole protein are correlated, and (iii) $`𝒞_i<0`$ shows that the fractions of neutral neighbors for residues $`i`$ and for the whole protein are anti-correlated.

### VI.7 Substitution process

Given an evolutionary trajectory $`\{x(𝐀_1),x(𝐀_2),\mathrm{}\}`$, the distribution of the number of substitutions taking place in a time $`t`$ can be computed by considering Eq. (5), where $`k`$, the number of attempted mutations, is a Poissonian variable of average value $`\mu t`$.

In order to handle the computation, we divide all values of $`x`$ in $`M`$ classes, choosing $`x_c`$ as representative value of all $`x`$’s belonging to class $`c`$. The number of operations needed to evaluate the substitution probability increases exponentially with the number of classes $`M`$. At the same time, the evaluation becomes more and more accurate as $`M`$ increases. We chose $`M=6`$ in our numerical computations as a reasonable compromise between accuracy and rapidity, checking that larger values of $`M`$ introduce only small changes.

Abbreviations: PDB: Protein Data Bank, SNC model: Structurally Constrained Neutral Model


format: text

text:
# Multiplicity, kinematics and rotation rates of very young brown dwarfs in Cha I

## 1. Introduction

Although a large number of brown dwarfs have been detected up to date, fundamental parameters have been measured only for a small subset. Consequently, there are still a lot of open questions, for example, the mechanism that leads to the formation of brown dwarfs is still poorly constrained.

We have studied a population of twelve very young bona fide and candidate brown dwarfs in the Cha I star forming region, Cha H$`\alpha `$ 1 to 12 (Comerón et al. 1999, 2000, Neuhäuser & Comerón 1998, 1999) by means of high-resolution spectroscopy ($`\lambda /\mathrm{\Delta }\lambda =\mathrm{40\hspace{0.17em}000}`$) with UVES at the 8.2 m Kueyen telescope of the VLT as well as by a photometric monitoring campaign carried out at the Danish 1.5 m telescope at ESO. These observations yielded radial and rotational velocities as well as rotational periods. The studied brown dwarfs are only a few million years old and their observation allows insights into the formation and early evolution of brown dwarfs.

## 2. UVES spectra I: Kinematic study

High-resolution spectroscopic observations were conducted with the echelle spectrograph UVES attached to the 8.2 m Kueyen telescope of the VLT. At least two spectra separated by a few weeks have been taken of nine of the twelve targets. The wavelength regime from 6600 Å to 10400 Å was covered with a spectral resolution of $`\lambda /\mathrm{\Delta }\lambda =\mathrm{40\hspace{0.17em}000}`$. Radial velocities have been measured by a cross-correlation of plenty of stellar lines of the object spectra against a template spectrum. Telluric lines have been used as wavelength reference. The achieved precision of the relative radial velocities range between 80 m s-1 and 600 m s-1, depending on the S/N of the individual spectra.

Recently, it has been proposed, that brown dwarfs are ejected protostars which have been cut off from the gas reservoir in the early accretion phase of the star formation process and could have therefore not accreted to stellar masses (Reipurth & Clarke 2001). A kinematic study of the sample (Joergens & Guenther 2001) showed that their radial velocity dispersion is relatively small (2.2 km s-1). It is significantly smaller than the radial velocity dispersion of the T Tauri stars in the field (3.6 km s-1) and slightly larger than that one of the surrounding molecular gas (1.2 km s-1). The distributions are displayed in Fig. 1. These results give suggestive evidence that there is no run-away brown dwarf among the studied sample. It cannot be ruled out that some of them have a larger 3D velocity dispersion. However, the studied brown dwarfs occupy a field of less than 12$`\times `$ 12 at a distance of 160 pc. Brown dwarfs born within this field and ejected during the early accretion phase in directions with a significant fraction perpendicular to the line-of-sight, would have flown out of the field a long time ago, given the age of the objects of 1–5 Myr. This is true even for the smallest ejection velocities of 2 km s-1 calculated by theorists (see also next paragraph). Therefore, the measurement of radial velocities is the very method to test if objects born in the field have significantly high velocities due to dynamical interactions during their formation process.

Very recent dynamical calculations (Bate et al. 2002, Sterzik & Durisen this volume) hint at rather small ejection velocities suggesting the possibility that the imprint of the ejection in the kinematics might not be an observable effect. They predict a velocity dispersion of ejected brown dwarfs of 2 km s-1 in 3D and 1.2 km s-1 in 1D. On the other hand, they predict that 10% of the brown dwarfs have a larger velocity than 5 km s-1 due to dynamical interactions. The velocities of the bona fide and candidate brown dwarfs in Cha I cover a total range of only 2.6 km s-1. Therefore, it is concluded that the ejection-model for the formation of brown dwarfs is not a likely formation mechanism for the brown dwarfs in Cha I.

## 3. UVES spectra II: RV Survey for companions

By means of time-resolved UVES spectra, a radial velocity survey for close companions to the targets was conducted. For each of the objects Cha H$`\alpha `$ 1, 2, 3, 5, 6, 7, 8 and 12 spectra have been taken in two nights separated by a few weeks. For Cha H$`\alpha `$ 4 spectra have been taken in 11 nights, separated by weeks to years.

We have found that the radial velocities for Cha H$`\alpha `$ 1, 2, 3, 5, 6, 7, 8 and 12 are constant within the measurements uncertainties, whereas Cha H$`\alpha `$ 4 shows small but significant variations (Fig. 2). Details will be published in a forthcoming paper (Joergens et al., in prep.).

An upper limit for the semiamplitude of a variation caused by a hypothetical companion was determined for each object by assuming that the total variability amplitude was recorded. Subsequently, upper mass limits for the hypothetical companions have been estimated. For this calculation an orbital separation of 0.1 AU was assumed. Primary masses from Comerón et al. (1999, 2000) have been used. Furthermore, circular orbits and an inclination of 90 were assumed. The assumption of $`i=90^{}`$ means that the derived mass limits for companions are upper limits for the minimum masses M$`{}_{2}{}^{}\mathrm{sin}i`$.

The upper limits for the semiamplitude of the ’constant’ objects (Fig. 2, left panel) range between 30 m s-1 for Cha H$`\alpha `$ 5 and 600 m s-1 for Cha H$`\alpha `$ 7 and 12. The deduced upper limits for M$`{}_{2}{}^{}\mathrm{sin}i`$ range between a tenth of a Jupiter mass (Cha H$`\alpha `$ 5) and 1–2 Jupiter masses (Cha H$`\alpha `$ 7 and 12). The small but significant variations of the radial velocity of Cha H$`\alpha `$ 4, if caused by a companion, correspond to a companion mass M$`{}_{2}{}^{}\mathrm{sin}i`$ of 0.8 MJup.

In summary, no evidence for substellar companions around the targets has been found, aside from variations of the radial velocity of Cha H$`\alpha `$ 4, which could be due to a companion of less than 1 MJup. There is the possibility that there are planetary mass companions with masses below one Jupiter mass since the presented observations were not sensitive enough to detect those. Furthermore, spectroscopic binaries in the sample might not have been detected due to the non-observation at the critical orbital phases. However, the overall picture that arises is that the multiplicity fraction of the studied brown dwarfs is rather low ($``$ 10%).

## 4. Photometric monitoring: Rotational periods

A photometric monitoring campaign has been carried out at the Danish 1.5 m telescope at ESO in order to study the time dependence of the brightness of the objects. Rotational periods have been measured by tracing modulations of the light curves due to magnetically induced surface spots. Rotational periods have been determined for the three brown dwarf candidates Cha H$`\alpha `$ 2, Cha H$`\alpha `$ 3 and Cha H$`\alpha `$ 6 of 2.8 d, 2.2 d and 3.4 d, respectively (Joergens et al. 2002a,b). See Fig. 3 and 4 for phase-folded light curves. The periods are consistent with measurements of rotational velocities $`v\mathrm{sin}i`$ from UVES spectra for the objects (8–26 km s-1). The observations show that brown dwarfs at an age of 1–5 Myr display surface spots like T Tauri stars and are moderately fast rotators in contrast to rapidly rotating old brown dwarfs consistent with them being in an early contracting stage. A comparison with rotational periods at 36 Myr (Eislöffel & Scholz 2001) indicates that most of the acceleration of brown dwarfs takes place in the first 30 million years or less of their lifetime (cf. Fig. 4, right panel).

## 5. Discussion and conclusions

The kinematic study of brown dwarfs in Cha I based on the measurement of precise radial velocities from high-resolution UVES spectra showed that their radial velocities have a dispersion of only 2.2 km s-1 and span a total range of only 2.6 km s-1. We therefore conclude that the ejection-model for the formation of brown dwarfs (Reipurth & Clarke 2001) is not a likely formation mechanism for the studied brown dwarfs since first calculations (Bate et al. 2002, Sterzik & Durisen 2002, this volume) predict that 10% of the brown dwarfs should have a larger velocity than 5 km s-1 due to the ejection.

The small binary fraction ($``$ 10%) found in the presented RV survey is in agreement with the result of a direct imaging survey for wide, low-mass companions to the same objects (Neuhäuser et al. (2002), Joergens et al. (2001), Neuhäuser et al. this volume). Combining both surveys, a wide range of possible companion separations has been covered. The exact separations depend on the companion masses. For example, for brown dwarf companions ($`<`$ 13 MJup) to the targets, separations $`<`$3 AU and between 50 and 1000 AU were covered. With more restricted separations ($`<`$0.1 AU and 300–1000 AU) the surveys were sensitive also to companion masses down to 1 MJup.

These results seem to be in contrast with the high multiplicity fraction observed for very young stars and may suggest that brown dwarfs form not by direct collapse of unstable cloud cores as stars. However, a significant comparison of the multiplicity fraction of stars and brown dwarfs at very young ages is still hampered by small-number statistics in the substellar regime. Such a study has to compare multiplicity fractions in a certain separation range, which has to agree for both samples. Furthermore, it is possible that the gravitational collapse of cloud cores of brown dwarf masses does not yield as much multiple systems as the collapse of clouds of solar masses. In addition, stars can have companions, which have only a tenth of the mass of the primary, whereas a companion of a brown dwarf with such a low mass ratio would be already a planet. However, there is still a lot to do for observers as well as for theorists in order to understand in which way objects of brown dwarf masses are formed.

Finally, we have determined rotational periods for the three brown dwarf candidates Cha H$`\alpha `$ 2, Cha H$`\alpha `$ 3 and Cha H$`\alpha `$ 6 of 2.2 to 3.4 days by means of a photometric monitoring campaign. These are the first rotational periods for very young brown dwarfs and among the first for brown dwarfs at all. They provide valuable data points in an as yet almost unexplored region of the substellar period-age diagram. A comparison of the determined rotational periods at the age of a few million years with rotational properties of older brown dwarfs ($`>`$36 Myr, Eislöffel & Scholz 2001) shows that most of the acceleration of brown dwarfs during the contraction phase takes place within the first 30 Myr or less. It is known that Cha H$`\alpha `$ 2 and 6 have optically thick disks (Comerón et al. 2000), therefore magnetic braking due to interactions with the disk may play a role for them. This is suggested by the fact, that among the three brown dwarf candidates with determined periods, the one without a detected disk, Cha H$`\alpha `$ 3, has the shortest period. If the interaction with the disk is responsible for the braking, the results from Eislöffel & Scholz (2001) indicate that brown dwarf inner disks have been dissipated at an age of 36 Myr. These limits for the time scale of disk dissipation for brown dwarfs are very similar to those for T Tauri stars, which are known to dissolve their inner disks within about the first 10 Myr (e.g. Calvet et al. 2000).

Further measurements of rotational periods for brown dwarfs are much-needed in order to complete the picture of angular momentum evolution in the substellar regime as well as that of rotationally induced phenomena, like dynamo activity and meteorological processes.

## References

Bailer-Jones C.A.L. & Mundt R. 2001, A&A 367, 218

Bate M.R., Bonnell I.A. & Bromm V. 2002, MNRAS 332, L65

Calvet N., Hartmann L. & Strom S.E. 2000, ’Evolution of Disk Accretion’, In: *Protostars and Planets IV*, ed. by Mannings V., Boss A.P., Russell S.S., Univ. of Arizona Press, Tucson, p. 377

Clarke F.J., Tinney C.G. & Covey K.R. 2002, MNRAS 332, 361

Comerón F., Rieke G.H. & Neuhäuser R. 1999, A&A 343, 477

Comerón F., Neuhäuser R. & Kaas A.A. 2000, A&A 359, 269

Eislöffel J. & Scholz A. 2001, ’Rotation and Atmospheres of Brown Dwarfs’. In: *Origins of stars and planets: The VLT view*, ESO workshop, Garching, Germany, April 24-27, 2001, in press

Joergens V., Guenther E., Neuhäuser R. et al. 2001 ’Multiplicity of young brown dwarfs in Cha I’. In: *Origins of stars and planets: The VLT view*, ESO workshop, Garching, Germany, April 24-27, 2001, in press, astro-ph/0106185

Joergens V. & Guenther E.W. 2001, A&A 379, L9

Joergens V., Fernandéz M., Neuhäuser R. & Guenther E.W. 2002a, Poster proceeding of the Star Spot Thinkshop, held in Potsdam, Germany in May 2002, in press, astro-ph/0206137

Joergens V., Fernandéz M. & Neuhäuser R. 2002b, A&A subm.

Martín E.L., Zapatero Osorio M.R. & Lehto H.J. 2001, ApJ 557, 822

Neuhäuser R. & Comerón F. 1998, Science 282, 83

Neuhäuser R. & Comerón F. 1999, A&A 350, 612

Neuhäuser R. & Brandner W., Alves J., Joergens V., Comerón F. 2002, A&A 384, 999

Reipurth B. & Clarke C. 2001, ApJ 122, 432


text:
# Observation of stochastic resonance in percolative Josephson media

## .1 1. Introduction

The phenomenon of *stochastic resonance* has been discussed in relation to diverse problems in nonlinear science, physics, chemistry and biology Hanggi . Generally speaking, stochastic resonance is the enhancement of the output signal-to-noise ratio caused by the injection of an optimal amount of noise into a periodically driven nonlinear system. This kind of behavior is often thought as counterintuitive, since here a stochastic force amplifies a small periodic signal. Its mechanism is usually explained in terms of motion of a particle in a double-well potential subjected to noise, in the presence of a time-periodic force. The periodic forcing leads to noise-enhanced transitions between the two wells and thus to an enhanced output of the forcing signal.

A clear example of nonlinear systems with few degrees of freedom is a superconducting loop with a Josephson junction, well known as a superconducting quantum interferometer (SQUID). With a proper choice of the size of the loop, this system undergoes bistable dynamics for magnetic flux trapped in the loop. There have already been experiments that reported operating SQUIDs under stochastic resonance conditions, both with external noise injection Hibbs  and with thermally generated intrinsic noise Lukens . The stochastic resonance effect can be considerably enhanced in a system of coupled bistable oscillators (see, e.g., Ref. Bulsara ). Therefore, it is interesting to study stochastic amplification for a Josephson media consisting of many superconducting contours with Josephson junctions.

Earlier we observed quantum interference effects in macroscopically inhomogeneous superconducting Sn-Ge thin-film composites near the percolation threshold PhysicaB . This system exhibits a considerable voltage noise under dc current bias and a rectification of ac current, which arise below the superconducting transition temperature. According to Ref. Ouboter , a dc voltage is observed when an ac current larger than the critical current passes through a system of two superconductors weakly connected by an asymmetric double point contact, i.e., the magnetic flux quantization induces the critical current oscillations and the respective voltage oscillations. We have argued PhysicaB , that the oscillatory dependence $`V_{dc}(H)`$ in Sn-Ge thin-film composites is related to quantum interference in randomly distributed asymmetric superconducting contours interrupted by Josephson *weak links*.

In Ref. PhysicaB  we reported measurements of the $`V_{dc}(H)`$ dependence for various orientations of the film relative to the field. The scale of the oscillatory structure in $`V_{dc}(H)`$ is inversely proportional to the cosine of the angle between the applied magnetic field and the normal to the sample plane. The emergence of the normal magnetic field component alone and also the antisymmetry of the oscillatory structure relative to $`H=0`$ indicate a quantum-interference origin of $`V_{dc}(H)`$. Moreover, it appears feasible to relate these active contours to the percolative cluster that has a well-known *fractal* structure. The existence of a wide and *self-similar* distribution of Josephson contour areas leads to a fractal character of the dependence $`V_{dc}(H)`$. We have suggested and verified a model for the origin of the $`1/f`$ voltage noise by a passive transformation of magnetic field oscillations with a fractal transfer function $`V_{dc}(H)`$ PhysicaB .

In the present paper, we study the noise-induced electrical response of granular Sn-Ge thin-film composites. We argue that a distributed network containing many superconducting loops with Josephson junctions may show a cooperative behavior as *stochastically resonating media*.

## .2 2. Experimental details and results

Josephson networks may occur naturally, e.g., in nonuniform superconducting materials such as granular thin films. We prepare granular Sn-Ge thin-film composites having monotonically varying structure by vacuum condensation of Sn on a long (60 mm) substrate along which a temperature gradient is created. Sn is deposited on the previously prepared 50 nm thick Ge layer. The thickness of the Sn layer is 60 nm. The metallic condensate is covered from the top with amorphous Ge. The structural change results in variation of the composite properties from metallic to insulating over the substrate. This crossover in properties is observed on a series consisting of 30 samples cut from different parts of the substrate. For the present investigations, we chose samples with properties near the percolation threshold, with a characteristic structure depicted in Fig. 1.

During the measurements, the samples were kept in exchange gas inside a superconducting solenoid. The electrical measurements were carried out according to the standard four-probe technique. A sinusoidal ac current of frequency $`f_I=100`$ kHz and amplitude $`I_{ac}=0.8`$ mA was produced by an HP3245A universal source connected to the current leads through a dc-decoupling transformer. Fast Fourier transformation spectra of the output voltage are measured by using an SR770 spectrum analyzer with a Blackman-Harris window function. We used the signal-to-noise ratio (SNR) as the major characteristic of stochastic resonance. The SNR was measured as the ratio of the voltage amplitude of the spectral line to the voltage noise level below it. The noise background in the signal bin is estimated by performing a linear fit to the peak clipped spectrum. The noise intensity (noise level) denotes the standard deviation $`\sigma _N`$ of the Gaussian white noise signal, which was supplied by the internal SR770 generator.

The transition of a sample into the superconducting state is smeared over 1.0 K with the center of the resistive transition at $`T_0=3.8`$ K. At temperatures below $`T_0`$ and with ac current $`I_{ac}`$ applied through the sample, we observed a rectified dc voltage $`V_{dc}`$, the magnitude of which oscillated as a function of the dc magnetic field $`H`$ applied perpendicular to the substrate (Fig. 2a). The amplitude and frequency of the current $`I_{ac}`$ did not significantly affect the general features of the $`V_{dc}(H)`$ dependence. The results could be always readily reproduced.

To observe the phenomenon of stochastic resonance, we study the rectified voltage dependence on magnetic field. The applied magnetic field consisted of three components (Fig. 2b): (i) a dc field $`H`$, which varied in the range between $`300`$ and $`+300`$ mOe, (ii) a small ac component with a frequency $`f_H`$ between 5 and 60 Hz and an amplitude $`H_{ac}=20`$ mOe, and (iii) Gaussian white noise $`H_{noise}`$ with the intensity $`\sigma _N`$ ranging up to 70 mOe. The Fourier spectra of voltage response are shown in Fig. 3 together with oscillograms of the input signal $`H_{ac}+H_{noise}`$. Figure 4 shows the dependence of the output SNR for the first harmonic of $`f_H`$ on the intensity of input noise $`H_{noise}`$. One can see that increasing the noise amplitude at first increases the SNR and then decreases it. Such maxima are rather characteristic for the phenomenon of stochastic resonance. Similar measurements taken at different magnetic fields and frequencies often showed multiple maxima such as those shown in Fig. 5.

## .3 3. Discussion

In summary, our experiments demonstrate the characteristic feature of the phenomenon of stochastic resonance, namely the nonmonotonic behavior of the SNR. At the optimum noise level the SNR increases up to 40. The presence of multiple maxima (Fig. 4 and 5) can be due to the effect of different Josephson contours in our structure, which is operated at the border of the percolation threshold.

We suppose that the nonmonotonic dependence of the SNR on frequency $`f_H`$ (Fig. 5) excludes other possible explanations (such as, e.g., a simple rectification effect due to a nonlinearity of the response) for the observed gain of a small input signal.

Detailed measurements taken at different frequencies, shown in Fig. 5 indicate, at least in some ranges of the dc magnetic field, the existence of parameter regions characterized by a significant gain for a relatively broadband signal. We interpret this behavior as a property of percolative Josephson media with a wide range of self-similar loops. The SNR gain in our system can be tuned to a desired operation frequency $`f_H`$ by changing the dc magnetic field H.

The nature of the stochastic resonance in the system studied can be related to the commonly known bistable oscillator behavior of the magnetic flux quantization loops. Moreover, in the presence of current bias $`I_{ac}`$ at relatively high frequency (at $`f_I`$ about 100 kHz) with amplitude larger than critical, our samples exhibit dynamical chaos. Such a regime is commonly characterized by a coexistence of multiple attractors in the phase space. Indeed, calculation of Lyapunov exponents from the time evolution of the voltage measured at constant current indicates presence of chaos in our system LT22 . In this case, the ”phase trajectory” of the system may stay long time in any of the attractors and perform irregular transitions between them. Synchronization of such intermittent transitions by a small input signal may lead to stochastic resonance as well Anischenko . Yet, these speculations require further investigations to be firmly justified.

This work was supported by the German-Ukrainian collaboration grant of the Bundesministerium für Bildung, Wissenschaft, Forschung und Technologie (BMBF) No. UKR-003-99.