 Research
 Open Access
 Published:
A stochastic differential equation model for transcriptional regulatory networks
BMC Bioinformatics volume 8, Article number: S4 (2007)
Abstract
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–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 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, cellcycle control factors and DNAbinding 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, cellcycle control factors and DNAbinding 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 nondegenerate 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 nonnegligible 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. Figure 2, Figure 3, and Figure 4 show comparative plots of the expression profiles for observed data and predictions from the stochastic differential equation model, with beta sigmoid and sigmoid prototypes of regulatory function. Two conclusions may be immediately drawn:

the SDE model can provide very good predictions of mRNA expression levels;

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 DNAbinding 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 relationships 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.
Models and methods
SDE model of timecontinuous 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 in time of the target gene amount of mRNA and the variation in time of the expression level of mRNA. Let \mathcal{R} be the set of potential regulators for the target gene. Denote by g_{ t }the function that models the transcription rate of the target gene at time t
gt : \mathcal{\text{P}}(\mathcal{R}) → \mathcal{R}_{+} (1)
where \mathcal{\text{P}}(\mathcal{R}) is the set of all possible subsets of \mathcal{R} and \mathcal{R}_{+} is the set of real positive numbers. Denote the real, positive mRNA degradation rate by λ.
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 intensity – assume without loss of generality that
X_{ t }= log(N_{ t }) (4)
Thus, the chain rule of the stochastic calculus applies (Itô formula) and the SDE obtained for X_{ t }yields
Local regulatory network
Consider an increasing sequence of temporal values
T = {t_{0} <t_{1} <...<t_{ n }} (6)
Let m be the cardinality of the set \mathcal{R} and let {X}_{t}^{i} be the mRNA expression level of the ith regulator from the set \mathcal{R}, measured at time t ∈ T. Denote
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 \mathcal{R}. 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
Figure 6 shows an example of beta sigmoid function shape. The parameter {t}_{s}^{i} corresponds to the point where the regulator expression level begins to increase. The maximal contribution of the regulator i is induced in the target gene at time {t}_{m}^{i}, when the mRNA expression level of the regulator attends its maximum, corresponding to the biological hypotheses. The beta sigmoid function degenerates after time 2{t}_{m}^{i}  {t}_{s}^{i} and becomes noninformative. For this reason we define the regulatory function as
where I_{ A }is the indicator function of the set A (I_{ A }(x) = 1 if x ∈ A and I_{ A }(x) = 0 if x ∉ A) and
is the sigmoid function; μ_{ i }and σ_{ i }are the mean and deviation of {\underset{\xaf}{X}}^{i}, the prototype of the regulatory function from [6].
The learning in the local network is driven by the SDE
where \tilde{c}_{0} = c_{0}  λ  σ^{2}/2. The network weights c_{1},...,c_{ m }carry information in both their magnitude and sign: positive values correspond to regulators with activation, and negative values correspond to repression.
Statistical analysis
For a given target gene, the aim of the statistical analysis is to extract from the time course mRNA levels

1.
the set of m regulators (model selection);

2.
their corresponding parameters σ and the set {\tilde{c}_{0}, c_{1},...,c_{ m }} of parameters estimation;
with the best fit with respect to Equation (15). The beta sigmoid as regulatory function adds supplementary parameters {t}_{s}^{i}, {t}_{m}^{i} and x_{ max }to the model. These parameters are estimated from the corresponding time course mRNA levels according to their definitions given in Equation (10) and employed in the computation of the estimators of σ and \underset{\xaf}{c}:
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) is considered in discrete form for each time interval [t_{ j }, t_{j+1}], j = {1, 2,...,n} that corresponds to time measurements. The estimators of σ and {\tilde{c}_{0}, c_{1},...,c_{ m }} are obtained maximizing the loglikelihood function log L (ML approach [13]) of the ndimensional random vector with elements
The computation of log L uses basic properties of Brownian Motion: the increments {W}_{{t}_{j+1}}{W}_{{t}_{j}} are pairwise independent and each increment is normally distributed, with zero mean and standard deviation given by \sqrt{{t}_{j+1}{t}_{j}}.
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\widehat{\mathrm{log}L} + 2(m + 1)
where \widehat{\mathrm{log}L} is the estimator of log L and is obtained from the functional invariance property of the maximum likelihood estimators \widehat{\sigma} and \widehat{\underset{\xaf}{c}}, i.e.,
= log L(\widehat{\sigma}, \widehat{\underset{\xaf}{c}}).
Let \mathscr{H} denote the set formed by a candidate pool of regulators of the target gene; denote by \mathscr{H} the cardinality of \mathscr{H}. Ideally, ML and AIC procedures shall be performed on each combination of regulators from \mathscr{H}. Since the number of all possible combinations of regulators is {2}^{\left\mathscr{H}\right}, 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 loglikelihood 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.
Availability
The method was implemented in R 2.2.1 (R Development Core Team, http://www.rproject.org/). The source code is available upon request.
Abbreviations
 SDE:

Stochastic Differential Equation
 AIC:

Akaike Information Criteria
 ML:

Maximum Likelihood
 QE:

Quadratic Error
References
Spellman PT, Sherlock G, Zhang MQ, Iyer V, Anders K, Eisen MB, Brown PO, Botstein D, Futcher B: Comprehensive identification of cell cycleregulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization. Mol Biol Cell 1998, 9: 3273–3297.
Cho RJ, Campbell MJ, Winzeler EA, Steinmetz L, Conway A, Wodicka L, Wolfsberg TG, Gabrielian AE, Landsman D, Lockhart DJ, Davis RW: A genome wide transcriptional analysis of the mitotic cell cycle. Mol Cell 1998, 2: 65–73. 10.1016/S10972765(00)801148
Chu S, DeRisi J, Eisen M, Mulholland J, Botstein D, Brown PO, Herskowitz I: The transcriptional program of sporulation in budding yeast. Science 1998, 282: 699–705. 10.1126/science.282.5389.699
Pilpel Y, Sudarsanam P, Church GM: Identifying regulatory networks by combinatorial analysis of promoter elements. Nat Genet 2001, 29: 153–159. 10.1038/ng724
Garvie CW, Wolberger C: Recognition of specific DNA sequences. Mol Cell 2001, 8: 937–946. 10.1016/S10972765(01)003926
Chen KC, Wang TY, Tseng HH, Huang CY, Kao CY: A stochastic differential equation model for quantifying transcriptional regulatory network in Saccharomyces cerevisiae. Bioinformatics 2005, 21: 2883–2890. 10.1093/bioinformatics/bti415
Kaern M, Elston TC, Blake WJ, Collins JJ: Stochasticity in gene expression: from theory to phenotypes. Nat Rev Genet 2005, 6: 451–464. 10.1038/nrg1615
Harbison CT, Gordon DB, Lee TI, Rinaldi NJ, Macisaac KD, Danford TW, Hannett NM, Tagne JB, Reynolds DB, Yoo J, Jennings EG, Zeitlinger J, Pokholok DK, Kellis M, Rolfe PA, Takusagawa KT, Lander ES, Gifford DK, Fraenkel E, Young RA: Transcriptional regulatory code of a eukaryotic genome. Nature 2004, 431: 99–104. 10.1038/nature02800
Transcriptional regulatory networks[http://www.csie.ntu.edu.tw/~b89x035/yeast]
Chen HC, Lee HC, Lin TY, Li WH, Chen BS: Quantitative characterization of the transcriptional regulatory network yeast cell cycle. Bioinformatics 2004, 20: 1914–1927. 10.1093/bioinformatics/bth178
Lee TI, Rinaldi NJ, Robert F, Odom DT, BarJoseph Z, Gerber GK, Hannett NM, Harbison CT, Thompson CM, Simon I, Zeitlinger J, Jennings EG, Murray HL, Gordon DB, Ren B, Wyrick JJ, Tagne JB, Volkert TL, Fraenkel E, Gifford DK, Young RA: Transcriptional Regulatory Network in Saccharomyces cerevisie. Science 2002, 298: 799–804. 10.1126/science.1075090
Identifying regulatory networks by combinatorial analysis of promoter elements[http://genetics.med.harvard.edu/~tpilpel/MotComb.html]
Casella G, Berger R: Statistical Inference. Belmont, CA: Duxbury Press; 2001.
Akaike H: A new look at the statistical model identification. IEEE Trans Autom Control 1974, AC19: 716–723. 10.1109/TAC.1974.1100705
Weisberg S: Applied linear regression. New York: John Wiley; 1985.
Acknowledgements
ACH is grateful for support from Laboratoire IGS CNRS Marseille where this work has been initiated. The authors thanks Dr. Karsten Suhre and Dr. Yves Vandenbrouck for useful discussions.
This article has been published as part of BMC Bioinformatics Volume 8, Supplement 5, 2007: Articles selected from posters presented at the Tenth Annual International Conference on Research in Computational Biology. The full contents of the supplement are available online at http://www.biomedcentral.com/14712105/8?issue=S5.
Author information
Authors and Affiliations
Corresponding author
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
ACH devised the model, performed statistical analysis and drafted the manuscript. MDQ implemented the model, collected the data and edited the manuscript. Both authors read and approved the final manuscript.
Electronic supplementary material
12859_2007_1917_MOESM1_ESM.xls
Additional File 1: Lists of strings (for the gene names) and numerical data (for their fitting parameters). The full table obtained from processing with the SDE beta sigmoid method the Spellman data [1] (XLS 2 MB)
Rights and permissions
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
About this article
Cite this article
ClimescuHaulica, A., Quirk, M.D. A stochastic differential equation model for transcriptional regulatory networks. BMC Bioinformatics 8, S4 (2007). https://doi.org/10.1186/147121058S5S4
Published:
DOI: https://doi.org/10.1186/147121058S5S4
Keywords
 mRNA Expression Level
 Regulatory Function
 Gene Expression Measurement
 Drift Term
 Transcriptional Regulatory Network