Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC3140286

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Definitions
- 3. Structural Identifiability Analysis
- 4. Practical Identifiability Analysis
- 5. Sensitivity-Based Identifiability Analysis
- 6. Application Examples
- 7. Discussion and Conclusion
- References

Authors

Related links

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/090757009PMCID: PMC3140286

NIHMSID: NIHMS293681

HONGYU MIAO: ude.retsehcor.cmru@oaim_uygnoh; XIAOHUA XIA: az.ca.pu.onitsop@aixx; ALAN S. PERELSON: vog.lnal@psa; HULIN WU: ude.retsehcor.tsb@uwh

See other articles in PMC that cite the published article.

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.

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.

A general dynamic system can be expressed as follows:

$$\stackrel{.}{\mathit{x}}(t)=\mathit{f}(t,\mathit{x}(t),\mathit{u}(t),\mathit{\theta}),$$

(2.1)

$$\mathit{y}(t)=\mathit{h}(\mathit{x}(t),\mathit{u}(t),\mathit{\theta}),$$

(2.2)

where ** x**(

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

Let ** θ** = (

$$\stackrel{.}{\mathit{x}}(t)=\mathit{f}(t,\mathit{x}(t),\mathit{u}(t),{\mathit{\theta}}_{c},{\mathit{\theta}}_{t}),$$

(2.3)

$$\mathit{y}(t)=\mathit{h}(\mathit{x}(t),\mathit{u}(t),{\mathit{\theta}}_{c},{\mathit{\theta}}_{t}).$$

(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].

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

Note that in the definition above, ** u**(

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:

Observable: Given an initial state **x**_{0} 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**(

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.

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*.

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

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**, **θ**_{1}) = **y**(**u**, **θ**_{2}) holds if and only if **θ**_{1} = **θ**_{2}.

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*_{0}*-identifiable*.

Locally strongly identifiable (**x**_{0}-identifiable): For an admissible input **u**(t) in the time range of interest [t_{0}, t_{1}] and a given initial state **x**_{0} = **x**(t_{0}) which is independent of **θ** and not an equilibrium point, if there exists an open set Θ^{0} within the parameter space Θ such that for any two different parameter vectors **θ**_{1}, **θ**_{2} Θ^{0}, the solutions **x**(t, **θ**, **u**) exists on [t_{0}, t_{0} + ε] (t_{0} < ε ≤ t_{1} − t_{0}) for both **θ**_{1} and **θ**_{2}, and **y**(t, **θ**_{1}, **x**_{0}, **u**(t)) ≠ **y**(t, **θ**_{2}, **x**_{0}, **u**(t)) on [t_{0}, t_{0} + ε], the system structure is said to be locally strongly identifiable (or **x**_{0}-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:

Structurally identifiable: Let
${C}_{\mathit{u}}^{N}[{t}_{0},{t}_{1}]$ denote the function space expanded by all input functions on [t_{0}, t_{1}] 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^{0} M, Θ^{0} Θ, and
${U}^{0}\subset {C}_{\mathit{u}}^{N}[{t}_{0},{t}_{1}]$ such that the system is locally strongly identifiable at **θ** given **u** for any **x**_{0} M^{0}, **θ** Θ^{0}, and **u** U^{0}.

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]:

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

$$\mathrm{\Phi}=\mathrm{\Phi}(\mathit{\theta},\mathit{u},\stackrel{.}{\mathit{u}},\dots ,{\mathit{u}}^{(k)},\mathit{y},\stackrel{.}{\mathit{y}},\dots ,{\mathit{y}}^{(k)}),\mathrm{\Phi}\in {R}^{q},$$

can be constructed after a finite number of steps of algebraic calculation or differentiation such that Φ = 0 and
$\mathit{det}{\scriptstyle \frac{\partial \mathrm{\Phi}}{\partial \mathit{\theta}}}\ne 0$ hold in the time range of interest [t_{0}, t_{1}] for any (**θ**, **x**_{0}, **u**) in an open and dense subset of
$\mathrm{\Theta}\times M\times {C}_{\mathit{u}}^{N}[{t}_{0},{t}_{1}]$, 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]:

Algebraically identifiable with known initial conditions: If a meromorphic function
$\mathrm{\Phi}=\mathrm{\Phi}(\mathit{\theta},{\mathit{x}}_{0},\mathit{u}({t}_{0}^{+}),\stackrel{.}{\mathit{u}}({t}_{0}^{+}),\dots ,{\mathit{u}}^{(k)}({t}_{0}^{+}),\mathit{y}({t}_{0}^{+}),\stackrel{.}{\mathit{y}}({t}_{0}^{+}),\dots ,{\mathit{y}}^{(k)}({t}_{0}^{+}))$, Φ 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
$\mathit{det}{\scriptstyle \frac{\partial \mathrm{\Phi}}{\partial \mathit{\theta}}}\ne 0$ hold for all (
$\mathit{\theta},{\mathit{x}}_{0},\mathit{u}({t}_{0}^{+}),\stackrel{.}{\mathit{u}}({t}_{0}^{+}),\dots ,{\mathit{u}}^{(k)}({t}_{0}^{+}),\mathit{y}({t}_{0}^{+}),\stackrel{.}{\mathit{y}}({t}_{0}^{+}),\dots ,{\mathit{y}}^{(k)}({t}_{0}^{+})$), where k is a positive integer, (
$\mathit{\theta},{\mathit{x}}_{0},\mathit{u}({t}_{0}^{+}),\stackrel{.}{\mathit{u}}({t}_{0}^{+}),\dots ,{\mathit{u}}^{(k)}({t}_{0}^{+})$) in an open and dense subset of Θ × M × U, (
$\mathit{u}({t}_{0}^{+}),\stackrel{.}{\mathit{u}}({t}_{0}^{+}),\dots ,{\mathit{u}}^{(k)}({t}_{0}^{+})$) and (
$\mathit{y}({t}_{0}^{+}),\stackrel{.}{\mathit{y}}({t}_{0}^{+}),\dots ,{\mathit{y}}^{(k)}({t}_{0}^{+})$) are the derivatives of **u** and **y** at time
${t}_{0}^{+}$, 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.

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.

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

$$\stackrel{.}{\mathit{x}}=\mathbf{A}(t,\mathit{x},{\mathit{\theta}}_{c})\mathit{x}+\mathit{u},$$

(3.1)

$$\mathit{y}=\mathbf{C}({\mathit{\theta}}_{c})\mathit{x},$$

(3.2)

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

$${\mathit{a}}_{k}({t}_{0})={\mathit{y}}^{(k)}({t}_{0}),$$

where *k* denotes the *k ^{th}* derivative of

$$\mathbf{C}\mathit{x}({t}_{0})={\mathit{a}}_{0}({t}_{0}),$$

(3.3)

$$\mathbf{C}\left[\sum _{i=1}^{k}\frac{(k-1)!}{(k-i)!(i-1)!}{\mathbf{A}}^{(k-i)}({t}_{0}){\mathit{x}}^{(i-1)}({t}_{0})+{\mathit{u}}^{(k-1)}({t}_{0})\right]={\mathit{a}}_{k}({t}_{0}),$$

(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

In nature, the power series expansion method is an approach to verify the *x*_{0} 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:

$$\stackrel{.}{\mathit{x}}=\mathbf{A}\mathit{x}+\mathbf{B}\mathit{u},$$

(3.5)

$$\mathit{y}=\mathbf{C}\mathit{x},$$

(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**^{−1} **AP** of **A** such that

$$\stackrel{.}{\mathit{x}}=({\mathbf{P}}^{-1}\mathbf{AP})\phantom{\rule{0.16667em}{0ex}}\mathit{x}+\mathbf{B}\mathit{u},$$

(3.7)

$$\mathit{y}=\mathbf{C}\mathit{x},$$

(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:

$$\stackrel{.}{\mathit{x}}=\mathit{f}(\mathit{x}(t,\mathit{\theta}),\mathit{\theta})+u(t)\mathit{g}\phantom{\rule{0.16667em}{0ex}}(\mathit{x}(t,\mathit{\theta}),\mathit{\theta}),$$

(3.9)

$$\mathit{y}=\mathit{h}(\mathit{x}(t,\mathit{\theta}),\mathit{\theta}).$$

(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.

Structural equivalence: Given two systems of the family in Eq. (3.9) and (3.10), if there exist two parameters **θ**_{1}, **θ**_{2} Θ such that for the same admissible input u(t), the solution of the two system exists for **θ**_{1} and **θ**_{2}, respectively, and the corresponding system outputs are the same, the system with parameter **θ**_{1} is said to be equivalent to the system with parameter **θ**_{2}, denoted by **θ**_{1} ~ **θ**_{2}.

Under the similarity transformation framework, the identifiability problem becomes a system equivalence problem: a system structure is identifiable if no equivalent systems exist for *θ*_{1}, *θ*_{2} Θ and *θ*_{1} ≠ *θ*_{2} [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.

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,

$$\mathit{y}(\mathit{u},{\mathit{\theta}}_{1})=\mathit{y}(\mathit{u},{\mathit{\theta}}_{2})\Rightarrow {\mathit{\theta}}_{1}={\mathit{\theta}}_{2}$$

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

$$\mathit{f}(\mathit{x},{\mathit{\theta}}_{1})=\mathit{f}(\mathit{x},{\mathit{\theta}}_{2})\Rightarrow {\mathit{\theta}}_{1}={\mathit{\theta}}_{2}$$

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

$$\stackrel{.}{x}(t)=\frac{{\mu}_{m}bl(t)x(t)}{{K}_{s}+bl(t)}-{K}_{d}x(t),$$

(3.11)

$$\stackrel{.}{l}(t)=-\frac{{\mu}_{m}l(t)x(t)}{Y({K}_{s}+bl(t))},$$

(3.12)

$$x(0)=0,l(0)=1.$$

(3.13)

Therefore, the right hand side function vector is

$$\mathit{f}(x,l,\mathit{\theta})=\left(\begin{array}{c}{\scriptstyle \frac{{\mu}_{m}\mathit{blx}}{{K}_{s}+bl}}-{K}_{d}x\\ -{\scriptstyle \frac{{\mu}_{m}lx}{Y({K}_{s}+bl)}}\end{array}\right),$$

and from ** f**(

$$\frac{{\mu}_{{m}_{1}}{b}_{1}lx}{{K}_{{s}_{1}}+{b}_{1}l}-{K}_{{d}_{1}}x=\frac{{\mu}_{{m}_{2}}{b}_{2}lx}{{K}_{{s}_{2}}+{b}_{2}l}-{K}_{{d}_{2}}x,$$

(3.14)

$$-\frac{{\mu}_{{m}_{1}}lx}{{Y}_{1}({K}_{{s}_{1}}+{b}_{1}l)}=-\frac{{\mu}_{{m}_{2}}lx}{{Y}_{2}({K}_{{s}_{2}}+{b}_{2}l)}.$$

(3.15)

Solving the two equations above, we have

$${K}_{{d}_{1}}={K}_{{d}_{2}},{\mu}_{{m}_{1}}={\mu}_{{m}_{2}},\frac{{b}_{1}}{{b}_{2}}=\frac{{K}_{{s}_{1}}}{{K}_{{s}_{2}}}=\frac{{Y}_{2}}{{Y}_{1}},$$

which indicates that only parameters (*K _{d}*,

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,

$$\nexists ({\mathit{\theta}}_{1},{\mathit{\theta}}_{2})\in {R}^{q}\times {R}^{q}\phantom{\rule{0.16667em}{0ex}}\mathit{such}\phantom{\rule{0.16667em}{0ex}}\mathit{that}\phantom{\rule{0.16667em}{0ex}}\mathit{y}(\mathit{u},{\mathit{\theta}}_{1})\equiv \mathit{y}(\mathit{u},{\mathit{\theta}}_{2})\phantom{\rule{0.16667em}{0ex}}\mathit{and}{\left|\right|{\mathit{\theta}}_{1}-{\mathit{\theta}}_{2}\left|\right|}_{\infty}\phantom{\rule{0.16667em}{0ex}}>\delta ,$$

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
, that is,

$$\underset{\_}{\mathbb{S}}\subset \mathbb{S}\subset \overline{\mathbb{S}}.$$

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

Walter et al. [112] thought that if
is empty, the model is identifiable; if
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
to help to improve mathematical models by reparameterizing unidentifiable parameters. Thus, the application of the direct test approach is very limited in practice.

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:

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,

$${\ddot{y}}_{1}{\stackrel{.}{y}}_{2}^{2}-5{\stackrel{.}{x}}_{1}{x}_{2}^{2}{\ddot{y}}_{2}-3{\theta}_{1}{y}_{1}{y}_{2}{u}_{1}+{\theta}_{2}{y}_{2}{\ddot{y}}^{2}$$

(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

$$\stackrel{.}{\mathit{x}}(t)-\mathit{f}(t,\mathit{x}(t),\mathit{u}(t),\mathit{\theta})=0,$$

(3.17)

$$\mathit{y}(t)-\mathit{h}(\mathit{x}(t),\mathit{u}(t),\mathit{\theta})=0.$$

(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*_{1}, *v*_{2}, …, *v _{r}*} denote the

Before we introduce more properties of {*v*_{1}, *v*_{2}, …, *v _{r}*}, some definitions and concepts on differential indeterminates and polynomials need to be described. First, the

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

$$v\prec \stackrel{.}{v}\phantom{\rule{0.16667em}{0ex}}(\forall v\in V),{v}_{1}\prec {v}_{2}\iff {\stackrel{.}{v}}_{1}\prec {\stackrel{.}{v}}_{2}\phantom{\rule{0.16667em}{0ex}}(\forall {v}_{1},{v}_{2}\in V)$$

where means ‘ranks lower than’.

Note that for the same indeterminate *v* *V*, the item with a higher degree ranks higher, e.g., *v* *v*^{2}. The following are two examples of ranking:

$$\mathit{u}\prec \stackrel{.}{\mathit{u}}\prec \cdots \prec \mathit{y}\prec \stackrel{.}{\mathit{y}}\prec \cdots \prec \mathit{\theta}\prec \stackrel{.}{\mathit{\theta}}\prec \cdots \prec \mathit{x}\prec \stackrel{.}{\mathit{x}}\prec \cdots ,$$

(3.19)

$$\mathit{u}\prec \mathit{y}\prec \mathit{\theta}\prec \mathit{x}\prec \stackrel{.}{\mathit{u}}\prec \stackrel{.}{\mathit{y}}\prec \stackrel{.}{\mathit{\theta}}\prec \stackrel{.}{\mathit{x}}\prec \cdots .$$

(3.20)

For a given ranking over a differential polynomial ring {*v*_{1}, *v*_{2}, …, *v _{r}*}, the highest ranking derivative in a differential polynomial

Partially reduced: For two differential polynomials P_{1}, P_{2} , let v_{P1} denote the leader of P_{1}, P_{2} is said to be partially reduced with respect to P_{1} if there exists no proper derivative of v_{P1} in P_{2}.

Reduced: For two differential polynomials P_{1}, P_{2} , let v_{P1} denote the leader of P_{1}, P_{2} is said to be reduced with respect to P_{1} if P_{2} is partially reduced with respect to P_{1} and the degree of v_{P1} in P_{2} is less than the degree of v_{P1} in P_{1}.

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

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*_{1}, *A*_{2}, …, *A _{r}*} and

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

- 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

Considering two differential polynomials *P*_{1} and *P*_{2}, assume the leader of *P*_{2} is *v*_{P2} and there exists a proper derivative of *v*_{P2}, e.g.,
${v}_{{P}_{2}}^{(k)}$ for *k* ≥ 1 in *P*_{1}, then *P*_{1} can be partially reduced by *P*_{2} as follows [49]. First, take derivatives of *P*_{2} up to the *k ^{th}* order

$$\begin{array}{c}{\stackrel{.}{P}}_{2}={S}_{{P}_{2}}{\stackrel{.}{v}}_{{P}_{2}}+{R}_{1},\\ {\ddot{P}}_{2}={S}_{{P}_{2}}{\ddot{v}}_{{P}_{2}}+{R}_{2},\\ \cdots \\ {P}_{2}^{(k)}={S}_{{P}_{2}}{v}_{{P}_{2}}^{(k)}+{R}_{k},\end{array}$$

where *S*_{P2} is the separant of *P*_{2}, and *R _{i}*(

$${S}_{{P}_{2}}^{r}{P}_{1}=Q{P}_{2}^{(k)}+\stackrel{\sim}{P},$$

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
${v}_{{P}_{2}}^{(k)}$ 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

$$\begin{array}{c}{A}_{1}(\mathit{u},\mathit{y}),\dots ,{A}_{m}(\mathit{u},\mathit{y}),\\ {B}_{1}(\mathit{u},\mathit{y},{\theta}_{1}),{B}_{2}(\mathit{u},\mathit{y},{\theta}_{1},{\theta}_{2}),\dots ,{B}_{n}(\mathit{u},\mathit{y},{\theta}_{1},\dots ,{\theta}_{q}),\\ {C}_{1}(\mathit{u},\mathit{y},\mathit{\theta},\mathit{x}),\dots ,{C}_{l}(\mathit{u},\mathit{y},\mathit{\theta},\mathit{x}),\end{array}$$

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*_{1}, *B*_{2}, …, *B _{n}*} to verify the structural identifiability:

Assume that no separant or initial of {B_{1}, B_{2}, …, 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].

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]:

Let Φ : R^{d+p+q} → R^{q} denote a function of model parameter **θ** R^{q}, system input **u** R^{n}, system output **y** R^{d} and their derivatives, that is,

$$\mathrm{\Phi}=\mathrm{\Phi}(\mathit{\theta},\mathit{u},\stackrel{.}{\mathit{u}},\dots ,{\mathit{u}}^{(k)},\mathit{y},\stackrel{.}{\mathit{y}},\dots ,{\mathit{y}}^{(k)}),$$

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
$({\mathit{\theta}}_{\ast},{\mathit{u}}_{\ast},{\stackrel{.}{\mathit{u}}}_{\ast},\dots ,{\mathit{u}}_{\ast}^{(k)},{\mathit{y}}_{\ast},{\stackrel{.}{\mathit{y}}}_{\ast},\dots ,{\mathit{y}}_{\ast}^{(k)})\in {R}^{d+p+q}$ such that

$$\mathrm{\Phi}({\mathit{\theta}}_{\ast},{\mathit{u}}_{\ast},{\stackrel{.}{\mathit{u}}}_{\ast},\dots ,{\mathit{u}}_{\ast}^{(k)},{\mathit{y}}_{\ast},{\stackrel{.}{\mathit{y}}}_{x},\dots ,{\mathit{y}}_{\ast}^{(k)})=0\phantom{\rule{0.16667em}{0ex}}\text{and}\left|{\scriptstyle \frac{\partial \mathrm{\Phi}}{\partial {\mathit{\theta}}_{\ast}}}\right|\ne 0.$$

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

$$\mathrm{\Phi}\approx \mathrm{\Phi}({\mathit{\theta}}_{\ast})+(\mathit{\theta}-{\mathit{\theta}}_{\ast}){\scriptstyle \frac{\partial \mathrm{\Phi}}{\partial {\mathit{\theta}}_{\ast}}};$$

since Φ(*θ** _{*}*) = 0 and
${\left({\scriptstyle \frac{\partial \mathrm{\Phi}}{\partial {\mathit{\theta}}_{\ast}}}\right)}^{-1}$ exists (i.e.,
$\left|{\scriptstyle \frac{\partial \mathrm{\Phi}}{\partial {\mathit{\theta}}_{\ast}}}\right|\ne 0$), a unique solution of

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*_{1}, *B*_{2}, …, *B _{n}*} in the characteristic set. However, a more rigorous approach is to verify whether
$\left|{\scriptstyle \frac{\partial \mathrm{\Phi}}{\partial \mathit{\theta}}}\right|\ne 0$ as suggested by the implicit function theorem method, where Φ is generated from

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]

$$\frac{dT}{dt}=\lambda -\rho T-\beta TV,\phantom{\rule{0.38889em}{0ex}}T(0)={T}_{0},$$

(3.21)

$$\frac{{dT}^{\ast}}{dt}=\beta TV-\delta {T}^{\ast},\phantom{\rule{0.38889em}{0ex}}{T}^{\ast}(0)={T}_{0}^{\ast},$$

(3.22)

$$\frac{dV}{dt}=N\delta {T}^{\ast}-cV,\phantom{\rule{0.38889em}{0ex}}V(0)={V}_{0},$$

(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}^{(3)}=(\frac{\stackrel{.}{V}}{V}-\rho -\beta V)[\ddot{V}+(\delta +c)\stackrel{.}{V}+\delta cV]+N\lambda \delta \beta V-\delta c\stackrel{.}{V}-(\delta +c)\ddot{V}.$$

(3.24)

Therefore, we can define

$$f=(\frac{\stackrel{.}{V}}{V}-\rho -\beta V)[\ddot{V}+(\delta +c)\stackrel{.}{V}+\delta cV]+N\lambda \delta \beta V-\delta c\stackrel{.}{V}-(\delta +c)\ddot{V}-{V}^{(3)},$$

(3.25)

$$\mathrm{\Phi}={\left[\begin{array}{ccccc}f& {f}^{(1)}& {f}^{(2)}& {f}^{(3)}& {f}^{(4)}\end{array}\right]}^{T},$$

(3.26)

such that Φ = 0 is automatically satisfied. If
${\scriptstyle \frac{\partial \mathrm{\Phi}}{\partial \mathit{\theta}}}\ne 0$, then *θ* = (*β*, *ρ*, *ν*, *μ*, *η*) can be identified according to the implicit function theorem, where *ν* = *δc*, *μ* = *δ* + *c*, and *η* = *Nλβδ*. That is, if

$$\mathit{Rank}\left[\frac{\partial \mathrm{\Phi}}{\partial \theta}\right]=\mathit{Rank}\left[\begin{array}{ccccc}{\scriptstyle \frac{\partial f}{\partial \beta}}& {\scriptstyle \frac{\partial f}{\partial \rho}}& {\scriptstyle \frac{\partial f}{\partial \eta}}& {\scriptstyle \frac{\partial f}{\partial \mu}}& {\scriptstyle \frac{\partial f}{\partial \nu}}\\ {\scriptstyle \frac{\partial {f}^{(1)}}{\partial \beta}}& {\scriptstyle \frac{\partial {f}^{(1)}}{\partial \rho}}& {\scriptstyle \frac{\partial {f}^{(1)}}{\partial \eta}}& {\scriptstyle \frac{\partial {f}^{(1)}}{\partial \mu}}& {\scriptstyle \frac{\partial {f}^{(1)}}{\partial \nu}}\\ {\scriptstyle \frac{\partial {f}^{(2)}}{\partial \beta}}& {\scriptstyle \frac{\partial {f}^{(2)}}{\partial \rho}}& {\scriptstyle \frac{\partial {f}^{(2)}}{\partial \eta}}& {\scriptstyle \frac{\partial {f}^{(2)}}{\partial \mu}}& {\scriptstyle \frac{\partial {f}^{(2)}}{\partial \nu}}\\ {\scriptstyle \frac{\partial {f}^{(3)}}{\partial \beta}}& {\scriptstyle \frac{\partial {f}^{(3)}}{\partial \rho}}& {\scriptstyle \frac{\partial {f}^{(3)}}{\partial \eta}}& {\scriptstyle \frac{\partial {f}^{(3)}}{\partial \mu}}& {\scriptstyle \frac{\partial {f}^{(3)}}{\partial \nu}}\\ {\scriptstyle \frac{\partial {f}^{(4)}}{\partial \beta}}& {\scriptstyle \frac{\partial {f}^{(4)}}{\partial \rho}}& {\scriptstyle \frac{\partial {f}^{(4)}}{\partial \eta}}& {\scriptstyle \frac{\partial {f}^{(4)}}{\partial \mu}}& {\scriptstyle \frac{\partial {f}^{(4)}}{\partial \nu}}\end{array}\right]=5,$$

(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

Note that with a high dimensional parameter space, like in the example that we are considering here, the matrix ${\scriptstyle \frac{\partial \mathrm{\Phi}}{\partial \theta}}$ can become very complicated to derive and its rank is difficult to evaluate. For example, one element in the matrix (3.27) is

$$\begin{array}{c}{\scriptstyle \frac{\partial {f}^{(4)}}{\partial \mu}}={V}^{(-5)}[\beta {V}^{(5)}{V}^{6}+10\beta \ddot{V}{V}^{(3)}{V}^{5}+5\beta \stackrel{.}{V}{V}^{(4)}{V}^{5}+\rho {V}^{(5)}{V}^{5}+{V}^{(6)}{V}^{5}\\ -6{({V}^{(3)})}^{2}{V}^{4}-8\ddot{V}{V}^{(4)}{V}^{4}-2\stackrel{.}{V}{V}^{(5)}{V}^{4}+12{\ddot{V}}^{3}{V}^{3}+44\stackrel{.}{V}\ddot{V}{V}^{(3)}{V}^{3}\\ +9{\stackrel{.}{V}}^{2}{V}^{(4)}{V}^{3}-78{\stackrel{.}{V}}^{2}{\ddot{V}}^{2}{V}^{2}-32{\stackrel{.}{V}}^{3}{V}^{(3)}{V}^{2}+84{\stackrel{.}{V}}^{4}\ddot{V}V-24{\stackrel{.}{V}}^{6}].\end{array}$$

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*^{(3)}) at five distinct time points *t*_{1}, ···, *t*_{5}. Denote the values of (*V*, , , *V*^{(3)}) at *t* = *t _{i}* as (
${V}_{i},{\stackrel{.}{V}}_{i},{\ddot{V}}_{i},{\stackrel{.}{V}}_{i}^{(3)}$) for

$${\mathrm{\Phi}}^{\ast}={[{f}_{1},{f}_{2},{f}_{3},{f}_{4},{f}_{5}]}^{T}=0.$$

(3.28)

If

$$\left|\frac{\partial {\mathrm{\Phi}}^{\ast}}{\partial \theta}\right|=\left|\begin{array}{ccccc}{\scriptstyle \frac{\partial {f}_{1}}{\partial \beta}}& {\scriptstyle \frac{\partial {f}_{1}}{\partial \rho}}& {\scriptstyle \frac{\partial {f}_{1}}{\partial \eta}}& {\scriptstyle \frac{\partial {f}_{1}}{\partial \mu}}& {\scriptstyle \frac{\partial {f}_{1}}{\partial \nu}}\\ {\scriptstyle \frac{\partial {f}_{2}}{\partial \beta}}& {\scriptstyle \frac{\partial {f}_{2}}{\partial \rho}}& {\scriptstyle \frac{\partial {f}_{2}}{\partial \eta}}& {\scriptstyle \frac{\partial {f}_{2}}{\partial \mu}}& {\scriptstyle \frac{\partial {f}_{2}}{\partial \nu}}\\ {\scriptstyle \frac{\partial {f}_{3}}{\partial \beta}}& {\scriptstyle \frac{\partial {f}_{3}}{\partial \rho}}& {\scriptstyle \frac{\partial {f}_{3}}{\partial \eta}}& {\scriptstyle \frac{\partial {f}_{3}}{\partial \mu}}& {\scriptstyle \frac{\partial {f}_{3}}{\partial \nu}}\\ {\scriptstyle \frac{\partial {f}_{4}}{\partial \beta}}& {\scriptstyle \frac{\partial {f}_{4}}{\partial \rho}}& {\scriptstyle \frac{\partial {f}_{4}}{\partial \eta}}& {\scriptstyle \frac{\partial {f}_{4}}{\partial \mu}}& {\scriptstyle \frac{\partial {f}_{4}}{\partial \nu}}\\ {\scriptstyle \frac{\partial {f}_{5}}{\partial \beta}}& {\scriptstyle \frac{\partial {f}_{5}}{\partial \rho}}& {\scriptstyle \frac{\partial {f}_{5}}{\partial \eta}}& {\scriptstyle \frac{\partial {f}_{5}}{\partial \mu}}& {\scriptstyle \frac{\partial {f}_{5}}{\partial \nu}}\end{array}\right|\ne 0,$$

(3.29)

by the implicit function theorem, there exists a unique solution of *θ*. Assuming that *β* ≠ 0, some algebra shows that the rank of (
${\scriptstyle \frac{\partial {\mathrm{\Phi}}^{\ast}}{\partial \mathit{\theta}}}$) is equal to the rank of

$$\mathrm{\sum}=\left[\begin{array}{ccccc}{V}_{1}({\ddot{V}}_{1}+\mu {\stackrel{.}{V}}_{1})& {\ddot{V}}_{1}+\mu {\stackrel{.}{V}}_{1}& {V}_{1}& {V}_{1}^{2}& {\stackrel{.}{V}}_{1}({V}_{1}^{-1}{\stackrel{.}{V}}_{1}-\rho -\beta {V}_{1})-{\ddot{V}}_{1}\\ {V}_{2}({\ddot{V}}_{2}+\mu {\stackrel{.}{V}}_{2})& {\ddot{V}}_{2}+\mu {\stackrel{.}{V}}_{2}& {V}_{2}& {V}_{2}^{2}& {\stackrel{.}{V}}_{2}({V}_{2}^{-1}{\stackrel{.}{V}}_{2}-\rho -\beta {V}_{2})-{\ddot{V}}_{2}\\ {V}_{3}({\ddot{V}}_{3}+\mu {\stackrel{.}{V}}_{3})& {\ddot{V}}_{3}+\mu {\stackrel{.}{V}}_{3}& {V}_{3}& {V}_{3}^{2}& {\stackrel{.}{V}}_{3}({V}_{3}^{-1}{\stackrel{.}{V}}_{3}-\rho -\beta {V}_{3})-{\ddot{V}}_{3}\\ {V}_{4}({\ddot{V}}_{4}+\mu {\stackrel{.}{V}}_{4})& {\ddot{V}}_{4}+\mu {\stackrel{.}{V}}_{4}& {V}_{4}& {V}_{4}^{2}& {\stackrel{.}{V}}_{4}({V}_{4}^{-1}{\stackrel{.}{V}}_{4}-\rho -\beta {V}_{4})-{\ddot{V}}_{4}\\ {V}_{5}({\ddot{V}}_{5}+\mu {\stackrel{.}{V}}_{5})& {\ddot{V}}_{5}+\mu {\stackrel{.}{V}}_{5}& {V}_{5}& {V}_{5}^{2}& {\stackrel{.}{V}}_{5}({V}_{5}^{-1}{\stackrel{.}{V}}_{5}-\rho -\beta {V}_{5})-{\ddot{V}}_{5}\end{array}\right].$$

(3.30)

As long as det(Σ) ≠ 0, we have
$({\scriptstyle \frac{\partial {\mathrm{\Phi}}^{\ast}}{\partial \mathit{\theta}}})\ne 0$. 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* ^{(3)} is not involved in the matrix Σ, *V* ^{(3)} should exist at any time point. For evaluating *V* ^{(3)} 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

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, 2*q* + 1 experimental observations are enough for identification of the *q* parameters if measurements are absolutely accurate. Sontag [100] explained that the number 2*q*+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 2*q* 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].

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:

$$\mathit{y}(t)=\mathit{h}(\mathit{x}(t),\mathit{u}(t),\mathit{\theta})=\mathit{\epsilon}(\mathit{t}),$$

(4.1)

where ** ε(t)** is measurement error with mean 0 and variance

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
*θ*_{0}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*θ*$$\mathit{ARE}=100\%\times \frac{1}{N}\sum _{i=1}^{N}\frac{\left|{\theta}_{0}^{(k)}-{\stackrel{\sim}{\theta}}_{i}^{(k)}\right|}{\left|{\theta}_{0}^{(k)}\right|},$$(4.2)where ${\theta}_{0}^{(k)}$ is the*k*-th element of*θ*_{0}and ${\stackrel{\sim}{\theta}}_{i}^{(k)}$ 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.

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 ** = [**_{1}, _{2}, ···, * _{q}*]

$$R=\left[\begin{array}{cccc}{r}_{11}({\widehat{\theta}}_{1},{\widehat{\theta}}_{1})& {r}_{12}({\widehat{\theta}}_{1},{\widehat{\theta}}_{2})& \cdots & {r}_{1q}({\widehat{\theta}}_{1},{\widehat{\theta}}_{q})\\ {r}_{21}({\widehat{\theta}}_{2},{\widehat{\theta}}_{1})& {r}_{22}({\widehat{\theta}}_{2},{\widehat{\theta}}_{2})& \cdots & {r}_{1q}({\widehat{\theta}}_{2},{\widehat{\theta}}_{q})\\ \vdots & \phantom{\rule{0.16667em}{0ex}}& \phantom{\rule{0.16667em}{0ex}}& \phantom{\rule{0.16667em}{0ex}}\\ {r}_{q1}({\widehat{\theta}}_{q},{\widehat{\theta}}_{1})& {r}_{q2}({\widehat{\theta}}_{q},{\widehat{\theta}}_{2})& \cdots & {r}_{qq}({\widehat{\theta}}_{q},{\widehat{\theta}}_{q})\end{array}\right],$$

(4.3)

where *r _{ij}* (

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, *σ*^{2}). In this case, for a general dynamic system (2.1) and (2.2), the Fisher information matrix can be given as

$$\mathit{FIM}=\sum _{i=1}^{N}{\left(\frac{\partial {\widehat{\mathit{y}}}_{i}}{\partial \widehat{\mathit{\theta}}}\right)}^{T}{\mathbf{V}}^{-1}\left(\frac{\partial {\widehat{\mathit{y}}}_{i}}{\partial \widehat{\mathit{\theta}}}\right),$$

(4.4)

where the subscript *i* denotes the *i ^{th}* time point of experimental observation,

$$\mathbf{C}={\mathit{FIM}}^{-1}.$$

(4.5)

Finally, the element *r _{ij}* of the correlation matrix can be defined as

$${r}_{ij}=\frac{{\mathbf{C}}_{ij}}{\sqrt{{\mathbf{C}}_{ii}{\mathbf{C}}_{jj}}},i\ne j,$$

(4.6)

$${r}_{ij}=1,i=j.$$

(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.

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,

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 **θ** Θ, **y**(**u**, **θ**) = **y**(**u**, **θ**^{*}) implies **θ** = **θ**^{*}.

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*_{1} ≤ *t*_{2} ≤ ··· ≤ *t _{N}*, then the

$${s}_{ij}({t}_{k})=\frac{\partial {y}_{i}({t}_{k},{\mathit{\theta}}^{\ast})}{\partial {\theta}_{j}},$$

(5.1)

where *y _{i}*(

$${\mathbf{S}}_{dN\times q}=\left[\begin{array}{ccc}{s}_{11}({t}_{1})& \cdots & {s}_{1q}({t}_{1})\\ \cdots & \ddots & \cdots \\ {s}_{d1}({t}_{1})& \dots & {s}_{dq}({t}_{1})\\ \vdots & \vdots & \vdots \\ {s}_{11}({t}_{N})& \cdots & {s}_{1q}({t}_{N})\\ \cdots & \ddots & \cdots \\ {s}_{d1}({t}_{N})& \cdots & {s}_{dq}({t}_{N})\end{array}\right].$$

(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].

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**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 *θ*^{*}

$$\begin{array}{l}{\mathit{y}}_{k}(\mathit{\theta})=\mathit{y}(\mathit{x}({t}_{k}),\mathit{u}({t}_{k}),\mathit{\theta})\\ \approx \mathit{y}(\mathit{x}({t}_{k}),\mathit{u}({t}_{k}),{\mathit{\theta}}^{\ast})+{\frac{\partial \mathit{y}(\mathit{x}({t}_{k}),\mathit{u}({t}_{k}),\mathit{\theta})}{\partial \mathit{\theta}}|}_{\mathit{\theta}={\mathit{\theta}}^{\ast}}\xb7(\mathit{\theta}-{\mathit{\theta}}^{\ast}),\end{array}$$

(5.3)

where *k* = 1, 2, …, *N* denotes the index of the measurement time points. Let *r** _{k}* denote the measurement at

$$\begin{array}{l}\mathit{RSS}(\mathrm{\Delta}\mathit{\theta})=\sum _{k=1}^{N}{\left[{\mathit{r}}_{k}-{\mathit{y}}_{k}({\mathit{\theta}}^{\ast})+{\frac{\partial \mathit{y}(\mathit{x}({t}_{k}),\mathit{u}({t}_{k}),\mathit{\theta})}{\partial \mathit{\theta}}|}_{\mathit{\theta}={\mathit{\theta}}^{\ast}}\xb7\mathrm{\Delta}\mathit{\theta}\right]}^{2}\\ =\sum _{k=1}^{N}{\left[{\frac{\partial \mathit{y}(\mathit{x}({t}_{k}),\mathit{u}({t}_{k}),\mathit{\theta})}{\partial \mathit{\theta}}|}_{\mathit{\theta}={\mathit{\theta}}^{\ast}}\xb7\mathrm{\Delta}\mathit{\theta}\right]}^{2},\end{array}$$

(5.4)

where *r** _{k}* −

$$\mathit{RSS}(\mathrm{\Delta}\mathit{\theta})={(\mathbf{S}\mathrm{\Delta}\mathit{\theta})}^{T}\xb7\mathbf{S}\mathrm{\Delta}\mathit{\theta},$$

(5.5)

where **S** is the sensitivity matrix defined in (5.2). Obviously, the minimum of *RSS*(Δ** θ**) is reached at

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

$$\mathit{corr}({\mathit{S}}_{\xb7i},{\mathit{S}}_{\xb7j})=\frac{\mathit{cov}({\mathit{S}}_{\xb7i},{\mathit{S}}_{\xb7j})}{\sigma ({\mathit{S}}_{\xb7i})\sigma ({\mathit{S}}_{\xb7j})},$$

(5.6)

where *S*_{·} * _{i}* (or

$${c}_{i}^{\mathit{tot}}=\sum _{j=1,j\ne i}^{q}\mid \mathit{corr}({\mathit{S}}_{\xb7i},{\mathit{S}}_{\xb7j})\mid \xb7I(\mid \mathit{corr}({\mathit{S}}_{\xb7i},{\mathit{S}}_{\xb7j})\mid \phantom{\rule{0.16667em}{0ex}}\ge 1-\delta ),$$

(5.7)

where *I* denotes the indicator function, and *δ* (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.

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],

$${s}_{k,j}^{i}=\frac{{\theta}_{j}}{{y}_{i}}\xb7\frac{\partial {y}_{i}({t}_{k},\mathit{\theta})}{\partial {\theta}_{j}}=\frac{\partial ln{y}_{i}({t}_{k},\mathit{\theta})}{\partial ln{\theta}_{j}},$$

(5.8)

where *i* {1, 2, …, *d*} denotes the index of system outputs, *j* {1, 2, …, *q*} the index of parameters and *k* {1, 2, …, *N*} the index of measurement time points. Then the normalized sensitivity matrix for each *y _{i}* is defined as

$${\mathbf{S}}_{N\times q}^{i}=\left[\begin{array}{ccc}{s}_{1,1}^{i}& \cdots & {s}_{1,q}^{i}\\ \vdots & \ddots & \vdots \\ {s}_{N,1}^{i}& \cdots & {s}_{N,q}^{i}\end{array}\right].$$

(5.9)

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

$$e(\mathrm{\Delta}{\theta}_{j})=\sum _{k=1}^{N}\sum _{i=1}^{d}{\left[\frac{{y}_{i}({t}_{k},{\mathit{\theta}}_{-j},{\theta}_{j}+\mathrm{\Delta}{\theta}_{j})-{y}_{i}({t}_{k},\mathit{\theta})}{{y}_{i}({t}_{k},\mathit{\theta})}\right]}^{2},$$

(5.10)

where *θ*_{−}* _{j}* denotes the parameter vector with the

$$os({\theta}_{j})=\frac{\partial e}{\partial ln{\theta}_{j}}=\sum _{k=1}^{N}\sum _{i=1}^{d}{\left[\frac{\partial ln{y}_{i}({t}_{k},\mathit{\theta})}{\partial ln{\theta}_{j}}\right]}^{2}=\sum _{k=1}^{N}\sum _{i=1}^{d}{({s}_{k,j}^{i})}^{2}.$$

(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**^{iT} **S*** ^{i}* are calculated to provide the information for model reduction. Let
${\lambda}_{j}^{i}$ denote the eigenvalues which are ordered non-decreasingly

$$\mid {\lambda}_{1}^{i}\mid \phantom{\rule{0.16667em}{0ex}}\le \phantom{\rule{0.16667em}{0ex}}\mid {\lambda}_{2}^{i}\mid \phantom{\rule{0.16667em}{0ex}}\le \cdots \le \phantom{\rule{0.16667em}{0ex}}\mid {\lambda}_{q}^{i}\mid .$$

(5.12)

Also, let the corresponding eigenvectors be denoted by

$${\mathbf{\Gamma}}^{i}=({\mathit{\gamma}}_{1}^{i},{\mathit{\gamma}}_{2}^{i},\cdots ,{\mathit{\gamma}}_{q}^{i})=\left[\begin{array}{ccc}{\gamma}_{1,1}^{i}& \cdots & {\gamma}_{1,q}^{i}\\ \vdots & \ddots & \vdots \\ {\gamma}_{q,1}^{i}& \cdots & {\gamma}_{q,q}^{i}\end{array}\right].$$

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:$$\begin{array}{l}{m}_{1}=arg\left(\underset{1\le j\le q}{max}\mid {\gamma}_{j,1}\mid \right),\\ {m}_{l}=arg\left(\underset{\begin{array}{c}1\le j\le q\\ j\ne {m}_{1},\cdots ,{m}_{l-1}\end{array}}{max}\mid {\gamma}_{j,l}\mid \right),\text{for}\phantom{\rule{0.16667em}{0ex}}l>1.\end{array}$$(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:^{i}$$\begin{array}{l}{m}_{1}=arg\left(\underset{1\le j\le q}{max}\sum _{h=1}^{q}{\gamma}_{j,h}^{2}\right),\\ {m}_{l}=arg\left(\underset{\begin{array}{c}1\le j\le q\\ j\ne {m}_{1},\cdots ,{m}_{l-1}\end{array}}{max}\sum _{h=1}^{q}{\gamma}_{j,h}^{2}\right),\text{for}\phantom{\rule{0.16667em}{0ex}}l>1.\end{array}$$(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:$$\begin{array}{l}{m}_{0}=arg\left(\underset{1\le j\le q}{max}\mid {\gamma}_{j,q}\mid \right),\\ {m}_{l}=arg\left(\underset{\begin{array}{c}1\le j\le q\\ j\ne {m}_{0},\cdots ,{m}_{l-1}\end{array}}{max}\mid {\gamma}_{j,q-l}\mid \right),\text{for}\phantom{\rule{0.16667em}{0ex}}l>0.\end{array}$$(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].

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 (

$${\mathit{S}}_{h}^{\perp}={\mathit{S}}_{h}-{\mathit{S}}_{h}^{\mathit{proj}}.$$

(5.16)

In the work of Yao et al. [120], the norm
$\left|\right|{\mathit{S}}_{h}^{\perp}\left|\right|$ was proposed to be the measurement of nearly linear dependency between the vector *S** _{h}* and the vector space

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.

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

$$\mathit{RSS}(\mathit{\theta})=\sum _{k=1}^{N}{[{\mathit{r}}_{k}-{\mathit{y}}_{k}(\mathit{\theta})]}^{2},$$

(5.17)

and let *λ*_{1} ≤ *λ*_{2} ≤ ··· ≤ *λ _{q}* denote the eigenvalues of

$$\mathit{RSS}({\mathit{\theta}}^{\ast}+\alpha {\mathit{\gamma}}_{j})\approx \mathit{RSS}({\mathit{\theta}}^{\ast})+\nabla \mathit{RSS}({\mathit{\theta}}^{\ast})\xb7\alpha {\mathit{\gamma}}_{j}+\frac{1}{2}{\alpha}^{2}{\mathit{\gamma}}_{j}^{T}\xb7{\mathbf{S}}^{T}\mathbf{S}\xb7{\mathit{\gamma}}_{j},$$

(5.18)

where *α* is an arbitrary small constant and *RSS*(*θ** ^{*}*) is the gradient of

$$m=arg\left(\underset{1\le h\le q}{max}(\mid {\mathit{\gamma}}_{j,h}\mid )\right).$$

(5.19)

In practice, *λ _{j}* is usually not exactly zero, therefore, a cut-off value

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

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.

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.

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

$$\frac{dT}{dt}=(\rho -{k}_{m}{T}_{m}-{k}_{w}{T}_{w}-{k}_{R}{T}_{m}w)T,$$

(6.1)

$$\frac{{dT}_{m}}{dt}=({\rho}_{m}+{k}_{m}T-{q}_{m}{T}_{w}){T}_{m}+0.25{k}_{R}{T}_{mw}T,$$

(6.2)

$$\frac{{dT}_{w}}{dt}=({\rho}_{w}+{k}_{w}T-{q}_{w}{T}_{m}){T}_{w}+0.25{k}_{R}{T}_{mw}T,$$

(6.3)

$$\frac{{dT}_{mw}}{dt}=({\rho}_{mw}+0.5{k}_{R}T){T}_{mw}+({q}_{m}+{q}_{w}){T}_{m}{T}_{w},$$

(6.4)

where *T*, *T _{m}*,

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}*,

$${y}_{1}=T,\phantom{\rule{0.38889em}{0ex}}{y}_{2}={T}_{m},\phantom{\rule{0.38889em}{0ex}}{y}_{3}={T}_{w},\phantom{\rule{0.38889em}{0ex}}{y}_{4}={T}_{mw}.$$

(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,

$${\stackrel{.}{y}}_{1}=\rho {y}_{1}-{k}_{m}{y}_{1}{y}_{2}-{k}_{w}{y}_{1}{y}_{3}-{k}_{R}{y}_{1}{y}_{4}.$$

(6.6)

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

$${\ddot{y}}_{1}=\rho {\stackrel{.}{y}}_{1}-{k}_{m}{({y}_{1}{y}_{2})}^{(1)}-{k}_{w}{({y}_{1}{y}_{3})}^{(1)}-{k}_{R}{({y}_{1}{y}_{4})}^{(1)},$$

(6.7)

$${y}_{1}^{(3)}=\rho {\ddot{y}}_{1}-{k}_{m}{({y}_{1}{y}_{2})}^{(2)}-{k}_{w}{({y}_{1}{y}_{3})}^{(2)}-{k}_{R}{({y}_{1}{y}_{4})}^{(2)},$$

(6.8)

$${y}_{1}^{(4)}=\rho {y}_{1}^{(3)}-{k}_{m}{({y}_{1}{y}_{2})}^{(3)}-{k}_{w}{({y}_{1}{y}_{3})}^{(3)}-{k}_{R}{({y}_{1}{y}_{4})}^{(3)},$$

(6.9)

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

$$\mathit{Rank}\left[\begin{array}{cccc}{y}_{1}& -{y}_{1}{y}_{2}& -{y}_{1}{y}_{3}& -{y}_{1}{y}_{4}\\ {\stackrel{.}{y}}_{1}& -{({y}_{1}{y}_{2})}^{(1)}& -{({y}_{1}{y}_{3})}^{(1)}& -{({y}_{1}{y}_{4})}^{(1)}\\ {\ddot{y}}_{1}& -{({y}_{1}{y}_{2})}^{(2)}& -{({y}_{1}{y}_{3})}^{(2)}& -{({y}_{1}{y}_{4})}^{(2)}\\ {y}_{1}^{(3)}& -{({y}_{1}{y}_{2})}^{(3)}& -{({y}_{1}{y}_{3})}^{(3)}& -{({y}_{1}{y}_{4})}^{(3)}\end{array}\right]=4.$$

(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*_{1} = *T* are needed to evaluate
${y}_{1}^{(4)}$, and at least four measurements of *y*_{2} = *T _{m}*,

$${\stackrel{.}{y}}_{2}-{k}_{m}{y}_{1}{y}_{2}-0.25{k}_{R}{y}_{1}{y}_{4}={\rho}_{m}{y}_{2}-{q}_{m}{y}_{2}{y}_{3}.$$

(6.11)

Similarly, by taking the higher derivatives of Eq. (6.11), the parameters (*ρ _{m}*,

$$\mathit{Rank}\left[\begin{array}{cc}{y}_{2}& -{y}_{2}{y}_{3}\\ {\stackrel{.}{y}}_{2}& -{({y}_{2}{y}_{3})}^{(1)}\end{array}\right]=2.$$

(6.12)

By the same token, the parameters (*ρ _{w}*,

$$\mathit{Rank}\left[\begin{array}{cc}{y}_{3}& -{y}_{2}{y}_{3}\\ {\stackrel{.}{y}}_{3}& -{({y}_{2}{y}_{3})}^{(1)}\end{array}\right]=2.$$

(6.13)

Finally, *ρ _{mw}* is identifiable if

$$\mathit{Rank}\phantom{\rule{0.16667em}{0ex}}[{y}_{4}]=1.$$

(6.14)

In summary, all parameters (*ρ*, *ρ _{m}*,

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

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

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]:

$$\frac{d}{dt}{T}_{U}(t)=\lambda -\rho {T}_{U}(t)-\eta (t){T}_{U}(t)V(t),$$

(6.15)

$$\frac{d}{dt}{T}_{I}(t)=\eta (t){T}_{U}(t)V(t)-\delta {T}_{I}(t),$$

(6.16)

$$\frac{d}{dt}V(t)=N\delta {T}_{I}(t)-cV(t),$$

(6.17)

where *T _{U}* is the concentration of uninfected target cells,

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*_{1}, *x*_{2} and *x*_{3} denote *T _{U}*,

$$\eta \prec {y}_{2}\prec {y}_{2}\prec \theta \prec {x}_{3}\prec {x}_{2}\prec {x}_{1},$$

(6.18)

where *θ* = [*λ*, *ρ*, *N*, *δ*, *c*]* ^{T}* is the vector of constant unknown parameters. We can eliminate

$${\ddot{y}}_{1}+(\rho +\delta ){\stackrel{.}{y}}_{1}+\delta \rho {y}_{1}-\delta \lambda +\eta (t){y}_{2}({\stackrel{.}{y}}_{1}+\delta {y}_{1}-\lambda )=0,$$

(6.19)

$${\ddot{y}}_{2}+(\delta +c){\stackrel{.}{y}}_{2}+\delta {cy}_{2}-\eta (t){y}_{2}(N\delta {y}_{1}-{\stackrel{.}{y}}_{2}-{cy}_{2})=0.$$

(6.20)

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

$$\eta (t)=\frac{{\ddot{y}}_{1}+(\rho +\delta ){\stackrel{.}{y}}_{1}+\delta \rho {y}_{1}-\delta \lambda}{-{y}_{2}({\stackrel{.}{y}}_{1}+\delta {y}_{1}-\lambda )},$$

(6.21)

or from Eq. (6.20) as

$$\eta (t)=\frac{{\ddot{y}}_{2}+(\delta +c){\stackrel{.}{y}}_{2}+\delta {cy}_{2}}{{y}_{2}(N\delta {y}_{1}-{\stackrel{.}{y}}_{2}-{cy}_{2})}.$$

(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

$$\begin{array}{l}{\ddot{y}}_{1}{y}_{2}{\stackrel{.}{y}}_{2}-{\stackrel{.}{y}}_{1}{y}_{2}{\ddot{y}}_{2}-\delta {y}_{1}{y}_{2}{\ddot{y}}_{2}+\lambda {y}_{2}{\ddot{y}}_{2}-(\delta +c){\stackrel{.}{y}}_{1}{y}_{2}{\stackrel{.}{y}}_{2}\\ +(\rho \delta +\rho +\delta -{\delta}^{2}-\delta c){y}_{1}{y}_{2}{\stackrel{.}{y}}_{2}+{cy}_{2}{\stackrel{.}{y}}_{2}+\rho {c\stackrel{.}{y}}_{1}{{y}_{2}}^{2}+(\rho \delta c-{\delta}^{2}c){y}_{1}{{y}_{2}}^{2}\\ -N\delta {y}_{1}{\ddot{y}}_{1}{y}_{2}+{c\ddot{y}}_{1}{{y}_{2}}^{2}-N\delta (\rho +\delta ){y}_{1}{\stackrel{.}{y}}_{1}{y}_{2}\\ -N{\delta}^{2}\rho {{y}_{1}}^{2}{y}_{2}+N{\delta}^{2}\lambda {y}_{1}{y}_{2}=0.\end{array}$$

(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,

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:

$$\frac{dT}{dt}=-\beta TV,$$

(6.24)

$$\frac{dI}{dt}=\beta TV-\delta I,$$

(6.25)

$$\frac{dV}{dt}=pI-cV,$$

(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_{50}*/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:

$${V}^{(3)}=\left[\ddot{V}+\delta cV+(\delta +c)\stackrel{.}{V}\right]({V}^{-1}\stackrel{.}{V}-\beta V)-\delta c\stackrel{.}{V}-(\delta +c)\ddot{V}.$$

(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.

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

$$\frac{dT}{dt}=-\beta TV,$$

(6.28)

$$\frac{{dI}_{1}}{dt}=\beta TV-{kI}_{1},$$

(6.29)

$$\frac{{dI}_{2}}{dt}={kI}_{1}-\delta {I}_{2},$$

(6.30)

$$\frac{dV}{dt}={pI}_{2}-cV,$$

(6.31)

where *I*_{1} is the number of latent infected epithelial cells that are not yet producing virus and *I*_{2} 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*_{1}, *I*_{2}, *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.

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.4*p* 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.

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

$$\mathbf{M}\stackrel{.}{\mathit{x}}(t)=\mathit{f}(t,\mathit{x}(t),\mathit{x}[\mathit{\tau}(t,\mathit{x})],\stackrel{.}{\mathit{x}}(t),\stackrel{.}{\mathit{x}}[\mathit{\tau}(t,\mathit{x})],\mathit{u}(t),\mathit{\theta}),$$

(7.1)

$$\mathit{y}(t)=\mathit{h}(\mathit{x}(t),\mathit{x}[\mathit{\tau}(t,\mathit{x})],\mathit{u}(t),\mathit{\theta}),$$

(7.2)

$$\mathit{x}({t}_{0})=\mathit{x}({t}_{0},\mathit{u}({t}_{0}),\mathit{\theta}),$$

(7.3)

where *t*_{0} is the starting value of the independent variable, ** x**(

$$\mathit{\tau}(t,\mathit{x})\le t,$$

(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

$$\text{rank}(\frac{\partial \mathit{h}}{\partial \mathit{x}})=d,$$

(7.5)

which implies that none of the *d* system outputs is trivial (e.g., linear combination of other outputs). Since ** τ**(

$$\mathit{x}(t)=\mathit{g}(t,\mathit{\theta}),\phantom{\rule{0.16667em}{0ex}}\text{if}\phantom{\rule{0.16667em}{0ex}}t\le {t}_{0}$$

(7.6)

where ** x**(

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.

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.

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,* 2*r* + 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]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |