Reverse Engineering Sparse Gene Regulatory Networks Using Cubature Kalman Filter and Compressed Sensing

Abstract

This paper proposes a novel algorithm for inferring gene regulatory networks which makes use of cubature Kalman filter (CKF) and Kalman filter (KF) techniques in conjunction with compressed sensing methods. The gene network is described using a state-space model. A nonlinear model for the evolution of gene expression is considered, while the gene expression data is assumed to follow a linear Gaussian model. The hidden states are estimated using CKF. The system parameters are modeled as a Gauss-Markov process and are estimated using compressed sensing-based KF. These parameters provide insight into the regulatory relations among the genes. The Cramér-Rao lower bound of the parameter estimates is calculated for the system model and used as a benchmark to assess the estimation accuracy. The proposed algorithm is evaluated rigorously using synthetic data in different scenarios which include different number of genes and varying number of sample points. In addition, the algorithm is tested on the DREAM4 in silico data sets as well as the in vivo data sets from IRMA network. The proposed algorithm shows superior performance in terms of accuracy, robustness, and scalability.

1. Introduction

Gene regulation is one of the most intriguing processes taking place in living cells. With hundreds of thousands of genes at their disposal, cells must decide which genes are to express at a particular time. As the cell development evolves, different needs and functions entail an efficient mechanism to turn the required genes on while leaving the others off. Cells can also activate new genes to respond effectively to environmental changes and perform specific roles. The knowledge of which gene triggers a particular genetic condition can help us ward off the potential harmful effects by switching that gene off. For instance, cancer can be controlled by deactivating the genes that cause it.

Gene expression is the process of generating functional gene products, for example, mRNA and protein. The level of gene functionality can be measured using microarrays or gene chips to produce the gene expression data [1]. More accurate estimation of gene expression is now possible using the RNA-Seq method. Intelligent use of such data can help improve our understanding of how the genes are interacting in a living organism [2–4]. Gene regulation is known to exhibit several modes; a couple of important ones include transcription regulation and posttranscription regulation [5]. While the theoretical applications of gene regulation are extremely promising, it requires a thorough understanding of this complex process. Different genes may cooperate to produce a particular reaction, while a gene may repress another gene as well. The potential benefits of gene regulation can only be reaped if a complete and accurate picture of genetic interactions is available. A network specifying different interconnections of genes can go a long way in understanding the gene regulation mechanism. The control and interaction of genes can be described through a gene regulatory network. Such a network depicts various interdependencies among genes where nodes of the network represent the genes, and the edges between them correspond to an interaction among them. The strength of these interactions represents the extent to which a gene is affected by other genes in the network. A key ingredient of this approach is an accurate and representative modeling of gene networks. Precise modeling of a regulatory network coupled with efficient inference and intervention algorithms can help in devising personalized medicines and cures for genetic diseases [6].

Various methods for gene network modeling have been proposed recently in the literature with varying degrees of sophistication [7–10]. These techniques can be broadly classified as static and dynamic modeling schemes. Static modeling includes the use of correlation, statistical independence for clustering [11–13], and information theoretic criteria [14–16]. On the other hand, dynamic models provide an insight into the temporal evolution of gene expressions and hence yield a more quantitative prediction on gene network behavior [17–20]. In order to incorporate the stochasticity of gene expressions, statistical techniques have been applied [13]. A rich literature is also available on the Bayesian modeling of gene networks [21–26]. Promoted in part by the Bayesian methods, the state-space approach is a popular technique to model the gene networks [27–33], whereby the hidden states can be estimated using the Kalman filter. In the case of nonlinear functions, the extended Kalman filter (EKF) and particle filter represent feasible approaches [33, 34]. However, the EKF relies on the first-order linear approximations of nonlinearities, while the particle filter may be computationally too complex. A comprehensive review of these methods can be found in [35].

In this paper, the gene network is modeled using a state-space approach, and the cubature Kalman filter (CKF) is used to estimate the hidden states of the nonlinear model [36, 37]. The gene expressions are assumed to evolve following a sigmoid squash function, whereas a linear function is considered for the expression data. The noise is assumed to be Gaussian for both the state evolution and gene expression measurements. As the gene network is assumed sparse, any simple mean square error minimization technique will not suffice for the estimation of static parameters. Therefore, a compressed sensing-based Kalman filter (CSKF) [38] is used in conjunction with CKF for reliable estimation of parameters. In case of statistical inference, it is essential to obtain some guarantees on the performance of estimators. In this regard, the Cramér-Rao lower bound (CRB) of the parameter estimates is used as a benchmarking index to assess the mean square error (MSE) performance of the proposed estimator which is evaluated here for a parameter vector. The performance of the proposed algorithm is tested on synthetically generated random Boolean networks in various scenarios. The algorithm is also tested using DREAM4 data sets and IRMA networks [39, 40].

The main contributions of this paper can be summarized as follows. (1)CKF is proposed for the estimation of states, and a compressed sensing-based Kalman filter is used for the estimation of system parameters. The genes are known to interact with few other genes only necessitating the use of sparsity constraint for more accurate estimation. The proposed algorithm carries out online estimation of parameters and is therefore computationally efficient and is particularly suitable for large gene networks.(2)The Cramér-Rao lower bound is calculated for the estimation of unknown parameters of the system. The performance of the proposed algorithm is compared to CRB. This comparison is significant as it shows room for improvement in the estimation of parameters. (3)The proposed algorithm is compared with the EKF algorithm. Using the false alarm errors, true connections, and Hamming distance as fidelity criteria, rigorous simulations are carried out to assess the performance of the algorithm with the increase in the number of samples. In addition, receiver operating characteristic (ROC) curves are plotted to evaluate the algorithms for different network sizes. It is observed that the proposed algorithm outperforms EKF in terms of accuracy and precision. The proposed algorithm is then applied to the DREAM4 10-gene and 100-gene data sets to assess the algorithm accuracy. The underlying gene network for the IRMA data sets is also inferred.

The rest of this paper is organized as follows. Section 2 describes the underlying system model for the gene expressions. The proposed CKF algorithm in combination with CSKF for gene network inference is formulated in Section 3. The derivation of CRB is shown in Section 4, and the simulation results and their interpretation are presented in Section 5. Finally, conclusions are drawn in Section 6.

2. System Model

Gene regulatory networks can be modeled as static or dynamical systems. In this work, state-space modeling is considered which is an instance of a dynamic modeling approach and can effectively cope with time variations. The states represent gene expressions, and their evolution in time, in general, can be expressed as 𝐱_𝑘=𝑔𝐱_𝑘−1+𝐰_𝑘𝑘=1,…,𝐾, (1) where 𝐾 is the total number of data points available, 𝐰_𝑘 is assumed to be a zero-mean Gaussian random variable with covariance 𝐐_𝑘=𝜎²_𝑤𝐈, and the function 𝑔(⋅) represents the regulatory relationship between the genes and is generally nonlinear. The microarray data is a set of noisy observations and is commonly expressed as a linear Gaussian model [41] 𝐲_𝑘=ℎ𝐱_𝑘+𝐯_𝑘, (2) where 𝐯_𝑘 is Gaussian-distributed random variable with zero mean and covariance 𝐒_𝑘=𝜎²_𝑣𝐈 and incorporates the uncertainty in the microarray experiments. In order to capture the gene interactions effectively, the following nonlinear state evolution model is assumed [33, 34]: 𝑥_𝑘,𝑛=^𝑁_𝑚=1𝑏_𝑛𝑚𝑓𝑥_{𝑘−1,𝑚}+𝑤_𝑘,𝑛,𝑘=1,…,𝐾,𝑛=1,…,𝑁, (3) where 𝑁 is the total number of genes in the network and 𝑓(⋅) is the sigmoid squash function 𝑓𝑥_{𝑘−1,𝑚}=11+𝑒^{−𝑥𝑘−1,𝑚}. (4) This particular choice for the nonlinear function ensures that the conditional distribution of the states remains Gaussian [41]. The multiplicative constants 𝑏_𝑛𝑚 quantify the positive or negative relations between various genes in the network. A positive value of  𝑏_𝑛𝑚 implies that the 𝑚th gene is activating the 𝑛th gene, whereas a negative value implies repression [34, 42]. The absolute value of these parameters indicates the strength of interaction.

The model given in (3) and (4) in the absence of any constraints may be unidentifiable and may result into overfitted solutions [43]. Assumptions on network structures are, therefore, necessary to obtain a connectivity matrix that agrees with the biological knowledge. In a gene regulatory network (GRN), the genes are known to interact with few other genes only. To this end, the coefficients 𝑏_𝑛𝑚s are estimated using sparsity constraints, as explained in the next section.

A discrete linear Gaussian model for the microarray data is considered which can be expressed at the 𝑘th time instant as [41] 𝐲_𝑘=𝐱_𝑘+𝐯_𝑘. (5) Stacking the unknown parameters together, the parameter vector to be estimated is 𝐛≜𝜙₁,𝜙₂,…,𝜙_𝑁, (6) where 𝜙_𝑛=[𝑏_𝑛1,…,𝑏_𝑛𝑁]. Plugging the values of states from (3) into (5), it follows that 𝐲_𝑘=𝐑_𝑘𝐛+𝐞_𝑘, (7) where 𝐑_𝑘≜𝐟_𝑘̃𝟎𝐟̃𝟎𝟎𝟎_𝑘𝐟̃𝟎𝟎𝟎𝟎⋱𝟎𝟎𝟎𝟎_𝑘𝐟̃_𝑘,≜𝑓𝑥_𝑘−1,1⋯𝑓𝑥_{𝑘−1,𝑁}. (8)(9) Thus, the gene network inference problem boils down to the estimation of system parameters 𝐛 using the observations 𝐲_𝑘, where the effective noise 𝐞_𝑘 is the sum of system and observation noises. The next section describes the proposed inference algorithm for sparse networks.

3. Method

In this section, the methodology proposed to infer the system parameters in (3) is described. The proposed cubature Kalman filter with sparsity constraints (CKFS) approach is succinctly illustrated in Figure 1. The specific details of this algorithm are as next presented.

Figure 1 

Block diagram of network inference methodology CKFS.

3.1. Cubature Kalman Filter

Kalman filter is a Bayesian filter which provides the optimal solution to a general linear state space inference problem depicted by (1) and (2) and assumes a recursive predictive-update process. The underlying assumption of Gaussianity for the predictive and the likelihood densities simplifies the Kalman filter algorithm to a two-step process, consisting of prediction and update of the mean and covariance of the hidden states. However, the presence of nonlinear functions in the state and measurement equations requires calculation of multidimensional integrals of the form nonlinear function  ×  Gaussian density [36], which in general is computationally prohibitive. Several solutions to this problem have been proposed including the EKF, which linearizes the nonlinear function by taking its first-order Taylor approximation, and the unscented Kalman filter (UKF), which approximates the probability density function (PDF) using a nonlinear transformation of the random variable. Recently, a new approach, CKF, has been proposed which evaluates the integrals numerically using spherical-radial cubature rules [36].

The next two subsections briefly explain the working of Bayesian filtering and the CKF solution for the nonlinear multidimensional integrals.

3.1.1. Time Update

Let the observations up to the time instant 𝑘 be denoted by 𝐝_𝑘; that is, 𝐝_𝑘≜[𝐲^𝑇₁,…,𝐲^𝑇_𝑘]^𝑇. In the prediction phase, also called the time update of the Bayesian filter, the mean and covariance of the Gaussian posterior density are computed as follows: 𝐱_{𝑘∣𝑘−1}=𝐸𝐟𝐱_𝑘−1∣𝐝_𝑘−1,𝐏_{𝑥𝑥,𝑘∣𝑘−1}=𝐸𝐟𝐱_𝑘𝐟^𝑇𝑥_𝑘−𝐱_{𝑘∣𝑘−1}𝐱^𝑇_{𝑘∣𝑘−1}+𝐐_𝑘−1, (10) where 𝐸 denotes the expectation operator and 𝐱_𝑘−1 is normally distributed with parameters (𝐱_{𝑘−1∣𝑘−1},𝐏_{𝑥𝑥,𝑘−1∣𝑘−1}). The third equality is a consequence of the zero-mean nature of Gaussian noise 𝐰 and its independence from 𝐝_𝑘. The estimates 𝐱_{𝑘−1∣𝑘−1} and 𝐏_{𝑥𝑥,𝑘−1∣𝑘−1} are assumed to be available from the previous iteration. Here, 𝐏_{𝑥𝑥,𝑘∣𝑘−1} is an estimate of the error covariance matrix.

3.1.2. Measurement Update

Since the measurement noise is also Gaussian, the likelihood density is given by 𝐲_𝑘−1∣𝐝_𝑘−1∶𝒩(𝐳_𝑘−1;𝐲_{𝑘∣𝑘−1},𝐏_{𝑥𝑥,𝑘∣𝑘−1}). As the measurements become available at the 𝑘th time instant, the mean and covariance of the likelihood density are calculated as follows: 𝐲_{𝑘∣𝑘−1}=𝐸𝐲_𝑘∣𝐝_𝑘−1,𝐏_{𝑦𝑦,𝑘∣𝑘−1}=𝐸𝐱_𝑘𝐱^𝑇_𝑘−𝐲_{𝑘∣𝑘−1}𝐲^𝑇_{𝑘∣𝑘−1}+𝐒_𝑘−1. (11) The updated posterior density, obtained from the conditional joint density of states, and the measurements can be expressed as 𝐱^𝑇_𝑘𝐲^𝑇_𝑘^𝑇𝐝_𝐤−𝟏𝐱∼𝒩^{𝑘∣𝑘−1}𝐲_{𝑘∣𝑘−1},𝐏_{𝑥𝑥,𝑘∣𝑘−1}𝐏_{𝑥𝑦,𝑘∣𝑘−1}𝐏^𝑇_{𝑥𝑦,𝑘∣𝑘−1}𝐏_{𝑦𝑦,𝑘∣𝑘−1}, (12) where 𝑃_{𝑥𝑦,𝑘∣𝑘−1}=𝐸𝐱_𝑘𝐱^𝑇_𝑘−𝐱_{𝑘∣𝑘−1}𝐲^𝑇_{𝑘∣𝑘−1} (13) is the cross-covariance matrix between the states and the measurements. Hence, the states and the corresponding error covariance matrix are updated by calculating the innovation 𝐳_𝑘−𝐳̂_{𝑘∣𝑘−1} and the Kalman gain 𝐊_𝐺,𝑖𝐱_𝑘∣𝑘=𝐱_{𝑘∣𝑘−1}+𝐊_𝐺,𝑘𝐲_𝑘−𝐲_{𝑘∣𝑘−1},𝑃_{𝑥𝑥,𝑘∣𝑘}=𝑃_{𝑥𝑥,𝑘∣𝑘−1}−𝐊_𝐺,𝑘𝑃_{𝑦𝑦,𝑘∣𝑘−1}𝐊^𝑇_𝐺,𝑘,𝐊_𝐺,𝑘=𝑃_{𝑥𝑦,𝑘∣𝑘−1}𝑃⁻¹_{𝑦𝑦,𝑘∣𝑘−1}. (14) The next subsection briefly describes the computation of high-dimensional integrals present in the equations above.

3.1.3. Computation of Integrals Using Spherical-Radial Cubature Points

In order to determine the expectations in (10), using a numerical integration method, a spherical-radial cubature rule is applied. This method calculates the cubature points 𝐗_{𝑗,𝑘−1∣𝑘−1} as follows [36]: 𝐗_{𝑗,𝑘−1∣𝑘−1}=𝐔_{𝑘−1∣𝑘−1}𝜁_𝑗+𝐱_{𝑘−1∣𝑘−1}, (15) where 𝜁_𝑗=ℓ/2[1]_𝑗, 𝑗=1,…,ℓ, ℓ=2𝑁 denotes the total number of cubature points and 𝐔_{𝑘−1∣𝑘−1} stands for the square root of the error covariance matrix; that is, 𝐏_{𝑥𝑥,𝑘−1∣𝑘−1}=𝐔_{𝑘−1∣𝑘−1}𝐔^𝑇_{𝑘−1∣𝑘−1}. (16) The cubature points are updated via the state equation 𝐗^∗_{𝑗,𝑘∣𝑘−1}=𝑔𝐗_{𝑗,𝑘−1∣𝑘−1}. (17) The propagated cubature points yield the state and error covariance estimates 𝐱_{𝑘∣𝑘−1}=1ℓ^ℓ_𝑗=1𝐗^∗_{𝑗,𝑘∣𝑘−1},𝐏_{𝑥𝑥,𝑘∣𝑘−1}=1ℓ^ℓ_𝑗=1𝐗^∗_{𝑗,𝑘∣𝑘−1}𝐗^∗𝑇_{𝑗,𝑘∣𝑘−1}−𝐱_{𝑘∣𝑘−1}𝐱^𝑇_{𝑘∣𝑘−1}+𝐐_𝑘−1. (18) The integrals in (11) and (14) can be evaluated in a similar manner. The next subsection explains the estimation of parameters in the system.

3.2. Estimation of Sparse Parameters Using Kalman Filter

The state estimates are obtained using the CKF as described in the previous subsection. In order to estimate the unknown parameters in the system model, one of the most commonly used methods involves stacking the parameters with the states and estimating them together. The estimation process performed in this manner is called joint estimation. Another method for the estimation of parameters consists of a two-step recursive process which is termed dual estimation. This process estimates the states in the first step, and with the assumption that states are known, parameters are estimated in the second step. These steps are repeated until the algorithm converges to the true values or until the amount of available observations is exhausted. This paper makes use of the latter technique.

The vector 𝐛 as defined in (6) is assumed to be evolving as a Gauss-Markov model. As discussed previously, the states are assumed to be known at this step. The system evolution equations can therefore be expressed as 𝐛_𝑘=𝐛_𝑘−1+𝜼_𝑘−1,𝐲_𝑘=𝐑_𝑘𝐛_𝑘+𝐞_𝑘, (19) where 𝜼_𝑘 stands for the i.i.d Gaussian noise and 𝐑_𝑘 is as defined in (8). It is observed that (19) is a system of linear equations with additive Gaussian noise, and therefore, the Kalman filter is the optimal choice for the estimation of parameter vector. The standard predict and update steps involved in Kalman filter are summarized as follows: 𝐛_{𝑘∣𝑘−1}=𝐛_{𝑘−1∣𝑘−1}+𝜼_𝑘,𝐏_{𝑏𝑏,𝑘∣𝑘−1}=𝐏_{𝑏𝑏,𝑘−1∣𝑘−1}+𝚺_𝜂𝑘,𝐮_𝑘=𝐲_𝑘−𝐑_𝑓𝑘𝐛_𝑘,𝐊_𝐺=𝐏_{𝑏𝑏,𝑘∣𝑘−1}𝐑^𝑇_𝑓𝑘𝐑_𝑓𝑘𝐏_{𝑏𝑏,𝑘∣𝑘−1}𝐑^𝑇_𝑓𝑘+𝜎²_𝑒𝐈⁻¹,𝐛_𝑘∣𝑘=𝐛_{𝑘∣𝑘−1}+𝐊_𝐺𝐮_𝑘,𝐏_{𝑏𝑏,𝑘∣𝑘}=𝐈−𝐊_𝐺𝐑_𝑓𝑘𝐏_{𝑏𝑏,𝑘∣𝑘−1}, (20) where 𝐊_𝐺 denotes the Kalman gain and 𝐏 represents the error covariance matrix.

The Kalman filter algorithm is based on an 𝑙₂-norm minimization criterion. As the gene networks are known to be highly sparse, the parameter vector is expected to have only a few nonzero values. A more accurate approach for estimating such a vector would be to introduce an additional constraint on its 𝑙₁-norm which is the core idea in compressed sensing [38, 44]. Such an 𝑙₁-norm constraint provides a unique solution to the underdetermined set of equations [45]. Therefore, instead of a simple 𝑙₂ norm minimization, the following constrained optimization problem is considered: 𝐛min_𝑘𝐛_𝑘−𝐛_𝑘²₂s.t.𝐛_𝑘≤𝜖. (21) The importance of this constraint can be judged by the fact that without it, the system would be rendered unidentifiable [43].

The problem (21) can be solved using a pseudomeasurement (PM) method which incorporates the inequality constraint (21) in the filtering process by assuming an artificial measurement ‖𝐛_𝑘‖₁−𝜖=0. This is concisely expressed as 0=𝐑𝐛_𝑘−𝜖,𝐑_𝜏=sign𝐛_𝜏(1),…,sign𝐛_𝜏(𝑁). (22) The value of the covariance matrix 𝚺_𝜖=𝜎²_𝜖𝐈 of the pseudonoise 𝜖 is selected in a similar manner as the process noise covariance in the EKF algorithm. However, it is found that large values of variances, that is, 𝜎²_𝜖≥100, prove sufficient in most cases [38]. Further details on selecting these parameters can be found in [38, 46]. The PM method solves (21) in a recursive manner for 𝐾_𝜏 iterations using the following steps: 𝐊^𝜏_𝐺=𝐏_𝜏𝐑^𝑇_𝜏𝐑_𝜏𝐏_𝜏𝐑^𝑇_𝜏+𝚺_𝜖⁻¹,𝐛_𝜏+1=𝐈−𝐊^𝜏_𝐺𝐑_𝜏𝐛_𝜏,𝐏_𝜏+1=𝐈−𝐊^𝜏_𝐺𝐑_𝜏𝐏_𝜏. (23) At each 𝑘th time instant, 𝐏_{𝑏𝑏,𝑘∣𝑘} and 𝐛_𝑘∣𝑘 obtained from (20) are considered as initial values; that is, 𝐛¹=𝐛_𝑘∣𝑘 and 𝐏₁=𝐏_{𝑏𝑏,𝑘∣𝑘} which is the error covariance matrix. The value of 𝐾_𝜏 is equal to the number of constraints, that is, the expected number of nonzero 𝐛_𝑚𝑛s in the system. Possible ways for calculating 𝐾_𝜏 include minimum description length (MDL) principle and Bayesian information criterion (BIC).

3.3. Inference Algorithm

The network inference algorithm is summarized in Algorithm 1. The algorithm consists of a recursive process which repeats itself for the number of observations present in the time-series data. For each time sample, the state estimate is obtained using the CKF, and the parameter estimate is obtained using the KF. Since the parameters are expected to be sparse, the estimates are then refined further using the CSKF algorithm. This iterative process results in a simple and accurate algorithm for gene network inference while considering a complex nonlinear model.

(1)    Input time series data set 𝐲. (2)   Initialize 𝐼,𝐾,𝜙0,𝐱0. (3)   for  𝑘=1,…,𝐾  do (4)    Find the state estimates using CKF following the time       and measurement update steps in Section 3. (5)    Estimate parameters 𝐛𝑘 from 𝐱𝑘 and 𝐲𝑘 using (20). (6)    for  𝜏=1,…,𝐾𝜏  do (7)     Update the parameters 𝐛𝑘 using (23). (8)    end  for (9)   end  for (10) return

Algorithm 1 

Network inference: CKFS.

4. Cramér-Rao Bound

The performance of an estimator can be judged by comparing it with theoretical lower bounds proposed in parameter estimation theory. The CRB establishes a lower bound on the MSE of an unbiased estimator [47]. In particular, the CRB states that the covariance matrix of the estimator 𝐛 is lower bounded by 𝔼𝐛−𝐛𝐛−𝐛^𝑇≽[𝐈(𝐛)]⁻¹, (24) where the matrix inequality ≽ is to be interpreted in the semidefinite sense and 𝐈(𝐛) is the Fisher information matrix (FIM) 𝐈(𝐛)()=𝔼𝜕ln𝑓𝐲∣𝐛()𝜕𝐛𝜕ln𝑓𝐲∣𝐛𝜕𝐛^𝑇. (25) The CRB for gene network inference can be calculated as follows. By stacking all the observations for 𝑘=1,…,𝐾, (7) can be written compactly in the matrix form 𝐲=𝐑𝐛+𝐞, (26) where 𝐲=[𝐲^𝑇₁,…,𝐲^𝑇_𝐾]^𝑇, 𝐑=[𝐑^𝑇₁,…,𝐑^𝑇_𝐾]^𝑇, and 𝐞=[𝐞^𝑇₁,…,𝐞^𝑇_𝐾]^𝑇. The PDF 𝑝(𝐲∣𝐛) is expressed as 𝑝()−()𝐲∣𝐛=𝐶exp𝐲−𝐑𝐛^𝑇(𝐲−𝐑𝐛)2𝜎²_𝑒, (27) where 𝐶 is a constant. The derivative of ln𝑝(𝐲∣𝐛) can be expressed as 𝜕ln𝑝(𝐲∣𝐛)𝜕𝜕𝐛=−()𝜕𝐛𝐲−𝐑𝐛^𝑇(𝐲−𝐑𝐛)𝜎²_𝑒=𝐑^𝑇𝐲−𝐑^𝑇𝐑𝐛𝜎²_𝑒. (28) It now follows that 𝜕ln𝑝(𝐲∣𝐛)()𝜕𝐛𝜕ln𝑝𝐲∣𝐛𝜕𝐛^𝑇=𝐑^𝑇(𝐲−𝐑𝐛)(𝐲−𝐑𝐛)^𝑇𝐑𝜎⁴_𝑒. (29) By taking the expectation of (29), the FIM in (25) is given by 𝐈(𝐛)=𝐑^𝑇𝐑𝜎²_𝑒. (30) The inverse of the FIM in (30) can be used to place a lower bound on the estimation error of the parameter vector 𝐛. Figure 2 shows the comparison of MSE of CKFS algorithm with CRB as a function of number of samples 𝐾 for one representative gene from the eight-gene network considered in Section 5.1. It is observed that the MSE of the estimated parameters decreases with increasing number of samples.

Figure 2 

Cramér-Rao bound on the estimation of parameters. The MSE for one of the representative 𝜃 is shown here for a network consisting of 8 vertices.

5. Results and Discussion

The simulation results of the CKFS algorithm are discussed in this section. The performance is first tested on synthetic data obtained from randomly generated Boolean networks under various scenarios and performance metrics. The algorithm is then assessed on the DREAM4 networks and the IRMA network.

5.1. Synthetic Data

Time-series data is produced from randomly generated Boolean networks using the system model (3) and (5). Two scenarios are considered for this purpose.

First, the comparison is performed by varying the number of sample size while keeping the network size fixed. The gene network consists of 8 genes and 20 vertices. In terms of network estimation, if the algorithm predicts an edge between two nodes which may not be present in reality, an error, referred to as false alarm error (F), is said to have occurred. Another situation is the indication of the absence of a vertex in the graph which in fact is present in the real network. This kind of error is termed missed detection (M). The summation of these two errors normalized over the total number of vertices in the network yields the Hamming distance. It is also important to consider the probability of predicting the true connections correctly which will be assessed by the true connections (T) metric. An algorithm with low Hamming distance and small false alarm error is particularly desirable as predicting an edge erroneously can be troublesome for biologists. True connections indicate the reliability of the predictions. Figure 3 illustrates the performance of the CKFS algorithm and that of the EKF algorithm proposed in [34] in terms of the metrics described above. It is important to mention here that the same system model is assumed by both CKFS and EKF algorithms for the purpose of this simulation. These metrics are the same as those used in [15]. The variances of both the system and measurement noises, 𝜎²_𝑤 and 𝜎²_𝑣, respectively, are taken to be 10⁻⁵ in all the simulations and are assumed to be known. It is noticed that EKF has a slightly lower false alarm rate when the number of samples is small; however, as the number of samples increases, CKFS yields a lower false alarm error. The Hamming distance for CKFS is also smaller than EKF indicating lesser cumulative error. True connections show a consistent behavior for the two algorithms when the number of samples is increased where CKFS is able to predict connections more accurately. These experiments show the superiority of CKFS in terms of lower error rate.

(a)

(b)

(c)

Figure 3 

(a), (b), and (c) False alarm errors, Hamming distance, and true connections. The synthetic networks consist of 8 vertices and 20 edges. The metric is normalized over the number of edges. CKFS gives lower error and predicts more true connections with the increase in the sample size of data.

To obtain a more rigorous evaluation, the performance of algorithms is then compared in a scenario which considers the sample size to be fixed and assumes networks of different sizes. The receiver operating characteristic (ROC) curves are plotted as performance measures. A higher area under the ROC curve (AUROC) shows more true positives for a given false positive, and therefore, indicates better classification. The performance of CKFS(𝑁,𝐸,𝐾) and EKF(𝑁,𝐸,𝐾) is shown in Figure 4, where 𝑁 stands for the number of nodes, 𝐸 represents the number of edges, and 𝐾 denotes the time points. It is observed that the CKFS exhibits superior performance than the EKF for networks of different sizes.

(a)

(b)

(c)

Figure 4 

ROC curves for the performance of CKFS and EKF using synthetic data. (𝑁,𝐸,𝐾) (a), (b), and (c) (5,10,20), (10,12,20), and (15,19,20). The area under the ROC curve for CKFS is more than that for EKF for various sized networks.

The complexity of the two algorithms is compared for synthetically generated networks with number of genes equal to 10, 20, 30, and 40. The sample size is kept to 50 time points for each of these networks, and the run time for EKF and CKFS algorithms is calculated as shown in Table 1. It is noted that EKF is faster for smaller network sizes, but as the network size increases, the run time gets much larger than that for CKFS. The main reason for this is that EKF [34] estimates the states and parameters by stacking them together which requires large-sized matrix multiplications at each iteration. The benefit associated with performing dual estimation, as in CKFS, is that the parameters are estimated separately from the states. Since the system is linear and one-to-one for parameters, inversion of much smaller matrices can be performed reducing the computational complexity of CKFS algorithm. CKFS is therefore particularly attractive for large-sized networks.

Table 1 

Run time in seconds for EKF and CKFS algorithms for varying network sizes for synthetically generated data. The number of sample points is fixed to 50.

5.2. DREAM4 Gene Networks

Several in silico networks have been produced in order to benchmark the performance of gene network inference algorithms. dialogue on reverse engineering assessment and methods (DREAM) in silico networks serve as one of the popular methods used for this purpose [39, 48]. In this section, the performance of the CKFS algorithm is evaluated using the 10-gene and 100-gene networks released online by the DREAM4 challenge. Five networks are produced using the known GRNs of Escherichia coli and Saccharomyces cerevisiae. The data sets for each of 10-gene network consists of 21 data points for five different perturbations. The inference is performed by using all the perturbations. The 100-gene network consists of data sets for ten perturbations. AUROC and area under the precision-recall curve (AUPR) are calculated for the five networks of both the data sets and shown in Tables 2 and 3, respectively. The quantities, precision and recall, are defined as 𝑃=𝑇/(𝑇+𝐹) and 𝑅=𝑇/(𝑇+𝑀), respectively. For comparison purposes, the values of the two quantities for time-series network identification (TSNI) algorithm that exploits ordinary differential equations are also given [39]. The CKFS algorithm is found to perform significantly better than the TSNI algorithm.

Table 2 

Area under the ROC curve (AUROC) and area under the PR curve (AUPR) for DREAM4 10-gene networks for the five different networks.

Table 3 

Area under the ROC curve (AUROC) and area under the PR curve (AUPR) for DREAM4 100-gene networks for the five different networks.

5.3. IRMA Gene Network

In addition to synthetic data, it is imperative to test the algorithms using real biological data. In this subsection, the performance of the CKFS algorithm is assessed using the in vivo reverse-engineering and modeling assessment (IRMA) network [40]. This network consists of five genes. Galactose activates the gene expression in the network, whereas glucose deactivates it. The cells are grown in the presence of galactose and then switched to glucose to obtain the switch-off data which represents the expressive samples at 21 time points. The switch-on data consists of 16 sample points and is obtained by growing the cells in a glucose medium and then changing to galactose. The system and measurement noise variances for the CKFS are assumed to be identical as in the previous simulations. Figure 5 shows the inferred network, the gold standard, and the network inferred using TSNI. It is observed that the CKFS algorithm succeeds to predict most of the interactions while giving lower false positives.

(a)

(b)

(c)

Figure 5 

The inferred IRMA networks. (a), (b), and (c) Gold standard, inferred network using CKFS, and inferred network using ODE [39, 40]. Black arrows indicate true connections, blue arrows indicate the edges that are correct, but their directions are reversed, and red arrows indicate false positives.

6. Conclusions

This paper presents a novel algorithm for inferring gene regulatory networks from time-series data. Gene regulation is assumed to follow a nonlinear state evolution model. The parameters of the system, which indicate the inhibitory or excitatory relationships between the genes, are estimated using compressed sensing-based Kalman filtering. The sparsity constraint on the parameters is crucial because the genes are known to interact with few other genes only. The use of CKF and the dual estimation of states and parameters renders the algorithm computationally efficient. The performance of CKFS is evaluated for synthetic data for different network sizes as well as varying sample points. ROC curves, Hamming distance, and true positives are used for comparing the accuracy of inferred network with EKF. It is observed that CKFS outperforms the EKF algorithm. In addition, CKFS gives advantages over EKF in terms of smaller run time for large networks. The Cramér-Rao lower bound is also determined for the parameters of the model and compared with the MSE performance of the proposed algorithm. Assessment using DREAM4 10-gene and 100-gene networks and IRMA network data corroborates the superior performance of CKFS. Future research directions include incorporating the estimation of model order in the network inference algorithm.

Acknowledgments

This work was supported by US National Science Foundation (NSF) Grant 0915444 and QNRF-NPRP Grant 09-874-3-235. The material in this paper was presented in part at the IEEE International Workshop on Genomic Signal Processing and Statistics (GENSIPS), San Antonio, TX, USA, December 2011.