An Extended Constrained Total Least-Squares Method for the

Article pubs.acs.org/IECR

An Extended Constrained Total Least-Squares Method for the Identification of Genetic Networks from Noisy Measurements Ugur Guner,‡,§ Hong Jang,†,§ Matthew J. Realff,‡ and Jay H. Lee*,† †

Department of Chemical and Biomolecular Engineering, Korea Advanced Institute of Science and Technology (KAIST), 291 Daehak-ro, Yuseong-gu, Daejeon 305-701, Republic of Korea ‡ School of Chemical and Biomolecular Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332, United States S Supporting Information *

ABSTRACT: We address the system identification problem of genetic networks using noisy and correlated time series data of gene expression level measurements. Least-squares (LS) is a commonly used method for the parameter estimation in the network reconstruction problems. The LS algorithm implicitly assumes that the measurement noise is confined only to the dependent variables. However, a discrete time model for the genetic network systems will lead to serially correlated noise terms that appear in both the dependent and independent variables. A constrained total least-squares algorithm (CTLS) used in signal and image processing applications showed significant improvements in such an estimation problem over the LS and total least-squares (TLS) methods. In this paper, we propose an extended CTLS algorithm that estimates parameters for all the dependent variables simultaneously, instead of estimating them separately for each dependent variable, as in the original CTLS algorithm. In addition, the CTLS algorithm is further generalized to assign weights to the error terms according to the variances or covariances of the measurement noise. We demonstrate its improved performance over the original CTLS method, as well as the commonly used LS and TLS methods on a widely adopted artificial genetic network example, under a variety of noise conditions.

1. INTRODUCTION It is becoming increasingly important to reconstruct biological networks in order to understand complex disease mechanisms, identify novel and safer protein targets for therapies, and design efficient drugs. Systems biology has emerged as a discipline to uncover the biological networks through measured data at the genetic activity level.1 Within it, computational methods for identifying the networks are important and have been growing in parallel with the availability of genomic data of various quality and quantity.2 At an abstract level, the components of biological systems can be simplified to a set of nodes that are connected to each other by links, with each link representing the associations between two components.3 Two important features of a network are its topology and dynamics. Topology refers to the connectivity of the nodes in a network. Dynamics involve the quantification of the connections and time course response of the node states. Many diseases involve interactions of complex networks at different time and length scales. It is important to understand the fundamentals of the networks to be able to identify novel targets for interventions that may help prevent or cure diseases. The most common approach to model the dynamics of a biological network is to view it as a network of gene products, typically mRNA and proteins, and to describe their rate of changes through a system of differential equations.4 Therefore, the modeling framework is that of a continuous time and deterministic dynamical system cast as ordinary differential equations (ODEs). In many studies using ODEs, the main underlying assumption is that the system is operating near equilibrium states, so that the dynamics can be approximated by linear differential equations.5−8 With this assumption, the © 2015 American Chemical Society

network connections are quantified as constant weights valid around the neighborhood of the steady state. An interaction matrix, the Jacobian, represents the constant weights as a matrix. This matrix reveals the nature and strength of the interactions between biological components. This approximation method has been applied successfully by several research groups to reverse-engineer biological networks.5−17 Measured data from most biological experiments exhibit substantial noise arising from errors in the measurement technique and intrinsic/extrinsic variations in the cells or organisms.18,19 Measurement noise can be addressed by improving experimental techniques such as fluorescent protein fusion for single-molecule sensitivity,20 or by replicating the data. In practice, however, improving data quality and/or quantity is often not feasible, because biological systems involve many highly fluctuating dynamics and repeating experiments can result in the inclusion of more biological variations.21 Therefore, it is crucial to develop analytical methods that allow robust identification of biological networks in the presence of a high, but poorly defined measurement noise. The majority of inference algorithms for discrete time models have focused on the least-squares (LS) regression.22 The LS regression implicitly assumes additive error terms to the equation, which means that errors are limited to the dependent variables only. A significant problem in biological networks from the regression standpoint is that both the independent and dependent variables can have high levels of noise. Received: Revised: Accepted: Published: 10583

April 15, 2015 August 17, 2015 October 6, 2015 October 6, 2015 DOI: 10.1021/acs.iecr.5b01418 Ind. Eng. Chem. Res. 2015, 54, 10583−10592

Article

Industrial & Engineering Chemistry Research

system, its dynamics move from the equilibrium so that the linearized first-order equation that is described by eq 2 may no longer be valid. Instead, a nonlinear model in the form of eq 1 may be required to explain the dynamics and therefore have to be used in the estimation. This will be a significant leap in complexity, compared to the formulation presented in Section 3. Our present study is limited to the cases of perturbations sufficiently small not to cause significant nonlinear effects. Since the time series measurement data are available at discrete time points, most inference methods are developed using discrete time equations similar to that given below:

Moreover, these noise terms are serially correlated with each other. The total least-squares (TLS) method has been proposed to treat the cases where both the independent and dependent variables are corrupted by noise.23−25 However, the TLS method does not address the serial correlation structure in the noise. More advanced inference algorithms are needed that can take into account this aspect of the noise component. To this end, in the context of biochemical network inference, Kim et al.26,27 proposed an application of the constrained total least-squares (CTLS) algorithm28,29 that is suited to the model formulation and noise structure of the network identification problem. In their CTLS approach for a multivariable network model, however, the parameter estimation for each dependent variable is calculated separately. Since errors propagate in time and across variables via the interaction matrix, additional improvements can be expected by estimating the parameters for all the dependent variables simultaneously. In this paper, an extended CTLS algorithm is presented that accounts for the noise correlation structure and estimates parameters for all the dependent variables simultaneously. In addition, we endow the formulation with the flexibility of assigning different weights to the error terms based on the variance/covariance information on the measurement error. We compared our method with the original CTLS method, as well as other common regression methods on a benchmark genetic network example containing four genes.22,26,30,31 A further improvement in accuracy is observed over the original CTLS method, as well as the other conventional methods.

xk̃ + 1 = Fdxk̃ + uk

where xk̃ ∈ n × 1 is the expression of genes at sample time point k, Fd ∈ n × n is the parameter matrix transformed to discrete time domain, and uk ∈ n × 1 is the perturbation at time point k. The aim of inference methods is to estimate the parameter values in Fd by using noisy and discrete time series data, and transform them to continuous time domain parameters in F, thereby reconstructing the target biological network. The transformation of Fd to F can be done analytically, using the following relation based on the zero-order hold discretization,26

F=

F= (1)

(4)

2 (Fd − I )(Fd + I )−1 Δt

(5)

where I ∈ n × n is an identity matrix. However, this formulation can also lead to problems when (Fd + I) is close to singular. In such cases, the Euler approximation can be employed, which is given as

where the ith element, x̂i(t), of column vector x(̂ t ) ∈ n × 1 denotes the expression activity (or mRNA concentration) of the ith gene in the network at time t, and the ith element ui(t) of column vector u(t ) ∈ p × 1 represents an external input to the ith gene, such as the perturbation to the expression level of the ith gene at time t. f(·) is a nonlinear function describing the n-dimensional dynamics of the gene network. The dynamics (eq 1) around a given equilibrium x0 can be approximated with the linear differential equations obtained through the first-order Taylor series expansion, dx(̃ t ) = Fx(̃ t ) + Bu(t ) dt

1 log Fd Δt

where Δt is the sampling time interval. When Fd is poorly estimated, log Fd can be a complex matrix. To avoid such an unphysical result, a bilinear transformation method from discrete to continuous domain can be employed. The bilinear transformation is given as9

2. GENETIC NETWORK SYSTEM DESCRIPTION In general, a genetic regulatory network for n genes is represented by the following ordinary differential equation (ODE): dx(̂ t ) = f (x(̂ t ), u(t )) dt

(3)

F=

1 (Fd − I ) Δt

(6)

which is highly sensitive to the choice of Δt. Gene expression data are usually subject to high levels of additive and multiplicative errors.8,32 In the context of the genetic network model, this can be represented as xk = xk̃ + ek

(2)

(7)

where xk ∈ n × 1 is the noise-corrupted observation decomposed into the true concentration value x̃k and measurement error term ek ∈ n × 1, both of which are assumed to be unknown. The latter term can be further decomposed as follows:

where x(̃ t ) ∈ n × 1 denotes a deviation from the equilibrium state, x̂(t) − x0. F ∈ n × n is the Jacobian matrix of f(x̂(t),u(t)) for x̂(t) at x 0 and represents interactions between the individual genes in the network around the equilibrium state. A negative element in F indicates an inhibition and a positive element in F indicates an activation between the gene pair. A zero element in F means the absence of interaction for the pair. Matrix B ∈ n × p is the effect of the external perturbations and usually lumped into u(t ) ∈ n × 1 for simplicity.5 Note that the perturbation u(t) should be chosen to be sufficiently small, to keep the system close to a target equilibrium state. If large perturbations are applied to the

ek = Λkxk + bk

(8) n×n

where the diagonal matrix Λk ∈  corresponds to multiplicative measurement error, which is usually defined as zeromean white noises, Λk , ii ≈ 5(0, σΛ 2). The term bk ∈ n × 1 stands for the drift noise, which is typically modeled as the Brownian motion or random walk, represented by the relation 10584

DOI: 10.1021/acs.iecr.5b01418 Ind. Eng. Chem. Res. 2015, 54, 10583−10592

Article

Industrial & Engineering Chemistry Research bk = bk − 1 + Wkxk Δt

θLS, i = (ATA)−1ATφi

(9) n×n

where each element of the diagonal matrix Wk ∈  is 2 26 assumed to be zero-mean white noise, Wk , ii ≈ 5(0, σW ). Using eqs 3 and 7, one can write the model for all components as xk + 1 − ek + 1 = Fd(xk − ek) + uk

The LS solution assumes that errors are present only in the dependent variables. If error matrix E1 is assumed to be zero, then the least-squares solution is unbiased, meaning that it is the true solution in an average sense. 3.2. Total Least-Squares Formulation. For the cases where the E1 matrix (containing the noise in the independent variables) is significant, the TLS method was developed. The TLS method minimizes both E1 and E2 in a Frobenius norm sense.23 By adding error terms, the standard form (eq 14) can be rewritten as shown below:

(10)

Equation 10 can be written for all time points k = 1, ..., m − 1 as the following matrix equation: X 2 − E2 = Fd(X1 − E1) + U

(11)

(A + ΔA)θi = φi + Δφi

In this equation, the measurement matrices are represented as X2 = [x2, ..., xm] ∈ n × (m − 1) and X1 = [x1, ..., xm−1] ∈ n × (m − 1). Similarly, the error matrices are E2 = [e2, ..., em] ∈ n × (m − 1) and E1 = [e1, ..., em−1] ∈ n × (m − 1), and the input matrix is U = [u1 , ..., um − 1] ∈ n × (m − 1). One can see from eq 11 that both the dependent and independent variables have error terms. Furthermore, E1 and E2 are serially correlated as they contain the same columns except for the first and last columns. Input vector uk can be time-varying perturbations to the gene’s expression.22 For simplicity, the input vector is assumed to be constant, u ∈ n × 1 for all time points. Then, eq 11 can be rewritten as X 2 − E2 = [Fd

⎧ ⎡ X1 ⎤ ⎡ E1 ⎤⎫ ⎪ ⎪ ⎥⎬ ⎢ ⎥−⎢ u]⎨ ⎪ ⎪ × − 1 ( m 1) × − 1 ( m 1) ⎢ ⎥ ⎥⎦⎭ ⎦ ⎣⎢ 0 ⎩⎣1

(19)

where ⎡ E 1 ⋯ E p ⎤T 1 1⎥ ΔA = −⎢ ⎢⎣ 0 p × p(m − 1) ⎥⎦

(20)

Δφi = {−[ E21 ⋯ E2p ]i }T

(21)

or, more compactly, as [ θiT −1](C + ΔC) = 0

(22)

where

(12)

where 1 is the row vector of ones of size (m − 1) and 01×(m−1) is the row vector of zeros of size (m − 1). The model described by eq 12 can be further extended for parallel experiment case p as follows: 1×(m−1)

C = [ A φi ]

(23)

ΔC = [Δ A Δφi ]

(24)

The TLS formulation searches for each ith row of the parameter matrix at a time by minimizing the Frobenius norm of the combined error in the both sides of the equation. This can be represented as

[ X 21 ⋯ X 2p ] − [ E21 ⋯ E2p ] ⎫ ⎧⎡ X1 ⋯ X p ⎤ ⎡ 1 E1 ⋯ E1p ⎤⎪ ⎪ 1 1 ⎥ 1 p ⎨⎢ ⎥ ⎢ ⎬ − = [ Fd u ⋯ u ] ⎢ ⎪⎣ Ip ⊗ 11 × (m − 1)⎥⎦ ⎢⎣ 0 p × p(m − 1) ⎥⎦⎪ ⎭ ⎩

(18)

min || ΔC ||F

(25)

θi

subject to eq 22. The TLS solution for this system becomes (13)

θTLS, i = (ATA − λ 2I )−1ATφi

where Ip ∈  is the identity matrix and the symbol “⊗” represents the Kronecker product. p×p

where λ is the smallest singular value of C. Compared to the least-squares solution in eq 18, the TLS solution in eq 25 has a correction term λ2 noise in the independent variables. More detailed derivation can be found elsewhere.24 One of the main assumptions in the TLS approach is that the error terms, E1 and E2, are independent, which is not true in our case. 3.3. Constrained Total Least-Squares Formulation. In the case of the genetic network system model, E1 and E2 are serially correlated. The structure of this correlation can be incorporated into the estimation problem by the CTLS approach.28 The CTLS method also estimates each ith row of the parameter matrix for each dependent variable separately one at a time, similarly as in the TLS method. It differs from TLS in its objective function, which is defined as

3. PARAMETER ESTIMATION METHODS 3.1. Least-Squares Formulation. To formulate a standard LS estimation method, eq 13 can be generally expressed as Aθi = φi

(14)

⎡ X 1 ⋯ X p ⎤T 1 1 ⎥ A=⎢ ⎢ I ⊗ 11 × (m − 1)⎥ ⎣ p ⎦

(15)

φi = {[ X 21 ⋯ X 2p ]i }T

(16)

θi = {[ Fd u1 ⋯ u p ]i }T

(17)

(26) 24

⎡θ ⎤ min[ θiT −1]CT(HθiHθTi )C ⎢ i ⎥ ⎣−1⎦ θi

where A and φi include the measurement data, and θi contains the parameters for the ith gene to be estimated. The [·]i term denotes the row vector that contains the ith row of matrix [·]. The LS solution for the ith row of the parameter matrix can be written as

(27)

where Hθi is defined as n

Hθi =

∑ θiGs − Gn+ p + 1 s=1

10585

(28) DOI: 10.1021/acs.iecr.5b01418 Ind. Eng. Chem. Res. 2015, 54, 10583−10592

Article

Industrial & Engineering Chemistry Research Gs = {Ip ⊗ [(Im − 1 ⊗ es)T 0(m − 1) × n ]}

where Q kj ∈ n × n is the weighting matrix corresponding to the inverse of the covariance matrix of the respective noise term. Since the multiplicative measurement error is involved in the error term of eq 8, the weighting matrix Qjk can be simply chosen as

∀ s = 1 , ..., n (29)

Gn + p + 1 = {Ip ⊗ [0(m − 1) × n (Im − 1 ⊗ ei)T ]}

(30)

where ei ∈ n × 1 is the vector of zeros where only the ith

⎡(x j )−2 0 ⎢ k ,1 ⎢ 0 (xkj,2)−2 Q kj = ⎢ ⎢ ⋮ ⋮ ⎢ ⎢ 0 0 ⎣

element is equal to 1. Ip ∈ p × p and Im − 1 ∈ (m − 1) × (m − 1) are the identity matrices. The objective described by eq 27 minimizes the Frobenius norm of only the minimal set of the errors, thereby considering the correlation between errors E1 and E2. (See ref 26 and its Supporting Information for details of the derivations of eqs 27−30.) 3.4. Extended Constrained Total Least-Squares Formulation. It is clear from the model in eq 11 that noise terms E1 and E2 are correlated as a function of full parameter matrix Fd, rather than each of its rows. In our extended CTLS (ECTLS), we address this problem by reformulating the CTLS optimization. For that, one can start by rewriting eq 10 as ekj+ 1 = xkj+ 1 + Fd(ekj − xkj) − u j

e3j = x3j + Fd(e 2j − x 2j) − u j

(33)

Fd u1 ,..., u p 1 e1 ,..., e1p

(40)

In this minimization problem, the decision variables are Fd, uj, and ej1 for j = 1, ..., p, and the input variables are the measurements xjk for j = 1, ..., p and k = 1, ..., m. Our approach to the above optimization involves an iteration loop with two steps. The first step assumes that parameter matrix Fd is given and calculates uj and ej1 for j = 1, ..., p satisfying the first-order necessary condition for optimality. The first-order necessary condition consists of following equations with the gradients for uj and ej1:

Substituting eq 32 into eq 33 yields e3j = x3j + (Fd)2 (e1j − x1j) − (Fd + I )u j

(34)

A similar relation for all time points, k = 1, ..., m − 1, can be obtained as follows: ekj+ 1 = xkj+ 1 + (Fd)k (e1j − x1j) − Sku j k−1

Sk = (Fd)

k−2

+ (Fd)

+ ··· + I

(35) (36)

p

p

p

∀ j = 1 , ..., p

j=1 k=1

p

(37)

(42)

∑ ((Fd)k )T Q kj+ 1(Fd)k

(43)

(44)

k=1 m−1

Z=

∑ ((Fd)k )T Q kj+ 1Sk

(45)

k=1 m−1

rj =

∑ ((Fd)k )T Q kj+ 1xkj+ 1

(46)

k=1

Also, solving eq 42 leads to the following set of expressions (see the Appendix for details): (u j)T = ((d j)T + (e1j − x1j)T Z)(H )−1

m−1

j=1 k=1

(41)

m−1

R=

(47)

m−1

∑ (e1j)T Q 1je1j + ∑ ∑ (ekj+ 1)T Q kj+ 1ekj+ 1 j=1

∂Etotal =0 ∂u j

m−1

Etotal is a function of errors at the initial time step for all the experiments, as well as the inputs and measurements at all the time points. The error terms for all the components and time points can be substituted by the corresponding expressions from eq 35. Hence, eji emerges as an independent decision variable in the problem formulation. Advantage of our formulation is that the error function (eq 37) can be easily modified to accommodate further information about the measurement error. Specifically, the objective function can be generalized to contain weighted sum of the error terms, in order to incorporate their variance or covariance information as follows: Etotal =

∀ j = 1 , ..., p

(e1j)T = ((x1j)T R + (u j)T ZT − (r j)T )(R + Q 1j)−1

∑ (e1j)T e1j + ∑ ∑ (ekj+ 1)T ekj+ 1 j=1

∂Etotal =0 ∂e1j

Solving eq 41 yields the following equations (see the Appendix for details):

An error function is defined the sum of squared error terms. It is represented for all time points k = 1, ..., m − 1 and experiments j = 1, ..., p are represented as follows: Etotal =

(39)

min Etotal (Fd , u1 , ..., u p , e11 , ..., e1p|x11 , ..., xm1 , x12 , ..., xmp)

where ejk is the error term at the kth time point of the jth experiment. Similarly, xjk is the measurement data at the kth time point of the jth experiment. One can write eq 31 for k = 1 and k = 2 as follows: (32)

0

Note that, in the standard case, where the noise variance information is missing, we can also assume that the weighting matrix is represented by a scalar times the identity matrix. The extended constrained total least-squares approach can be simply seen as an unconstrained minimization of Etotal.

(31)

e 2j = x 2j + Fd(e1j − x1j) − u j

⎤ ⎥ ⎥ ⋯ 0 ⎥ ⋱ ⋮ ⎥ ⎥ ⋯ (xkj, n)−2 ⎥⎦ ⋯

j

d = (38)

∑ (Sk)T Q kj+ 1xkj+ 1 k=1

10586

(48) DOI: 10.1021/acs.iecr.5b01418 Ind. Eng. Chem. Res. 2015, 54, 10583−10592

Article


⎡ 1 + A (x ̂ (t )/K )n24 ⎤ ⎛ x ̂ (t ) ⎞ dx 2̂ (t ) 24 4 24a d 2 = V2s⎢ ⎟ n24 ⎥ − V 2 ⎜ dt ⎣ 1 + (x4̂ (t )/K 24a) ⎦ ⎝ k 2d + x 2̂ (t ) ⎠

m−1

H=

∑ (Sk)T Q kj+ 1Sk

(49)

k=1

(55)

By substituting eq 43 into eq 47 and rearranging the equation for uj, we can get an analytical solution for uj that satisfies the first-order necessary condition for optimality as follows:

⎡ ⎤ dx3̂ (t ) 1 + A32 (x 2̂ (t )/K32a)n32 = V3s⎢ n32 n31 ⎥ dt ⎣ (1 + (x 2̂ (t )/K32a) )(1 + (x1̂ (t )/K31i) ) ⎦ ⎛ x3̂ (t ) ⎞ − V 3d⎜ ⎟ ⎝ k 3d + x3̂ (t ) ⎠

(u j)T = ((d j)T Z −1(Q 1j + R ) − (x1j)T Q 1j − (r j)T )(HZ −1(Q 1j + R ) − ZT)−1

(50)

⎡ 1 + A (x ̂ (t )/K )n43 ⎤ ⎛ x4̂ (t ) ⎞ dx4̂ (t ) 43 3 43a d = V4s ⎢ ⎟ n43 ⎥ − V 4 ⎜ dt 1 ( x ( t )/ K ) + ⎣ ⎦ ⎝ k4d + x4̂ (t ) ⎠ 3̂ 43a

ej1

Then, an analytical solution for can be obtained from eq 43 by substituting solution (50). To check the convexity of the first step problem, the Hessian matrix for the objective function (eq 37) is calculated as (see the Appendix for details) 2 ⎤ ⎡ ∂ 2E total ∂ Etotal ⎥ ⎢ j j j 2 ∂e1 ∂u ⎥ ⎡ 2(R + Q 1j)T − 2Z ⎤ ⎢ (∂e1 ) ⎥ ⎥=⎢ ⎢ 2 T⎥ ⎢ ∂ Etotal ∂ 2Etotal ⎥ ⎢⎣ − 2ZT 2H ⎦ ⎢ j j j 2 ⎥ u e ∂ ∂ u ( ) ∂ ⎦ ⎣ 1

(57)

where x̂i(t) is the concentration of mRNAi for i = 1, ..., 4. Vsi and Vdi correspond to the maximum transcription and degradation rates, respectively, and Kij is the Michaelis constants. Parameter values for this model are given in Table 1. Table 1. Parameter Values for the Differential Equations eqs 54−57) Describing the Artificial Four-Gene Network System

∀ j = 1 , ..., p

Maximum Transcription Rates [nM/h]

(51)

parameter

All the eigenvalues of the Hessian matrix are positive, which indicates the first optimization problem is strictly convex. The second step is an optimization that searches for parameter matrix Fd minimizing the objective function in (eq 37) with the exception of the first term, for constant ej1 and uj for j = 1, ..., p given from the first step: p

min ∑ Fd

m−1

∑

(ekj+ 1)T Q kj+ 1ekj+ 1 (52)

j=1 k=1

In our implementation, we used the interior-point method in fmincon in the MATLAB Optimization Toolbox, which can handle a nonlinear multivariable objective function and constraints. The termination criterion for the iteration between the two steps can be set as ||(Fd)l − (Fd)l − 1 ||2 ||(Fd)l − 1||2

parameter

value

Vs1 Vs2 Vs3 Vs4 Michaelis Constants

5 3.5 3 4 [nM]

Vd1 200 Vd2 500 Vd3 150 Vd4 500 Other Constants

parameter

value

parameter

value

K14a K24a K32a K43a K12i K31i K1d K2d K3d K4d

1.6 1.6 1.5 0.15 0.5 0.7 30 60 10 50

A14 A24 A32 A43 n12 n14 n24 n31 n32 n43

4 4 5 2 1 2 2 1 2 2

Equilibrium state values for this system are given by x1,0 = 0.4920, x2,0 = 0.6052, x3,0 = 0.1866, and x4,0 = 0.6514. It is wellknown that a proper estimation performance requires a sufficient perturbation of the system. In this test case, timeseries data are obtained from constant parameter perturbation experiments. Four parallel simulated experiments are performed. In each experiment, one of the four VSi values are perturbed by 100% in the negative direction. Within each experiment, measurements are taken at a uniform sample time of 36 s for several time points varying from 3 to 21. Multiplicative measurement noise terms are added to the data in a manner shown in eq 8 (σΛ2 = 0.02).26 The corresponding true Jacobian matrix is obtained as follows:

where ϵ is taken as a sufficiently small number.

4. SIMULATION CASE STUDY AND DISCUSSION We applied our algorithm on an artificial four-gene network example as a simulation case study. This genetic network is modeled by nonlinear ordinary differential equations that are given in the Supporting Information of ref 30 It is used as a test case to evaluate the parameter estimation performance in several previous studies.8,22 Kim et al.26 also applied the proposed CTLS algorithm on this example under various conditions such as different noise types and quantities of data. The differential equations for the network are given as

⎡−6.45 −2.92 0 2.54 ⎤ ⎢ ⎥ −8.17 0 0 3.93 ⎥ F=⎢ ⎢ −2.31 2.80 − 14.46 0 ⎥ ⎢ ⎥ ⎣ 0 0 10.22 −9.74 ⎦

⎤ ⎡ 1 + A14 (x4̂ (t )/K14a)n14 dx1̂ (t ) = V1s⎢ n14 n12 ⎥ dt ⎣ (1 + (x4̂ (t )/K14a) )(1 + (x 2̂ (t )/K12i) ) ⎦ −

value

Maximum Degradation Rates [nM/h]

≤ϵ (53)

⎛

x1̂ (t ) ⎞ ⎟ ⎝ k1d + x1̂ (t ) ⎠

V1d⎜

(56)

(58)

In the parameter estimation, the true Jacobian matrix as well as the perturbation inputs are assumed to be unknown. The

(54) 10587


Article


The second error definition is based on the Frobenius norm of the difference between the true and estimated Jacobian matrices.

main objective is to estimate the Jacobian matrix with the generated dataset, using the proposed methods. The estimated parameter matrices from LS, TLS, CTLS, and two ECTLSs for data points of 6, averaged over the 1000 data sets are shown in Table 2. From the results, it can be seen that the estimated parameter matrix from ECTLS with the weighting parameter is the closest to the true Jacobian matrix.

εF = || F ̂ − F ||

The bilinear transformation is used to convert Fd to F, as given in eq 5. However, this transformation can lead to singularity problems. Therefore, the condition number is checked for each solution and the Euler approximation is employed (see eq 6) whenever the solution is close to singularity. Furthermore, the minimization problem solved by this algorithm is nonlinear and generally is nonconvex. Hence, the quality of the solution is dependent on the initial guesses for Fd. In this case, the solution provided by the LS algorithm can be used as an initial guess in the minimization problem. However, in some cases, it may happen that the LS solution does not provide a good initial guess and, as a result, the minimization algorithm may diverge. For such test cases, a boundary constraint on Fd can be applied as26

Table 2. Estimated Parameter Matrices from LS, TLS, CTLS, and Two ECTLSs for Data Points of 6 (Averaged over the 1000 Datasets) estimated parameter matrix, F̂

method

LS

⎡− 28.03 − 5.96 − 2.94 4.55⎤ ⎢ ⎥ 4.96 5.03⎥ ⎢ − 3.23 − 23.42 ⎢ − 3.23 2.63 − 21.07 0.32 ⎥ ⎢ ⎥ ⎣ 2.75 0.86 20.36 − 21.40 ⎦

TLS

⎡− 7.45 − 4.11 4.62 2.49 ⎤ ⎢ ⎥ 2.03 − 9.16 16.08 1.61⎥ ⎢ ⎢− 2.45 2.12 − 14.36 − 0.02 ⎥ ⎢ ⎥ ⎣ 3.47 3.22 5.50 − 11.72 ⎦

CTLS

⎡− 7.52 − 4.05 0.08 2.80 ⎤ ⎥ ⎢ − 0.86 − 8.88 2.35 3.46 ⎥ ⎢ ⎢ − 2.63 2.14 − 14.63 0.08 ⎥ ⎥ ⎢ ⎣ 0.80 − 0.81 21.31 − 10.67 ⎦

ECTLS, Qjk = I

⎡− 7.37 − 3.98 0.54 2.70 ⎤ ⎢ ⎥ 1.54 3.47 ⎥ ⎢− 0.59 − 8.66 ⎢− 2.60 2.12 − 14.65 0.10 ⎥ ⎢ ⎥ ⎣ 0.46 − 0.50 16.03 − 10.14 ⎦

ECTLS

⎡− 6.28 − 3.92 0.68 2.56 ⎤ ⎥ ⎢ 1.10 3.39 ⎥ ⎢− 0.38 − 7.95 ⎢− 2.55 2.13 − 14.45 0.04 ⎥ ⎥ ⎢ ⎣ 0.33 − 0.57 15.25 − 9.57 ⎦

−h|Fd ,LS| ≤ Fd ≤ h|Fd ,LS|

1 N1

n

n

∑ ∑ |αij| + i=1 j=1

1 N2

n

n

∑ ∑ |βij| i=1 j=1

(59)

⎧ f̂ − f ij ⎪ ⎪ ij for fij ≠ 0 αij = ⎨ fij ⎪ ⎪ ⎩ 0 otherwise

(60)

⎧ 0 for f ≠ 0 ij ⎪ βij = ⎨ ⎪ fiĵ otherwise ⎩

(61)

(63)

where h is asumed to have a value of 2 in this work. The parameter ϵ for the termination is set to 10−6. In the case study, all the methods are compared based on the two error criteria and the computation time (3.40 GHz CPU, 8GB RAM, MATLAB R2014b ver.) for different numbers of samples used, ranging from 3 to 21. CTLS refers to the original constrained total least-squares method. Our reformulation of CTLS is represented as ECTLS. ECTLS is further classified as ECTLS-1 and ECTLS-2. ECTLS-1 is the simpler version of the algorithm assuming the weighting matrix is defined as a scalar times the identity matrix (Qjk = I), while ECTLS-2 is the version considering the multiplicative measurement error with the weighting matrix (eq 39). The results are averaged over 1000 estimation trials done with Monte Carlo simulation data, and they are tabulated in Table 3. As it can be seen, for small numbers of data, the TLS method has a larger error, compared to the LS method. This is due to the minimum requirement of data for TLS.23 The ECTLS consistently gives the most accurate results, compared with all other regression methods, including CTLS. In particular, ECTLS-2 shows superior results to ECTLS-1 between the two ECTLS versions. The ECTLS-2 method reduces the εM error by an average of 6% and improves εF by an average of 10%, compared to CTLS. We also observed an average reduction of 11% and 22% for the standard deviations of εM and εF errors, respectively, for ECTLS, in comparison to the CTLS. It should be noted that the accuracy of the estimation increases as the number of data points used increases. In the case where only three data points are used, all the methods perform identically, since the number of equations is equal to that of the unknowns, which results in a single unique solution. However, the proposed ECTLS algorithm requires longer computation time, which increases rapidly as the number of samples in each experiment increases. The second step in the ECTLS algorithm takes a major part of the total computation time, as it is implemented to reformulate the optimization problem of eq 54. In particular, reconstructing the objective function involves p × (m − 1) multiplication operations of the unknown parameter matrix. In addition, the iteration of the

To test and compare the performance of the proposed ECTLS method with the original CTLS method, as well as with the LS and TLS methods, two different measures are adopted to quantify the deviation of the estimated Jacobian from the true one. The first error criterion is defined as ϵM =

(62)

where fiĵ and f ij are the ith row and jth column elements for matrices F̂ and F respectively. F is the true Jacobian matrix and F̂ is the estimated Jacobian obtained by applying one of the LS, TLS, CTLS, and ECTLS methods and then transforming the result to the continuous time domain using eqs 4−6. N1 and N2 terms are the number of nonzero and zero elements in the true Jacobian matrix. 10588


Article

Industrial & Engineering Chemistry Research Table 3. Estimation Results for the Four-Gene Network Example with Only White Multiplicative Noise εM method

mean

LS TLS CTLS ECTLS-1 ECTLS-2

94.46 94.46 94.46 94.46 94.42


16.18 63.85 14.83 12.94 12.64


7.47 9.95 5.34 5.21 5.14


4.89 5.53 3.01 2.99 2.95


3.24 3.14 1.73 1.68 1.65

Table 4. Estimation Results for the Four-Gene Network Example with Both Drift and White Multiplicative Noise and Nine Data Points

εF STD

mean

STD

3 Samples/Experiments 35.43 358.19 116.82 35.43 358.19 116.82 35.43 358.19 116.82 35.43 358.19 116.82 35.42 358.11 116.78 6 Samples/Experiments 5.13 66.98 17.38 231.23 480.90 3246.97 5.70 62.45 27.44 4.60 51.53 15.69 4.44 49.77 14.62 9 Samples/Experiments 2.32 32.89 8.55 7.85 49.65 104.84 1.95 23.02 7.67 1.82 21.69 6.48 1.77 21.16 6.17 12 Samples/Experiments 1.59 22.56 6.09 2.03 26.43 10.84 1.11 13.98 4.74 1.05 13.54 4.32 1.03 13.12 4.16 21 Samples/Experiments 0.95 15.71 3.98 1.11 16.06 6.23 0.55 8.64 2.13 0.55 8.45 2.15 0.54 7.86 1.95

εM

mean computation time (s)

method

0.0002 0.0002 0.0317 0.0153 0.0194


0.0002 0.0004 0.9612 3.3134 3.3415


0.0002 0.0004 0.8983 5.0435 4.4186


0.0002 0.0005 0.9015 6.0620 8.1012


0.0002 0.0007 1.1479 20.593 32.565

mean

εF STD

mean

Strength of Drift Noise = 2.0 12.04 4.29 51.33 35.22 190.56 180.85 10.66 5.58 46.07 9.41 3.40 38.19 9.12 3.27 36.71 Strength of Drift Noise = 1.0 8.71 2.87 37.47 12.12 7.31 60.92 6.75 2.51 28.44 6.50 2.31 26.63 6.35 2.24 25.76 Strength of Drift Noise = 0.1 7.52 2.38 32.88 9.78 6.56 47.18 5.33 2.05 23.03 5.24 1.92 21.71 5.15 1.88 21.12 Strength of Drift Noise = 0.05 7.43 2.37 32.73 9.76 4.31 45.65 5.38 1.93 22.82 5.22 1.78 21.70 5.14 1.75 21.08

STD 16.24 743.20 27.43 12.64 11.80 10.22 82.91 9.96 8.31 7.96 8.84 47.19 8.22 6.98 6.69 8.56 25.65 7.97 6.51 6.27

simulation data points are also tabulated in Tables S1 and S2 in the Supporting Information. It is obvious from Table 4 that, as the strength of the drift noise decreases, all the averaged error indices decreases steadily. Furthermore, the results in Table S1 and S2 are improved for all of the methods, compared to the results in Table 4, since more information is provided in the dataset for the parameter estimation. Most importantly, the proposed ECTLS algorithm shows significantly improved accuracy compared with all the other methods. Once again, ECTLS-2 is confirmed to be the most accurate method for the case having both drift and white multiplicative noise. Especially, the εM and εF are reduced by an average of 7% and 11% (13% and 23% for STD), respectively, compared to CTLS.

loop amplifies the number of matrix multiplications resulting in higher computation time. These aspects of the current algorithm implementation leave significant room for improvement in the speed. Even though the accuracy of ECTLS is always higher than CTLS, as shown in Table 3, the benefit of the rigorous error modeling in the ECTLS algorithm gets reduced as more and more information is contained in the dataset. The difference in performance of the two algorithms decreases as the number of data points increases, while the computation time of ECTLS is higher and increases more rapidly than CTLS. For many biological systems, however, datasets having sufficiently rich information is seldom accessible. In this context, ECTLS can be recommended to obtain more-accurate estimates (e.g., the cases of 6−9 data points in Table 3). To evaluate the estimation performance in the presence of the drift type noise, we applied the drift noise model in eq 9. For different strengths of the drift noise, covariance σW2 is changed gradually from 0.05 to 2.0. The averaged results for 1000 Monte Carlo simulation data points are tabulated in Table 4. The number of sample points in the dataset is kept constant at 9 and all the other conditions remain the same, as in Table 3. The effect of drift noise on the parameter estimation is further investigated for the increased number of sample points in the dataset from 9 to 12 and 21, while keeping all the other parameters same. The averaged results for 1000 Monte Carlo

5. CONCLUSION In this paper, we addressed the problem of biological network identification from noisy measurements. It has been shown that the resulting estimation model has noise terms in both the dependent and independent variables and they are serially correlated. The LS method is the most popular approach for such an estimation problem but the inherent assumption in the LS method is that noise is serially independent and confined to the dependent variables only. The TLS method is capable of taking errors in independent variables into account, but with no reference to possible correlations in the noise terms. The CTLS algorithm is a further improvement on TLS that can incorporate the correlation structure in the noise terms of the independent and dependent variables. Therefore, the CTLS method is well-suited to the estimation problem at hand. The existing study on the application of the CTLS method in the context of biological networks has reported significantly 10589


Article

Industrial & Engineering Chemistry Research m−1

∂{∑k = 1 (xkj+ 1)T Q kj+ 1(Fd)k (e1j − x1j) + (e1j − x1j)T ((Fd)k )T Q kj+ 1xkj+ 1}

improved estimation results, but was is formulated considering the errors in each dependent variable independently, leaving further room for improvement. In this work, we introduced an ECTLS method that simultaneously estimates the errors in all the dependent variables and across time. The ECTLS method can also incorporate the variance/covariance information in the error terms through weight parameters in the objective function. We demonstrated its improved performance over the existing CTLS method and other estimation methods on a four-gene network under various conditions, including different numbers of sample points and levels of noise. Although the CTLS methods seem to improve the parameter estimation significantly over the existing methods, the error levels are still high, despite reasonable noise levels in the data. The situation could be further improved by using network connectivity data based on prior knowledge, combined with optimal experimental design to achieve parameter estimation of high accuracy.

∂e1j m−1

=2

k=1

= 2(r j)T

2(e1j)T Q 1j + 2(e1j)T R − 2(x1j)T R − 2(u j)T ZT + 2(r j)T = 0 (A7)

⎡

(e1j)T

+

−

− x1j)

m−1

∂u j

+ (xkj+ 1)T Q kj+ 1(Fd)k (e1j − x1j)

m−1

+ (e1j − x1j)T ((Fd)k )T Q kj+ 1xkj+ 1 − (xkj+ 1)T Q kj+ 1Sku j ⎤ − (u j)T (Sk)T Q kj+ 1xkj+ 1 + (u j)T (Sk)T Q kj+ 1Sku j}⎥ ⎥⎦

= −2

∂Etotal =0 ∂e1j

= −2(d j)T

(A1)

m−1

∂u j m−1

= − 2(e1j − x1j)T =

(A2)

− 2(e1j

−

x1j)T Z

(A11)

m−1

∂{∑k = 1 (u j)T (Sk)T Q kj+ 1Sku j} ∂u j

m−1

= 2(u j)T

∑ (Sk)T Q kj+ 1Sk k=1

j T

= 2(u ) H (A12)

= 2(e1j)T Q 1j

The results from eqs A10−A12 are summed, and the summation is substituted into eq A9,

(A3)

−2(d j)T − 2(e1j − x1j)T Z + 2(u j)T H = 0

m−1

(A13)

Solving eq A9 for (uj)T yields

∂e1j m−1

m−1

∑ ((Fd)k )T Q kj+1(Fd)k k=1

− 2(x1j)T

(u j)T = ((d j)T + (e1j − x1j)T Z)(H )−1

∑ ((Fd)k )T Q kj+1(Fd)k k=1

= 2(e1j)T R − 2(x1j)T R

(A14)

To check the second-order optimality condition, a Hessian matrix of the following is needed:

(A4)

2 ⎡ ∂ 2E ⎤ total ∂ Etotal ⎥ ⎢ j j ⎢ (∂e1j)2 ∂e1 ∂u ⎥ ⎢ 2 ⎥ ⎢ ∂ Etotal ∂ 2Etotal ⎥ ⎢ j j j 2 ⎥ ⎣ ∂u ∂e1 (∂u ) ⎦

m−1

∂{∑k = 1 − (u j)T (Sk)T Q kj+ 1(Fd)k (e1j − x1j) − (e1j − x1j)T ((Fd)k )T Q kj+ 1Sku j} ∂e1j m−1

∑ (Sk)T Q kj+ 1(Fd)k k=1

= − 2(u j)T ZT

∑ ((Fd)k )T Q kj+1Sk k=1

∂{∑k = 1 (e1j − x1j)T ((Fd)k )T Q kj+ 1(Fd)k (e1j − x1j)}

= − 2(u j)T

(A10)

∂{∑k = 1 − (u j)T (Sk)T Q kj+ 1(Fd)k (e1j − x1j) − (e1j − x1j)T ((Fd)k )T Q kj+ 1Sku j}

This expression is calculated by taking the derivative of each term on the right-hand side of the expression for Etotal. The following expressions show the derivative of each term one by one:

= 2(e1j)T

∑ (xkj+ 1)T Q kj+ 1Sk k=1

The first step searches for matrix Fd that minimizes Etotal, with respect to ej1 and uj. The first-order necessary condition for optimality is related with two gradients for ej1 and uj. The equation for the first gradient is represented as

∂e1j

(A9)

∂{∑k = 1 −(xkj+ 1)T Q kj+ 1Sku j − (u j)T (Sk)T Q kj+ 1xkj+ 1}

− (u j)T (Sk)T Q kj+ 1(Fd)k (e1j − x1j) − (e1j − x1j)T ((Fd)k )T Q kj+ 1Sku j

∂{(e1j)T Q 1je1j}

(A8)

The derivation is done in similar steps from eqs A3−A6. The following expressions show the derivative of each term one by one:

k=1

x1j)T ((Fd)k )T Q kj+ 1(Fd)k (e1j

+ (u j)T ZT − (r j)T )(R + Q 1j)−1

∂Etotal =0 ∂u j

m−1

(e1j

=

yields

((x1j)T R

For the second gradient, the equation is represented as

∑ ⎢(e1j)T Q 1je1j + ∑ {(xkj+ 1)T Q kj+ 1xkj+ 1 ⎣ j=1 ⎢

(ej1)T

Solving for

APPENDIX By substituting eq 35 into eq 38, the total error becomes p

(A6)

The derivative of the total error, with respect to the initial error, is the summation of the derivative of its terms. Adding the expressions derived in eqs A3−A6, one can rewrite eq A2 as follows:

■

Etotal =

∑ (xkj+ 1)T Q kj+ 1(Fd)k

(A5) 10590

(A15) DOI: 10.1021/acs.iecr.5b01418 Ind. Eng. Chem. Res. 2015, 54, 10583−10592

Article


(6) Gardner, T. S.; Di Bernardo, D.; Lorenz, D.; Collins, J. J. Inferring genetic networks and identifying compound mode of action via expression profiling. Science 2003, 301 (5629), 102−105. (7) Tegner, J.; Yeung, M. S.; Hasty, J.; Collins, J. J. Reverse engineering gene networks: Integrating genetic perturbations with dynamical modeling. Proc. Natl. Acad. Sci. U. S. A. 2003, 100 (10), 5944−5949. (8) Sontag, E.; Kiyatkin, A.; Kholodenko, B. N. Inferring dynamic architecture of cellular networks using time series of gene expression, protein and metabolite data. Bioinformatics 2004, 20 (12), 1877−1886. (9) Bansal, M.; Della Gatta, G.; Di Bernardo, D. Inference of gene regulatory networks and compound mode of action from time course gene expression profiles. Bioinformatics 2006, 22 (7), 815−822. (10) Brynildsen, M. P.; Tran, L. M.; Liao, J. C. A Gibbs sampler for the identification of gene expression and network connectivity consistency. Bioinformatics 2006, 22 (24), 3040−3046. (11) Zavlanos, M.; Pappas, G. J.; Julius, A.; Boyd, S. Genetic network identification using convex programming. IET Syst. Biol. 2009, 3 (3), 155−166. (12) Zavlanos, M. M.; Julius, A. A.; Boyd, S. P.; Pappas, G. J. Inferring stable genetic networks from steady-state data. Automatica 2011, 47 (6), 1113−1122. (13) Ö ksüz, M.; Sadıkoğlu, H.; Ç akır, T. Sparsity as cellular objective to infer directed metabolic networks from steady-state metabolome data: a theoretical analysis. PLoS One 2013, 8 (12), e84505. (14) Wu, S.; Liu, Z.-P.; Qiu, X.; Wu, H., High-dimensional ordinary differential equation models for reconstructing genome-wide dynamic regulatory networks. In Topics in Applied Statistics; Springer: Berlin, Heidelberg, Germany, 2013; pp 173−190. (15) Studham, M. E.; Tjärnberg, A.; Nordling, T. E.; Nelander, S.; Sonnhammer, E. L. Functional association networks as priors for gene regulatory network inference. Bioinformatics 2014, 30 (12), i130−i138. (16) Hecker, M.; Lambeck, S.; Toepfer, S.; Van Someren, E.; Guthke, R. Gene regulatory network inference: data integration in dynamic modelsA review. BioSystems 2009, 96 (1), 86−103. (17) di Bernardo, D.; Thompson, M. J.; Gardner, T. S.; Chobot, S. E.; Eastwood, E. L.; Wojtovich, A. P.; Elliott, S. J.; Schaus, S. E.; Collins, J. J. Chemogenomic profiling on a genome-wide scale using reverse-engineered gene networks. Nat. Biotechnol. 2005, 23 (3), 377− 383. (18) Elowitz, M. B.; Levine, A. J.; Siggia, E. D.; Swain, P. S. Stochastic gene expression in a single cell. Science 2002, 297 (5584), 1183−1186. (19) Raj, A.; Rifkin, S. A.; Andersen, E.; van Oudenaarden, A. Variability in gene expression underlies incomplete penetrance. Nature 2010, 463 (7283), 913−918. (20) Taniguchi, Y.; Choi, P. J.; Li, G.-W.; Chen, H.; Babu, M.; Hearn, J.; Emili, A.; Xie, X. S. Quantifying E. coli proteome and transcriptome with single-molecule sensitivity in single cells. Science 2010, 329 (5991), 533−538. (21) Raj, A.; van Oudenaarden, A. Nature, nurture, or chance: stochastic gene expression and its consequences. Cell 2008, 135 (2), 216−226. (22) Schmidt, H.; Cho, K. H.; Jacobsen, E. W. Identification of small scale biochemical networks based on general type system perturbations. FEBS J. 2005, 272 (9), 2141−2151. (23) Golub, G. H.; Van Loan, C. F. An analysis of the total least squares problem. SIAM J. Numer. Anal. 1980, 17 (6), 883−893. (24) Van Huffel, S.; Vandewalle, J. The Total Least Squares Problem: Computational Aspects and Analysis; Society for Industrial and Applied Mathematics (SIAM): Philadelphia, PA, 1991; Vol. 9. (25) Jacklin, N.; Ding, Z.; Chen, W.; Chang, C. Noniterative convex optimization methods for network component analysis. IEEE/ACM Trans. Comput. Biol. Bioinf. 2012, 9 (5), 1472−1481. (26) Kim, J.; Bates, D. G.; Postlethwaite, I.; Heslop-Harrison, P.; Cho, K.-H. Least-squares methods for identifying biochemical regulatory networks from noisy measurements. BMC Bioinf. 2007, 8 (1), 8. (27) Montefusco, F.; Cosentino, C.; Kim, J.; Amato, F.; Bates, D. G. Reverse engineering partially-known interaction networks from noisy

Each element in the matrix (eq A15) can be obtained by differentiating the left side of eq A7 (or eq A13), with respect to ej1 (or uj), ∂ 2Etotal (∂e1j)2

= 2(R + Q 1j)T

(A16)

∂ 2Etotal = −2Z ∂e1j∂u j

(A17)

∂ 2Etotal T j j = − 2Z ∂u ∂e1

(A18)

∂ 2Etotal (∂u j)2

= 2HT

(A19)

Then, the Hessian matrix is given as 2 ⎡ ∂ 2E ⎤ total ∂ Etotal ⎥ ⎢ j j ⎢ (∂e1j)2 ∂e1 ∂u ⎥ ⎡ 2(R + Q 1j)T −2Z ⎤ ⎥ ⎢ 2 ⎥=⎢ 2 ⎢ ⎥ T T ⎢ ∂ Etotal ∂ Etotal ⎥ ⎣ −2Z 2H ⎦ ⎢ j j j 2 ⎥ ⎣ ∂u ∂e1 (∂u ) ⎦

(A20)

which is positive-definite.

■

ASSOCIATED CONTENT

S Supporting Information *

This material is available free of charge via the Internet at http://pubs.acs.org/. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.5b01418. Two tables, including additional estimation results for drift noise cases (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*Tel.: +82-42-350-3926. Fax: +82-42-350-3910. E-mail: [email protected]. Author Contributions §

These authors contributed equally.

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the Advanced Biomass R&D Center (ABC) of Global Frontier Project funded by the Ministry of Science, ICT and Future Planning (No. 20110031354).

■

REFERENCES

(1) Kitano, H. Systems biology: A brief overview. Science 2002, 295 (5560), 1662−1664. (2) Kitano, H. Computational systems biology. Nature 2002, 420 (6912), 206−210. (3) Barabasi, A.-L.; Oltvai, Z. N. Network biology: Understanding the cell’s functional organization. Nat. Rev. Genet. 2004, 5 (2), 101−113. (4) De Jong, H. Modeling and simulation of genetic regulatory systems: A literature review. J. Comput. Biol. 2002, 9 (1), 67−103. (5) Yeung, M. S.; Tegnér, J.; Collins, J. J. Reverse engineering gene networks using singular value decomposition and robust regression. Proc. Natl. Acad. Sci. U. S. A. 2002, 99 (9), 6163−6168. 10591


Article

Industrial & Engineering Chemistry Research data. In Proceedings of the 18th IFAC World Congress, Milano, Italy, Aug. 28−Sept. 2, 2011; International Federation of Automatic Control: Laxenburg, Austria, 2011; pp 11679−11684. (28) Abatzoglou, T. J.; Mendel, J. M. Constrained Total Least Squares. In ICASSP 87: Proceedings, 1987 International Conference on Acoustics, Speech, and Signal Processing, Dallas, TX, April 6−9, 1987; IEEE: New York, 1987; pp 1485−1488. (29) Abatzoglou, T. J.; Mendel, J. M.; Harada, G. A. The constrained total least squares technique and its applications to harmonic superresolution. Signal Process., IEEE Trans. 1991, 39 (5), 1070−1087. (30) Kholodenko, B. N.; Kiyatkin, A.; Bruggeman, F. J.; Sontag, E.; Westerhoff, H. V.; Hoek, J. B. Untangling the wires: a strategy to trace functional interactions in signaling and gene networks. Proc. Natl. Acad. Sci. U. S. A. 2002, 99 (20), 12841−12846. (31) Stigter, J. D.; Molenaar, J. Network inference via adaptive optimal design. BMC Res. Notes 2012, 5 (1), 518. (32) Ideker, T.; Thorsson, V.; Siegel, A. F.; Hood, L. E. Testing for differentially-expressed genes by maximum-likelihood analysis of microarray data. J. Comput. Biol. 2000, 7 (6), 805−817.

10592


An Extended Constrained Total Least-Squares Method for the

Recommend Documents