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

**|**HHS Author Manuscripts**|**PMC2885820

Formats

Article sections

- Abstract
- 1. Introduction
- 2. A Multi-Stage Smoothing-Based (MSSB) Approach
- 3. The Spline-Enhanced Nonlinear Least Squares (SNLS) Approach
- 4. Simulation Studies
- 5. Viral Dynamic Parameter Estimation for AIDS Patients
- 6. Discussion and Conclusion
- References

Authors

Related links

Ann Appl Stat. Author manuscript; available in PMC 2010 June 15.

Published in final edited form as:

Ann Appl Stat. 2010 March 1; 4(1): 460–483.

doi: 10.1214/09-AOAS290PMCID: PMC2885820

NIHMSID: NIHMS175615

Department of Biostatistics and Computational Biology, University of Rochester School of Medicine and Dentistry, 601 Elmwood Avenue, Box 630, Rochester, New York 14642, Email: ude.retsehcor.tsb@gnailh

See other articles in PMC that cite the published article.

Modeling viral dynamics in HIV/AIDS studies has resulted in deep understanding of pathogenesis of HIV infection from which novel antiviral treatment guidance and strategies have been derived. Viral dynamics models based on nonlinear differential equations have been proposed and well developed over the past few decades. However, it is quite challenging to use experimental or clinical data to estimate the unknown parameters (both constant and time-varying parameters) in complex nonlinear differential equation models. Therefore, investigators usually fix some parameter values, from the literature or by experience, to obtain only parameter estimates of interest from clinical or experimental data. However, when such prior information is not available, it is desirable to determine all the parameter estimates from data. In this paper, we intend to combine the newly developed approaches, a multi-stage smoothing-based (MSSB) method and the spline-enhanced nonlinear least squares (SNLS) approach, to estimate all HIV viral dynamic parameters in a nonlinear differential equation model. In particular, to the best of our knowledge, this is the first attempt to propose a comparatively thorough procedure, accounting for both efficiency and accuracy, to rigorously estimate all key kinetic parameters in a nonlinear differential equation model of HIV dynamics from clinical data. These parameters include the proliferation rate and death rate of uninfected HIV-targeted cells, the average number of virions produced by an infected cell, and the infection rate which is related to the antiviral treatment effect and is time-varying. To validate the estimation methods, we verified the identifiability of the HIV viral dynamic model and performed simulation studies. We applied the proposed techniques to estimate the key HIV viral dynamic parameters for two individual AIDS patients treated with antiretroviral therapies. We demonstrate that HIV viral dynamics can be well characterized and quantified for individual patients. As a result, personalized treatment decision based on viral dynamic models is possible.

The AIDS pandemic continues to pose a serious threat to public health worldwide. AIDS has killed an estimated 2.9 million people with close to 4.3 million people newly infected with HIV in 2006 alone. The total number of people living with HIV has reached its highest level (WHO web site). Although highly active antiretroviral therapy (HAART) regimens are effective in suppressing plasma HIV RNA levels (viral load) below the limit of detection, many patients may fail HAART treatments due to drug resistance, lack of potency, poor drug adherence, pharmacokinetic problems, and adverse effects. In addition, the complexity of regimens and lack of full understanding of the pathogenesis of HIV infection also pose great challenges to AIDS researchers. Over the past 2 decades, many mathematicians and statisticians have developed mechanism-based models and statistical approaches to assist in understanding HIV pathogenesis and have made significant contributions in this area (Ho et al., 1995; Wei et al, 1995, Perelson et al., 1996; Perelson et al., 1997; Wu et al., 1999). In particular, differential equation models have been widely used in describing dynamics and interactions of HIV and the immune system. Some survey on these models can be found in Perelson and Nelson (1999), Nowak and May (2000), and Tan and Wu (2005).

For the mechanism-based models of HIV infection, one critical question is how to use experimental or clinical data to estimate the parameters in the nonlinear differential equation models which do not have closed-form solutions in most cases. Researchers have made a substantial effort to get an approximate closed-form solution under various assumptions, and then use the standard regression approach to estimate the dynamic parameters in the models (Ho et al., 1995; Wei et al, 1995, Perelson et al., 1996; Perelson et al., 1997; Wu et al., 1999; Wu and Ding, 1999; Wu, 2005). But these approximations and assumptions may not always hold, in particular, for patients undergoing long-term treatment.

In this paper, we consider the following HIV dynamic model for patients under long term treatments (e.g., Perelson and Nelson 1999),

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

(1.1)

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

(1.2)

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

(1.3)

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

In practice, some parameters can be fixed from previous studies and only the remaining parameters are needed to be estimated. However, when such prior knowledge is not available, it is important to estimate all viral dynamic parameters, (λ, *ρ, N, δ, c, η*(*t*))^{T}, for each individual patient from the clinical measurements of *V*(*t*) and *T*(*t*), since the estimated dynamic parameters may be used to guide clinical decisions for individualized treatment. Although AIDS investigators have tried to estimate some of these dynamic parameters based on viral load data (Ho et al., 1995; Wei et al, 1995, Perelson et al., 1996; Perelson et al., 1997; Wu et al., 1999), none of them have successfully estimated all these dynamic parameters directly for an individual patient. Huang and Wu (2006) and Huang, Liu and Wu (2006) have made an effort to use the Bayesian method to estimate all these parameters, but strong priors are required for most parameters in order to make all parameters identifiable. In the work of Xia (2003), Filter, Xia and Gray (2005), Gray et al. (2005), Ouattara, Mhawej and Moog (2008), all model parameters were estimated for HIV dynamic models with only constant coefficients. In this paper, to the best of our knowledge, this effort represents the first attempt to simultaneously estimate all the constant and time-varying HIV viral dynamic parameters for individual patients using both viral load and total CD4+ T cell count data.

In this paper we propose two estimation approaches. The first one is a multi-stage smoothing-based (MSSB) approach. We derive the direct relationships between the unknown dynamic parameters and measurement variables from the original differential equation models, and then we employ these relationships to formulate regression models for the unknown parameters. These regression models involve not only the time function of measurement variables, but also their derivatives. We propose using a nonparametric smoothing method, such as the local polynomial smoothing, to obtain the measurement variables and their derivatives which are substituted into the regression models to estimate the unknown dynamic parameters (including the time-varying parameter) in the second step. This approach avoids directly solving the ODEs by numerical methods and is computationally efficient. The second approach that we propose is to use a spline-enhanced nonlinear least squares (SNLS) method. This approach is to use splines to approximate the time-varying parameter so that the original ODE model becomes a model only containing constant parameters. The standard NLS method can be employed to estimate the unknown dynamic parameters and spline coefficients. This approach needs to use a numerical method to repeatedly solve the nonlinear ODEs for the high-dimensional NLS optimization problem which is computationally challenging. Some cutting-edge computational techniques are necessary to solve this problem. Note that, the MSSB approach is computationally efficient but the derived estimates are quite rough; however, these estimates can be used as the initial values for the SNLS method so that the two approaches can be efficiently combined for practical applications.

The rest of this paper is organized as follows. In Section 2, we briefly discuss the identifiability of the HIV dynamic model, then propose the multi-stage smoothing-based (MSSB) approach. In Section 3, we introduce the spline-enhanced nonlinear least squares (SNLS) method. In particular, a hybrid optimization technique combining a gradient method and a global optimization algorithm will be used to tackle the high-dimensional NLS optimization problem. The simulation studies are presented to evaluate the performance of the proposed approaches in Section 4. We apply the proposed methods to estimate HIV viral dynamic parameters (including the time-varying infection rate) from data for two AIDS patients in Section 5. We demonstrate that some of these dynamic parameters are estimated from clinical data for the first time. We conclude the paper with some discussions in Section 6. The details of the proposed hybrid optimization algorithm are given in the Appendix.

Before introducing the estimation methods, it is critical to determine whether all unknown model parameters can be uniquely and reliably identified from a given system input and measurable output. The technique to answer this question is called identifiability analysis. There are mainly two types of identifiability problems: structural (mathematical) identifiability and practical (statistical) identifiability. Structural identifiability analysis techniques use model structure information only (without knowledge of experimental observations) to determine whether all unknown parameters are uniquely identifiable. Therefore, it is also called the prior identifiability analysis. In contrast, practical identifiability analysis techniques should be applied after model fitting to verify the reliability of estimates, so it is called the posterior identifiability analysis (Rodriguez-Fernandez, Egea and Banga 2006). A thorough discussion and review of identifiability techniques and theories is out of the scope of this study, for more details, the interested reader is referred to Ritt (1950), Bellman and Åström (1970), Kolchin (1973), Pohjanpalo (1978), Cobelli et al. (1979), Walter (1987), Vajda et al. (1989), Ollivier (1990), Chappel and Godfrey (1992), Ljung and Glad (1994), Audoly et al. (2001), Xia and Moog (2003), Jeffrey and Xia (2005), Wu et al. (2008), and Miao et al. (2008). In this study, we employed the differential algebra approach (Ljung and Glad, 1994) and it can be easily verified that all constant and time-varying parameters in models (1.1) - (1.3) are structurally identifiable.

Parameter estimation for ODE models has been investigated using the least squares principle by mathematicians (Hemker, 1972; Bard, 1974; Li, Osborne and Prvan, 2005), computer scientists (Varah, 1982), and chemical engineers (Ogunnaike and Ray, 1994; Poyton et al., 2006). Mathematicians have focused on the development of efficient and stable algorithms to solve the least squares problem. Recently statisticians have started to develop various statistical methods to estimate dynamic parameters in ODE models. For example, Putter et al. (2002), Huang and Wu (2006), and Huang, Liu and Wu (2006) have developed hierarchical Bayesian approaches to estimate dynamic parameters in HIV dynamic models for longitudinal data. Li et al. (2002) proposed a spline-based approach to estimate time-varying parameters in ODE models. Ramsay (1996) proposed a technique named principal differential analysis (PDA) for estimation of differential equation models (see a comprehensive survey in Ramsay and Silverman, 2005). Recently Ramsay et al. (2007) applied a penalized spline method to estimate the constant dynamic parameters in ODE models. Chen and Wu (2008, 2009) and Liang and Wu (2008) proposed a two-step smoothing-based approach to estimate both constant and time-varying parameters in ODE models separately.

In this section, we adopt the ideas from Liang and Wu (2008) and Chen and Wu (2008, 2009) to use the smoothing-based approach to estimate all dynamic parameters, including both constant and time-varying parameters, in model (1.1) - (1.3) in three stages. Stage I is to smooth the noisy data to estimate the state variables and their derivatives; Stage II is to estimate constant dynamic parameters (λ, *ρ, c*); Stage III is proposed to estimate both constant dynamic parameters, *δ* and *N*, and the time-varying parameter *η*(*t*). We introduce these three stages in details in the following subsections.

In this subsection we briefly describe how we use the local polynomial regression technique to estimate the functions *T*(*t*) and *V*(*t*) and their derivatives. Consider a general situation that a state variable, *X*(*t*), is observed at *n* time points, (*Y*_{1}, …, *Y _{n}*), i.e.,

$${\sum _{i=1}^{n}\left\{{Y}_{i}-\sum _{j=0}^{p}{\zeta}_{j}{\left({t}_{i}-t\right)}^{j}\right\}}^{2}{K}_{h}\left({t}_{i}-t\right),$$

where *K*(·) is a symmetric kernel function, *K _{h}*(·) =

$${\widehat{X}}^{\left(q\right)}\left(t\right)={\xi}_{\left(q+1\right)}^{\mathrm{T}}{({\mathbf{T}}_{\mathrm{1+}q,t}^{\mathrm{T}}{\mathbf{W}}_{t}{\mathbf{T}}_{1+q,t})}^{-1}{\mathbf{T}}_{1+q,t}^{\mathrm{T}}{\mathbf{W}}_{t}Y,\phantom{\rule{0.2em}{0ex}}\mathrm{for}\phantom{\rule{0.2em}{0ex}}q=0,1,2,$$

where *ξ** _{q}* is the (

We apply the local polynomial smoothing technique to the viral load data, *V*(*t*) and total CD4+ T cell count data, *T*(*t*) = *T _{U}*(

In this subsection we apply the approach proposed by Liang and Wu (2008) to estimate constant dynamic parameters (λ, *ρ, c*) in models (1.1)-(1.3). Notice that we can combine equations (1.1) and (1.2), and obtain

$$\frac{d}{dt}\left[{T}_{U}\left(t\right)+{T}_{I}\left(t\right)\right]=\mathrm{\lambda}-\rho {T}_{U}\left(t\right)-\delta {T}_{I}\left(t\right).$$

Recalling that *T*(*t*) = *T _{I}*(

$$\frac{d}{dt}T\left(t\right)=\mathrm{\lambda}-\rho \left[T\left(t\right)-{T}_{I}\left(t\right)\right]-\delta {T}_{I}\left(t\right).$$

From the above equation, we obtain

$${T}_{I}=\frac{-\mathrm{\lambda}}{\rho -\delta}+\frac{-\rho}{\underset{6}{\rho -\delta}}T+\frac{1}{\rho -\delta}{T}^{\prime},$$

where *T*′ = *dT*(*t*)/*dt*. Substituting this expression into equation (1.3) and letting
${\alpha}_{0}=\frac{N\delta \mathrm{\lambda}}{\rho -\delta},{\alpha}_{1}=\frac{N\delta \rho}{\rho -\delta}\phantom{\rule{0.1em}{0ex}}\mathrm{and}\phantom{\rule{0.1em}{0ex}}{\alpha}_{2}=\frac{N\delta}{\rho -\delta}$,we have

$${V}^{\prime}\left(t\right)={\alpha}_{0}+{\alpha}_{1}T\left(t\right)+{\alpha}_{2}{T}^{\prime}\left(t\right)-cV\left(t\right),$$

(2.1)

where *V*(*t*) and *T*(*t*) for *t* = *t*_{1}, *t*_{2}, …, *t _{n}* are the measurements of viral load and CD4+ T cell count from clinical studies which are subject to measurement errors.

Let (*t*), (*t*), ′(*t*) and ′(*t*) be the estimates of *V*(*t*), *T*(*t*), *V*′(*t*) and *T*′(*t*) from Subsection 2.1, respectively, and substitute these estimates in (2.1). We then obtain a linear regression model (Liang and Wu 2008),

$$\widehat{{V}^{\prime}}\left(t\right)={\alpha}_{0}+{\alpha}_{1}\widehat{T}\left(t\right)+{\alpha}_{2}{\widehat{T}}^{\prime}\left(t\right)-c\widehat{V}\left(t\right)+\epsilon \left(t\right),$$

(2.2)

where *ε*(*t*) includes all substitution errors. Using the least squares regression technique, we obtain the estimates of *α*_{0},*α*_{1}, and *α*_{2} which are termed as the pseudoleast squares (PsLS) estimates in Liang and Wu (2008). Notice that, from the above derivation, we have the relationship λ = −*α*_{0}/*α*_{2} and *ρ* = *α*_{1}/*α*_{2}. Thus, we can obtain the estimates of *c*, λ and *ρ* from model (2.2).

Write *α* = (*α*_{0},*α*_{1},*α*_{2},*c*) and Λ = diag(1, ,′, −). Let
${\mu}_{l}\left(K\right)={\int}_{-1}^{1}{z}^{l}K\left(z\right)dz$ for *l* = 0, 1, 2. Using the arguments similar to those in Liang and Wu (2008), we can establish asymptotics for the smoothing-based estimators of *α*.

The Delta-method can then be used for establishment of the asymptotic properties of the estimators for λ and *ρ*. Notice that, when the mean of *ε*(*t*) is zero, the least squares estimates of *α*_{0}, *α*_{1}, and *α*_{2} are unbiased. However, the substitution error *ε*(*t*) in model (2.2) is not mean zero, instead its mean is in the order of *h*^{2}, and variance is of the order (*nh*)^{−1} which goes to zero when *n* → ∞. Thus, *ε*(*t*) is different from a standard measurement error with mean zero and constant variance,

Now we propose how to estimate the remaining constant parameters, *N* and *δ*, and the time-varying parameter, *η*(*t*). Recalling (1.2), we have

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

(2.3)

while from (1.3), we can obtain

$$\frac{d}{dt}{T}_{I}\left(t\right)=\frac{{V}^{\left(2\right)}\left(t\right)+c{V}^{\prime}\left(t\right)}{N\delta}.$$

(2.4)

A combination of (2.3) and (2.4) deduces

$$\begin{array}{rl}{V}^{\left(2\right)}\left(t\right)+c{V}^{\prime}\left(t\right)=\eta \left(t\right)N\delta T\left(t\right)V\left(t\right)-& \eta \left(t\right)\left\{{V}^{\prime}\left(t\right)V\left(t\right)+c{V}^{2}\left(t\right)\right\}\\ -& \delta \left\{{V}^{\prime}\left(t\right)+cV\left(t\right)\right\}.\end{array}$$

(2.5)

Let *Z*(*t*) = ^{(2)}(*t*)+*ĉ*′(*t*), *U*_{1}(*t*) = −{′(*t*)(*t*)+*ĉ*^{2}(*t*)}, *U*_{2}(*t*) = − {′(*t*)+ *ĉ*(*t*)}, *U*_{3}(*t*) = (*t*)(*t*) in which all functions and parameters have been estimated from Stages I and II. Thus, from (2.5), we can formulate a semiparametric time-varying coefficient regression model (Wu and Zhang, 2006):

$$Z\left(t\right)={U}_{1}\left(t\right)\delta +{U}_{2}\left(t\right)\eta \left(t\right)+{U}_{3}\left(t\right)N\delta \eta \left(t\right)+{\epsilon}^{\ast}\left(t\right),$$

(2.6)

where *ε**(*t*) includes all substitution errors. Note that *δ* and *Nδ* are constant parameters, while *η*(*t*) is an unknown time-varying parameter. To fit this model, we can approximate the time-varying parameter *η*(*t*) by the local polynomial approach, a basis spline method or other nonparametric techniques (Wu and Zhang, 2006). The basis spline approach is straightforward and simple if the function *η*(*t*) does not fluctuate dramatically as in our case of HIV dynamic models. Thus, we select to use the B-spline approximation here, i.e.,

$$\eta \left(t\right)\approx \sum _{j=1}^{s}{a}_{j}{b}_{j},k\left(t\right)$$

(2.7)

where *a _{j}* are constant B-spline coefficients, and

$$\begin{array}{c}Z\left(t\right)={U}_{1}\left(t\right)\delta +\sum _{j=1}^{s}\{{b}_{j,k}\left(t\right){U}_{2}\left(t\right)\}{a}_{j}\\ +\sum _{j=1}^{s}\{{b}_{j,k}\left(t\right){U}_{3}\left(t\right)\}N\delta {a}_{j}+{\epsilon}^{\ast}\left(t\right).\end{array}$$

(2.8)

This becomes a standard linear regression model and we can easily estimate the unknown constant parameters *δ, a _{j}*, and

Also notice that, Wu and Ding (1999) suggested that fitting viral dynamic data to different models for different time periods may result in better estimates of viral dynamic parameters. For example, we may fit a viral dynamic model with a closed-form solution to the early (first week) viral load data to obtain a better estimate of parameters *δ* and *c* as proposed by Perelson et al. (1996). These estimates can then be substituted into the regression models in Stage II and III to obtain the estimates of other parameters. We will adopt this alternative strategy in our real data analysis.

In summary, we are able to use a multi-stage approach to derive the estimates of all dynamic parameters in an HIV dynamic model. This approach only involves parametric and nonparametric/semiparametric regressions and the implementation is straightforward. The numerical evaluations of ODEs are avoided and the initial values of the state variables are not required. However, there are two issues should be addressed for this approach. First, in this study, we focused on models (1.1) - (1.3) and the three stages were therefore particularly proposed for the specific model. However, the MSSB approach itself is a general method for estimating parameters in ODE models; for general guidelines in use of the MSSB approach, the interested reader is referred to Liang and Wu (2008). Second, although the estimators based on the MSSB approach have some attractive asymptotic properties, some limitations exist for this approach. In particular, it needs to estimate the derivatives of state variables which is sensitive to measurement noise when the data are sparse. This may result in biased estimates of parameters and measurement errors. In order to improve the estimates, we propose the spline-enhanced least squares approach in the next section.

In this section, we introduce a spline-enhanced nonlinear least squares (SNLS) method to refine parameter estimates for model (1.1) - (1.3). The basic idea is to approximate the time-varying parameter *η*(*t*) by the B-spline approach. Then model (1.1) - (1.3) becomes the standard ODEs with constant parameters only. We can use the standard NLS approach to estimate the constant dynamic parameters and B-spline coefficients by numerically evaluating the ODEs repeatedly. This method is computationally intensive. Li et al. (2002) used a similar idea to approximate the unknown time-varying parameter in a pharmacokinetic ODE model.

Different types of splines can be constructed based on different basis functions, such as the well-known piecewise polynomial splines and basis splines. Note that, for an arbitrary spline function of a specific degree and smoothness over a given domain partition, a linear combination of basis splines of the same degree and smoothness over the same partition can always be found to represent this spline function (de Boor, 1978). Therefore, B-splines can be employed to approximate the time-varying parameter *η*(*t*) in this study without loss of generality.

Similar to the approximation (2.7) in the previous section, *η*(*t*) can be approximated by a B-spline of order *k* with *s* control points (de Boor 1978), i.e.,
$\eta \left(t\right)\approx {\sum}_{j=1}^{s}{a}_{j}{b}_{j,k}\left(t\right)$. In addition, it should be noted that once the number and positions of control points are determined, the number and positions of knots are also automatically determined at the knot average sites (de Boor, 1978). Thus, our model Eq. (1.1) - (1.3) can be approximated by

$$\frac{d}{dt}{T}_{U}\left(t\right)=\mathrm{\lambda}-\rho {T}_{U}\left(t\right)-\left\{\sum _{j=1}^{s}{a}_{j}{b}_{j,k}\left(t\right)\right\}{T}_{U}\left(t\right)V\left(t\right),$$

(3.1)

$$\frac{d}{dt}{T}_{I}\left(t\right)=\left\{\sum _{j=1}^{s}{a}_{j}{b}_{j,k}\left(t\right)\right\}{T}_{U}\left(t\right)V\left(t\right)-\delta {T}_{I}\left(t\right),$$

(3.2)

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

(3.3)

where *θ* = (λ, *ρ, N, δ, c, a*_{1}, *a*_{2}, …, *a _{s}*)

$${Y}_{1i}=T\left({t}_{i}\right)\phantom{\rule{0.2em}{0ex}}+{\in}_{1i},i=1,2,\cdots ,{n}_{T}$$

(3.4)

$${Y}_{2j}=V({t}_{j})\phantom{\rule{0.2em}{0ex}}+{\in}_{2j},j=1,2,\cdots ,{n}_{V},$$

(3.5)

where *ε*_{1}* _{i}* and

$$\mathit{RSS}\left(\theta \right)=\sum _{i=1}^{{n}_{T}}{\{{Y}_{1i}-T({t}_{i},\theta )\}}^{2}+\sum _{j=1}^{{n}_{V}}{\{{Y}_{2j}-V({t}_{j},\theta )\}}^{2},$$

(3.6)

where *n _{T}* is the total number of CD4+ T cell measurements and

For the SNLS approach, a critical step is to minimize a multimodal objective function (3.6) over a high-dimensional parameter space. In practice, it is challenging to find the global solution to such problems if the parameter space is high-dimensional (as always the case in many biomedical problems), and/or the parameter values are of different orders of magnitude, and/or the objective function is multimodal or even not smooth with noisy data. Thus, it is critical to develop an efficient and stable optimization algorithm.

There are three main categories of optimization methods: direct search methods, gradient methods, and global optimization methods. However, both the direct search methods (e.g., the Simplex method) and the gradient methods (e.g, the Levenberg-Marquardt method and the Gauss-Newton method) can be easily trapped by local minima and even just fail for our problems (Miao et al., 2008). For details and application of such methods in ODE models, the reader is referred to Nocedal and Wright (1999) and Englezos and Kalogerakis (2001). Therefore, the global optimization methods are more suitable for the parameter estimation problem for ODE models. Moles, Banga and Keller (2004) compared the performance and computational cost of seven global optimization methods, including the differential evolution method (Storn and Price 1997). Their results suggest that the differential evolution method outperforms the other six methods with a reasonable computational cost. Unfortunately, existing global optimization methods are very computationally intensive. Improved performance can be achieved by combining gradient methods and global optimization methods, called hybrid methods. A hybrid method based on the scatter search and sequential quadratic programming (SQP) has been proposed by Rodriguez-Fernandez, Egea and Banga (2006), who showed that the hybrid scatter search method is much faster than the differential evolution method for a simple HIV ODE model. However, our preliminary work did not show improved performance of the scatter search method without SQP with respect to the differential evolution method in terms of computational cost and convergence rate. Thus, the performance improvement of the hybrid scatter search method is mainly due to the incorporation of the SQP local optimization method. In addition, Vrugt and Robinson (2007) suggested that a significant improvement of efficiency can be achieved by combining multiple global optimization algorithms. Here we combined the differential evolution and the scatter search method and incorporated with the SQP local optimization method (e.g., SOLNP by Ye (1987)) for parameter estimation of ODE models. We present the details of the proposed algorithm in the Appendix.

Another critical problem that needs to be addressed in order to fit model (3.1) - (3.3) is to determine the order *k* and the number and positions of control points *s* for the spline approximation. In general, we assume that *η*(*t*) is a first-order continuous function of time. Therefore, B-splines of order 2 (piecewise straight lines) should not be considered. Also, since the high order B-spline approximation (e.g., *k* ≥ 5) may introduce unnecessary violent local oscillations (called Runge’s phenomenon) (Runge 1901), we consider up to 4* ^{th}* order B-splines in this study. As to the positions of the control points, to account for highly-skewed data, we select the control points’ positions such that they are equally-spaced in the logarithm time scale. We can use AIC, BIC and AICc criteria (Akaike, 1973; Schwarz, 1978; Burnham and Anderson, 2004) to determine the order and the knots, i.e.,

$$\mathit{AIC}=-2\phantom{\rule{0.1em}{0ex}}\mathrm{In}\phantom{\rule{0.1em}{0ex}}L+2K,$$

(3.7)

$$\mathit{BIC}=-2\phantom{\rule{0.1em}{0ex}}\mathrm{In}\phantom{\rule{0.1em}{0ex}}L+K\phantom{\rule{0.1em}{0ex}}\mathrm{In}\phantom{\rule{0.1em}{0ex}}\left(N\right),$$

(3.8)

$$AI{C}_{c}=\mathit{AIC}+\frac{2K\left(K+1\right)}{N-K-1},$$

(3.9)

where *L* denotes the likelihood function, *N* the total number of observations, and *K* the number of unknown parameters. Under the normality assumption of measurement errors, these model selection criteria can be rewritten as

$$\mathit{AIC}=N\phantom{\rule{0.1em}{0ex}}\mathrm{In}\phantom{\rule{0.1em}{0ex}}\left(\frac{\mathit{RSS}}{N}\right)+2K,$$

(3.10)

$$\mathit{BIC}=N\phantom{\rule{0.1em}{0ex}}\mathrm{In}\phantom{\rule{0.1em}{0ex}}\left(\frac{\mathit{RSS}}{N}\right)+K\phantom{\rule{0.1em}{0ex}}\mathrm{In}\phantom{\rule{0.1em}{0ex}}\left(N\right),$$

(3.11)

$$AI{C}_{c}=N\phantom{\rule{0.1em}{0ex}}\mathrm{In}\phantom{\rule{0.1em}{0ex}}\left(\frac{\mathit{RSS}}{N}\right)+\frac{2NK}{N-K-1},$$

(3.12)

where *RSS* is the residual of sum squares obtained from the NLS model fitting. It should also be noticed that, as the number of spline knots goes to large, the associated parameter space will be high-dimensional. In this case, the novel and efficient optimization methods such as the hybrid optimization algorithms are necessary to locate the global minima of the NLS objective function. In addition, it is necessary to pre-determine a good and informative search range for each of the unknown parameters in order to ease the computational burden of the proposed optimization algorithm. Thus, it is necessary and a good strategy to combine the proposed two approaches to estimate the dynamic parameters in a complex dynamic system. That is, first the MSSB approach introduced in the previous section can be used to obtain a rough estimate and the search range for the unknown parameters, and then the SNLS approach can be used to refine the estimates.

In this section, we evaluate the performance of the MSSB and the SNLS approaches based on Eq. (1.1) - (1.3) by Monte Carlo simulations. An ad hoc bandwidth selection procedure is used for selecting a proper bandwidth *h* for MSSB. See Liang and Wu (2008) for more detailed discussions. The following parameter values were used to generate simulation data: *T _{U}*(0) = 600,

To evaluate the performance of the MSSB and SNLS approaches for smaller samples sizes and larger variances that are similar to the actual clinical data evaluated in the next section, we also performed the simulations for the number of measurements *n* = *n _{T}* =

$$\mathit{ARE}=\frac{1}{N}\sum _{j=1}^{N}\frac{\mid \theta -{\widehat{\theta}}_{j}\mid}{\mid \theta \mid}\times 100\%,$$

where * _{j}* is the estimate of parameter

The simulation results are reported in Table 1. From these results, we can see that the SNLS approach significantly improves the performance of the MSSB approach in all simulation cases as we expected. When the number of measurements is large (e.g., *n* = *n _{T}* =

We applied the proposed MSSB and SNLS approaches to estimate viral dynamic parameters from data for individual AIDS patients based on model Eq. (1.1) - (1.3). Two HIV-1 infected patients were treated with a four-drug antiretroviral regimen in combination with an immune-based therapy. Frequent viral load measurements were scheduled at baseline and after initiating the combination treatment: 13 measurements during the first day, 14 measurements from day 2 to week 2, and then one measurement at each of the following weeks, 4, 8, 12, 14, 20, 24, 28, 32, 36, 40, 44, 48, 52, 56, 64, 74 and 76, respectively. Total CD4+ T cell counts were monitored at weeks 2, 4 and monthly thereafter.

As suggested in Section 2, first we applied the nonlinear regression model in Perelson et al. (1996) to fit the viral load data for the first week to obtain the estimates of *δ* and *c*. The estimation results for the two patients are: Patient I, *δ* = 1:09 and *c* = 2:46; and Patient II, *δ* = 0:43 and *c* = 3:78. Then we applied the proposed MSSB and SNLS approaches to estimate all other viral dynamic parameters including the time-varying parameter *η*(*t*) in model (1.1) - (1.3). The log_{10} transformation of the data was used in order to stabilize the variance and computational algorithms.

A rough estimate and a reasonable search range for each of the unknown parameters and the initial values of the state variables were obtained using the MSSB method. Then we applied the SNLS approach to refine the estimates. For the SNLS approach, we selected the smoothing parameters using the AIC, BIC, and AICc criteria. We considered all the combinations of two spline orders (3 and 4) and 8 different numbers of control points (from 3 to 10), and the results are reported in Table 2 and Table 3 for Patient I and Patient II, respectively. Note that as a practical guideline, if the number of unknown parameters exceeds *n*/40 (where *n* is the number of measurements), the AICc instead of AIC should be used. For our clinical data, *n* is about 40, and the number of unknown parameters varies between 8 and 15 for different models (the unknown initial conditions were also considered as unknown parameters), which is much larger than *n*/40 = 40/40 = 1. So the AICc is more appropriate for our case. In fact, the AICc converges to AIC as the number of measurements gets larger, thus the AICc is often suggested to be employed regardless of the number of measurements (Burnham and Anderson, 2004). Therefore, our model selection is mainly based on AICc although the AIC and BIC scores are also reported in Table 2 and Table 3 for comparisons. From Table 2, we can see that both AICc and BIC selected the best spline model for *η*(*t*) as order 3 and 5 control points for Patient I although the AIC selected a different model. From Table 3, we can see that the AICc selected the same spline model (order 3 with 5 control points) for Patient II and both AIC and BIC ranked this model as the second best model. Based on above discussions, we are confident that the spline model with order 3 and 5 control points is the best approximation to the time-varying parameter *η*(*t*). Thus, we will mainly report and discuss the results from this model.

Model Selection for Patient I: the time-varying parameter η(t) is approximated using splines of order 3 or 4 with the number of knots from 3 to 10

Model Selection for Patient II: the time-varying parameter η(t) is approximated using splines of order 3 or 4 with the number of knots from 3 to 10

Fig. 1 shows the data and model fitting results from the best SNLS approach for the two patients. The fitted curves are reasonably well. Table 4 reports the estimation results from both MSSB and SNLS methods for the two patients. The MSSB approach provided good search ranges of unknown parameters for the SNLS approach. From the SNLS estimates, the proliferation rates of uninfected CD4+ T cells are λ = 397 and 45 cells per day, respectively for Patient I and II; the death rates are *ρ* = 0.49 and 0.10 per day which correspond to the half-life of 1.4 and 6.9 days for the two patients, respectively; the numbers of virions produced by each of the infected cells are 265 and 1114 per cell for Patient I and II, respectively.

The estimated values and associated 95% confidence intervals (CI) of viral dynamic parameters for two AIDS patients. Perelson’s model (Perelson et al., 1996) was first fitted to the viral load data within the first week to obtain estimates of **...**

Note that the nonlinear regression estimates of the death rate of infected cells, *δ*, for the two patients are 1.09 and 0.43 per day with the corresponding half-life of 0.64 and 1.61 days, respectively, which indicates a higher death rate of infected cells compared to the uninfected cells for both patients. The estimates of viral clearance rate (*c*) were 2.46 and 3.78 per day with the corresponding half-life of 0.28 and 0.18 days for the two patients, respectively. These results are similar to those from the previous studies (Perelson et al. 1996, 1997) and biologically reasonable. However, to our knowledge, this is the first time the estimates of the proliferation rate and death rate of uninfected CD4+ T cells (λ and *ρ*) and the number of virions produced by an infected cell (*N*) have been obtained directly from the clinical data of AIDS patients. From the estimation results in Table 4, we can see that the death rate of infected cells is 2 to 4 fold higher than that of uninfected cells for Patient I and II, respectively. We can also see that the difference in the estimated dynamic parameters between the two patients is large. This is consistent with the argument that the between-subject variation of viral dynamics in AIDS patients is large (Wu and Ding, 1999; Wu et al., 1999). Thus, the personalized treatment is necessary for AIDS patients.

The estimated time-varying parameter, the infection rate *η*(*t*) for Patient I and II are plotted in Fig. 2. From the estimates, we can see that the infection rates for the two patients have similar patterns which are primarily due to similar patterns of viral load and CD4+ T cell counts for the two patients. It seems that there was an initial fluctuation in the infection rate at the beginning of treatment, and then the infection rate was stabilized after one or two weeks. We can also see that the estimated infection rate is not zero, which may suggest the imperfectness of the long-term treatment and the development of drug-resistant mutants. This is inconsistent with the perfect treatment assumption in Perelson et al. (1996, 1997) although it may be valid in a short period of time after initiating the treatment.

We have proposed two approaches to identify all dynamic parameters including both constant and time-varying parameters in an HIV viral dynamic model which is characterized by a set of nonlinear differential equations. This is a very challenging problem in the history of HIV dynamic studies. The proposed multi-stage smoothing-based (MSSB) approach is straightforward and easy to implement since it does not require numerically solving the ODEs and the initial values of the state variables are not needed. However, one limitation of this approach is that it requires the estimates of derivatives of the state variables, which are usually poor when the data are sparse. This may result in poor estimates of unknown dynamic parameters. The second method that we proposed is the spline-enhanced nonlinear least squares (SNLS) approach. This method is more accurate in estimating the unknown parameters, but it requires numerical evaluations of ODEs repeatedly in the optimization procedure and the convergence of the computational algorithm is problematic when the parameter space is high-dimensional. We have proposed a hybrid optimization technique to deal with the non-convergence and the local solution problems, but the computational cost is high and a good search range for each of the unknown parameters is required. In practical implementation, we propose to combine the two approaches, i.e., first using the MSSB approach to obtain rough estimates and possible ranges for the unknown parameters and then using the SNLS approach to refine the estimates. We have applied this strategy to estimate all dynamic parameters from data for two AIDS patients.

It is very important to estimate all dynamic parameters in HIV dynamic models directly from clinical data when prior knowledge is not available or fixing some parameters will significantly affect the estimates of other parameters in a dynamic model. To our knowledge, this is the first attempt to propose a procedure to estimate both constant and time-varying parameters in the proposed HIV dynamic model directly from clinical data. Although Huang, Liu and Wu (2006) have attempted to do so using a Bayesian approach, strong priors were used for most of the dynamic parameters and a parametric form for the time-varying parameter was employed. Notice that the time-varying parameter *η*(*t*) is a function of antiviral treatment effect in the HIV dynamic model (Huang et al., 2006). It is very important to estimate this time-varying parameter for each AIDS patient individually in order to better design the treatment strategy and personalize the treatment for each individual. In addition, to better assess the treatment efficacy, the time-varying infection rate *η*(*t*) can be further modeled as *η*(*t*) = *η _{pre}* × (1 −

We used both viral load and CD4+ T cell count data in order to identify all dynamic parameters in the HIV dynamic model (1.1)-(1.3). Usually AIDS investigators believe that CD4+ T cell count data are very noisy with a large variation due to both measurement errors and natural variations of CD4+ T cell counts. Thus, it is important to improve the data quality in order to get more reliable estimates of dynamic parameters although it is beyond statisticians’ control. Our study also suggests that the frequent measurements, in particular the viral load measurements during the early stage after initiating an antiviral therapy, are important to identify some constant dynamic parameters such as *c* and *δ*, and long-term monitoring of viral load and CD4+ T cell counts is necessary to estimate other parameters, in particular, the time-varying parameter.

This research was partially supported by the NIAID/NIH grants AI62247, AI059773, AI50020, and AI055290. Also, Liang’s research was partially supported by the NSF grant DMS0806097, and Miao’s research was partially supported by the Provost’s Multidisciplinary Award and the DCFAR Mentoring Award from University of Rochester. The authors thank an Associate Editor and referees for their constructive comments.

The differential evolution (DE) algorithm (Storn and Price, 1997) is a typical evolutionary algorithm that searches the optimum by inheritance, mutation, selection, and crossover of parent populations (e.g., vectors of possible parameter values in the parameter space). The initial population is randomly generated from a uniform distribution within the search range to cover the whole region, denoted by *x _{i},G*(

$${\nu}_{i,G+1}={x}_{{r}_{1},G}+F({x}_{{r}_{2},G}-{x}_{{r}_{3},G})$$

(A.1)

where *v _{i,G}*

The scatter search method was first proposed by Glover (1977) and extended later by Laguna and Marti (2003, 2005). The scatter search method is not a genetic algorithm; instead, it locates the optimum by tracking the search history and balancing visiting frequency within different segments of the search region using sophisticated strategies. Here we give a brief introduction to this method. For more details, the reader is referred to Laguna and Marti (2003) and Rodriguez-Fernandez, Egea and Banga (2006).

Let *l _{i}* and

$$\mathbf{F}=\left[\begin{array}{cccc}{f}_{\mathrm{11}}& {f}_{12}& \cdots & {f}_{1m}\\ {f}_{21}& {f}_{22}& \cdots & {f}_{2m}\\ \vdots & \vdots & \vdots & \vdots \\ {f}_{q1}& {f}_{q2}& \cdots & {f}_{qm}\end{array}\right].$$

(A.2)

Then we can calculate the probability of historic visits to *s _{ij}* as

$${p}_{ij}=\frac{\frac{1}{{f}_{ij}}}{\sum _{k=1}^{m}\frac{1}{{f}_{ik}}},i=1,2,\cdots ,q;j=1,2,\dots ,m.$$

(A.3)

Since the scatter search method is population-based, the first population needs to be generated to start the algorithm. Let *N _{G}* denotes the total number of parameter vectors in the first population, then

$${z}_{i}\le \sum _{j=1}^{k}{p}_{ij},k=1,2,\cdots ,m.$$

Note that after each new parameter vector is generated, the matrix **F** is also updated.

Once the first population is generated, the parameter vectors are divided into two categories: *elite vectors* and *diverse vectors*. Let *n _{e}* denote the number of all elite vectors, which is usually a small fixed number (e.g., 20). Then

- Type 1:
*p*_{1}=*x*^{(1)}*− d* - Type 2:
*p*_{2}=*x*^{(1)}+*d* - Type 3:
*p*_{3}=*x*^{(2)}+*d*

where *d* = *r*^{T} · (*x*^{(2)} − *x*^{(1)}) with all the components of *r* generated from an uniform distribution on [0, 1]. If both *x*^{(1)} and *x*^{(2)} are in the first half elite vectors in terms of fitness, one new vector of type 1, one of type 3 and two of type 2 are generated; if only *x*^{(1)} belongs to the first half elite vectors in terms of fitness, one new vector of each type is generated; if both *x*^{(1)} and *x*^{(2)} belong to the second half elite vectors in terms of farthest distance, then one new vector of type 2 and one of either type 1 or 3 are generated. The fitness of these new generated vectors is then compared to that of the elite vectors. The new vectors with smaller objective function values (that is, better fit) will replace the elite vectors with larger objective function values. This procedure is repeated until no new vectors can replace any elite vectors. Now the algorithm can either stop or continue by regenerating diverse vectors.

^{*}Supported in part by NIAID/NIH grants AI055290, AI50020, RR024160, and AI078498.

1. Akaike H. Information theory as an extension of the maximum likelihood principle. In: Petrov BN, Csaki F, editors. Second International Symposiumon Information Theory. Budapest: Akademiai Kiado; 1973. pp. 267–281.

2. 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(1):55–65. [PubMed]

3. Bard Y. Nonlinear Parameter Estimation. London: Academic; 1974.

4. Bellman R, Åström KJ. On structural identifiability. Mathematical Biosciences. 1970;7:329–339.

5. Burnham KP, Anderson DR. Multimodel inference: understanding AIC and BIC in model selection. Sociological Methods Research. 2004;33:261.

6. Carra’Ferro G, Gerdt VP. Improved Kolchin-Ritt algorithm. Programming and Computer Software. 2003;29(2):83–87.

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

8. Chen J, Wu H. Efficient local estimation for time-varying coefficients in deterministic dynamic models with applications to HIV-1 dynamics. Journal of the American Statistical Association. 2008;103:369–384.

9. Chen J, Wu H. Estimation of time-varying parameters in deterministic dynamic models with application to HIV infections. Statistica Sinica. 2009 in press.

10. Cobelli C, Lepschy A, Jacur R. Identifiability of compartmental systems and related structural properties. Mathematical Biosciences. 1979;44:1–18.

11. De Boor C. A Practical Guide to Splines. New York: Springer-Verlag; 1978.

12. Englezos P, Kalogerakis N. Applied Parameter Estimation for Chemical Engineers. New York: Marcel-Dekker; 2001.

13. Fan J, Gijbels I. Local Polynomial Modeling and its Applications. London: Chapman and Hall; 1996.

14. Filter RA, Xia X, Gray CM. Dynamic HIV/AIDS parameter estimation with appliation to a vaccine readiness study in Southern Africa. IEEE Transactions on Biomedical Engineering. 2005;52(5):284–291. [PubMed]

15. Glover F. Heuristics for integer programming using surrogate constraints. Decision Sciences. 1977;8:156–166.

16. Gray CM, Williamson C, Bredell H, Puren A, Xia X, Filter R, Zijenah L, Cao H, Morris L, Vardas E, Colvin M, Gray G, McIntyre J, Musonda R, Allen S, Katzenstein D, Mbizo M, Kumwenda N, Taha T, Karim SA, Flores J, Sheppard HW. Viral dynamics and CD4+ T celll counts in subtype C human immunodeficiency virus type 1-infected individuals from Southern Africa. Aids Research and Human Retroviruses. 2005;21(4):285–291. [PubMed]

17. Hemker PW. Numerical methods for differential equations in system simulation and in parameter estimation. In: Hemker HC, Hess B, editors. Analysis and Simulation of Biochemical Systems. 1972. pp. 59–80.

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

19. Huang Y, Wu H. A Bayesian approach for estimating antiviral efficacy in HIV dynamic models. Journal of Applied Statistics. 2006;33:155–174.

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

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

22. Joshi M, Seidel-Morgenstern A, Kremling A. Exploiting the bootstrap method for quantifying parameter confidence intervals in dynamic systems. Metabolic Engineering. 2006;8:447–455. [PubMed]

23. Kolchin E. Differential Algebra and Algebraic Groups. New York: Academic Press; 1973.

24. Laguna M, Marti R. Scatter Search: Methodology and Implementations in C. Boston: Kluwer Academic Publishers; 2003.

25. Laguna M, Marti R. Experimental testing of advanced scatter search designs for global optimization of multimodal functions. Journal of Global Optimization. 2005;33:235–255.

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

27. Li Z, Osborne M, Prvan T. Parameter estimation in ordinary differential equations. IMA Journal of Numerical Analysis. 2005;25:264–285.

28. Liang H, Wu H. Parameter estimation for differential equation models using a framework of measurement error in regression models. Journal of the American Statistical Association. 2008;103:1570–1583. [PMC free article] [PubMed]

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

30. Miao H, Dykes C, Demeter LM, Cavenaugh J, Park SY, Perelson AS, Wu H. Modeling and estimation of kinetic parameters and replicative fitness of HIV-1 from flow-cytometry-based growth competition experiments. Bulletin of Mathematical Biology. 2008;70:1749–1771. [PubMed]

31. 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(1):292–300. [PMC free article] [PubMed]

32. Moles CG, Banga JR, Keller K. Solving nonconvex climate control problems: pitfalls and algorithm performances. Appl Soft Comput. 2004;5:35–44.

33. Nocedal J, Wright SJ. Numerical Optimization. New York: Springer Verlag; 1999.

34. Nowak MA, May RM. Virus Dynamics: Mathematical Principles of Immunology and Virology. Oxford: Oxford University Press; 2000.

35. Ogunnaike BA, Ray WH. Process Dynamics, Modeling, And Control. New York: Oxford University Press; 1994.

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

37. Ouattara DA, Mhawej MJ, Moog CH. Clinical tests of therapeutical failures based on mathematical modeling of the HIV infection. Joint special issue of IEEE transactions on circuits and systems and IEEE transactions on Automatic Control (Special issue on systems biology) 2008;53:230–241.

38. Perelson AS, Neumann AU, Markowitz M, Leonard JM, Ho DD. HIV-1 dynamics in vivo: virion clearance rate, infected cell life-span, and viral generation time. Science. 1996;271:1582–1586. [PubMed]

39. Perelson AS, Essunger P, Cao YZ, Vesanen M, Hurley A, Saksela K, Markowitz M, Ho DD. Decay characteristics of HIV-1-infected compartments during combination therapy. Nature. 1997;387:188–191. [PubMed]

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

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

42. Poyton AA, Varziri MS, Mcauley KB, Mclellan PJ, Ramsay JO. Parameter estimation in continuous-time dynamic models using principal differential analysis. Computer and Chemical Engineering. 2006;30:698–708.

43. Putter H, Heisterkamp SH, Lange JMA, Wolf F. A Bayesian approach to parameter estimation in HIV dynamic models. Statistics in Medicine. 2002;21:2199–2214. [PubMed]

44. Ramsay JO. Principal differential analysis: data reduction by differential operators. Journal of the Royal Statistical Society, Series B. 1996;58:495–508.

45. Ramsay JO, Silverman BW. Functional Data Analysis. 2. New York: Springer; 2005.

46. Ramsay JO, Hooker G, Campbell D, Cao J. Parameter estimation for differential equations: a generalized smoothing approach (with discussions) Journal of the Royal Statistical Society, Series B. 2007;69:741–796.

47. Ritt JF. Differential Algebra. Providence, RI: Amer Math Soc; 1950.

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

49. Runge C. Über empirische Funktionen und die Interpolation zwischen äquidistanten Ordinaten. Zeitschrift für Mathematik und Physik. 1901;46:224–243.

50. Schwarz G. Estimating the dimensions of a model. The Annals of Statistics. 1978;6:461–464.

51. Storn R, Price K. Differential evolution - a simple and efficient heuristic for global optimization over continuous spaces. Journal of Global Optimization. 1997;11:341–359.

52. Tan WY, Wu H. Deterministic and Stochastic Models of AIDS Epidemics and HIV Infections with Intervention. Singapore: World Scientific; 2005.

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

54. Varah JM. A spline least squares method for numerical parameter estimation in differential equations. SIAM Journal on Scientific Computing. 1982;3:131–141.

55. Vrugt JA, Robinson BA. Improved evolutionary optimization from genetically adaptive multimethod search. PNAS. 2007;104:708–711. [PubMed]

56. Walter E. Identifiability of Parameteric Models. Oxford: Pergamon Press; 1987.

57. Wei X, Ghosh SK, Taylor ME, Johnson VA, Emini EA, Deutsch P, Lifson JD, Bonhoeffer S, Nowak MA, Hahn BH, Saag MS, Shaw GM. Viral dynamics in human immunodeficiency virus type 1 infection. Nature. 1995;373:117–122. [PubMed]

58. Wu H, Ding A. Population HIV-1 dynamics in vivo: applicable models and inferential tools for virological data from AIDS clinical trials. Biometrics. 1999;55:410–418. [PubMed]

59. Wu H, Kuritzkes DR, Mcclernon DR, Kessler H, Connick E, Landay A, Spear G, Heath-Chiozzi M, Rousseau F, Fox L, Spritzler J, Leonard JM, Lederman MM. Characterization of viral dynamics in human immunodeficiency virus type 1-infected patients treated with combination antiretroviral therapy: relationships to host factors, cellular restoration and virological endpoints. Journal of Infectious Diseases. 1999;179:799–807. [PubMed]

60. Wu H. Statistical methods for HIV dynamic studies in AIDS clinical trials. Statistical Methods in Medical Research. 2005;14:171–192. [PubMed]

61. Wu H, Zhu H, Miao H, Perelson AS. Parameter identifiability and estimation of HIV/AIDS dynamic models. Bulletin of Mathematical Biology. 2008;70:785–799. [PubMed]

62. Wu H, Zhang JT. Nonparametric Regression Methods for Longitudinal Data Analysis. New Jersey: Wiley; 2006.

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

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

65. Ye Y. PhD thesis. Department of ESS, Stanford University; 1987. Interior Algorithms for Linear, Quadratic and Linearly Constrained Non-linear Programming.

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