DNA microarray analysis has emerged as a leading technology to enhance our understanding of gene regulation and function in cellular mechanism controls on a genomic scale. This technology has advanced to unravel the genetic machinery of biological rhythms by collecting massive gene-expression data in a time course. Here, we present a statistical model for clustering periodic patterns of gene expression in terms of different transcriptional profiles. The model incorporates biologically meaningful Fourier series approximations of gene periodic expression into a mixture-model-based likelihood function, thus producing results that are likely to be closer to biological relevance, as compared to those from existing models. Also because the structures of the time-dependent means and covariance matrix are modeled, the new approach displays increased statistical power and precision of parameter estimation. The approach was used to reanalyze a real example with 800 periodically expressed transcriptional genes in yeast, leading to the identification of 13 distinct patterns of gene-expression cycles. The model proposed can be useful for characterizing the complex biological effects of gene expression and generate testable hypotheses about the workings of developmental systems in a more precise quantitative way.

---

ALMOST all of life's phenomena populating the earth can be described in terms of periodic rhythms that result from the planet's rotation and orbit around the sun (GOLDBETER 2002). Cell division (MITCHISON 2003), circadian rhythms (CROSTHWAITE 2004; PROLO et al. 2005), morphogenesis of periodic structures such as somites in vertebrates (DALE et al. 2003), and the complex life cycles of some microorganisms (LAKIN-THOMAS and BRODY 2004; ROVERY et al. 2005) are all excellent representatives of biological rhythms. From a mechanistic perspective, rhythmic behavior arises in genetic and metabolic networks as a result of nonlinearities associated with various modes of cellular regulation.

The complexity and inherent periodicity of rhythmic processes can be well described by mathematical models. For example, mathematical models of increasing complexity for the genetic regulatory network producing circadian rhythms in the fly Drosophila predict the occurrence of sustained circadian oscillations of the limit cycle type (LELOUP et al. 1999). When incorporating the effect of light, the models account for phase shifting of the rhythm by light pulses and for entrainment by light–dark cycles. The models also provide an explanation for the long-term suppression of circadian rhythms by a single pulse of light. Stochastic simulations were developed to test the robustness of circadian oscillations with respect to molecular noise (GONZE et al. 2002, 2003).

The recent development of gene-expression technologies has offered biologists a unique opportunity to more closely study the mechanisms that control a particular biological rhythm. PANDA et al. (2002) employed high-density oligonucleotide arrays to trace gene expression in mouse tissue samples taken every 4 hr during two complete circadian cycles. They also identified clusters of circadian-regulated genes among >7000 genes with a cosine wave-fitting algorithm (HARMER et al. 2000). About 650 cycling transcripts were detected to be under circadian regulation specific to either the suprachiasmatic nuclei or the liver. These studies allowed the authors to understand how the major oscillator in the suprachiasmatic nuclei and the liver regulates behavioral and physiological rhythms in the whole organism.

To associate the profile of gene expression with a physiological function of interest, it is crucial to cluster the types of gene expression on the basis of their periodic patterns. Statistical modeling and algorithms play a central role in cataloging dynamic gene-expression profiles. While various computational models have been developed for gene clustering based on static microarray data (EISEN et al. 1998; GHOSH and CHINNAIYAN 2002; MCLACHLAN et al. 2002; RAMONI et al. 2002; FAN and REN 2006), considerable attention has been paid to methodological derivations for detecting temporal patterns of gene expression in a time course based on functional principal component analysis or mixture model analysis (HOLTER et al. 2001; QIAN et al. 2001; BAR-JOSEPH et al. 2003; LUAN and LI 2003; PARK et al. 2003; WAKEFIELD et al. 2003; ERNST et al. 2005; STOREY et al. 2005; MA et al. 2006; NG et al. 2006; INOUE et al. 2007). In particular, LUAN and LI (2004) proposed a statistical model for characterizing different types of periodic rhythms of gene expression. These methods show some unique utilization to identify genes with varying expression profiles over time, providing new tools to elucidate a comprehensive picture of the life process.

A common feature of these clustering methods is that they model time-dependent gene-expression profiles using nonparametric approaches, such as cubic splines. Although nonparametric approaches are statistically flexible and computationally fast, they often do not provide rigorous biological interpretations of results even if the gene-expression data analyzed contain some biologically meaningful patterns. This article was motivated by the fundamental idea that key features of many biological processes can be described by parsimonious mathematical equations. For example, FRANK (1926) used the Fourier series to approximate periodic and quasi-periodic biological phenomena. ATTINGER et al. (1966) provided a detailed mathematical formulation of this approximation for a biological rhythmic system through both theoretical and experimental approaches. Fourier series approximation as an analytical tool (PRIESTLEY 1981) has been widely applied to study the mechanisms and patterns of biological rhythmicity, including the cyclic organization of preterm and term neonates during the neonatal period (BEGUM et al. 2006) and pharmacodynamics (MAGER and ABERNETHY 2007).

Fourier approximation has also been used to model periodic gene expression, leading to the detection of periodic signals in various organisms including yeast and human cells (SPELLMAN et al. 1998; WICHERT et al. 2004; KIM et al. 2006). GLYNN et al. (2006) proposed a Lomb–Scargle periodogram approach based on the fast Fourier transform to model unevenly spaced gene-expression time series and then characterize periodic patterns of gene expression by using a multiple-hypothesis testing procedure with a controlled false discovery rate. Our statistical model presented in this article integrates the Fourier series approximation into a mixture-model framework to mathematically cluster microarray genes on the basis of their distinct patterns of periodic expression. The advantage of this mixture-based model includes its solid statistical foundation of testing the number of gene clusters. Furthermore, through the implementation of the Fourier series approximation, it is possible to test several biologically meaningful characteristics such as sharp peaks in the ordinary periodograms calculated from the Fourier transform of the time series (DURBIN 1967). We use a published example for cell cycle-regulated genes in Saccharomyces cerevisiae (SPELLMAN et al. 1998) to validate the new model. The statistical behavior of the model is examined through simulation studies.

Mixture model:

Finite mixture models (MCLACHLAN and PEEL 2000) have been widely used to model the distributions of a variety of random phenomena. Multivariate normality is generally assumed for multivariate data of a continuous nature. This multivariate normal mixture model is employed to detect different patterns in gene-expression profiles.

Assume that n genes are measured at multiple time points. Our model is able to consider unevenly spaced time intervals and different measurement schedules for each gene. Let be the T_i-dimensional gene-expression data for gene i. Suppose there are J components in the mixture model. This means that any one of the genes (i) is assumed to arise from one (and only one) of the J possible periodic patterns of expression, the distribution of whose expression data is expressed as the J-component mixture probability density function; i.e.,

(1)

where is the mixture proportions that are nonnegative and sum to unity; u_i = (u_i1,..., u_iJ) contains the component- (or pattern-) specific mean vector for gene i; and contains residual variances and covariances among the T_i time points for gene i that are common for all gene-expression patterns. The probability density function of the jth gene-expression pattern or cluster, , is assumed to be multivariate normal with mean vector and common covariance matrix .

From the above analysis, different gene-expression patterns can be detected by estimating the number of mixture components, the time-dependent means of each component, and the covariance matrix. By comparing the differences of the mean vector among different components, it is possible to test how response curves differ among different components and, further, to address fundamentally important biological issues regarding the interplay between gene expression and biological rhythms.

Modeling the mean curve:

If the genes studied are periodically regulated, their time-dependent expression can be accurately approximated by a Fourier series approximation (SPELLMAN et al. 1998). Fourier series approximations can assess periodicity, so we can identify the genes whose RNA levels varied periodically within the cell cycle and, further, find the associated amplitudes and phases of these cycles by estimating the mathematical parameters in the Fourier series approximation. A general form of the Fourier signal is given as

(2)

where a₀ is the average value of s_K(t) and the other a_k and b_k coefficients are the amplitude coefficients that determine the times at which the gene achieves peak and trough expression levels, respectively, and is the period of the signal of gene expression.

In practice, the Fourier series of Equation 2 can be approximated by the first K + 1 term. By analyzing the error of the Fourier series approximation, expressed as the difference between the signal and the (K + 1)-term series,

it is found that, when a third-order approximation is used (i.e., K = 3), the unused terms (k = 4,..., ) from the series together are only 3% of the signal. Thus, the time-dependent expression value of a gene can be adequately modeled by a Fourier series approximation of the first three orders. Let denote the vector containing Fourier parameters of various orders for the genes with pattern j, which is specified as

(3)

The mean expression value of gene i when its pattern is j at time point t_i can be approximated by these Fourier series parameters; i.e.,

Modeling the covariance structure:

A number of approaches can be used to model the covariance structure of serial measurements. A commonly used approach for structuring the covariance is the first-order autoregressive [AR(1)] model (DAVIDIAN and GILTINAN 1995; VERBEKE and MOLENBERGHS 2000). One advantage of using the AR(1) model is that it provides a general expression for calculating the determinant and inverse of the matrix for any number of time points measured. But it assumes variance stationarity and correlation stationarity; i.e., the residual variance at different time points is the same, expressed as ², and the correlation between two different time points for gene i, and , decreases exponentially in with time lag, expressed as corr(.

To remove the heteroscedastic problem of the residual variance, which violates a basic assumption of the simple AR(1) model, two approaches can be used. The first approach is to model the residual variance by a parametric function of time, as proposed by PLETCHER and GEYER (1999). But this approach needs to implement additional parameters for characterizing the time-dependent change of the variance. The second approach is to embed CARROLL and RUPPERT's (1984) transform-both-sides (TBS) model into the growth-incorporated finite mixture model (WU et al. 2004), which does not need any more parameters. Both empirical analyses with real examples and computer simulations suggest that the TBS-based model can increase the precision of parameter estimation and computational efficiency. Furthermore, the TBS model preserves original biological means of the curve parameters although statistical analyses are based on transformed data. In this study, we used the TBS-based AR(1) model, with the structuring parameters arrayed in .

Computational algorithm:

The EM algorithm and Nelder–Mead simplex algorithm were used to estimate the unknown parameters . The observed data ) are regarded as being incomplete. Let z_ij be a missing variable, defined as 1 if y_i arises from the jth component of the mixture model (i = 1,..., n; j = 1,..., J), and write z_i = (z_i1,..., z_iJ)'. The variables z₁,..., z_n are independent and obey a multinomial distribution consisting of one category among J possibilities with probabilities (₁,..., _J). Then, we have

The complete data log-likelihood is

In the E step, z_ij is replaced by the posterior probability (P_ij) of the ith gene that belongs to the jth pattern. This is equal to the conditional expectation (E[z_ij | y_i]) of the complete data log-likelihood log L_c(), given the observed data y, and is computed as

(4)

In the M step, we intend to choose the value of that maximizes . For the derivation process of the estimates of the unknown parameters in the M step, see the APPENDIX. The E and M steps are iterated repeatedly until the parameter estimates converge to finally obtain the maximum-likelihood estimates (MLEs) of the unknown parameters.

Model selection:

Testing for the number of components in a mixture is an important but very difficult problem that has not been completely resolved. One of the leading selection methods is the Akaike information criterion (AIC) (AKAIKE 1974). This is designed to be an approximately unbiased estimator of the expected Kullback–Leibler information of a fitted model. The minimum AIC produces a selected model that is close to the best possible choice. Within the context of the mixture model, the number of components, k, is chosen to minimize

where N(k) is the number of independent parameters within . By the same token, the Bayesian information criterion (BIC) (SCHWARZ 1978) is also available, expressed as

The selected model is the one with the smallest BIC. There is no clear consensus on which criterion is best to use, although the empirical work of FRALEY and RAFTERY (1998) seems to favor the BIC. Since information criteria penalize models with additional parameters, the AIC and BIC model order selection criteria are based on parsimony. Note that since the BIC imposes a greater penalty for additional parameters than does the AIC, the BIC always provides a model with a number of parameters not greater than that chosen by the AIC.

Hypothesis testing:

The significance of overall differences in transcriptional expression profile among different groups of microarray genes is tested by formulating the following hypotheses:

(5)

The log-likelihood ratio (LR) test statistic is then calculated by

(6)

where the tilde and circumflex stand for the MLEs of the unknown parameters under the null and alternative hypotheses, respectively. The critical threshold for claiming distinguishable expression patterns can be determined on the basis of simulation studies. The null hypothesis means that no different patterns of periodic expression exist among the genes studied, whereas the alternative hypothesis states that at least two different patterns can be identified. Under the null hypothesis, time-dependent expression data for n genes are simulated by assuming that the data follow a multivariate normal distribution with mean vector u and covariance matrix . According to Equation 3, individual elements in u are approximated by the Fourier series function of a particular order. A set of parameters that describes the shape of the Fourier series function can be taken from their estimates obtained from real data analyses. Similarly, the structure of is modeled by AR(1) with the two underlying parameters ( and ²) taken from the estimates.

The functional clustering model described in the main text is then used to analyze the expression data simulated under the condition of no distinct groups, provide the estimates of the curve parameters and covariance-structuring parameters under the null and alternative hypotheses, and calculate the LR test statistic. This procedure is repeated 1000 times, leading to 1000 LR values. The empirical distribution of the LR test statistic over 1000 replicates is then examined. The 95th percentile of the empirical distribution is then regarded as the critical threshold for claiming the existence of distinct patterns of periodic gene expression.

A worked example:

SPELLMAN et al. (1998) reported the results of 800 periodically expressed transcriptional genes in the genome of yeast (S. cerevisiae). DNA microarrays were used to analyze mRNA levels for six yeast strains in cell cultures that have been synchronized by three independent methods, -factor arrest, elutriation, and arrest of a cdc15 temperature-sensitive mutant. Each method produces populations of yeast cells synchronized in terms of their phase in the cell cycle.

As described by SPELLMAN et al. (1998), RNA was extracted from each of the samples and from a control sample (unsynchronized cells). cDNAs were labeled with Cy5 fluor (red) for synchronized samples and Cy3 fluor (green) for the controls. Mixed labeled control and experimental cDNAs were hybridized to individual microarrays containing all 6178 yeast genes. The average fluorescence intensity for each fluor within each array spot was recorded. The data of gene expression were measured by normalized fluorescence ratio log₂(Cy5/Cy3) at different time points. All the microarray data reported in SPELLMAN et al. (1998) are given at http://cellcycle-www.stanford.edu.

To validate the usefulness of the model proposed in this article, we analyze time-dependent gene-expression data derived from the cdc15 experiment. The data contain 800 genes, 632 of which have complete data during all 24 time points. The remaining 168 genes contain missing values at some time points. All the 800 genes are analyzed simultaneously.

When the Fourier series approximation is used to cluster the periodic patterns of gene expression, two issues should be determined in the following sequence. First, what is the best order for Fourier series function that explains the time-dependent data? Second, what is the optimal number of components for the mixture model that each correspond to a different expression pattern? Model selection criteria, AIC and BIC, were used to determine the best Fourier series order and best number of mixture components for Spellman et al.'s data. The two criteria provide similar results (Figure 1), although our analysis is mostly based on BIC. It seems that a higher order of Fourier series can better fit time-dependent data than a lower order, reflecting the complexity of dynamic changes of gene expression. A lower order of Fourier series tends to detect a smaller number of gene-expression patterns than a higher order. For example, the first order detects 13 components, whereas the second and third orders detect 18 and >20 components, respectively. This makes sense because a higher order of Fourier series function has more power to discern subtle differences in gene-expression profiles. When closely looking at the BIC curves (Figure 1), the first order displays a dramatic decrease when the number of clusters is 6–8, whereas a dramatic decrease for the second and the third order occurs at 12–14 and 18–20 clusters, respectively.

View larger version (15K): In this window In a new window Download PPT slide   FIGURE 1.— Component number-dependent AIC and BIC values of model fitting by a Fourier series function of order 1–3 for 800 genes collected from the yeast genome.

 
To illustrate different periodic patterns of gene-expression profiles concordant with cell cycles, we used the first-order Fourier series to detect 13 patterns whose profiles (Figure 2) were drawn with the Fourier parameters estimated from the proposed model (Table 1). These 13 patterns differ dramatically in the overall shape of curves as defined by parameter sets . On the basis of these estimates, a number of hypothesis tests can be made about the developmental patterns of gene expression. Table 1 also provides the estimates of the proportions of mixture components. The 13 patterns occur at different frequencies among the observed genes.

View larger version (16K): In this window In a new window Download PPT slide   FIGURE 2.— Thirteen periodic patterns of gene-expression profiles approximated by a first-order Fourier series function for 800 genes collected from the yeast genome.

 

View this table: In this window In a new window   TABLE 1 MLEs of the Fourier parameters for 13 different periodic patterns of gene expression among 800 genes collected from the yeast genome based on the first-order approximation

 

Computer simulation:

We performed Monte Carlo simulation studies to investigate statistical properties of the functional clustering model proposed. A total of 1000 genes were simulated whose time-dependent expression was measured at 24 equally spaced time points. All the genes were sorted into three distinct patterns with varying proportions. The simulated gene-expression profiles follow an arbitrary form of periodic function. As has been mathematically clear, a periodic function can be approximated by a Fourier series function. The Fourier parameters used for our simulation were assigned by values that are within their spaces according to SPELLMAN et al.'s (1998) data. The time-dependent covariance matrix of gene expression was structured by the AR(1) model. Different residual variances were used in the simulation to examine the effect of residual errors on parameter estimation.

The optimal number of components for the simulated data was determined by calculating AIC and BIC values. As shown in Figure 3, the model can correctly estimate the number of components. On the basis of results from 1000 simulation replicates, the model can provide reasonably accurate and precise estimates of all Fourier parameters (Tables 2 and 3). The precision of parameter estimation depends on the proportion of a gene-expression pattern being better for a more frequent than for a less frequent pattern. As expected, increasing residual variance will reduce the estimation precision of parameters (Table 2 vs. Table 3). To show the robustness of the model, an additional simulation study based on a second-order Fourier series approximation was performed. The results suggest that all parameters can be reasonably estimated even if the number of parameters being estimated is increased (Tables 4 and 5).

View larger version (16K): In this window In a new window Download PPT slide   FIGURE 3.— Component number-dependent AIC and BIC values of model fitting by a first- and a second-order Fourier series function for 1000 simulated genes under different residual variances.

 

View this table: In this window In a new window   TABLE 2 MLEs of the first-order Fourier parameters for three periodic patterns of gene expression among 1000 simulated genes with a given set of values (a0, a1, b1, T) under the residual variance of 2 = 0.3

 

View this table: In this window In a new window   TABLE 3 MLEs of the first-order Fourier parameters for three periodic patterns of gene expression among 1000 simulated genes with a given set of values (a0, a1, b1, T) under the residual variance of 2 = 2.0

 

View this table: In this window In a new window   TABLE 4 MLEs of the second-order Fourier parameters for three periodic patterns of gene expression among 1000 simulated genes with a given set of values (a0, a1, b1, T) under the residual variance of 2 = 0.3

 

View this table: In this window In a new window   TABLE 5 MLEs of the second-order Fourier parameters for three periodic patterns of gene expression among 1000 simulated genes with a given set of values (a0, a1, b1, T) under the residual variance of 2 = 2.6

 
NG et al. (2006) proposed a random-effect mixture model for clustering gene-expression profiles through the incorporation of covariate information. When the covariate is time, Ng et al.'s model functions as ours does. However, these two models are different in three aspects. First, our model allows gene expression measured at unevenly spaced time intervals and gene-specific differences in measurement schedule, although these issues can be incorporated into Ng et al.'s model through extensive modifications. Second, we derived a closed form for the estimation of all the Fourier parameters for each gene cluster within the EM framework, whereas Ng et al. estimated the period of the expression cycle for the mean Fourier curve of all genes by using the least-squares estimation approach. Third, our model is able and flexible enough to consider Fourier series approximation of any arbitrary order. We conducted an additional simulation experiment to compare the results from our and Ng et al.'s models. The data from 1000 time-dependent gene expressions were simulated by assuming that these follow a multivariate normal distribution with gene cluster-specific mean vectors each fitted by a group of Fourier parameters and covariance matrix structured by the AR(1) model. Both our and Ng et al.'s models can correctly estimate the number of gene clusters, i.e., three in this simulated example (Table 6). While our model is able to provide reasonably accurate estimates of the Fourier parameters, Ng et al.'s model was quite biased for many parameter estimates. This comparative analysis suggests that our model will perform better than Ng et al.'s model when gene-expression profiles follow Fourier series approximations.

View this table: In this window In a new window   TABLE 6 Comparisons between the results from our and NG et al.'s (2006) models for clustering time series gene-expression profiles that are approximated by the first-order Fourier series

 

The use of microarray gene-expression technologies to understand developmental questions has received considerable attention in recent years (PANDA et al. 2002). Statistical approaches for analyzing time-dependent gene expression data have been proposed (HOLTER et al. 2001; QIAN et al. 2001; BAR-JOSEPH et al. 2003; LUAN and LI 2003; PARK et al. 2003; WAKEFIELD et al. 2003; ERNST et al. 2005; STOREY et al. 2005; MA et al. 2006; NG et al. 2006; INOUE et al. 2007), but most of them are limited in both biological and statistical aspects. First, these approaches mostly based on a clustering analysis were not implemented with biological principles of gene expression that are related to a life process. For this reason, the results obtained from these approaches may not be biologically relevant and, thus, may be less useful for deciphering the developmental machinery of gene expression. The model proposed in this article integrates mathematical aspects of periodic gene expression into the analytical framework, thereby allowing for the interplay between gene expression and development.

Second, many existing approaches to clustering genes of a similar expression pattern on the basis of a similarity measure have not considered the autocorrelation of time series data and, therefore, fail to remove systematical measurement errors. Although some authors implemented time-dependent correlations into their models (e.g., LUAN and LI 2003; NG et al. 2006), biological meanings of gene expression were not well considered. In the model proposed, the statistical principle of functional data analysis has been embedded into the model by structuring the time-dependent covariance matrix. This, on the one hand, de-noises repeated measurement errors and increases the effectiveness of the model and, on the other hand, enhances the model's power due to a reduced number of parameters being estimated. As an illustration, we used a simple AR(1) model to approximate the covariance structure. Other models, such as a structured antedependence model (ZIMMERMAN and NÚÑEZ-ANTÓN 2001), can also be incorporated (see JAFFRÉZIC et al. 2003; ZHAO et al. 2005), allowing the choice of an optimal model for structuring the covariance matrix. Finally, the mixture model-based approach allows for the estimation of the frequencies of various patterns of gene expression and the calculation of the posterior probability of each gene that belongs to a particular pattern.

The mixture-based approach incorporated by Fourier series approximation is a promising technique for detecting periodic gene-expression patterns. A major advantage of this approach lies in its remarkable flexibility to ask and address fundamental biological questions at the interplay between gene-expression and developmental patterns. Several important hypotheses can be made from this approach, including those about the differences of gene expression in Fourier curve shapes, curve features, and duration of gene expression based on individual Fourier parameters that describe biological characteristics of periodic cycles. For example, the peak-to-trough ratio, a_m/b_m, reflects the amplitude of expression profile and can be tested for its differences among the gene groups detected. If the mean curve is modeled with the Fourier series of order one, i.e.,

the hypothesis can be expressed as

The slope of the gene expression profile may change with time, which suggests the occurrence of gene expression x time interaction effects during a time course. The differentiation of u(t) with respect to time t represents a slope of gene expression. If the slopes at a particular time point t* are different between the curves of different gene groups, this means that significant gene expression x time interaction occurs between this time point and next. The test for gene expression x time interaction can be formulated with the hypotheses

The effect of gene expression x time interaction can be examined during a given time course.

The new approach was used to analyze a real data set for periodic gene expression. The results from this approach suggest that it would be useful for the identification of gene clusters in terms of their periodic expression patterns. Through simulation studies, this approach has proved to provide reasonable accuracy and precision of parameter estimation and can be directly used to analyze a real data set of periodic gene expression. This approach can be modified or extended in the following areas. First, the clustering and estimation of different gene-expression profiles depends on the precise estimation of covariance functions. FAN et al. (2007b) proposed a semiparametric approach for modeling the covariance structure, which has been shown to be particularly powerful for functional data collected at irregular and subject-specific time points. The incorporation of Fan et al.'s approach into our functional clustering model is expected to improve its power for gene clustering. Second, when repeated measurement includes a high number of time points, the structuring of the covariance matrix may be quickly problematic. A handful of statistical models for dimension reduction proposed by J. Fan and his group (FAN et al. 2007a; FAN and LV 2008) can be incorporated into our model, in a hope to increase the tractability of high-dimensional data. With these and other modifications, the approach for gene clustering presented in this article could be useful for addressing some development-relevant questions in genetic control of complex biological processes. The computer code for the approach proposed in this article is available at statgen.ufl.edu.

In what follows, we derive the log-likelihood equations for estimating the unknown vector . The log-likelihood of parameters constructed on the basis of the mixture model is expressed as

and the posterior probability with which the ith gene belongs to the jth pattern is defined by Equation 4. Since , we have

By setting it equal to zero, we have

(A1)

By plugging in into the right side of Equation A1, we have

(A2)

Assuming that the model is implemented by a first-order Fourier series approximation, unknown Fourier parameters are specified as , where c_j = (a_0j, a_1j, b_1j). We have

Note that

Let

and then we have

To solve c_j, we conduct the differentiation as

which leads to

We further have

If is invertable, we have

(A3)

Note that can be written as

If we write , where

we have

and

By solving it for ², we have

so

(A4)

For the AR(1) model, we have

and

Let

Then, we have

and

Therefore, we have

By solving the above log-likelihood equation, the MLE of can be obtained as

(A5)

The conditional expectation estimated with Equation 4 in the E step is used to solve the unknown parameters with log-likelihood Equations A2–A5 in the M step. But since it is impossible to derive a closed form for , we implement the simplex algorithm to estimate this parameter in the M step. The E and M steps are repeated until the estimates are stable.

We thank the two anonymous referees for their constructive comments on the manuscript, which led to significant improvement of its presentation. Part of this work was conducted when R. Wu spent his sabbatical leave at Princeton University. The preparation of this manuscript was partially supported by National Science Foundation grant no. 0540745 to R. Wu and the Brain Korea 21 project in 2007.

¹ These authors contributed equally to this work.

AKAIKE, H., 1974 A new look at the statistical model identification. IEEE Trans. Automat. Contr. 19: 716–723.[CrossRef]

ATTINGER, E. O., A. ANNE and D. A. MCDONALD, 1966 Use of Fourier series for the analysis of biological systems. Biophys. J. 6: 291–304.[Medline]

BAR-JOSEPH, Z., G. K. GERBER, D. K. GIFFORD, T. S. JAAKKOLA and I. SIMON, 2003 Continuous representations of time-series gene expression data. J. Comput. Biol. 10: 341–356.[CrossRef][Medline]

BEGUM, E. A., B. MOTOKI, O. MAKOTO, Y. HATSUMI, K. MASATOSHI et al., 2006 Emergence of physiological rhythmicity in term and preterm neonates in a neonatal intensive care unit. J. Circadian Rhythms 4: 11.[CrossRef][Medline]

CARROLL, R. J., and D. RUPPERT, 1984 Power-transformations when fitting theoretical models to data. J. Am. Stat. Assoc. 79: 321–328.[CrossRef]

CROSTHWAITE, S. K., 2004 Circadian clocks and natural antisense. RNA FEBS Lett. 567: 49–54.

DALE, J. K., M. MAROTO, M. L. DEQUEANT, P. MALAPERT, M. MCGREW et al., 2003 Periodic notch inhibition by lunatic fringe underlies the chick segmentation clock. Nature 421: 275–278.[CrossRef][Medline]

DAVIDIAN, M., and D. M. GILTINAN, 1995 Nonlinear Models for Repeated Measurement Data. Chapman & Hall, London.

DURBIN, J., 1967 Tests of serial independence based on the cumulated periodogram. Bull. Int. Stat. Inst. 42: 1039–1049.

EISEN, M. B., P. T. SPELLMAN, P. O. BROWN and D. BOTSTEIN, 1998 Cluster analysis and display of genome-wide expression patterns. Proc. Natl. Acad. Sci. USA 95: 14863–14868.[Abstract/Free Full Text]

ERNST, J., G. J. NAU and Z. BAR-JOSEPH, 2005 Clsutering short time series gene expression data. Bioinformatics 21: i159–i168.[Abstract]

FAN, J., and J. LV, 2008 Sure independence screening for ultra-high dimensional feature space (with discussion). J. R. Stat. Soc. Ser. B 70: 849–911.

FAN, J., and Y. REN, 2006 Statistical analysis of DNA microarray data in cancer research. Clin. Cancer Res. 12: 4469–4473.[Abstract/Free Full Text]

FAN, J., Y. FAN and J. LV, 2007a High dimensional covariance matrix estimation using a factor model. J. Econom (in press).

FAN, J., T. HUANG and R. Z. LI, 2007b Analysis of longitudinal data with semiparametric estimation of covariance function. J. Am. Stat. Assoc. 35: 632–641.

FRALEY, C., and A. E. RAFTERY, 1998 How many clusters? Which clustering method? Answers via model-based cluster analysis. Technical Report 329. Department of Statistics, University of Washington, Seattle.

FRANK, O., 1926 Die theorie der pulswellen. Zool. Biol. 85: 91–130.

GHOSH, D., and A. M. CHINNAIYAN, 2002 Mixture modeling of gene expression data from microarray experiments. Bioinformatics 18: 275–286.[Abstract/Free Full Text]

GLYNN, E. F., J. CHEN and A. R. MUSHEGIAN, 2006 Detecting periodic patterns in unevenly spaced gene expression time series using Lomb-Scargle periodograms. Bioinformatics 22: 310–316.[Abstract/Free Full Text]

GOLDBETER, A., 2002 Computational approaches to cellular rhythms. Nature 420: 238–245.[CrossRef][Medline]

GONZE, D., J. HALLOY and A. GOLDBETER, 2002 Stochastic versus deterministic models for circadian rhythms. J. Biol. Phys. 28: 637–653.[CrossRef]

GONZE, D., J. HALLOY, J. C. LELOUP and A. GOLDBETER, 2003 Stochastic models for circadian rhythms: influence of molecular noise on periodic and chaotic behavior. C. R. Biol. 326: 189–203.[CrossRef][Medline]

HARMER, S. L., J. B. HOGENESCH, M. STRAUME, H.-S. CHANG, B. HAN et al., 2000 Orchestrated transcription of key pathways in Arabidopsis by the circadian clock. Science 290: 2110–2113.[Abstract/Free Full Text]

HOLTER, N. S., A. MARITAN and M. CIEPLAK, 2001 Dynamic modeling of gene expression data. Proc. Natl. Acad. Sci. USA 98: 1693–1698.[Abstract/Free Full Text]

INOUE, L. Y., M. NEIRA, C. NELSON, M. GLEAVE and R. ETZIONI, 2007 Cluster-based network model for time-course gene expression data. Biostatistics 8: 507–525.[Abstract/Free Full Text]

JAFFRÉZIC, F., R. THOMPSON and W. G. HILL, 2003 Structured antedependence models for genetic analysis of multivariate repeated measures in quantitative traits. Genet. Res. 82: 55–65.[CrossRef][Medline]

KIM, B.-R., R. C. LITTELL and R. L. WU, 2006 Clustering the periodic pattern of gene expression using Fourier series approximations. Curr. Genomics 7: 197–203.[CrossRef]

LAKIN-THOMAS, P. L., and S. BRODY, 2004 Circadian rhythms in microorganisms: new complexities. Annu. Rev. Microbiol. 58: 489–519.[CrossRef][Medline]

LELOUP, J.-C., D. GONZE and A. GOLDBETER, 1999 Limit cycle models for circadian rhythms based on transcriptional regulation in Drosophila and Neurospora. J. Biol. Rhythms 14: 433–448.[Abstract/Free Full Text]

LUAN, Y., and H. LI, 2003 Clustering of time-course gene expression data using a mixed-effects model with B-splines. Bioinformatics 19: 474–482.[Abstract/Free Full Text]

LUAN, Y., and H. LI, 2004 Model-based methods for identifying periodically expressed genes based on time course microarray gene expression data. Bioinformatics 20: 332–339.[Abstract/Free Full Text]

MA, P., C. I. CASTILLO-DAVIS, W. ZHONG and J. S. LIU, 2006 A data-driven clustering method for time course gene expression data. Nucleic Acids Res. 34: 1261–1269.[Abstract/Free Full Text]

MAGER, D. E., and D. R. ABERNETHY, 2007 Use of wavelet and fast Fourier transforms in pharmacodynamics. J. Pharmacol. Exp. Ther. 321: 423–430.

MCLACHLAN, G., and D. PEEL, 2000 Finite Mixture Models. John Wiley & Sons, New York.

MCLACHLAN, G. J., R. W. BEAN and D. PEEL, 2002 A mixture model-based approach to the clustering of microarray expression data. Bioinformatics 18: 414–422.

MITCHISON, J. M., 2003 Growth during the cell cycle. Int. Rev. Cytol. 226: 165–258.[Medline]

NG, S. K., G. J. MCLACHLAN, K. WANG, L. B.-T. JONES and S.-W. NG, 2006 A mixture model with random-effects components for clustering correlated gene-expression profiles. Bioinformatics 22: 1745–1752.[Abstract/Free Full Text]

PANDA, S., T. K. SATO, A. M. CASTRUCCI, M. D. ROLLAG, W. J. DEGRIP et al., 2002 Melanopsin (Opn4) requirement for normal light-induced circadian phase shifting. Science 298: 2213–2216.[Abstract/Free Full Text]

PARK, T., S. G. YI, S. LEE, S. Y. LEE, D. H. YOO et al., 2003 Statistical tests for identifying differentially expressed genes in time-course microarray experiments. Bioinformatics 19: 694–703.[Abstract/Free Full Text]

PLETCHER, S. D., and C. J. GEYER, 1999 The genetic analysis of age-dependent traits: modeling a character process. Genetics 153: 825–835.[Abstract/Free Full Text]

PRIESTLEY, M. B., 1981 Spectral Analysis and Time Series. Academic Press, San Diego.

PROLO, L. M., J. S. TAKAHASHI and E. D. HERZOG, 2005 Circadian rhythm generation and entrainment in astrocytes. J. Neurosci. 25: 404–408.[Abstract/Free Full Text]

QIAN, J., B. STENGER, C. A. WILSON, J. LIN, R. JANSEN et al., 2001 Partslist: a web-based system for dynamically ranking protein folds based on disparate attributes, including whole-genome expression and interaction information. Nucleic Acids Res. 29: 1750–1764.[Abstract/Free Full Text]

RAMONI, M. F., P. SEBASTIANI and I. S. KOHANE, 2002 Cluster analysis of gene expression dynamics. Proc. Natl. Acad. Sci. USA 99: 9121–9126.[Abstract/Free Full Text]

ROVERY, C., M. V. LA, S. ROBINEAU, K. MATSUMOTO, P. RENESTO et al., 2005 Preliminary transcriptional analsysis of spoT gene family and of membrane proteins in Rickettsia conorii and Rickettsia felis. Ann. NY Acad. Sci. 1063: 79–82.[CrossRef][Medline]

SCHWARZ, G., 1978 Estimating the dimension of a model. Ann. Stat. 6: 461–464.[CrossRef]

SPELLMAN, P. T., S. SHERLOCK, M. Q. ZHANG, V. R. IYER, K. ANERS et al., 1998 Comprehensive identification of cell cycle-regulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization. Mol. Biol. Cell 9: 3273–3297.[Abstract/Free Full Text]

STOREY, J. D., W. XIAO, J. T. LEEK, R. G. TOMPKINS and R. W. DAVIS, 2005 Significant analysis of time course microarray experiments. Proc. Natl. Acad. Sci. USA 102: 12837–12842.[Abstract/Free Full Text]

VERBEKE, G., and G. MOLENBERGHS, 2000 Linear Mixed Models for Longitudinal Data. Springer-Verlag, New York.

WAKEFIELD, J., C. ZHOU and S. G. SELF, 2003 Modelling gene expression data over time: curve clustering with informative prior distributions, pp. 721–732 in Bayesian Statistics 7, Proceedings of the Seventh Valencia International Meeting, edited by J. M. BERNARDO, M. J. BAYARRI, J. O. BERGER, A. P. DAWID, D. HECKERMAN, A. F. M. SMITH and M. WEST. Oxford University Press, London/New York/Oxford.

WICHERT, S., K. FOKIANOS and K. STRIMMER, 2004 Identifying periodically expressed transcripts in microarray time series data. Bioinformatics 20: 5–20.[Abstract/Free Full Text]

WU, R. L., C.-X. MA, M. LIN, Z. H. WANG and G. CASELLA, 2004 Functional mapping of quantitative trait loci underlying growth trajectories using a transform-both-sides logistic model. Biometrics 60: 729–738.[CrossRef][Medline]

ZHAO, W., Y. Q. CHEN, G. CASELLA, J. M. CHEVERUD and R. L. WU, 2005 A nonstationary model for functional mapping of complex traits. Bioinformatics 21: 2469–2477.[Abstract/Free Full Text]

ZIMMERMAN, D. L., and V. NÚÑEZ-ANTÓN, 2001 Parametric modeling of growth curve data: an overview (with discussion). Test 10: 1–73.[CrossRef]

Communicating editor: C. HALEY

International Access Link

Online ISSN: 1943-2361  Print ISSN: 0016-6731
Copyright 2009 by the Genetics Society of America
phone: 412-268-1812   fax: 412-268-1813   email: genetics-gsa{at}andrew.cmu.edu