Several other models have been proposed to explain the Zipf’s law-like distribution of mRNA abundance. Without exception, those models were based on the dynamics of mRNA synthesis and degradation in a cell, rather than evolutionary processes. Furusawa and Kaneko proposed a model which asserts that the distribution originates from the balance between the effect of up- regulating genes and down-regulating genes in the intracellular gene regulation network [33]. Kuznetosov assumed that the dynamics can be described by the birth-death stochastic process [34]. Similarly, Ueda et al. formulated this process as geometrical Brownian motion [37]. According to the formulation of any of those dynamic models, it follows that the expression level of every gene would fluctuate within the full range, namely from 0.0 copy/ cell to the maximum value permitted, virtually at random. To avoid this consequence, those models must make an additional assumption that each gene has its own dynamic range of expression levels. This course of thinking would lead to the basic idea of the evolutionary model, where the gene expression level is determined by its genome. Materials and Methods Estimation of mRNA Abundance Based on EST Frequencies We downloaded the UniGene data files of Homo sapiens (Build #210), Mus musculus (Build #170), Gallus gallus (Build #39), and X. laevis (Build #82) from the NCBI FTP site. With reference to the property table of cDNA libraries provided by the BodyMap-Xs database (http://bodymap.jp) [68], we selected cDNA libraries having the largest number of EST sequences among the libraries which were originated from pathologically normal tissues of adult organisms, and which did not suffer from experimental procedures such as normalization, subtraction, selection, and full-length cDNA enrichment. We estimated the mRNA abundance by dividing the number of EST sequences belonging to each UniGene cluster in each cDNA library by the total number of EST sequences in the cDNA library. The list of the dbEST library IDs we used in the analysis is as follows: liver 6989 (human), 2484 (mouse), 11,222 (chicken), 5540 (Xenopus); brain 18318 (human), 12634 (mouse), 15560 (chicken), 8910 (Xenopus); testis 18476 (human), 11976 (mouse), 15563 (chicken), 12882 (Xenopus); kidney 18374 (human), 7268 (mouse), 11220 (chicken), 11985 (Xenopus). Monte Carlo Simulation (1) Simple Stochastic Process For the examination of the convergence to the Zipf’s law-like distribution, we assumed that the mRNA abundance at the initial state (t = 0) was fi(0) = M/N for each gene (i = 1,2, …, N). At the transition to the next generation, the mRNA abundance of each gene was altered by mutation following fi’(t+1) = fi(t) (1+kit), where kit was the normally distributed random variable. Then we normalized , N the mRNA abundance to fi00ðtz1Þ~M|f 0iðtz1Þ P f 0jðtz1Þ j~1 in order to grant P f 00jðtz1Þ~M at this time. For each gene i, if j fi0(t+1),L then the value was set to fi0(t+1) = L. Then we assigned fi(t+1) = fi0(t+1); that is, fi(t+1) is the mRNA abundance of gene i in the next generation. For the time development of gene expression divergence, we assumed that the distribution of mRNA abundance at the initial state (t = 0) followed the Zipf’s law-like distribution generated by the Monte Carlo simulation described above. We executed two independent simulations described above during the same number of generations (T), and we calculated the gene expression divergence defined by Khaitovich et al. between those two Gene Expression Evolution hypothetical species [22]. The computer program for the Monte- Carlo simulation was developed in our laboratory and written in Java version 1.6.0. The fitting procedure was executed using Mathematica version 6.0. Monte Carlo Simulation (2) Random Mating Population with Natural Selection Using a random mating, finite, and constant size (Ne = 100) population of diploid organisms, we assumed that every locus was inherited independently. We assumed that a given number of genes (n) were subject to natural selection, and that the other genes (N-n) were not subject to natural selection. For the intensity of selection (defined as the reciprocal of the proportion of individuals selected from a population to be used as parents) to be the same across different values of n, we used the following selection schemes. We assumed that at the initial state (t = 0), the mRNA abundance followed the Zipf’s law-like distribution generated by the refined neutral model. The initial distribution was identical for all individuals. Here, f1(0).f2(0).….fN(0) were satisfied. The mRNA abundance of each gene in each individual evolved by the process described above. The model parameters were: L = 1.0, M = 300,000, and N = 20,000. In the optimum abundance model, we assumed that the optimum values distribution was given as 10X+1, where X was a random variable having a skewed normal distribution with parameter a = 5 (Figure 5). We generated this random number sequence using an sn (the skew-normal and skew-t distributions) package of R statistics system version 2.8.0. The optimum values satisfied the inequalities h1.h2.….hn, where hi was the optimum value of gene i. The relative fitness of an individual was given by the reciprocal of the sum of Euclidian distances of each gene, from the optimum values to the realized expression levels. Then a given number of individuals having the largest relative fitness survived to propagate the next generation. In the relative abundance model, we assumed the optimum relationships to be f1.f2.….fn. We tested the inequalities for all combinations of genes (C(n,2) = n!/((n-2)!2!) times), and we counted the number of inequalities satisfied in each individual. We used this count as the score of the individual, and a given number of individuals having the largest score survived to propagate the next generation. Since huge numbers of inequalities must be tested for all individuals in every generation, this simulation was very time consuming when we assume that n = N = 20,000. Therefore we confined the number of genes undergoing natural selection to n = 500. The average pair-wise variance was calculated on the n genes undergoing the natural selection. For those simulations, we used Endeavor Pro7000 (Epson) Intel Core i7 3.20 GHz 64 cores and SPARC Enterprise M8000 (Fujitsu) SPARC64 VII 2.52 GHz 632 cores Generalized Chi-Square Goodness-of-Fit Test The consistency of the model and the data were tested by the generalized chi-square goodness-of-fit test. A common problem of the standard Pearson’s chi-square test is that because its power depends on sample size, small and unimportant departures from a specified reference distribution may be detected with large samples in general. The generalized chi-square goodness-of-fit test assumes that the observed distribution may contain a small discrepancy from the reference distribution that is eventually detected as the sample size increases [45]. The formula for the hypothesis (H0) of the generalized chi-square test is H0: d = d0, d0.0, where d is the discrepancy measure defined as d~ X k ðpi{piÞ2, pi is the i-th cell pi probability given by the reference distriibution pi~pið1zeiÞ, where PLoS ONE | www.plosone.org 9 November 2009 | Volume 4 | Issue 11 | e7943