Exact reconstruction of gene regulatory networks using compressive sensing

Background We consider the problem of reconstructing a gene regulatory network structure from limited time series gene expression data, without any a priori knowledge of connectivity. We assume that the network is sparse, meaning the connectivity among genes is much less than full connectivity. We develop a method for network reconstruction based on compressive sensing, which takes advantage of the network’s sparseness. Results For the case in which all genes are accessible for measurement, and there is no measurement noise, we show that our method can be used to exactly reconstruct the network. For the more general problem, in which hidden genes exist and all measurements are contaminated by noise, we show that our method leads to reliable reconstruction. In both cases, coherence of the model is used to assess the ability to reconstruct the network and to design new experiments. We demonstrate that it is possible to use the coherence distribution to guide biological experiment design effectively. By collecting a more informative dataset, the proposed method helps reduce the cost of experiments. For each problem, a set of numerical examples is presented. Conclusions The method provides a guarantee on how well the inferred graph structure represents the underlying system, reveals deficiencies in the data and model, and suggests experimental directions to remedy the deficiencies. Electronic supplementary material The online version of this article (doi:10.1186/s12859-014-0400-4) contains supplementary material, which is available to authorized users.


Introduction
Mathematical modeling of biological signaling pathways can provide an intuitive understanding of their behavior [1][2][3]. However, since typically only incomplete knowledge of the network structure exists and the system dynamics is known to be sufficiently complex, the challenge has become to show that the identified networks and corresponding mathematical models are enough to adequately represent the underlying system. In the last years, many data-driven mathematical tools have been developed and applied to reconstruct graph representations of gene regulatory networks (GRNs) from data. These include Bayesian networks, regression, correlation, mutual information and system-based approaches [4][5][6][7][8][9][10]. Also, these approaches either focus on static or on time series data. The latter approach has the advantage of being able to identify dynamic relationships between genes.
However, data-driven reconstruction of the network structure itself remains in general a difficult problem; nonlinearities in the system dynamics and measurement noise make this problem even more challenging. For linear time invariant (LTI) systems, there exist necessary and sufficient conditions for network reconstruction [11]. However, for time-varying or nonlinear systems, there has not been as yet any statistical guarantee on how well the inferred model represents the underlying system [12][13][14][15]. The recent work [16] addresses the problem of data-driven network reconstruction, together with measurement noise and unmodelled nonlinear dynamics, yet this work points out that these complications impose a limit on the reconstruction, and with strong nonlinear terms the method fails. Additionally, identifying whether important nodes in the graph structure are missing, how many are missing, and where these nodes are in the interconnection structure remains a challenging problem.
Moreover, in order to continue to have an impact in systems biology, identification of the graph topology from data should be able to reveal deficiencies in the model and suggest new experimental directions [17]. For example, Steinke et al. [18] presented a Bayesian method for identifying a gene regulatory network from micro-array measurements in perturbation experiments and showed how to use optimal design to minimize the number of measurements.
Since biological regulatory networks are known to be sparse [19][20][21][22][23], meaning that most genes interact with only a small number of genes compared with the total number in the network, many methods [12][13][14][15][24][25][26][27] take advantage of the sparsity. The methods typically use l 1 -norm optimization, which leads to a sparse representation of the network and improves the ability to find the actual network structure. Even though many methods [14,15] show that the reconstruction results are fairly good, the methods cannot guarantee exact recovery. This stems from the fact that the so-called incoherence condition is typically not satisfied for the matrix in a linear measurement model Y = q, where q is the signal to be reconstructed, Y is the measurement, and is known as the sensing matrix. In this paper, we construct Y and from time series protein or gene expression data with candidate basis functions or kinetic features, and q reflects the (unknown) underlying GRN structure to be reconstructed. Roughly, incoherence is a measure of the correlation between columns of the sensing matrix. Since the incoherence condition of the sensing matrix provides a metric of performance, this is one of the motivating factors for the use of compressive sensing (CS) [28] in GRN reconstruction. CS is a signal processing technique for efficiently acquiring and reconstructing a signal by taking advantage of the signal's sparsity and allowing the entire signal to be determined from relatively few measurements under a certain condition, i.e., it requires that the incoherence condition to be satisfied. In the Human Epidermal Growth Factor Receptor2 (HER2) positive breast cancer signaling pathway that we studied in [29,30], time series data sets consist of only 8 time point measurements of 20 protein signals, and we would like to use this limited data to identify a graph structure which could have 20 × 20 or 400 edges.
In this paper, we consider reconstruction of biochemical reaction networks (protein-protein interaction) or GRNs (including mRNA, transcription factors). Biological processes in cells are properly performed by gene regulation, signal transduction and interactions between proteins and thus, the network structure that we consider is an abstraction of the system's biochemical reaction dynamics, describing the manifold ways in which one substance affects all the others to which it is connected. Thus, we are interested in directed graph representations of signaling pathways and propose models of biochemical reaction networks or GRNs as a set of basis functions or features that best explain a set of time series observations such as protein expression or gene expression time series data. We develop a new algorithm for GRN reconstruction based on CS. First, we focus on sparse graph structures using limited time series data with all nodes accessible and no measurement noise. We test the network reconstruction algorithm on a simulated biological pathway in which the structure is known a priori. We explain the limitation of the proposed method's performance when the dataset naturally has high coherence and propose a way to overcome this limitation by designing additional effective experiments. Thus, one big part of the whole study is the method for proposing new experiments. As systems biology matures, experimental study design becomes more important in systems biology rather than simply building computational models of experimental data. Thus, the proposed method can be useful for experimental design. By revealing deficiencies in the data and model, the method can suggest experimental directions to remedy the deficiencies. By doing this, we can minimize the number of experiments and reduce the cost of experiments.
Next, the proposed algorithm is extended to a more general problem: we consider partially corrupted data with data inconsistencies due to model mismatch, and measurement noise affecting all the data. Typically, data inconsistencies may result from missing nodes in the model; or in some cases arbitrary data corruption may result from human errors such as mislabeling or the improper use of markers or antibodies. The question is whether one can still recover the graph structure reliably under these conditions. Inspired by a robust error correction method [31], the exact recovery of the graph structure can be guaranteed under suitable conditions on the node dynamics, provided that hidden nodes can affect relatively few nodes in the graph structure. Also, a set of numerical examples is provided to demonstrate the method, including some from an RPPA (Reverse Phase Protein Array) dataset [32] collected from HER2 positive breast cancer cell lines.
In this paper, the main contributions are the following.
• The CS framework uses the coherence of the sensing matrix as a performance index, which allows us to assess and optimize mathematically network reconstruction.
• Coherence also provides a guideline for optimizing experiment design for network reconstruction. • By utilizing an error correction method in conjunction with the CS framework, network reconstruction may be performed even when there are hidden nodes and measurement noise.

Overview: compressive sensing
Consider measurements Y ∈ R m of a signal q ∈ R n : where ∈ R m×n is called the sensing matrix.
One key question [33] is how many measurements m are needed to exactly recover the original signal q from : • If m > n and is a full rank matrix, then the problem is overdetermined. If m = n and is a full rank matrix, the problem is determined and may be solved uniquely for q. • If m < n, the problem is underdetermined even if has full rank. We can restrict q ∈ R n to the subspace which satisfies Y = q. However, q cannot be determined uniquely.
For the underdetermined case, the least squares solution q * = arg min q q l 2 = * ( * ) −1 Y is typically used as the "best guess" in many applications. However, if q is known to be sparse, meaning that many of its components are zero, one might expect that fewer than n measurements are needed to recover it. It is thus of interest to obtain a good estimator for underdetermined problems in which q is assumed to be s-sparse, meaning that s of the elements of q are nonzero. In principle, the theory of compressive sensing (CS) asserts that the number of measurements needed to recover q ∈ R n is a number proportional to the compressed size of the signal (s), rather than the uncompressed size (n) [34]. To be able to recover q, CS relies on two properties: sparsity, which pertains to the signals of interest, and incoherence, which pertains to the sensing matrix. Proposition 1. [33]. Suppose that any 2s columns of an m × n matrix are linearly independent (this is a reasonable assumption if m ≥ 2s). Then, any s-sparse signal q ∈ R n can be reconstructed uniquely from q.
The proof [33] of the above proposition also shows how to reconstruct an s-sparse signal q ∈ R n from the measurement Y = q where q is the unique sparsest solution to Y = q: q * = arg min q q l 0 subject to Y = q (2) and q l 0 := n i=1 I(q i = 0) is the cardinality of q. However, the l 0 -minimization is computationally intractable (NP-hard in general). Recent breakthroughs enable approximating l 0 -optimization by using l 1minimization which is a convex optimization problem and can be solved in a simple but effective way by linear programming: The l 1 -minimization in Equation (3) requires mild oversampling, more specifically, m ≥ c · μ( , ) 2 · s log n for some positive constant c where can be decomposed as a product of a sparsity basis , and an orthogonal measurement system [28]. We have m measurement in the domain (sensing modality) under bases (signal model) and μ represents the coherence defined as follows: Definition 1. [28]. The coherence of a matrix is given by where j and k denote columns of . In plain English, the coherence measures that the largest correlation between any two columns in .
Several theoretical results [35] ensure that the l 1minimization guarantees exact recovery whenever the sensing matrix is sufficiently incoherent. For example, we can say that is incoherent if μ is small. Coherence is a key property in the compressed sensing framework because if two columns are closely correlated, it would be very difficult to distinguish whether the influence from the components in the sparse signal comes from one or the other (recall that the measurement Y(= q) is a linear combination of each columns of the sensing matrix with the components in q as coefficients). Also, numerical experiments suggest that in practice, most s-sparse signals are in fact recovered exactly once m ≥ 4s [28]. Here, "exact recovery" means that we find the sparsest solution (q) such that Y = q.
Therefore, if the sensing matrix satisfies the incoherence condition (m ≥ c · μ( , ) 2 · s log n, or in practice, m ≥ 4s), a sufficiently sparse signal q can be exactly recovered from the limited dataset without any prior knowledge of the number of nonzero elements, their locations, and their values. On the other hand, if the condition is not satisfied, exact recovery cannot be guaranteed [17,34]. However, it is possible to use the property of coherence to guide biological experiment design, basically to collect a more informative dataset. As we will discuss in this paper, this can be done by inhibiting or stimulating certain genes to manipulate the gene expression.

CS can help reconstruct GRNs
In graph theory, a digraph can be represented by G = (V , E) where V and E represent nodes and edges respectively. For GRNs, each node represents a gene and each edge represents an influence map which models how genes affect each other. For example, the interactions could be how genes inhibit or stimulate each other. Since the connectivities of GRNs are typically unknown, often the best we can do is to select a set of possible candidate functions encoding possible unknown connectivities between genes.
In this section, we formulate a data-driven network identification problem into CS framework: first, we define a dynamical model of gene regulatory network. Then, we encode system dynamics into the sensing matrix ( ) and denote unknown connectivities between genes by q, a signal to be recovered.

A dynamical model of gene regulatory networks (GRNs)
We consider a dynamical system described by: where x ∈ R n denotes the concentrations of the ratelimiting species which can be measured in experiments; ∈ R n is a vector whose elements are the change in concentrations of the n species over time; f(·) : R n → R n represents biochemical reactions, such as those governed by mass action kinetics, Michaelis-Menten, or Hill kinetics. Thus, f(·) can include functions of known form such as product of monomials, monotonically increasing or decreasing Hill functions, simple linear terms, and constant terms [24]. u ∈ R n denotes the control input which could represent inhibitions and stimulations; For example, if we consider protein-protein interaction, u represents druginduced perturbation such as inhibition and stimulation. For a GRN, we consider the use of siRNA as a control input to the system for gene knockdown; g(·) : R n h → R n represents influence from hidden nodes x h ∈ R n h , which cannot be measured in experiments; n h is the number of hidden nodes and unknown; and w ∈ R n represents energy-bounded process noise or measurement noise. Here x,ẋ and u is assumed to be known where x can be measured in experiments andẋ may not be measured directly in experiments but we could calculate these quantities by interpolating x and using numerical derivatives. In [36], the authors pointed out that although data on time derivative can be difficult to obtain especially in the presence of noise, it is possible to estimate the gene expressions relatively accurately by repeating measurement with careful instrumentation and statistics [37].
Since we do not know whether important nodes in the gene regulatory network are missing, how many missing nodes there are, and how they affect system dynamics (i.e., x h , n h and g(x h )), we denote a vector of (unknown) influence from hidden nodes' dynamics by v( g(x h )); Also, without loss of generality, since we know u, we define y as follows:

Formulating a dynamical system as a GRN
The nonlinear function f(x) can be decomposed into a linear sum of scalar basis functions, f b,i (x) ∈ R, where we select the set of possible candidate basis functions: where N is the number of possible candidate basis functions and q ij are unknown parameters which reflect underlying structure, i.e., influence of f b,i (x) on the j-th state (ẋ j ). Typically, we may choose a larger set than necessary, and allow the CS method to indicate the importance of each function, as we shall describe. Thus the Equation (5) can be written as follows: where is the vector field which includes possible candidate basis functions. In this way, any biochemical reactions can be represented by a linear map S q and a function F b (x) where S q encodes the underlying graph structure and F b (x) includes all possible candidate functions that could be included in the biochemical reactions.
In practice, we can construct F b (x) by selecting the most commonly used candidate basis functions to model GRNs or protein-protein interaction, for example, all monomials, binomials, other combinations and Hill function.

Example 1. Consider the simple nonlinear ordinary differential equations (ODEs):
where x i denotes the concentration of the i-th species,ẋ i is the change in concentration of the i-th species, the γ i denotes protein decay rate, the k i denotes the maximum promoter/inhibitor strength. Here, there is no input (u = 0), no hidden node (v = 0) and no process noise (w = 0). Also, n act represents positively cooperative binding (activation) and n ihb represents negative cooperative binding (inhibition). The set of ODEs corresponds to a topology where gene 1 is activated by gene 2, gene 2 is inhibited by gene 3, and gene 3 is inhibited by gene 1 as shown in Figure 1(A). We can write Equation (8) as follows: We can also consider a version of Equation (8) in which there exists a hidden node (x h ) affecting (x 1 ) as shown in Figure 1(B), as well as process noise.

Formulating GRN into the CS framework
Suppose the time series data are sampled from a real experimental system at discrete time points t k . By taking the transpose of both sides of Equation (7), we obtain where q = S q = [q 1 q 2 . . . q n ] ∈ R N×n . Assuming that M successive data points are sampled, then define: This leads to n independent equations: where ∈ R M represents the M successive data points, ∈ R M×N consists of N possible candidate bases which are functions of given time series data x and q i ∈ R N represents the unknown influence map corresponding to the i-th species. Since a biochemical reaction network is typically sparse, as a consequence, q i is sparse and we have N M for because we assume the limited time series data and may choose a larger set of basis functions than necessary.
Although we formulate n independent linear regression problems in Equation (11), we consider n independent  (11) together by stacking y i as follows:

equations in Equation
Now, Equation (12) is in the form of the CS formulation in Equation (1) where Y ∈ R n·M is the measurement, ∈ R n·M×n·N is the sensing matrix which consists of basis functions for the given time series data, v [v 1 ; v 2 ; . . . ; v n ] ∈ R n·M represents possibly large corruption from hidden nodes and w [w 1 ; w 2 ; . . . ; w n ] ∈ R n·M represents process or measurement noise.
In this paper, we make the following assumptions: Assumption 1. We consider two cases: • (ideal case) We assume that we can measure all states x (i.e., there is no hidden node (x h ), thus v = 0) and there is no measurement noise (w = 0). Also, there are enough columns in in Equation (11) to represent the underlying system in Equation (6).
) is assumed to be sparse, and process noise (w = 0). Here, we can also consider the case where the columns (also known as the dictionaries in the CS literature [28]) of may not be able to represent the underlying system. In this case, we consider the influence from the missing dictionaries as v.
We first consider the ideal case for simplicity in explaining the main results, and then extend the proposed method to the more general case.

Formulating GRN identification problem into CS
Many existing algorithms [16,24] consider the n independent linear regression problems in Equation (11) separately. Since the columns of matrix are composed of time series data as in Equation (10), it is difficult to a priori guarantee low correlation and sometimes even suffers rank deficiency. Also, Pan et al. [24] pointed out that correlation between the columns of is usually high (μ( ) is close to 1).
Intuitively, if two columns of the sensing matrix are highly correlated, it is hard to distinguish the corresponding components in the sparse signal q (note that the measurement Y is a linear combination of each column of the sensing matrix with the components in q as coefficients). In order to deal with high coherence in Equation (11), many method combines CS with different techniques such as Bayesian formulation [24], Kalman filter [25], and Granger causality [26]. Also, it is a well known problem in the lasso formulation of network inference: if there is high coherence in the sensing matrix, one can use an elastic net which combines the l 1 and l 2 norms. Although each reconstruction result [24][25][26] might be the optimal solution in the sense of its formulation (i.e., maximum likelihood), the identified graph may not represent the underlying GRN. In other words, if the data set is not informative enough to fully explore the underlying system, while the identified graph structure based on the given data set may be an optimal solution of the particular optimization problem, it may not represent the true system.
In this paper, our goal is to get the smallest data informative enough to recover the underlying graph structure exactly. Since the proposed method maintains the CS framework by reducing coherence of the sensing matrix, the method is fundamentally different than any other methods which make use of different techniques in conjunction with the l 1 optimization [24][25][26] which leads to a sparse representation of the network. Hence, we use all the properties of CS in order to access the ability to exactly reconstruct the underlying graph structure, reveal deficiencies in the data and model, and design new experiments to remedy the deficiencies if necessary.
While maintaining the CS framework, in order to deal with high coherence, we formulate Equation (12) instead of Equation (11) since we have strongly uncorrelated columns in . In other words, since has many independent columns, we have more degrees of freedom to reduce coherence of . We will show that by using a transformation, the components of the sensing matrix can be made more uniformly distributed so that we could reduce coherence.
Moreover, since each q i has different degrees of sparsity in general, if we consider n independent equations in Equation (11), should satisfy the incoherence condition M ≥ 2 max i q i 0 stated in Proposition 1.
On the other hand, only needs to satisfy the condition M · n ≥ 2 · n i=1 q i l 0 . Since the averaged sparsity = n i=1 q i l 0 n ≤ max i q i l 0 is smaller than max i q i l 0 , we can reduce the required number of samples (M). Also, in case of rank deficiency, we can simply remove the corresponding rows (note that the rank deficiency is more likely caused by the row since N M).

Optimal design of sensing matrix Reducing coherence by transformation
In order to reduce coherence, first we rearrange with respect to a spatial information, and then we consider the transformation in order to reduce coherence (see Method: Rearranging the sensing matrix for mathematical details): where ⊗ represents the Kronecker product a ,¯ ∈ R M·n×N·n represents the rearranged sensing matrix multiplied by the transformation s j , andz(j),v andw are defined in Equation (20). We want to find s has full row rank (14) where ψ s ij ∈ R n represents the j-th column of s i . In this paper, we propose a heuristic approach and a novel way to find s i by solving the optimization iteratively to reduce coherence (see Method: Section 'Randomly chosen matrix ' and 'Finding the optimal transformation' for details).

Example 2.
(reducing coherence for a linear system) consider a simple linear system (ẋ = Ax) where n = 5, N = 5, M = 4, s = 10 (note that n · M ≥ 2s), and the elements of A are randomly chosen such that there is no isolated node. Figure 2 shows the sensing matrix for both (top left) and the transformed sensing matrix (top right, denoted by ). By reducing coherence, the components of the transformed sensing matrix (top right) are more uniformly distributed and the coherence is reduced by up to 0.6 although μ( ) is close to 1 (bottom left) in Figure 2(A). Also, Figure 2(B) shows the result of the inferred graph structure based on given time series data without any a priori information where the x-axis represents indices of the influence map (i.e., the 1 st , 2 nd , . . . , n th rows of influence map S q ; note that for a linear system, S q = A).
Here, there are 5 states (n = 5) in a linear system so the influence map A has 25 elements. Although L1 and L2 b norm minimizations fail to recover the exact signal, CS in Equation (17) (11) without reducing coherence and (L2) solves Equation (12) with l 2 -regularization. Example 3. (reducing coherence for a nonlinear system, n = 3, N = 9, M = 5, s = 6) consider simple nonlinear ODEs as follows:

(see Recovery of gene regulatory networks for details) recovers the exact signal (bottom right) by reducing coherence of the sensing matrix. Note that (L1) solves the n independent equations in Equation
2, n act = 4 (activation Hill coefficient) and n ihb = 4 (inhibition Hill coefficient). The set of ODEs corresponds to a topology (x 2 → x 1 → x 3 x 2 ) as shown in Figure 3(A). Figure 3(B) shows that we can reduce the coherence by up to 0.8 and (C) shows that only CS recovers the exact graph structure, and L2 regularization does not encourage sparsity but distributes the coefficients to be more similar to each other.

Designing effective experiments
Consider the case where the sensing matrix is not incoherent. If the coherence condition (M · n > cμ(¯ ) 2 · s log(N · n)) is not satisfied, exact recovery cannot be guaranteed [34]. We use the transformation ( s i or P = P P in Equation (24); see Method section 'Reducing coherence by transformation') in order to reduce the coherence but obviously, sometimes we might have inherent limits to how much the coherence can be reduced. There are possible reasons: • Since we solve the relaxed problem in Equation (24) iteratively, P might be sub-optimal. • If the time series data of two different gene expressions, x i and x j are highly correlated, it might be difficult to reduce coherence. In this case, we need to design a new experiment to remedy deficiencies in the data.
As we mentioned, the incoherence of the sensing matrix can be used not only as a good metric to guarantee exact recovery but also as a guideline for designing new experiments. For example, from the coherence distribution, we can identify which columns of the sensing matrix have high coherence, i.e., f b,i (x) and f b,j (x) in Equation (6). Intuitively, in order to reduce ambiguities from the highly   To show the effectiveness of the new experiment, we design two different experiment sets and compare the reconstruction results with each other; for Exp#2, we perturb x 3 and for Exp#3, we perturb x 2 . As we mentioned earlier, intuitively, we expect that Exp#3 to be a more informative experiment to identify the graph structure since we would like to reduce the coherence between x 2 and x 4 . As we expected, in Exp#3, we can reduce the coherence more than that of Exp#2 as shown in Figure 4, and recover the exact show that x 2 and x 4 cause high coherence. Thus, in order to reduce this high coherence, we should perturb either x 2 or x 4 . We design two different experiment sets; for Exp#2, we perturb x 3 and for Exp#3, we perturb x 2 . As we expected, Exp#3 is a more effective experiment to reduce the coherence between x 2 and x 4 and we recover the exact graph structure as shown in Figure 5. On the other hand, the coherence of Exp#2 remains almost the same as that of Exp#1 and we fail to recover the exact graph structure in Figure 5.
graph structure as shown in Figure 5. On the other hand, the coherence of Exp#2 remains almost the same as that of Exp#1 and we fail to recover the exact graph structure. This numerical example illustrates that by using the property of coherence, we can guide biological experiment design more effectively.

Example 5.
(limitation of reducing coherence by , nonlinear dynamics, n = 3, N = 9, M = 5, s = 9) consider the following set of ODEs: where (γ 1 , γ 2 , γ 3 ) = (−0.25, −0.23, −0.26), k +12 = 1.2, k +13 = 1.25, k −21 = 2.8, k −23 = 2.1, k −31 = 2.7, k −32 = 1.8, n act = 4 (activation Hill coefficient) and n ihb = 4 (inhibition Hill coefficient). The corresponding topology is shown in Figure 6(B) (note that we intentionally choose the symmetric structure and similar parameters). The reconstruction error using Exp#1 data is shown in Figure 6(A) (left, bottom) and the reconstruction error illustrates difficulties of resolving ambiguities from x 2 and x 3 . This can be captured by the coherence distribution of the sensing matrix based on the Exp#1 dataset; the correlation between the columns corresponding to x 1 and x 2 is close to the correlation between the columns corresponding to x 1 and x 3 . Based on the coherence distribution, we design two trials; for Exp#2, we perturb x 1 and for Exp#3, we perturb x 2 . As we expected, Exp#2 is not an effective experiment in terms of information. On the other hand, by using Exp#3, we can reduce both the maximum coherence and the averaged coherence, and reconstruct the exact graph structure as shown in Figure 6(A).
Both Examples 4 and 5 illustrate that if the transformed sensing matrix is not incoherent enough to guarantee exact recovery, we can design a new experiment based on the distribution of coherence. Also, we show that the proposed experiment can help to reduce coherence more and thus reconstruct the exact graph structure. For a fair comparison, we use the same number of time points as M here. However, in practice, we can also stack all the experimental data sets together if we assume that the linear map S q where the subscript Z i ,¯ i ,v i ,w i represents the i-th experiment. As the number of measurements increase (M), one may be able to reduce the coherence. However, one can reduce the coherence only if the additional measurements provide us more useful information. As a trivial example, one could stack exactly the same data on top of the first, and increase M to 2M, however the coherence is exactly the same as that of the original dataset.

Recovery of gene regulatory networks
In this section, we present reconstruction of the exact graph structure and show how the condition for exact recovery will be used. First, we consider the ideal case where there are no hidden nodes and no measurement noise. Second, we extend the ideal case to the more general case.

Reconstructing gene regulatory networks (ideal case)
In Equation (13), q represents the s-sparse network structure which we want to reconstruct from the  min q l 1 s.t.Z =˜ q (CS) (17) whereZ PZ,˜ P¯ and P is the optimal transformation in Equation (24).

Proposition 2.
If the sensing matrix˜ constructed from time series data, multiplied by the optimal transformation, P, has 2s linearly independent columns, then any s-sparse network structure q can be reconstructed uniquely from Z =˜ q.
The requirement of 2s linearly independent columns in Proposition 2 may be translated to an incoherence condition on the sensing matrix. That is, if the unknown s-sparse signal q is reasonably sparse, it is possible to recover q under the incoherence condition on the sensing matrix. Although the sensing matrix consists of redundant dictionaries, the coherence of the sensing matrix can be reduced. In a heuristic way, we multiply the redundant dictionaries k,: by a randomly chosen matrix s k at each time step k and iterate this step until the coherence is decreased. Or, we can find the optimal transformation P(or P) in Equation (24) to reduce the coherence. In the previous numerical examples, we illustrated that the coherence of the sensing matrix is decreased by transformation and showed the exact reconstruction of graph structure.

Example 6. (Statistics)
Here, we compare the success rate of the proposed method with other methods such as L1 and L2. Figure 7 shows statistics of 50 trials for a simple linear case (for each trial, we randomly generate the influence map). Here, we count the number of successes of each method when any of the methods recover the exact structure. By reducing coherence, we can improve the success rate as shown in Figure 7(A). Also, L2-regularization does not encourage sparsity but distributes the coefficients to be more similar to each other as shown in Figure 7(B).

Graph reconstruction with hidden nodes
The main contributions of the proposed method in the previous section is the conversion of the problem of inferring graph structure into the CS framework. Then, we demonstrate that one could recover sparse graph structures from only a few measurements. However, for practical use, the proposed method needs to be able to deal with both sparsely corrupted signals (v = g(x h )) and measurement noise (w) in Equation (5).
In general, the assumption of accessibility or observability of all nodes [17] is not satisfied. Thus, we focus on the case in which the hidden node affects observable nodes directly as shown in Figure 1(B). Also, without loss of generality, the hidden node dynamics could be any arbitrary dynamic model. Or, even if there is no hidden node, a small portion of the biological experiment dataset could be in practice contaminated by large error resulting from, for example, mislabeling, or improper use of markers or antibodies. Moreover, all biological datasets are contaminated by at least a small amount of noise from measurement devices. Therefore, the proposed method should be robust. We note this goes beyond the results in [11] due to the consideration of hidden node dynamics with measurement noise.
Here, the question is whether it is still possible to reconstruct the graph structure reliably when measurements are corrupted. Since hidden nodes and measurement noise are considered, the number of time points is assumed to be greater than that of the previous case [17] (i.e., no corruption and no measurement noise). Thus, the number of rows of the sensing matrix is assumed to be greater than the number of columns. If the number of the time points M(< N) is limited, then we can stackz(j) with different s k or including different dataset as described in Equation (16).
In CS literature [31], two decoding strategies for recovering the signal from a corrupted measurement are introduced, where the corruption includes both a possible sparse vector of large errors and a vector of small error affecting all the entries. It is shown that two decoding schemes allow the recovery of the signal with nearly the same accuracy as if no sparse large errors occurred. Our contribution is converting the problem of inferring the graph structure with hidden nodes into the highly robust error correction method framework [31] (see Method: Section 'Two-step refinements') and showing how this can improve the reliability of reconstruction.

Example 7.
(arbitrary corruption with no measurement noise) Consider nonlinear ODEs as follows: where h act , h ihb represents Hill functions for activation and inhibition respectively, and v i represents arbitrary corruption shown in Figure 8 (at each time step, we choose card(v(j)) = 1). The magnitude of v i are about 50% of the magnitude ofẋ. Since we consider arbitrary corruption, we need more time points (M > N). By using two-step refinements, first we estimate sparse large corruption v as shown in Figure 8(B) (top) and then, we reconstruct q (bottom).
In practice, a specific node is corrupted by a hidden node and a small portion of the dataset can be largely corrupted by human error. Also, since we choose the set of possible candidate basis functions of the sensing matrix in Equation (6), the columns of the sensing matrix may not be able to represent the underlying system (i.e., missing dictionaries). Then, we can consider the influence from these missing dictionaries as v. (18) with different parameters and consider sparse large corruption v and small magnitude noise w (1% of the magnitude ofẋ). Figure 9 shows the time series data and reconstruction result.

Geometric view
In Equation (17), since we assume all nodes are accessible and perfect measurement (meaning that there is no hidden node, v = 0 and no measurement noise, w = 0), we can solve Equation (27) directly without filtering out the unmodelled dynamics in Equation (26) (i.e., Z =¯ q in Equation (27)). If there exist hidden nodes or measurement noise, we can still provide an unambiguous indication of the existence of these corruption (e = 0).
The intuition is that Z(=¯ q +¯ e) can be decomposed as the superposition of an arbitrary element in V (= Z − QQ * e) and of an element in V ⊥ (= QQ * e) as shown in Figure 10. In other words, Z can be decomposed as the superposition of modelled dynamics and anomalies caused by hidden node or unmodelled dynamics. This geometric view enables us to understand how we could reveal deficiencies in our model: •Z = Q * Z = 0: there is no hidden node =⇒ Z ∈ R(¯ ) •Z = 0: Z cannot be represented by¯ q so there might be hidden nodes or our dictionaries in the sensing matrix¯ are not sufficient to represent Z (revealing deficiencies in our model or dictionaries).

HER2 Overexpressed breast cancer
We apply the proposed algorithm to study a breast cancer signaling pathway by reconstructing the graph structure using an RPPA dataset [32] as shown in Figure 11 (see also Additional file 1: Figure S1, S2, and S3 for details: each figure presents the RPPA dataset and the result of graph reconstruction compared with L1, L2optimization). Here, we choose small networks which are composed of 3 nodes and known to be sparsely connected, i.e., PI3K → PDK → Akt and PDK → Akt → mTOR in order to satisfy our assumption such that the influence on observable nodes from a hidden node should be sparse (i.e., v is sparse). The graph structures identified by the proposed method are consistent with the current understanding of the networks, whereas those found using L1-and L2-optimizations fail to reconstruct the known structure as shown in Figure 12. Also, an abstract model of the breast cancer signal pathway proposed by M. Moasser [38] is considered, as shown in Additional file 1: Figure S3(B) where PHLPP isoforms are a pair of protein phosphatases, PHLPP1 and PHLPP2, which are important regulators of Akt serinethreonine kinases (Akt1, Akt2, Akt3) and conventional protein kinase C (PKC) isoforms. PHLPP may act as a tumor suppressor in several types of cancer due to its ability to block growth factor-induced signaling in cancer cells [39]. PHLPP dephosphorylates Ser473 (the hydrophobic motif ) in Akt, thus partially inactivating the kinase [40]. Unfortunately, in our RPPA dataset, we do not have PHLPP so we simply consider three nodes (Akt pT308 , Akt pS473 and mTOR). Figure 11(right) shows the result of the proposed method using the RPPA dataset. The reconstructed graph structure matches up to the known structure (Additional file 1: Figure S3(B)). Specifically, our result can capture the partial inactivating characteristics of PHLPP (i.e., mTOR(→ PHLPP) Akt pS473 ).

Discussion
Many network inference algorithms use l 1 -norm optimization to find the actual network structure, but they cannot guarantee exact recovery. In this paper, we propose a novel approach to reconstruct GRNs based on compressive sensing with a dynamic system model. We demonstrate that the incoherence condition is essential for exact recovery. This is not properly considered in inference methods based on l 1 -norm optimization. Also, we illustrate how the incoherence condition can be used to design new experiments effectively. Finally, we consider a more general setting, in which hidden genes exist and all measurements are contaminated by noise, and we show that the proposed method leads to reliable reconstruction.
We compare all the results with l 1 -and l 2 -norm optimization, since l 1 -norm optimization is widely used in network inference and many inference methods combine the l 1 -and l 2 -norms (for example, elastic net). We highlight the improvement by the proposed method compared with these commonly used methods. Our use of the method to propose new experiments, and its extension Figure 10 Geometric view of two-step refinement. A geometric view of two consecutive l 1 -norm optimizations. By using two-step refinements, first we estimate sparse large corruption v and then, we reconstruct q. For the ideal situation with no unobservable species (genes) or noise, we directly reconstruct q without estimating v.
to a more general setting goes beyond the results in any network inference methods.
From the lessons learned from the numerical study using both synthetic datasets and real RPPA datasets of HER2 overexpressed breast cancer signaling pathways, the proposed method can be applied to sub-networks in order to satisfy our assumption about the number of hidden nodes relative to the network size, and has great potential for reconstructing GRN by designing effective experiments. Note that the proposed method is not limited to the network size itself but rather by this assumption and the incoherence condition. In the ideal case with no unobservable species (genes) and no noise, the method can be applied to entire networks if the incoherence condition is satisfied (i.e., n · M ≥ 2s, practically n · M ≥ 4s) where n is the number of nodes, M is the number of time points and s represents the sparsity of the network structure. For example, suppose we only have one time point measurement (M = 1) and then, the proposed method cannot be applicable for any network (n > 2 where n is the number of nodes) since it is impossible to satisfy the incoherence condition. As an example, if we consider n nodes which are simply connected (i.e., s = n − 1), then n · M < 2s. With a reasonable number of time points (M), the incoherence condition could be satisfied. Similarly, in the practical case with unobservable species (genes) and measurement noise, in order to guarantee the exact recovery of the graph structure, we should satisfy Assumption 1 for a given network (hidden node only affects relatively few nodes in the given network) and the incoherence condition. Suppose we choose a small sub-network where a hidden node or unobservable node affects all nodes in the chosen sub-network: then, it is impossible to infer the exact structure of the sub-network. However, if a hidden node only affects relatively few nodes in the given network, we can still reconstruct the exact structure of the sub-network by inferring the hidden node influence first. Therefore, the proposed method is not limited to the network size itself but rather by Assumption 1 and the incoherence condition. In practice, we can work on a sub-network at a time and integrate the identified sub-networks.
Finally, we evaluate performance with respect to the size of the problem under typical parameters (the number of genes, the number of kinetic features and connectivity) in (Additional file 1: Table S1) which provides us a brief guideline of application, e.g., for an n-genes network with certain kinetic features, how many time points/how much resolution are appropriate for the reconstruction.
In this paper, we do not consider any a priori knowledge of connectivity. However, since we may have partial information of connectivity, we can also use prior knowledge of connectivity in the form of known nonzero elements in q allowing fewer experiments to be performed. Also, here we only consider protein expression or gene expression levels but we can also consider many other Figure 11 HER2+ overexpressed breast cancer. CS reconstruction result using Reverse Phase Protein Array data (SKBR3 cell line, Serum [32]). See Additional file 1: Figure S1, S2 and S3 for further details (red: activation, blue: inhibition). Here, we choose small networks which are composed of 3 nodes and known to be sparsely connected in order to satisfy our assumption such that the influence on observable nodes from a hidden node should be sparse. The identified graph structures are consistent with the current understanding of the networks. types of useful information that can be incorporated into the network reconstruction process, such as sequence motifs and direct binding measurements (e.g., ChIP-chip and ChIP-seq). Since it has been shown that ChIP-seq signals of Histone modification are more correlated with transcription factor motifs at promoter sites in comparison to RNA level, time series Histone modification ChIP-seq could provide a more reliable inference of GRNs in comparison to method based on expression level. We are interested in continuing this research direction for large-scale networks.

Conclusion
We proposed a method for reconstructing sparse graph structures based on time series gene expression data without any a priori information. We demonstrated that the proposed method can reconstruct graph structure reliably. Also, we illustrated that coherence in the sensing matrix can be used as a guideline for designing effective experiments.
Second, the proposed method is extended to the cases in which dynamics is corrupted by hidden nodes and the measurement is corrupted by human error in addition to the measurement noise. Using a two-step refinement procedure, we demonstrate good performance for the reconstruction of graph structure. A set of numerical examples is implemented to illustrate the method and its performance. Also, a biological example of HER2 overexpressed breast cancer using an RPPA dataset is studied. We are currently applying our method to recover the HER2 signaling pathway, where a significant part of the network is currently unknown.

Rearranging the sensing matrix
Define z(j) as a vector of each component of y i at the j-th time point: where j,: represents the j-th row of . Consider s j ∈ R n×n and multiply Equation (19)

Randomly chosen matrix
The optimization problem in Equation (14) is not trivial because of the constraint. One simple and heuristic approach is that we select s i by (normalized) randomly chosen matrix with independent identically distributed (i.i.d.) random variable where¯ has full row rank, calculate μ(¯ ) and run this with several times to reduce coherence. Since the randomly chosen matrix, s i , spreads out the component of i,: uniformly, we can reduce coherence.

Finding the optimal transformation
This heuristic approach may not be enough to reduce the coherence. Consider a nonsingular matrix P ∈ R M·n×M·n and˜ = P¯ where¯ is constructed by heuristic way, i.e., randomly chosen matrix.
where P P P = P ∈ S M·n is positive definite and i denote the i-th column of˜ . Note that if P P = I, We can also combine · ∞ and · 1 to reduce the coherence. Note that for · ∞ , we minimize the maxim coherence and for · 1 , we minimize the sum of the all possible combinations of the columns of¯ , i.e.,μ ij . Thus, if certain bases are highly correlated, P or P makes the components of the sensing matrix spread out enough to differentiate the influences from those bases. Since we ignore the denominator, the optimal solution of Equation (24) may be suboptimal. Thus, we can also combine heuristic way and Equation (24) iteratively to reduce coherence of the sensing matrix in practice.

Two-step refinements Sparse large corruption with no measurement noise
Recall Equation (13) where¯ diag s 1 , s 2 , . . . , s M and v(j) represents sparse large corruption at the j-th time point, that could result from the existence of hidden nodes. We assume that the influence from hidden nodes is sparse and unknown (i.e., v(j) is assumed to be sparse). In other words, hidden nodes can affect only a few nodes' dynamics (intuitively, if hidden nodes affect all nodes, there is no way to reconstruct the graph structure). Then, we consider the reconstruction of graph structure q from the corrupted signal Z.
Since we assume the number of rows of the sensing matrix¯ is greater than the number of columns (M · n > N · n), we consider Q * which annihilates the sensing matrix¯ on the left (Q * ¯ = 0) where Q * is any (M · n − N · n) × M · n matrix whose kernel is the range of in R M·n (rank(Q * ) + nullity(Q * ) = M · n): Then, the following two-step optimization problem enables us to compute q: • Filter unmodelled dynamics out from the measurement: • Reconstruct q : min q l 1 s.t. Z −¯ v =¯ q orZ = P(Z−¯ v) = P¯ q =˜ q (P in Equation (23)) If we could somehow get an accurate estimatev from Equation (26), Equation (27) represents the problem of reconstructing the graph structure q. The intuition is that Z(=¯ q+¯ v) can be decomposed as the superposition of modelled dynamics and anomalies caused by hidden node or unmodelled dynamics.
The two step convex optimization problems in Equation (26) and Equation (27) are l 1 -norm optimization problems in CS. Thus, if the sensing matrix Q * ¯ and¯ (or P¯ ) satisfy the incoherence condition, signals v and q can be recovered exactly [17,34]. Here, there are many possible choices of Q * but we have to choose Q * to satisfy the incoherence condition for the exact recovery of v [17]. To choose such a Q * , we observe that¯ can be denoted as follows using Singular Value Decomposition (SVD): where U 1 ∈ R M·n×r , U 2 ∈ R M·n×(M·n−r) , ∈ R r×r , V 1 ∈ R N·n×r , V 2 ∈ R N·n×(N·n−r) and r is the rank of¯ . Suppose we choose Q * such that Q * = U 2 . Then: can be used as a tuning matrix for satisfying the incoherence condition. A geometric view in Results and discussion: Section 'Geometric view'. enables us to understand how we could reveal deficiencies in our model.

Sparse large corruption with measurement noise
While considering influence from hidden nodes is interesting, it still may not be realistic to assume that except for hidden nodes, one is able to measure the node dynamics with infinite precision. A better model would assume that there is measurement noise. Consider the problem of recovering the graph structure q from the vector Z which is corrupted by measurement noisew in Equation (13): . . .
where e(j) = v(j) + w(j), w is Gaussian noise N (0, σ ) assumed to be bounded w l 2 ≤ . In general, we can consider any corruption decomposed into sparse large error v and small magnitude error w [31]. Then, modified two-step refinements can be applied as follows: where the parameters 1 , 2 above depend on the magnitude of the small errors , which can be determined as in [31].

Endnotes
a If A is an m × n matrix and B is a p × q matrix, then the Kronecker product product A ⊗ B is the mp × nq block matrix: