Hybrid Taguchi-Differential Evolution Algorithm for Parameter Estimation of Differential Equation Models with Application to HIV Dynamics

Abstract

This work emphasizes solving the problem of parameter estimation for a human immunodeficiency virus (HIV) dynamical model by using an improved differential evolution, which is called the hybrid Taguchi-differential evolution (HTDE). The HTDE, used to estimate parameters of an HIV dynamical model, can provide robust optimal solutions. In this work, the HTDE approach is effectively applied to solve the problem of parameter estimation for an HIV dynamical model and is also compared with the traditional differential evolution (DE) approach and the numerical methods presented in the literature. An illustrative example shows that the proposed HTDE gives an effective and robust way for obtaining optimal solution, and can get better results than the traditional DE approach and the numerical methods presented in the literature for an HIV dynamical model.

1. Introduction

Recently, the mathematical modeling of the epidemiology and immunology dynamics of human immunodeficiency virus (HIV) has been proven to be valuable in understanding the HIV pathogenesis (see, e.g., [1–15] and the references therein). Many of the proposed approaches only work under the assumption that the parameters of the HIV dynamical models are known in advance. However, in the real world, it may be difficult to determine the parameters because of the complexity of HIV dynamical models. Therefore, in recent years, it has become a hot issue for estimating the parameters of the HIV dynamical models (see, e.g., [16–25] and the references therein). On the other hand, for solving the problem of parameter estimation in the HIV dynamical models, the methods of estimation such as the Bayesian approach [17, 21], the online estimation algorithm [18], the Monte Carlo technique [20], the numerical method with combined Adomian/Alienor approach [22], the multiple time point method [23], the nonlinear least squares method [24], and the multistage smoothing-based method and spline-enhanced nonlinear least squares approach [25] have been presented in the literature. In particular, Miao et al. [24] introduced the differential evolution (DE), a global optimization algorithm, to solve the nonlinear least squares problem to avoid the local minima and convergence difficulty in nonlinear optimization.

The DE [26] has received considerable attention concerning its potential as a new optimization technique for complex nonlinear problems and has been successfully used in various areas (see, e.g., [27–34] and the references therein). The main specific feature of the DE as an optimization method is its implicit parallelism, which is a result of the evolution and the hereditary-like process. The advantages of using a DE for solving the optimization problems are its global solution finding property, simple but powerful search capability, easy-to-understand concept, compact structure, having only a few control parameters, ease of use, and high convergence characteristics [35, 36].

The DEs are designed by using the concept of evolutionary algorithms, such as multipoint searching, mutation operation crossover operation and selection operation, but crossover operation based on the random process is not a systematic reasoning way for breeding better offspring (or trial individual vectors). Therefore, in order to seek the optimal breeding in the DEs such that the efficiency of the DEs for estimating all the parameters in the HIV dynamical model [37] can be further promoted, the purpose of this paper is to propose a new DE approach by introducing a systematic reasoning. This proposed new DE approach is named the hybrid Taguchi-differential evolution (HTDE). The HTDE combines the DE [26] with the Taguchi method [38, 39]. The Taguchi method, a robust design approach, uses many ideas from statistical experimental design, where some of the factors (or individuals) are related, for evaluating and implementing improvements in products, processes, and equipment. Factors (or individuals) are called related when the desirable experimental region of some factors (or individuals) depends on the level settings of other factors (or individuals). Two major tools used in the Taguchi method are (i) signal-to-noise ratio (SNR) which measures quality, and (ii) orthogonal arrays which are used to study many design parameters simultaneously [40, 41]. In the HTDE, the Taguchi method is to provide a new systematic crossover operation to replace the original crossover operation of DE. Then, the systematic reasoning ability of the Taguchi-method-based crossover operation is used to breed better individuals in order to generate the representative individuals to be the new potential offspring (or trial vectors). So, the Taguchi-method-based crossover operation can enhance DEs, such that the HTDE can be robust, statistically sound, and quickly convergent. Thus, the HTDE can be an efficient optimization algorithm to deal with the problem of parameter estimation for an HIV dynamical model, which is also regarded an optimization problem.

This paper is organized as follows. Section 2 describes the HTDE for solving the parameter estimation problem of an HIV dynamical model. In Section 3, we evaluate the efficiency of the proposed HTDE by comparing our results with those obtained from the traditional DE approach and the numerical methods of Manseur et al. [22] respectively, by using the same example given by Manseur et al. [22]. Finally, Section 4 offers some conclusions.

2. Parameter Estimation of an HIV Dynamical Model with HTDE Approach

In this paper, the following three-dimensional model of HIV [22, 37] described by 𝑥̇₁(𝑡)=𝑆−𝑑𝑥₁(𝑡)−𝛽𝑥₁(𝑡)𝑥₃(𝑡), (2.1a)𝑥̇₂(𝑡)=𝛽𝑥₁(𝑡)𝑥₃(𝑡)−𝜇₁𝑥₂(𝑡), (2.1b)𝑥̇₃(𝑡)=𝑘𝑥₂(𝑡)−𝜇₂𝑥₃(𝑡), (2.1c) is considered, where 𝑥₁(𝑡) denotes amount (quantity) of healthy CD4+ T cells, 𝑥₂(𝑡) denotes the infected CD4+ T cells, 𝑥₃(𝑡) denotes the viral load, and the positive constants 𝑆,𝑑,𝛽,𝜇₁,𝑘, and 𝜇₂ denote the system parameters for an HIV dynamical model described as the rate of production of the healthy cells, the death rate of the health calls, the infection rate of healthy cells CD4 by virus HIV, the death rate of infected cells, the rate of production of free virus, and the death rate of free virus, respectively. Equation (2.1a) denotes the population dynamics of the healthy cells. In the presence of HIV, the healthy cells interact with the virus and its reproduction rate decreases according to the term −𝛽𝑥₁(𝑡)𝑥₃(𝑡). Equation (2.1b) denotes the population dynamics of the infected cells. The growth of infected cells is proportional to the amount of healthy cells infected by virus and will be discounted by the amount of cells in destruction −𝜇₁𝑥₂(𝑡). Equation (2.1c) represents the dynamics of the concentration of free virus. The free virus increases in proportion to the infected cells 𝑘𝑥₂(𝑡), and depend on this natural decline rate −𝜇₂𝑥₃(𝑡).

Remark 2.1. In clinical practice, it is difficult to measure the amount of healthy CD4+ T cells. The measurement of total CD4+ T cells in blood may not be representative of the total body pool of CD4+ T cells that are involved in HIV infection. In addition, the variation of total CD4+ T cells counts in blood is large [23]. Therefore, in this paper, for evaluating the performance of the proposed HTDE approach to an HIV dynamical model given by Manseur et al. [22], and comparing the results obtained from our proposed HTDE approach with those from the traditional DE approach [26] and the numerical methods given by Manseur et al. [22]. By using the same example, we assume the measurements of viral load and T cells (either the amount of healthy CD4+ T cells or the total number of CD4+ T cells) are available.
When estimating the system parameters, suppose the structures of the system is known in advance, and thus the HIV dynamical model can be described as follows: 𝑥̇₁(𝑡)=𝑆−𝑑𝑥₁(𝑡)−𝛽𝑥₁(𝑡)𝑥₃(𝑡), (2.2a)𝑥̇₂(𝑡)=𝛽𝑥₁(𝑡)𝑥₃(𝑡)−𝜇₁𝑥₂(𝑡), (2.2b)𝑥̇₃(𝑡)=𝑘𝑥₂(𝑡)−𝜇₂𝑥₃(𝑡), (2.2c) which is expected to match the actual HIV dynamical model (2.1a)–(2.1c), where the positive constants 𝑆,𝑑,𝛽,𝜇₁,𝑘, and 𝜇₂ are the estimated parameters of an HIV dynamical model, and 𝑥_𝑖(𝑡)(𝑖=1,2,3) denote the state variables of an HIV dynamical model obtained from the estimated parameters {𝑆,𝑑,𝛽,𝜇₁,𝑘,𝜇₂}.
The parameter estimation problem of an HIV dynamical model considered here is to search for the optimal parameters 𝑆,𝑑,𝛽,𝜇₁,𝑘, and 𝜇₂ such that the performance index 𝐽=³_𝑖=1^𝑇_𝑛=0(𝑥_𝑖(𝑛)−𝑥_𝑖(𝑛))² (2.3) is minimized, where 𝑛 is referred to as the sampling time point, 𝑇 is the total number of sampling, and after a period of transient process, a set of state variables is selected as the initial state vector [𝑥₁(0),𝑥₂(0),𝑥₃(0)] for the parameter estimation. This is equivalent to the optimization problem: minimize𝐽=𝑓(𝑆,𝑑,𝛽,𝜇₁,𝑘,𝜇₂)≡𝑓(𝑎₁,𝑎₂,…,𝑎₆) (2.4) subject to 𝑥_𝑖≤𝑎_𝑖≤𝑥_𝑖 (𝑖=1,2,…,6), and 𝑎_𝑖 are the positive real numbers, where 𝑥_𝑖 and 𝑥_𝑖 are the lower and upper bounds of 𝑎_𝑖 given from the practical consideration, respectively. That is, the estimated parameters {𝑎₁,𝑎₂,…,𝑎₆} are generated at random from within user-defined bounds [𝑥_𝑗,𝑥_𝑗]. Due to there often being multiple local optima in the landscape of the performance index 𝐽, it is difficult to obtain the global optimal parameters. Therefore, the proposed HTDE approach will be applied to search for the global or close-to-global optimal parameters 𝑎_𝑖(𝑖=1,2,…,6) such that the performance index 𝐽 is minimized. Therefore, the parameter estimation problem of an HIV dynamical model can become an optimization problem (2.4), where (2.4) is a nonlinear function with the continuous variables.

Remark 2.2. It is obvious that if one set of estimated parameters {𝑎₁,𝑎₂,…,𝑎₆} is specified, then 𝑥_𝑖(𝑡)(𝑖=1,2,3) can be determined from the HIV dynamical model (2.2a)–(2.2c), and thus the value of the performance index (2.3) corresponding to this set of {𝑎₁,𝑎₂,…,𝑎₆} can be calculated. Given another set of estimated parameters {𝑎₁,𝑎₂,…,𝑎₆}, there obtains another value of the performance index (2.3). That is, the value of the performance index (2.3) is actually dependent on the set of estimated parameters {𝑎₁,𝑎₂,…,𝑎₆}. Therefore; it is assumed that there is no measurement error in this paper
The HTDE combines the traditional DE [26] with the Taguchi method [38, 39]. In the HTDE, the Taguchi method is inserted between the crossover operations of a DE. Then, by using two major tools (signal-to-noise ratio and orthogonal arrays) of the Taguchi method, the systematic reasoning ability of the Taguchi method is incorporated in the crossover operations to systematically select the better genes to achieve the crossover operations, and consequently enhance the DEs [32]. The detailed steps of the HTDE are described as follows. The detailed description of the Taguchi method can be found in the books presented by Taguchi et al. [38] and Wu [39].

Detailed Steps: HTDE
Step 1 (parameter setting). The population size 𝑝_𝑠, the number of generations 𝑁_𝑔, the positive real mutation scaling factor 𝐹, and the user-defined bounds [𝑥_𝑖,𝑥_𝑖](𝑖=1,2,…,6).Step 2 (initialization). In the real coding representation, each individual is encoded as a vector of floating point numbers, with the same length as the vector of decision parameters. For convenience and simplicity, an individual 𝑋^𝐺_𝑖 is denoted as 𝑋^𝐺_𝑖=[𝑥^𝐺_𝑖1,𝑥^𝐺_𝑖2,𝑥^𝐺_𝑖3,𝑥^𝐺_𝑖4,𝑥^𝐺_𝑖5,𝑥^𝐺_𝑖6], (2.5) where 𝐺 denotes the current generation, 𝑥^𝐺_𝑖𝑗(𝑗=1,2,…,6) denote the estimated parameters (i.e.,{𝑎₁,𝑎₂,…,𝑎₆}),𝑖=1,2,…,𝑝_𝑠. The initial values of 𝑥^𝐺_𝑖𝑗 are chosen at random from within user-defined bounds [𝑥_𝑗,𝑥_𝑗].Step 3. Generate a uniformly distributed random number 𝛼, where 𝛼∈[0,1].Step 4. Let 𝑥^𝐺_𝑖𝑗=𝑥_𝑗+𝛼(𝑥_𝑗−𝑥_𝑗), where 𝑥^𝐺_𝑖𝑗 denotes the 𝑗th element of 𝑋^𝐺_𝑖 in (2.5). Repeat six times and produce an individual vector [𝑥^𝐺_𝑖1,𝑥^𝐺_𝑖2,𝑥^𝐺_𝑖3,𝑥^𝐺_𝑖4,𝑥^𝐺_𝑖5,𝑥^𝐺_𝑖6].Step 5. Repeat the above two steps 𝑝_𝑠 times and produce 𝑝_𝑠 initial feasible individual vectors 𝑋^𝐺_𝑖(𝑖=1,2,…,𝑝_𝑠).Step 6. Initialize the iteration index 𝐼=1.Step 7. Set the initial target index 𝑖=1.Step 8 (mutation operation). In every generation, each individual vector 𝑋^𝐺_𝑖(𝑖=1,2,…,𝑝_𝑠) in the population becomes a target vector. For each target vector 𝑋^𝐺_𝑖, DE applies a differential mutation operation to generate a mutated individual 𝑉_𝑖^𝐺+1, according to 𝑉_𝑖^𝐺+1=𝑋^𝐺_𝑗+𝐹(𝑋^𝐺_𝑛−𝑋^𝐺_𝑞), (2.6) where 𝑋^𝐺_𝑗,𝑋^𝐺_𝑛, and 𝑋^𝐺_𝑞 are randomly selected from the population such that 𝑗,𝑛, and 𝑞 belong to {1,2,…,𝑝_𝑠} and 𝑖≠𝑗≠𝑛≠𝑞;𝐹 is a mutation scaling factor which controls the evolution rate. In the event that mutation causes an element 𝑣^𝐺+1_𝑖𝑗 in the vector 𝑉_𝑖^𝐺+1=[𝑣^𝐺+1_𝑖1,𝑣^𝐺+1_𝑖2,𝑣^𝐺+1_𝑖3,𝑣^𝐺+1_𝑖4,𝑣^𝐺+1_𝑖5,𝑣^𝐺+1_𝑖6] (i.e., {𝑎₁,𝑎₂,…,𝑎₆}) to shift outside the allowable interval [𝑥_𝑗,𝑥_𝑗], then 𝑣^𝐺+1_𝑖𝑗 is set to 𝑥_𝑗 if 𝑣^𝐺+1_𝑖𝑗<𝑥_𝑗, and 𝑥_𝑗 if 𝑣^𝐺+1_𝑖𝑗>𝑥_𝑗. The mutated vector 𝑉_𝑖^𝐺+1 will be used in the following Taguchi-method-based crossover operation as a donor vector for breeding a better offspring or trial vector 𝑈_𝑖^𝐺+1.Detailed Steps: Taguchi-Method-Based Crossover Operation
Substep 8.1. Set 𝑚=1. Take the three columns of the orthogonal array 𝐿₄(2³) to allocate the three factors 𝑋^𝐺_𝑗,𝑋^𝐺_𝑛, and 𝑋^𝐺_𝑞 with two levels (individuals), respectively.Substep 8.2. Calculate the performance index 𝐽 with (2.3) and the signal-to-noise ratio with smaller-the-better characteristic for the new individual 𝑉_𝑖^𝐺+1 using (2.6).Substep 8.3. If 𝑚>4, then go to Substep 8.4. Otherwise, 𝑚=𝑚+1, return to Substep 8.2, and continue through Substep 8.3.Substep 8.4. Calculate the effects of the various factors.Substep 8.5. One optimal trial vector 𝑈_𝑖^𝐺+1 is generated based on the results from Substep 8.4.Step 9. Add one to the target index 𝑖. If 𝑖<𝑝_𝑠, then go to Step 8.Step 10. Add one to the iteration index 𝐼. If 𝐼<𝑁_𝑔, then go to Step 7.Step 11. Stop and display the optimal individual vector (i.e., the optimal parameters {𝑎₁,𝑎₂,…,𝑎₆}) and its performance index 𝐽.

3. Results and Comparisons

In this section, we evaluated the performance of the proposed HTDE approach to an HIV dynamical model given by Manseur et al. [22] and compared the results obtained from our proposed HTDE approach with those from the traditional DE approach [26], and the numerical methods with combined Adomian/Alienor approach and the classical optimization method with Levenberg-Marquardt approach given by Manseur et al. [22] by using the same example. Besides, for the sake of fair comparisons between our proposed HTDE approach and the traditional DE approach, the same evolutionary environments used in this paper are the population size 𝑝_𝑠=100, the same maximum number of generations 𝑁_𝑔=100, the mutation scaling factor 𝐹=0.8, and the same search interval [𝑥_𝑗,𝑥_𝑗]=(0,2](𝑗=1,2,…,6) for the estimated parameters {𝑎₁,𝑎₂,…,𝑎₆}, respectively.

The actual parameters of an HIV dynamical model with initial condition [𝑥₁(0)𝑥₂(0)𝑥₃(0)]=[10.20.8] that will be estimated are assumed to be 𝑆=1,𝑑=0.8,𝛽=1,𝜇₁=0.8,𝑘=1, and 𝜇₂=0.01078 [22]. In the following, we will apply the proposed HTDE to find the optimal parameters {𝑆,𝑑,𝛽,𝜇₁,𝑘,𝜇₂} (i.e., {𝑎₁,𝑎₂,…,𝑎₆}) such that the performance index 𝐽 (2.3) is minimized, where the sampling time is 0.1 second, and the total sampling number 𝑇 is 10.

For the HIV dynamical models (2.1a)–(2.1c) and (2.2a)–(2.2c), after using the proposed HTDE to execute 100 independent runs, we can get the best values, the mean values, and the standard deviation of the obtained data for the estimated parameters, the performance index and the computational time with Intel(R) Core(TM)2 Duo CPU 2.80 GHz and 2.00 GB RAM as shown in Tables 1 and 2, respectively. On the other hand, if we employ the traditional DE to execute 100 independent runs, we can also get the best values, the mean values, the standard deviation of the obtained data for the estimated parameters, the performance index and the computational time with Intel(R) Core(TM)2 Duo CPU 2.80 GHz and 2.00 GB RAM as shown in Tables 1 and 2, respectively. Besides, the estimated parameters and the performance indices are obtained from two transformations by using the combined Adomian/Alienor method and a class optimization method with the Levenberg-Marquardt method [22], respectively, are also shown in Table 1. Figures 1, 2, 3, 4, 5, 6 and 7 show the average convergence results of performance index 𝐽 and estimated parameters {𝑎₁,𝑎₂,…,𝑎₆} with respect to the number of generations in 100 independent runs by using the proposed HTDE and the traditional DE, respectively, for the HIV dynamical model. The responses of the 𝑥₁(𝑡),𝑥₂(𝑡), and 𝑥₃(𝑡) for the HIV dynamical model with the mean values of the estimated parameters in 100 independent runs obtained from the proposed HTDE and the traditional DE are, respectively, shown in Figures 8, 9, and 10.

Table 1 

Comparison of results for the estimated parameters and the performance index with 100 independent runs for the HIV dynamical model.

Table 2 

Comparison of results for the computational time (in minute) with 100 independent runs for the HIV dynamical model.

Figure 1 

Average convergence results of performance index in 100 independent runs by using the HTDE and the DE, respectively, for the HIV dynamical model.

Figure 2 

Average convergence results of 100 independent runs by using the HTDE and the DE, respectively, for the estimated parameter 𝑎1 of the HIV dynamical model.

Figure 3 

Average convergence results of 100 independent runs by using the HTDE and the DE, respectively, for the estimated parameter 𝑎2 of the HIV dynamical model.

Figure 4 

Average convergence results of 100 independent runs by using the HTDE and the DE, respectively, for the estimated parameter 𝑎3 of the HIV dynamical model.

Figure 5 

Average convergence results of 100 independent runs by using the HTDE and the DE, respectively, for the estimated parameter 𝑎4 of the HIV dynamical model.

Figure 6 

Average convergence results of 100 independent runs by using the HTDE and the DE, respectively, for the estimated parameter 𝑎5 of the HIV dynamical model.

Figure 7 

Average convergence results of 100 independent runs by using the HTDE and the DE, respectively, for the estimated parameter 𝑎6 of the HIV dynamical model.

Figure 8 

Responses of 𝑥1(𝑡) by using the HTDE and the DE, respectively, for the HIV dynamical model.

Figure 9 

Responses of 𝑥2(𝑡) by using the HTDE and the DE, respectively, for the HIV dynamical model.

Figure 10 

Responses of 𝑥3(𝑡) by using the HTDE and the DE, respectively, for the HIV dynamical model.

From Table 1 and Figures 1–10, the following results can be observed: (i) the proposed HTDE can obtain better results including the best values of estimated parameters and performance index, and the mean values of estimated parameters and performance index than the traditional DE and the numerical methods presented by Manseur et al. [22]; (ii) the proposed HTDE gives a smaller standard deviation of the obtained data for the estimated parameters and the performance index than the traditional DE, so the proposed HTDE has a stable solution quality; (iii) the average convergence results for the proposed HTDE are better than those of the traditional DE in 100 independent runs regarding the estimated parameters and the performance index as well as the proposed HTDE has a better solution convergence quality; (iv) although the computational time of the proposed HTDE is about four times the traditional DE, the better estimated parameters and performance index are concerned rather than the computational time in the problem of parameter estimation for an HIV dynamical model. Thus, it can be concluded that the proposed HTDE approach can give a more effective and robust way for finding the actual parameters than the traditional DE approach [26] and the numerical methods presented by Manseur et al. [22] for the HIV dynamical model.

4. Conclusions

An HTDE approach has been presented in this paper based on the Taguchi-method-based crossover operation for solving the problem of parameter estimation for an HIV dynamical model under the minimization of a performance index (2.3). The HTDE can be easily employed to find the system parameters of an HIV dynamical model. An illustrative example regarding an HIV dynamical model has shown that the proposed HTDE approach is effective and robust to estimate the system parameters of an HIV dynamical model, and outperforms the traditional DE approach [26] and the numerical methods presented by Manseur et al. [22]. Although this HTDE approach is proposed under the framework of HIV dynamical model, it should be generally applicable to any optimization problem.

Acknowledgment

This work was in part supported by a Grant from the Chi-Mei Medical Center and Kaohsiung University Research Foundation (99CM-KMU-01) and the National Science Council, Taiwan, Republic of China, under Grant number NSC 99-2320-B-037-026-MY2. The authors also thank the editor and the reviewers for their suggestions and comments.