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.

1 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 Markov Chain Monte Carlo (MCMC) methods for the dynamic Bayesian network model [6] and the ordinary differential equation model [7], as well as the Kalman filtering methods for the state-space model [4,8] and the ordinary differential equation model [3]. Some survey papers can be found in [9-12].

Due to the ‘stochastic’ nature of the gene expression, the Kalman filtering approach based on the state-space model is one of the most competitive methods for inferring the GRN. The Kalman filter is optimal for linear Gaussian systems. However, the GRN is generally highly nonlinear. Hence, advanced filtering methods for nonlinear dynamic systems should be considered. The extended Kalman filter (EKF) is probably the most widely used nonlinear filter which uses the first-order Taylor series expansion to linearize the nonlinear model. However, the accuracy of the EKF is low when the system is highly nonlinear or contains large uncertainty. The point-based 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.

2 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 = [ x1,,…,x] is the state vector and x denotes the gene expression level of the i-th gene at time k.  is a nonlinear function that characterizes the regulatory relationship among the genes.  is the state noise and it is assumed to follow a Gaussian distribution with mean and covariance matrix , i.e., .

Following [8], we use the following nonlinear function in the state Equation (1):

with

and

In (2),  is the regulatory coefficient matrix with a denoting the regulation coefficient from gene j to gene i. Note that a positive coefficient a indicates that gene j activates gene i and a negative a indicates that gene j represses gene i. In (4), μ is a parameter. Note that  and μ 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]. 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  (·) is some nonlinear function,  is the measurement noise, which is assumed to follow the Gaussian distribution with mean  and covariance matrix , i.e., . For example, if the noise corrupted expression levels are observed, then () = .

3 Network inference using point-based Gaussian approximation filters

3.1 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 .

The optimal Bayesian filter is composed of two steps: prediction and filtering. Specifically, given the prior pdf at time k - 1, the predicted conditional pdf 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 ,  has a Gaussian distribution, i.e., . The second assumption is that (,) are jointly Gaussian given .

It then follows from the second assumption that given ,  has a Gaussian distribution, i.e., Using (1) and the first assumption, we have the predicted mean and covariance  given respectively by

and

where denotes the multivariate Gaussian pdf with mean and covariance .

Then, following the second assumption, given ,  is Gaussian distributed, i.e., Using the conditional property of the multivariate Gaussian distribution, the filtered mean and covariance  are given respectively by

with

3.2 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 can be used to approximate the integral

then the general Gaussian-type integral can be approximated by

where =  and  can be obtained by Cholesky decomposition or singular value decomposition (SVD).

Using (17), we can then approximate (8) and (9) as follows:

and

where is the transformed quadrature point obtained from the covariance decomposition, i.e.,

Similarly, we can approximate (12), (14) and (15) as follows:

where 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 third-degree quadrature filters, such as the UKF and the third-degree cubature Kalman filter (CKF3), when the system is highly nonlinear or contains large uncertainty [16]. In this paper, we consider the UKF, CKF3, and the fifth-degree cubature Kalman filter (CKF5). Other filters such as the central difference filter [14] and divided difference filter [28] can also be used. The CKF5 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.

In the unscented Kalman filter (UKF), we have N = 2 n + 1 where n is the dimension of . The quadrature points and the corresponding weights are given respectively by

and

where κ is a tunable parameter, and  is the i-th n-dimensional unit vector in which the i-th element is 1 and other elements are 0.

The left-hand side of (16) can be rewritten as

Consider the integral

By letting  = r with  = 1 and , I() can be rewritten in the spherical-radial coordinate system as

where , and is the spherical surface measure or the area element on U.

Note that (31) contains two types of integrals: the radial integral and the spherical integral .

If the radial rule can be approximated by

and the spherical integral can be approximated by

then (31) 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 is given by

with

Moreover, the set of points is given by

3.3 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  in (2) and the parameters  = [ μ1,…,μ] in (4). The problem of network inference then becomes to estimate  and . To do that, we incorporate the unknown parameters  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  and .

Specifically, we denote and the augmented state vector . Then, the augmented state equation can be written as

Note that  and  can be obtained from , and with , where  denotes an m × m all-zero matrix.

In the remainder of this paper, we assume that the noisy gene expression levels are observed. Therefore, the augmented measurement equation becomes

where , denotes an n × (n2 + n) all zeros matrix.

The point-based Gaussian approximation filters can then be used to obtain the estimate of the augmented state, , 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) become

which are the same as the filtering updates for Kalman filters.

4 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. And to accommodate range constraints, we employ PDF-truncated Gaussian approximation filters.

4.1 Optimization-based approach for sparsity constraints

Note that under the Gaussian assumption, the state estimation of the Kalman filter is equivalently given by the solution to the following optimization problem [30,31]

To incorporate the prior information of the GRN, (46) is modified as

where 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 , i.e.,

Therefore, in this case, where c1 and c2 are constants. And, (47) can be rewritten as

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 L1 - norm constraint to the GRN so that the regulatory coefficient matrix  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] where e = 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 should incur if e = 1. Thus, (47) can be rewritten as

Note that as in [20], here we do not force a = 0 corresponding to e = 1 but rather use an L1 - 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 entry-wise product operation of two matrices.

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.

For convenience, we consider a general optimization problem of the form

where and is a smooth function. Note that if , then (51) becomes (49); and if , then (51) becomes (50). Note that .

The solution to (51) can be iteratively obtained by solving a sequence of optimization problems. As in Newton’s method, the Taylor series expansion of around the solution at the t-th iteration is given by

where ∇L is the gradient of L and α is such that α mimics the Hessian ∇2L. Then, is given by [32]

The equivalent form of (53) is given by [32]

The solution to (54) is given by [32], where

is the soft thresholding function with sign() and being component-wise operators.

Finally, the iterative procedure for solving (51) is given by

And the iteration stops when the following condition is met

where ε is a given small number.

4.2 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 vector

The PDF truncation method [31] can be employed to incorporate the above range constraint into the Gaussian approximation filters, by converting the updated mean and covariance  to the pseudo mean and covariance which are then used in the next prediction and filtering steps.

We next briefly outline the PDF truncation procedure. We use and to denote the mean and covariance after the first i constraints have been enforced. Initially, we set and . Consider the following transformation

where and are obtained from the Jordan canonical decomposition and is obtained by using the Gram-Schmidt orthogonalization and it satisfies [33]

Then, the upper bound is transformed to [33]

Similarly, the lower bound is transformed to

The constraint requires that the first element of  lies between and . Hence, only the truncated PDF of the first element of  is considered and it is given by [33]

Then, the mean and variance of the first element of  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 by

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 by and .

5 Numerical results

5.1 Synthetic network

In this section, a synthetic network that contains eight genes is used to test the performance of the EKF, the UKF, the CKF3, the CKF5, and their corresponding filters incorporating the prior information. Forty data points are collected to infer the structure of the network. The system noise and measurement noise are assumed to be Gaussian distributed with means  and covariances and  = 0.01 8, respectively. The connection coefficient matrix is given by

and μ = 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 μ 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].

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.

Table 1

Comparison of UKF with different

	True positive rate						False positive rate			Positive predictive rate
Filters	Min		Max		Avg		Min	Max	Avg	Min	Max	Avg
UKF(κ = -5)	0.7576		0.9355		0.8472		0.5000	0.7647	0.5955	0.5094	0.6279	0.5824
UKF(κ = -2)	0.7576		0.9355		0.8406		0.5161	0.7647	0.5933	0.5094	0.6279	0.5825
UKF(κ = 0)	0.7576		0.9375		0.8426		0.5161	0.7647	0.5918	0.5094	0.6364	0.5840
UKF(κ = 2)	0.7576		0.9375		0.8407		0.5152	0.7353	0.5895	0.5098	0.6279	0.5841
UKF(κ = 5)	0.7576		0.9063		0.8394		0.5161	0.7353	0.5933	0.5192	0.6279	0.5821

The inference results of the EKF, the UKF, the CKF3, and the CKF5 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.

Table 2

Comparison of different filters

		True positives #			False positives #			True negatives #			False negatives #
Filters	Min	Max	Avg	Min	Max	Avg	Min	Max	Avg	Min	Max	Avg
EKF	2	17	10.60	23	44	36.4	2	15	7.08	2	24	9.92
UKF	25	29	26.80	16	26	19.28	8	16	13.06	2	8	4.86
CKF3	25	30	26.74	16	26	19.10	8	15	13.14	2	8	5.02
CKF5	25	29	26.64	16	26	19.24	8	16	13.08	1	8	5.04
	True positive rate						False positive rate			Positive predictive rate
Filters	Min		Max		Avg		Min	Max	Avg	Min	Max	Avg
EKF	0.0769		0.8667		0.5224		0.6053	0.9545	0.8358	0.0800	0.3208	0.2231
UKF	0.7576		0.9355		0.8472		0.5	0.7576	0.5955	0.5094	0.6279	0.5824
CKF3	0.7576		0.9375		0.8426		0.5161	0.7647	0.5918	0.5094	0.6364	0.5840
CKF5	0.7576		0.9667		0.8417		0.5000	0.7647	0.5946	0.5094	0.6279	0.5814

Based on the above tests, in the rest of the paper, only the UKF is used.

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) with λ = 2 is shown in Table 3. It can be seen that the average TP # and TN # of the UKF are both higher than those of the UKF. In addition, the average FP # and FN # of the UKF are lower than those of the UKF. Hence, the UKF predicts more correct links and nonlinks than the UKF. Moreover, the UKF produces less incorrect links and missed links than the UKF. The average TPR and the precision of the UKF are higher than those of the UKF. In addition, the average FPR of the UKF is lower than that of the UKF. Hence, by using the existing network information, the inference accuracy can be improved.

Table 3

Inferred results of the conventional filter and filters incorporating the prior information

		True positives #			False positives #			True negatives #			False negatives #
Filters	Min	Max	Avg	Min	Max	Avg	Min	Max	Avg	Min	Max	Avg
UKF	25	29	26.80	16	26	19.28	8	16	13.06	2	8	4.86
UKFp1	25	29	27.34	14	19	16.52	13	18	15.72	2	8	4.42
UKFp2	23	26	24.16	13	16	13.86	16	18	17.20	7	10	8.78
UKFp3	25	29	26.70	12	24	17.50	9	19	14.50	3	8	5.30
	True positive rate						False positive rate			Positive predictive rate
Filters	Min		Max		Avg		Min	Max	Avg	Min	Max	Avg
UKF	0.7576		0.9355		0.8472		0.5	0.7647	0.5955	0.5094	0.6279	0.5824
UKFp1	0.7576		0.9355		0.8614		0.4375	0.5935	0.5121	0.5778	0.6744	0.6239
UKFp2	0.6970		0.7879		0.7335		0.4194	0.5000	0.4462	0.5897	0.6667	0.6355
UKFp3	0.7576		0.9063		0.8348		0.3871	0.7273	0.5463	0.5294	0.6923	0.6049

The performance of UKF with different λ is shown in Table 4. It is seen that the performance of UKF and UKF is close when λ is small since only small regulation is imposed on the solution. When λ is large, the difference between the UKF and UKF is large. In particular, the UKF provides sparser solution than the UKF when λ is large. It can be seen from Table 4, the average FPR of UKF decreases with the increasing of λ. The average TPR of UKF, however, does not increase monotonically with the increasing of λ. The average PPR of UKF increases with the increasing of λ. Hence, roughly speaking, the UKF is better than the UKF when large λ is used.

Table 4

Comparison of UKF using different

	True positive rate						False positive rate			Positive predictive rate
Filters	Min		Max		Avg		Min	Max	Avg	Min	Max	Avg
UKFp1 (λ  = 0.1)	0.7576		0.9355		0.8484		0.5000	0.7647	0.5900	0.5094	0.6279	0.5850
UKFp1 (λ = 0.5)	0.7576		0.9677		0.8535		0.4688	0.7647	0.5696	0.5094	0.6512	0.5948
UKFp1 (λ = 1)	0.7576		0.9355		0.8614		0.4375	0.5935	0.5121	0.5778	0.6744	0.6239
UKFp1 (λ = 5)	0.7500		0.9355		0.8439		0.3548	0.5455	0.4672	0.5814	0.7105	0.6456
UKFp1 (λ = 10)	0.7273		0.9063		0.8217		0.3226	0.4848	0.4156	0.6190	0.7368	0.6695

To consider the strength of the links, rather than setting it to 1, e (in the indicator matrix ) is set to different values. Large e 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 . To compare the performance of with UKF, for , the values of the second row in Equation (79) is multiplied by 5. The performance of using different λ is given in Table 5. It can be seen from Tables 4 and 5 that the performance of and UKF is close when λ is small, e.g., λ = 0.1. In addition, the average TPR and FPR of is smaller than that of UKF for all tested λ except for λ = 0.1. Hence, PPR is used to evaluate the performance of and UKF. Although the average PPR of and UKF is close when the λ is large, e.g., λ = 10, the average PPR of is consistently higher than that of UKF. The results indicate that the inference accuracy of and UKF are close when λ is very small or very large. The inference accuracy of outperforms UKF when the appropriate strength of the link and parameter λ are used.

Table 5

Effect of strength of the links using different

	True positive rate						False positive rate			Positive predictive rate
Filters	Min		Max		Avg		Min	Max	Avg	Min	Max	Avg
UKF_p~1 (λ = 0.1)	0.7576		0.9677		0.8484		0.4688	0.7647	0.5713	0.5094	0.6512	0.5929
UKF_p~1 (λ = 0.5)	0.7576		0.9333		0.8468		0.4516	0.7059	0.5422	0.5385	0.6512	0.6057
UKF_p~1 (λ = 1)	0.7500		0.9032		0.8221		0.3750	0.5758	0.4953	0.5814	0.6842	0.6257
UKF_p~1 (λ = 5)	0.7273		0.8750		0.8220		0.3548	0.5000	0.4169	0.6098	0.7179	0.6684
UKF_p~1 (λ = 10)	0.7500		0.8750		0.8214		0.3226	0.5000	0.4143	0.6098	0.7368	0.6696

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 to denote the UKF incorporating this false prior knowledge. In Table 6, the performance of with different λ is shown. It can be seen from Tables 4 and 6 that the average TPR of is smaller than that of UKF when λ is small, e.g., λ = 0.1,0.5. In addition, the average FPR of is larger than that of UKF when λ is large, e.g., λ = 5,10. Moreover, although the average PPR of is close to that of UKF when λ is small, the average PPR of is consistently lower than that of UKF. Hence, as expected, the results indicate that the false prior knowledge will lead to worse inference result.

Table 6

Effect of false prior information using different

	True positive rate						False positive rate			Positive predictive rate
Filters	Min		Max		Avg		Min	Max	Avg	Min	Max	Avg
UKF_p¯1 (λ = 0.1)	0.7576		0.9667		0.8491		0.5000	0.7647	0.5933	0.5094	0.6279	0.5835
UKF_p¯1 (λ = 0.5)	0.7576		0.9355		0.8535		0.4839	0.7647	0.5962	0.5094	0.6279	0.5836
UKF_p¯1 (λ = 1)	0.7576		0.9333		0.8572		0.4839	0.7059	0.6001	0.5200	0.6279	0.5830
UKF_p¯1 (λ = 5)	0.6970		0.8125		0.7546		0.4194	0.5938	0.5000	0.5682	0.6486	0.6062
UKF_p¯1 (λ = 10)	0.5758		0.7576		0.6810		0.3226	0.5000	0.4066	0.5676	0.7059	0.6369

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.

As shown in Table 3, the average TP # and FP # of UKF are lower than those of UKF and the average TN # and FN # of UKF are higher than those of UKF. Hence, UKF produces less links, including correct and incorrect ones. In addition, UKF 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 is lower than that of UKF and the average precision of UKF is higher than that of UKF. Hence, by incorporating LASSO, the inference accuracy is improved.

A representative inference result of UKF 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 does not.

Figure 1

Inferred regulations of UKF and true regulations.

Figure 2

Inferred regulations of UKF and true regulations.

The performance UKF with different λ is shown in Table 7. It is seen that the performance of UKF and UKF is close when λ is small since only small regulation is imposed on the solution. When λ is large, the difference between UKF and UKF is large. The average TPR and FPR of UKF decrease with the increasing of the λ. The average PPR does not increase monotonicallly with the increasing of λ. Generally speaking, for different λ, UKF is more sensitive than that of UKF. Although the performance of UKF depends on λ, the average PPR of UKF is consistently higher than that of UKF. Hence, roughly speaking, UKF has better performance than UKF.

Table 7

Comparison of UKF using different

	True positive rate						False positive rate			Positive predictive rate
Filters	Min		Max		Avg		Min	Max	Avg	Min	Max	Avg
UKFp2 (λ = 0.1)	0.7576		0.9355		0.8304		0.4839	0.6970	0.5699	0.5306	0.6512	0.5914
UKFp2 (λ = 0.5)	0.6970		0.8710		0.7750		0.4194	0.5758	0.4902	0.5682	0.6585	0.6198
UKFp2 (λ = 1)	0.6970		0.7879		0.7335		0.4194	0.5000	0.4462	0.5897	0.6667	0.6355
UKFp2 (λ = 5)	0.4545		0.6667		0.5501		0.3226	0.4516	0.3791	0.5714	0.6471	0.6064
UKFp2 (λ = 10)	0.4545		0.5455		0.4800		0.2903	0.3871	0.3523	0.5556	0.6538	0.5920

Incorporating the range constraint The existing network information can be used to provide the rough range constraint of . A tight constraint is forced on the regulation coefficient a 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. As shown in Table 3, the average FPR of UKF is lower than that of UKF and the average precision of UKF is higher than that of UKF.

5.2 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.’

Figure 3

Pathway model of the five genes in yeast protein synthesis network.

The GRN is inferred by the UKF and UKF. The predicted gene expressions using parameters estimated by UKF 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 (P1 (1 ≤ i ≤ 5)) are shown in Figure 5. It can be seen that the filter converges since the variance P1 approaches zero. The results for other regulatory coefficients are similar and not shown here. The evaluation of the inferred GRN by UKF and UKF is shown in Table 8.

Figure 4

True gene expression and model output.

Figure 5

Variance of regulatory coefficients.

Table 8

Inferred results of the UKF and UKF

Filters	True positives #	False positives #	True negatives #	False negatives #
UKF	1	7	14	3
UKFp2	2	3	18	2
Filters	TPR		FPR	Precision
UKF	0.25		0.3333	0.1250
UKFp2	0.5000		0.1429	0.4000

By incorporating the sparsity constraint, UKF provides much better inference results than UKF. As shown in Table 8, the TP # and TN # of UKF 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 is lower than that of UKF and the TPR and the precision of UKF is higher than that of UKF.

6 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.

Competing interests

All authors declare that they have no competing interests.