A Computationally Efficient Algorithm for Testing ... - ACS Publications

Ind. Eng. Chem. Res. 2010, 49, 6125–6134

6125

A Computationally Efficient Algorithm for Testing the Identifiability of Polynomial Systems with Applications to Biological Systems Amos Ben-Zvi* Department of Chemical and Materials Engineering, UniVersity of Alberta, Edmonton, Alberta, Canada

Mathematical models of physical systems often contain many more parameters than can be estimated from observations. It is useful to check whether such models are structurally identifiable before choosing an experimental plan for parameter estimation. Many identifiability tests have been proposed in the literature. However, the proposed methods are computationally complex, or even intractable, for systems containing more than a handful of states or parameters. In this work, a linear algebra-based approach for testing ordinary differential equations and index-one differential algebraic equation systems with linear output maps for local structural identifiability is presented. The proposed method is computationally efficient, as it does not require repeated differentiation of the model equations. Furthermore, the proposed model can be used by experimenters to determine which set of measurements should be made in order to estimate specific parameters within the model. The effectiveness of the proposed approach is demonstrated by testing the identifiability of a 15-state model describing NF-κB regulation. The model is shown to be identifiable even if process measurements are related to the state variables by proportionality constants (themselves unknown). Furthermore, the proposed procedure is used to compute a table that relates to each unknown parameter or parameter set a set of state variables whose observation is sufficient for parameter estimation. Introduction Mathematical models of physical systems typically contain more parameters than can be estimated from available data. Lack of parameter estimability may be due to the quality and quantity of data available. Alternatively, the model may be structurally unidentifiable. In this latter case, no additional data (regardless of quality or quantity) is sufficient to uniquely identify every parameter in the model. Structural identifiability implies that there exists some set of tests which can be used to uniquely estimate every parameter. Lack of structural identifiability implies that the effect of some parameters cannot be ascertained from the model predictions regardless of experimental design. These unidentifiable (and hence superfluous) parameters can be removed from the model (or fixed at an arbitrary value) without affecting the predictive ability of the model. Furthermore, if these unidentifiable parameters are not removed from the model, then the parameter estimation problem is ill-posed. The issue of structural identifiability has been well studied in the literature.1 Domains of applications include: models of HIV dynamics,2 biotechnology,3 reactor systems,4 water treatment,5 metabolics,6 aerospace,7 and photophysics.8 Identifiability is assessed under ideal conditions (e.g., noiseless measurements, perfect input actuation). Identifiability is a necessary, but not sufficient, condition for experimental estimability which is the ability to estimate parameters from real-world data. The class of mathematical models considered in this work include sets of coupled, nonlinear ordinary differential equations (ODEs), and index-one differential algebraic equation (DAE) systems. Several methods have been proposed in the literature testing an ODE model for structural identifiability. These include a system-wide isomorphism approach,9,10 a power series-based approach,11 a functional generating series-based approach,12 a geometry-based approach,13 and a linearization-based approach.14 Methods that are suitable for testing the identifiability * To whom correspondence should be addressed. E-mail: abenzvi@ ualberta.ca. Tel.: 1-780-492-7651. Fax: 1-780-492-2881.

of DAE systems include the differential algebraic approach2,15,16 and a linearization-based approach.17 With the exception of the state isomorphism approach, each of the proposed approaches relies, in practice, on repeated differentiation of the model equations. As a result, the equations and computations required for testing a model for identifiability often become intractable for systems containing more than a handful of states or parameters. This issue has not been ameliorated by the increase in computing power. As recently as 2008, Jimenez-Horneo et al.3 compared several approaches for identifiability of a six-state, nine-parameter model and concluded that “in some cases computational complexity is so high that the problem is intractable, even with powerful software”. Similarly, in their work, Fujarewicz et al.18 state that “unfortunately, the system of differential equations analyzed in this paper is highly complicated and the so-called similarity transformation cannot be performed easily”. The need for a “computationally efficient” test for identifiability was highlighted in a recently published manuscript describing a novel numerical method for testing identifiability of models over a predefined interval.19 In this manuscript, the authors state that “the proposed technique only works for low dimensional problems, since the proposed algorithm takes an exponential amount of time in terms of the dimension”. It should be noted, however, that the method proposed by Lagrange et al.19 is applicable to nonpolynmial systems, whereas the method proposed in this work is restricted to polynomial systems. Linearization-based approaches14,17 provide a computationally efficient test for identifiability. However, these approaches provide only sufficient conditions for identifiability. That is, even if the linearized system is found to be unidentifiable, the original nonlinear system may be identifiable. While it has been shown that a linearization-based test can be used to treat unidentifiable nonlinear systems,4 this is only possible if the results obtained from the linearized system analysis are combined with numerical simulations and a heuristic-based approach for parameter reduction.

10.1021/ie9018512  2010 American Chemical Society Published on Web 06/02/2010

6126

Ind. Eng. Chem. Res., Vol. 49, No. 13, 2010

In this work, a computationally efficient algorithm for testing the identifiability of ODE and index-one DAE systems is presented. The approach relies only on linear algebra to generate a set of constraints on the set of possible model parameters. Furthermore, the proposed approach does not rely on a linearization of the model. Rather, it relies on describing a set of polynomials (representing the model equations) as a set of linear functions (in monomials of the model states and outputs). As a result, the approach is restricted to models that are polynomial in the inputs, states, and outputs. The proposed approach yields only sufficient conditions for identifiability. However, the proposed approach is also constructive in the following sense: it can be used to associate with every parameter (sub)set, a set of outputs whose measurement is sufficient for identifiability of the chosen parameter set. This property is useful for cases where only a subset of the model parameters are of interest to the experimenter. Furthermore, in most practical applications, there is a limit to the number and quality of measurements that are possible. By associating a subset of the measurements with a subset of the parameters, the proposed procedure allows experimenters to focus on observations that are critical for the estimation of key parameters. Background System Description. The class of systems under consideration in this work is the set of observable nonlinear ODE and index-one DAE systems relating the states x(t) ∈ Rn, inputs u(t) ∈ Rs, and parameters θ ∈ P ⊆ Rp in the following manner: F(x˙(t), x(t), u(t), θ) ) 0

(1)

where the mapping F is polynomial in x˙, x, u, and θ with ∂F/∂x˙ full rank on a dense subset of the state space. In the case of ODE systems, ∂F/∂x˙ ) In×n and the mapping F can be written as F(x˙(t), x(t), u(t), θ) ) x˙(t) - f(x(t), u(t), θ). In this work, it will be assumed that the measured variables y(t) ) C(θ)x(t) are related to the states by an invertible matrix C(θ). While this assumption of linearity is somewhat restrictive, it is often the case that the relation between the measured variables and the system states can be approximated by a static gain over a sufficiently small interval (in the states).18 The more restrictive assumption that the matrix C(θ) is invertible for θ ∈ P can be relaxed postanalysis. That is, under the proposed framework, one would assume that all states are observable and then determine which observations can be omitted while maintaining identifiability. The restriction of System 1 to being polynomial in the parameters is not as restrictive as one might imagine. The class of systems described by eq 1 include ordinary differential equation systems of the form x˙(t) ) f(x(t), u(t), θ) where the mapping f contains a ratio of polynomials in θ including, for example, Monod,20 Haldane,21 Caperon-Meyer,22 and Lehman23 kinetics. Physical systems that have been modeled using equations that are polynomial in the parameters include distillation systems and isothermal reactors,24 the human hypothalamic-pituitary-adrenal axis,25 and modeling of viral and immune system dynamics.26 Furthermore, as pointed out by Ljung and Glad,27 it is possible to represent some nonpolynomial functions (e.g., x ) sin(y)) in the form of polynomial equations in the variables and their time-derivatives. To study the structural properties of the model, it is assumed that the initial condition of System 1 is arbitrary (e.g., not necessarily zero). This assumption allows identifiability to be assessed independent of the controllability of the system. A

method for assessing the identifiability of nonlinear systems with fixed initial conditions was proposed by Saccomani et al.16 Let System 1 be defined over a time interval u ) [to, tf] with tf > to. The input of System 1 can be viewed as a sufficiently differentiable mapping u:u f Rs. As a result, for any specified input, System 1 can be seen as generating a parameter-dependent mapping (or behavior) M(θ): (u, x(0)) f y which assigns to each input trajectory and initial condition a unique output trajectory defined over a time interval u. The model inputoutput behavior can be used to define local identifiability. Local Identifiability. A model is locally structurally identifiable if there is a locally unique (bijective) relationship between the model parameters θ and the mapping M(θ)( · , x(0)) which describes how an input signal, u, is transformed by the system to produce an output signal, y. Formally, System 1 is defined as locally identifiable about θo ∈ P if for any other point θ˜ ∈ P with |θo - θ˜ | < ε, for an arbitrarily small ε, the following holds (on all but a countable subset of P): M(θo) ) M(θ˜ ) ⇒ θo ) θ˜ This definition of local identifiability is similar to the one used by Ljung.28 Generally, the mapping M cannot be solved for analytically. Rather, a set of equations is typically obtained which characterize M. This set is used to compute constraints on the set of possible parameter values. For linear time invariant (LTI) ODE systems, the mapping M can be characterized by the coefficient in the transfer function relating the input and initial condition to the output.29 For nonlinear ODE systems, the mapping M can be characterized by the coefficients in a Taylor or a power series expansion of y as a function of time.3,30 For polynomial DAE systems, M can be characterized by coefficients in a characteristic set of the differential polynomials which make up System 1.31 The existence of this characteristic set has been shown by Ritt,32 who also provides an algorithm for its construction. As a practical matter, the computation of the characteristic set for a set of polynomials is computationally intensive for more than a handful of variables (e.g., inputs, outputs, and states). For differential polynomials, this computational complexity is even higher due to repeated symbolic differentiation. In this work, a method for obtaining sufficient conditions for local identifiably of System 1 is introduced. The approach is based on the idea that one does not need to completely compute the characteristic set of a system to determine identifiability. Linear Equations. Consider a linear system of equations g(ξ) ) A(ξ)φ

(2)

where ξ ∈ RR is a set of known quantities, g(ξ) ∈ Rβ, A(ξ) is an β × γ matrix, and φ ∈ Rγ. System 2 does not have a unique solution for φ if the rank of A is less than γ. If the rank of A is less than γ then only a (possibly empty) subset of the elements in φ can be computed from eq 2. The subset of φ which can be uniquely determined from eq 2 can be obtained by computing a row-reduced echelon matrix that is equivalent to A(ξ).33 Specifically, let the matrix Q(ξ) define a (finite) set of elementary row operations such that the product matrix Q(ξ)A(ξ) is in rowreduced echelon form. Then the solutions to eq 2 are also the solutions to Q(ξ)g(ξ) ) Q(ξ)A(ξ)φ. Furthermore, if the vector φ is defined by a mapping φ ) Φ(θ), where θ ∈ Rp are the system parameters, then by the inverse function theorem, there exists (locally about a nominal value θ0 ∈ Rp) an inverse mapping Φ-1: φ f θ, if and only if


6127

the rank of the Jacobian ∂Φ/∂θ|θ)θ0 is equal to p. Let M ) A(ξ)(∂Φ/∂θ)|θ)θ0. A necessary and sufficient condition for the local existence of a unique solution for θ in eq 2 with φ ) Φ(θ) is that

where A(ξ) is an n × γ matrix. Letting g(ξ) ) -a0(ξ), eq 5 can be written as

rank[M ] ) p

Along solutions of System 1, F(ξ) ) 0 (by definition) and eq 6 is equivalent to

34

(3)

This is because M ) ∂(A(ξ)oΦ)/∂θ can be seen as the Jacobian of the composition mapping A(ξ)oΦ with respect to θ. The rank condition in eq 3 implies that A(ξ)oΦ has an inverse locally about θo which can be used to compute θ from known quantities. This condition in eq 3 forms the basis for the proposed identifiability test. The Proposed Method The proposed method is divided into two parts. First, the model equations are placed in the form of eq 2. Second, the rank of M is tested using symbolic computation software. Note that this step is not computationally intensive for a large class of systems, as it requires only an implementation of linear elimination algorithms (e.g., Gaussian elimination). The algorithm provides sufficient conditions. If rank(M ) ) p, then the system is concluded to be locally identifiable. If rank(M ) < p, then the results of the proposed test are indeterminate. Furthermore, the null space of M can be used to determine which parameters may not be identifiable.4 The first step is to convert the model equations to the form of eq 2. This is done by first finding the state variables from the model equations using the relation x(t) ) C-1y(t). The next step is to transform the mapping F into the form of eq 2. Let ξ ) (y, y˙, u) be a set of directly measured or known quantities. The mapping F in System 1 is polynomial in θ. By definition, this implies that F can be written as γ

F(ξ, θ) ) a0(ξ) +

∑ a (ξ)Φ (θ) γ

i

(4)

i)1

for some non-negative integer γ, where ai ∈ Rn are coefficients of the function Φi(θ), which is a rational function of monomials in θ. That is, the function F is linear in ratios of monomials of θ. Note that, under redefinition of Φi and γ, the coefficient a0 can always be made nonzero (as a function of ξ). For example, if a0 is identically zero, one may redefine Φi for i ) 1, 2, · · · , γ as Φγ ) 1 and Φi f Φi/Φγ, and redefine aγ as a0 and γ as γ - 1. As an example, consider the one-state system given by F ) (x + θ1)x˙2 - x3θ2 + θ3

θ43

+

ξ12θ1 θ42

-

ξ23θ2 θ43

+ θ3

Recalling that the output mapping is assumed invertible (i.e., θ4 * 0), the mapping F can be written as ξ12ξ2 + ξ12θ1θ4 - ξ23θ2 + θ3θ43 which is in the form of eq 4 with a0(ξ) ) ξ12ξ2, a1(ξ) ) ξ12, a2(ξ) ) -ξ23 and a3(ξ) ) 1. The linearity of F in Φ(θ) implies that F can be written as F(ξ, θ) ) a0(ξ) + A(ξ)Φ(θ)

g(ξ) ) A(ξ)Φ(θ)

(6)

(7)

with every entry in the mapping g nonzero as a function of ξ. Equation 7 characterizes the behavior, M, of System 1. As shown in the previous section, eq 7 has a unique solution (at a particular time t ∈ u ) if and only if the rank condition in eq 3 holds. As a result, eq 3 provides sufficient conditions for the identifiability of System 1. However, eq 7 does not provide necessary conditions for identifiability, because the relations in eq 7 must hold at eVery time t ∈ u. More generally, eq 7 imposes a set of n constraints on the set of possible parameter values (for every t ∈ u ). However, the system may have more parameters than states. As a result, it may be necessary to evaluate eq 7 at several points in time to generate additional constraints on the set of possible parameter values. Evaluating eq 7 at different times t1 < t2 < · · · < tζ all in u yields the following equation g˜(ξ) ) A˜(ξ)Φ(θ)

(8)

with g˜ ) [g(ξ(t1))T, g(ξ(t2))T, ..., g(ξ(tζ))T]T T T T T A˜(ξ) ) [A(ξ(t1)) , A(ξ(t2)) , ..., A(ξ(tζ)) ]

Equation 8 has, locally, a unique solution to θ if and only if there exists some choice of input such that

(

rank

)

∂A˜(ξ)oΦ )p ∂θ

(9)

Equation 8 can be used to determine a lower bound for the number of times that must be used to evaluate identifiability. Given that A(ξ) is a β × γ matrix and therefore A˜(ξ) is a ζβ × γ matrix. A necessary condition for the rank of A˜(ξ) to be greater than p is that ζβ g p. By assumption β > 0 is a positive dimension and therefore this condition can be rewritten as ζ g pβ-1; that is, it is not possible to assess identifiability if ζ < pβ-1. As an example, consider the trivial one-state system y˙(t) ) θ1y(t) + θ2u(t)

with y ) θ4x. The mapping F can be written in terms of the known quantities ξ1 ) y˙, ξ2 ) y, and ξ3 ) u as ξ12ξ2

F(ξ) ) -g(ξ) + A(ξ)Φ(θ)

(5)

This system can be put in the form of eq 7 with g(ξ) ) y˙(t), A ) [y(t), u(t)], and Φ(θ) ) θ ) [θ1, θ2]T. In this example, the number of independent equations describing the system dynamics is β ) 1 and the number of parameters is p ) 2. As a result, the system must be evaluated at two or more points in time to obtain conditions for identifiability of both θ1 and θ2. Consider t1 ∈ T and t2 ∈ T with t1 * t2. Equation 7, evaluated at t1 and t2, can be written in the form of eq 8 with g˜ ) [y˙(t1), y˙(t2)]T, Φ(θ) ) θ, and A˜(ξ) )

[

y(t1) u(t1) y(t2) u(t2)

]

The matrix A˜(ξ) in this example is square and therefore there exists a unique solution for θ if and only if det (A˜(ξ)) * 0 for some choice of input. It is important to note that the determinant

6128


function is a polynomial in the matrix entries and, as a result, the condition det (A˜(ξ)) ) 0 implies that there must be algebraic constraints on the process inputs and outputs (as functions of time). As the proposed method is restricted to ODE and indexone DAE systems with arbitrary initial conditions, no algebraic constraints exist on the system outputs other than the ones defined by eq 3. More generally even for the nonsquare case, rank deficiency of A˜(ξ) implies that any submatrix consisting of p rows of A˜(ξ) must have a zero determinant. The requirement to evaluate eq 3 at several points in time to obtain eq 8 is made in order to avoid having to generate additional constraints by differentiating eq 3. That is, instead of choosing two points in time at which to evaluate eq 3, one may, alternatively, choose to differentiate eq 3 to obtain the additional constraints ∂g ∂A dξ dξ ) (ξ) · Φ(θ) (ξ) · ∂ξ dt ∂ξ dt

(10)

where (∂A/∂ξ)(ξ) is a tensor of appropriate dimensions. However, the differentiation of eq 3 typically results in very large expressions even for very simple models.17 This differentiation can instead by approximated by choosing two time values t1 and t1 + δ for an arbitrarily small δ > 0 such that both t1 and t1 + δ are in u. The resulting set of equations can be rearranged to obtain g(ξ(t1 + δ)) - g(ξ(t1)) (A(ξ(t1 + δ)) - A(ξ(t1))) Φ(θ) ) δ δ (11) Taking the limit as δ f 0 in eq 11 one recovers eq 10. Obtaining additional algebraic constraints by evaluating eq 8 at different time values can therefore be used to recover the results that one would obtain by taking time derivatives, while greatly reducing the related computational cost. The proposed procedure, however, still maintains the inherent limitations of a derivative-based approach. Namely, the conditions for identifiability are only sufficient. One can never rule out identifiability because it is always possible, by considering additional points in time, that identifiability may be verified. However, the ability to select many different times at which to evaluate eq 7 overcomes the limitation of some expansion-based approaches where, typically, the initial time is used. This is because the initial time may be uninformative if, for example, the timederivatives of some states are zero at the initial time. In the proposed method, it is assumed that the system outputs are related to states by an invertible matrix C(θ). However, assumption of invertibility for C(θ) can be relaxed. Consider a partitioning of the vector of outputs into a vector of observed variables ya ) [y1, y2, · · · , yk] and a set of unobserved variables yb ) [yk+1, yk+2, · · · , yn], for some k < n. Without loss of generality, one may choose ya and yb such that

[ ] [ ]

ya(t) Ca(θ) y(t) ) ) x(t) yb(t) Cb(θ) Consider a subset ξa ) [ya, u] of the elements in ξ. Without loss of generality, let ξ ) [ξa, ξb], with ξb ) yb. Consider a situation where there exists an invertible matrix H(ξ) such that the equation H(ξ)g˜(ξ) ) H(ξ)A˜(ξ)Φ(θ) can be partitioned as

gã(ξa) ) Aã(ξa)Φ(θ) g˜b(ξa, ξb) ) A˜(ξa, ξb)Φ(θ) Then the observations ya associated with ξa are sufficient for identifiability if

(

rank

)

∂Aã(ξ)oΦ )p ∂θ

(12)

That is, if eq 12 holds, then it is possible to measure the effect of every parameter uniquely only using the known variables ξa. Stated equivalently, the outputs yb are not necessary for identifiabilty. As a result, the assumption of an invertible matrix C(θ) can be replaced with a nonsquare matrix Ca(θ) so that ya ) Ca(θ)x. The proposed algorithm has several key features. First, it is computationally efficient as it requires only application of row reduction to a matrix containing polynomials of the inputs, outputs, and the first time-derivative of the outputs. Furthermore, in many practical applications (including the example in this work) the computation of a large number of derivatives is not necessary to obtain conditions for identifiability. This row reduction can be efficiently accomplished using symbolic manipulation tools such as Maple (Maplesoft, 2010). Second, the proposed test is local due to the local scope of the application of the inverse function theorem. Finally, if the condition in eq 8 does not hold, then the set of solutions Θ ) {θ ∈ Rp|g˜(ξ) ) A˜(ξ)Φ(θ)} is the set of all possible parameter values. Numerical Example. Nuclear factor κB (NF-κB) is a protein and a transcription factor that plays an important role in the body’s immune response, including the regulation of cellular response to stresses such as viral and bacterial infection. A mathematical model of the regulatory module of NF-κB has been published.35 This model includes the two-compartment kinetics of the activators IκB (IKK) and NF-κB, the two inhibitors A20 and IκBR, along with their complexes. The model proposed in ref 35 is made up of a set of 15 ordinary differential equations (ODEs). The model contains 29 parameters. Only a subset of these parameters possess experimentally derived estimates that have been published in the literature. In Fujarewicz et al.18 an adjoint sensitivity-based approach for fitting mathematical models is presented. This approach is applied to a model of the NF-κB regulatory module.35 The model fitting is done by comparing model predictions with observations made using electrophoretic mobility shift assay (EMSA). EMSA is a technique for determining the extent of protein-DNA interaction and, indirectly, the concentration of different proteins and their complexes. Furthermore, the relationship between protein (or complex) concentration and the measured signal is approximated, for a single experiment, by yi ) wixi where yi is the quantification of the measured signal i, xi is the concentration of protein i, and wi is a proportionality constant (to be estimated) for i ∈ {1, 2, · · · , 15}. After fitting the data and obtaining parameter estimates for the unknown parameters and the proportionality constants w1 to w15, Fujarewicz et al.18 note the following: “The preliminary conclusion is that there is a large set in the parameter space resulting in similar behavior of the model and similar value of the performance index.... It may be possible that not all parameters can be estimated from the data or only some functions of the parameters can be estimated. Unfortunately the system of differential equations analyzed in this paper is highly


complicated and the so-called similarity transformation cannot be performed easily”. Furthermore, it is reported that the fitted parameters obtained during cross-validation differ by as much as 2 orders of magnitude from the nominal fitted values. The ability to estimate every parameter value accurately depends, in practice, on the quantity and quality of data that is available. Sensitivity analysis36 is the tool most commonly used to assess the effect of parameters on model predictions. Sensitivity analysis has previously been applied to analyze models of the NF-κB signaling pathways.37,38 Parameter sensitivity, however, only provides a local (in the parameter space) approximation of the effect of parameters on the model predictions. Furthermore, the results obtained from sensitivity analysis are specific to the initial conditions chosen for simulation and the level of noise in the measured signal. Sensitivity is, therefore, not a structural property of the model, but a property of a specific experimental design applied to the model. The following two questions regarding the NF-κB model are treated in this work: (1) Identifiability. Can the value of every unknown parameter be independently identified, at least under ideal conditions? (2) Choice of Measurements. Which observations are necessary to estimate which parameter? In this example, it will be shown that unknown parameters in the NF-κB model can be estimated (including the unknown proportionality constants, as necessary). However, estimation of every unknown parameter requires measurement of 13 of the 15 states. A table showing which observations are sufficient for the estimation of each parameter is constructed. The approach taken to analyze the NF-κB model is computationally efficient and can be readily applied to different sets or subsets of the unknown parameters. Procedure for Testing the Identifiability of the NF-KB Model. In this work the NF-κB model is analyzed using the proposed approach. The proposed approach is, however, augmented by first partitioning the full system in to a set of subsystems whose identifiability can be checked independently. By analyzing these subsystems, the number of possible unidentifiable parameters can be reduced. It has previously been shown that partitioning a system into a set of independent subsystems can be used to reduce the number of parameters which have to be tested for identifiability.39 Once the number of parameters which have to be tested for identifiability is reduced, the proposed approach is applied. The system equations are put in the form of eq 2 and then the rank condition in eq 3 is checked. This procedure is applied to the system for two cases. In the first case, the system states are assumed to be directly measurable. In the second case, the states are assumed to be measured via observation gains w1 to w15, themselves unknown parameters, so that yi ) wixi for i ) 1, 2, · · · , 15. It is shown that the observation gains can be estimated independently of the other unknown parameters. This implies that the system with unknown observation gains is also identifiable. Finally, a table is constructed that allows experimenters to decide which observations are critical for the identification of each parameter. A list of unnecessary, necessary, and sufficient observations can be constructed for each parameter or parameter set to be identified. This information may be critical for systems where obtaining experimental data is either costly or time-consuming. The NF-KB Model Equations. The equations for the NF-κB regulation network are given by18 x˙ ) f(x, u)

(13)

6129

Table 1. Definition of State Variables state

description

x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15

cytoplasmic concentration of IKK kinase in neutral state cytoplasmic concentration of IKK kinase in active state cytoplasmic concentration of IKK kinase in inactive state concentration of (IKKa | IκBR) complexes concentration of (IKKa | IκBR | NFκB) complexes concentration of free cytoplasmic NF-κB protein concentration of free nuclear NF-κB protein concentration of protein A20 concentration of A20 transcript concentration of free cytoplasmic IκBR protein concentration of free nuclear IκBR protein concentration of IκBR transcript cytoplasmic concentration of (IκBR | NFκB) complexes nuclear concentration of (IκBR | NFκB) complexes concentration control gene transcripts

where f ) [f1, · · · , f15]T is given by f1 ) kprod - kdegx1 - k1ux1

(14)

f2 ) k1ux1 - k3x2 - k2ux2x8 - kdegx2 - a2x2x10 + t1x4 a3x2x13 + t2x5 (15) f3 ) k3x3 + k2ux2x8

(16)

f4 ) a2x2x10 - t1x4

(17)

f5 ) a3x2x13 - t2x5

(18)

f6 ) c6ax13 - a1x6x10 + t2x5 - i1x6

(19)

f7 ) i1kVx6 - a1x7x11

(20)

f8 ) c4x9 - c5x8

(21)

f9 ) c2 + c1x7 - c3x9

(22)

f10 ) -a2x2x10 - a1x6x10 + c4ax12 - c5ax10 - i1ax10 + e1ax11 (23) f11 ) -a1x7x11 + i1akVx10 - e1ax14

(24)

f12 ) c2a + c1ax7 - c3ax12

(25)

f13 ) a1x6x10 - c6ax13 - a3x2x13 + e2ax14

(26)

f14 ) a1x7x11 - e2akVx14

(27)

f15 ) c2c + c1cx7 - c3cx15

(28)

where the states x1 to x15 are given in Table 1. A detailed description of how the model equations are derived is available in Lipniacki et al.35 The parameters in the NF-κB model are divided into two sets: P1 ) {kprod, kdeg, k1, k3, k2, e2a, i1, i1a, c5, c4a} P2 ) {a1, a2, a3, t1, t2, c1a, c2a, c3a, c5a, c6a, c1 c2, c3, c4, kV, e1a, c1c, c2c, c3c} As was done in Fujarewicz et al.,18 the parameters in the set P1 are to be estimated from experimental data, while the parameters in P2 are to be treated as known. Nominal values for all parameters in P1 and P2 are listed in Lipniacki et al.35 The input u in the NF-κB system takes on values in the set {0, 1} with 1 and 0 implying the external stimulation of the signaling pathway, and lack thereof, respectively.

6130


Identification of Subsystems. The NF-κB system can be divided into three subsystems. Each subsystem is a subset equations which can be analyzed independently. In this section, the information that may be obtained from these subsystems is quantified. First, the case where x1 is known exactly is discussed. Second, the case where only a gain-dependent measurement yi ) wixi for some unknown gain wi ∈ R, with i ∈ {1, 2, · · · , 15} is discussed. The first subsystem is given by eq 14 when u ) 0. Note that the time-solution of x1 does not depend on any of the other states in the model. If u ) 0 then measurement of x1, under idealized conditions, can be used to determine both kprod and kdeg. This is because the derivative x˙1, when plotted against x1, will form a linear graph with an intercept of kprod and a slope of -kdeg. This approach for computing kprod and kdeg is not useful in practice, because it requires the estimation of x˙1 from noisy data. However, the fact that kprod and kdeg can be identified under ideal conditions implies that each parameter has a unique contribution to the dynamics of x1, and thus neither parameter can be removed outright from the model. Furthermore, this analysis shows that there exists an experimental design where measurements of x1 alone are sufficient to estimate kprod and kdeg. The second subsystem is eq 14 when u ) 1. Plotting the derivative x˙1 versus x1 for this subsystem will yield a graph with a slope of -(kdeg + k1). As a result, using information from both subsystems (i.e., u ) 0 and u ) 1) each of kdeg, kprod, and k1 can be estimated. If only a gain-dependent measurement y1 ) w1x1 is available, then the two subsystems become y˙1 ) kprodw1 - kdegy1 and y˙1 ) kprodw1 - (kdeg + k1)y1 respectively, and, as a result, only kdeg and k1 can be uniquely determined from the two subsystems. The third proper subsystem of the NF-κB system is given by eq 16 when u ) 0. In this case, a measurement of x3 can be used to compute k3. If y3 ) x3w3 is measured, the subsystem becomes y˙3 ) k3y3 which implies that k3 can be estimated regardless of the value of w3. Using the three subsystems identified in this section, it has been shown that four parameters kprod, k1, kdeg, and k3 can be identified if exact values of x1 and x3 are known as a function of time. Furthermore, if gain-dependent measurements are available, then three of the parameters, k1, kdeg, and k3, are identifiable. As a result, one can reduce the set of possibly unidentifiable parameters in the NF-κB system to the set P3 ) P1 - {k1, kdeg, k3} ) {kprod, k2, e2a, i1, i1a, c5, c4a} Identifiability Analysis of the Main System. Using the method proposed in this work, the identifiability of the NF-κB system can be checked in a computationally efficient manner using linear algebra. Specifically, the equations that make up NF-κB are linear in the elements of the unknown parameters set P3. This implies that the function f(x, u) can be rewritten as f(x, u) ) A(ξ)θ + b(x, u)

where θ ) [kprod, k2, e2a, i1, i1a, c5, c4a]T, ξ ) [x˙, x, u], and the entries in both the matrix A(ξ) ) (∂f(x, u)/∂θ)|ξ)[x˙, x, u] and the vector b ) f(x, u) - (∂f(x, u)/∂θ)θ are not functions of the unknown parameters. Equation 13 can be rewritten A(ξ)θ ) g(ξ)

(29)

where g(ξ) ) x˙ - b(x, u), and the matrix A(ξ) can be rowreduced by premultiplying eq 29 with a matrix, Q, of elementary row operations giving

where

[

QA(ξ)θ ) Q(g(ξ))

1 0 0 0 0 0 0 0 0 0 0 0 ux2x8 0 -x6 0 0 0 0 0 0 0 0 0 0 0 x14 0 x12 0 0 0 0 -x10 0 QA(ξ) ) -x 0 0 0 0 0 0 8 0 0 0 0 0 0 kVx12 0 0 0 0 0 0 0 l l l l l l l 0 0 0 0 0 0 0 and

[

]

(30)

ξ)[x˙,x,u]

kdegx1 + k1ux1 - x˙1 -k3x3 - x˙3 -c6ax13 + a1x6x10 - t2x5 - x˙6 Qg(ξ) ) -a1x6x10 + t6ax13 + a3x2x13 - x˙13 a2x2x10 + a1x6x10 + c5ax10 - e1ax11 - x˙10 -c4x9 - x˙8 l

]

ξ)[x˙,x,u]

The identifiability of the NF-κB system can be evaluated using eq 9 with the time of sampling chosen so that none of u(t), x2(t), x8(t), x14(t), x10(t), and x12(t) are identically zero. Using the notation of eq 9 with ξ )[x˙(t), x(t), u(t)], A˜(ξ) ) QA(ξ) and Φ ) θ. In this case,

(

rank

)

∂A˜(ξ)oΦ ) rank(QA(ξ)) ) 7 ∂θ

and the system is identifiable. It is important to note that the procedure proposed in this work is not for parameter estimation. Rather it is for deciding whether parameters can, even in theory, be identified from experimental data. More specifically, the analysis proposed in this work allows the determination of which states should be measured in order to identify specific parameters. For example, the first row in eq 30 can be written as kprod ) kdegx1 + k1ux1 - x˙1 which implies that given knowledge of x1(t) * 0, as well as x˙1(t) and u(t), one may in theory compute kprod. Of course, typically, one would not rely on a single observation nor on an estimation of the time-derivative of the states to estimate the system parameters. However, the parameter kprod is structurally identifiable if one is able to measure the state x1. Similarly, k2, e2a, i1, c5, c4a appear alone in the second, third, fourth, sixth, and seventh lines of eq 30, respectively. As a result, k2 can be determined, for u * 0, if x2, x8, and x3 are measured, e2a can be determined if x5, x6, x10, and x13 are measured, i1 can be


6131

Identifiability Analysis of the Main System with Unknown Measurement Gain. The relationship between measured values (e.g, using blotting technique) and concentrations (i.e., state values) is typically nonlinear. However, over a short-range of concentrations, this relationship may be approximated by the linear relation yi ) wixi, where wi for i ) 1, 2, ..., 15 are unknown parameters.18 Generally, the magnitude of the these measurement gains (i.e., wi) depends on the technique that is used. Specifically, in the case of Western Blots or electrophoretic mobility shift assays, the gains are typically normalized within a blot so that the maximum response is scaled to unity.18 Letting θi )(wi)-1 for i ) 1, 2, ..., 15, System 13 can be rewritten as ˆf(y˙, y, u) ) 0 with ˆf ) [fˆ1, ..., ˆf15], where ˆfi for i ) 1, 2, ..., 15 are given by ˆf1 ) -y˙1θ1 + kprod - kdegθ1y1 - k1uθ1y1

(31)

ˆf2 ) -y˙2θ2 + k1uθ1y1 - k3θ2y2 - k2uθ2θ8y2y8 - kdegθ2y2 a2θ2θ10y2y10 + t1θ4y4 - a3θ2θ13y2y13 + t2θ5y5 (32)

Figure 1. A matrix showing the dependence between the states and unknown parameters.

determined if x2, x14, x6, x10, and x13 are measured, c5 can be determined if x9 and x8 are measured, and c4a if x12, x2, x6, x7, x10 and x11 are measured. Finally, i1a only appears in the fifth row of eq 30, and can be estimated if c4a can be estimated and x11 is measured. The situation is summarized in Figure 1. The dark squares in the ijth entry of the table in Figure 1 indicate that the state xi can be used to estimate the parameter in column j. Figure 1 can be used to determine which states are the most critical for parameter estimation. The rows corresponding to states x6 and x10 both have four dark squares indicating that they can be used for estimating four different parameters. The states x4 and x15 are the least useful states for parameter estimation because neither is necessary for the estimation of any parameter. As a result, the system would be identifiable even if states x4 and x15 were not measured. That is, even though the assumption of the output gain matrix C(θ) being invertible was used in the algorithm, the results obtained for this system suggest that not every state need be measured for identifiability. Another use of Figure 1 is that it can be used to determine which parameters are identifiable given a set of measured states. For example, if only x2, x3, x8, and x9 are measured, then only the parameters k2 and c5 can be uniquely identified. As a second example consider the set of measured variables reported by Lipniacki et al.35 In this case the measured variables are the free nuclear NF-κB concentration (x7), total cytoplasmic IKBR concentration (x10 + x13), A20 transcript concentration (x9), total IKK cytoplasmic concentration (x1 + x2 + x3), and active IKK concentration (x2). According to Figure 1, even if all of these variables are measured exactly and sufficiently often, it is not possible to guarantee a unique estimate of any of the model parameters. That is, no column corresponding to model parameters in Figure 1 has dark squares only in rows corresponding to observed states. For example, for the identification of parameter k2, in the second column, one would need to measure both x2 and x3. However, as only the sum x1 + x3 is known, it is not possible to compute x3 independently of x1 and therefore k2 is not identifiable in this scenario.

ˆf3 ) -y˙3θ3 + k3θ3y3 + k2uθ2θ8y2y8

(33)

ˆf4 ) -y˙4θ4 + a2θ2θ10y2y10 - t1θ4y4

(34)

ˆf5 ) -y˙5θ5 + a3θ2θ13y2y13 - t2θ5y5

(35)

ˆf6 ) -y˙6θ6 + c6aθ13y13 - a1θ6θ10y6y10 + t2θ5y5 - i1θ6y6 (36) ˆf7 ) -y˙7θ7 + i1kVθ6y6 - a1θ7θ11y7y11

(37)

ˆf8 ) -y˙8θ8 + c4θ9y9 - c5θ8y8

(38)

ˆf9 ) -y˙9θ9 + c2 + c1θ7y7 - c3θ9y9

(39)

ˆf10 ) -y˙10θ10 - a2θ2θ10y2y10 - a1θ6θ10y6y10 + c4aθ12y12 c5aθ10y10 - i1aθ10y10 + e1aθ11y11 (40) ˆf11 ) -y˙11θ11 - a1θ7θ11y7y11 - i1akVθ10y10 - e1aθ14y14 (41) ˆf12 ) -y˙12θ12 + c2a + c1aθ7y7 - c3aθ12y12

(42)

ˆf13 ) -y˙13θ13 + a1θ6θ10y6y10 - c6aθ13y13 - a3θ2θ13y2y13 + e2aθ14y14 (43) ˆf14 ) -y˙14θ14 + a1θ7θ11y7y11 - e2akVθ14y14

(44)

ˆf15 ) -y˙15θ15 + c2c + c1cθ7y7 - c3cθ15y15

(45)

The set of unknown parameters for this case is given by P4 ) {kprod, k2, e2a, i1, i1a, c5, c4a, θ1, θ2, ..., θ15} While eqs 31-45 are not linear in the unknown parameters, they are polynomial functions of the unknown parameters. As a result, they can be treated as linear functions of products of the elements in P4 (i.e., monomials in the unknown parameters). The set of (possibly dependent) monomials in the elements of P4 which appear in eqs 31-45 is given by Pm ) {kprod, θ1, ..., θ15, θ2θ8, θ7θ11, θ10θ2, θ13θ2, θ14e2a, θ8c5, θ10i1a, θ12c4a, θ6i1, θ10θ6} Defining

6132


θˆ 1 ) kprod θˆ 2 ) θ2θ8k2 θˆ ) θ θ

θˆ 7 ) θ8c5 θˆ 8 ) θ10i1a θˆ ) θ c

θˆ 4 ) θ10θ2 θˆ 5 ) θ13θ2 θˆ ) θ e

θˆ 10 ) θ6i1 θˆ 11 ) θ10θ6

3

6

7 11

9

ˆ in Equation 46 Scheme 1. The Row Reduce Matrix Q

12 4a

14 2a

and redefining θ ) [θˆ , θ1, ..., θ15], eqs 31-45 can be rewritten in the form of eq 29 with A(ξ) ) (∂fˆ/∂p)(y˙, y, u) and g(ξ) ) F(y˙, y, u) - A(ξ)θ both independent of any unknown parameters. As in the previous section, the matrix A can be row-reduced ˆ (shown in Scheme 1) and the set of 15 using a matrix Q equations is given by ˆ A(ξ)θ ) Q ˆ q(ξ) Q

(46)

Note that the rank of A in this case is at most 15, which is less than the number of unknown parameters, and as a result, one cannot use eq 46 to deduce identifiability. However, by ˆ A(ξ) and b(ξ) at different values in time, one evaluating Q can obtain additional constraints on the set of possible parameters. For example, the last row of eq 46, given by -c1c(-c3ay12 + y˙12) c1c w12 + (-c3cy15 + y˙15)w15 ) c - c2c c1a c2a 1a (47) -1

-1

contains two unknown parameters θ12 ) w12 and θ15 ) w15 . If one chooses to evaluate eq 47 at two different times t1 and t2, the general form given by eq 8 will contain the following two rows: -c1c(-c3ay12(t1) + y˙12(t1)) w12 + (-c3cy15(t1) + y˙15(t1))w15 ) c1a c1c c - c2c c2a 1a -c1c(-c3ay12(t2) + y˙12(t2)) w12 + (-c3cy15(t2) + y˙15(t2))w15 ) c1a (-c3cy15(t2) + y˙15(t2)) The gains w12 and w15, can be computed by choosing t1 and t2 so that the following condition holds:

[

]

-c1c(-c3ay12(t1) + y˙12(t1)) (-c3cy15(t1) + y˙15(t1)) c1a *0 det -c1c(-c3ay12(t2) + y˙12(t2)) (-c3cy15(t2) + y˙15(t2)) c1a (48) The determinant function in eq 48 is a algebraic relation in the variables y12, y˙12, y15 and y˙15. If one assumes that there are no implied algebraic constraints in the output variables and their time-derivatives then the determinant in eq 48 will not, in general, be identically zero as a function of time. As a result, it is possible to choose times t1 and t2 so that θ12 and θ15 are identifiable. Using a similar argument, the secondlast (i.e., 14th) row of eq 46 can be used to show that θ9 and θ12 are identifiable for an appropriate choice of sampling times. This process can be repeated with each row of eq 46 to show that each of the coefficients θ1-θ15 can be estimated independently of the unknown parameters. A list of the

Table 2. List of Estimated Coefficients and Known Quantities Necessary for Estimation for the Last Seven Rows of eq 46 row

estimated coefficients

known quantities used for estimation

15 14 13 12 11 10 9

θ12, θ15 θ9, θ12 θ7, θ9 θ5, θ6, θ13, θ14 θ1, θ2, θ3, θ4 θ10, θ11, θ13, θ14 θ8

y12, y15 y9, y12 y7, y9 θ7, y5, y6, y13, y14 θ5, y1, y2, y3, y4 θ5, θ6, θ7, y13, y14 θ3, y2, y3, y8

parameters estimated and outputs used in each row is summarized in Table 2. The fact that the system gains θ1-θ15 (i.e., the observation gains) are all identifiable independently of the unknown parameter set P3 (i.e., the system parameters considered in the case where y ) x) implies that all parameters in the system, including the observation gains w, are identifiable. The calculations involved in the identifiability analysis presented in this section are not computationally expensive. If one were to try and use a linearization-based approach, the calculations would be more complex. This is because one would have to linearize the system about some parameterdependent steady state value. This steady-state value can be computed by setting ˆf(y, u) ) 0 and solving for y. However, the mapping ˆf is polynomial, and, as a result, there may be several steady state solutions for y. To use the linearizationbased approach one would have to determine which steady state solutions for y are physically realizable (for some nominal parameter value) and, if more than one steady state solution is realizable, one would have to consider several discrete scenarios. The method proposed by Ljung and Glad27 could be applied to this system. However, this approach requires the implementation of Ritt’s algorithm which is generally more computationally intensive than row reduction. Indeed, in their work, Ljung and Glad state “Honestly, it must be added that the actual complexity of this algorithm increases very rapidly with the size of the problem. Our


current implementations have difficulties when the sum of the number of parameters and the number of states exceeds 10 or so.” The authors go on to state that “The identification question can be answered without a complete characteristic set. It is an important challenge to find practical implementations that are more efficient than our current ones.” An alternative to the linearization or differential-algebra based approaches is the isomorphism-based approach.9 However, it has been shown that the isomorphism-based approach is not tractable for this system.18 Finally, one may employ a series-based identifiability test.40 However, this approach would require the manipulation of equations describing the time-derivatives of the model equations which contain many terms because ˆf is polynomial in the system states. Using the proposed framework, it has been shown that the NF-κB model is identifiable from some input and set of initial conditions. However, this does not guarantee that all parameters can be estimated using any initial condition or any data set. As a result, the analysis presented in this work is a good “first step” toward parameter estimation and experimental design because it allows an experimenter to determine whether any parameters can be removed outright from the model. Furthermore, it allows specific outputs and parameters to be paired for estimation. The proposed approach however, does not address the issue of how many or what type of experiments should be performed to accurately estimate all parameters. Conclusions Mathematical models of physical systems often have more parameters than can be estimated from experimental observations. A system is called identifiable if there exist some set of experiments which can, under ideal conditions, be used to uniquely estimate each parameter value. Local identifiability implies that the parameters are identifiable in a neighborhood of some nominal value. Obtaining necessary and sufficient conditions for identifiability of a given model is often computationally intensive. In this work, a computationally efficient approach for testing nonlinear LTI ODE and index-1 DAE polynomial models for identifiability is proposed. The approach consists of two steps. First, the model equations are posed as a set of linear equations in monomials of the system parameters. Second, the linear model is tested for identifiability using a rank condition. The proposed approach allows a table to be constructed which can be used to identify which process variables must be measured in order to estimate only a subset of the unknown parameters. As an illustrative example, the identifiability of a NF-κB system model is analyzed. It is shown that the model is identifiable if all states are measured. A subset of states is identified whose measurement can guarantee identifiability. A table is constructed to illustrate the relationship between measured variables and individual parameter identifiability. Direct observation of process states in the NF-κB system is not generally possible. The relationship between the measured values and state variables is approximated using a set of observation gains which must be experimentally determined. In this work, it is shown that these observation gains can be estimated independently of the other unknown system parameters, and as a result, the system (including the observation gains) is identifiable. The proposed approach provides sufficient conditions for ODE and index-one DAE systems whose equations can be expressed as polynomials in the outputs and unknown

6133

parameters. The approach is only effective on a dense subset of the parameter space (e.g., all of the output gains cannot be identically zero), and assumes that the model states are related to the measured outputs by a matrix. The main advantage of the proposed approach is that it is computationally tractable for systems which may be difficult to test for identifiability using other means. While the approach does not identify an ideal or optimal set of experiments for parameter estimation, it is a first step toward a proper experimental design which can be used for parameter estimation. Acknowledgment The author acknowledges the financial support of the Canadian Natural Science and Engineering Research Council. Literature Cited (1) Walter, E.; Lecourtier, Y. Math. Biosci. 1981, 56, 1–25. (2) Xia, X.; Moog, C. H. IEEE Trans. Autom. Control 2003, 48, 330– 336. (3) Jimenez-Hornero, J. E.; Santos-Duenas, I. M.; Garcia Garcia, I. Math. Biosci. 2008, 116, 154–162. (4) Ben-Zvi, A.; McLellan, P. J.; McAuley, K. B. AIChE J. 2004, 50, 2493–2501. (5) Deng, Y.; Chen, J. Water Sci. Technol. 2006, 53, 93–99. (6) Nikerel, I. E.; van Winden, W. A.; Verheijen, J. T.; Heijnen, J. J. Metab. Eng. 2009, 11, 20–30. (7) Denis-Vidal, L.; Jauberthie, C.; Joly-Blanchard, G. IEEE Trans. Autom. Control 2006, 51, 154–158. (8) Boens, N.; Ameloot, M. Int. J. Quantum Chem. 2006, 106, 300– 315. (9) Vajda, S.; Godfrey, K. R.; Rabitz, H. Math. Biosci. 1989, 93, 217– 248. (10) Peeters, R. L. M.; Hanzon, B. Automatica 2005, 41, 513–529. (11) Pohjanpalo, H. Math. Biosci. 1978, 41, 21–33. (12) Lecourtier, Y.; Lamnabhi-Lagarrigue, F.; Walter, E. Volterra and generating power series approach to identifiability testing. In Identifiability of Parametric Models; Walter, E., Ed.; Pergamon: Oxford, 1987; pp 5066. (13) Tunali, E. T.; Tarn, T. IEEE Trans. Autom. Control 1987, AC-32, 146–154. (14) Grewal, M.; Glover, K. IEEE Trans. Autom. Control 1976, AC21, 833–837. (15) Vajda, S. IEEE Trans. Autom. Control 1987, AC-32, 182–184. (16) Saccomani, M. P.; Audoly, S.; D’Angio´, L. A new differential algebra algorithm to test the identifiability of nonlinear systems with given initial conditions. Proc. 40th Conf. Dec. Control 2001, 3108–3113. (17) Ben-Zvi, A.; McLellan, P. J.; McAuley, K. B. Can. J. Chem. Eng. 2006, 84, 590–596. (18) Fujarewicz, K.; Kimmel, M.; Lipniacki, T.; Sierniak, A. IEEE/ACM Trans. Comput. Biol. Bioinform. 2007, 4, 322–335. (19) Lagrange, S.; Delanoue, N.; Jaulin, L. Automatica 2008, 44, 2959– 2962. (20) Monod, J. Annu. ReV. Microbiol. 1949, 3, 371–394. (21) Marcos, N. I.; Guay, M.; Dochain, D.; Zhang, T. J. Process Control 2004, 14, 217–328. (22) Caperon, J.; Meyer, J. Deep Sea Res. Oceanogr. 1972, 19b, 619– 632. (23) Lehman, J. T.; Botkin, D. B.; Likens, G. E. Limnol. Oceanogr. 1975, 20, 343–378. (24) Bequette, B. W. Process Dynamics Modeling, Analysis, and Simulation; Prentice Hall: Upper Saddle River, NJ, 1998. (25) Ben-Zvi, A.; Vernon, S.; G., B. PLoS Comput. Biol. 2009, 5, 1– 10. (26) Perelson, A. S. Nat. ReV. Immunol. 2002, 2, 28–36. (27) Ljung, L.; Glad, T. Automatica 1994, 30, 265–276. (28) Ljung, L. System Identification, Theory for the User, 2nd ed.; Prentice Hall: Upper Saddle River, NJ, 1999. (29) Bellman, R.; Astrom, K. J. Math. Biosci. 1970, 7, 329–339. (30) Walter, E.; Pronzato, L. Identifiabilities and nonlinearities in nonlinear systems. Model. Estimat. 1995, 111–143. (31) Audoly, S.; Bellu, G.; D’Angio´, L.; Saccomani, M. P.; Cobelli, C. IEEE Trans. Biomed. Eng. 2001, 48, 55–65.

6134


(32) Ritt, J. F. Differential Algebra; American Mathematical Society: New York, 1950. (33) Hoffman, K.; Kunze, R. Linear Algebra, 2nd ed.; Prentice-Hall: New York, 1971. (34) Chicone, C. Ordinary Differential Equations with Applications; Springer-Verlag: New York, 1999. (35) Lipniacki, T.; Paszek, P.; Brasier, A. R.; Luxon, B.; Kimmel, M. J. Theor. Biol. 2004, 228, 195–215. (36) Bard, Y. Nonlinear Parameter Estimation; Academic Press: New York, 1974. (37) Ihekwaba, A.; Broomhead, D.; Grimley, R.; Benson, N.; White, M.; Kell, D. Syst. Biol. 2005, 152, 153–160.

(38) Cho, K. H.; Kolch, W.; Wolkenhauer, O. Simulation 2003, 79, 726– 739. (39) Jayasankar, B.; Ben-Zvi, A.; Huang, B. Comput. Chem. Eng. 2009, 33, 484–492. (40) Walter, E. Lecture Notes in Biomathematics; Springer-Verlag: Berlin, 1982.

ReceiVed for reView November 23, 2009 ReVised manuscript receiVed March 18, 2010 Accepted April 26, 2010 IE9018512

A Computationally Efficient Algorithm for Testing ... - ACS Publications

Recommend Documents