Gene Expression Evolution Figure 1. Log (rank) versus log (frequency) plot for mRNA abundance. (A) The mRNA abundance of human (red), mouse (green), chicken (blue), and X. laevis (light blue). (B) Result of the Monte Carlo simulation, L = 0.0. (C) Result of the Monte Carlo simulation, L = 1.0. At the initial state, the mRNA abundance of all genes is set to M/N. Other model parameters of (B) and (C) are: M = 300,000, N = 20,000, and sk = 0.05. The line shows y = 0.1/x doi:10.1371/journal.pone.0007943.g001 value of t to make b be near 1.0 concurrently. Therefore this model is inconsistent with observations. We refined this previous model by adding a few cytological boundary conditions. First, we assume that the total number of mRNA molecules (M) in a given cell is constant throughout a specific evolutionary process. Secondly, the number of expressed genes (N) in a given cell type is constant throughout this evolutionary process. In typical mammalian cells, M and N are estimated to be 300,000–500,000 and 10,000–30,000, respectively [42–44]. This assumption implies that the mRNA abundance in a cell is determined by the relative, instead of absolute, affinity of the gene regulation proteins and the cis-elements among the genome, since those proteins would ‘‘choose’’ their target from a collection of cis-elements in the gene expression process. This corresponds to the normalization step in this model (see Materials and Methods). It should be noted that by assuming the maximum capacity of mRNA in a cell (M), a limitless increase in the variance of the mRNA abundance by the Gibrat process leads most genes to decrease their mRNA abundances without limit. That is, most genes will lose their gene expression ability eventually (Figure 1 and Figure S1). Obviously this seems unrealistic and contradicts the assumption of the constancy of N. To avoid this, we make a third assumption that an individual who loses gene expression of functional genes will be eliminated from the population by negative selection. We formulate this assumption by prohibiting the mRNA abundance of each gene from being lower than a certain limit L; in the formula, fiðtÞ§L for any i and t (Materials and Methods). Some genes may have a different lower limit L from other genes according to their function; however, we ignore these as minor cases for the purpose of modeling gene expression evolution. We next investigated the long-run behavior of the refined model by a series of Monte Carlo simulations. The simulations indicate that the process approaches the steady-state, Zipf’s law-like distribution regardless of the initial state distribution when model parameters are set to the typical value of mammalian cells (M = 300,000, N = 20,000), and the lower limit of expression level L is set to 1.0 [copy/cell] (Figure 1C). The consistency of the model (with L = 1.0, M = 300,000, and N = 10,000) and data has been tested by the generalized chi-square goodness-of-fit test [45]. We use the EST-based gene expression profiles of livers of human, mouse, chicken, and Xenopus laevis, since liver tissue is composed primarily of a homogeneous population of differentiated cells. Estimating the discrepancy measure d from human liver EST profiles (d = 0.0268), the P-values of the generalized goodness-of-fit tests of human, mouse, chicken, and X. laevis were 1.0 (n = 2,384), 0.997 (n = 762), 0.794 (n = 935), and 0.999 (n = 1304), respectively (df = 4). Even with a rather conservative assumption that the relative error in the EST profiles is only 3% for each bin (d = 9.061024), P-values of human, mouse, chicken, and X. laevis were 0.68, 0.46, 0.0025, and 0.10, respectively (see Materials and Methods). In conclusion, the model is not rejected by the data. The relatively small P-value of the chicken data in the latter case can be readily explained by the measurement error of EST-based profiling. We also show by simulation that the standard deviation of the random coefficient kit has little effect on the form of the steady- state distribution, as expected from the central limit theorem, but that it affects the evolutionary rate. That is, the number of generations required for convergence increases as the standard deviation of kit decreases (Figure S2). This result implies that the exponent b is solely determined by the model parameters L, M, and N. Therefore, we examined the relationship between the values of L, M, N, and the resultant exponent b of the steady-state distribution (Figure 2). The result of the Monte Carlo simulation shows that exponent b is much more susceptible to L but hardly PLoS ONE | www.plosone.org 3 November 2009 | Volume 4 | Issue 11 | e7943