Quantitative inference of dynamic regulatory pathways via microarray data

Biological phenomena at different organismic levels have revealed some sophisticated systematic architectures of cellular and physiological activities implicitly. These architectures were built upon the biochemical processes before the emergence of proteome and transcriptome [1–3]. Under the molecular machinery, the biochemical processes are mostly interpreted as frameworks of connectivity between biochemical compounds and proteins, which are synthesized from genes to function as transcription factors binding to regulatory sites of other genes, as enzymes catalyzing metabolic reactions, or as components of signal transduction pathways [4–6]. This implies that, in order to understand the molecular mechanism of genes in the control of intracellular or intercellular processes, the scope should be broadened from DNA sequences coding for proteins to the systems of genetic regulatory pathways determining which genes are expressed, when and where in the organism and to which extent [7]. In the experience of engineering field, the systematic architecture and dynamic model could investigate the characteristics of signaling regulatory pathways [8]. Therefore, how to construct the dynamic model of a signaling pathway from the system structure point of view might be the first key to the door of system biology. Most biological phenomena directly or indirectly influenced by genes such as metabolism, stress response, and cell cycle are well studied on the molecular basis. Thus, identification of a signal transduction pathway could be traced back to the genetic regulatory level. The rapid advances of genome sequencing and DNA microarray technology make possible the quantitative analysis of signaling pathway besides the qualitative analysis. More particularly, the embedded time-course feature of microarray data would promote the system analysis of signal regulatory pathways as well, which is very mature in the field of engineering.

In addition to northern blots and reverse transcription-polymerase chain reaction (RT-PCR), which study a small number of genes in a single assay, the transcriptome analysis has, via DNA microarray technology [9], managed to achieve high-throughput monitoring of the almost genome-wide mRNA expression levels in living cells or tissues. Two types of available microarrays, the spotted cDNA and in situ synthesized oligonucletide [10] chips, which permit the spatiotemporal expression levels of genes to be rapidly measured in a massively parallel way, are used in different experimental requirements and stocked in the databases on net, such as Stanford Microarray Database (SMD) [11], Gene Expression Omnibus(GEO) [12] in NCBI, and ArrayExpress [13] in EBI. Microarray experiments are now routinely used to collect large-scale time series data that facilitate quantitative genetic regulatory analysis while qualitative discussion is the traditional thinking [14–17].

Several analytic methods have been proposed to infer genetic interrelations from gene expression data. In the coarse-scale approach of clustering, the underlying conjecture is that co-expression is indicative of the co-regulation, thus clustering may identify genes that have similar functions or are involved in the related biological processes. The most widely used method is the unsupervised hierarchical clustering [18]. This approach has an increasing number of nested classes by similarity measurement and resembles a phylogenetic classification. If we know the number of clusters in advance, the k-means clustering [19] could assign gene elements into a fixed number k of clusters in a way to minimize the overall Pearson or Euclidean distances of each member internally in the same cluster. Other algorithms such as the neural-network-based self-organizing maps (SOM) [20], singular value decomposition (SVD) or principal component analysis (PCA) [21], and fuzzy clustering methods [19] also have their own advantages and limitations. Alternative supervised clustering algorithm of support vector machine [22], which uses prior biological information of cluster for training, would enhance the accuracy of clustering. However, the nature of clustering algorithms apparently cannot uncover the causal interactions between genes just by grouping. Regarding the causality of pathways, the clustering analysis needs to cooperate with sequence motif detection [23]. It is also important to note that models using clustering analysis are static and thus can not describe the dynamic evolution of gene expression, even in the type of time-course microarray data.

A statistical model of Bayesian network [24] was proposed to model genetic regulatory networks. Basically, the technique uses a probabilistic score to evaluate the networks with respect to the expression data and searches for the network with the optimal score. The dynamic Bayesian network [25] was proposed to learn the network structure and parameters by maximizing the posterior probability via Bayes rule of prior probability and marginal likelihood. Another algorithm of Boolean networks [26] can also be employed to model the dynamic evolution of gene expression, where the state of a gene can be simplified to being either active (on, 1) or inactive (off, 0). The probabilistic nature of Bayesian networks is capable of handling noise inherent in both the biological processes and the microarray experiments. This makes Bayesian networks superior to Boolean networks, which are deterministic in nature. The validity of dynamic Bayesian networks is evaluated by the sensitivity-specificity score ratios [25], which depend on the training size, the degree of accuracy of prior assumption. A genetic regulatory network based on the first order differential equation with given decay rates was discussed in [27].

In this study, the dynamic system approach could be employed to model how a target gene's expression profile is regulated by its upstream regulatory genes from the system causality point of view. Then, with the causal dynamic model, the upstream regulatory function can be extracted from the expression profile of the target gene by the optimal estimation method, i.e. maximum likelihood estimation. Since merely the second-order differential equation is employed to model the dynamic evolution of the target gene, only a few parameters need to be estimated. Furthermore, the derived regulatory function is closely related to the causal upstream information of the pathway and will create a basis for inferring the regulatory pathway from the system biology point of view.

In either eukaryote or prokaryote, signaling regulatory pathways are considered as responses to the physiological activities or the deviation from homeostasis, which would affect the normal states of an organism. Among these signaling regulatory pathways, cell cycle [17] is one of the most conspicuous features of life which plays an important role in growth and cellular differentiation in all organisms. In plants, the stress-induced pathways [28] are very important to survivability under the abiotic environmental treatment such as drought, salinity and cold[29]. If these critical pathways can be identified from quantitative analysis in silico, the defect of biological processes would be predicted and corrected before hand. Our aim is to construct signaling regulatory pathways quantitatively by the system inference approach with a dynamic model and microarray data.

In this study, a second-order differential equation, which has been widely used to model many physical dynamic systems with good characteristics, is proposed to model the time-profile evolutional behavior of a target gene. The regulatory function is taken as the driving input of the dynamic equation of the target gene. Using the dynamic equation and microarray data, we first extract the regulatory function for each target gene. According to the extracted regulatory function, we deduce their upstream regulators to trace back upstream signaling pathways. Then, upstream regulatory genes are taken as target genes to trace back their upstream regulatory genes. Iteratively, we can construct the whole regulatory pathway to the genome wide using the dynamic regulatory model and microarray data from the system biology point of view. Finally, we give some independent validation of our approach by repeating the analysis with randomly reshuffling the time order of microarray data and see if the proposed pathways are destroyed.

We have applied our dynamic system approach to two genetic regulatory pathways with microarray data sets publicly available on net [15, 30]. One is the circadian regulatory pathway in Arabidopsis thaliana [31, 32], and the other is the metabolic shift pathway from fermentation to respiration in yeast Saccharomyces cerevisiae [33]. The circadian system is an essential signaling pathway that allows organisms to adjust cellular and physiological processes in anticipation of periodic changes of light in the environment [34–38]. According to the synchronously dynamic evolution of microarray data, we have successively identified the core signaling transduction from light receptors to the endogenous biological clock [39, 40], which is coupled to control the correlatively physiological activity with paces on a daily basis. On the other hand, the diauxic shift [41] from the exhausted fermentable sugar of anaerobic metabolism to aerobic growth is correlated with widespread changes in the expression of genes involved in fundamental cellular processes such as carbon metabolism, protein synthesis, and carbohydrate storage. [28, 31, 32, 42–47] The architecture of the signaling pathway correlative to glycolysis or gluconeogenesis during the diauxic shift is properly built up. With the dynamic system approach, not only the regulatory abilities between causal genes could be derived, but also the delays of regulatory activity are specified. These quantitative characteristics will help determine the intrinsic frameworks of connectivity in the above interesting pathways from the system biology point of view.

The proposed methods in this study would be divided into four steps. In the first step, a dynamic model using the second-order differential equation is developed to describe the expression profile data as output and the regulatory function as input to denote the implicit characteristics of each gene with some parameters. With the help of the second-order dynamic model, we would then extract the upstream regulatory function from the expression profile of the target gene using the optimal estimation method. In the third step, the regulatory function estimated will help seek the correlative regulatory signals from the upstream paths. Iteratively, we can reconstruct the whole signaling regulatory pathway by linking up the upstream regulatory paths. Finally, some biological filters using available biological knowledge are employed to prune the constructed signaling regulatory pathway to improve the accuracy of the proposed method.

I. Dynamic system description of signaling regulatory model

The second-order differential equation is well used in the description of dynamic system evolved from the causality of gene regulatory function. Let X _i (t) denote the expression profile of the i-th gene at time point t. The following second-order differential equation is proposed to model the expression level of the i-th gene,

where G _i (t) is the upstream regulatory function to influence the expression profile X _i (t) of the i-th gene while a _i , and b _i are the parameters that characterize the dynamic inherent property of the gene like degradation and oscillation, and ε _i (t) is the noise of current microarray data or the residue of the model. In general, the second-order differential equation has been widely used to model dynamic systems to characterize efficiently the dynamic properties of damping and resonance of systems in physics and engineering.

Obviously, the clue of upstream regulatory pathways is in G _i (t). Thus, the first step is to detect the upstream regulatory function G _i (t) from both dynamic equation in (1.1) and microarray data. However, to detect the input regulatory function G _i (t) from both equation (1.1) and microarray data directly is not easy. In this situation, a Fourier decomposition technique is employed to decompose G _i (t) as a synthesis of some harmonic sinusoid functions so that the signal detection problem of G _i (t) is reduced to a simple parameter estimation problem.

Accordingly, we can decompose G _i (t) by the following Fourier series,

Then the detection of G _i (t) becomes how to estimate the Fourier coefficients of α _n and β _n , which are the magnitudes of different harmonics of cos(nωt) and sin(nωt), for n = 0,..., N in equation (1.2), respectively. In science and engineering, the Fourier series has been widely employed to synthesize any continuous functions with finite energy. The estimation of α _n , β _n and the detection of G _i (t) in equation (1.2) are given in Methods in the sequel.

As a result of parameter estimation in Methods, the detection

of regulatory function G _i (t) could be derived as follows,

Since the input regulatory function G _i (t) of a target gene is usually due to the transcriptional binding or some physical interactions from the upstream regulatory genes, in the following, we would trace back to the corresponding regulatory genes from input regulatory function of the target gene.

II. Inference of the regulatory pathway via

Apparently the input regulatory function G _i (t) in equation (1.1) contains the driving information for the target gene's expression from the upstream regulatory genes. The identified regulatory function from equation (1.3) could be interpreted as the regulatory connectivity through transcriptional binding or protein-protein interaction imposed on the i-th target gene. Nevertheless, the expression data of protein type which should be considered directly in practice are by now unavailable and unreliable to trace back upstream regulatory genes. Instead, the expression data on mRNA level which is now widely available from microarray assays would make tracing back the upstream regulatory pathway possible under proper assumptions. All along the paper we assume that the expression levels of mRNA transcripts are proportional to the actual number of corresponding proteins in the cell. This assumption is indeed a strong approximation since post-transcription is known to play a very important role in down regulating the number of the transcription factor in the cell.

Before the inference of upstream regulatory genes, it is reasonable to confine the effect of the regulatory genes on the regulated target gene. The saturated activity of expression level reveals that the regulatory ability cannot extend unlimitedly. The sigmoid function is often chosen to express the nonlinear saturation with proper parameters. Here, we apply the sigmoid transformation to represent the 'on' and 'off activities of the regulatory genes on binding or not to motifs of the target gene. So the regulatory signal

shown below with the parameter set of θ _j = {γ, M _j , τ _j } is the sigmoid transformation of X _j (t), the expression profile of the j-th regulatory gene.

where γ is the transition rate, M _j is the mean expression of the j-th regulatory gene's profile, and τ _j is the corresponding signal transduction delay.

The delay activity should be considered in order to describe the signal transduction delay τ _j from the j-th regulatory gene to the target gene. The delay τ _j would be computed by statistical correlation between the regulatory signal transformed from the j-th regulatory gene and the identified regulatory function of the target gene. The delay τ _j is determined by the following maximum correlation criterion,

where r _τ is the correlation between and under variable delay τ. If there are many τ to achieve the maximum correlation in (2.2), then only the smallest one is chosen.

Using the correlation method, we trace back R _i regulatory genes whose regulatory signals are most correlated with the regulatory function of the i th target gene, i.e. choose R _i genes with maximum correlation but with smaller τ _j in (2.2). The determination of number R _i will be discussed later. Then, we construct the regulatory pathway by tracing back R _i regulatory genes from the identified regulatory function of the target gene as the following kinetic relationship,

where c _ij is the pathway kinetic parameters from the regulatory gene j to the target gene i, R _i are the searched upstream regulatory genes, the constant c_{i 0}represents the basal level to denote the regulatory function other than upstream regulatory genes, and e _i (t) is the residue of the model.

Furthermore, to estimate the pathway kinetic parameters c _ij , equation (2.3) for m time points should be written in the following regression form,

where

, .

We assume that each element in the error vector, e _i (t _k ), k = {1,..., m}, is an independent random variable with a normal distribution with zero mean and variance σ². By maximum likelihood parameters estimation method (see Methods), the estimates of σ² and Ω _i are given as follows, which is solved as

and

It should be noted that with the combination of biological knowledge about the transcriptional factors, protein phosphorylation, post-transcriptional and specific enzyme regulation of target genes, lots of putative and verified genes correlated with the target genes are pruned by this biological filter for the more efficient and accurate searching of R _i upstream regulatory genes in equation (2.3). For example, suppose the expression profile of gene j has a high correlation with regulation function G _i (t) of target gene i. However, if gene j is not a transcription factor, protein phosphorylation, post-transcription or specific enzyme of target gene i, it will be deleted from the candidates of R _i upstream regulatory genes because it may be only a co-expressed gene with the target gene corregulated by the other gene. On the contrary, a verified regulatory gene should be recruited into the candidates even with small correlation with G _i (t).

Finally, we take the well-known Akaike Information Criterion (AIC) into account for determining the number R _i of regulatory signal [42],

The first term in AIC is the residual variance and the second term R _i is the number of regulatory genes. AIC includes both the estimated residual variance and model complexity in one statistic, which decreases as σ² decreases and increases as R _i increases. AIC has been widely employed to determine the complexity of system modeling science and engineering [42]. The optimal number R _i of the upstream regulatory genes will be determined by the minimization of the AIC value in equation(2.7).

Now, for the selected target genes in the interesting pathway, we could search for the optimal R _i upstream regulatory genes by AIC in equation (2.7) after the biological filtering and determine their pathway kinetic parameters c _ij of regulatory signal by equation (2.6). After biological filter pruning, if the number of candidates of regulatory genes is still less than R _i determined by AIC, then some genes, which are highly correlative to G _i (t) but not of transcription factors or signaling proteins of target gene i, should be recruit into candidates to uncover regulatory relationships that were not suspected to be connected. After the combination of equations (2.3) and (1.1), the whole regulatory pathway is obtained as

for i = {l, 2,..., L}, and L is the number of target genes in the pathway. The sub-paths related to the i-th target gene in the interesting pathway could be detected by the inference algorithm. Then, it is natural that the whole regulatory pathway would be constructed by the links of all the sub-paths. We also outline the whole flowchart of our dynamic inferring algorithm as shown in Figure 1 for an overview.

Microarray expression analysis by the dynamic system approach offers an opportunity to generate functional regulation interpretation on the genome-wide scale. The crucial ontology behind using dynamic system techniques is that the causality between gene expression profiles could be identified according to the differential equation underlying a dynamic system. Therefore, because the microarray data were harvested with time progression, the simultaneously varied gene expressions implicated in a genetic regulatory system would be detected to infer the regulatory pathways in spite of the versatile interactions such as transcriptional control, protein phosphorylation, or specific enzyme regulation.

The clustering method answers the problem of what is the functional catalogue of a specific gene by the identification of resembling patterns of gene expressions. Similarly, the co-regulations of upstream genes in our method also imply their concurrent functions. In contrast to the clustering algorithm, the causality of time-course data has been smoothly drawn by our dynamic method. The Bayesian networks were used merely for forward probabilistic estimation with time transition lacking in the feedback linkages. This unidirectional problem would not happen in our algorithm. Owing to the quantitative regulatory abilities of our model, we have a greater diversity of regulatory influence than the Boolean networks, which are deterministic with merely two states.

In our dynamic system approach, we not only can link qualitatively the upstream genes to the downstream ones iteratively, but also indicate quantitatively their regulatory relationships, including the regulatory abilities and the activation delays. In terms of the regulatory abilities, the comparison between the upstream regulatory genes of a target gene can inspire us to ask which one is significant biologically and whether it is a positive or negative influence on the investigated gene. Moreover, the speculation of activation delays benefits the empirical reference by providing us when the upstream regulatory genes might interact with their target genes. Since any gene can be considered as a target gene to trace back its upstream regulators, these regulators are then considered as target genes to trace back their upstream regulators. Iteratively, the genetic regulatory pathway (or network) can be constructed to the genome-wide. According to the qualitative and quantitative features imbedded, two regulatory pathway examples are characterized as in Figure 5 and Figure 6 for the identification of the proposed method. In addition, using the Akaike Information Criterion (AIC), a proper number of regulatory genes would be affirmed. As a result, many links overlap with well-known regulatory and signaling pathways in the previous literature and several putative ones are also found. Furthermore, the activation or repression relationships inferred via the microarray data would distinctly uncover the overall effect of regulatory interactions among casual genes in pathways on the transcriptional level.

In the two pathways under investigation, we have a more detailed understanding about the regulatory interactions among phytochrome, crytochrome and biological clock in the circadian regulatory system. On the other hand, the sophisticated knowledge of the metabolic pathway after the diauxic shift can be unfolded properly in our analysis. Furthermore, the independent validation of our approach is also given by randomly reshuffling the time order of microarray experiment. We found that the proposed pathways in Figures 5 and 6 are all destroyed as shown in Figures 7 and 8, respectively. The successful analysis of these two pathways implies the development of a valid and high-throughput method. All of the programs have been released [see Additional file 1]

There are some shortcomings in our study. First, although the time-course microarray data are available, its lower samplings will distort the real changes of gene expressions, especially for quick dynamic evolution. A more sampling experiment with respect to the intrinsic turnover rate is expected to have more precise analysis. Secondly, a regulatory gene with larger activation delay would not be recognized because the less activation delay criterion is used, but this might be overcome by properly relaxing the criterion. Thirdly, activation profiles under the proteome should be highly correlated with the transcriptional profiles to elevate the interpretation of our system model. In general, the synchronous time-course microarray assay is more suitable to underlie the transcriptional binding among causal genes, but an inference of physical interactions in the post-transcriptional level also has sufficient feasibility in our study.

In the near future, the most pressing task is to investigate our presumed paths in the laboratory. As the pathway construction algorithms are further developed, we expect this system approach to have immense impact in elucidating the underlying molecular mechanisms of pathways in a variety of organisms, especially after the maturation of the protein chips. Ultimately, we envision that biologists will perform routine pathway inference to seek some novel regulations and to identify the evolutionarily conserved links.