Models of Gene Expression Evolution with Natural Selection In the previous sections, we discussed neutral models of gene expression evolution in which evolutionary changes in mRNA abundances are assumed to be irrelevant to the functions of the proteins. However, since it is evident that adequate regulation of the mRNA abundance is a prerequisite for efficient protein function, it is necessary to extend the model to account for the effect of natural selection which optimizes the mRNA abundance. As we mentioned previously, the distributions of mRNA abundance in cells are similar in a variety of species, including mammals, insects, plants, and lower eukaryotes. It seems that this distribution is not easily explained from a natural selection model of gene expression evolution, because these species having different sets of genes are adapted to different environments and different ways of life. If those genes have their own optimum mRNA abundance (mi) according to their functions, then the observed distribution of mRNA abundance in cells would differ among the various species. There is no reason that each of the optimum values {mi} (i = 1,…,N) follows Zipf’s law-like distribu- tion. Rather, {mi} can take any shape, even one far from the Zipf’s law-like distribution depending on the functions of the proteins, since distantly related species, such as insects and plants, have different sets of genes and are adapted to different environments with different ways of life. The problem here is how these two seemingly contradictory phenomena can be explained simulta- neously from a unified model. First, we examined whether the steady-state distribution can be altered, assuming the optimum mRNA abundance for each gene. Since Zipf’s law-like distributions are observed in a variety of phenomena, one may expect that this type of distribution is so robust that it is not altered by additional factors such as natural selection. To examine this hypothesis, we undertook the Monte Carlo simulation with random mating populations with finite population size comprised of sexually reproducing diploid organisms. We assume that the mRNA abundance of each locus is determined by the average effect of two genes on the locus, and that fitness of each individual is determined by the mRNA abundance of a given set of genes. At present, there is insufficient information on ‘‘the distribution of the optimum mRNA abundances’’ for simulation. We can only conjecture that there is no reason for a variety of species to have an identical optimum value distribution. Therefore, we assume that the optimum value distribution has a slightly distorted shape from the Zipf’s law-like distribution with the optimum value distribution given as 10X+1, where X is a random variable having a skewed normal distribution (Figure 5). Then we examined whether the steady-state distribution converges to the Zipf’s law-like distribution. The Monte Carlo simulation demonstrated that the steady-state distribution is readily altered by the effect of natural selection, although only when a small number of genes (about 0–100) is assumed to be at their optimum did the process converge into a Zipf’s law-like distribution (Figure 5). It should be noted that in these simulations the selection intensity (defined as the reciprocal of the proportion of individuals selected from a population to be used as parents) was the same across different values of n, the number of loci undergoing the natural selection (Materials and Methods). This result showed that the Zipf’s law-like distribution was not robust against the effect of natural selection under the absolute optimum value model. This suggests that it is not easy to explain the universality of the Zipf’s law-like distribution under the absolute optimum value model without assuming that the optimum value distribution itself follows Zipf’s law. However, Gene Expression Evolution there is no reason to think that this additional assumption holds. Rather, if there exists an optimum value distribution, it seems natural that distantly related species would have different optimum value distributions as mentioned above. Bedford and Hartl (2009) [55] reported that, in several Drosophila species, average pairwise variance in the mRNA abundance as they defined it was saturated at around 0.3. They pointed out that the variance is expected to equal 1.0 under the neutral model of gene expression. Indeed, the expectation of our refined neutral model is also 1.0. For an explanation of this observation, Bedford and Hartl (2009) [55] proposed an evolutionary model based on the Ornstein-Uhlenbeck (OU) process described as following a stochastic differential equation: dxiðtÞ~ðmi{xiðtÞÞlidtzdidW ðtÞ, where, xi(t) = log2fi(t) and fi(t) is the mRNA abundance of gene i at time t, and mi is the log of optimum mRNA abundance of gene i. W is the standard Brownian motion (also called a Wiener process), and li and di are positive constants corresponding to the strength of natural selection and random genetic drift acting on gene i, respectively. In other words, Bedford and Hartl argued that evolutionary change in the mRNA abundance can be described by the balance between the drift of Brownian motion and natural selection pulling the mRNA abundance toward its optimum. However, since the Bedford-Hartl model assumed an optimum value mi (i = 1, …, N) of the mRNA abundance for every gene, it is not easy to explain the Zipf’s law-like distribution of mRNA abundance in cells from this model. In our simulation scheme, a selection intensity of 1.67 is necessary for constraining the Bedford’s average pairwise variance to around 0.3; however, the steady-state distribution of mRNA abundance departs from the Zipf’s law-like distribution under this selection intensity (Figure 6). In order to reconcile this apparent inconsistency, we introduce another model of gene expression evolution. In this model, we assume that the genes have no intrinsic optimum mRNA abundance such as mi. Instead, the mRNA abundance suitable for their function is determined by their relative relationship to the other genes, such as fi.fj. Such relationships would reflect the whole structure of gene regulatory pathways, and would thus be as complex. Here we assume that the relationships are f1.f2.….fn, for simplicity. The Monte Carlo simulation demonstrated that this model can explain both the Bedford’s average pairwise variance and the Zipf’s law-like distribution simultaneously (Materials and Methods, Figure 7). The behavior of the relative abundance model can be summarized as follows: the distribution of mRNA abundance in a cell is determined both by random genetic drift and by cytological constraints, and the rank of mRNA abundance of each gene is optimized by natural selection. Discussion In this study, we introduce a novel, explicitly defined model of gene expression evolution. Our model assumes that (1) the mRNA abundance is affected by mutations proportional to the abundance before mutation; (2) the number of expressed genes and total number of mRNA molecules in a cell type is nearly constant throughout a specific evolutionary process; (3) the mRNA abundance of each gene has a lower limit near 1.0 copy/cell; and (4) the mRNA abundance of each gene (fi) suitable for its function is determined by its relative relationship to other genes, i.e., fi.fj. PLoS ONE | www.plosone.org 6 November 2009 | Volume 4 | Issue 11 | e7943