Should We Expect Substitution Rate to Depend on Population Size?

Corresponding author: Joshua L. Cherry, Department of Human Genetics, University of Utah, 15 N 2030 E RM 2100, Salt Lake City, UT 84112-5330., cherry{at}genetics.utah.edu (E-mail).

Communicating editor: R. R. HUDSON

	ABSTRACT

TOP ABSTRACT THE MODEL NUMERICAL RESULTS DISCUSSION APPENDIX 1 LITERATURE CITED

The rate of nucleotide substitution is generally believed to be a decreasing function of effective population size, at least for nonsynonymous substitutions. This view was originally based on consideration of slightly deleterious mutations with a fixed distribution of selection coefficients. A realistic model must include the occurrence and fixation of some advantageous mutations that compensate for the loss of fitness due to deleterious substitutions. Some such models, such as so-called "fixed" models, also predict a population size effect on substitution rate. An alternative model, presented here, predicts the near absence of a population size effect on substitution rate. This model is based on concave log-fitness functions and a fixed distribution of mutational effects on the selectively important trait. Simulations of an instance of the model confirm the approximate insensitivity of the substitution rate to population size. Although much experimental evidence has been claimed to support the existence of a population size effect, the body of evidence as a whole is equivocal, and much of the evidence that is supposed to demonstrate such an effect would also suggest that it is very small. Perhaps the proposed model applies well to some genes and not so well to others, and genes therefore vary with regard to the population size effect.

THE probability of fixation of a mutant allele with a given selection coefficient (s) depends on the effective population size (N_e). Specifically, the ratio of the probability of fixation of a newly arising allele to that for a neutral allele is a function of the product N_es. If we define S to be equal to 2N_es for a haploid population, or 4N_es for a diploid population, and if |s| is small, then this ratio is approximately equal to

(1)

(KIMURA 1962 ). This is often approximated for deleterious mutations by the statement that alleles with |s| less than some threshold, such as 1/N_e, are "effectively neutral."

This dependence on effective population size, along with assumptions of a constant distribution of s among mutants and a constant mutation rate, has been used to derive a dependence of substitution rate (k) on N_e. The exact form of the relationship between N_e and substitution rate will depend upon the mutational spectrum. OHTA 1977 considered an exponential distribution for -s, and KIMURA 1979 considered a gamma distribution with ß = . With an additional assumption, they derived the approximate relationships that k 1/N_e and k , respectively.

These treatments considered deleterious mutations, which presumably are much more common than advantageous ones. The occurrence of only deleterious mutations, and the fixation of some fraction of them, would lead to a continuous decline in fitness, which seems unrealistic. For this reason, it is generally acknowledged that there must be at least occasional fixation of advantageous mutations. The view that emerges is that the mutational spectrum includes both deleterious and advantageous alleles, with the former being much more common than the latter. While the average selection coefficient is negative among mutants, the action of selection is such that it is zero among substitutions ("accepted" mutations). A sort of balance exists between the tendency of mutations to be deleterious and the higher probability of fixation of more advantageous alleles.

There is a problem with this view if we consider the mutational spectrum of s to be fixed even as substitutions occur. The problem is that the average value of s among substitutions will be a strictly increasing function of N_e. This means that a "balance" will be achieved for only a single value of N_e. Clearly, real populations will not usually have precisely this critical N_e. For smaller populations, one would again predict a continuous decline in fitness, albeit a somewhat slower one. For larger populations, one would expect a continuous increase in fitness, which also seems unrealistic. The assumptions of this model must be incorrect. It must be that as fitness decreases (e.g., because of a decrease in N_e), some factor increases the average of s among substitutions and eventually stops the decline in fitness and establishes a steady state. Similarly, this same factor must eventually stop the increase in fitness when N_e is increased. This factor must exert its effect through a dependence of the distribution of selection coefficients on the fitness of the parental allele.

One way in which the mutational spectrum of s might change with fitness is that as more deleterious mutations are accumulated, more ways to improve the sequence and fewer ways to make it worse are available. If we simplistically imagine that each site in a molecule can be either optimal or not optimal, and that fitness is a function of the number of optimal sites, then it follows that the less fit an allele, the higher the ratio of advantageous to disadvantageous mutations. This kind of model may be a good description of the evolution of codon usage, where for a small population or small selection coefficients (due, e.g., to a low level of gene expression) a state of mutational equilibrium will be approached. Such models can, like models with only deleterious mutations, lead to a population size effect on substitution rate. "Fixed" models (OHTA and TACHIDA 1990 ), in which the distribution of fitness among mutant alleles is independent of the fitness of the parent allele, can be thought of as rather extreme instances of this type of model.

I present here an alternative to such models that may more realistically describe the evolution of protein sequences. This model, like those discussed above, involves changes to the mutational distribution of selection coefficients as alleles of different fitness become fixed. However, the reason for these changes is different, as is the nature of the changes. This model, unlike those discussed above, predicts the near absence of a population size effect on substitution rate. I present results of calculations and computer simulations for an instance of this model. These results confirm that substitution rate is approximately independent of population size. I then discuss the model in light of experimental evidence and certain theoretical concerns.

	THE MODEL

TOP ABSTRACT THE MODEL NUMERICAL RESULTS DISCUSSION APPENDIX 1 LITERATURE CITED

The model that is the main topic of this article was inspired by work on the evolutionary implications of concave fitness functions (HARTL et al. 1985 ), although concavity is not a requirement (see below). Suppose that the distribution of mutants is constant, not with regard to s, but with regard to some parameter x of which fitness is an increasing function. For example, x might be some measure of enzyme activity. HARTL et al. 1985 considered the fitness function

(2)

where x is a quantity proportional to the catalytic efficiency of the enzyme. This relationship comes from metabolic control analysis (KACSER and BURNS 1973 ) and has the same form as the Michaelis-Menten equation, although it describes a different phenomenon (it does not relate reaction velocity to substrate concentration, and its horizontal asymptote does not correspond to saturation of an enzyme). This fitness function is plotted in Figure 1. Suppose that the distribution of x among mutants is constant for all values of x (I call such a parameter an equimutable parameter). The distribution of s will then change with x for two reasons. First of all, as x decreases, the fitness curve becomes steeper, leading to larger changes in fitness for a given change in x. Second, as fitness decreases with x, any change in fitness will constitute a larger fractional change in fitness, and hence a larger magnitude of s. Because we are interested in fractional fitness changes, a logarithmic scale for fitness is appropriate. Figure 2 illustrates graphically how, when the logarithm of fitness is a concave function of x, the distribution of selection coefficients broadens as fitness decreases. Because of the change in the distribution of s with x, a steady state can be reached at some point along this curve for a range of N_e. At this steady state, the steepness of the log-fitness function is such that the distribution of S leads to a balance between upward and downward substitutions (the average effect of substitution will be zero, but there will not necessarily be equal numbers of upward and downward substitutions). There will of course be stochastic variations about this steady state. A negative variation will lead to larger selection coefficients, and hence a compensatory upward trend. Conversely, an upward variation will lead to a downward trend. Thus, the steady state has a form of stability.

View larger version (10K): In this window In a new window Download PPT slide   Figure 1. The fitness function x/(1 + x). This function approaches a horizontal asymptote, which is indicated by the dashed line.

View larger version (20K): In this window In a new window Download PPT slide   Figure 2. The effect of a concave log-fitness function on the distribution of selection coefficients. For parental alleles of different fitness, identical distributions of mutational effects on the quantitative trait correspond to very different distributions of effects on the log of fitness. For the nearly neutral mutations of interest, effects on the log of fitness are nearly equal to selection coefficients. The point with lower fitness and a broader distribution of selection coefficients is a steady-state point for a smaller population than is the point with higher fitness.

Consider a sudden decrease in N_e for a population at such a steady state. Immediately, |S| for any given mutation will decrease, despite s remaining unchanged. Due to the narrower distribution of S, the value of x will tend to move downward. As x moves downward, the mutational distribution of s broadens. Eventually, if the new N_e is not too small, a point will be reached at which the broadening of the distribution of s will compensate for the smaller value of N_e to yield an equally broad distribution of N_es or S, and a steady state will be achieved for the new population size. Linear approximation of the fitness function at the equilibrium points implies that the distribution of S will be identical at steady state for the old and new population sizes (steady state is achieved when the slope of log-fitness is approximately proportional to 1/N_e, so the distribution of N_es is approximately invariant). Identity of the distribution of S implies that the same fraction of mutations will be fixed. If mutation rates are constant, this means equal rates of substitution regardless of population size, contrary to the view that substitution rate is a decreasing function of N_e. A change in N_e leads to a shift along the fitness curve that approximately compensates for the change in N_e by a change in the distribution of s.

A more formal development of the above is given in the Appendix 1 There it is shown that the requirement for the existence of a stable steady state for some range of N_e is that the logarithm of the fitness function is concave. It is also shown that concavity of the fitness function itself is a sufficient, but not a necessary, condition for concavity of the log-fitness function. It is noteworthy that for a simple linear fitness function, the logarithm of the fitness function is concave. A constant distribution of selection coefficients corresponds to an exponential fitness function.

This model bears formal similarity to certain models of gene interactions invoked to explain the cessation of Muller's ratchet. KONDRASHOV 1994 discussed synergistic epistasis as a force that slows the ratchet. In this context "epistasis" means deviation from multiplicativeness of fitness across loci, which corresponds to deviation from linearity of the log-fitness function. "Synergistic epistasis" corresponds to concavity of the log-fitness function. WAGNER and GABRIEL 1990 propose a role for compensatory (advantageous) mutations in halting Muller's ratchet in the absence of recombination. The assumptions of their model lead to an increasing sensitivity of the log of fitness to mutation as the genotype becomes less fit, which can be interpreted as epistasis.

AKASHI 1996 discussed a model in which the fitness function is concave and the fraction of deleterious mutations increases as a fitness optimum is approached. He noted that on this model the decline in fitness due to a decrease in population size would eventually cease at a new equilibrium, despite an excess of deleterious over advantageous mutations. Because it assumes that the ratio of deleterious to advantageous mutations changes with fitness, Akashi's model would predict a population size effect on substitution rate. However, this effect might be small, as discussed below.

	NUMERICAL RESULTS

TOP ABSTRACT THE MODEL NUMERICAL RESULTS DISCUSSION APPENDIX 1 LITERATURE CITED

Specification of an instance of the model:
A more concrete discussion will be based on the fitness function discussed above. However, rather than assuming that enzyme activity is equimutable, I think it more realistic to assume that the distribution of fractional change of activity is constant and, therefore, that the logarithm of activity is equimutable (one way to justify this assumption is to assume equimutability of a certain free energy of interaction, as catalytic efficiency is an exponential function of this free energy). This assumption is not an important part of the general model, which applies to a large class of fitness functions. It does illustrate that it may be a transformation of the parameter usually measured, rather than the parameter itself, that is equimutable. It also eliminates the possibility of fitness becoming negative.

For the purpose of discussion I assume a simplistic mutational spectrum in which x can either go up by a fixed amount or down by this same amount. This amount is chosen to be such that deleterious mutations correspond to 10% loss of activity. Each allele has an integral value x associated with it, and the interpretation of this value is that enzyme activity is proportional to (1/0.9)^x. Fitness, then, is given by

(3)

with x now interpreted as a logarithmic measure of enzyme activity. I further assume that 1% of all mutations increase x by one, with the remaining 99% decreasing it by one.

Several things can be done with this model and the results compared. First I calculate, on the basis of the linear approximation discussed above, the relative substitution rate at steady state. Second, using the actual fitness function, I obtain numerical solutions concerning the steady state for various population sizes. Finally, I present the results of simulations of the model, again for various population sizes.

Numerical calculations for the equilibrium point:
With linear approximation of the fitness function, the steady-state fraction of accepted mutations can be calculated on the basis of the distribution of x without regard to N_e or the fitness function (see Appendix 1). For the particular mutational distribution being discussed, the relative substitution rate is predicted to be ~0.0938.

With the actual fitness function in hand, one can compute the position of the steady state, along with the relative rate of substitution at this steady state, given a value of N_e. This value for the rate of substitution does not take into account the existence of stochastic variations around the steady state. The situation around a steady state is illustrated graphically in Figure 3, where the rates of upward and downward substitution are plotted as functions of x for a particular value of N_e. The steady state occurs where these two curves cross. Also plotted is the sum of these rates, which gives the overall relative rate of substitution as a function of x. It can be seen from this plot that small fluctuations around this steady state do not radically change the rate of substitution and, furthermore, that negative and positive fluctuations have somewhat compensatory effects on the overall rate. The steady-state values of the relative substitution rate for several effective population sizes, calculated numerically as explained in the Appendix 1 are shown in Table 1 and Figure 4. These values are not too different from the value 0.0938 obtained through linear approximation. Furthermore, the values for different N_e are even closer to one another, reflecting a certain type of similarity in the deviations from linearity at different points along the curve.

View larger version (14K): In this window In a new window Download PPT slide   Figure 3. The steady state for Ne = 100. The increasing and decreasing dashed curves represent the calculated rates of fixation of deleterious and advantageous mutations, respectively, as a fraction of the total mutation rate. The point at which these two curves cross is the steady-state point. The solid curve represents the total fixation rate and is the sum of the other two functions. Its value at the steady-state point gives the steady-state rate of substitution, relative to the mutation rate.

View larger version (11K): In this window In a new window Download PPT slide   Figure 4. Relative substitution rate, as a function of Ne, estimated in various ways. The lines represent calculations based on the linear approximation (dashed line) and the specific fitness function (solid line). The plotted points represent simulation results with mutation rate equal to 0.001 (crosses) and 0.01 (circles).

Computer simulations:
To check these results and to take account of realities such as fluctuations about the steady state, I have run simulations of this model. In the simulation, each haploid individual in the (n + 1)th generation is parented by an individual chosen from the nth generation. The parent is chosen randomly, with the probability of parenting being proportional to the relative fitness of the individual in generation n. Mutation is applied with some probability and with an effect on x in accordance with the distribution specified above. The population is an ideal one for which N_e = N. In each run of the simulation, the population is initialized with an allele value near its equilibrium point, as determined by prior simulation. After 30,000 generations, a count of substitutions is begun, and the simulation run for another 100,000 (µ = 0.01) or 1,000,000 (µ = 0.001) generations. The number of substitutions accumulated in these 100,000 or 1,000,000 generations is then determined for a single individual in the population. Repeated runs were performed for each population size.

The results of these simulations are shown in Table 1 and Figure 4. It can be seen that there is only a slight dependence of substitution rate on population size. The most important thing to note about this small dependence is that its direction is opposite to that predicted by OHTA 1977 and KIMURA 1979 . The increase in substitution rate with population size is presumably the result of increased polymorphism in the larger populations. Equation 1 applies to a newly arising mutant in an otherwise monomorphic population. In reality, the simultaneous existence of many alleles will lead to deviations from theory based on Equation 1. This explanation accounts for the increasing deviations of the simulations from the steady-state approximation as population size increases. It also accounts for the increase in the deviation from theory with higher mutation rate. A similar effect was seen by GILLESPIE 1994 in his simulations of certain models, and he gave it the same explanation. The effect is probably best understood as a dependence on the product of N_e and mutation rate. Real populations often have N_e much larger than those simulated here. However, they also have much smaller mutation rates.

GILLESPIE 1994 simulated a model with a concave log-fitness function (for sufficiently large values of the quantitative trait). However, he used a distribution of mutational effects that is symmetric about the parent allele value, rather than one for which the mean effect of mutation is negative. Presumably this is why he obtained the result that the rate of substitution was indistinguishable from the mutation rate, as is the case for strictly neutral mutations. While this result is technically a case of population size independence, it does not illustrate the fact that there can be N_e-independence even when the substitution rate is significantly lower than the mutation rate because of selection.

	DISCUSSION

TOP ABSTRACT THE MODEL NUMERICAL RESULTS DISCUSSION APPENDIX 1 LITERATURE CITED

Qualitative differences among allowable fitness functions:
The particular fitness functions discussed here are merely examples. Some ways in which fitness functions may differ qualitatively, and yet still be compatible with the proposed model, are noteworthy.

HARTL et al. 1985 , in discussing how Equation 2 can lead to the evolution of selective neutrality, emphasized that populations would be found on the plateau of fitness. While this may be true of real populations, the ideas developed here are equally applicable to populations whose steady-state points are far from the plateau for whatever reason. In fact, a concave function need not have a plateau, and the argument for N_e-independence holds for unbounded log-fitness functions, such as ln x.

For some fitness functions, such as the one used in the simulations described here (Equation 3), the derivative of the log-fitness function has an upper bound. This means that there must be a sufficiently large effective population size for a steady state to exist. In populations below this minimum size, which will depend on the mutational distribution, the gene will experience an indefinite downward slide in fitness, eventually being lost. The logarithms of other fitness functions, such as Equation 2, have derivatives that grow unbounded as x decreases.

The log-fitness functions specified here have been concave everywhere. A log-fitness function may be concave in only some places and yet still support stable steady states for a range of N_e. In fact a function that is alternately concave and convex might have multiple steady states for a given N_e.

The ratio of deleterious to advantageous substitutions:
The simplistic distribution of mutational effects discussed here has the property that the effects of advantageous and deleterious mutations are of identical fixed sizes. With such a mutational spectrum, the steady-state rates of advantageous and disadvantageous substitutions are necessarily equal. With other mutational spectra, however, this need not be the case. It might be that advantageous substitutions are rarer than deleterious ones, but are of correspondingly larger average effect, or vice versa. For advantageous mutations, the mean effect size must be greater among substitutions than among mutations because of the action of selection (assuming some variation in effect size). Similarly, the magnitude of the mean effect size of deleterious substitutions must be lower than that for mutations. If the mean effect size is the same for advantageous and deleterious mutations, then the mean effect size of advantageous substitutions will be greater than that for deleterious substitutions, and there will be correspondingly fewer advantageous than deleterious substitution events. While there is no good reason to make this assumption of equal means, different mutational distributions can lead to higher as well as lower fractions of deleterious substitutions. These considerations would seem to contradict GILLESPIE's (1994) claim that "there are no biologically realistic models in which most of the substitutions of mutations of very small effect are deleterious."

Population size, generation time, and the molecular clock:
The claim that substitution rate is a decreasing function of N_e has been invoked to explain the near constancy of the "molecular clock," measured in years, despite differences in generation time among lineages (OHTA 1972 , OHTA 1977 ). Longer generation times, it is argued, are found among organisms with larger body sizes and smaller population numbers. The increase in rate of substitution due to smaller N_e offsets the decreased number of generations per unit time. The model presented here, because it implies approximate independence from N_e, is incompatible with such an explanation. This fact might be taken as an argument against the model. However, it is uncertain whether this explanation of constancy is required. There is some evidence that there is a generation time effect, even for nonsynonymous substitutions (LI et al. 1987 ). Also, other explanations of constancy, such as equal numbers of cell generations per unit time (WILSON et al. 1977 ), have been put forth. Furthermore, this explanation, which depends upon the reciprocal nature of two unknown relationships, may create more problems than it solves. There are presumably many sets of species with roughly equal generation times but very different population sizes. The existence of a population size effect would imply different rates of substitution among such groups of organisms. In fact the estimates of NEI and GRAUR 1984 indicate a large amount of such variation for actual population size, although their estimates of effective population size show much less spread.

It would be enlightening to learn the rates of substitution in closely related organisms with similar lifestyles and generation times but radically different long-term population sizes. OHTA's (1972) model would predict a lower rate in the organisms with a larger population size, at least for nonsynonymous substitutions. The model proposed here would predict nearly identical rates, provided that the mutation rate per generation is the same in all lineages.

Claims of evidence for a population size effect:
Some experimental evidence has been claimed to support a population size effect on amino acid substitution rate. OHTA 1995 estimated the synonymous and nonsynonymous substitution rates for 49 genes in three mammalian lineages and compared the ratios of synonymous to nonsynonymous rates (d_S/d_N). She concluded that d_S/d_N is greater for the rodent lineage than the primate lineage and claimed that this difference is due to the larger population size of rodents. Ohta emphasized the existence of the difference in d_S/d_N, but did not discuss the size of the effect. Her numbers indicate that d_S/d_N for the rodent line is ~1.5 times that for the primate line. This is a rather modest effect, especially considering issues such as the potential pitfalls of correcting for multiple substitutions (NEI and GOJOBORI 1986 ). EASTEAL and COLLET 1994 present evidence that the apparently higher d_S/d_N for the "rodent lineage" is actually due to a high value for an interior branch of the tree that is not on the line leading to rodents. Although their study was based on a different (and smaller) set of genes, it casts serious doubt on a population size effect as an explanation of Ohta's results.

MORAN 1996 reports a larger difference in d_S/d_N between lineages, also invoking a difference in N_e to explain the effect. In this case endosymbiotic bacteria (Buchnera) are compared to enteric bacteria (Escherichia coli and Salmonella typhimurium). The endosymbionts are found to have smaller d_S/d_N for the genes encoding the enzymes of tryptophan biosynthesis. It is interesting to compare these numbers to those for mammalian genes. The d_S/d_N ratios for the 80 enzymes among WOLFE and SHARP's (1993) collection of rodent genes (median = 13.0) are comparable to those for the enteric trp genes (median = 17.5). (Mammalian genes for nonenzymes often have lower, but sometimes much higher, d_S/d_N.) The values for the endosymbionts stand out as exceptionally low (median = 2.2), despite the fact that N_e is likely larger for them than for mammalian species. The similarity of values between the enteric bacteria and mammals might be taken as evidence against a strong population size effect, and the lower values for Buchnera are considered to be the result of some other factor.

The comparison of rates in enteric bacteria and mammals need not be confined to these few enteric genes. WOLFE and SHARP 1993 compared their estimates of rates for rodent genes to estimates of enteric rates (SHARP 1991 ). The data indicate that d_S/d_N is higher for enterics than for rodents, but only by a factor of 3.4. This ratio is tiny compared to the difference between the effective population sizes of enteric bacteria and rodents. It is also much too small to compensate for the difference in generation times between these groups. Furthermore, as Wolfe and Sharp point out, the difference in d_S/d_N may be due to the inclusion of different types of genes in the two data sets. Most notably, the enteric data set contains mostly enzymes. The enzymes in the rodent data set have higher d_S/d_N than nonenzymes. Consideration of only the enzymatic subset of the rodent data reduces the ratio of d_S/d_N values to 2.3.

Insensitivity of the substitution rate to the fitness function:
If the proposed model is a mixed bag for constancy of substitution rates across lineages, perhaps it provides an explanation for another kind of constancy. It is sometimes remarked that the amino acid substitution rate is rather variable among different proteins, especially by comparison to the rate of silent substitution. From one point of view, though, this variation is surprisingly small. LI and GRAUR 1991 presented estimates of substitution rates for 36 mammalian genes. Excluding some exceptionally slowly changing genes such as histones and actin, the nonsynonymous substitution rates of most of the genes listed (30/36) were within an order of magnitude of 10^-9 per site per year. Different proteins are of radically different importance to fitness, some being essential to life, while others are subtle adaptive refinements. Even among essential proteins, one might expect large variations in the sensitivity of fitness to mutation. From the perspective that stronger selection leads to lower substitution rate, one might expect tremendous variation in substitution rate among proteins. From this point of view, the clustering of substitution rates around 10^-9 per site per year is remarkable.

The proposed model can account for this lack of variation. In the context of this model, differences in "strength of selection" might be interpreted as differences in fitness function. By an argument based on Equation A1, and similar to that for N_e-independence, it can be shown that the substitution rate is approximately independent of fitness function. This is so because the steady state will be reached at approximately the same slope of log-fitness, regardless of the choice of fitness function. The differences that do exist among proteins would be explained by differences in the distribution of x among mutants, including differences in the ratio of favorable to unfavorable mutations.

Of course another explanation of the similarity of substitution rates is that the majority of mutations are either extremely deleterious or extremely close to neutrality. However, such a model also implies N_e-independence. Indeed it is exactly such a model that the "nearly neutral theory" was intended to replace (OHTA 1972 ).

The insensitivity of the substitution rate to changes in the fitness function may help to explain rate constancy across lineages. The selective importance of a given gene is likely to vary among organisms. In a conventional model, this variation would be reflected in differences in the distribution of selection coefficients of mutants. Such differences would have the same effect as population size differences on the all-important product N_es and would be a source of substitution rate variation among lineages for the gene in question. In contrast, the model that I propose would predict no such variation in the steady-state rate among lineages.

Variance of the substitution process:
Neutralist theories, in their simplest forms, predict that substitution is a Poisson process and hence has a variance equal to its mean. Neutral theories have been criticized on the grounds that amino acid sequence data indicate a variance much greater than the mean (GILLESPIE 1984 , GILLESPIE 1991 ). As shown in Table 1, the simulations yielded estimates of index of dispersion (ratio of variance to mean) that are somewhat higher than one. This is not surprising because substitutions are not independent. Rather, a substitution event changes the predominant fitness value and hence the future rate of substitution. This inflation of variance is, in GILLESPIE's (1991) terminology, a "residual" rather than a "lineage" effect. While the estimates of index of dispersion from the simulations are not as high as those observed for real proteins, other combinations of fitness functions and mutational distributions might yield higher values.

Furthermore, many neutralist explanations of the overdispersed clock are in some way compatible with the proposed model. TAKAHATA 1987 discussed how a fluctuating population size could lead to overdispersion of the clock, because bursts of substitution would occur when the population is small. The present model does not predict higher steady-state substitution rates in small populations, but it does predict bursts of substitution when the population size is suddenly changed. These bursts can occur when the population size is increased as well as decreased and might account for variance among lineages with different population size histories. TAKAHATA 1987 also discussed the "fluctuating neutral space" model, in which successive substitutions affect the substitution rate by changing the level of selective constraint. This model of Takahata's predicts a residual variance effect. It shares with the model presented here a lack of constancy of the distribution of selection coefficients. However, Takahata seemed to be concerned with differences in the distribution of selection coefficients even among alleles of the same fitness. It would be unrealistic to deny the existence of such differences. In the context of the model presented here, one could allow for differences in the mutational spectrum of x among parent alleles with the same value of x and still maintain that the value of x for an allele is not predictive of its mutational spectrum. Such a state of affairs might be described as "weak equimutability," with "strong equimutability" implying the additional (and unrealistic) constraint of equal distributions of x for all alleles.

Plausibility of the model:
The question of how realistic the proposed model is may be divided into two parts. First, there is the question of how well a model in which fitness is a function of a single variable captures reality. Second, the assumption that there is an equimutable parameter should be examined.

The notion that fitness is a fixed function of only enzyme activity (as defined above) is simplistic. This model is meant to describe the effects of changes in catalytic efficiency due to changes in amino acid sequence. Changes to the level of protein production, due to regulatory mutations or to regulation itself, are not well modeled. Such changes will have an effect on fitness due in part to the cost of protein production, so that enzyme activity is not the sole determinant of fitness. Furthermore, Equation 2 is based on the assumption that the activities of other enzymes in a biochemical pathway are fixed. In reality, these activities too will change mutationally. Complications similar to these may exist for quantitative traits other than enzyme activity. These complications raise difficult issues such as what it is that ultimately limits fitness. Despite such complications, concave log-fitness functions may have great value as models of reality. The model that I am proposing is no more simplistic, in these and other regards, than many other models for the distribution of selection coefficients, including those used to derive effects of population size on substitution rate.

The argument presented here for an approximately constant substitution rate regardless of N_e depends upon the existence of an equimutable parameter. Of course the distribution of mutational effects on x may be systematically different for parent alleles with different values of x. In particular, less fit alleles may, as discussed above, have more of a tendency to be improved, and less to be harmed, by mutation. While in some cases there will exist a transformation that yields an approximately equimutable parameter, many types of nonconstancy cannot be "transformed away." For example, the appropriate type of transformation will not change the ratio of favorable to unfavorable mutations, so if this quantity changes, as in the model discussed by AKASHI 1996 , there can be no equimutable parameter. It seems plausible that deviations from equimutability are small over the range of x values needed to accommodate a range of N_e. In any case, concave log-fitness functions, which merely reflect a kind of diminishing returns, may be sufficiently common that concavity usually contributes to the attainment of a steady state, even when deviations from equimutability also contribute. In such a situation, the population size effect will at least be mitigated by the contribution of the changing slope of the log-fitness function. Genes might differ greatly with regard to the relative contributions of the two effects and therefore with regard to the degree of population size effect. In the limit as concavity becomes the dominating force in the attainment of a steady state, the population size effect disappears.

	ACKNOWLEDGMENTS

I thank Jon Seger for discussions and advice. This work was supported in part by National Institutes of Health grant 5 P50 HG-00199-07.

Manuscript received November 21, 1997; Accepted for publication June 10, 1998.

	APPENDIX 1

TOP ABSTRACT THE MODEL NUMERICAL RESULTS DISCUSSION APPENDIX 1 LITERATURE CITED

The model and its consequences:
Consider a locus that exerts its effect on fitness solely through a parameter x that it controls. Further assume that x is equimutable, as defined above. For small x, the change in fitness can be approximated by w = w'(x)x. However, we are interested in s, which is the ratio w/w, so the important relationship is

(A1)

The distribution of s, according to this approximation, will vary with x only by scaling. The requirement for the existence of a stable steady state (for at least some range of N_e) is that the distribution of s broadens as x decreases. Equation A1 shows that what is required is concavity of ln w(x). It is easily shown (see below) that concavity of w(x) is a sufficient, but not a necessary, condition for concavity of ln w(x), where w(x) is positive. The steady states for various values of N_e are achieved at values of w'(x)/w(x), such that w'(x)/w(x)1/N_e, and hence the distribution of

is approximately constant, given that the distribution of x is constant. Thus, because the fraction of accepted mutations depends only on the distribution of S, the rate of substitution is independent of N_e at steady state.

Concavity of fitness and log-fitness functions:
Suppose that f(x) is a concave function (f''(x) < 0) and is positive (as a fitness function must be, unless it is zero). Suppose that g(x) = ln f(x). Then

and

The squared terms must be nonnegative, and the concavity and positiveness of f guarantee that f(x)f''(x) is negative, so g''(x) < 0. Thus, concavity of the fitness function implies concavity of the log-fitness function.

For the nonnecessity of the concavity of f for that of its logarithm, consider such examples as f(x) = x³, a convex function for which ln f(x) = 3 ln x is concave, or even f(x) = x, a function with zero second derivative for which ln f(x) = ln x is concave. With regard to the latter function, note that a constant distribution of s among mutants does not correspond to a linear fitness function, but rather to an exponential fitness function and hence a linear log-fitness function.

Predictions of steady states and substitution rates:
The condition for an equilibrium is that the mean effect of substitutions on the quantitative trait is zero (the concavity condition is needed for the stability of this equilibrium). For the particular model used in the simulations, the requirement is that 0.01 () = 0.99 () , where S_d and S_u are the values of S for down and up mutations, respectively. The left- and right-hand sides of this equation are the rates, relative to the total mutation rate, of favorable and unfavorable substitutions, respectively. The total relative rate of substitution at any point, including the equilibrium, is given by the sum of the upward and downward rates, namely 0.01 () + 0.99 () .

Under linear approximation, S_d = -S_u, so letting S = S_u the equilibrium condition becomes 0.01 () = 0.99 () . The solution of this equation is S = ln 99 (more generally for this type of model, S = ln r, where r is the ratio of the rate of down to up mutation). The relative rate of substitution for this equilibrium value of S is ~0.0938.

Given the fitness function (in this case Equation 3), it is possible to do a similar calculation that uses actual values for selection coefficients, rather than relying on a linear approximation. Again, the equilibrium condition is 0.01 () = 0.99 () , but now selection coefficients are obtained as functions of the parental x value, which involve the fitness function. The selection coefficient for an advantageous mutation is given by and that for a deleterious one by . For any specified value of N_e, S_d and S_u are expressed as functions of x, and the resulting equilibrium equation is solved for x (here done numerically using Matlab's fsolve function). The values S_d and S_u for this x are then used to calculate the steady-state relative substitution rate for the specified N_e. This process is illustrated graphically in Figure 3.

	LITERATURE CITED

TOP ABSTRACT THE MODEL NUMERICAL RESULTS DISCUSSION APPENDIX 1 LITERATURE CITED

AKASHI, H., 1996  Molecular evolution between Drosophila melanogaster and D. simulans: reduced codon bias, faster rates of amino acid substitution, and larger proteins in D. melanogaster. Genetics 144:1297-1307[Abstract].

EASTEAL, S. and C. COLLET, 1994  Consistent variation in amino-acid substitution rate, despite uniformity of mutation rate: protein evolution in mammals is not neutral. Mol. Biol. Evol. 11:643-647[Abstract].

GILLESPIE, J. H., 1984  Molecular evolution over the mutational landscape. Evolution 38:1116-1129.

GILLESPIE, J. H., 1991 The Causes of Molecular Evolution. Oxford University Press, New York.

GILLESPIE, J. H., 1994  Substitution processes in molecular evolution. III. Deleterious alleles. Genetics 138:943-952[Abstract].

HARTL, D. L., D. E. DYKHUIZEN, and A. M. DEAN, 1985  Limits of adaptation: the evolution of selective neutrality. Genetics 111:655-675[Abstract/Free Full Text].

KACSER, H. and J. A. BURNS, 1973  The control of flux. Symp. Soc. Exp. Biol. 32:65-104.

KIMURA, M., 1962  On the probability of fixation of mutant genes in a population. Genetics 47:713-719[Free Full Text].

KIMURA, M., 1979  Model of effectively neutral mutations in which selective constraint is incorporated. Proc. Natl. Acad. Sci. USA 76:3440-3444[Abstract/Free Full Text].

KONDRASHOV, A. S., 1994  Muller's ratchet under epistatic selection. Genetics 136:1469-1473[Abstract].

LI, W.-H., and D. GRAUR, 1991 Fundamentals of Molecular Evolution. Sinauer Associates, Sunderland, MA.

LI, W.-H., M. TANIMURA, and P. M. SHARP, 1987  An evaluation of the molecular clock hypothesis using mammalian DNA sequences. J. Mol. Evol. 25:330-342[Medline].

MORAN, N. A., 1996  Accelerated evolution and Muller's rachet in endosymbiotic bacteria. Proc. Natl. Acad. Sci. USA 93:2873-2878[Abstract/Free Full Text].

NEI, M. and T. GOJOBORI, 1986  Simple methods for estimating the numbers of synonymous and nonsynonymous nucleotide substitutions. Mol. Biol. Evol. 3(5):418-426[Abstract].

NEI, M. and D. GRAUR, 1984  Extent of protein polymorphism and the neutral mutation theory. Evol. Biol. 17:73-118.

OHTA, T., 1972  Evolutionary rate of cistrons and DNA divergence. J. Mol. Evol. 1:150-157.

OHTA, T., 1977 Extension of the neutral mutation drift hypothesis, pp. 148–167 in Molecular Evolution and Polymorphism, edited by M. KIMURA. National Institute of Genetics, Mishima, Japan.

OHTA, T., 1995  Synonymous and nonsynonymous substitutions in mammalian genes and the nearly neutral theory. J. Mol. Evol. 40:56-63[Medline].

OHTA, T. and H. TACHIDA, 1990  Theoretical study of near neutrality. I. Heterozygosity and rate of mutant substitution. Genetics 126:219-229[Abstract].

SHARP, P. M., 1991  Determinants of DNA sequence divergence between Escherichia coli and Salmonella typhimurium: codon usage, map position, and concerted evolution. J. Mol. Evol. 33:23-33[Medline].

TAKAHATA, N., 1987  On the overdispersed molecular clock. Genetics 116:169-179[Abstract/Free Full Text].

WAGNER, G. P. and W. GABRIEL, 1990  Quantitative variation in finite parthenogenetic populations: What stops Muller's ratchet in the absence of recombination? Evolution 44:715-731.

WILSON, A. C., S. S. CARLSON, and T. J. WHITE, 1977  Biochemical evolution. Annu. Rev. Biochem. 46:573-639[Medline].

WOLFE, K. H. and P. M. SHARP, 1993  Mammalian gene evolution: nucleotide sequence divergence between mouse and rat. J. Mol. Evol. 37:441-456[Medline].