 Methodology Article
 Open Access
 Published:
Genomic data assimilation using a higher moment filtering technique for restoration of gene regulatory networks
BMC Systems Biologyvolume 9, Article number: 14 (2015)
Abstract
Background
As a result of recent advances in biotechnology, many findings related to intracellular systems have been published, e.g., transcription factor (TF) information. Although we can reproduce biological systems by incorporating such findings and describing their dynamics as mathematical equations, simulation results can be inconsistent with data from biological observations if there are inaccurate or unknown parts in the constructed system. For the completion of such systems, relationships among genes have been inferred through several computational approaches, which typically apply several abstractions, e.g., linearization, to handle the heavy computational cost in evaluating biological systems. However, since these approximations can generate false regulations, computational methods that can infer regulatory relationships based on less abstract models incorporating existing knowledge have been strongly required.
Results
We propose a new data assimilation algorithm that utilizes a simple nonlinear regulatory model and a state space representation to infer gene regulatory networks (GRNs) using timecourse observation data. For the estimation of the hidden state variables and the parameter values, we developed a novel method termed a higher moment ensemble particle filter (HMEnPF) that can retain first four moments of the conditional distributions through filtering steps. Starting from the original model, e.g., derived from the literature, the proposed algorithm can sequentially evaluate candidate models, which are generated by partially changing the current best model, to find the model that can best predict the data. For the performance evaluation, we generated six synthetic data based on two real biological networks and evaluated effectiveness of the proposed algorithm by improving the networks inferred by previous methods. We then applied timecourse observation data of rat skeletal muscle stimulated with corticosteroid. Since a corticosteroid pharmacogenomic pathway, its kinetic/dynamics and TF candidate genes have been partially elucidated, we incorporated these findings and inferred an extended pathway of rat pharmacogenomics.
Conclusions
Through the simulation study, the proposed algorithm outperformed previous methods and successfully improved the regulatory structure inferred by the previous methods. Furthermore, the proposed algorithm could extend a corticosteroid related pathway, which has been partially elucidated, with incorporating several information sources.
Background
Gene regulatory networks (GRNs) are fundamental for sustaining complex biological systems in cells. Although a comprehensive understanding of intracellular systems is still far from complete, many findings regarding intracellular systems have been published as a result of recent technological advances in biotechnology, e.g., microarray and ChipSeq. By combining these findings, we can construct biological simulation models in which the dynamics of biomolecules are described by mathematical equations, e.g., the MichaelisMenten model [1] and Ssystem [2]. However, simulation results may not match results from biological observations due to inaccurate or missing information about intracellular systems.
In oder to infer unknown parts of biological systems, there exist roughly two major approaches, i.e., simulation modelbased and statistical approaches. In constructing biological simulation models, regulatory relationships among biomolecules are collected from the literature. To represent the regulatory systems, mathematical equations, often differential equations [14], are given to describe the dynamic behavior of the involved biomolecules. The parameter values of these simulation models have been estimated to be consistent with the data by some computational methodologies. Several methods have been proposed to infer regulatory structures [5,6], to reproduce the dynamic behavior of biological systems recorded in the literature [710], and to improve published pathways so that they are consistent with the data [11,12]. However, differential equationbased approaches are computationally intensive when updating parameter values and simulation results simultaneously. Therefore, they cannot be applied to more than several genes when much of the regulatory structure is unknown.
A statistical approach using more abstracted models, e.g., Bayesian networks [1316] and the state space model [1721], have been successfully applied to infer the structure of transcriptional regulation using data from biological observations. Whereas purely datadriven methods need to explore a large model space, some studies have further incorporated other information, e.g., literaturerecorded pathways and TFs information [2226]. In contrast, these approximations can generate false regulations; there is a tradeoff relationship between accuracy and computational ease. To overcome the problem, methods to improve and deconvolve networks, which are inferred by some computational approaches, utilizing less abstract models to better predict the data have been also proposed recently [2729]. In following the direction, we should apply models that can maximally emulate the nonlinear dynamics of gene regulatory networks and establish a method for estimating the parameter values that maximize the ability to predict the data.
For this purpose, we proposed a novel data assimilation algorithm utilizing a simple nonlinear model, termed the combinatorial transcription model [5,30], and a state space representation [31,32], to infer GRNs by restoring networks that is inferred by some GRNs inference methods or derived from the literature. Since the nonlinearity results in generating nonGaussian conditional distributions of the hidden state variables, we applied the unscented Kalman filter (UKF) [3335] that can efficiently calculate the approximated conditional distributions as Gaussian distributions [36]. However, UKF cannot satisfy the requirements for estimating accurate parameter values of the model; thus, the first four moments of the conditional distributions of the hidden states should be retained. To address this problem, we developed a novel method, termed a highermoment ensemble particle filter (HMEnPF), which can retain the first two moments and the third and fourth central moments throughout the prediction, filtering, and smoothing steps. Starting from an original network, which is derived from the literature or some GRNs inference methods, the proposed algorithm using HMEnPF improves the network based on the nonlinear state space model. Furthermore, the combinatorial transcription model was extended so that the model can handle additional biomolecules such as drugs.
To show the effectiveness of the proposed algorithm, we first prepared synthetic timecourse data and compared the proposed algorithm to GeneNet [37,38] based on an empirical graphical Gaussian model (GGM), G1DBN [39] based on dynamic Bayesian networks using first order conditional dependencies, and the previous algorithm using UKF only [36]. For this comparison, six synthetic data with 30 timepoints were generated based on a WNT5A [40] and a yeastcellcycle network [41]. As an application example, we used the timecourse microarray data after stimulating rat skeletal muscle with corticosteroid, which were downloaded from the GEO database (GSE490). For this experiment, we also utilized corticosteroid pharmacogenomics [42,43], a previously defined regulatory structure for rat skeletal muscle [44], TF information from ITFP (Integrated Transcription Factor Platform) [45] and the original network inferred by G1DBN. Consequently, we proposed candidate pathways for an extension of corticosteroidrelated pathways.
Methods
State space representation of combinatorial transcription model
Let x _{ i }(t) be the abundance of the ith (i=1,…,p) gene as a function of time t. As a gene regulatory system, we assume that x _{ i }(t) is controlled by its synthesis and degradation processes, and that the quantity of synthesis is regulated by the other genes as described by
where f _{ i } is a function of the regulatory effect on the ith gene by other genes, x(t)=(x _{1}(t),…,x _{ p }(t))^{′}, $\boldsymbol {\theta }_{f_{i}}$ is a vector of tuning parameters for f _{ i }, u _{ i } is a synthesis coefficient, d _{ i } is a degradation coefficient and v(t) is a random system noise. Here, (·)^{′} stands for transposition. f _{ i } is often considered as a hill function, such as the MichaelisMenten model [1].
Since the estimation of parameter values maximizing prediction ability is a computationally heavy task when using differential equations, difference equations have been typically utilized to analyze biological systems [4,5,17,18,20,21,46]. The impact of such substitution was discussed previously [3,4]. In this study, we handle a simple nonlinear difference equation based on the combinatorial transcription model [5,30,36] described by
where x _{ i,t } is the amount of the ith gene at time t, a _{ i,j } is an individual effect by the jth gene on the ith gene, b _{ i,(j,k)} is a combinatorial effect of the jth and kth genes on the ith gene and $\mathcal A_{i}$ is an active set of genes regulating the ith gene. Since this model is a simple extension of a linear model to express a combinatorial effect by two different genes, b _{ j,j } is not considered.
Under the framework of data assimilation, in order to combine the simulation results with the observed experimental data, we apply a state space representation of Eq. (2) given by
where x _{ t }=(x _{1,t },…,x _{ p,t })^{′}, A=(a _{1},…,a _{ p })^{′}∈R ^{p×p} is a linear effect matrix, a _{ i }=(a _{ i,1},…,a _{ i,p })^{′} (i=1,…,p), $B=(\boldsymbol {b}_{1}, \ldots, \boldsymbol {b}_{p})'\in R^{p \times p^{2}}$ is a combinatorial effect matrix, b _{ i }=(b _{ i,(1,1)},…,b _{ i,(1,p)},b _{ i,(2,1)},…,b _{ i,(p,p)})^{′} (i=1,…,p), vec is a transformation function $(R^{p \times p} \rightarrow R^{p^{2}})$ , u=(u _{1},…,u _{ p })^{′}, and v _{ t }∼N(0,Q) and w _{ t }∼N(0,R) are system and observational noises with diagonal covariance matrices, respectively. We define an entire set of time points $\mathcal T=\{1, \ldots, T\}$ and the observed time set $\mathcal {T}_{\textit {obs}}$ ( $\mathcal {T}_{\textit {obs}}\subset \mathcal {T}$ ), and consider $\mathcal {T}_{\textit {obs}}=\mathcal {T}$ in the following for simplicity. Note that A and B should be sparse matrices, and we also consider an active set of elements ${\mathcal B_{i}}(i=1, \ldots, p)$ , which are sets of nonzero columns in the ith row of B.
Incorporation of biomolecules affecting biological systems
Although the regulatory system of Eqs. (3) and (4) can only represent dynamic regulation among genes, other biomolecules, such as drugs, can affect the regulatory system. To address these cases, we further consider a term representing the concentration of other biomolecules as represented by
where d _{ t } is an Mdimensional vector containing the concentration of the biomolecules at the tth time point, G=(g _{1},…,g _{ p })^{′} is an p×M matrix and g _{ i }=(g _{ i,1},…,g _{ i,M })^{′} (i=1,…,p) is an Mdimensional vector representing their regulatory effects on the ith gene. As with ${\mathcal A}_{i}$ and ${\mathcal B}_{i}$ , we consider an active set of elements ${\mathcal G}_{i}$ for the ith row of the drug effect G. A conceptual view of Eq. (5) is illustrated in Figure 1. In using Eqs. (4) and (5), we try to infer the regulatory structure among genes and estimate the values of θ={A, B, G, u, Q, R, μ _{0}}.
A highermoment ensemble particle filter
Let Y _{ t } be {y _{1},…,y _{ t }}. In estimating the parameter values and calculating the likelihood of Eqs. (4) and (5), conditional probability densities p(x _{ t }Y _{ t−1}), p(x _{ t }Y _{ t }) and p(x _{ t }Y _{ T }) can be nonGaussian forms. Thus, since these probability densities can be analytically intractable, we applied a type of Monte Carlo approach termed ensemble approximation. In this approach, for example, a probability density p(x _{ t }) is approximated by
where $\boldsymbol {x}_{t}^{(n)}$ is the nth sample from p(x _{ t }), N is the number of samples and δ is a Dirac delta function. A sample $\boldsymbol {x}_{t}^{(n)}$ and a set of samples $\{\boldsymbol {x}_{t}^{(n)}\}$ are called particle and ensemble, respectively. Previously, many types of ensemble approximation methods have been developed to obtain conditional distributions of the hidden state variables in nonlinear state space models, e.g., the ensemble Kalman filter (EnKF) [47] and the particle filter (PF) [48,49]. Here, our requirements for this study are the followings; (i) the number of particles is not reduced through filtering steps since p and the dimension of θ can be too high for the resampling procedure and (ii) the third and fourth moments of probability densities of the hidden states should be kept to optimize θ as explained in the next subsection. To satisfy these requirements, we extended a method termed the ensemble particle filter (EnPF) [50,51], which can retain the first two moments through filtering steps, and developed a novel method termed a highermoment ensemble particle filter (HMEnPF) that can additionally retain third and fourth central moments without reducing particles. The procedure of the proposed method is explained below.
Prediction step
In this step, we attempt to calculate p(x _{ t+1}Y _{ t }) after obtaining p(x _{ t }Y _{ t }) (t=0,…,T−1). Let $\boldsymbol {x}_{tt}^{(n)}$ be a sample from a conditional probability density p(x _{ t }Y _{ t }). Initially, generate particles $\boldsymbol {x}_{00}^{(n)} \sim p(\boldsymbol {x}_{0})$ for n=1,…,N. Then, for t=0,…,T−1,

1.
Generate particles $\boldsymbol {v}_{t}^{(n)}\sim N(0, Q)$ for n=1,…,N.

2.
Calculate $\boldsymbol {x}_{t+1t}^{(n)}$ by applying $\boldsymbol {x}_{tt}^{(n)}$ and $\boldsymbol {v}_{t}^{(n)}$ to Eq. (5) for n=1,…,N.
Filtering step
In this step, we attempt to calculate p(x _{ t+1}Y _{ t+1}) after obtaining p(x _{ t+1}Y _{ t }) (t=0,…,T−1). This step consists of the following three substeps termed “Particle Filter Step”, “Ensemble Kalman Filter Step” and “Merging Step”.
At the tth ( $t \in {\mathcal T}_{\textit {obs}}$ ) time step,

1.
“Particle Filter Step” is firstly executed to obtain $\left \{\hat {\boldsymbol {x}}_{tt}^{(n)}\right \}$ that is according to the theoretically exact conditional probability density p(x _{ t }y _{ t }) as follows.

(a)
Resample $\hat {\boldsymbol {x}}_{tt}^{(n)}$ according to
$$\begin{array}{*{20}l} p(\boldsymbol{x}_{t}\boldsymbol{y}_{t}) &= \frac{1}{\sum_{\dot{n}=1}^{N}p\left(\boldsymbol{y}_{t}\boldsymbol{x}_{tt1}^{(\dot{n})}\right)}\\ &~~~~~ \times \sum_{n=1}^{N}p\left(\boldsymbol{y}_{t}\boldsymbol{x}_{tt1}^{(n)}\right)\delta\left(\boldsymbol{x}_{t}\boldsymbol{x}_{tt1}^{(n)}\right). \end{array} $$((7)) 
(b)
Calculate the first and second moments $\mu _{tt} = E\left [\left \{\hat {\boldsymbol {x}}_{tt}^{(n)}\right \}\right ]$ and $V_{tt}=Var\left [\left \{\hat {\boldsymbol {x}}_{tt}^{(n)}\right \}\right ]$ , respectively.

(c)
Standardize $\hat {\boldsymbol {x}}_{tt}^{(n)}$ as
$$\begin{array}{*{20}l} \hat{\boldsymbol{z}}_{tt}^{(n)} = V_{tt}^{\frac{1}{2}} \cdot \left(\hat{\boldsymbol{x}}_{tt}^{(n)}  \boldsymbol{\mu}_{tt}\right). \end{array} $$((8)) 
(d)
Calculate the third and fourth central moments $\hat {\boldsymbol {m}}_{tt}^{(3)}=E\left [\left \{\hat {\boldsymbol {z}}_{tt}^{(n)}\right \}^{3}\right ]$ and $\hat {\boldsymbol {m}}_{tt}^{(4)} =E\left [\left \{\hat {\boldsymbol {z}}_{tt}^{(n)}\right \}^{4}\right ]$ , respectively.

(a)

2.
“Ensemble Kalman Filter Step” is secondly executed to obtain $\left \{\tilde {\boldsymbol {x}}_{tt}^{(n)}\right \}$ that is calculated under the Gaussian assumption with regard to p(x _{ t }y _{ t }) as follows.

(a)
Generate particles $\boldsymbol {w}_{t}^{(n)}\sim N(0, R)$ for n=1,…,N.

(b)
Calculate Kalman gain
$$\begin{array}{*{20}l} K_{t} = V_{tt1}\left(V_{tt1} + R_{t}\right)^{1}, \end{array} $$((9))where $V_{tt1} = Var\left [\left \{\boldsymbol {x}_{tt1}^{(n)}\right \}\right ]$ and $R_{t} = Var\left [\left \{\boldsymbol {w}_{t}^{(n)}\right \}\right ]$ .

(c)
Calculate $\tilde {\boldsymbol {x}}_{tt}^{(n)}$ as
$$\begin{array}{*{20}l} \tilde{\boldsymbol{x}}_{tt}^{(n)} = \boldsymbol{x}_{tt1}^{(n)} + K_{t} \left(\boldsymbol{y}_{t}  \boldsymbol{x}_{tt1}^{(n)} + \boldsymbol{w}_{t}^{(n)}\right). \end{array} $$((10)) 
(d)
Calculate the first and second moments $\tilde {\mu }_{tt} = E\left [\left \{\tilde {\boldsymbol {x}}_{tt}^{(n)}\right \}\right ]$ and $\tilde {V}_{tt}=Var\left [\left \{\tilde {\boldsymbol {x}}_{tt}^{(n)}\right \}\right ]$ , respectively.

(e)
Standardize $\tilde {\boldsymbol {x}}_{tt}^{(n)}$ as
$$\begin{array}{*{20}l} \tilde{\boldsymbol{z}}_{tt}^{(n)} = \tilde{V}_{tt}^{\frac{1}{2}} \cdot \left(\tilde{\boldsymbol{x}}_{tt}^{(n)} \tilde{\boldsymbol{\mu}}_{tt}\right). \end{array} $$((11)) 
(f)
Calculate the third and fourth central moments $\tilde {\boldsymbol {m}}_{tt}^{(3)}=E\left [\left \{\tilde {\boldsymbol {z}}_{tt}^{(n)}\right \}^{3}\right ]$ and $\tilde {\boldsymbol {m}}_{tt}^{(4)} =E\left [\left \{\tilde {\boldsymbol {z}}_{tt}^{(n)}\right \}^{4}\right ]$ , respectively.

(a)

3.
“Merging Step” is finally executed to obtain $\left \{\boldsymbol {x}_{tt}^{(n)}\right \}$ of which the first, second, third central and fourth central moments match to those of $\left \{\hat {\boldsymbol {x}}_{tt}^{(n)}\right \}$ . Here, we consider a standardization function S(γ,α,β) that transforms a normal random vector γ into a normalized random vector x whose the third and fourth central moments are α and β, respectively. From a previous study [52], we have S(γ,α,β) and S _{ inv }(x,α,β) that transforms x to γ as explained in the Additional file 1. Then, we obtained $\boldsymbol {x}_{tt}^{(n)}$ as
$$\begin{array}{*{20}l} \boldsymbol{x}_{tt}^{(n)} &= \hat{V}_{tt}^{\frac{1}{2}}S\left(\boldsymbol{z}_{t}^{(n)}, \hat{\boldsymbol{m}}_{tt}^{(3)}, \hat{\boldsymbol{m}}_{tt}^{(4)}\right) + \hat{\boldsymbol{\mu}}_{tt}, \end{array} $$((12))$$\begin{array}{*{20}l} \boldsymbol{z}_{t}^{(n)} &= S_{inv}\left(\tilde{\boldsymbol{z}}_{tt}^{(n)}, \tilde{\boldsymbol{m}}_{tt}^{(3)}, \tilde{\boldsymbol{m}}_{tt}^{(4)}\right). \end{array} $$((13))
Smoothing step
The smoothing step used for calculating x _{ ts } (s>t) was essentially the same as the filtering step. The details of the smoothing step can be found in the Additional file 2.
Parameter estimation using EMalgorithm
Let X _{ T }={x _{0}, …, x _{ T }} be the set of state variables. The loglikelihood of observation data is given by
where p(x _{0}) is a probability density of Ndimensional Gaussian distributions N(μ _{0},Σ _{0}), p(x _{ t }x _{ t−1}) and p(y _{ t }x _{ t }) can be probability densities of Ndimensional nonGaussian distributions obtained by ensemble approximation.
In this study, we estimate the values of θ by maximizing Eq. (14) using the EMalgorithm. Thus, the conditional expectation of the joint loglikelihood of the complete data (X _{ T },Y _{ T }) at the lth iteration
is iteratively maximized with respect to θ until the convergence is attained. More details are included in the Additional file 3.
Network restoration algorithm
We consider an algorithm to explore the best model by sequentially evaluating candidate models generated from the current best model ${\mathcal M}_{\textit {current}}$ by partially modifying the regulation. Briefly, given the original model ${\mathcal M}_{\textit {original}}$ , we attempt to sequentially create and evaluate candidates that are generated by adding, deleting and replacing regulatory components of ${\mathcal M}_{\textit {current}}$ until the best model is no longer updated. A conceptual view is illustrated in Figure 2.
Due to the heavy computational cost to evaluate the model by HMEnPF, we proposed a novel algorithm for reconstructing GRNs with combining UKF and HMEnPF as described in Algorithms 1, 2, 3, 4 and 5 and illustrated in Figure 3. Compared to EnKF (the computational task of EnKF is included in HMEnPF), the computational cost for UKF in prediction, filtering, and smoothing steps are roughly $\frac {2p+1}{N}$ , $\frac {1}{N}$ and $\frac {2}{N\cdot \mathcal {T}_{\textit {obs}}}$ , respectively. The theoretical details of UKF for the combinatorial transcription model were discussed previously [36]. Briefly, the proposed algorithm first calculate e _{ a }, e _{ b } and e _{ g } explained in the Additional file 4 for all candidate models, next evaluate the top r _{1} candidates for each row by UKF and then evaluate the r _{2} top candidates by HMEnPF. Note that, when the systems include G, regulations by the drugs are inferred in the same way as A in Algorithms 1, 2, 3, 4 and 5. In Results and discussion section, we set {r _{1},r _{2},a d d _{ max },d e l _{ max }}={5,5,+∞,+∞}.
Results and discussion
Comparison using synthetic data
To show the effectiveness of the proposed algorithm, we prepared synthetic timecourse gene expression data based on the synthetic networks, WNT5A [40] and a yeast cellcycle [41], as illustrated in Figures 4 and 5, respectively. For each network and three different system noises, we generated five timecourses ( ${\mathcal T}=\{1,2,\ldots,30\}$ ) by using Eqs. (3) and (4); thus, six sets of five timecourses were prepared. The values of the parameters A, B and u in Eq. (3) were determined between 1 and 1, the system noise v _{ t } and observational noise w _{ t } were generated according to Gaussian distributions with a mean 0 and three variances 0.01, 0.05 and 0.1, and that with a mean 0 and a variance 0.3, respectively. For the original networks to be improved by the proposed algorithm, we utilized GeneNet [37,38] based on an empirical graphical Gaussian model (GGM) and G1DBN [39] based on dynamic Bayesian networks using first order conditional dependencies. After the restoration, the original and improved networks were evaluated by true positive (TP), false positive (FP), true negative (TN), false negative (FN), precision rate $\left (PR=\frac {\text {TP}}{\text {TP}+\text {FP}}\right)$ , recall rate $\left (RR={\frac {\text {TP}}{\text {{TP}+\text {FN}}}}\right)$ and Fmeasure $\left (=\frac {\text {2PR}\cdot \text {RR}}{\text {PR}+\text {RR}}\right)$ . Note that, since GeneNet infers undirected regulations among genes, we compared its results to the undirected true networks. In addition, since a directed network is required for the original network, we transformed the undirected network inferred by GeneNet as follows; (i) a true directed regulation was set when an inferred undirected regulation was correct, and (ii) a false directed regulation of which direction was randomly selected was set when an inferred undirected regulation was incorrect. Here, to clarify the significance of HMEnPF, we also showed the results of the previous algorithm using UKF only [36]. These results are summarized in Tables 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 and 12.
The results indicate that the proposed algorithm using HMEnPF and only UKF could outperform G1DBN and GeneNet, and the proposed algorithm showed better performance than that of using UKF only. This concludes that retaining higher moment information can improve the accuracy of approximation and estimate correct parameter values. Additionally, we recognized that the performance of the proposed algorithm strongly depends on the accuracy of the original network. Thus, to obtain better results, we should carefully construct original networks or select inference methods for creating the original network. Note that the Jar file of the proposed algorithm is available at: http://sunflower.kuicr.kyotou.ac.jp/~thasegw/, and the synthetic data, the parameter values and the original networks are in the Additional file 5.
Inference using real data
As an application example, we analyzed microarray timecourse gene expression data from rat skeletal muscle [42,43]. The microarray data were downloaded from the GEO database (GSE490). The timecourse gene expression data was measured at 0, 0.25, 0.5, 0.75, 1, 2, 4, 5, 5.5, 7, 8, 12, 18, 30, 48, and 72 [h] (16 time points) after stimulation of corticosteroid, but we removed data at 48 and 72 [h] (steady state profiles) for computational efficiency. The data at time 0 represent controls (untreated). There were two, three, or four replicated observations for each time point.
Because corticosteroid pharmacokinetics/dynamics in skeletal muscle have been modeled based on differential equations [43], the timedependent concentration of corticosteroid in nucleus in rat skeletal muscle can be obtained for d _{ t } as explained in the Additional file 6. Furthermore, corticosteroid catabolic/anabolic processes in rat skeletal muscle have been partially established [44]; thus, we handled gene (i) TFs, Trim63, Akt1, Akt2, Mstn, Mtor, Irs1, and (ii) nonTFs, Akt3, Anxa3, Bcat2, Bnip3, Foxo1, Igf1, Igf1r, Pik3c3, Pik3cb, Pik3cd, Pik3c2g, Rheb, Slc2a4 with their regulatory relationships. Additionally, we handled genes (iii) TFs, Cebpb, Cebpd, Gpam, Srebf1 and (iv) nonTFs, Rxrg, Scarb1, Scd, Gpd2, Mapk6, Ace, Ptpn1, Ptprf, Edn1, Agtr1a, Ppard, Hmgcs2, Serpine1, Il6r, Mapk14, Ucp3 and Pdk4 that are suggested as corticosteroid related genes [42]. Note that the microarray (GSE490) does not include three genes in the original pathway [44], Redd1, Bcaa and Klf15. In summary, we handled the concentration of corticosteroid in nucleus, these 40 genes (shown in Table 13) and an original network that was inferred by G1DBN with regulatory relationships among (i) and (ii). Note that TFs information was derived from ITFP (Integrated Transcription Factor Platform) [45].
Consequently, we obtained the improved network as illustrated in Figure 6. A purple circle, blues circles, and green circles represent corticosterid, TF candidates and nonTF candidates, respectively. In the center of this figure, there exist corticosteroid regulations to several TF and nonTF genes and regulatory effects transmit to down stream genes of TF candidates genes. In addition, there exist some interesting findings. At first, genes included in ‘response to insulin stimulus (GO:0032868)’ and ‘insulin receptor binding (GO:0005158)’, ‘Igf1’, ‘Akt1’, ‘Akt2’, ‘Srebf1’, ‘Ptprf’, ‘Mtor’ and ‘Ptpn1’, construct a regulatory component in the bottom right of this figure. Including ‘Cebpd’ and ‘Cebpb’, which are assumed to be candidate genes for insulinrelated transcription factors and selected as hub genes, functional relationships between corticosteroid and insulinrelated functions were reported [53]. On the other hand, ‘Irs1’, ‘Bcat2’, ‘Edn1’, ‘Ucp3’, ‘Pdk4’, ‘Mstn’, ‘Foxo1’ and ‘Rxrg’ that are involved in ‘positive regulation of metabolic process (GO:0009893)’ and ‘fatty acid metabolic process (GO:0006631)’ build the other regulatory process. Since some combinatorial regulations were inferred, it is conceivable that higher moment approximation can affect the estimation results beyond linear models.
Conclusions
In this paper, we developed a novel approach to restore original GRNs to be consistent with timecourse mRNA expression data based on the combinatorial transcription model. Since we applied a state space representation with the nonlinear system equation in the context of data assimilation, the conditional distributions of the hidden variables can be nonGaussian distributions. In contrast to the previous approaches using particle filter, Gaussian approximation and regressionbased solutions, our proposed approach, HMEnPF, can retain the first, second, third central and fourth central moments through filtering steps to estimate near optimal parameter values by the EMalgorithm.
According to the comparison results using six synthetic data based on the real biological pathways, the proposed algorithm successfully explored better models than the previous methods, G1DBN and GeneNet, considering linear relevance. Moreover, the proposed algorithm using HMEnPF outperformed that of using UKF. This concludes that HMEnPF retaining parts of higher moment information can improve the accuracy of the estimation of the parameter values beyond unscented approximation (that cannot retain any moment through filtering steps based on Gaussian approximation). Through the experimental results, we also observed that the performance of the restoration algorithm strongly depends on the original network, which was prepared by literature information or some GRNs inference methods. Thus, one of significant points is to select methods to infer the original network. On the other hand, the proposed method has some limitations. For example, we require timecourse data in which the number of time points should be more than 10 or so. Moreover, due to its heavy computational costs, the calculation for more than 20−30 genes without TF information can be infeasible.
As an application example, we prepared corticosteroid pharmacogenomic pathways in rat skeletal muscle that have been investigated and established a part of regulatory relationships and related genes. Additionally, the intracellular concentration of corticosteroid that directly/indirectly affects gene expression can be obtained by the previously developed differential equations and TF information for rat genes can also be utilized. In summary, we inferred the regulatory relationships among 40 genes and corticosteroid with fixing the established pathways and restricting that only TF candidates can regulate other genes. G1DBN was employed to construct the original model for the proposed algorithm. Consequently, several combinatorial regulations and regulations by corticosteroid were inferred by extending the original network. Since previous linear models may not be able to infer these regulations, the proposed algorithm can be valuable to restore inferred and literaturebased networks to be consistent with the data.
References
 1
Savageau MA. Biochemical systems analysis: II. The steadystate solutions for an npool system using a powerlaw approximation. J Theor Biol. 1969; 25(3):370–9.
 2
Savageau MA, Voit EO. Recasting nonlinear differential equations as ssystems: a canonical nonlinear form. Math Biosci. 1987; 87(1):83–115.
 3
Elowitz MB, Leibler S. A synthetic oscillatory network of transcriptional regulators. Nature. 2000; 403(6767):335–8.
 4
de Jong H. Modeling and simulation of genetic regulatory systems: A literature review. J Comput Biol. 2002; 9(1):67–103.
 5
Opper M, Sanguinetti G. Learning combinatorial transcriptional dynamics from gene expression data. Bioinformatics. 2010; 26(13):1623–9.
 6
Henderson J, Michailidis G. Network reconstruction using nonparametric additive ode models. PLoS ONE. 2014; 9(4):94003.
 7
Koh CHH, Nagasaki M, Saito A, Wong L, Miyano S. DA 1.0: parameter estimation of biological pathways using data assimilation approach. Bioinformatics. 2010; 26(14):1794–6.
 8
Matsuno H, Nagasaki M, Miyano S. Hybrid petri net based modeling for biological pathway simulation. Nat Comput. 2011; 10:1099–120.
 9
Ramsay JO, Hooker G, Campbell D, Cao J. Parameter estimation for differential equations: a generalized smoothing approach. J R Stat Soc: Ser B (Stat Methodology). 2007; 69(5):741–96.
 10
Quach M, Brunel N, d’AlcheBuc F. Estimating parameters and hidden variables in nonlinear statespace models based on odes for biological networks inference. Bioinformatics. 2007; 23(23):3209–16.
 11
Hasegawa T, Yamaguchi R, Nagasaki M, Imoto S, Miyano S. Comprehensive pharmacogenomic pathway screening by data assimilation. In: Proceedings of the 7th International Conference on Bioinformatics Research and Applications. ISBRA’11. Berlin, Heidelberg: Springer: 2011. p. 160–171.
 12
Hasegawa T, Nagasaki M, Yamaguchi R, Imoto S, Miyano S. An efficient method of exploring simulation models by assimilating literature and biological observational data. Biosystems. 2014; 121(0):54–66.
 13
Friedman J, Hastie T, Tibshirani R. Sparse inverse covariance estimation with the graphical lasso. Biostatistics. 2007; 9(3):432–41.
 14
Kim S, Imoto S, Miyano S. Dynamic bayesian network and nonparametric regression for nonlinear modeling of gene networks from time series gene expression data. Biosystems. 2004; 75(13):57–65.
 15
Young W, Raftery A, Yeung K. Fast Bayesian inference for gene regulatory networks using ScanBMA. BMC Syst Biol. 2014; 8(1):47.
 16
Zacher B, Abnaof K, Gade S, Younesi E, Tresch A, Fröhlich H. Joint Bayesian inference of conditionspecific miRNA and transcription factor activities from combined gene and microRNA expression data. Bioinformatics. 2012; 28(13):1714–20.
 17
Barenco M, Tomescu D, Brewer D, Callard R, Stark J, Hubank M. Ranked prediction of p53 targets using hidden variable dynamic modeling. Genome Biol. 2006; 7(3):25.
 18
Beal MJ, Falciani F, Ghahramani Z, Rangel C, Wild DL. A bayesian approach to reconstructing genetic regulatory networks with hidden factors. Bioinformatics. 2005; 21:349–56.
 19
Hasegawa T, Yamaguchi R, Nagasaki M, Miyano S, Imoto S. Inference of gene regulatory networks incorporating multisource biological knowledge via a state space model with l1 regularization. PLoS ONE. 2014; 9(8):105942.
 20
Hirose O, Yoshida R, Imoto S, Yamaguchi R, Higuchi T, CharnockJones DS, et al. Statistical inference of transcriptional modulebased gene networks from time course gene expression profiles by using state space models. Bioinformatics. 2008; 24:932–42.
 21
Rangel C, Angus J, Ghahramani Z, Lioumi M, Sotheran E, Gaiba A, et al. Modeling tcell activation using gene expression profiling and statespace models. Bioinformatics. 2004; 20:1361–72.
 22
Sabatti C, James GM. Bayesian sparse hidden components analysis for transcription regulation networks. Bioinformatics. 2006; 22(6):739–46.
 23
Asif HMS, Sanguinetti G. Largescale learning of combinatorial transcriptional dynamics from gene expression. Bioinformatics. 2011; 27(9):1277–83.
 24
Eduati F, De Las Rivas J, Di Camillo B, Toffolo G, SaezRodriguez J. Integrating literatureconstrained and datadriven inference of signalling networks. Bioinformatics. 2012; 28(18):2311–7.
 25
do Rego TG, Roider HG, de Carvalho FAT, Costa IG. Inferring epigenetic and transcriptional regulation during blood cell development with a mixture of sparse linear models. Bioinformatics. 2012; 28(18):2297–303.
 26
Tian Y, Zhang B, Hoffman E, Clarke R, Zhang Z, Shih IM, et al. Knowledgefused differential dependency network models for detecting significant rewiring in biological networks. BMC Syst. Biol. 2014; 8(1):87.
 27
Barzel B, Barabási ALL. Network link prediction by global silencing of indirect correlations. Nat Biotechnol. 2013; 31(8):720–5.
 28
Feizi S, Marbach D, Medard M, Kellis M. Network deconvolution as a general method to distinguish direct dependencies in networks. Nat Biotechnol. 2013; 31(8):726–33.
 29
Nakajima N, Tamura T, Yamanishi Y, Horimoto K, Akutsu T. Network completion using dynamic programming and leastsquares fitting. Sci World J. 2012; 2012:1–8.
 30
Wang W, Cherry JM, Nochomovitz Y, Jolly E, Botstein D, Li H. Inference of combinatorial regulation in yeast transcriptional networks: A case study of sporulation. Proc Nat Acad Sci USA. 2005; 102(6):1998–2003.
 31
Kalman RE. A New Approach to Linear Filtering and Prediction Problems. Trans ASME  J Basic Eng. 1960; 82(Series D):35–45.
 32
Shumway RH, Stoffer DS. An approach to time series smoothing and forecasting using the em algorithm. J Time Ser Anal. 1982; 3(4):253–64.
 33
Julier SJ, Uhlmann JK. A new extension of the kalman filter to nonlinear systems. In: Proc. of AeroSense: The 11th Int. Symp. on Aerospace/Defense Sensing, Simulations and Controls: 1997. p. 182–193.
 34
Julier SJ, Uhlmann JK. Unscented filtering and nonlinear estimation. Proc IEEE. 2004; 92(3):401–22.
 35
Chow SM, Ferrer E, Nesselroade JR. An unscented kalman filter approach to the estimation of nonlinear dynamical systems models. Multivariate Behavioral Res. 2007; 42(2):283–321.
 36
Hasegawa T, Mori T, Yamaguchi R, Imoto S, Miyano S, Akutsu T. An efficient data assimilation schema for restoration and extension of gene regulatory networks using timecourse observation data. J Comput Biol. 2014; 21(11):785–98.
 37
Schäfer J, Strimmer K. An empirical bayes approach to inferring largescale gene association networks. Bioinformatics. 2005; 21(6):754–64.
 38
OpgenRhein R, Strimmer K. From correlation to causation networks: a simple approximate learning algorithm and its application to highdimensional plant gene expression data. BMC Syst Biol. 2007; 1(1):37.
 39
Lébre S. Inferring dynamic genetic networks with low order independencies. Stat App Genet Mol Biol. 2009; 8(1):1–38.
 40
Kim S, Li H, Dougherty ER, Cao N, Chen Y, Bittner M, et al. Can markov chain models mimic biological regulation. J Biol Syst. 2002; 10(4):337–28093357.
 41
Kanehisa M, Goto S, Sato Y, Furumichi M, Tanabe M. Kegg for integration and interpretation of largescale molecular data sets. Nucleic Acids Res. 2012; 40(D1):109–14.
 42
Almon RR, DuBois DC, Jin JY, Jusko WJ. Temporal profiling of the transcriptional basis for the development of corticosteroidinduced insulin resistance in rat muscle. J Endocrinol. 2005; 184(1):219–32.
 43
Yao Z, Hoffman EP, Ghimbovschi S, DuBois DC, Almon RR, Jusko WJ. Mathematical modeling of corticosteroid pharmacogenomics in rat muscle following acute and chronic methylprednisolone dosing. Mol Pharm. 2008; 5(2):328–39.
 44
Shimizu N, Yoshikawa N, Ito N, Maruyama T, Suzuki Y, Takeda SI, et al. Crosstalk between Glucocorticoid Receptor and Nutritional Sensor mTOR in Skeletal Muscle. Cell Metab. 2011; 13(2):170–82.
 45
Zheng G, Tu K, Yang Q, Xiong Y, Wei C, Xie L, et al. Itfp: an integrated platform of mammalian transcription factors. Bioinformatics. 2008; 24(20):2416–7.
 46
Greenfield A, Hafemeister C, Bonneau R. Robust datadriven incorporation of prior knowledge into the inference of dynamic regulatory networks. Bioinformatics. 2013; 29(8):1060–7.
 47
Evensen G. Sequential data assimilation with a nonlinear quasigeostrophic model using monte carlo methods to forecast error statistics. J Geophys Res. 1994; 99:10143–62.
 48
Gordon NJ, Salmond DJ, Smith AFM. Novel approach to nonlinear/nongaussian bayesian state estimation. IEEE Proc F, Radar Signal Process. 1993; 140(2):107–13.
 49
Kitagawa G. Monte Carlo Filter and Smoother for NonGaussian Nonlinear State Space Models. J Comput Graphical Stat. 1996; 5(1):1–25.
 50
Anderson LJ, Anderson LS. A monte carlo implementation of the nonlinear filtering problem to produce ensemble assimilations and forecasts. Monthly Weather Rev. 1999; 127(12):2741–58.
 51
Pham DT. Stochastic methods for sequential data assimilation in strongly nonlinear systems. Monthly Weather Rev. 2001; 129(5):1194–207.
 52
Zhao Y, Lu Z. Fourthmoment standardization for structural reliability assessment. J Struct Eng. 2007; 133(7):916–24.
 53
Foti D, Iuliano R, Chiefari E, Brunetti A. A nucleoprotein complex containing sp1, c/ebpb, and hmgiy controls human insulin receptor gene transcription. Mol Cell Biol. 2003; 23(8):2720–32.
Acknowledgements
The supercomputing resource was provided by Human Genome Center, the Institute of Medical Science, the University of Tokyo (http://sc.hgc.jp/shirokane.html).
This work was partly supported by GrantinAid for JSPS Fellows Number 249639.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
The work presented here was carried out in collaboration between all authors. TH conceived and designed the study, and wrote the manuscript. TM assisted in constructing biological models and preparing the manuscript. RY and TS provided statistical expertise and careful manuscript review. SM and SI provided valuable advises in developing the proposed method from the point of view of statistics and bioinformatics. TA supervised the whole project. All authors read and approved the final manuscript.
Additional files
Additional file 1
The standardization function for HMEnPF. The standardization function and its inverse function for HMEnPF are described in the file.
Additional file 2
The procedure of the smoothing step in HMEnPF. The procedure of the smoothing step in HMEnPF is described in the file.
Additional file 3
The detailed solution of the EMalgorithm for the estimation of the parameter values. The detailed solution of the E and Msteps in the EMalgorithm for HMEnPF are described in the file.
Additional file 4
The functions measuring the effectiveness of regulatory edges. The functions measuring the effectiveness of regulatory edges e _{ a }, e _{ b } and e _{ g } are described in the file.
Additional file 5
The synthetic data from WNT5A and a yeast cellcycle networks. The six synthetic data from WNT5A and a yeast cellcycle networks, their parameter values and the original networks are included in the file.
Additional file 6
The corticosteroid pharmacokinetics/dynamics in rat skeletal muscle. The differential equations of rat pharmacokinetics/dynamics are described in the file.
Rights and permissions
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Gene regulatory networks
 Time series analysis
 Systems biology
 Data assimilation
 Monte Carlo