Gene Expression Evolution Figure 5. Distribution of mRNA abundance under the optimum mRNA abundance model. (A) Optimum value distribution dassumed in the Monte Carlo simulation. (B–F) Steady-state distributions of mRNA abundance where n genes are subject to natural selection. Red points represent genes undergoing natural selection. Blue points are neutral genes. Model parameters are: L = 1.0, M = 300,000, N = 20,000, and sk = 0.01. Population size (Ne) is 100. t = 50,000. The selection intensity is 1.67. doi:10.1371/journal.pone.0007943.g005 A series of Monte Carlo simulations demonstrated that this model can explain the following three phenomena simultaneously: (i) the clock-like accumulation of gene expression divergence reported by Khaitovich (2004) [22]; (ii) the diminished upper limit of the average pairwise variance reported by Bedford (2009) [55]; (iii) the Zipf’s law-like distribution of mRNA abundance in a wide variety of species. Although the effect of cytological constraints on the time development of gene expression divergence has been predicted by Lemos et al (2005) [23], an explicit formulation of cytological constraints has not been given in previous works. The present work introduces a formulation of cytological constraints (assump- tions 2 and 3) and a scheme of natural selection optimizing the mRNA abundance (assumption 4), based on the relationship between the evolutionary change in the mRNA abundance and the distribution of mRNA abundance in a cell. One feature of this model is that locus-specific gene expression divergence should not increase linearly with generation number; on the contrary, it behaves like a random walk. Therefore the distribution of locus-specific gene expression divergence should exhibit a rather different and complex shape depending on the divergence time. On the other hand, Khaitovich’s gene expression divergence was defined by the average value over all loci comprising the transcriptome, instead of being locus-specific. In this case, the divergence should increase linearly with generation time; that is, a clock-like accumulation of gene expression divergence. This is an essential difference between our model and an infinite site model of neutral polymorphisms of nucleotide sequences in which locus-specific divergence accumulates propor- tionally with generation time. Although our model currently focuses on describing the evolutionary process of a single cell for simplicity, this model might also be applicable for describing the gene expression evolution of multicellular organisms. The results of genome-wide gene expression profiling of higher eukaryotes, including humans, indicate that most genes can be clearly classified into two well- PLoS ONE | www.plosone.org 7 November 2009 | Volume 4 | Issue 11 | e7943