becomes sufficiently large, and that this might confuse the analysis. Therefore they pointed out that it is important to identify an adequate time scale for discriminating the effect of natural selection from that of random genetic drift of selectively neutral variations. In order to address this problem, it is crucial to construct a neutral model which can predict the long-run behavior of evolutionary changes in mRNA abundance. A neutral model of gene expression evolution was first proposed by Khaitovich and coworkers (2004) [22]. They constructed this model based on the observation that gene expression divergence increases proportionally with divergence time in lineages of primates and rodents, which is termed as a ‘‘clock-like’’ accumu- lation of gene expression divergence. This observation can be explained from the assumption that mutations cause changes in relative amounts of expression levels irrespective of gene function [22,26]. However, since those studies were confined mainly with relatively short terms of gene expression evolution, such as between humans and chimpanzees, the long-run behavior of the neutral model of gene expression evolution has not been well studied. To investigate long-run behavior, in this study we focused on a property of the distribution of mRNA abundance. As soon as genome-wide gene expression profiling techniques were devel- oped, it was revealed that there is a common tendency in the distribution of mRNA abundance: a few genes are expressed intensely and most genes are expressed at quite low levels. It is now known that this distribution can be well described by Zipf’s law [31] (or its mathematical equivalent, called the power law) from vertebrates to lower eukaryotes [32–37]. This law states that there is a relationship between the mRNA abundance (f [copy/cell]) and its abundance rank (r) represented by f !r{b. The exponent b is the absolute value of the slope in a logarithmic rank-frequency plot. It is remarkable that the value of exponent b is near 1.0 in most normal tissues. (More specifically, the exponent b is near 1.0 in normal tissues composed primarily of a homogeneous population of differentiated cells, such as liver and muscle. In other normal tissues composed of a mixture of different types of cells, including brain, testis, and kidney, the exponent b tends to be slightly lower than 1.0, as expected [35]). Here it should be noted that if evolutionary forces affect the mRNA abundance of each gene, they inevitably affect the distribution of it. Moreover, the distribution of mRNA abundance can be expected to be generated by the long-term effect of evolutionary forces on gene expression regulation. Therefore, if most evolutionary changes in the mRNA abundance can be explained by the neutral hypothesis, then the above two mentioned phenomena, namely the clock-like accumulation of gene expression divergence and Zipf’s law of the transcriptome, should be explained with the same neutral model. Indeed, Ogasawara et al. (2003) [35] suggested that Zipf’s law of the transcriptome might originate from a process which is quite similar to the neutral model introduced by Khaitovich et al. (2004) [22]. However, as pointed out in this study, the previous model cannot explain the uniformity of the exponent b, which suggests that long-run behavior of previous models is not consistent with observations. In this study we propose refined models which can explain the two phenomena simultaneously by adding a few cytological constraints to the previous models. Results A Neutral Model of Gene Expression Evolution and the Genesis of Zipf’s Law-Like Distribution in mRNA Abundance Consider the mRNA abundance fi(t) (number of mRNA molecules in a cell) of gene i in a given cell type at generation t. Gene Expression Evolution We regard fi(t) as a heritable quantitative trait whose value is determined solely by the genetic effect, and we do not take account of the fluctuation caused by physiological and environmental factors. We assume that mRNA abundance fi(t) is affected by mutations with proportionality to the abundance before mutation. This can be expressed as follows: fiðtÞ~ð1zkitÞfiðt{1Þ where kit is a mutually independent and identically distributed (iid) random variable with mean 0 and variance s2, and is also statistically independent of fi(t). This assumption means that the probability of evolutionary change, from say 1.0 copy/cell to 2.0 copies/cell, is equal to the probability of change from 10 copies/cell to 20 copies/cell, and not to the probability from 10 copies/cell to 11 copies/cell. The neutral model of gene expression evolution introduced by Khaitovich et al. (2004) [22] is also based on essentially the same idea. In this study, we refer to the model assuming that evolutionary changes in mRNA abundance are irrelevant to the function of the proteins as neutral model of gene expression evolution. On the other hand, we refer to the model assuming that the range of evolutionary changes is confined depending on the function of the protein as the natural selection model of gene expression evolution. We use the terminology from previous studies [22,26], but it should be noted that the neutral model of gene expression evolution and Kimura’s neutral model of molecular evolution are different concepts. Kimura’s neutral theory asserts that the majority of DNA sequence polymorphisms observed within a species has no effect on the fitness, because the majority of harmful mutations are eliminated from the population by negative selection. Therefore it is obvious that both the natural selection model and the neutral model of gene expression evolution do not conflict with Kimura’s neutral model of molecular evolution. Assuming that the absolute value of kit is small compared with 1, we approximate from the Taylor expansion of lnð1zxÞ~xzO x2 that ln fi(t)~ln fið0Þz t X kij j~1 Since we assume that kij (j = 1, 2,…, t) is iid random variables, according to the central limit theorem, P kij will be approxi- mately a normal distribution with mean 0 and variance sk2 = s 2 t when t??. Therefore, ln fiðtÞ is asymptotically normally distributed and hence fi(t) is lognormally distributed. In general, a positive random variable Z is said to be lognormally distributed with two parameters, mean m and variance s2, if X ~lnZ is normally distributed with them. It is known that when s2 is sufficiently large, the lognormal distribution would appear almost linear on both a log-log plot of the probability density function (power-law plot) and on a Zipfian plot (Figure 1). This is the model of the genesis of Zipf’s law of the transcriptome suggested by Ogasawara [35], and this process is generally known as the Gibrat process or the law of proportionate effect which was first introduced as an explanation of firm size distribution [38–41]. However, since the variance of the resultant distribution of ln fiðtÞ is given by s2t, it increases with generation t without limit. Consequently, the exponent b also increases with t without limit (Figure 1 and Figure S1). Obviously this situation is unrealistic because there is no reason for many species to have an appropriate PLoS ONE | www.plosone.org 2 November 2009 | Volume 4 | Issue 11 | e7943