Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
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-AOAS290
PMCID: PMC2885820



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.

Keywords and phrases: HIV viral dynamics, Ordinary differential equation, time-varying parameter, parameter identifiability, differential algebra, hybrid optimization, local polynomial smoothing, semiparametric regression

1. Introduction

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




where TU is the concentration of uninfected target CD4+ T cells, TI the concentration of infected cells, V the viral load, λ the proliferation rate of uninfected target cells, ρ the death rate of uninfected target cells, η(t) the time-varying infection rate depending on antiviral drug efficacy, δ the death rate of infected cells, c the clearance rate of free virions, and N the number of virions produced by a single infected cell on average. In this system, TU(t), TI(t) and V(t) are state variables, and (λ, ρ, N,δ, c,η(t))T are unknown dynamic parameters. Notice that here we do not distinguish the infected cells by various sub-populations such as productively infected cells, latently infected cells and long-lived infected cells (Perelson et al. 1996, 1997), since we intend to model the viral dynamics for an HIV-infected patient under long-term treatment. We also use a time-varying parameter η(t) to model the infection rate since the infection rate may change nonparametrically due to the variation in treatment effect over time. Generally speaking, all kinetic parameters in this model could be time-varying (not only η(t)); however, it is usually unfeasible to do so due to the limited data collected in the clinical studies. Thus, this model provides a flexible although simple approach for studying long-term viral dynamics. In clinical trials or clinical practice, viral load, V(t), and total CD4+ T cell count, T(t) = TU(t)+TI(t), are closely monitored and measured over time.

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.

2. A Multi-Stage Smoothing-Based (MSSB) Approach

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.

2.1. Stage I: Local polynomial estimation of the state variables

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, (Y1, …, Yn), i.e., Yi = X(ti)+ ei for i = 1, …, n, where (e1, e2, …, en) are independent measurement errors with mean zero. Assume that the fourth derivative of X(t) exists. For each given time point t0, we approximate the function X(ti) locally by a pth-order polynomial,, X(ti)X(t0)+(tit0)X(1)(t0)++X(p)(t0)(tit0)p/p!j=0pζj(t0)(tit0)j, for ti, i = 1, …, n; in a neighborhood of the point t0, where ζj(t0) = X(j)(t0) for j = 0, 1, …, p. Following the local polynomial fitting (Fan and Gijbels 1996), the estimators X(ν) (t) of X(ν)(t) (ν = 0, 1, 2 in our case) can be obtained by minimizing the locally weighted least-squares criterion,


where K(·) is a symmetric kernel function, Kh(·) = K(·/h)/h, and h is a proper bandwidth. Let Y = (Y1, …, Yn)T, denote Tp,t as an n ×(p + 1) design matrix whose ith row is (1, tit, …,(tit)p)T, and Wt an n × n diagonal matrix of kernel weights; i.e., diag{Kh(t1t), …, Kh(tnt)}. Assuming that the matrix Tp,tTWtTp,t is not singular, a standard weighted least squares theory leads to the solution, ζ^=(Tp,tTWtTp,t)1Tp,tTWtY, where p = 1, 2, 3. We use the local linear regression to estimate X(t), the local quadratic regression to estimate X′ (t), and the local cubic regression to estimate X(2)(t). Consequently, using the above notations, the estimators X(q)(t) can be expressed as


where ξq is the (q + 2) × 1 vector having 1 in the (q + 1) entry and zeros in the other entries.

We apply the local polynomial smoothing technique to the viral load data, V(t) and total CD4+ T cell count data, T(t) = TU(t) + TI(t) to obtain the estimates of V(t), T(t), V′ (t) and T′ (t) that will be used in Stage II, and the estimate of V(2)(t) that will also be used in Stage III.

2.2. Stage II: PsLS estimates

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


Recalling that T(t) = TI(t) + TU(t) and substituting TU(t) = T(t) − TI(t) in the above equation, we obtain


From the above equation, we obtain


where T′ = dT(t)/dt. Substituting this expression into equation (1.3) and letting α0=Nδλρδ,α1=Nδρρδandα2=Nδρδ,we have


where V(t) and T(t) for t = t1, t2, …, tn are the measurements of viral load and CD4+ T cell count from clinical studies which are subject to measurement errors.

Let V(t), T(t), V′(t) and T′(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),


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, T,T′, −V). Let μl(K)=11zlK(z)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 [alpha] 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 h2, 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,

2.3. Stage III: Semiparametric regression for estimation of both constant and time-varying parameters

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


while from (1.3), we can obtain


A combination of (2.3) and (2.4) deduces


Let Z(t) = V(2)(t)+ĉV′(t), U1(t) = −{V′(t)V(t)+ĉV2(t)}, U2(t) = − {V′(t)+ ĉV(t)}, U3(t) = T(t)V(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):


where ε*(t) includes all substitution errors. Note that δ and 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.,


where aj are constant B-spline coefficients, and bj,k(t) the basis functions of order k. For more details about B-splines such as construction of higher order basis functions via recurrence relations, the reader is referred to de Boor (1978). Thus, model (2.6) can be approximately written as


This becomes a standard linear regression model and we can easily estimate the unknown constant parameters δ, aj, and Nδaj(j = 1, 2, …, s) from which we can derive the estimates of δ, N, and η(t). Note that AIC, BIC or AICc can be used to determine the order of the B-splines in our numerical data analysis. See detailed discussion on this in the next section.

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.

3. The Spline-Enhanced Nonlinear Least Squares (SNLS) Approach

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.

3.1. Spline Approximation and Nonlinear Least Squares

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., η(t)j=1sajbj,k(t). 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




where θ = (λ, ρ, N, δ, c, a1, a2, …, as)T are unknown constant parameters. Note that we have measurements of total CD4+ T cell counts, T(ti) = TU(ti) + TI(ti), and viral load, V(ti), i.e., the measurement model can be written as



where ε1i and ε2j are assumed to be mean zero with constant variances and independent. Then the standard NLS estimator can be obtained by minimizing the objective function


where nT is the total number of CD4+ T cell measurements and nV is the total number of viral load observations; and T(ti) and V(tj) are numerical solutions to Eq. (3.1) - (3.3) using the Runge-Kutta method. However, if ε1i and ε2j are correlated, the weighted (generalized) NLS method can be used. To stabilize the variance and the computational algorithms, the log-transformation of the data is usually used in practice. More generally, the weighted NLS approach can be used to more efficiently estimate the unknown parameters if the weights for different terms in the objective function are known. Usually the Fisher-Information-Matrix can be used to obtain the confidence intervals for the unknown parameters, but the bootstrap approach is more precise although it is more computationally intensive (Joshi, Seidel-Morgenstern, Kremling, 2006). We will use the bootstrap method in our real data analysis.

3.2. Hybrid Optimization and Spline Parameter Selection

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 4th 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.,




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




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.

4. Simulation Studies

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: TU(0) = 600, TI(0) = 30, V(0) = 105, λ = 36, ρ = 0.108, N = 1000, δ = 0.5, c = 3, and η(t) = 9 × 10−5× {1 0.9 cos(πt/ 1000)}. Let T = TU + TI denote the total number of infected and uninfected CD4+ T cells and V denote the viral load, the measurement model (3.4)-(3.5) was used to simulate the observation data with measurement noise, where ε1i and ε2j were assumed to be independent and generated from normal distributions with mean zero and constant variances σ12 and σ22, respectively. Eq. (1.1)-(1.3) were numerically solved within the time range [0, 20] (days) using the Runge-Kutta method, and solutions were output for equally-spaced time intervals of 0.1 nd 0.2 which correspond to the number of measurements 200 and 100, respectively. The measurement errors σ12 = 20, 30, 40 and σ22 = 100, 150, 200 were added to the numerical results of the ODE model according to the measurement model (3.4) - (3.5), respectively.

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 = nT = nV = 30 and 50, and σ12 = 202, 302, 402 and σ22 = 502, 752, 1002. To evaluate the performance of the estimation methods, we use the average relative estimation error (ARE) which is defined as


where [theta w/ hat]j is the estimate of parameter θ from the jth simulation data set and N is the total number of simulation runs. We applied the MSSB and SNLS approaches to 500 simulated data sets for each simulation scenario. We used the AICc and determined that the time-varying parameter η(t) was approximated by a spline of order 2 with 3 knots in our simulation studies.

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 = nT = nV = 100 or 200) and variances are small, the SNLS approach performs extremely well, but the AREs of the MSSB estimates of some parameters such as ρ and δ are very large. When the number of measurements is small (similar to our real data analysis in the next section, i.e., n = nT = nV = 30 or 50) and the variances are large, the performance of the SNLS estimates is still reasonably good, but the MSSB approach performs poorly for all parameter estimates. Similarly, the estimate of time-varying parameter η(t) from the SNLS approach outperforms that of the MSSB approach. Thus, the simulation results and our experience suggest that it is a good strategy to use the MSSB approach to obtain rough search ranges for unknown parameters, and then use the SNLS approach to refine the estimates.

Table 1
Simulation study: a comparison of the MSSB and SNLS approaches. The average relative error (ARE) was calculated based on 500 simulation runs. The time-varying parameter η(t) is approximated by a 3-knots spline of order 2

5. Viral Dynamic Parameter Estimation for AIDS Patients

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

Table 2
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
Table 3
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.

Fig 1
The SNLS model fitting for two AIDS patients
Table 4
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.

Fig 2
The estimated infection rate η(t) (solid line) and its bootstrap 95% confidence interval (dashed line)

6. Discussion and Conclusion

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 −ηtreat(t)), where ηpre denotes the infection rate at the pre-treatment equilibrium, and ηtreat(t) the time-varying infection rate after treatment. By comparing ηpre with η(t), it will be easy to quantify the treatment effects on infection rate. However, since ηpre is strongly correlated with ηtreat(t), ηpre has to be fixed to estimate ηtreat(t). Unfortunately, the pre-treatment equilibrium data were not collected in this study such that we could not determine the value of ηpre and therefore to do the comparison. Finally, the treatment could also affect parameter N; however, limited by the clinical measurements, we simplified our model without considering parameter N as time-varying.

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 xi,G(i = 1, 2, …, NP), where NP is the number of members in this generation. The subsequent population inherits and mutates by randomly mixing the previous generation with certain weights. Storn and Price (1997) recommended a mutation method


where vi,G+1 is the ith member of generation (G + 1), xr1,G, xr2,G and xr3,G the members of the previous generation G, and F > 0 the amplification factor. Integers r1, r2 and r3 are randomly chosen from {1, 2, …, NP } which are mutually different from each other and different from i. Also, to increase diversity, the mutated member vi,G+1 exchanges its components with xi,G with a given probability, the crossover ratio. For this purpose, a number within [0, 1] is generated for each component of vi,G+1 by a uniform random number generator. If this number is greater than the crossover ratio, then the component of xi,G is kept. Finally, the best member in the mutated population is selected by comparing the values of objective function. The convergence rate of the differential evolution method depends on specific mutation strategies used, i.e., the amplification factor and the crossover ratio. This method has been successfully employed in previous studies (Miao et al., 2008; Miao et al., 2009). For more details about the DE algorithm, the reader is referred to Storn and Price (1997).

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 li and ui denote the lower and upper boundary of the ith component of the parameter vector, respectively. Then the region between li and ui can be divided into m segments (e.g., m = 4), which are usually of equal length. Let sij (j = 1, 2, …, m) denote the jth segment of the search region [li, ui ]. If a parameter value candidate of θi falls in sij, we call it a visit to sij. Let fij denote the total number of visits to sij, then a visiting history matrix F can be formed


Then we can calculate the probability of historic visits to sij as


Since the scatter search method is population-based, the first population needs to be generated to start the algorithm. Let NG denotes the total number of parameter vectors in the first population, then m parameter vectors are generated first such that the ith component of the jth vector falls into sij. In this way, all elements of matrix F become one. To generate the rest (NGm) parameter vectors in the first population, a random vector z = (z1, z2, …, zq)T is generated for each parameter vector with zi(i = 1, 2, …, q) following a uniform distribution on [0, 1]. Then the ith component of the parameter vector to be generated falls into sik for the smallest k satisfying


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 ne denote the number of all elite vectors, which is usually a small fixed number (e.g., 20). Then ne/2 elite vectors are chosen based on the smallest objective function values; the rest half of the elite vectors are chosen to be farthest from all the first half elite vectors in the sense of Euclidean distances. Thus, both fitness and diversity are considered in the construction of the elite vectors. Then the new parameter vectors can be generated by combination of the elite vectors. Let x(1) and x(2) denote two different elite vectors and x(1) has a smaller objective function value than x(2). Then three types of combinations can be employed to generate new parameter vectors

  • Type 1: p1 = x(1)− d
  • Type 2: p2 = x(1) + d
  • Type 3: p3 = x(2) + d

where d = rT · (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.