Gene regulatory network inference by point-based Gaussian approximation filters incorporating the prior information

The extended Kalman filter (EKF) has been applied to inferring gene regulatory networks. However, it is well known that the EKF becomes less accurate when the system exhibits high nonlinearity. In addition, certain prior information about the gene regulatory network exists in practice, and no systematic approach has been developed to incorporate such prior information into the Kalman-type filter for inferring the structure of the gene regulatory network. In this paper, an inference framework based on point-based Gaussian approximation filters that can exploit the prior information is developed to solve the gene regulatory network inference problem. Different point-based Gaussian approximation filters, including the unscented Kalman filter (UKF), the third-degree cubature Kalman filter (CKF3), and the fifth-degree cubature Kalman filter (CKF5) are employed. Several types of network prior information, including the existing network structure information, sparsity assumption, and the range constraint of parameters, are considered, and the corresponding filters incorporating the prior information are developed. Experiments on a synthetic network of eight genes and the yeast protein synthesis network of five genes are carried out to demonstrate the performance of the proposed framework. The results show that the proposed methods provide more accurate inference results than existing methods, such as the EKF and the traditional UKF.


Introduction
Inferring gene regulatory network (GRN) has become one of the most important missions in system biology. Genome-wide expression data is widely used due to the development of several high-throughput experimental technologies. The gene regulatory network can be inferred from a number of gene expression samples taken over a period of time. Modeling of GRN is required before its structure can be inferred. Common dynamical modeling methods of GRN include Bayesian networks [1], Boolean networks [2], ordinary differential equations [3], state-space models [4,5], and so on. Various approaches based on different models have been used to infer the network from observed gene expression data, such as the http://bsb.eurasipjournals.com/content/2013/1 /16 Gaussian approximation filters have been recently proposed to improve the performance of the EKF, which employ various quadrature rules to compute the integrals involved in the exact Bayesian estimation. Many filters fall into this category, such as the unscented Kalman filter (UKF) [13], the Gauss-Hermite quadrature filter [14], the cubature Kalman filter (CKF) [15], and the sparse-grid quadrature filter [16]. Besides the point-based Gaussian approximation filters, the particle filter has drawn much attention recently [17]. The particle filter uses random particles with weights to represent the probability density function (pdf ) in the Bayesian estimation and provides better estimation result than the EKF. The main problem of the particle filter is that the computational complexity is high, and therefore, it is hard to use for high-dimensional problems, such as the problem considered in this paper. The EKF and the particle filter have been used for the inference of GRN [4,8,18]. In this paper, we consider the point-based Gaussian approximation filters. Our main objective is to provide a framework of incorporating network prior information into the filters. For example, some gene regulations may be known [19] from literature and the inference accuracy of GRN can be improved by incorporating the known regulations of the GRN [20]. Integration of the prior knowledge or constraints with the GRN inference algorithm has been introduced to improve the inference result. The DNA motif sequence in gene promoter regions is incorporated in [21] while modeling of transcription factor interactions is incorporated in [22]. As mentioned in [20], experimentally determined physical interactions can be obtained. In addition, the sparsity constraint is frequently used in the inference of the GRN. To the best of the authors' knowledge, the most related work in incorporating the prior information in Bayesian filters is [8]. In that work, rather than directly getting the inference results from the filter, an optimization method is used. In particular, a cost function is used in which the sparsity constraint is enforced. However, the cost function in [8] does not consider the uncertainty of the state in the filtering. That cost function in fact is not coupled well with the filtering algorithm. In addition, it did not consider other kinds of prior information. In this paper, we propose a new framework that incorporates the prior information effectively in the filtering algorithm by solving a constrained optimization problem. Efficient recursive algorithms are provided to solve the associated optimization problem.
The remainder of this paper is organized as follows. In Section 2, the modeling of gene regulatory network is introduced. The point-based Gaussian approximation filters are briefly introduced in Section 3. The proposed new filtering framework is described in Section 4. In Section 5, experimental results are provided. Finally, concluding remarks are given in Section 6.

State-space modeling of gene regulatory network
The GRN can be described by a graph in which genes are viewed as nodes and edges depict causal relations between genes. The structure of GRN reveals the mechanisms of biological cells. Analyzing the structure of GRN will pave the way for curing various diseases [23]. The learning of GRN has drawn much attention recently due to the availability the microarray data. By analyzing collected gene expression levels over a period of time, one can identify various regulatory relations between different genes. To facilitate the analysis of the GRN, modeling of GRN is required. Different models can be used, such as Bayesian networks [1], Boolean networks [2], ordinary differential equation [3], and state-space model [4,5]. The state-space model has been widely used because it incorporates noise and can make use of computationally efficient filtering algorithms [5]. Thus, we also use the state-space modeling of GRN in this paper. Under the discrete-time state-space modeling of the gene regulatory networks, the network evolution from time k to time k − 1 can be described by where x k =[x 1,k , . . . , x n,k ] T is the state vector and x i,k denotes the gene expression level of the i-th gene at time k. f is a nonlinear function that characterizes the regulatory relationship among the genes. v k is the state noise and it is assumed to follow a Gaussian distribution with mean 0 and covariance matrix Q k , i.e., v k ∼ N (0, Q k ).
Following [8], we use the following nonlinear function in the state Equation (1): with and In (2), A is the regulatory coefficient matrix with a ij denoting the regulation coefficient from gene j to gene i. Note that a positive coefficient a ij indicates that gene j activates gene i and a negative a ij indicates that gene j represses gene i. In (4), μ i is a parameter. Note that A and μ i are unknown parameters. The discrete-time nonlinear stochastic dynamic system [24] shown in Eqs. (1)- (3) have been successfully used to describe the GRN [4,8]. Equation (4) is also called Sigmoid function which is frequently used since it is consistent with the fact that all concentrations get saturated at some point in time [25]. http://bsb.eurasipjournals.com/content/2013/1/16 The Sigmoid function has been used in modeling GRN to verify various methods, such as artificial neural network [26], simulated annealing and clustering algorithm [27], extended Kalman filter [4], particle filter [8], and Genetic programming and Kalman filtering [25].
For the measurement model, we consider the following general nonlinear observation equation where h(·) is some nonlinear function, n k is the measurement noise, which is assumed to follow the Gaussian distribution with mean 0 and covariance matrix R k , i.e., n k ∼ N (0, R k ). For example, if the noise corrupted expression levels are observed, then h(x) = x.

Gaussian approximation filters
In this section, the framework of point-based Gaussian approximation filters for the state-space dynamic model is briefly reviewed. We consider the state-space model consisting of the state Equation (1) and the measurement Equation (5). We denote y k [y 1 , . . . , y k ]. The optimal Bayesian filter is composed of two steps: prediction and filtering. Specifically, given the prior pdf p(x k−1 |y k−1 ) at time k − 1, the predicted conditional pdf p(x k |y k−1 ) is given by After the measurement at time k becomes available, the filtered pdf is given by The pdf recursions in (6) and (7) are in general computationally intractable unless the system is linear and Gaussian. The Gaussian approximation filters approximate (6) and (7) by invoking Gaussian assumptions. Specifically, the first assumption is that given y k−1 , The second assumption is that (x k , y k ) are jointly Gaussian given y k−1 .
It then follows from the second assumption that given y k−1 , x k has a Gaussian distribution, i.e., x k |y k−1 ∼ N (x k|k−1 , P k|k−1 ). Using (1) and the first assumption, we have the predicted meanx k|k−1 and covariance P k|k−1 given respectively bŷ and where φ x;x, P denotes the multivariate Gaussian pdf with meanx and covariance P. Then, following the second assumption, given y k = [ y k−1 , y k ], x k is Gaussian distributed, i.e., x k |y k ∼ N (x k|k , P k|k ). Using the conditional property of the multivariate Gaussian distribution, the filtered meanx k|k and covariance P k|k are given respectively bŷ and witĥ

Point-based Gaussian approximation filters
The integrals in (8), (9), (12), (14) and (15) are Gaussian type that can be efficiently approximated by various quadrature methods. Specifically, if a set of weighted points (γ i , w i ), i = 1, . . . , N can be used to approximate the integral where P = SS T and S can be obtained by Cholesky decomposition or singular value decomposition (SVD).
Using (17), we can then approximate (8) and (9) as follows: and where ξ k−1,i is the transformed quadrature point obtained from the covariance decomposition, i.e., Similarly, we can approximate (12), (14) and (15) as follows: whereξ k,i is the transformed point obtained from the decomposition of the predicted covariance, i.e., Various numerical rules can be used to form the approximation in (16), which lead to different Gaussian approximation filters. In particular, the unscented transformation, the Gauss-Hermite quadrature rule, and the sparse-grid quadrature rules are used in the unscented Kalman filter (UKF), the Gauss-Hermite quadrature Kalman filter (GHQF), and the sparse-grid quadrature filter (SGQF), respectively.
Recently, the fifth-degree quadrature filter has been proposed and shown to be more accurate than the thirddegree quadrature filters, such as the UKF and the thirddegree cubature Kalman filter (CKF 3 ), when the system is highly nonlinear or contains large uncertainty [16]. In this paper, we consider the UKF, CKF 3 , and the fifth-degree cubature Kalman filter (CKF 5 ). Other filters such as the central difference filter [14] and divided difference filter [28] can also be used. The CKF 5 is based on Mysovskikh's method which uses fewer point than the fifth-degree quadrature filter in [16]. In the following, different numerical rules used in (16) are briefly summarized.

Unscented transform
In the unscented Kalman filter (UKF), we have N = 2n+1 where n is the dimension of x. The quadrature points and the corresponding weights are given respectively by and where κ is a tunable parameter, and e i is the i-th ndimensional unit vector in which the i-th element is 1 and other elements are 0.

Cubature rules
The left-hand side of (16) can be rewritten as Consider the integral By letting x = rs with s T s = 1 and r = √ x T x, I(h) can be rewritten in the spherical-radial coordinate system as where U n = {s ∈ R n : s = 1}, and σ (·) is the spherical surface measure or the area element on U n . Note that (31) contains two types of integrals: the radial integral ∞ 0 h r (r) r n−1 exp −r 2 dr and the spherical integral U n h s (s)dσ (s).
If the radial rule can be approximated by and the spherical integral can be approximated by A third-degree cubature rule to approximate (29) is obtained by using the third-degree spherical rule and radial rule [15]: Remark: The third-degree cubature rule is identical to the unscented transformation with κ = 0.
To construct the fifth-degree cubature rule, the Mysovskikh's method [29] and the moment matching method [16] are used to provide the fifth-degree spherical rule and radial rule, respectively. The final fifth-degree cubature rule is given by where the point s (i) 1 is given by Moreover, the set of points {s (i) 2 } is given by

Augmented state-space model for network inference
In the state-space model for gene regulatory networks described in Section 3.2, the underlying network structure is characterized by the n × n regulatory coefficient matrix A in (2) and the parameters μ = [μ 1 , . . . , μ n ] in (4). The problem of network inference then becomes to estimate A and μ. To do that, we incorporate the unknown parameters A and μ into the state vector to obtain an augmented state-space model, and then apply the point-based Gaussian approximation filters to estimate the space vector and thereby obtaining the estimates of A and μ. Specifically, we denote θ = [a 11 , a 12 , · · · , a 1n , · · · , a nn , μ 1 , · · · , μ n ] T and the augmented state vectorx k = Then, the augmented state equation can be written as In the remainder of this paper, we assume that the noisy gene expression levels are observed. Therefore, the augmented measurement equation becomes where B = [I n , O n×(n 2 +n) ], O n×(n 2 +n) denotes an n × (n 2 + n) all zeros matrix. The point-based Gaussian approximation filters can then be used to obtain the estimate of the augmented state,x k , from which the estimates of the unknown network parameters, i.e.,Â andμ can then be obtained.
Note that since the measurement Equation (41) is linear, the filtering Equations (10, 11) becomê and which are the same as the filtering updates for Kalman filters.

Incorporating prior information
In practice, some prior knowledge on the underlying GRN is typically available. In this section, we outline approaches to incorporating such prior knowledge into the point-based Gaussian approximation filters for network inference. In particular, we consider two types of prior information, namely, sparsity constraints and range constraints on the network. For networks with sparsity constraints, we incorporate an iterative thresholding procedure into the Gaussian approximation filters. http://bsb.eurasipjournals.com/content/2013/1/16 And to accommodate range constraints, we employ PDFtruncated Gaussian approximation filters.

The optimization formulations
Note that under the Gaussian assumption, the state estimationx k|k of the Kalman filter is equivalently given by the solution to the following optimization problem [30,31] x k|k = arg min with To incorporate the prior information of the GRN, (46) is modified as where J p (x) is a penalty function associated with the prior information and λ is a tunable parameter that regulates the tightness of the penalty. For example, in gene regulatory networks, each gene only interacts with a few genes [20]. To capture such a sparsity constraint, a Laplace prior distribution can be used for the connection coefficient matrix A, i.e., Therefore, in this case, J p (x) =− log p(A) = c 1 A 1 + c 2 where c 1 and c 2 are constants. And, (47) can be rewritten asJ Note that (49) can also be interpreted as the result of applying the least squares shrinkage selection operator (LASSO) to (47). The LASSO adds an L 1 -norm constraint to the GRN so that the regulatory coefficient matrix A tends to be sparse with many zero elements.
As another example, if some known regulatory relationship exists, then it should be taken into account to improve the estimation accuracy. Specifically, define an n×n indicator matrix E =[e i,j ] where e ij = 1 indicates that there is a lack of regulation from gene j to gene i. Then, similar to the use of LASSO, a penalty on a ij should incur if e ij = 1. Thus, (47) can be rewritten as Note that as in [20], here we do not force a ij = 0 corresponding to e ij = 1 but rather use an L 1 -norm penalty. The advantage of such an approach is that it allows the algorithm to pick different structures but more likely to pick the edges without penalties. 'o' denotes the entrywise product operation of two matrices.

Iterative thresholding algorithm
Solving the optimization problems in (49) and (50) is not straightforward since |a| is non-differentiable at a = 0. In the following, an efficient solver called the iterative thresholding algorithm is introduced.

PDF truncation method for range constraints
If the range constraints on the regulatory coefficients are available, the inference accuracy can be improved by enforcing the constraints in the Gaussian approximation filters.
In particular, assume that we impose the following range constraints on the state vectorx The PDF truncation method [31] can be employed to incorporate the above range constraint into the Gaussian approximation filters, by converting the updated mean x k|k and covariance P k|k to the pseudo meanx t k|k and covariance P t k|k which are then used in the next prediction and filtering steps.
We next briefly outline the PDF truncation procedure. We usex t k|k,i and P t k|k,i to denote the mean and covariance after the first i constraints have been enforced. Initially, we setx t k|k,0 =x k|k and P t k|k,0 = P k|k . Consider the following transformation where S i and D i are obtained from the Jordan canonical decomposition S i D i S T i = P t k|k,i and G i is obtained by using the Gram-Schmidt orthogonalization and it satisfies [33] Then, the upper bound e T ix ≤ d i is transformed to [33] Similarly, the lower bound e T ix ≥ c i is transformed to The constraint requires that the first element of z k,i lies betweenc i andd i . Hence, only the truncated PDF of the first element of z k,i is considered and it is given by [33] Then, the mean and variance of the first element of z k,i after imposing the i-th constraint are given respectively by Thus, the mean and covariance of the transformed state vector after imposing the i-th constraint are given respectively byz By taking the inverse transform of (63), we then get After imposing all n constraints, the final constrained state estimate and covariance at time k are given respectively byx t k|k x t k|k,n and P t k|k P t k|k,n .

Synthetic network
In this section, a synthetic network that contains eight genes is used to test the performance of the EKF, the and μ i = 2, i = 1, · · · , 8. For the filter, each coefficient inÂ is initialized from a Gaussian distribution with mean 0 and variance 0.2. Moreover, the coefficient μ i is initialized from a Gaussian distribution with mean 1.5 and variance 0.2. The system state is initialized using the first measurement.
The metric used to evaluate the inferred GRN is the true positive rate (TPR), the false positive rate (FPR), and the positive predictive value (PPV). They are given by [34] where the number of true positives (TP#) denotes the number of links correctly predicted by the inference algorithm; the number of false positives (FP#) denotes the number of incorrectly predicted links; the number of true negatives (TN#) denotes the number of correctly predicted nonlinks; and the number of false negatives (FN#) denotes the number of missed links by the inference algorithm [8].

Comparison of the EKF with point-based Gaussian approximation filters
The UKF with different parameter κ is tested. The simulation results based on 50 Monte Carlo runs are shown in Table 1. It can be seen that UKFs with κ = 0, 2, 5 have slight better performance than UKFs with κ = −5, −2.
One possible reason is that the weights of all sigma points used in the UKF are all positive when κ ≥ 0. In general, all positive weights will guarantee better stability of the filtering algorithm. However, it should be emphasized that, in this specific example, there is no big difference between UKFs with different κ. In addition, the objective of this paper was to investigate the proposed filter incorporating the prior information. Hence, the UKF is used to denote UKF with κ = 3 − n and compare with the filters incorporating the prior information.
The inference results of the EKF, the UKF, the CKF 3 , and the CKF 5 are summarized in Table 2, all results are based on 50 Monte Carlo runs. It can be seen that all point-based Gaussian approximation filters have better performance than the EKF since the average(avg) FPR is lower and the average TPR and precision are higher than that of the EKF. Although the CKFs exhibit slightly better  filtering performance than the UKF in some runs, they are comparable in terms of TPR, FPR, and PPV. Based on the above tests, in the rest of the paper, only the UKF is used.

Comparison of the UKF and the UKF incorporating the prior information
As mentioned above, the UKF is used as a typical filter to compare the performance with and without the prior information.
Incorporating existing network information The following prior existing network information is assumed to be known: 1) gene1, gene5, and gene7 have little possibility to regulate gene2; 2) gene2, gene3, gene8 have little possibility to regulate gene7. Hence, the indicator matrix in (50) is given by The comparison of the UKF and the UKF incorporating the existing network information (denoted by UKF p1 ) with λ = 2 is shown in Table 3. It can be seen that the average TP# and TN# of the UKF p1 are both higher than those of the UKF. In addition, the average FP# and FN# of the UKF p1 are lower than those of the UKF. Hence, the UKF p1 predicts more correct links and nonlinks than the UKF. Moreover, the UKF p1 produces less incorrect links and missed links than the UKF. The average TPR and the precision of the UKF p1 are higher than those of the UKF. In addition, the average FPR of the UKF p1 is lower than that of the UKF. Hence, by using the existing network information, the inference accuracy can be improved.
The performance of UKF p1 with different λ is shown in Table 4. It is seen that the performance of UKF p1 and UKF is close when λ is small since only small regulation is imposed on the solution. When λ is large, the difference between the UKF p1 and UKF is large. In particular, the UKF p1 provides sparser solution than the UKF when λ is large. It can be seen from Table 4, the average FPR of UKF p1 decreases with the increasing of λ. The average TPR of UKF p1 , however, does not increase monotonically with the increasing of λ. The average PPR of UKF p1 increases with the increasing of λ. Hence, roughly speaking, the UKF p1 is better than the UKF when large λ is used.
To consider the strength of the links, rather than setting it to 1, e ij (in the indicator matrix E) is set to different values. Large e ij is used if the strength of the link from gene j to gene i is strong. For convenience, the UKF considering the strength of links is denoted as UKFp 1 . To compare the performance of UKFp 1 with UKF p1 , for UKFp 1 , the values of the second row in Equation (79) is multiplied by 5. The performance of UKFp 1 using different λ is given in Table 5. It can be seen from Tables 4 and 5 that the performance of UKFp 1 and UKF p1 is close when λ is small, e.g., λ = 0.1. In addition, the average TPR and FPR of UKFp 1 is smaller To consider the effect of false prior knowledge, the second row of the indicator matrix in Equation (79) is changed to [0, 1, 1, 1, 0, 1, 0, 1], which conflicts with the truth. For convenience, we use UKFp 1 to denote the UKF incorporating this false prior knowledge. In Table 6, the performance of UKFp 1 with different λ is shown. It can be seen from Tables 4 and 6 that the average TPR of UKFp 1 is smaller than that of UKF p1 when λ is small, e.g., λ = 0.1, 0.5. In addition, the average FPR of UKFp 1 is larger than that of UKF p1 when λ is large, e.g., λ = 5, 10. Moreover, although the average PPR of UKFp 1 is close to that of UKF p1 when λ is small, the average PPR of UKFp 1 is consistently lower than that of UKF p1 . Hence, as expected, the results indicate that the false prior knowledge will lead to worse inference result.

Incorporating LASSO
The problem setup is the same as before except that the LASSO rather than the existing network information is used. The UKF incorporating LASSO is denoted as UKF p2 .
As shown in Table 3, the average TP# and FP# of UKF p2 are lower than those of UKF and the average TN# and FN# of UKF p2 are higher than those of UKF. Hence, UKF p2 produces less links, including correct and incorrect ones. In addition, UKF p2 produces more nonlinks and missed links. It is consistent with the fact that the LASSO tends to provide a sparse solution. It can be seen from Table 3 that the average FPR of UKF p2 is lower than that of UKF and the average precision of UKF p2 is higher than that of UKF. Hence, by incorporating LASSO, the inference accuracy is improved.
A representative inference result of UKF p2 and the true regulations are shown in Figure 1. For comparison, the inference result of UKF and the true regulations are shown in Figure 2. By comparing Figure 2 and Figure 1, it can be seen that UKF falsely predicts the nonlinks from gene1 to gene2, from gene3 to gene6, from gene4 to gene8, from gene5 to gene2, and from gene6 to gene4 while UKF p2 does not.
The performance UKF p2 with different λ is shown in Table 7. It is seen that the performance of UKF p2 and UKF is close when λ is small since only small regulation is imposed on the solution. When λ is large, the difference between UKF p2 and UKF is large. The average TPR and FPR of UKF p2 decrease with the increasing of the λ. The average PPR does not increase monotonicallly with the increasing of λ. Generally speaking, for different λ, UKF p2 is more sensitive than that of UKF p1 . Although the performance of UKF p2 depends on λ, the average PPR of UKF p2 is consistently higher than that of UKF. Hence, roughly speaking, UKF p2 has better performance than UKF.

Incorporating the range constraint
The existing network information can be used to provide the rough range constraint ofx. A tight constraint is forced on the regulation coefficient a ij when there is a small regulation possibility from genej to genei and a loose constraint is forced on the regulation coefficient with no prior information. In the simulation, for the coefficients corresponding to the zero elements in (79), the lower bound and the upper bound are set as −10 and 10, respectively. For the coefficients corresponding to the nonzero elements in (79), the lower bound and the upper bound are set as −0.1 and 0.1, respectively. The UKF incorporating the range constraint is denoted as UKF p3 . As shown in Table 3, the average FPR of UKF p3 is lower than that of UKF and the average precision of UKF p3 is higher than that of UKF.

Yeast protein synthesis network
In this section, time-series gene expression data of the yeast protein synthesis network is used. Five genes (HAP1, CYB2, CYC7,CYT1, and COX5A) of the yeast protein synthesis network are considered and 17 data points which can be found in [35] are collected. The regulation relationship between them has been revealed by the biological experiment and shown in Figure 3. The dashed arrow in Figure 3 denotes 'repression' and the solid arrow denotes 'activation.' The GRN is inferred by the UKF and UKF p2 . The predicted gene expressions using parameters estimated by UKF p2 and the true measured gene expressions are shown in Figure 4. It can be seen that the model output fits the measured data well. The variances of the regulatory coefficients of HAP1 (P 1i (1 ≤ i ≤ 5)) are shown in Figure 5. It can be seen that the filter converges since the variance P 1i approaches zero. The results for other regulatory coefficients are similar and not shown here. The evaluation of the inferred GRN by UKF and UKF p2 is shown in Table 8.
By incorporating the sparsity constraint, UKF p2 provides much better inference results than UKF. As shown in Table 8, the TP# and TN# of UKF p2 are higher than those of UKF and the FP# and FN# are lower than those of UKF. In addition, it can be seen from Table 8, the FPR of UKF p2 is lower than that of UKF and the TPR and the precision of UKF p2 is higher than that of UKF.

Conclusions
In this paper, we have proposed a framework of employing the point-based Gaussian approximation filters which incorporates the prior knowledge to infer the gene regulatory network (GRN) based on the gene expression data. The performance of the proposed framework is tested by a synthetic network and the yeast protein synthesis network. Numerical results show that the inference accuracy of the GRN by the proposed point-based Gaussian approximation filter incorporating the prior information is higher than using the traditional filters without incorporating prior knowledge. The proposed method works for small-and medium-size GRNs due to the computational complexity considerations. It remains a future research topic how to adapt the proposed inference framework to handle large GRNs at reasonable computational complexity.