BACKGROUND
a major goal of systems biology is to integrate biological functions of individual genes in terms of their interactions.
time course gene expression profiling, which can capture the global transcriptional responses to signals during a biological process of interest, offers a major data source to achieve this goal  <cit> .
in network modeling of gene expression data, assessing pair-wise relationships is often a starting point.
in early days, correlation coefficient  <cit> , euclidean distance, as well as their variations, such as partial correlations, empirical bayes and bootstrap methods  <cit> , were used.
they are effective for computing direction free linear dependence when the data are independent.
networks constructed this way are essentially co-expression networks.
while having the appeal of being simple and intuitive, correlation metrics have limitations when applied to time course data.
they assume independence of the order of the data points, while in reality the data at each time step depend on the previous time points.
ignoring the inter-time point dependence not only loses sensitivity toward detecting interactions but could also lead to erroneous predictions  <cit> .
significant phase shift in the timing of expression changes have also been observed for highly associated genes  <cit> .
some studies tried to identify the phase lag directly by shifting gene expression time series with respect to each other until the optimal alignment is reached.
for instance, qian et al proposed a local clustering approach based on optimal pair-wise alignment through dynamic programming  <cit> ; schmitt et al used the pearson's correlation  <cit> , balasubramaniyan et al used the spearman rank correlation  <cit> , pereda et al used cross correlation  <cit> , to compute the maximum time-lagged similarity between two transcript profiles, and utilized the results to identify clusters.
the degree of lag varies widely in different gene pairs, and these approaches need multiple runs to find the lag that best aligns each pair.
the performance of the alignment depends on whether the lags are close to integer numbers of the sampling steps of the experiment.
more sophisticated methods were also developed.
aach and church implemented both simple and interpolative time warping based on dynamic programming to identify an optimal alignment of two gene expression time series  <cit> .
expanding this approach, liu and müller proposed a non-parametric time-synchronized iterative mean updating technique to construct modes of temporal structure in gene expression profiles  <cit> .
bar-joseph et al.
<cit>  developed an approach to align temporal data sets using piecewise spline fitting, extracting shift and stretch parameters for each data set.
butte et al.
<cit>  utilized digital signal-processing tools, including power spectral densities, coherence, transfer gain, and phase shift, to find pair-wise gene associations based on periodically expressed time-invariant gene profiles.
more recently, a hidden markov model based approach was utilized to infer the timing in gene expression changes under different experimental conditions  <cit> .
linear and non-linear multivariate analysis and signal processing techniques were also introduced to analyze time series microarray data  <cit> .
several studies used pair wise mutual information to infer interactions and regulatory relationships between genes  <cit> .
this method assumes a fixed time delay, which might not be true across different experimental conditions.
in frequency domain time series analysis, causality and interrelationship among the components can be studied using coherence and partial coherence.
graphical models based on such analysis have been studied by butte et al  <cit>  and salvador et al  <cit> .
however, albo et al  <cit>  showed that partial coherence-based causality measures are sensitive to measurement noise.
apparently, more studies are needed to fully utilize the dynamics underlying the temporal gene expression pattern, and to better understand the complex spatial-temporal architecture of transcriptome.
recently, increasing evidence, including those from the advancement of single-cell time course gene expression profiling technologies  <cit> , suggest that like other complex dynamic systems in nature, synchrony through oscillation and frequency modulation is a general strategy for an organism to coordinate the transcription of multiple target genes in responses to external signals  <cit> .
examples include the p53-mdm <dig> feedback loop  <cit> , the nf-κb signaling pathway  <cit> , and calcium responsive pathways  <cit> .
these further emphasize the need of new methods to study and utilize the dynamics.
the oscillations in gene expression, like other oscillations in biological systems  <cit> , are most often pulsatile or relaxed oscillations rather than harmonic, thus calling for mathematical methods rooted from phase space analysis  <cit> .
in this study, we investigate the potential of network inference using the phase locking analysis technique  <cit> .
this approach is based on the following concepts originating from nonlinear dynamics  <cit> : if two time series interact with each other, there will be a process of rhythmic adjustment resulting from the interaction, leading to phase locking.
phase-locked oscillators progress through their trajectories in phase space at the same pace , or rational ratios with respect to each other .
conversely, such phase locking phenomenon can be utilized to infer interaction between two dynamic systems, even for weak interactions  <cit> .
recently kim et al clustered genes of synchronized oscillatory pattern  during yeast cell cycle, and observed that genes in the same cluster were closely associated, as evidenced by the sharing of go terms and biogrid interactions  <cit> .
in this study we will apply the phase locking analysis to the stanford yeast cell cycle data  <cit> , and examine the phase locking  between transcription factors and targets, between gene pairs with prior evidence of other types of interaction, and between cell cycle genes.
based on the results, we will propose a new network inference approach utilizing the phase locking index, and examine the modular structure of the networks constructed and the biological themes shared by genes in the network modules.
RESULTS
phase locking of interacting genes
distribution of λ
we first examined the distribution of phase locking index λ in all four datasets of he stanford yeast cell cycle study  <cit> , the histograms are given in figure  <dig>  in each case the distribution is close to normal , indicating that the way we define the threshold λ for significant locking,  <dig> standard deviations above mean, is reasonable.
in all histograms there is a hint of a high-λ tail, likely contributed from gene pairs that are phase locked.
phase locking between the cell cycle regulating tfs and their targets
following the original gene expression study of cell cycles  <cit> , several groups have investigated yeast transcription binding using the chip-chip technology  <cit> .
these data provide useful information of which genes are potentially transcription regulation targets of each tf.
we have obtained the data from simon et al, where the promoter binding by the  <dig> known cell cycle regulating tfs were studied  <cit> .
both chip-chip and microarray data are noisy, and we found no direct quantitative dependence of λ or r on the binding p-value .
however, we obtained interesting results when using the target and non-target control groups as benchmarks.
the receiver operating characteristic  curve was used to determine whether or not λ can potentially discriminate targets from non-targets.
roc is a graphical plot of the sensitivity versus 1-specificity, namely the fraction of true positives versus the fraction of false positives, as the discrimination threshold of a classifier is varied.
the area under curve  reflects the performance.
the roc of a random classifier would be a 45° line with auc =  <dig> .
in figure  <dig> the roc plot for each tf in alpha factor arrest experiment is given.
overall, λ is significantly better than either a random classifier or r , suggesting that expression levels from tf and target genes are likely to exhibit phase locking.
in contrast, the expression levels of tf and their targets seem to correlate little.
in fact, using r does not fare any better than a random classifier .
the results from other datasets are similar, and the auc of roc plots from all four datasets are summarized in figure  <dig>  the list of target genes that exhibit significant phase locking with the tf is given in additional file  <dig>
phase locking of biogrid gene pairs
to further investigate the potential of phase locking analysis in network inference, we examined phase locking between gene pairs that have evidence of other types of interaction according to biogrid http://www.thebiogrid.org.
biogrid is a freely accessible database of physical  and genetic interactions, curated from high-throughput data and literature  <cit> .
of all possible gene pairs in our data sets, ~ <dig>  are annotated in biogrid.
we constructed five bootstrapping  <cit>  sets that consisted of the same number of biogrid gene pairs, randomly sampled from all possible pairs.
the distribution of phase locking index was examined in each group.
we found that the distribution for the biogrid pairs is skewed toward higher λ than the bootstrapping sets .
the odds ratio  of the enrichment of phase locked pairs at a given threshold value λ  <dig>  when there is prior evidence of interaction according to biogrid, can be calculated by:
odds ratio=pp
where p is the probability of a biogrid gene pair having a significant phase locking λ >λ  <dig> , and p is such probability of a bootstrapping gene pair.
the or of the four datasets at different λ  <dig> is given in figure  <dig>  evidently, it increases with more stringent λ  <dig>  and is significantly greater than  <dig> above the cutoff values of λ that we used in this study .
these results further validate the idea that λ can be utilized to identify interacting gene pairs.
phase locking among cell cycle genes
we anticipate the expression of the  <dig> cell cycle genes to show significant coordination, as they should all respond to the master signals that drive the cell cycle, and form an interaction network that together regulate the cell cycle progression.
we found that indeed they exhibit phase locking among themselves more often than with random genes or than within random genes.
on the other hand, when evaluated by the correlation coefficient r of their expression levels, cell cycle genes are not different from random genes.
to illustrate the point we present in figure  <dig> colormaps of the adjacency matrix a from the alpha factor arrest dataset, where the non-zero elements are presented by a black pixel.
according to λ, in the subregion of cell cycle genes, there are significantly more black pixels than in the subregion of random genes, or in the subregion of cross interactions.
to give a more quantitative description in figure  <dig> we give the network degrees of all genes in the network.
evidently, based on phase locking, the cell cycle genes form a highly connected subnetwork, whilst random genes are sparsely connected.
in sharp contrast, according to r, there is literally no difference between the cell cycle genes and the random genes, both with moderate connectivity.
results from the other three data sets are similar.
note that in each case the total number of genes is less than their purported values , this is because we have removed genes with more than one missing values in their time series measurements.
phase lags and m:n phase locking
in general we found that the phase lag Ψ <dig> between locked pairs varies widely, and most of them are not close to  <dig>  this is consistent with observations made by others  <cit> .
figure  <dig> gives the distribution of the relative phase Ψ <dig> as defined in equation  <dig> of phase locked tf and target pairs.
they span the whole range of , with no consistent pattern.
this is not surprising as it is the protein of the tf that interact directly with the target gene transcription, not the expression of the tf transcript  <cit> .
indeed, it has been observed that response delay varies widely in gene expression regulation  <cit> .
additional file  <dig> presents some examples of phase locked biogrid gene pairs with significant phase lags, where the correlations are low, even with time-lagged correlations.
more examples and movies of phase lagged gene pairs are available at our website http://zen.dom.uab.edu:8080/phase/demo/.
the universal existence of phase lags and their variation suggest the advantage of a method like the phase locking analysis that can automatically account for them.
the pair-wise alignment approaches will need to align each pair individually thus adding an extra step in network modeling.
in this study we included the consideration of m:n phase locking, up to max{m, n} =  <dig>  we found that there were a significant number of gene pairs exhibiting higher than 1: <dig> locking.
the number breakdowns are listed in additional file  <dig>  in the cell cycle gene analysis, on average there are  <dig>  gene pairs in each data set that exhibit significant phase locking.
out of these pairs, ~ <dig> % are from higher order locking, with ~half contributed from the 1: <dig> and 2: <dig> locking.
in the tf-target pairs, the proportion of higher order locking is significantly higher at ~10%.
on average, in each dataset there are  <dig> significantly phase locked tf-target pairs.
<dig> are due to higher than 1: <dig> locking, and ~70% of them are contributed from the 1: <dig> or 2: <dig> locking.
the number of phase locked pairs, versus max{m, n}, is given in figure  <dig>  interestingly, an exponential dependence is evident.
including higher order phase locking significantly improved the performance of λ as a classifier to detect tf targets.
in figure  <dig> we also included results that only considered 1: <dig> locking.
clearly, the performance is not as good as when higher order locking is included , though it is still better than r or a random classifier .
the basis of higher order phase locking is the different intrinsic dynamic time scales of the different time series.
figure  <dig> presents some examples of order  <dig> locked biogrid gene pairs.
the corresponding pubmed id that contains evidence of their interaction is also given.
more examples of high order locked pairs from the cell cycle genes, and tf-target pairs are given in additional files  <dig> & <dig>  in these examples, the values of λ <dig> , and r are low, and interaction would have been missed if using them as the metric.
demo movies of higher order phase locking are available at our website http://zen.dom.uab.edu:8080/phase/demo/.
the wide range of phase lags and the significant proportion of higher order locking emphasize the need of a rigorous method like the phase locking analysis that can automatically take care of them.
agreement between the  <dig> datasets
as observed by others, we see large variation among the  <dig> datasets.
however, when we examine the phase locked pairs, we found that there is significant concordance between the four datasets, with p <  <dig>  .
the venn diagram is given in figure 10a.
interestingly, we find that often for gene pairs that show consistent phase locking in multiple datasets, their correlation do not follow the trend.
figure 10b is an example, where the tf-target pair swi5-ash <dig>  <cit> , shows high λ consistently in the alpha, cdc <dig> and elu experiments, but their expressions correlate in no dataset.
network inference using phase locking
in our study so far, we have demonstrated that phase locking in expression changes is a good indicator of interaction.
it is therefore natural to utilize it to construct gene interaction networks.
highly connected genes tend to be essential genes
we first examined the network defined by the adjacency matrix given in equation  <dig>  analysis of the yeast protein-protein interaction network in the past has revealed that highly connected genes are more likely to be essential for survival  <cit> .
in a recent study of co-expression network by zhang et al  <cit> , strong positive correlations between the network degree and the functional essentiality of genes were also observed.
in our phase locked networks, we also observed a positive correlation between the network degree ki of a gene and its likelihood of being essential.
the relationship is nicely captured in figure 11a.
all genes were ranked according to their ki and partitioned into  <dig> equal sized bins.
the proportion of essential genes in each bin was then plotted against the mean ki.
a linear dependence is evident .
we also examined the degree distribution of essential versus non-essential genes, given in figure 11b.
compared to the non-essential ones, essential genes are significantly skewed toward having higher numbers of phase locked partners .
genes in network modules have focused biological themes
the structure of the phase locked network is further studied through hierarchical clustering using the topological overlapping matrix given in equation  <dig>  as depicted in figure  <dig> modular organization is evident.
genes in the grey regions did not belong to any modules .
most of the modules  contain higher numbers of essential genes than the scattered set .
based on the theory of phase locking and our results so far, one anticipates genes from the same module to be highly interactive among themselves.
this is indeed true for all the modules when we annotate them with biogrid information .
in addition, ontological analysis of each identified module was performed using gostat  <cit> , to examine the biological themes and the functional relationship of the module members.
we found that genes in the same module tend to be involved in a set of focused signaling pathways.
at a stringent cutoff p =  <dig> , genes in the  <dig> modules share a total of  <dig> go categories, and  <dig> at p =  <dig> , whilst random gene sets of the same sizes and the scattered gene set do not share any go category even at a very loose cut off p =  <dig> .
in table  <dig> we list the go biological processes shared by genes in each module at p <  <dig> .
a common theme is evident that most of these processes are critical to dna replication and cell cycle regulation.
for instance, genes in the turquoise module likely form a cell cycle regulatory module, while the yellow module is likely associated with rna processing.
discussion
the importance of phase space information
synchrony through oscillation is a common, maybe the most efficient, way to coordinate regulation in complex nonlinear dynamic systems.
it is also a ubiquitous phenomenon in biological systems, where oscillations are observed in all organisms across a wide range of temporal and spatial scales, and are believed to play an important role in maintaining homeostasis and delivering encoded information  <cit> .
examples include the synchronized oscillations in interneuron networks, pulsatile endocrine hormone secretion, circadian oscillators, somite segmentation, and innumerable others.
higher order phase locking occurs frequently, reflecting the multi-stability of complex systems, and is believed important to function.
for instance, in the study of cardiorespiratory synchronization, when plotting the instantaneous respiratory phase at the occurrence of a heartbeat versus time, schäfer and co-workers found n: <dig> synchronization between heart and respiration  <cit> .
in a study with anesthetized rats, stefanovska et al.
further observed lengthy synchronization epochs, and transitions from one ratio to another.
they suggested that such transitions might be useful in monitoring depth of anesthesia  <cit> .
increasing evidence suggests that, as in many other complex systems in nature, oscillation and frequency modulation is also a general strategy for an organism to coordinate multi-gene responses to external signals  <cit> .
it is found that the transcription factor activities, rather than levels of transcription factor expression, mediate transcriptional regulations  <cit> .
in the negative feedback loop between the tumor suppressor p <dig> and the oncogene mdm <dig>  p <dig> is expressed in a series of discrete pulses after dna damage, leading to oscillations in mdm <dig>  <cit> .
the amplitude of the oscillations was much more variable than the period, suggesting strong temporal regulation.
in the nf-κb signaling pathway, nf-κb  localization showed asynchronous oscillations following cell stimulation that decreased in frequency with increased iκbα transcription.
transcription of target genes depended on its oscillation persistence, and thus the functional consequences of nf-κb signaling likely depend on temporal characteristics of the oscillations  <cit> .
in yeast cells it has been shown for several calcium stress responsive tfs  that calcium concentration controls the frequency, but not the duration, of their oscillatory localization bursts, and the oscillation propagates to the expression of downstream genes.
it has been argued that such frequency modulation of localization bursts ensures proportional expression of multiple target genes across a wide dynamic range of expression levels  <cit> .
these facts all imply the importance to study synchrony in expression oscillation, to understand the information encoded and the underlying interaction/regulation mechanisms.
data from these studies  <cit>  also indicate that oscillations in gene expression, like most other oscillations in biological systems, are often pulsatile rather than harmonic.
therefore, mathematical methods rooted from phase space analysis are desirable.
the latter can potentially lead to new efficient network modeling algorithms, and help to understand the complex spatial-temporal architecture of transcriptome.
advantages and limitations of the phase locking analysis
additional to the theoretical appeals, we believe that phase locking analysis has several advantages in clustering genes of similar patterns and in network modeling.
firstly, compared to approaches that primarily rely on the similarities in the amplitude domain patterns, phase locking utilizes the dynamics underlying the temporal pattern, which is more robust against noise.
this is particularly appealing in network modeling of gene expression data, as they usually contain high noise.
also the transcript abundance measurements often reflect a compressed, even altered representation of the true expression changes due to technical issues  <cit> .
these can significantly mask the true patterns in amplitude changes.
in contrast, phase locking analysis, which focuses on the timing of the changes, will be less affected by the noise and the technical issues in the microarray gene expression study.
in fact, it is known that noisy coupled nonlinear dynamic systems may synchronize in phase whilst their amplitudes remain uncorrelated  <cit> .
in our analysis, we have seen ample examples where interacting gene pairs  exhibit obvious phase locking, but have very low correlation .
secondly, in phase locking analysis, the phase lags in gene expression changes between different genes are automatically accounted for, and the performance is not affected by the amount of lag.
on the other hand, the performance of alignment approaches depends on whether the lag is close to an integer number of the time steps of the experiment, and they need to be adjusted for each pair as the lags of different gene pairs vary greatly.
thirdly, phase locking does not require the two time series to have the same dynamic time scale, or the same frequency.
it is known that some pathways or gene groups in a cell respond to external signals at a much faster time scale than others  <cit> .
high order m:n phase locking analysis can take care of such interacting gene pairs, whilst they would be missed by the alignment method.
a limitation of the phase locking analysis is its reliance on the temporal spectrum to accurately derive the instantaneous phase, which could significantly affect results when either the number of sampling points or the sampling frequency is too low.
note that higher sampling frequencies are needed to detect high order phase locking.
we have carried out a set of numerical simulations and observed significant deterioration in performance when the number of sampling points is reduced to ~ <dig> or lower .
this limitation is not unique to the phase locking analysis.
all pair-wise alignment approaches suffer from the same limitations as they all rely on adequate sample size to make a good assessment of whether the expression patterns of the pair are similar or not.
there is also a caveat with the application of phase locking analysis to network modeling.
gene interactions or phase locking occurs inside each cell.
high throughput time series gene expression studies commonly measure a population of cells all at one time, effectively averaging the expressions of each gene across the whole cell population.
in processes where there is high synchrony in the whole cell population, such as the cell cycle study presented here, phase locking between gene expression changes that occur inside each cell is preserved at population level, and can be detected from population measurements.
in other biological processes, where cellular heterogeneity plays a key role, information of the signaling dynamics and phase locking inside each cell could be lost in population-level measurement.
again, in such situations, the performance of the co-expression alignment approaches to detect interaction will also be affected.
lastly, the recent advancements in single-cell techniques has enabled the generation of time-series gene expression measurements in a large number of individual living cells  <cit> .
we believe that phase locking analysis will be particularly suitable for such data.
the dynamic information at the level of individual living cells will be critical to unravel how a genetic network operates at the systems level.
CONCLUSIONS
a major challenge in systems biology is to reconstruct gene networks that are involved in basic cellular processes, and to understand how alterations to the network affect functions and consequently phenotypes.
interactions between genes can result in expression amplitude variations as well as temporal patterns.
therefore, network inference utilizing temporal domain information deserves more attention.
in this study, we investigated the potential of the phase locking analysis in network modeling of time course gene expression data.
we demonstrated that interacting gene pairs, including transcription regulation interaction, protein-protein, or genetic interaction, are more likely to exhibit phase locking in their gene expression changes, and vice versa.
among the phase locked pairs, up to ~10% are contributed from higher order locking, and the relative phase difference spans across the whole range of .
based on these findings, we constructed interaction networks and revealed that genes with higher network degrees are more likely to be essential genes.
utilizing the phase-locking index based topological overlapping matrix, we further investigated the modular structures in the network.
we showed that genes forming network modules are more likely to be essential genes than scattered genes in the network, and members of the same module tend to be involved in the same biological functions and processes.
in view of the importance of the frequency domain signal in transcription regulation, we believe that the phase locking analysis can potentially lead to new network modeling approaches and help to understand the dynamic designs of the intracellular signaling networks.

 a major goal of systems biology is to integrate biological functions of individual genes in terms of their interactions. time course gene expression profiling, which can capture the global transcriptional responses to signals during a biological process of interest, offers a major data source to achieve this goal. in network modeling of gene expression data, assessing pair-wise relationships is often a starting point. in early days, correlation coefficient, euclidean distance, as well as their variations, such as partial correlations, empirical bayes and bootstrap methods, were used. they are effective for computing direction free linear dependence when the data are independent. networks constructed this way are essentially co-expression networks. while having the appeal of being simple and intuitive, correlation metrics have limitations when applied to time course data. they assume independence of the order of the data points, while in reality the data at each time step depend on the previous time points. ignoring the inter-time point dependence not only loses sensitivity toward detecting interactions but could also lead to erroneous predictions. significant phase shift in the timing of expression changes have also been observed for highly associated genes. some studies tried to identify the phase lag directly by shifting gene expression time series with respect to each other until the optimal alignment is reached. for instance, qian et al proposed a local clustering approach based on optimal pair-wise alignment through dynamic programming; schmitt et al used the pearson's correlation, balasubramaniyan et al used the spearman rank correlation, pereda et al used cross correlation, to compute the maximum time-lagged similarity between two transcript profiles, and utilized the results to identify clusters. the degree of lag varies widely in different gene pairs, and these approaches need multiple runs to find the lag that best aligns each pair. the performance of the alignment depends on whether the lags are close to integer numbers of the sampling steps of the experiment. more sophisticated methods were also developed. aach and church implemented both simple and interpolative time warping based on dynamic programming to identify an optimal alignment of two gene expression time series. expanding this approach, liu and müller proposed a non-parametric time-synchronized iterative mean updating technique to construct modes of temporal structure in gene expression profiles. bar-joseph et al. developed an approach to align temporal data sets using piecewise spline fitting, extracting shift and stretch parameters for each data set. butte et al. utilized digital signal-processing tools, including power spectral densities, coherence, transfer gain, and phase shift, to find pair-wise gene associations based on periodically expressed time-invariant gene profiles. more recently, a hidden markov model based approach was utilized to infer the timing in gene expression changes under different experimental conditions. linear and non-linear multivariate analysis and signal processing techniques were also introduced to analyze time series microarray data. several studies used pair wise mutual information to infer interactions and regulatory relationships between genes. this method assumes a fixed time delay, which might not be true across different experimental conditions. in frequency domain time series analysis, causality and interrelationship among the components can be studied using coherence and partial coherence. graphical models based on such analysis have been studied by butte et al and salvador et al. however, albo et al showed that partial coherence-based causality measures are sensitive to measurement noise. apparently, more studies are needed to fully utilize the dynamics underlying the temporal gene expression pattern, and to better understand the complex spatial-temporal architecture of transcriptome. recently, increasing evidence, including those from the advancement of single-cell time course gene expression profiling technologies, suggest that like other complex dynamic systems in nature, synchrony through oscillation and frequency modulation is a general strategy for an organism to coordinate the transcription of multiple target genes in responses to external signals. examples include the p53-mdm feedback loop, the nf-κb signaling pathway, and calcium responsive pathways. these further emphasize the need of new methods to study and utilize the dynamics. the oscillations in gene expression, like other oscillations in biological systems, are most often pulsatile or relaxed oscillations rather than harmonic, thus calling for mathematical methods rooted from phase space analysis. in this study, we investigate the potential of network inference using the phase locking analysis technique. this approach is based on the following concepts originating from nonlinear dynamics: if two time series interact with each other, there will be a process of rhythmic adjustment resulting from the interaction, leading to phase locking. phase-locked oscillators progress through their trajectories in phase space at the same pace, or rational ratios with respect to each other. conversely, such phase locking phenomenon can be utilized to infer interaction between two dynamic systems, even for weak interactions. recently kim et al clustered genes of synchronized oscillatory pattern during yeast cell cycle, and observed that genes in the same cluster were closely associated, as evidenced by the sharing of go terms and biogrid interactions. in this study we will apply the phase locking analysis to the stanford yeast cell cycle data, and examine the phase locking between transcription factors and targets, between gene pairs with prior evidence of other types of interaction, and between cell cycle genes. based on the results, we will propose a new network inference approach utilizing the phase locking index, and examine the modular structure of the networks constructed and the biological themes shared by genes in the network modules.  phase locking of interacting genes distribution of λ we first examined the distribution of phase locking index λ in all four datasets of he stanford yeast cell cycle study, the histograms are given in figure in each case the distribution is close to normal, indicating that the way we define the threshold λ for significant locking, standard deviations above mean, is reasonable. in all histograms there is a hint of a high-λ tail, likely contributed from gene pairs that are phase locked. phase locking between the cell cycle regulating tfs and their targets following the original gene expression study of cell cycles, several groups have investigated yeast transcription binding using the chip-chip technology. these data provide useful information of which genes are potentially transcription regulation targets of each tf. we have obtained the data from simon et al, where the promoter binding by the known cell cycle regulating tfs were studied. both chip-chip and microarray data are noisy, and we found no direct quantitative dependence of λ or r on the binding p-value. however, we obtained interesting results when using the target and non-target control groups as benchmarks. the receiver operating characteristic curve was used to determine whether or not λ can potentially discriminate targets from non-targets. roc is a graphical plot of the sensitivity versus 1-specificity, namely the fraction of true positives versus the fraction of false positives, as the discrimination threshold of a classifier is varied. the area under curve reflects the performance. the roc of a random classifier would be a 45° line with auc =. in figure the roc plot for each tf in alpha factor arrest experiment is given. overall, λ is significantly better than either a random classifier or r, suggesting that expression levels from tf and target genes are likely to exhibit phase locking. in contrast, the expression levels of tf and their targets seem to correlate little. in fact, using r does not fare any better than a random classifier. the results from other datasets are similar, and the auc of roc plots from all four datasets are summarized in figure the list of target genes that exhibit significant phase locking with the tf is given in additional file phase locking of biogrid gene pairs to further investigate the potential of phase locking analysis in network inference, we examined phase locking between gene pairs that have evidence of other types of interaction according to biogrid http://www.thebiogrid.org. biogrid is a freely accessible database of physical and genetic interactions, curated from high-throughput data and literature. of all possible gene pairs in our data sets, ~ are annotated in biogrid. we constructed five bootstrapping sets that consisted of the same number of biogrid gene pairs, randomly sampled from all possible pairs. the distribution of phase locking index was examined in each group. we found that the distribution for the biogrid pairs is skewed toward higher λ than the bootstrapping sets. the odds ratio of the enrichment of phase locked pairs at a given threshold value λ when there is prior evidence of interaction according to biogrid, can be calculated by: odds ratio=pp where p is the probability of a biogrid gene pair having a significant phase locking λ >λ, and p is such probability of a bootstrapping gene pair. the or of the four datasets at different λ is given in figure evidently, it increases with more stringent λ and is significantly greater than above the cutoff values of λ that we used in this study. these results further validate the idea that λ can be utilized to identify interacting gene pairs. phase locking among cell cycle genes we anticipate the expression of the cell cycle genes to show significant coordination, as they should all respond to the master signals that drive the cell cycle, and form an interaction network that together regulate the cell cycle progression. we found that indeed they exhibit phase locking among themselves more often than with random genes or than within random genes. on the other hand, when evaluated by the correlation coefficient r of their expression levels, cell cycle genes are not different from random genes. to illustrate the point we present in figure colormaps of the adjacency matrix a from the alpha factor arrest dataset, where the non-zero elements are presented by a black pixel. according to λ, in the subregion of cell cycle genes, there are significantly more black pixels than in the subregion of random genes, or in the subregion of cross interactions. to give a more quantitative description in figure we give the network degrees of all genes in the network. evidently, based on phase locking, the cell cycle genes form a highly connected subnetwork, whilst random genes are sparsely connected. in sharp contrast, according to r, there is literally no difference between the cell cycle genes and the random genes, both with moderate connectivity. results from the other three data sets are similar. note that in each case the total number of genes is less than their purported values, this is because we have removed genes with more than one missing values in their time series measurements. phase lags and m:n phase locking in general we found that the phase lag Ψ between locked pairs varies widely, and most of them are not close to this is consistent with observations made by others. figure gives the distribution of the relative phase Ψ as defined in equation of phase locked tf and target pairs. they span the whole range of, with no consistent pattern. this is not surprising as it is the protein of the tf that interact directly with the target gene transcription, not the expression of the tf transcript. indeed, it has been observed that response delay varies widely in gene expression regulation. additional file presents some examples of phase locked biogrid gene pairs with significant phase lags, where the correlations are low, even with time-lagged correlations. more examples and movies of phase lagged gene pairs are available at our website http://zen.dom.uab.edu:8080/phase/demo/. the universal existence of phase lags and their variation suggest the advantage of a method like the phase locking analysis that can automatically account for them. the pair-wise alignment approaches will need to align each pair individually thus adding an extra step in network modeling. in this study we included the consideration of m:n phase locking, up to max{m, n} = we found that there were a significant number of gene pairs exhibiting higher than 1: locking. the number breakdowns are listed in additional file in the cell cycle gene analysis, on average there are gene pairs in each data set that exhibit significant phase locking. out of these pairs, ~ % are from higher order locking, with ~half contributed from the 1: and 2: locking. in the tf-target pairs, the proportion of higher order locking is significantly higher at ~10%. on average, in each dataset there are significantly phase locked tf-target pairs. are due to higher than 1: locking, and ~70% of them are contributed from the 1: or 2: locking. the number of phase locked pairs, versus max{m, n}, is given in figure interestingly, an exponential dependence is evident. including higher order phase locking significantly improved the performance of λ as a classifier to detect tf targets. in figure we also included results that only considered 1: locking. clearly, the performance is not as good as when higher order locking is included, though it is still better than r or a random classifier. the basis of higher order phase locking is the different intrinsic dynamic time scales of the different time series. figure presents some examples of order locked biogrid gene pairs. the corresponding pubmed id that contains evidence of their interaction is also given. more examples of high order locked pairs from the cell cycle genes, and tf-target pairs are given in additional files & in these examples, the values of λ, and r are low, and interaction would have been missed if using them as the metric. demo movies of higher order phase locking are available at our website http://zen.dom.uab.edu:8080/phase/demo/. the wide range of phase lags and the significant proportion of higher order locking emphasize the need of a rigorous method like the phase locking analysis that can automatically take care of them. agreement between the datasets as observed by others, we see large variation among the datasets. however, when we examine the phase locked pairs, we found that there is significant concordance between the four datasets, with p <. the venn diagram is given in figure 10a. interestingly, we find that often for gene pairs that show consistent phase locking in multiple datasets, their correlation do not follow the trend. figure 10b is an example, where the tf-target pair swi5-ash, shows high λ consistently in the alpha, cdc and elu experiments, but their expressions correlate in no dataset. network inference using phase locking in our study so far, we have demonstrated that phase locking in expression changes is a good indicator of interaction. it is therefore natural to utilize it to construct gene interaction networks. highly connected genes tend to be essential genes we first examined the network defined by the adjacency matrix given in equation analysis of the yeast protein-protein interaction network in the past has revealed that highly connected genes are more likely to be essential for survival. in a recent study of co-expression network by zhang et al, strong positive correlations between the network degree and the functional essentiality of genes were also observed. in our phase locked networks, we also observed a positive correlation between the network degree ki of a gene and its likelihood of being essential. the relationship is nicely captured in figure 11a. all genes were ranked according to their ki and partitioned into equal sized bins. the proportion of essential genes in each bin was then plotted against the mean ki. a linear dependence is evident. we also examined the degree distribution of essential versus non-essential genes, given in figure 11b. compared to the non-essential ones, essential genes are significantly skewed toward having higher numbers of phase locked partners. genes in network modules have focused biological themes the structure of the phase locked network is further studied through hierarchical clustering using the topological overlapping matrix given in equation as depicted in figure modular organization is evident. genes in the grey regions did not belong to any modules. most of the modules contain higher numbers of essential genes than the scattered set. based on the theory of phase locking and our results so far, one anticipates genes from the same module to be highly interactive among themselves. this is indeed true for all the modules when we annotate them with biogrid information. in addition, ontological analysis of each identified module was performed using gostat, to examine the biological themes and the functional relationship of the module members. we found that genes in the same module tend to be involved in a set of focused signaling pathways. at a stringent cutoff p =, genes in the modules share a total of go categories, and at p =, whilst random gene sets of the same sizes and the scattered gene set do not share any go category even at a very loose cut off p =. in table we list the go biological processes shared by genes in each module at p <. a common theme is evident that most of these processes are critical to dna replication and cell cycle regulation. for instance, genes in the turquoise module likely form a cell cycle regulatory module, while the yellow module is likely associated with rna processing. discussion the importance of phase space information synchrony through oscillation is a common, maybe the most efficient, way to coordinate regulation in complex nonlinear dynamic systems. it is also a ubiquitous phenomenon in biological systems, where oscillations are observed in all organisms across a wide range of temporal and spatial scales, and are believed to play an important role in maintaining homeostasis and delivering encoded information. examples include the synchronized oscillations in interneuron networks, pulsatile endocrine hormone secretion, circadian oscillators, somite segmentation, and innumerable others. higher order phase locking occurs frequently, reflecting the multi-stability of complex systems, and is believed important to function. for instance, in the study of cardiorespiratory synchronization, when plotting the instantaneous respiratory phase at the occurrence of a heartbeat versus time, schäfer and co-workers found n: synchronization between heart and respiration. in a study with anesthetized rats, stefanovska et al. further observed lengthy synchronization epochs, and transitions from one ratio to another. they suggested that such transitions might be useful in monitoring depth of anesthesia. increasing evidence suggests that, as in many other complex systems in nature, oscillation and frequency modulation is also a general strategy for an organism to coordinate multi-gene responses to external signals. it is found that the transcription factor activities, rather than levels of transcription factor expression, mediate transcriptional regulations. in the negative feedback loop between the tumor suppressor p and the oncogene mdm p is expressed in a series of discrete pulses after dna damage, leading to oscillations in mdm. the amplitude of the oscillations was much more variable than the period, suggesting strong temporal regulation. in the nf-κb signaling pathway, nf-κb localization showed asynchronous oscillations following cell stimulation that decreased in frequency with increased iκbα transcription. transcription of target genes depended on its oscillation persistence, and thus the functional consequences of nf-κb signaling likely depend on temporal characteristics of the oscillations. in yeast cells it has been shown for several calcium stress responsive tfs that calcium concentration controls the frequency, but not the duration, of their oscillatory localization bursts, and the oscillation propagates to the expression of downstream genes. it has been argued that such frequency modulation of localization bursts ensures proportional expression of multiple target genes across a wide dynamic range of expression levels. these facts all imply the importance to study synchrony in expression oscillation, to understand the information encoded and the underlying interaction/regulation mechanisms. data from these studies also indicate that oscillations in gene expression, like most other oscillations in biological systems, are often pulsatile rather than harmonic. therefore, mathematical methods rooted from phase space analysis are desirable. the latter can potentially lead to new efficient network modeling algorithms, and help to understand the complex spatial-temporal architecture of transcriptome. advantages and limitations of the phase locking analysis additional to the theoretical appeals, we believe that phase locking analysis has several advantages in clustering genes of similar patterns and in network modeling. firstly, compared to approaches that primarily rely on the similarities in the amplitude domain patterns, phase locking utilizes the dynamics underlying the temporal pattern, which is more robust against noise. this is particularly appealing in network modeling of gene expression data, as they usually contain high noise. also the transcript abundance measurements often reflect a compressed, even altered representation of the true expression changes due to technical issues. these can significantly mask the true patterns in amplitude changes. in contrast, phase locking analysis, which focuses on the timing of the changes, will be less affected by the noise and the technical issues in the microarray gene expression study. in fact, it is known that noisy coupled nonlinear dynamic systems may synchronize in phase whilst their amplitudes remain uncorrelated. in our analysis, we have seen ample examples where interacting gene pairs exhibit obvious phase locking, but have very low correlation. secondly, in phase locking analysis, the phase lags in gene expression changes between different genes are automatically accounted for, and the performance is not affected by the amount of lag. on the other hand, the performance of alignment approaches depends on whether the lag is close to an integer number of the time steps of the experiment, and they need to be adjusted for each pair as the lags of different gene pairs vary greatly. thirdly, phase locking does not require the two time series to have the same dynamic time scale, or the same frequency. it is known that some pathways or gene groups in a cell respond to external signals at a much faster time scale than others. high order m:n phase locking analysis can take care of such interacting gene pairs, whilst they would be missed by the alignment method. a limitation of the phase locking analysis is its reliance on the temporal spectrum to accurately derive the instantaneous phase, which could significantly affect results when either the number of sampling points or the sampling frequency is too low. note that higher sampling frequencies are needed to detect high order phase locking. we have carried out a set of numerical simulations and observed significant deterioration in performance when the number of sampling points is reduced to ~ or lower. this limitation is not unique to the phase locking analysis. all pair-wise alignment approaches suffer from the same limitations as they all rely on adequate sample size to make a good assessment of whether the expression patterns of the pair are similar or not. there is also a caveat with the application of phase locking analysis to network modeling. gene interactions or phase locking occurs inside each cell. high throughput time series gene expression studies commonly measure a population of cells all at one time, effectively averaging the expressions of each gene across the whole cell population. in processes where there is high synchrony in the whole cell population, such as the cell cycle study presented here, phase locking between gene expression changes that occur inside each cell is preserved at population level, and can be detected from population measurements. in other biological processes, where cellular heterogeneity plays a key role, information of the signaling dynamics and phase locking inside each cell could be lost in population-level measurement. again, in such situations, the performance of the co-expression alignment approaches to detect interaction will also be affected. lastly, the recent advancements in single-cell techniques has enabled the generation of time-series gene expression measurements in a large number of individual living cells. we believe that phase locking analysis will be particularly suitable for such data. the dynamic information at the level of individual living cells will be critical to unravel how a genetic network operates at the systems level.  a major challenge in systems biology is to reconstruct gene networks that are involved in basic cellular processes, and to understand how alterations to the network affect functions and consequently phenotypes. interactions between genes can result in expression amplitude variations as well as temporal patterns. therefore, network inference utilizing temporal domain information deserves more attention. in this study, we investigated the potential of the phase locking analysis in network modeling of time course gene expression data. we demonstrated that interacting gene pairs, including transcription regulation interaction, protein-protein, or genetic interaction, are more likely to exhibit phase locking in their gene expression changes, and vice versa. among the phase locked pairs, up to ~10% are contributed from higher order locking, and the relative phase difference spans across the whole range of. based on these findings, we constructed interaction networks and revealed that genes with higher network degrees are more likely to be essential genes. utilizing the phase-locking index based topological overlapping matrix, we further investigated the modular structures in the network. we showed that genes forming network modules are more likely to be essential genes than scattered genes in the network, and members of the same module tend to be involved in the same biological functions and processes. in view of the importance of the frequency domain signal in transcription regulation, we believe that the phase locking analysis can potentially lead to new network modeling approaches and help to understand the dynamic designs of the intracellular signaling networks. 
in nonlinear dynamic systems, synchrony through oscillation and frequency modulation is a general control strategy to coordinate multiple modules in response to external signals. increasing evidence suggests that frequency modulation is also common in transcription regulation. interestingly, many of the phase locked gene pairs exhibit higher order than 1: locking, and significant phase lags with respect to each other. we show that it is efficient at identifying network modules of focused biological themes that are important to cell cycle regulation. using the yeast cell cycle data, we show that significant phase locking exists between transcription factors and their targets, between gene pairs with prior evidence of physical or genetic interactions, and among cell cycle genes. in addition, it can automatically address issues of arbitrary time lags or different dynamic time scales in different genes, without the need for alignment. when compared with simple correlation we found that the phase locking metric can identify gene pairs that interact with each other more efficiently. conversely, the synchrony information can be utilized to infer interaction. in this study, we investigate the potential of phase locking analysis, a technique to study the synchrony patterns, in the transcription network modeling of time course gene expression data. based on these findings we propose a new phase locking metric for network reconstruction using time course gene expression data.
 in nonlinear dynamic systems, synchrony through oscillation and frequency modulation is a general control strategy to coordinate multiple modules in response to external signals. conversely, the synchrony information can be utilized to infer interaction. increasing evidence suggests that frequency modulation is also common in transcription regulation.  in this study, we investigate the potential of phase locking analysis, a technique to study the synchrony patterns, in the transcription network modeling of time course gene expression data. using the yeast cell cycle data, we show that significant phase locking exists between transcription factors and their targets, between gene pairs with prior evidence of physical or genetic interactions, and among cell cycle genes. when compared with simple correlation we found that the phase locking metric can identify gene pairs that interact with each other more efficiently. in addition, it can automatically address issues of arbitrary time lags or different dynamic time scales in different genes, without the need for alignment. interestingly, many of the phase locked gene pairs exhibit higher order than 1: locking, and significant phase lags with respect to each other. based on these findings we propose a new phase locking metric for network reconstruction using time course gene expression data. we show that it is efficient at identifying network modules of focused biological themes that are important to cell cycle regulation.  our result demonstrates the potential of phase locking analysis in transcription network modeling. it also suggests the importance of understanding the dynamics underlying the gene expression patterns. 