ON IDENTIFIABILITY OF NONLINEAR ODE MODELS AND APPLICATIONS IN VIRAL DYNAMICS

doi:10.1137/090757009

SIAM Rev Soc Ind Appl Math. Author manuscript; available in PMC 2011 July 20.

Published in final edited form as:

SIAM Rev Soc Ind Appl Math. 2011 January 1; 53(1): 3–39.

doi: 10.1137/090757009

PMCID: PMC3140286

NIHMSID: NIHMS293681

ON IDENTIFIABILITY OF NONLINEAR ODE MODELS AND APPLICATIONS IN VIRAL DYNAMICS

HONGYU MIAO,^* XIAOHUA XIA,^† ALAN S. PERELSON,^‡ and HULIN WU^§

^* Department of Biostatistics and Computational Biology, University of Rochester School of Medicine and Dentistry, 601 Elmwood Avenue, Box 630, Rochester, New York 14642, USA

^† Department of Electrical, Electronic and Computer Engineering, University of Pretoria, Lynnwood Road, Pretoria 0002, South Africa

^‡ Theoretical Biology and Biophysics Group, MS-K710, Los Alamos National Laboratory, Los Alamos, NM 87545, USA

^§ Department of Biostatistics and Computational Biology, University of Rochester School of Medicine and Dentistry, 601 Elmwood Avenue, Box 630, Rochester, New York 14642, USA

HONGYU MIAO: hongyu_miao/at/urmc.rochester.edu; XIAOHUA XIA: xxia/at/postino.up.ac.za; ALAN S. PERELSON: asp/at/lanl.gov; HULIN WU: hwu/at/bst.rochester.edu

See other articles in PMC that cite the published article.

Abstract

Ordinary differential equations (ODE) are a powerful tool for modeling dynamic processes with wide applications in a variety of scientific fields. Over the last 2 decades, ODEs have also emerged as a prevailing tool in various biomedical research fields, especially in infectious disease modeling. In practice, it is important and necessary to determine unknown parameters in ODE models based on experimental data. Identifiability analysis is the first step in determing unknown parameters in ODE models and such analysis techniques for nonlinear ODE models are still under development. In this article, we review identifiability analysis methodologies for nonlinear ODE models developed in the past one to two decades, including structural identifiability analysis, practical identifiability analysis and sensitivity-based identifiability analysis. Some advanced topics and ongoing research are also briefly reviewed. Finally, some examples from modeling viral dynamics of HIV, influenza and hepatitis viruses are given to illustrate how to apply these identifiability analysis methods in practice.

Keywords: ODE modeling, structural identifiability, practical identifiability, sensitivity-based identifiability, viral dynamics

1. Introduction

Ordinary differential equation (ODE) models have been widely used to model physical phenomena, engineering systems, economic behavior, and biomedical processes. In particular, ODE models have recently played a prominent role in describing both the within host dynamics and epidemics of infectious diseases and other complex biomedical processes (e.g., [2, 15, 59, 74, 75, 77]). Great attention has been paid to the so-called forward problem or simulation problem, i.e., predicting and simulating the results of measurements or output variables for a given system with given parameters. However, less effort has been devoted to the inverse problem, i.e., using the measurements of some state or output variables to estimate the parameters that characterize the system, especially for nonlinear ODE models without closed-form solutions.

In reality, before rigorous parameter estimation methods can be applied to an ODE model to formally estimate the model parameters based on experimental data, a serious barrier to overcome is how to verify whether the model parameters are identifiable based on the measurements of output variables or their functions when the ODE model does not have a closed-form solution. Further questions include, if not all model parameters are identifiable, are a subset of parameters identifiable? How many measurements, at which time points, are necessary to identify the identifiable parameters? To answer these questions, identifiability analysis should be done before tackling the inverse problem.

The literature on ODE identifiability analysis is found in journals from a variety of scientific fields such as mathematics, biomedical modeling, engineering, and statistics; in addition, various techniques and methodologies from these disciplines are employed in ODE identifiability studies. Therefore, it is useful to have a comprehensive review on these methods and approaches, and their further applications, e.g., in experimental design [32, 70, 71, 96]. In this paper, we review identifiability methods with a focus on nonlinear ODE models for which close-form solutions are not available. In Section 2, we review various definitions related to identifiability analysis. We review various techniques for structural identifiability analysis in Section 3. In addition to theoretical (structural) identifiability, it is also important to evaluate practical identifiability when experimental data are contaminated with measurement noise. This will be reviewed in Section 4. Sensitivity analysis is widely used in mathematical modeling to evaluate how sensitive output variables are to parameter values and input variables. Some sensitivity analysis techniques can also be used to evaluate parameter identifiability in ODE models, as will be reviewed in Section 5. We illustrate identifiability techniques using examples from infectious disease modeling in Section 6. We conclude this paper with some discussions and summary in Section 7.

2. Definitions

A general dynamic system can be expressed as follows:

(2.1)

(2.2)

where x(t) [set membership]

R^m is a vector of state variables (or dependent variables), y(t) [set membership]

R^d the measurement or output vector, u(t) [set membership]

R^p the known system input vector, and θ [set membership]

R^q the parameter vector. The system given by Eq. (2.1) is an ordinary differential equation model (ODE model). For the inverse problems, θ is unknown and has to be estimated based on experimental data. There are three common scenarios for θ:

constant parameters only;
time-varying parameters only;
a mixture of both constant and time-varying parameters.

Let θ = (θ_c, θ_t), where θ_c denotes the constant unknown parameters and θ_t denotes the time-varying unknown parameters. Now Eq. (2.1) and (2.2) can be re-written in the form:

(2.3)

(2.4)

Before we introduce definitions of identifiability, we review three important and useful concepts developed in control theory: reachable, controllable and observable [60, 111].

Definition 2.1

Reachable: For a certain initial state x₀ of interest, a state x₁ is said to be reachable if there exists a trajectory x(t) starting from x₀ which can achieve x₁ in a finite time given an admissible system input u(t).

Note that in the definition above, u(t) is called an admissible input (or admissible control) if it satisfies all system constraints at any time of interest and a solution of the dynamic system exists. The existence of such an admissible u(t) leads to the definition of controllability.

Definition 2.2

Controllable: If there exists an admissible u(t) which can transfer an initial state of interest to a target state in a finite time, the dynamic system is said to be controllable.

Controllability is an important property of a dynamic system since it indicates whether a system will respond to a certain input and behave as expected. One important application of this concept is stabilization of dynamic systems. In biomedical research, the dynamic system could be a virus, such as HIV, infecting a human and the system input could be antiretroviral therapy; how to control or stabilize the virus via intentionally-designed intervention strategies is still an interesting and challenging research topic [1, 22, 55, 58, 80, 81, 126].

To better understand dynamic system structure and behavior, it is also necessary to obtain measurements of the system output variables. However, we may not be able to directly measure the state variables; instead, we may be able to measure output variables which are functions of system input and state variables, as specified in Eq. (2.2) or (2.4). If it is possible to determine the system state from system output measurements, the system is observable according to the following definition:

Definition 2.3

Observable: Given an initial state x₀ and an admissible control u(t), if the current system state x(t) can be determined only through the system output y(t) in a finite time, the system is said to be observable.

In the definition above, it is usually assumed that the system output y(t) can be measured without error. The three definitions introduced so far describe the relationships among four basic factors of a dynamic system: initial state, input, current state and output. There are also other interesting concepts and definitions originating in control theory for connecting these four basic factors of dynamic systems (e.g., [23, 38, 50, 111]).

The concepts discussed above have been introduced for systems with known parameters. However, these concepts, especially controllability and observability, are also directly related to system (parameter) identifiability. A system which is controllable and observable has strong connections among input, state and output variables, and such strong connections may indicate that the system is identifiable. The reader is referred to [4, 19, 20, 98, 99, 101, 110, 123] for further discussions.

Definition 2.4

Identifiable: The dynamic system given by Eqs. (2.1) and (2.2) is identifiable if θ can be uniquely determined from the given system input u(t) and the measurable system output y(t); otherwise, it is said to be unidentifiable.

Limited system inputs may not result in system outputs with sufficient information for uniquely determining system parameters. Thus, it is also necessary to have informative input signals in order to identify system parameters. This idea was formalized by introducing the concept of a persistently exciting input [63, 64, 65]. Simply speaking, an input is said to be persistently exciting if enough information on the output variables can be generated from the input to identify system parameters in the sense that all estimates of system parameters converge to their true values in a finite time [47]. The assumption of persistently exciting inputs is a prerequisite for structural identifiability analysis [10] as will be discussed in the next section.

Ljung and Glad [64] introduced two important concepts, globally identifiable and locally identifiable.

Definition 2.5

Globally identifiable: A system structure is said to be globally identifiable if for any admissible input u(t) and any two parameter vectors θ₁ and θ₂ in the parameter space Θ, y(u, θ₁) = y(u, θ₂) holds if and only if θ₁ = θ₂.

Definition 2.6

Locally identifiable: A system structure is said to be locally identifiable if for any θ within an open neighborhood of some point θ_* in the parameter space, y(u, θ₁) = y(u, θ₂) holds if and only if θ₁ = θ₂.

Both definitions use the concept of one-to-one mapping between parameters and system input/output. With the development of various identifiability analysis techniques, more specific definitions of identifiability have been introduced by a number of authors [7, 10, 21, 56, 64, 114, 119]. For instance, Tunali and Tarn [104] introduced a definition of identifiability when an initial state is given, which was termed locally strongly identifiable. A similar concept was introduced in [56], called x₀-identifiable.

Definition 2.7

Locally strongly identifiable (x₀-identifiable): For an admissible input u(t) in the time range of interest [t₀, t₁] and a given initial state x₀ = x(t₀) which is independent of θ and not an equilibrium point, if there exists an open set Θ⁰ within the parameter space Θ such that for any two different parameter vectors θ₁, θ₂ [set membership]

Θ⁰, the solutions x(t, θ, u) exists on [t₀, t₀ + ε] (t₀ < ε ≤ t₁ − t₀) for both θ₁ and θ₂, and y(t, θ₁, x₀, u(t)) ≠ y(t, θ₂, x₀, u(t)) on [t₀, t₀ + ε], the system structure is said to be locally strongly identifiable (or x₀-identifiable).

We notice that this definition is specific for differential equation systems. But it is stringent with respect to the initial state. More generally, Xia and Moog [119] introduced structurally identifiability as follows:

Definition 2.8

Structurally identifiable: Let equation M5

denote the function space expanded by all input functions on [t₀, t₁] which are differentiable up to the order N, and let M denote an open set of initial system states. The system structure is said to be structurally identifiable if there exist open and dense subsets M⁰ [subset or is implied by]

M, Θ⁰

Θ, and

such that the system is locally strongly identifiable at θ given u for any x₀ [set membership]

M⁰, θ

Θ⁰, and u

U⁰.

This definition is also interchangeably called geometrical identifiable [104, 56]. Besides these identifiability definitions based on one-to-one mappings between system parameters and system input-output, Glad [37] and Ljung and Glad [64] introduced a definition of identifiability based on the algebraic equations consisting of the system input and output, which was called algebraically identifiable. This definition is directly related to identifiability analysis techniques [29, 64, 118, 119]:

Definition 2.9

Algebraically identifiable: Based on algebraic equations of system state, input and output, if a meromorphic function

can be constructed after a finite number of steps of algebraic calculation or differentiation such that Φ = 0 and equation M8

hold in the time range of interest [t₀, t₁] for any (θ, x₀, u) in an open and dense subset of equation M9

, where k is a positive integer, , …, u^(k) the derivatives of u, and , …, y^(k) the derivatives of y, the system structure is said to be algebraically identifiable.

Similarly algebraically identifiability with known initial conditions can be defined as follows [119]:

Definition 2.10

Algebraically identifiable with known initial conditions: If a meromorphic function equation M10

, Φ

R^q, can be constructed from algebraic equations of system state, input and output after a finite number of steps of algebraic calculation or differentiation such that Φ = 0 and equation M11

hold for all ( equation M12

), where k is a positive integer, ( equation M13

) in an open and dense subset of Θ × M × U, ( equation M14

) and (

) are the derivatives of u and y at time equation M16

, respectively, the system structure is said to be algebraically identifiable with known initial conditions.

A number of studies have considered system identifiability given initial conditions [26, 37, 64, 83, 104, 119] and reported that known initial conditions can help to identify more parameters. Particularly, the work of Wu et al. [117] clarified that, giving initial conditions is equivalent to have more observations on system output such that parameter estimation reliability can be improved, especially for dynamic systems sensitive to initial conditions.

3. Structural Identifiability Analysis

In this section, we will review structural identifiability methods in details. We will also discuss the advantages and disadvantages of these methods in practical applications in order to help the practitioners choose the appropriate approach for specific problems. Furthermore, we will discuss the minimum number of observations obtained via structural identifiability analysis to uniquely determine all identifiable parameters, keeping in mind that this number could be much higher for real problems due to the presence of measurement error or model uncertainty.

The concept of structural identifiability was first introduced by Bellman and Åström [10]. As suggested by its name, the corresponding techniques verify system identifiability by exploring the system structure (that is, the model itself). Early structural identifiability analysis techniques were developed from control theories in 1970s for linear models, especially compartmental models. For instance, Bellman and Åström [10] proposed an analysis technique for linear ODE models based on Laplace transforms. Later, the method of power series expansion was proposed by Pohjanpalo [84] and the similarity transformation method was proposed by Walter and Lecourtier [113] for linear ODE models. These methods have been well summarized in [3] and [51]. In this paper, we focus on identifiability methods for nonlinear ODE models instead of linear models.

Some of the approaches for linear ODE models such as the similarity transformation method have been extended by Vajda and Rabitz [107], Vajda et al. [106] and Chappel and Godfrey [16] to nonlinear ODE models. However, the extension only works for a limited number of simple nonlinear problems [7]. For general nonlinear models, new techniques are needed. A simple approach for this purpose, called direct test, was proposed by Denis-Vidal and Joly-Blanchard [25] and Walter et al. [112]. The basic idea of this approach is to directly use the identifiability definition to verify parameter identifiability, either analytically [25] or numerically [112]. However, the analytical direct test is not suitable for high-dimensional problems and the numerical direct test also has some limitations due to the use of a cut-off value.

Under the framework of differential algebra [91], new methods and algorithms have also been developed to target identifiability of general nonlinear models [14, 64, 78]. The differential algebra approach can utilize the power of symbolic computations, which requires much less human intervention. Since the differential algebra method was introduced to investigate the structural identifiability problem [64, 78] in the early 1990s, it has been successfully applied to nonlinear differential equation models, including models with time-varying parameters [7]. Ljung and Glad [64] summarized three conditions under which the system structure is globally identifiable, locally identifiable or not identifiable, respectively; however, to verify the three conditions, rigorous mathematical theories need to be further developed.

Xia and Moog [119] proposed another method based on the Implicit Function Theorem. By taking derivatives of observable system outputs with respect to independent variables (e.g., time), all latent variables (unobservable system state variables) can be eliminated after algebraic calculations and a finite number of equations consisting of known system inputs, observable system outputs and unknown parameters can be formed. Then a matrix (called the identification matrix) consisting of the partial derivatives of these equations with respect to unknown parameters (and usually their derivatives with respect to independent variables) can be formed. If the identification matrix is non-singular, this system is identifiable. This method has the advantages of theoretical and practical simplicity and has been successfully applied to HIV dynamic models of up to 6-dimensions [70, 119]. However, this method requires high order derivatives, thus the matrix can easily become very complicated and the singularity of the matrix becomes difficult to verify. Wu et al. [117] further extended this method by considering multiple time points instead of high order derivatives to overcome the disadvantages of Xia’s method. The methods based on the Implicit Function Theorem can be applied alone to verify system identifiability and they can also be employed to verify the three conditions in the differential algebra approach. However, for dynamic models with time-varying parameters, the singularity of the identification matrix is difficult to evaluate and no reliable conclusion can be easily drawn. Therefore, in practice, the differential algebra approach and the Implicit Function Theorem approach may have to be combined to solve a problem. In addition, if initial conditions are unknown, the correlation between unknown initial conditions and other model parameters can not be verified by structural identifiability analysis unless such unknown initial conditions explicitly appear on the right-hand side of Eq. (2.1).

Before moving onto the technical details, it is also helpful to mention it here that structural identifiability analysis methods are not widely used in practice, yet, due to either the computational complexity or the lack of mature computer implementations.

3.1. Power Series Expansion and Similarity Transformation

Grewal and Glover [39] studied the identifiability problem for nonlinear ODE models by considering local linearization of nonlinear systems. However, “the linearized system being identifiable” is only a necessary condition for “the nonlinear system being identifiable” instead of a sufficient condition. Therefore, the local linear approximation cannot answer the question completely. Pohjanpalo [84] proposed another approach called power series expansion to better handle nonlinear problems.

For the power series expansion method, the function f in Eq. (2.3) is assumed to be infinitely differentiable with respect to time t, u and x in the time range of interest [t₀, t₁]; and the same assumption is needed for x, y and u with respect to time, and for h with respect to x. Such assumptions are necessary because the power series expansion may require derivatives of arbitrary orders. The nonlinear system considered by Pohjanpalo [84] is of the following form:

(3.1)

(3.2)

which is very restrictive. Consider the derivatives of system output y at time t₀

where k denotes the k^th derivative of y. Therefore, the system input and output can be connected by their derivatives with respect to time at t₀:

(3.3)

(3.4)

where k = 1, …, ∞. Since the derivatives of y are theoretically observable, they are considered as known. Therefore, an infinite number of equations can be obtained from Eq. (3.3) and (3.4) to solve for θ_c simultaneously. If the solution is unique, then the system structure is (locally) identifiable.

In nature, the power series expansion method is an approach to verify the x₀ identifiability (or local strong identifiability). However, this method has a serious drawback: high order derivatives are needed for a high dimensional problem and the resulting equations can easily become too complicated to solve. This disadvantage has prevented this method from becoming popular in practice.

Walter and Lecourtier [113] initially proposed the similarity transformation method for linear ODE models. The system concerned here is in the form:

(3.5)

(3.6)

where A, B and C are matrices of constant coefficients. The basic idea of this method is to find the similar matrix S = P⁻¹ AP of A such that

(3.7)

(3.8)

where P is a non-singular matrix. It is straight forward to show that if the only possible similar transformation of A is P = I, the system is uniquely and globally identifiable; if a finite number of P ≠ I can be found, the system is locally identifiable (or nonuniquely identifiable); otherwise, the system is unidentifiable.

Vajda et al. [106] and Vajda and Rabitz [107] extended the work of Walter and Lecourtier [113] and proposed the similarity transformation method to tackle the nonlinear ODE systems by making use of the local state isomorphism theorem [44, 50]. The nonlinear system considered in Vajda et al. [106] is of the following form:

(3.9)

(3.10)

Note that although this system is a single-input system, the conclusion based on this system can be generalized to multiple-input systems. It is necessary to introduce the definition of structural equivalence before we further introduce the similarity transformation method.

Definition 3.1

Structural equivalence: Given two systems of the family in Eq. (3.9) and (3.10), if there exist two parameters θ₁, θ₂ [set membership]

Θ such that for the same admissible input u(t), the solution of the two system exists for θ₁ and θ₂, respectively, and the corresponding system outputs are the same, the system with parameter θ₁ is said to be equivalent to the system with parameter θ₂, denoted by θ₁ ~ θ₂.

Under the similarity transformation framework, the identifiability problem becomes a system equivalence problem: a system structure is identifiable if no equivalent systems exist for θ₁, θ₂ [set membership]

Θ and θ₁ ≠ θ₂ [106].

Knowledge about Lie algebra is needed to better understand the similarity transformation method; however, Lie algebra itself is a very rich topic, which will not be introduced in detail here. The interested reader is referred to [36]. Based on the work of Hermann and Krener [44], Vajda et al. [106] eventually proposed the similarity transformation approach to verify the global identifiability, for which a set of partial differential equations need to be formed and solved.

In summary, before the similarity transformation method can be applied, it is required that the system is both controllable and observable. Further more, a set of partial differential equations need to be generated and solved [106] to verify the system identifiability. These two disadvantages make the similarity transformation method not feasible for general nonlinear systems in practice.

3.2. Direct Test

Recall the definition of global (or local) identifiability, the key is to verify whether the same system output will result in a unique set of parameter values. That is,

should be satisfied either globally or locally if the model is identifiable. Based on this sufficient condition, Denis-Vidal and Joly-Blanchard [25] proposed to verify the identifiability of uncontrolled and autonomous systems by directly comparing the right hand side function f in Eq. (2.1). Note that here f = f(x(t), θ) does not explicitly depend on t and u(t) for uncontrolled and autonomous systems. Therefore, the problem becomes whether

can hold globally or locally. Denis-Vidal and Joly-Blanchard [25] used the following model for quantifying microbial growth [46] to illustrate the basic idea

(3.11)

(3.12)

(3.13)

Therefore, the right hand side function vector is

and from f(x, l, θ₁) = f (x, l, θ₂) we have

(3.14)

(3.15)

Solving the two equations above, we have

which indicates that only parameters (K_d, μ_m) are identifiable but the rest are not.

Although the analytical direct test approach described above is conceptually simple, it usually requires advanced mathematical skills to obtain analytic solutions and hence is difficult to be applied in practice. If a certain number of state variables are not measured, such latent variables have to be eliminated first (e.g., by taking higher order derivatives) in order to use the analytical direct test approach. It may be necessary to employ computer algebra tools, instead of algebraic manipulations by hand, for complicated models as suggested by Raksanyi et al. [87]. However, computer algebraic computation can easily become unfeasible for complicated nonlinear ODE models and Walter et al. [112] illustrated that the conclusions drawn from the analytical direct test approach can be misleading for certain types of models. Instead, Walter et al. [112] proposed the numerical direct test approach. For a model to be identifiable, Walter et al. [112] considered the following conditions which should be satisfied in practice,

where δ is a positive cut-off value chosen by the user. Techniques such as the algorithm SIVIA and forward-backward contractor for constraint satisfaction problems (CSP) were employed to find the inner and outer approximations of the solution set An external file that holds a picture, illustration, etc.
Object name is nihms293681ig1.jpg

An external file that holds a picture, illustration, etc.
Object name is nihms293681ig1.jpg

, that is,

For details of the algorithms for solving CSP (called interval arithmetic), see [53].

Walter et al. [112] thought that if An external file that holds a picture, illustration, etc.
Object name is nihms293681ig2.jpg

is empty, the model is identifiable; if An external file that holds a picture, illustration, etc.
Object name is nihms293681ig3.jpg

is not empty, then the model is not identifiable. However, the choice of the cut-off value δ is arbitrary, which seriously restricts the application of the numerical direct test method. In the parameter space of some models, there may exist a continuous and flat hypersurface on which the objective function (a function to be minimized for an optimization problem which evaluates how good a solution is, e.g., the residual sum of squares) has the same minimum value, which suggests the unidentifiability of the model. Under such circumstances, the numerical direct test approach, however, may still misleadingly conclude that the model is identifiable. In addition, it is difficult to verify which parameters are identifiable and which are not by using the numerical direct test method. Therefore, no useful information can be derived from An external file that holds a picture, illustration, etc.
Object name is nihms293681ig3.jpg

to help to improve mathematical models by reparameterizing unidentifiable parameters. Thus, the application of the direct test approach is very limited in practice.

3.3. Differential Algebra

The methods discussed in the previous subsections are difficult to apply to general nonlinear systems due to the difficulties in developing sufficient or necessary conditions for system identifiability and solving the resulting equations corresponding to such conditions. Also, rigorous training and advanced mathematical skills are required to use these methods. Is it possible to leave such tedious algebraic calculations to a computer instead of a human? The idea has motivated the development of methods under the framework of differential algebra [91] and have yielded some promising results [7, 14, 64].

Compared to other methods, the differential algebra approach has the following advantages: well-established theories, feasibility to general nonlinear dynamic systems, and availability of several computational algorithms (e.g., [14, 49, 57, 91]) and software packages (e.g., diffalg in Maple^© by Hubert [49], DAISY by Bellu et al. [12]). Theories and algorithms developed in abstract algebra and computer algebra are very helpful to understand differential algebra. For details of abstract algebra and computer algebra, the interested reader is referred to [28] and [72]. For details of differential algebra, the interested reader is referred to [14, 57, 64, 78, 91]. Here we only review some important concepts, theories and algorithms of differential algebra.

The first important concept is that of a differential polynomial. Here we give the definition for general dynamic systems:

Definition 3.2

Differential polynomial: If an expression is constructed from variables t, x, u and y, parameter θ = (θ_c, θ_t) and constants using only the operations of addition, subtraction, multiplication, constant positive whole number exponents, and constant positive whole number derivatives, it is called a differential polynomial.

For example,

(3.16)

is a valid differential polynomial. Note that for the problems considered in this paper, the derivatives in the definition above are with respect to time t only.

Now the dynamic system in Eqs. (2.1) and (2.2) can be rewritten as

(3.17)

(3.18)

If the left-hand side of Eqs. (3.17) and (3.18) is in the form of a differential polynomial after necessary algebraic computation or transformation, the structural identifiability of this system can be investigated in the differential algebra framework.

Clearly, an infinite number of differential polynomial equations can be formed by adding, scaling, multiplying, and differentiating both sides of Eqs. (3.17) and (3.18). It can be easily proved that the solution to Eqs. (3.17) and (3.18) is also the solution to all those generated equations. Therefore, the structure identifiability of Eqs. (3.17) and (3.18) can be investigated from those infinite number of generated equations. Let

{v₁, v₂, …, v_r} denote the differential polynomial ring with the differential indeterminates v₁, v₂, …, v_r [91]. For the dynamic systems under consideration, v_i [set membership]

V, i = 1, 2, …, r can be any component of x, y, u and θ, and

{v₁, v₂, …, v_r} is the set of the infinite number of generated differential polynomials. As mentioned above, the derivative on the ring

{v₁, v₂, …, v_r} is with respect to time t only, such a ring is called an ordinary differential ring.

Before we introduce more properties of

{v₁, v₂, …, v_r}, some definitions and concepts on differential indeterminates and polynomials need to be described. First, the order of a differential indeterminate is defined as the order of the derivative of that indeterminate and the degree of a differential indeterminate is defined as the exponent of that indeterminate. For example, in the first term equation M42

in Eq. (3.16), the order of y₂ is one and the degree of ₂ is two. To compare multiple differential polynomials, ranking needs to be defined [64, 91]:

Definition 3.3

Ranking: A total ordering of all the indeterminates and their derivatives is called a ranking if

where

means ‘ranks lower than’.

Note that for the same indeterminate v [set membership]

V, the item with a higher degree ranks higher, e.g., v [precedes]

v². The following are two examples of ranking:

(3.19)

(3.20)

For a given ranking over a differential polynomial ring

{v₁, v₂, …, v_r}, the highest ranking derivative in a differential polynomial P [set membership]

{v₁, v₂, …, v_r} is called the leader of P. Therefore, to rank two differential polynomials P₁ and P₂, the leaders of P₁ and P₂ are compared first, then the second highest ranked derivatives are compared if the leaders of P₁ and P₂ rank the same, and so on. To generalize this ranking concept to differential polynomials, the concepts of partially reduced and reduced were also introduced [64, 91]:

Definition 3.4

Partially reduced: For two differential polynomials P₁, P₂ [set membership]

, let v_P₁ denote the leader of P₁, P₂ is said to be partially reduced with respect to P₁ if there exists no proper derivative of v_P₁ in P₂.

Definition 3.5

Reduced: For two differential polynomials P₁, P₂ [set membership]

, let v_P₁ denote the leader of P₁, P₂ is said to be reduced with respect to P₁ if P₂ is partially reduced with respect to P₁ and the degree of v_P₁ in P₂ is less than the degree of v_P₁ in P₁.

With the definitions above, an autoreduced set can be introduced as follows:

Definition 3.6

Autoreduced set: A differential polynomial set is said to be an autoreduced set if any two differential polynomials in this set are reduced with respect to each other.

Autoreduced sets can also be ranked [64]. For two autoreduced sets A = {A₁, A₂, …, A_r} and B = {B₁, B₂, …, B_s}, A ranks lower than B if there exists an integer k, 1 ≤ k ≤ min(r, s) such that rank A_i = rank B_i (i = 1, 2, …, k − 1) and A_k [precedes]

B_k. Consider Eqs. (3.17) and (3.18) again, since an infinite number of differential polynomials can be generated with admissible operations, an infinite number of autoreduced sets can also be generated. Among these autoreduced sets, the set ranking the lowest is the most important and is called the characteristic set:

Definition 3.7

Characteristic set: Among all the autoreduced sets formed from a finite number of differential polynomials, the set ranking the lowest is called a characteristic set.

We now explain the relationship between the characteristic set and structure identifiability. The identifiability problem is to verify whether θ can be uniquely determined from differential polynomials with indeterminates u, y and θ only. Obviously, there exist an infinite number of sets of differential polynomials that can be employed to perform the identifiability analysis. However, the characteristic set has been proved to be the ‘best’ set among all such sets [91], where the word ‘best’ means the lowest rank. In summary, the characteristic set has the following properties:

Differential polynomials in the characteristic set satisfy Eqs. (3.17) and (3.18);
Differential polynomials in the characteristic set are in the simplest form possible;
Differential polynomials in the characteristic set have the exact information as in Eqs. (3.17) and (3.18) to verify system identifiability.

A number of algorithms have been developed to find the characteristic set, e.g., the Ritt algorithm [91], the Ritt-Kolchin algorithm [57], and the improved Ritt-Kolchin algorithm [14]. The implementation of such algorithms can be found in the diffgrob2 package [67] or the diffalg package [49]. The basic idea of these algorithms is to eliminate the higher ranking variables such as x so that differential polynomials with indeterminates u, y and θ can be obtained via symbolic computations. The key operation in the elimination process is called pseudo-division. Before we discuss the details of the pseudo-division, more definitions and notations need to be introduced. First, we call the coefficient of the highest power of the leader the initial. In addition, for a differential polynomial P and its leader v_P, we call the initial of Ṗ the separant of P, denoted by S_P. For example, using the ranking (3.19) combined with x₁ [precedes]

₁

···

x₂

₂

···, the initial of (3.16) is (−5₁ÿ₂) and the separant is (−10₁x₂ÿ₂).

Considering two differential polynomials P₁ and P₂, assume the leader of P₂ is v_P₂ and there exists a proper derivative of v_P₂, e.g., equation M46

for k ≥ 1 in P₁, then P₁ can be partially reduced by P₂ as follows [49]. First, take derivatives of P₂ up to the k^th order

where S_P₂ is the separant of P₂, and R_i(i = 1, …, k) the rest terms. Second, from the last equation above, equation M48

can be expressed in terms of equation M49

, S_P₂ and R_k, none of which contains equation M50

. Finally, substitute the expression of equation M51

into P₁ and then P₁ can be re-written as

where r is an integer, and Q and are differential polynomials. Furthermore, Q can be called a pseudo-quotient; and contains no proper derivatives of equation M53

and can be called a pseudo-remainder. The procedure described above is called pseudo-division. In this way, a set of differential polynomials can be reduced to generate an autoreduced set, and eventually, the characteristic set. More details of these computational algorithms can be found in [14, 49, 57]. However, we notice that the algorithms to find the characteristic set are still under development and the existing software packages do not always work well.

Ljung and Glad [64] concluded that each differential polynomial in the characteristic set can be in one of the following three forms if the ranking (3.19) is employed

where A, B and C are differential polynomials and the subindex m, n and l denote the number of differential polynomials. Ljung and Glad [64] proved a theorem about {B₁, B₂, …, B_n} to verify the structural identifiability:

Theorem 3.8

Assume that no separant or initial of {B₁, B₂, …, B_n} is zero,

If for some B_i, 1 ≤ i ≤ n, in the characteristic set one has B_i = _i, then the model structure is not identifiable;
If all B_i, 1 ≤ i ≤ n, in the characteristic set are of order zero and degree one in θ_i, and there exists non-degenerate solution (y, u, θ_*, x) for some θ_*, the model structure is globally identifiable at θ_*;
If all B_i, 1 ≤ i ≤ n, in the characteristic set are of order zero in θ_i, and some B_j is of degree greater than one in θ_j, and there exists non-degenerate solution (y, u, θ_*, x) for some θ_*, the model structure is locally identifiable at θ_*.

Although the results in Ljung and Glad [64] are for time-invariant parameters, they can be easily extended to time-varying parameter cases by treating θ_t as state variables or system inputs [7].

3.4. Implicit Function Theorem

Another approach based on the implicit function theorem was proposed by Xia and Moog [119]. For a general introduction of the implicit function theorem, the reader is referred to [73]. The theorem for identifiability analysis based on the implicit function theorem is given as follows [119]:

Theorem 3.9

Let Φ : R^d+p+q → R^q denote a function of model parameter θ [set membership]

R^q, system input u [set membership]

Rⁿ, system output y [set membership]

R^d and their derivatives, that is,

where k is a non-negative integer. Assume that Φ has continuous partial derivatives with respect to θ. The system structure is said to be locally identifiable at θ_* if there exists a point equation M56

such that

We can easily prove this theorem by considering the Taylor expansion of Φ at θ_*

since Φ(θ_*) = 0 and equation M59

exists (i.e., equation M60

), a unique solution of θ exists and the system is locally identifiable at θ_*.

Carefully examining this theorem, we find that it is the same as the algebraical identifiability definition 2.9. The implicit function theorem method can be employed alone to verify system identifiability, and it can also be used as a supplement to the differential algebra method. Theorem 3.8 suggests verifying system identifiability by examining specific forms of differential polynomials B = {B₁, B₂, …, B_n} in the characteristic set. However, a more rigorous approach is to verify whether equation M61

as suggested by the implicit function theorem method, where Φ is generated from B = {B₁, B₂, …, B_n}.

Xia and Moog [119] and Jeffrey and Xia [56] proposed a method to generate the function Φ by taking higher order derivatives of system output y to eliminate all latent (unobservable) state variables. For example, consider the following HIV model [81]

(3.21)

(3.22)

(3.23)

where T is the number of uninfected T cells susceptible to infection (target cells), T^* the number of infected T cells, and V the viral load. Using the method of Xia and Moog [119], we can take the 3rd-order derivative of V to obtain

(3.24)

Therefore, we can define

(3.25)

(3.26)

such that Φ = 0 is automatically satisfied. If equation M68

, then θ = (β, ρ, ν, μ, η) can be identified according to the implicit function theorem, where ν = δc, μ = δ + c, and η = Nλβδ. That is, if

(3.27)

we can identify the five parameters (β, ρ, ν, μ, η) in the model; furthermore, in the original model, N and λ cannot be identified separately. The identification function Φ involves the 7^th order derivative of V which requires at least 8 measurements of V to evaluate. Such information, a by-product of identifiability analysis, is useful to design new experiments.

Note that with a high dimensional parameter space, like in the example that we are considering here, the matrix equation M70

can become very complicated to derive and its rank is difficult to evaluate. For example, one element in the matrix (3.27) is

Thus, it is not easy to evaluate the rank of the above matrix. To avoid such complexity and evaluation of the high order derivatives, an alternative method for formulating the identification functions Φ(·) was proposed by Wu et al. [117]. Suppose we have the quantities (V, , , V⁽³⁾) at five distinct time points t₁, ···, t₅. Denote the values of (V, , , V⁽³⁾) at t = t_i as ( equation M72

) for i = 1, ···, 5. Let equation M73

, where f is specified in Eq. (3.25). Then we have

(3.28)

(3.29)

by the implicit function theorem, there exists a unique solution of θ. Assuming that β ≠ 0, some algebra shows that the rank of ( equation M76

) is equal to the rank of

(3.30)

As long as det(Σ) ≠ 0, we have equation M78

. Note that in Eq. (3.30), the matrix Σ also involves unknown parameters (μ, ρ, β); thus, to numerically determine Σ’s rank, nominal values of these parameters (i.e., obtained from literature) are needed.

Note that although V ⁽³⁾ is not involved in the matrix Σ, V ⁽³⁾ should exist at any time point. For evaluating V ⁽³⁾ at one time point, at least four measurements of V are needed. In order to form the five identification equations, at least eight measurements of V are necessary. This conclusion is consistent with that of the method proposed by Xia and Moog [119]. Note that this model is more likely to be locally identifiable than globally identifiable since Σ also contains unknown parameters. Compared to the method of Xia and Moog [119], the method of Wu et al. [117] is less computationally intensive and easier to implement since only the lower-order derivatives (the 3^rd or lower order derivatives in our case) of V need to be evaluated.

In this section, as a by-product of structural identifiability analysis, we illustrate how to calculate the minimum number of observation points for parameter identifiability, which is an important issue in experimental design. Thorough and in depth discussions of experimental design for dynamic systems are beyond the scope of this document, so we only briefly discuss the influence of the observation timing on identifiability here. Sontag [100] reported a very general and simple conclusion: for any ODE model with q unknown constant parameters, 2q + 1 experimental observations are enough for identification of the q parameters if measurements are absolutely accurate. Sontag [100] explained that the number 2q+1 is frequently met in geometry and dynamical system theory and it is the embedding dimension in the Whitney’s theorem [115, 116] for abstract manifolds or the embedding dimension of q-dimensional attractors [102]. However, an intuitive explanation for the number 2q is that each parameter in the righthand side of an ODE model is used to quantify the change of the state variables (like a slope) and at least two data points are needed to determine a slope. As to the influence of the observation times on identifiability, the work of Sontag [100] indicated that increasing the number of data points will not help to identify more parameters of a dynamic model once the model enters its steady state and produces only flat responses all the time. In principle, data points collected to capture violent nonlinear behavior of the dynamic system will be more informative for determining parameter values, as suggested in [103].

4. Practical Identifiability Analysis

Before we introduce various techniques of practical identifiability analysis, note that structural identifiability analysis provides a theoretical ground for practical identifiability analysis. If the structural analysis suggests that a model is not theoretically identifiable, the practical analysis is not necessary since theoretical unidentifiability must imply practical unidentifiability. Thus, only theoretically identifiable models need further practical identifiability analysis.

Structural identifiability analysis can be done without any actual experimental observation, so it is also called prior identifiability analysis. There are two basic assumptions upon which structural identifiability analysis heavily rely: model structures are absolutely accurate and measurements are exact (no measurement errors). However, these two assumptions are clearly not valid in practice. For instance, in biomedical research, both model uncertainty and measurement error are usually large. Therefore, even when structural identifiability analysis suggests that model parameters can be uniquely identified, the estimates of model parameters may still be unreliable. Thus, it is necessary to evaluate whether structurally identifiable parameters can be reliably estimated with acceptable accuracy from noisy data. This is so-called practical identifiability analysis or posterior identifiability analysis. It is strongly recommended to perform both structural and practical identifiability analyses in practice to insure the reliability of parameter estimation. For the rest of this section, we assume our measurement or output model to have measurement errors as follows:

(4.1)

where ε(t) is measurement error with mean 0 and variance σ²(t).

4.1. Monte Carlo Simulation

The history of Monte Carlo simulations can be traced back to the work of Metropolis and Ulam [69]. As implied by its name, this method is a sampling technique using random numbers and probability distributions. More specifically, the Monte Carlo simulation method define possible inputs first (e.g., measurement noise level), then randomly generate inputs according to certain probability distributions (e.g., normal distribution with zero mean), and then use the inputs to do certain calculations (e.g., add random errors to data and fit model to the simulated noisy data), and finally aggregate individual computation results (e.g., the average error in parameter estimates). It is not only useful for practical identifiability analysis but also helpful for experimental design. Monte Carlo simulation is very popular and widely used to assess the performance of statistical estimation methods in the statistical literature.

Once parameters or a subset of parameters of a model are determined to be theoretically (structurally) identifiable, one can use Monte Carlo simulations to evaluate whether the theoretically identifiable parameters can be reliably estimated with acceptable accuracy from noisy data. Obviously, in order to evaluate the practical (statistical) identifiability, statistical estimation methods, such as the least squares approach, need to be readily available. However, statistical parameter estimation for nonlinear ODE models is beyond the scope of this review paper and will not be reviewed here.

Monte Carlo simulations allow us to simulate various scenarios with different numbers of observations at different levels of noise or measurement error for different experimental designs although such designs may not be feasible for practical experiments. The simulated data can be used to evaluate whether model parameters or a subset of parameters can be reliably estimated under different conditions. In general, a Monte Carlo simulation procedure can be outlined as follows:

Determine the nominal parameter values θ₀ for simulation studies, which can be obtained by fitting the model to experimental data if available. Otherwise it can be obtained from the literature or other resources.
Use the nominal parameter values to numerically solve the ODE model to get the solution of the output or measurement variables at the experimental design time points.
Generate N sets (e.g., 1000) of simulated data from the output or measurement model (4.1) with a given measurement error level.
Fit the ODE model to each of the N simulated data sets to obtain parameter estimate _i, i = 1, 2, …, N.
Calculate the average relative estimation error (ARE) for each element of θ as

(4.2)
where is the k-th element of θ₀ and is the k-th element of _i.

The ARE can be used to assess whether each of the parameter estimates is acceptable or not. For a very small measurement error, the parameter estimates should be close to the true values and the ARE should be close to 0. When the measurement error increases, the ARE of the parameter estimates will also increase. However, the ARE for some of the parameter estimates may increase significantly and some others may just increase a little. However, for a reasonable or practical level of measurement error, if the ARE of a parameter estimate is unacceptably high, we claim that this parameter is not practically or statistically identifiable. In practical applications, some parameters may not be sensitive to measurement errors and can always be well estimated, some other parameters may be quite sensitive to measurement errors and their AREs are large even for a small measurement error, and at the same time, some parameters may be in the middle ground [70]. In addition, there is no clear cut rule on how high the ARE need to be before they are claimed “unacceptable” for a particular problem. Thus, the practical identifiability relies on the underlying problem and judgement of the investigators. Also notice that various statistical estimation approaches can be employed to obtain the parameter estimates, and the ARE may depend on the estimation methods.

Monte Carlo simulations can also be used to design better practical experiments. Different designs for different sample sizes under different conditions can be evaluated using the AREs. The best design and the necessary sample size can be determined based on the Monte Carlo simulation results. We will illustrate the application of this method in § 6.1.

4.2. Correlation Matrix

Although the Monte Carlo simulation approach is easy to understand and simple to implement, the associated computational cost is high since a large number of model fits to data need to be performed. Rodriguez-Fernandez et al. [92, 93] proposed an alternative approach for practical identifiability analysis of ODE models by examining the correlations between model parameters. This requires much less computation and is relatively simple if measurement errors follow an identical and independent distribution (i.i.d.).

The idea behind this approach is straightforward. Assume that the parameter estimate = [₁, ₂, ···, _q]^T has been obtained after fitting a model to experimental data. The correlation matrix of the parameter estimates can then be calculated based on the Fisher information matrix (FIM) [33, 109] in the following form:

(4.3)

where r_ij (i, j = 1, 2, …, q and −1 ≤ r_ij ≤ 1) is the correlation coefficient between parameter estimates _i and _j. If there exists a strong positive correlation between parameter estimates _i and _j, that is, the correlation coefficient r_ij is close to 1, parameters θ_i and θ_j are said to be practically undistinguishable. A strong correlation between two parameters indicates that one parameter strongly depends on another parameter and these two parameters cannot be separately estimated.

A derivation of the expression for the correlation matrix was provided by Rodriguez-Fernandez et al. [92]. For simplicity, the measurement errors were assumed to be uncorrelated and follow an identical normal distribution with mean zero, that is, N(0, σ²). In this case, for a general dynamic system (2.1) and (2.2), the Fisher information matrix can be given as

(4.4)

where the subscript i denotes the i^th time point of experimental observation, N the total number of observations, ŷ_i the model approximation of observation, the model parameter estimate, and V a known positive definite matrix of weights on variances. It can be proved that the covariance matrix C is equal to the inverse of FIM according to the Cramèr-Rao theorem [89], that is

(4.5)

Finally, the element r_ij of the correlation matrix can be defined as

(4.6)

(4.7)

Guedj et al. [40] tackled the practical identifiability problem of HIV dynamic models. They developed their approach under the framework of maximum likelihood estimation instead of least squares estimation, but their results are still based on the Fisher information matrix and the idea is similar to that in [92].

A limitation of the correlation matrix approach is that it requires not only the parameters but also their correlation matrix to be reliably estimated. This may be a problem for an model with most parameters unidentifiable since the correlation matrix estimate may strongly depend on the parameter estimates. If any two parameters are not distinguishable, the parameter estimates and their correlation matrix estimate may be poor. In addition, the correlation matrix approach only allows you to check whether any pair of parameters are distinguishable or not; to evaluate correlations between more than two parameters, the sensitivity-based identifiability analysis techniques (e.g., eigen-decomposition of the sensitivity matrix) should be considered, as described in the next section.

5. Sensitivity-Based Identifiability Analysis

Sensitivity analysis itself is a rich topic. The interested reader is referred to the comprehensive survey by Saltelli et al. [94] and Cacuci [13]. Sensitivity analysis (SA) is often used to assess the variation of system output induced by different input factors including model parameters. The sensitivity analysis idea can also be used to evaluate the identifiability of unknown parameters.

Sensitivity-based identifiability analysis is similar to the structural analysis approach in the sense that both approaches do not require actual experimental data (although the sensitivity-based method could require the number and locations of measurement time points, see details below), and both approaches assume that measurements are precise without error. However, the sensitivity-based method does not directly use the model structure information, which is a critical difference between the structural and the sensitivity-based approaches. The sensitivity-based method is similar to the practical analysis approach in the sense that both methods require pre-specified parameter values (either nominal or actual estimates), and both need to know the number and locations of measurement time points. However, the sensitivity-based method is different from the practical analysis approach in the sense that the sensitivity-based method does not take measurement error into account. Thus, the sensitivity-based method is a technique between the structural (theoretical) identifiability and practical identifiability analysis. We review such methods in this section.

A nominal parameter value is required for the sensitivity-based approach. Thus, parameter identifiability is evaluated with respect to a specific point in the parameter space by sensitivity-based methods. For this reason, the concept of at-a-point identifiability was introduced by Ljung and Glad [64] and Quaiser and Mönnigmann [86] as follows,

Definition 5.1

Globally at-a-point identifiable: Let θ^* denote a fixed point in the parameter space Θ. A system is said to be globally at-a-point identifiable if for any admissible input u(t) and any parameter vector θ [set membership]

Θ, y(u, θ) = y(u, θ^*) implies θ = θ^*.

Definition 5.2

Locally at-a-point identifiable: Let θ^* denote a fixed point in the parameter space Θ. A system is said to be locally at-a-point identifiable if for any admissible input u(t) and any parameter vector θ within an open neighborhood of θ^*, y(u, θ) = y(u, θ^*) implies θ = θ^*.

The sensitivity-based identifiability analysis techniques reviewed in this section only examine at-a-point identifiability. The sensitivity of measurable system responses with respect to parameter values is used to assess the identifiability of unknown parameters by these methods. More specifically, assume that the locations and the number of time points at which the system responses or state variables will be measured have been given, denoted by t₁ ≤ t₂ ≤ ··· ≤ t_N, then the sensitivity coefficient at each time point t_k (k = 1, 2, …, N) for a given nominal parameter vector θ^* is defined as

(5.1)

where y_i(i = 1, 2, …, d) denotes the i^th component of y (y [set membership]

R^d) and θ_j the j^th(j = 1, 2, …, q) component of θ (θ [set membership]

R^q). Thus, the sensitivity matrix for all time points is defined as

(5.2)

A number of identifiability analysis techniques have been developed based on this sensitivity matrix. Simply speaking, the larger the sensitivity coefficients are, the more notable the system responses or measurable state variables are with respect to the changes of parameters. In that sense, a parameter is likely to be identifiable if the system output is highly sensitive to a small perturbation of this parameter; otherwise, the parameter is likely to be unidentifiable. In addition, if there exists a strong correlation between any two parameters, these two parameters are very likely to be indistinguishable from each other. Such parameter dependence can also be evaluated by examining the dependence of the sensitivity matrix columns. We review four typical methods along this line in detail: the correlation method [51, 93, 122], the principle component analysis (PCA) method [24, 34, 54], the orthogonal method [35, 120, 121], and the eigenvalue method [85, 86, 96, 108].

5.1. Correlation Method

The correlation method was first proposed by Jacquez and Greif [51]. The method was originally developed to study identifiability for linear compartment models and the derivation was given for a single output system only (that is, y [set membership]

R). However, this method is not limited to linear models and single output systems.

Consider the first order Taylor expansion of the system output near the pre-specified nominal parameter vector θ^*

(5.3)

where k = 1, 2, …, N denotes the index of the measurement time points. Let r_k denote the measurement at t_k without errors and Δθ = θ − θ^*, then the residual sum of squares between the exact measurements and the linear approximation is

(5.4)

where r_k − y_k(θ^*) = 0 based on our assumptions. Finally, we can rewrite Eq. (5.4) in terms of the sensitivity matrix,

(5.5)

where S is the sensitivity matrix defined in (5.2). Obviously, the minimum of RSS(Δθ) is reached at S^T S · Δθ = 0. If S^T S is of full rank, the unique solution of S^T S · Δθ = 0 is = θ^*, which indicates that the model parameters θ are locally identifiable at θ^*. If S^T S is singular, then there exists at least one non-trivial solution ≠ θ^* such that the model parameters are not identifiable at θ^*. Two important issues should be noticed: first, a similar expression can be derived under the framework of maximum likelihood estimation [40], so the derivation is not limited to ordinary least squares; secondly, only local identifiability can be inferred based on the rank of S^T S since the linear approximation Eq. (5.3) is used [86].

It is also desirable to determine which parameters are not identifiable if S^T S is not of full rank. For this purpose, the correlations between parameters can be calculated based on the sensitivity matrix. More specifically, if we examine the columns of the sensitivity matrix defined in Eq. (5.2), it is clear that each column is the sensitivity of the system responses at all time points with respect to one particular parameter. Thus, the sample correlation of two columns is an estimate of the correlation between two parameters that can be calculated as

(5.6)

where S _· _i (or S _· _j) denotes the i^th (or j^th) column of the sensitivity matrix S, cov(S _· _i, S _· _j) the sample covariance between S _· _i and S _· _j, and σ (S _· _i) and σ (S _· _j) the sample standard deviations of S _· _i and S _· _j, respectively. If the calculated correlation coefficient between any two parameters is close to one, these two parameters are not distinguishable. However, such a decision always involves a pair of parameters. Is it possible to determine which one in this pair is more problematic and should be fixed or removed from the model? Quaiser and Mönnigmann [86] proposed the concept of total correlation for this question,

(5.7)

where I denotes the indicator function, and δ [set membership]

(0, 1) the cut-off value specified by the user. The parameter with the highest total correlation is the first candidate to be fixed or removed from the model.

If we compare the correlation method with the correlation matrix method for practical identifiability analysis described in the previous section [92], we can find that the Fisher information matrix and the sensitivity matrix are somehow similar; however, the ways of calculating correlations are different.

5.2. Tuning Importance Method and Principle Component Analysis

This category of methods have been developed to reduce model complexity by discarding nonsignificant parameters, such as the tuning importance method [24, 95, 105] and principle component analysis(PCA) [24, 34, 54]. More specifically, both methods rank all parameters first, and then these parameters are determined as identifiable or unidentifiable according to their ranks.

One important and interesting feature of the tuning importance and PCA methods is that they are based on the normalized sensitivity matrix, which is different from the sensitivity matrix defined in Eq. (5.2). To construct the normalized sensitivity matrix, the dimensionless sensitivity coefficient was defined as follows [24, 30],

(5.8)

where i

{1, 2, …, d} denotes the index of system outputs, j [set membership]

{1, 2, …, q} the index of parameters and k [set membership]

{1, 2, …, N} the index of measurement time points. Then the normalized sensitivity matrix for each y_i is defined as

(5.9)

For the tuning importance method, the following objective function was introduced by Seigneur et al. [95] and Turányi [105],

(5.10)

where θ ₋_j denotes the parameter vector with the j^th component removed. By following the same procedure as the correlation method in the last subsection, the overall sensitivity can be obtained and expressed in terms of the dimensionless sensitivity coefficients

(5.11)

Thus, the larger the overall sensitivity of one parameter, the more sensitive the system response is with respect to small perturbations of this parameter. Since all the parameters can be ranked according to their overall sensitivities, the parameters ranking the lowest are candidates to be unidentifiable and to be discarded.

For the PCA method, the eigenvalues and eigenvectors of the matrix S^{i^T} Sⁱ are calculated to provide the information for model reduction. Let equation M99

denote the eigenvalues which are ordered non-decreasingly

(5.12)

Also, let the corresponding eigenvectors be denoted by

Three strategies examining the eigenvalues and eigenvectors were proposed to rank all the parameters by Jolliffe [54]:

Starting with the eigenvector corresponding to the smallest absolute eigenvalue, loop over each eigenvector to locate the component with the maximum absolute value and the corresponding parameter at the maximum component location is marked as unidentifiable and to be removed if it has not be marked before. This procedure is summarized as follows:
(5.13)
Starting with the eigenvector corresponding to the smallest absolute eigenvalue again, loop over each row of matrix Γⁱ and calculate the sum of squares of all components in each row. The parameter corresponding to the location of the row with the largest sum of squares is determined as unidentifiable for removal. This procedure is summarized as follows:
(5.14)
Starting with the eigenvector corresponding to the largest absolute eigenvalue, loop over each eigenvector component to locate the largest one and mark the parameter if it is not marked before. The marked parameter is not selected as an unidentifiable parameter immediately; instead, a rank is assigned to this parameter. Eventually, all parameters are ranked and the parameters ranking the lowest are thought unidentifiable. This procedure is summarized as follows:
(5.15)

Froemel [34] proposed a simple strategy to integrate the rankings from all the three strategies described above. Further details are given in [34] and [86].

5.3. Orthogonal Method

The orthogonal method was proposed by Yao et al. [120]. The basic idea of this approach is to examine the (nearly) linear dependencies of columns of the sensitivity matrix S defined in Eq. (5.2). Thus, both the sensitivity of system response with respect to parameter values and the dependency between parameters regarding their effects on the system responses can be simultaneously evaluated to determine a set of identifiable parameters.

Unlike the correlation method, the orthogonal method does not calculate the correlation between different columns of S. Instead, the perpendicular distance of one column to the vector space spanned by other columns is calculated as a measurement of the linear dependency. It is an iterative procedure. More specifically, at the first iteration, the column of S with the largest sum of squares is removed from S and selected as the first element of an empty set S_I. At the (j + 1)^th (j [set membership]

{1, ··· q − 1}) iteration, j columns have been removed from S and selected into S_I, and a vector space spanned by all the columns in S_I is denoted by V_{S_I}. For column S_h still in S, the orthogonal projection equation M105

of this column on the vector space V_{S_I} is calculated and the perpendicular is then obtained as

(5.16)

In the work of Yao et al. [120], the norm equation M107

was proposed to be the measurement of nearly linear dependency between the vector S_h and the vector space V since the shorter the distance is, the larger the dependence. For a norm equation M108

which is nearly zero, the corresponding column S_h is thought to be nearly linearly dependent and thus is not identifiable. At the first iteration, the column of S with the largest norm is selected into S_I. However, another reasonable alternative criterion is to consider the angle between S_h and equation M109

(that is,

) since this criterion can select the best candidate even if the norms of different columns are orders of magnitude different. Finally, at each iteration, a pre-specified cut-off value δ will be used for the perpendicular distance or the angle of all columns in S. Once the distance or angle of one column is smaller than δ, this column is thought to be linearly dependent on V_{S_I}; therefore, the column with the largest equation M111

is removed from S and selected into the set S_I. This procedure is repeated until S becomes empty. The vectors in S_I determine which parameters are identifiable, which is the primary goal of the orthogonal method (i.e., find the identifiable rather than unidentifiable parameters).

Since the cut-off value δ in this method is an arbitrary value specified by users, the number of unidentifiable parameters strongly depends on the selection of δ. Due to this problem, Quaiser and Mönnigmann [86] proposed to rank all the parameters based on the values of norms or angles instead of simply dividing them into identifiable or unidentifiable groups.

5.4. Eigenvalue Method

The eigenvalue method was first proposed by Vajda et al. [108] and then further developed by Quaiser et al. [85], Schittkowski [96], and Quaiser and Mönnigmann [86]. This approach is based on the properties of eigenvalues and eigenvectors of S^T S, where S is the sensitivity matrix defined in Eq. (5.2). To illustrate this method, consider the residual sum of squares between system outputs and experimental measurements

(5.17)

and let λ₁ ≤ λ₂ ≤ ··· ≤ λ_q denote the eigenvalues of S^T S in a non-decreasing order and (γ₁, γ₂, ···, γ_q) the corresponding normalized eigenvectors. Note that since S^T S is symmetric and positive semi-definite, all its eigenvalues are real and non-negative. Given a nominal parameter vector θ^*, the Taylor expansion of RSS at θ^* along an eigenvector is approximately

(5.18)

where α is an arbitrary small constant and [nabla]

RSS(θ^*) is the gradient of RSS at θ^*. Although [nabla]

RSS(θ^*) is not necessarily zero if r_k is not an exact measurement, since θ^* is a nominal parameter vector that may not minimize RSS, the second order term equation M114

can become zero if the eigenvalue λ_j corresponding to θ_j is equal to zero since S^T S · γ_j = λ_jγ_j and equation M115

. That is, along the direction of θ_j with λ_j = 0, the change of RSS is expected to be nearly zero. The selection criterion for unidentifiable parameters is given by

(5.19)

In practice, λ_j is usually not exactly zero, therefore, a cut-off value δ needs to be specified. For a detailed implementation algorithm, the interested reader is referred to [86].

The four sensitivity-based identifiability analysis methods described above were also reviewed and compared in Quaiser and Mönnigmann [86]. In general, all four approaches are applicable to general ODE models; however, the eigenvalue method and the orthogonal method are better designed to globally evaluate and compare the influences of parameter values on system outputs such that these two methods outperform the correlation method and the principle component analysis method. In addition, both the eigenvalue method and the orthogonal method are easy to implement. Note that, the method proposed in [88] made the assumption that the nominal parameter vector θ^* minimizes the objective function, which should be interpreted as r_k = y_k(θ^*) as in (5.4) for the correlation method. In addition, it should be mentioned that in practice, if the sensitivity matrix is of full rank but with eigenvalues of different orders of magnitude, the parameters corresponding to the smallest eigenvalues are theoretically identifiable but likely to be practically unidentifiable. Under such a circumstance, the sensitivity-based approaches are still useful in the sense of determining practically unidentifiable parameters.

Finally, it is also interesting to combine the sensitivity-based approaches with the practical identifiability methods introduced in Section 4 to study the identifiability of a dynamic system. Notice that the sensitivity-based approaches do not require statistical estimation of unknown parameters, which can be done before the practical identifiability analysis.

6. Application Examples

In this section, we illustrate the applications of both structural and practical identifiability analysis techniques through examples in modeling viral dynamics. We summarize the identifiability analysis results for popular models of HIV infection, influenza infection, and hepatitis B and C virus dynamics.

6.1. HIV Model with Constant Parameters

Miao et al. [70] proposed the following model to describe a growth competition assay to quantify HIV replication fitness:

(6.1)

(6.2)

(6.3)

(6.4)

where T, T_m, T_w and T_mw are the numbers of uninfected cells, cells infected by mutant viruses, cells infected by wild-type viruses, and cells infected by both mutant and wild-type viruses (dual-infection). Let (λ, λ_m, λ_w, λ_mw) represent the proliferation rates of T, T_m, T_w and T_mw and (δ, δ_m, δ_w, δ_mw) the death rates of T, T_m, T_w and T_mw, respectively. Then ρ = λ − δ, ρ_m = λ_m − δ_m, ρ_w = λ_w − δ_w and ρ_mw = λ_mw − δ_mw are the net growth rates of T, T_m, T_w and T_mw, which are the differences between the corresponding proliferation rates and death rates. Parameters (k_m, k_w, k_R) are infection rates of mutant virus, wild-type virus and recombinant virus respectively, and q_m and q_w are dual infection rates.

For this example, the implicit function method of Xia and Moog [119] is employed to investigate the structural identifiability. Since all state variables (T, T_m, T_w, T_mw) are experimentally measurable, the outputs of the system are

(6.5)

By taking derivatives of one of the four equations in this HIV viral fitness model, the structural identifiability of the model can be evaluated. To demonstrate this, here we start with the first equation,

(6.6)

By taking higher orders (up to the 4th order) of derivatives on both sides of Eq. (6.6), we get

(6.7)

(6.8)

(6.9)

When Eqs. (6.7)–(6.9) are written in matrix form, we find that the parameters are identifiable if

(6.10)

Note that the rank of this matrix can be evaluated numerically if the analytical form is not available, and nominal parameter values are not needed for this case since the matrix above involves no parameters. Since the left hand side of Eq. (6.9) has a derivative of order 4, at least five measurements of y₁ = T are needed to evaluate equation M127

, and at least four measurements of y₂ = T_m, y₃ = T_w, and y₄ = T_mw are needed to evaluate their derivatives of order 3. Since k_m and k_R can be identified from Eqs. (6.7)–(6.9) if Eq. (6.10) holds, Eq. (6.2) can be rewritten as

(6.11)

Similarly, by taking the higher derivatives of Eq. (6.11), the parameters (ρ_m, q_m) are identifiable if

(6.12)

By the same token, the parameters (ρ_w, q_w) are identifiable if

(6.13)

Finally, ρ_mw is identifiable if

(6.14)

In summary, all parameters (ρ, ρ_m, ρ_w, ρ_mw, k_m, k_w, k_R, q_m, q_w) are structurally identifiable if at least five measurements of y₁ = T and four measurements of y₂ = T_m, y₃ = T_w, and y₄ = T_mw are available and if all coefficient matrices (6.10), (6.12), (6.13) and (6.14) are of full rank at least for some local time points.

Monte Carlo simulations were also performed to evaluate the practical identifiability of this HIV viral fitness model by Miao et al. [70]. The AREs of all 9 parameters for three measurement error levels (0%, 5%, and 30%) are duplicated in Table 6.1. Note that in this simulation study, we assumed that there were 1000 replicates of data at each time point although this may not be feasible in practical experiments. But this will help us to evaluate whether the unknown parameters are practically identifiable when the sample size of the data (with noise) is large enough. We can see that when there is no measurement error (σ = 0%), all the 9 parameters can be well identified (the maximum ARE is 0.4%), which confirms our theoretical identifiability analysis results. This also indicates that the parameter estimation method is good and the parameter estimates converges to the true parameter values when the sample size is large enough and the measurement error is small enough. However, when the measurement error increases to 5% and 30%, the ARE of parameter ρ_mw rapidly increases to 556% and 2062%, respectively. The ARE of ρ_w also increases to 39% and 201%, while the ARE of k_R increases to 28% and 106%, respectively. The AREs of parameters (q_w, ρ_m) are reasonable for the case of small measurement error (σ = 5%), but increase to 49% and 59% for the large measurement error case (σ= 30%), respectively. The AREs for other four parameters (ρ, k_m, k_w, q_m) are reasonable for all cases.

Table 6.1

Practical identifiability analysis by Monte Carlo simulations. The average relative error (ARE) is calculated based on 1000 simulation runs; 1000 replicates at each time point are generated for each simulated data set (Table 6 in Miao et al. [70]).

To further investigate the practical identifiability of unknown parameters under practical experimental conditions, we performed more simulations for different numbers of time points and different numbers of replicates at each time point. The simulation results are given in Table 6.2. One can see that the ARE of parameter ρ_mw ranges from 410% to 2130%. This, combined with the results in Table 6.1, indicates that the ρ_mw is practically unidentifiable. Considering the practical case of 9 time points and 9 replicates for each time point, the ARE of parameter ρ_w is 86%, which indicate that it may be difficult to accurately identify the parameter ρ_w unless the sample size is unrealistically large (say, 100 replicates for each time point). For parameter k_R, the AREs are also large (ranging from 62% to 108%) for practical cases (the number of replicates is 3, 6, or 9). For parameters (ρ_m, q_w), the AREs are reasonable (ranging from 22% to 38%) for most reasonable sample sizes, thus (ρ_m, q_w) can be considered as reasonably identifiable. The parameters (ρ, k_m, k_w, q_m) are very well identified (the AREs ranging from 3% to 22%) in all cases.

Table 6.2

Practical identifiability analysis using Monte Carlo simulations based on 1000 simulation runs with the measurement error level σ = 1.5% (Table 7 in Miao et al. [70]).

6.2. HIV Model with Constant and Time Varying Parameters

In this section, we consider another dynamic system that is widely used to describe HIV dynamics in HIV-infected patients with antiretroviral treatment [17, 48, 82]:

(6.15)

(6.16)

(6.17)

where T_U is the concentration of uninfected target cells, T_I the concentration of infected cells, V (t) the viral load, λ the source rate of uninfected T cells, ρ the death rate of uninfected T cells, η(t) the time-varying infection rate which is a function antiviral treatment efficacy, δ the death rate of infected cells, c the clearance rate of free virions, and N the average number of virions produced by a single infected cell over its lifetime. T_U (t), T_I(t) and V (t) are state variables, and (λ, ρ, N, δ, c, η(t))^T are unknown dynamic parameters.

The differential algebra approach for structural identifiability analysis (Section 3.3) requires one to eliminate the latent (unobservable) state variables from the dynamic equation in order to evaluate the identifiability. The concept of ranking is introduced such that computer algorithms can be designed to eliminate variables or their derivatives with higher rank. For notation simplicity, let x₁, x₂ and x₃ denote T_U, T_I and V, respectively. In the dynamic model (6.15)–(6.17), we can measure viral load (x₃ = V) and total CD4+ T cell counts (x₁ + x₂ = T_U + T_I). Let y₁ and y₂ denote the measurable variables x₁ + x₂ and x₃, respectively. We adopt the ranking

(6.18)

where θ = [λ, ρ, N, δ, c]^T is the vector of constant unknown parameters. We can eliminate x₁, x₂ and x₃ using the ranking (6.18) to obtain

(6.19)

(6.20)

Note η (t) can be expressed in terms of measurable state variables and other unknown constant parameters either from Eq. (6.19) as

(6.21)

or from Eq. (6.20) as

(6.22)

Thus, the time-varying parameter η (t) is identifiable if all the constant parameters are identifiable.

To verify the identifiability of constant parameters θ, Eq. (6.21) and Eq. (6.22) can be combined to obtain

(6.23)

The above equation just involves measurable state variables and constant parameters. Eq. (6.23) is of order 0 and of degree > 1 in θ, so Eq. (6.23) satisfies the third situation in Theorem 3.8 [64] in Section 3.3), thus θ = (λ, ρ, N, δ, c)^T is locally identifiable. Consequently, η (t) is also locally identifiable. In addition, the identifiability of θ can also be easily verified using the implicit function method based on Eq. (6.23). The identifiability of other similar HIV dynamic models have been studied in [56, 117, 119].

6.3. Influenza A Virus Infection

The purpose of this section is to illustrate possible problems if identifiability analysis is ignored by considering influenza infection in humans, an important infectious disease. Baccam et al. [8] proposed a target cell-limited model for influenza A virus infection:

(6.24)

(6.25)

(6.26)

where T is the number of uninfected target cells (epithelial cells), I is the number of productively infected cells, and V is the infectious viral titer expressed in TCID₅₀/ml which is the only state variable to be measured. Since this is a low-dimension nonlinear dynamic system, the implicit function method can be employed.

Considering the case that only V can be measured (and thus, for example, the initial number of target cells T(0) is not known), we can derive the follow equation from Eq. (6.24)–(6.26) by eliminating the unmeasurable state variables:

(6.27)

Obviously, only the parameters (β, δ, c) can be identified and the minimum number of required measurements of V is 6, and the parameter p is not identifiable in this case. Similarly, we consider the cases in which both I and V are measured, both T and V are measured, or all three state variables are measured. For the cases that any two or more state variables are measured, all four parameters (β, δ, c, p) are structurally (theoretically) identifiable. We summarize the structural identifiability analysis results for all cases in Table 6.3.

Table 6.3

Structural identifiability of the target cell-limited influenza model in Baccam et al. [8].

Baccam et al. [8] further proposed another target cell-limited influenza model with delayed virus production as follows:

(6.28)

(6.29)

(6.30)

(6.31)

where I₁ is the number of latent infected epithelial cells that are not yet producing virus and I₂ the number of productively infected epithelial cells. Again, we can use the implicit function theorem method to investigate the structural (theoretical) identifiability of this model. We summarize the identifiability analysis results in Table 6.4. For this model, if only V is measured, four parameters (β, δ, c, k) are identifiable and parameter p is not identifiable. However, if any two or more of the four state variables (T, I₁, I₂, V) are measured, all five parameters (β, δ, c, k, p) are theoretically identifiable. We also summarize the minimum number of required measurements for each of the state variables in Table 6.4.

Table 6.4

Structural identifiability of the target cell-limited influenza model with delayed virus production in Baccam et al. [8].

In the paper by Baccam et al. [8], only the virus titers were measured at 8 time points during days 1–8 of infection from six patients, but some of these measurements were below detection. According to the identifiability analysis in Table 6.3 and Table 6.4p is not identifiable. To fit the four dimensional model (6.28)–(6.31), all 8 data points need to be used, otherwise more parameters may be unidentifiable. However, since the identification equation (6.27) does not involve the unidentifiable parameter p, one may fix p, which does not affect the estimates of other parameters. Since the identifiability analysis was not considered, the estimates of kinetic parameters in Baccam et al. [8] should be interpreted with caution since Baccam et al. [8] fixed T(0) to avoid the identifiability problem and T(0)’s value was not directly from data.

Influenza infection in human is a very complex problem, therefore, much more complicated models were proposed; however, the problem is that such work usually over-parameterizes the model and ignores parameter identifiability, which makes it difficult to directly such models to data. For example, Hancioglu et al. [43] proposed a 10-equation model for influenza A virus infection (details not shown). We can show that, to verify the identifiability of all the 27 parameters in that model, almost all the 10 state variables need to be measured, which is nearly impossible to do in practice due to technical and ethical limitation. In summary, when fitting model to data, models with parameter identifiability verified should be considered.

7. Discussion and Conclusion

Ordinary differential equations are an important tool for quantifying a dynamic process in many scientific fields, and recently it has been widely used in modeling biomedical processes, in particular for modeling infectious diseases and viral dynamics. It is critical to estimate the unknown kinetic parameters in ODE models from experimental data in biomedical applications. However, it is not trivial and not apparent to know whether the unknown parameters in general nonlinear ODE models are identifiable based on the experimental data. Thus, an identifiability analysis is a prerequisite before any statistical method is applied to estimate the unknown parameters from the experimental data.

Three main categories of identifiability techniques have been developed for general ODE models. The first is structural (theoretical) identifiability analysis, which can be used to evaluate whether all parameters can be theoretically identified by manipulating the model structure. Two assumptions are needed for such analysis: 1) the model structure is absolutely accurate; and 2) measurement is exact (no measurement error). Although these two assumptions are not realistic in practice, it is still necessary to study theoretical identifiability. The second type of identifiability analysis is practical identifiability analysis in which both model uncertainty and practical measurement errors are considered. Structural identifiability analysis can be done before experiments for data collection are designed. In fact, the structural identifiability analysis can provide useful information, such as the minimum number of measurements at distinct time points, for experimental design. If an ODE model turns out to be unidentifiable or only a subset of model parameters are identifiable via structural identifiability analysis, the model may need to be modified or some of the parameters may need to be fixed before statistical methods are applied to estimate the unknown parameters. Otherwise, statistical estimates may not be reliable. Even though some parameter estimates can be obtained from an unidentifiable model, the estimates may be a local estimates or an arbitrary set of estimates that can overfit the observation data. If the structural identifiability analysis confirms that an ODE model is globally or locally identifiable, practical identifiability analysis should be done to check the reliability and sensitivity of estimates to measurement errors and model uncertainty. Based on the results of practical identifiability analysis, a model can be further re-fined by model selection techniques [71]. Practical identifiability analysis can also be used to better design future experiments. The third type of identifiability analysis technique is based on the sensitivity matrix. Similar to structural identifiability analysis, sensitivity-based methods do not require experimental observations and cannot account for model uncertainty and measurement errors. Like practical identifiability analysis, sensitivity-based methods also require at least one nominal value of parameters. Note that so far it is still difficult to do structural identifiability analysis for high-dimensional ODEs or complicated ODEs. In this case, the practical identifiability analysis may not be reliable since the structural (theoretical) identifiability of the model is unknown. The class of sensitivity-based methods, which is a technique between the structural (theoretical) identifiability and practical identifiability analysis, can be used in such case.

Besides the identifiability analysis techniques for ODEs, the identifiability analysis of delay differential equation (DDE) models should also be discussed. A general form of a DDE system is given as follows

(7.1)

(7.2)

(7.3)

where t₀ is the starting value of the independent variable, x(t) [set membership]

R^m is a vector of state variables, y(t) [set membership]

R^d the measurable system output vector, u(t) [set membership]

R^p the known system input vector, θ [set membership]

R^q the parameter vector, M a constant coefficient matrix (or mass matrix), and τ(t, x) a vector of delay functions. It is required that

(7.4)

that is, the value of delay functions should always be smaller than or equal to the current time which is reasonable since the future value is unknown yet. Another important assumption is that

(7.5)

which implies that none of the d system outputs is trivial (e.g., linear combination of other outputs). Since τ(t, x) ≤ t, it is necessary to know the value of x(t) when t ≤ t₀, i.e., the history function:

(7.6)

where x(t), t ≤ t₀ is only a function of time and parameters. It is obvious that the DDE system described here is much more complicated than an ODE system, and currently it is still impossible to numerically solve a very general DDE model. However, for some relatively simple DDE models, a number of numerical methods have been proposed and implemented. For details of such algorithms, the reader is referred to the work by Ascher and Petzold [6], Bellen and Zennaro [11], Guglielmi and Hairer [41, 42], and Shampine and Thompson [97]. Particularly, Guglielmi and Hairer [41, 42] proposed and implemented a comparatively general solver (called Radau IIA) for DDE models, which is recommended for practical applications due to its efficiency and stability. For examples of DDE modeling of HIV infection, the reader is referred to [90] and [76]. A number of independent studies have tackled the identifiability problem for DDE models [4, 5, 27, 31, 66, 68, 79, 124, 125]. However, most of these previous works deal with very simple and specific DDE models (e.g., [68]) and the generalizability of these results are limited due to a lack of understanding of the important feature of DDE models: the propagation of discontinuity at t₀ or in the history functions from lower to higher order derivatives of state variables if there is any. Such a feature makes DDE models easily become bifurcated or just unsolvable [11, 41, 42]. The identifiability conclusions based on model structure manipulation are not reliable unless analytical solutions of DDE models can be obtained and analyzed. However, it is surprising that almost all the existing work tried to tackle the identifiability problem by manipulating model structures (e.g., [5]). Generally, for complicated systems as described in Eqs. (7.1) to (7.3), it is extremely difficult to manipulate the model structures to study identifiability problems. Thus, the methodologies for DDE model identifiability are still in its infancy and promising approaches are likely to be numerical methods such as the practical or sensitivity-based methods (e.g., [9]), although it may require development and reliable realization of DDE numerical algorithms.

In addition, it should be mentioned that some of the identifiability techniques such as the differential algebra method can be extended to study the identifiability of partial differential equation (PDE) models (e.g., [49]), but the identifiability analysis for more complicated models such as PDEs or stochastic differential equations is beyond the scope of this paper.

Finally, after the identifiability analysis is done, statistical estimation methods should be used to estimate the unknown parameters in the model. The practical identifiability analysis also requires reliable parameter estimation methods be available. Recently statistical estimation methods for ODE models have attracted a great deal of attention from statisticians. Some novel and efficient estimation methods particular for nonlinear ODE models have been published in the statistical literature [17, 18, 40, 48, 61, 62, 88]. Besides the standard least squares approach [70, 71], more reliable and computationally efficient estimation methods and their theoretical foundations have been developed [17, 18, 62, 88]. However, the topic of statistical estimation methods for ODE, DDE and PDE models is beyond the scope of this paper.

Acknowledgments

This work was partially supported by NIAID/NIH research grants AI055290, AI50020, AI28433, AI078498, RR06555, University of Rochester Provost Award, and University of Rochester DCFAR (P30AI078498) Mentoring Award.

References

1. Adams BM, Banks HT, Kwon HD, Tran HT. Dynamic multidrug therapies for HIV: Optimal and STI control approaches. Mathe Biosci Eng. 2004;1:223–241. [PubMed]

2. Aluru S. Haendbook of Computational Molecular Biology (Chapman & All/Crc Computer and Information Science Series) Chapman & Hall/CRC; Boca Raton, FL: 2005.

3. Anderson DH. Compartmental modeling and tracer kinetics. In: Levin S, editor. Lecture Notes in Biomathematics. Vol. 50. Springer; Berlin: 1983.

4. Anguelova M. PHD thesis. Department of Mathematical Sciences, Chalmers University of Technology and Göteborg University; Göteborg, Sweden: 2007. Observability and identifiability of nonlinear systems with applications in biology.

5. Anguelova M, Wennberg B. State elimination and identifiability of the delay parameter for nonlinear time-delay systems. Automatica. 2008;44:1373–1378.

6. Ascher UM, Petzold LR. The numerical solution of delay-differential-algebraic equations of retarded and neutral type. SIAM J Numer Anal. 1995;32:1635–1657.

7. Audoly S, Bellu G, D’Angio L, Saccomani MP, Cobelli C. Global identifiability of nonlinear models of biological systems. IEEE Trans Biomed Eng. 2001;48:55–65. [PubMed]

8. Baccam P, Beauchemin C, Macken CA, Hayden FG, Perelson AS. Kinetics of influenza A virus infection in humans. J Virol. 2006;80:7590–7599. [PMC free article] [PubMed]

9. Banks HT, Bortz DM. A parameter sensitivity methodology in the context of HIV delay equation models. J Math Biol. 2005;50:607–625. [PubMed]

10. Bellman R, Åström KJ. On structural Identifiability. Math Biosci. 1970;7:329–339.

11. Bellen A, Zennaro M. Numerical solution of delay differential equations. Oxford University Press; Oxford: 2003.

12. Bellu G, Saccomani MP, Audoly S, D’Angi L. DAISY: A new software tool to test global identifiability of biological and physiological systems. Comput Methods Programs Biomed. 2007;88:52–61. [PMC free article] [PubMed]

13. Cacuci DG. Sensitivity and Uncertainty Analysis: Theory. Vol. 1. Chapman Hall/CRC; Boca Raton, FL: 2003.

14. Carra’Ferro G, Gerdt VP. Improved Kolchin-Ritt algorithm. Program Comput Soft. 2003;29:83–87.

15. Castillo-Chavez C, Blower S, vande Driessche P, Kirschner D, Yakubu AA. Mathematical approaches for emerging and reemerging infectious diseases: models, methods, and theory. Springer; New York: 2002.

16. Chappel MJ, Godfrey KR. Structural identifiability of the parameters of a nonlinear batch reactor model. Math Biosci. 1992;108:245–251. [PubMed]

17. Chen J, Wu H. Estimation of time-varying parameters in deterministic dynamic models with application to HIV infections. Statistica Sinica. 2008a;18:987–1006.

18. Chen J, Wu H. Efficient local estimation for time-varying coefficients in deterministic dynamic models with applications to HIV-1 dynamics. J Amer Stat Assoc. 2008b;103:369–384.

19. Cobelli C, Romanin-Jacur G. Controllability, observability and structural identifiability of multi input and multi output biological compartmental systems. IEEE Trans Biomed Eng. 1976;23:93–100. [PubMed]

20. Cobelli C, Polo A, Romanin-Jacur G. A computer program for the analysis of controllability, observability and structural identifiability of biological compartmental systems. Comput Programs Biomed. 1977;7:21–36. [PubMed]

21. Cobelli C, Lepschy A, Jacur R. Identifiability of compartmental systems and related structural proerties. Math Biosci. 1979;44:1–18.

22. Coffin JM. HIV Population Dynamics in Vivo: Implications for Genetic Variation, Pathogenesis, and Therapy. Science. 1995;267:483–489. [PubMed]

23. Davidson EJ. Connectability and structural controllability of composite systems. Automatica. 1997;13:109–123.

24. Degenring D, Froemel C, Dikta G, Takors R. Sensitivity analysis for the reduction of complex metabolism models. J Proc Cont. 2004;14:729–745.

25. Denis-Vidal L, Joly-Blanchard G. An easy to check criterion for (un)identifiability of uncontrolled systems and its applications. IEEE Trans Auto Cont. 2000;45:768–771.

26. Denis-Vidal L, Joly-Blanchard G, Noiret C. Some effective approaches to check the identifiabiity of uncontrolled nonlinear systems. Math Comput Simul. 2001;57:35–44.

27. Denis-Vidal L, Jauberthie C, Joly-Blanchard G. Identifiability of a nonlinear delayed-differential aerospace model. IEEE Trans Auto Cont. 2006;51:154–158.

28. Deskins WE. Abstract algebra. Dover: 1996.

29. Diop S, Fliess M. On nonlinear observability. Proc. of the first Europ. Cont. Conf; Paris, Hermes. 1991. pp. 152–157.

30. Dougherty EP, Hwang JT, Rabitz H. Further developments and applications of the Green’s function method of sensitivity analysis in chemical kinetics. J Chem Phy. 1979;71:1794–1808.

31. Ferretti G, Maffezzoni C, Scattolini R. On the identifiability of the time delay with least squares methods. Automatica. 1996;32:449–453.

32. Filter RA, Xia X, Gray CM. Dynamic HIV/AIDS parameter estimation with application to a vaccine readiness study in Southern Africa. IEEE Trans Biomed Eng. 2005;52:784–791. [PubMed]

33. Frieden BR. Science from Fisher Information. Cambridge University Press; New York: 2004.

34. Froemel C. Master thesis. Fachhochschule Aachen, Abteilung Juelich, Fachbereich Physikalische Technik, Studienrichtung Technomathematik; 2003. Parameterreduktion in stoffwechselmodellen mit methoden der statistik.

35. Gadkar KG, Gunawan R, Doyle FJ. Iterative approach to model identification of biological networks. BMC Bioinformatics. 2005;6:155–175. [PMC free article] [PubMed]

36. Gilmore R. Lie groups, Lie algebras, and some of their applications. Dover: 2006.

37. Glad ST. Solvability of differential algebraic equations and inequalities: an algorithm. Europ. Cont. Conf; Brussels. 1997.

38. Grantham WJ, Vincent TL. Modern control systems analysis and design. Wiley; New York: 1993.

39. Grewal M, Glover K. Identifiability of linear and nonlinear dynamical systems. IEEE Trans Autom Cont. 1976;21:833–837.

40. Guedj J, Thiebaut R, Commenges D., 3 Practical identifiability of HIV dynamics models. Bull Math Biol. 2007;69:2493–2513. [PubMed]

41. Guglielmi N, Hairer E. Implementing Radau IIA methods for stiff delay differential equations. Computing. 2001;67:1–12.

42. Guglielmi N, Hairer E. Computing breaking points in implicit delay differential equations. Adv Comput Math. 2008;29:229–247.

43. Hancioglu B, Swigon D, Clermont G. A dynamical model of human immune response to influenza A virus infection. J Theor Biol. 2007;246:70–86. [PubMed]

44. Hermann R, Krener AJ. Nonlinear controllability and observability. IEEE Trans Autom Cont. 1977;22:728–740.

45. Ho DD, Neumann AU, Perelson AS, Chen W, Leonard JM, Markowitz M. Rapid turnover of plasma virions and CD4 lymphocytes in HIV-1 infection. Nature. 1995;373:123–126. [PubMed]

46. Holmberg A. On the practical identifiability of microbial growth models incorporating Michaelis-Menten type nonlinearities. Math Biosci. 1982;62:23–43.

47. Hsu K, Novara C, Vincent T, Milanese M, Poolla K. Parametric and nonparametric curve fitting. Automatica. 2006;42:1869–1873.

48. Huang Y, Liu D, Wu H. Hierarchical Bayesian methods for estimation of parameters in a longitudinal HIV dynamic system. Biometrics. 2006;62:413–423. [PMC free article] [PubMed]

49. Hubert E. Essential components of algebraic differential equations. J Symb Comput. 1999;28:657–680.

50. Isidori A. Nonlinear control systems. 3. Springer; New York: 1995.

51. Jacquez JA. Compartmental analysis in biology and medicine. 2. The University of Michigan Press; Ann Arbor, Michigan: 1985.

52. Jacquez JA, Greif P. Numerical parameter identifiability and estimability: Integrating identifiability, estimability, and optimal sampling design. Math Biosci. 1985;77:201–227.

53. Jaulin L, Kieffer M, Didrit O, Walter E. Applied Interval Analysis. Springer-Verlag; London: 2001.

54. Jolliffe IT. Discarding variables in a principal component analyssis. I: artificial data. J Royal Stat Soc. 1972;21:160–172.

55. Jeffrey AM, Xia X, Craig IK. When to initiate HIV therapy: a control theoretic approach. IEEE Trans Biomed Eng. 2003;50:1213–1220. [PubMed]

56. Jeffrey AM, Xia X. Identifiability of HIV/AIDS model. In: Tan WY, Wu H, editors. Deterministic and stochastic models of AIDS epidemics and HIV infections with intervention. World Scientific; 2005.

57. Kolchin E. Differential algebra and algebraic groups. Academic Press; New York: 1973.

58. Krakovska O, Wahl LM. Optimal drug treatment regimens for HIV depend on adherence. J Theoret Biol. 2007;246:499–509. [PubMed]

59. Lauffenburger DA, Linderman JL. Receptors: models for binding, trafficking, and signaling. Oxford University Press; Oxford: 1993.

60. Levine WS. The control handbook. CRC Press; 1996.

61. Li L, Brown MB, Lee KH, Gupta S. Estimation and inference for s spline-enhanced population pharmacokinetic model. Biometrics. 2002;58:601–611. [PubMed]

62. Liang H, Wu H. Parameter estimation for differential equation models using a framework of measurement error in regression model. J Amer Stat Assoc. 2008;103:1570–1583. [PMC free article] [PubMed]

63. Ljung L. Convergence analysis of parametric identification methods. IEEE Trans Autom Cont. 1978;23:770–783.

64. Ljung L, Glad T. On global identifiability for arbitrary model parametrizations. Automatica. 1994;30:265–276.

65. Ljung L. System identification: theory for the user. 2. Prentice Hall; New Jersey: 1999.

66. Lunel SMV. Parameter identifiability of differential delay equations. Int J Adapt Cont Sig Proc. 2001;15:655–678.

67. Mansfield EL, Clarkson PA. Applications of the differential algebra package diffgrob2 to classical symmetries of differential equations. J Symb Comput. 1997;23:517–533.

68. Merino JA, Biasi JD, Plusquellec Y, Houin G. Local identifiability for two and three-compartment pahrmacokinetic models with time-lags. Med Eng Phy. 1998;20:261–268. [PubMed]

69. Metropolis N, Ulam S. The Monte Carlo Method. J Amer Stat Assoc. 1949;44:335–341. [PubMed]

70. Miao H, Dykes C, Demeter LM, Perelson AS, Wu H. Modeling and estimation of kinetic parameters and replication fitness of HIV-1 from Flow-Cytometry-based growth competition experiments. Bull Math Biol. 2008;70:1749–1771. [PubMed]

71. Miao H, Dykes C, Demeter LM, Wu H. Differential equation modeling of HIV viral fitness experiments: Model identification, model Selection, and multi-model inference. Biometrics. 2009;65:292–300. [PMC free article] [PubMed]

72. Mishra B. Algorithmic algebra. Springer; New York: 1993.

73. Munkres JR. Analysis on Manifolds. Addison-Wesley; 1991.

74. Murray JD. Mathematical biology I: An introduction. Springer; New York: 2002.

75. Murray JD. Mathematical biology II: Spatial models and biomedical applications. Springer; New York: 2003.

76. Nelson PW, Perelson AS. Mathematical analysis of delay differential equation models of HIV-1 infection. Math Biosci. 2002;179:73–94. [PubMed]

77. Nowak MA, May RM. Virus dynamics: mathematical principles of immunology and virology. Oxford University Press; Oxford: 2000.

78. Ollivier F. PhD Thesis. école Polytechnique; Paris, France: 1990. Le problème de l’identifiabilité globale: étude thé orique, méthodes effectives et bornes de complexité

79. Orlov Y, Belkoura L, Richar JP, Dambrine M. On identifiability of linear time-delay systems. IEEE Trans Auto Cont. 2002;47:1319–1324.

80. Perelson AS, Essunger P, Cao Y, et al. Decay characteristics of HIV-1-Infected compartments during combination therapy. Nature. 1997;387:188–191. [PubMed]

81. Perelson AS. Modelling viral and immune system dynamics. Nat Rev Immunol. 2002;2:28–36. [PubMed]

82. Perelson AS, Nelson PW. Mathematical analysis of HIV-1 dynamics in vivo. SIAM Rev. 1999;41:3–44.

83. Pia Saccomani M, Audoly S, D’Angi L. Parameter identifiability of nonlinear systems: the role of initial conditions. Automatica. 2003;39:619–632.

84. Pohjanpalo H. System identifiability based on the power series expansion of the solution. Math Biosci. 1978;41:21–33.

85. Quaiser T, Marquardt W, Mönnigmann M. Local identifiability analysis of large signaling pathway models. In: Allgoewer F, Reuss M, editors. 2nd Conference Foundation of Systems Biology in Engineering, Proceedings Plenary and Contributed Papers. Fraunhofer IRB Verlag; Stuttgart, Germany: 2007. pp. 465–470.

86. Quaiser T, Mönnigmann M. Systematic identifiability testing for unambiguous mechanistic modeling - application to JAK-STAT, MAP kinase, and NF-κB signaling pathway models. BMC Syst Biol. 2009;3:50–71. [PMC free article] [PubMed]

87. Raksanyi A, Lecourtier Y, Walter E, Venot E. Identifiability and distinguishability testing in computer algebra. Math Biosci. 1985;77:245–266.

88. Ramsay JO, Hooker G, Campbell D, Cao J. Parameter estimation for differential equations: A generalized smoothing approach (with discussions) J Roy Stat Soc B. 2007;69:741–796.

89. Rao CR. Information and the accuracy attainable in the estimation of statistical parameters. Bull Calcutta Math Soc. 1945;37:81–91.

90. Rebecca VC, Ruan S. A delay-differential equation model of HIV infection of CD4+ T-cells. Math Biosci. 2000;165:27–39. [PubMed]

91. Ritt JF. Differential Algebra. American Mathematical Society; Providence, RI: 1950.

92. Rodriguez-Fernandez M, Egea JA, Banga JR. Novel metaheuristic for parameter estimation in nonlinear dynamic biological systems. BMC Bioinformatics. 2006a;7:483–500. [PMC free article] [PubMed]

93. Rodriguez-Fernandez M, Mendes P, Banga JR. A hybrid approach for efficient and robust parameter estimation in biochemical pathways. Biosystems. 2006b;83:248–265. [PubMed]

94. Saltelli A, Chan K, Scott M. Wiley Series in Probability and Statistics. John Wiley & Sons; New York: 2000. Sensitivity analysis.

95. Seigneur C, Stephanopoulos G, Carr RW. Dynamic sensitivity analysis of chemical reaction system: a variational method. Chem Eng Sci. 1982;37:845–853.

96. Schittkowski K. Experimental design tools for ordinary and algebraic differential equations. Ind Eng Chem Res. 2007;46:9137–9147.

97. Shampine LF, Thompson S. Solving DDEs in MATLAB. App Num Math. 2001;37:441–458.

98. Sontag ED. On the observability of polynomial systems. SIAM J Cont Opt. 1979;17:139–151.

99. Sontag ED. Spaces of observables in nonlinear control. in Proc Int Congress of Mathematicians, Birkhäuser Verlag, Basel. 1995;2:1532–1545.

100. Sontag ED. For differential equations with r parameters, 2r + 1, experiments are enough for identification. J Nonlin Sci. 2002;12:553–583.

101. Sontag ED, Wang Y, Megretski A. Input classes for identification of bilinear systems. IEEE Trans Autom Cont. 2009;54:195–207.

102. Takens F. Detecting strange attractors in turbulence. In: Rand DA, Young LS, editors. Proc Symp Dyna Syst Turbulence. Springer-Verlag; Berlin: 1981.

103. Thomaseth K, Cobelli C. Generalized sensitivity functions in physiological system identification. Ann Biomed Eng. 1999;27:607–616. [PubMed]

104. Tunali T, Tarn TJ. New results for identifiability of nonlinear systems. IEEE Trans Autom Contorl. 1987;32:146–154.

105. Turányi T. Sensitivity analysis of complex kinetic systems. Tools and applications. J Math Chem. 1990;5:203–248.

106. Vajda S, Godfrey K, Rabitz H. Similarity transformation approach to identifiability analysis of nonlinear compartmental models. Math Biosci. 1989a;93:217–248. [PubMed]

107. Vajda S, Rabitz H. State isomorphism approach to global identifiability of nonlinear systems. IEEE Trans Autom Cont. 1989;34:220–223.

108. Vajda S, Rabitz H, Walter E, Lecourtier Y. Qualitative and quantitative identifiability analysis of nonlinear chemical kinetiv-models. Chem Eng Commun. 1989b;83:191–219.

109. Van Trees HL. Detection, estimation, and modulation theory, Part I. Wiley; New York: 1968.

110. Vidal R, Chiuso A, Soatto S. Observability and identification of jump linear systems. Proceedings of 41st IEEE conference on decision and control; 2002. pp. 3614–3619.

111. Vincent TL, Grantham WJ. Nonlinear and optimal control systems. Wiley-Interscience; New York: 1997.

112. Walter E, Braems I, Jaulin L, Kieffer M. Guaranteed numerical computation as an alternative to computer algebra for testing models for identifiability. Lecture Notes in Computer Science. 2004:124–131.

113. Walter E, Lecourtier Y. Unidentifiable compartmental models: what to do? Math Biosci. 1981;56:1–25.

114. Walter E. Identifiability of parameteric models. Pergamon Press; Oxford: 1987.

115. Whitney H. Differentiable manifolds. Ann Math. 1936;37:645–680.

116. Whitney H. The self-intersections of a smooth n-manifold in 2n-space. Ann Math. 1944;45:220–246.

117. Wu H, Zhu H, Miao H, Perelson AS. Identifiability and statistical estimation of dynamic parameters in HIV/AIDS dynamic models. B Math Biol. 2008;70:785–799. [PubMed]

118. Xia X. Estimation of HIV/AIDS parameters. Automatica. 2003;39:1983–1988.

119. Xia X, Moog CH. Identifiability of nonlinear systems with applications to HIV/AIDS models. IEEE Trans Autom Contorl. 2003;48:330–336.

120. Yao KZ, Shaw BM, Kou B, McAuley KB, Bacon DW. Modeling ethylene/butene copolymerization with multi-site catalysts: Parameter estimability and experimental design. Polym React Eng. 2003;11:563–588.

121. Yue H, Brown M, Knowles J, Wang H, Broomhead DS, Kell DB. Insights into the behaviour of systems biology models from dynamic sensitivity and identifiablity analysis: a case study of an NF-kappaB signalling pathway. Mol BioSyst. 2006;2:640–649. [PubMed]

122. Zak DE, Gonye GE, Schwaber JS, Doyle FJ. Importance of input perturbations and stochastic gene expression in the reverse engineering of genetic regulatory networks: insights from an identifiability analysis of an in silico network. Genome Res. 2003;13:2396–2405. [PubMed]

123. Zazworsky RM, Knudsen HK. Comments on “Controllability, observability and structural identifiability of multi input and multi output biological compartmental systems. IEEE Trans Biomed Eng. 1977;24:495–496. [PubMed]

124. Zhang J, Xia X, Moog CH. Parameter identifiability of nonlinear systems with time-delay. IEEE Trans Auto Cont. 2006;51:371–375.

125. Zhang J, Xia X. Identifiability problems of time-delay HIV models. 17th IFAC World Congress; Seoul, Korea. 2008.

126. Zurakowski R, Teel AR. A model predictive control based scheduling method for HIV therapy. J Theoret Biol. 2006;238:368–382. [PubMed]