A stochastic differential equation model for transcriptional regulatory networks

Background This work explores the quantitative characteristics of the local transcriptional regulatory network based on the availability of time dependent gene expression data sets. The dynamics of the gene expression level are fitted via a stochastic differential equation model, yielding a set of specific regulators and their contribution. Results We show that a beta sigmoid function that keeps track of temporal parameters is a novel prototype of a regulatory function, with the effect of improving the performance of the profile prediction. The stochastic differential equation model follows well the dynamic of the gene expression levels. Conclusion When adapted to biological hypotheses and combined with a promoter analysis, the method proposed here leads to improved models of the transcriptional regulatory networks.


Background
The production of independent sets of time courses of microarray data [1][2][3], obtained for the most studied eukaryotic organism Saccharomyces cerevisiae, improved the knowledge on the relationship between genes through the transcriptional process in the cell. The mechanism of the gene expression regulation is not entirely known, yet progress has been made by combining in silico approaches with the analysis of experimental data. In particular, contributions from a qualitative analysis realized by the recognition of specific promoter sequences, binding sites, and transcription factors are enhanced by quantitative studies obtained from microarray gene expression data [4,5]. The transcriptional regulatory network, built from thousands of genes, has a dynamical nature: the transcriptional program adapts itself to organismal development through the cell cycle, or as a response to changes in environment. In a systemic view the network architecture is potentially established by a qualitative analysis while quantitative methods address the main dynamical aspects -the network switches and the level of its parameters. This type of information may be obtained by processing from The Tenth Annual International Conference on Research in Computational Biology Venice, Italy. 2-5 April 2006 gene expression data that keep track of the variations in the experimental conditions and temporal modifications suited for the understanding of a particular transcriptional behavior.
A mathematical model for the processing of time dependent gene expression data has been sought to describe the dynamical aspects of regulation and to estimate the level of contribution for each transcriptional regulator in a succession of events. In this work we strengthen the model proposed in [6], by means of a novel pattern for the regulatory function. This model uses a SDE to describe the dynamics of the target mRNA expression level that reflects the actual knowledge about the stochasticity in gene expression, in a biological framework [7]. The drift term of the SDE depends on the regulatory rate of the target gene. The noise term is modeled by a Brownian Motion process which accounts for the superposition of small random factors that arise dynamically. The regulation rate is obtained as a linear combination of the regulatory functions of specific elements of the network. We propose a beta sigmoid function as the prototype of the regulatory function, designed to keep track of the local temporal patterns of the target gene regulators.
Our analysis shows that the utilization of the beta sigmoid function enhances the results in [6] where sigmoid functions were considered. The comparison was made by applying the model to the same test data set as in [6], given by gene expression measurements of the mRNA levels of 6178 S. cerevisiae ORFs at 18 time points under the α factor synchronization method [1]. A candidate pool of potential regulators was constructed by joining transcription factors, cell-cycle control factors and DNA-binding transcriptional regulators as found in the literature [1,8]. We performed the same statistical analysis from [6] based on the maximum likelihood principle for the estimation of the model parameters. The AIC strategy was used for the selection of the best fitting combination of the pool regulators. With the addition of beta sigmoid pattern, the SDE model renders good prediction results even in the case of the previously worst fitted genes obtained by [6]. The procedure proposed herein may be well suited to quantify transcriptional regulatory networks, provided it is tailored to the characteristics of the input data set.

Results
The model was evaluated on the data set from [1] that provides gene expression measurements of the mRNA levels of 6178 S. cerevisiae ORF s at 18 time points under the α factor synchronization method. To compare our results with those from [6], we used the data set of 216 potential regulator listed at [9], constructed by joining transcription factors, cell-cycle control factors and DNA-binding transcriptional regulators described in the literature [1,8,10]. This set has been created with respect to the regulation of the cell cycle process. There are about 800 genes identified to be involved in the cell cycle of the budding yeast [1] and we performed our analysis on the entire data set. This fact did not carry any methodological artifact because the target genes are processed independently. The advantage is that good prediction results may imply new hypothesis about the regulators of particular genes. The output of our analysis is bipartite. For each gene we provide 1. the parameters of the goodness of fit: log likelihood (log L), AIC and QE of the predicted mRNA levels with respect to the observed values 2. the corresponding regulators with their regulatory effect expressed by the local network weights; positive weights correspond to activator genes and negative weights correspond to repressor genes.
The full output of our analysis is presented in Additional file 1. The beta sigmoid function was non-degenerate for 72% of the expression values. Tables 1 and 2 show that the use of the beta sigmoid model for the regulatory function improves the fitting parameters for 15 of the 20 genes depicted in [6] as worst fitted. The prediction of five genes (YGR269W, FIT3, HSH49, ASH1 and ATS1) shows a substantial amelioration. Over the entire data set, the novel model of regulatory function improved the prediction of 29% of the gene expression profiles. The distribution of the improvement is presented in the histogram in Figure  1 computed for the difference between the quadratic error of the predicted mRNA levels with the sigmoid model and the quadratic error of the predicted mRNA levels with the beta sigmoid model. We note that there is a non-negligible number of genes with their expression level better fitted by the beta sigmoid regulatory pattern. The prediction results for a selection of 10 genes in this category are shown in Table 3 and Table 4. • there exist genes for which the SDE model with beta sigmoid regulatory function gives a better prediction than the SDE model with sigmoid regulatory function.
A good quality of fitting of a particular gene allows the consideration of the regulators associated by the model for further investigation such as DNA-binding sites or promoter architectures. The quadratic error of prediction with beta sigmoid regulatory function is less than 0.5 for 1885 genes from the entire data set (see Additional file 1).

Discussion and conclusions
The global view of the regulatory network is a cascade model, with genes regulating genes regulating other genes at their turn [11]. The SDE model [6] revisited here addresses the network local connections, i.e. the strict neighborhood of one target gene. The drift term of the stochastic differential equation is given by the regulation rate which quantifies the local network architecture by a linear combination of regulatory functions of regulating genes. The choice of the regulatory function pattern is a central aspect of the model, since the fitting of the gene expression profiles is very sensitive with respect to the drift term of the SDE. This model has the ability to extract from a given set of potential regulators those that fit the target gene expression profile.
The prototype of regulatory function introduced in [6] has a sigmoid pattern, built on the statistical characteristics of mRNA expression levels -see Equation (14). By keeping track of the temporal pattern of regulation, we show that the prediction of target gene expression profiles is improved for 29% of genes tested. We propose a prototype of regulatory function supported by a beta sigmoid model, built on temporal parameters extracted from the expression profiles of the regulators -see Equation (13). The SDE method relies on the assumption that the best fit of the target expression profiles is informative for the identification of the regulators and of their contribution. Thus, for the study of a specific set of target genes, our prototype of regulatory function may give more accurate results and provides a switch for the model proposed in [6]. Conceptually the beta sigmoid model has the advantage to correspond to the biological process of regulation: the temporal window of the peak defined by the shape of the beta sigmoid function reflects features of the regulation mechanism.
The regulation of gene expression in eukaryotes is a complex phenomenon and various particularities from one type of gene to another may occur. Hence the regulatory pattern may vary from gene to gene [4]. This fact is revealed in our result which shows that there are genes for which we can choose the best model between the beta sigmoid and the sigmoid pattern while for other genes neither of them fits the data. Before reaching this conclusion one has to be aware about the limitation induced from the selection of the set of potential regulators since incomplete information at this level may deteriorate the results.
Further research on more complex and explicit regulatory functions are foreseen from the availability of data sets and studies on various experimental condition for the budding yeast (sporulation [3], diauxic shift, heat and cold shock, treatment with DTT, pheromone and DNAdamaging agents [12]). In this framework a challenging task could be the study of the existence of possible rela- Fitting parameters obtained with the beta sigmoid model for the set of the 20 worst fitted genes with the sigmoid regulatory model; the asterisk * marks an improvement of the fitting with respect to the results from Table 2. The prediction of five genes (YGR269W, FIT3, HSH49, ASH1 and ATS1) shows a significant improvement. The fitting parameters from the sigmoid model of regulatory functions of the genes from Table 3. Genes fitted by the SDE model with beta sigmoid as regulatory function. The set of the worst fitted 20 genes by the sigmoid model, sorted in the increasing order of the log-likelihood.
tionships between the type of regulation pattern and the gene specificity.
This work provides a second implementation of the algorithm based on the SDE model, enlarged with a new type of regulatory pattern. The predictions from the algorithm may be improved with better strategies for the selection of the candidate pool of regulators. Moreover, the algorithm is a potential tool for the investigation of the interactions between the regulators of a target gene, modeled with a drift term defined by a non linear combination of regulatory functions.
This study shows that the SDE framework constitutes a reliable tool for the analysis of the transcriptional regulatory networks, provided it is completed with a validation of the identified regulators by a promoter analysis.

SDE model of time-continuous gene expression data
Let T denote a discrete set that corresponds to the time instants of the gene expression measurements. Consider two stochastic processes defined for a given target gene, (N t ) t∈T and (X t ) t∈T that model, respectively, the variation The model proposed in [6] assumes that from time t to Δt the transcription and degradation process are given by where (W t ) t∈T is a Brownian Motion process that models the random error and σ is a positive scaling parameter.
Consider infinitesimal time intervals, that is Δt → 0; from this it follows that the relation in Equation (2) becomes a stochastic differential equation Since N t is proportional with the signal intensity S t , and X t = log(S t -B) -where B is the background intensityassume without loss of generality that Comparative plot between the observed and the predicted values of mRNA expression levels of gene YALO61W The regulatory network is represented locally, in the neighborhood of the target gene, as a superposition of regulatory elements pictured in Figure 5. The local network relationship is modeled by the regulatory rate function, built from the observable information, i.e. the regulators mRNA expression levels, as: where F i denotes the regulatory functions of the potential regulators from . The constants c 0 , c 1 ,...,c m are the learning parameters of the network; they modulate the network behavior and carry information about the local regulatory process.

Beta sigmoid pattern of regulation
The regulatory function is a key element of the model and fits the quantitative pattern with a specific regulator that acts on the mRNA expression of the target gene.
Our work revealed a prototype of the regulatory function based on the beta sigmoid function, given by

Regulator functions
Learning network constants genes Regulator Noise .
. . For the evaluation of the impact of the beta sigmoid regulatory function model, the network weights are estimated from gene expression data with the standard statistical procedure described in detail in [6]. Equation (15)  The criteria used for the selection of the regulators is AIC [14]. Between any two combinations of regulators, the best combination is that for which the AIC of the regulators has the smallest value. The computation of AIC follows from AIC = -2 + 2(m + 1) where is the estimator of log L and is obtained from the functional invariance property of the maximum likelihood estimators and , i.e., = log L( , ).
Let denote the set formed by a candidate pool of regulators of the target gene; denote by | | the cardinality of . Ideally, ML and AIC procedures shall be performed on each combination of regulators from . Since the number of all possible combinations of regulators is , an enumeration algorithm for those sets will explode quickly. The heuristic procedure used is the forward selection strategy [15]. At first the regulator with the biggest log-likelihood with respect to the target gene is selected. A new regulator is added if it will increase the AIC more than any other single regulator outside the current combination. The actual implementation stops for a combination of maximum 10 regulators, exactly as done in [6]. Under these conditions the performance of the algorithm we propose is expressed by an order of magnitude equal to O(nm 2 ). In practice this is a slight enhancement compared to the algorithm proposed in [6] for which the order of magnitude equals O(n 2 m 2 ) -since for actual experimental data the number of time courses n is quite small. The difference in the performance of the two algorithms comes from the fact that the search of the maximum is less costly than the computation of the statistical parameters for a data set.