Between-host HIV model: stability analysis and solution using memetic computing

: HIV poses a great threat to humanity for two major reasons. First it attacks the immunity system of the body and second, it is epidemic in nature. Mathematical models of HIV have been instrumental in understanding and controlling the infection. In this paper, we solve the between host epidemic model of HIV, described by nonlinear coupled differential equations, by using memetic computing. Under this model, the sexually active population is divided into four classes and we investigate the transfer of individuals from one class to another. The solution consists of Bernstein polynomials whose parameters have been optimized by using differential evolution as global and sequential quadratic programming as local optimizer. Our second contribution is the stability analysis of this model. The disease-free equilibrium is stable while endemic equilibrium is unstable within the practical range of the values of parameters.


Introduction
HIV infection is considered among the world's biggest health problems affecting millions of people around the world. HIV infects various body cells in the immune system, particularly CD 4 + T cells. Without any treatment, this infection leads to the failure of the host immune system and ultimately to death [1]. The spread of infection is considered to be a vital part of any infectious disease where the dynamics of the disease primarily depends on the rate of transmission from an infectious individual to a susceptible individual [2]. The deadly and epidemic nature of HIV has several biological characteristics and has attracted epidemiologists, mathematicians and biologists resulting in various mathematical models that have proven to be valuable for understanding the evolution of the disease, the epidemiological patterns of HIV, especially the mechanism associated with the spread of the disease [3].
The HIV transmission models are used to specify the dynamics of the spread of the disease by a system of equations. In these equations, the transition rates between defined states are described quantitatively [4]. To describe the HIV transmission and simulate the natural history of infection and disease process, these models use various biological and behavioral variables. By using the valid models and accurate datasets, they are capable of reliably projecting the infection rates [5].
To describe the HIV dynamics, three broad classes of mathematical models [6], namely, immunological, * Correspondence: musharif.ahmed@riphah.edu.pk This work is licensed under a Creative Commons Attribution 4.0 International License. epidemiological, and immunoepidemiological models have been developed. At individual level, the within-host viral dynamics are described as immunological models, while at the population level, transmission dynamics is explained using between-host transmission and are categorized as epidemiological models. The dynamics of population-level HIV transmission as a function of immune viral responses at individual level are described using the immunoepidemiological models [7][8][9].
The epidemiological models of HIV infections categorize each individual as susceptible or infected. Infected individuals are assumed to transmit HIV virus to susceptible individuals with the same transmission rate in the duration of disease [10,11]. It is also studied that host biological factors such as age, sex, genetic susceptibility, time since infection, and status of immune system may cause variation in infectiousness of infected individuals [12]. Various models belonging to these three categories have been under continuous investigation and modification.
These models are described using coupled nonlinear differential equations and their exact solutions do not exist. As a consequence, researchers attempted to come up with nonanalytical solutions by employing various numerical techniques, which are largely deterministic in nature. However, the quest for faster, more stable, and more accurate algorithms is still in progress. Recently, the abundance of the computing power has motivated researchers to use the evolutionary algorithms for the solution of various research problems in engineering and medicine.
In the last two decades, a variety of approaches have been proposed to solve various models of HIV infection. For example, to solve the immunological models, different techniques have been used, such as homotopy analysis method [13], quasilinearization-Lagrangian method [14], the shifted boubaker Lagrangian approach [15], Laplace-Adomian decomposition method [16], and Legendre wavelet method [17]. The epidemiological models are studied both from stability analysis perspective, as well as, their numerical solutions [1,[18][19][20][21]. In this paper, we propose Memetic computing approach to solve the HIV epidemic model. The use of memetic computing is a recent phenomenon which has become popular in the last decade [22][23][24][25][26]. The cost function which is an error function is not at all convex. It has many local minima and one global minimum. Such a situation demands the use of algorithms like genetic algorithm or differential evolution to search for a global minimum as they do not get stuck in a local minimum. This global search method is then integrated with a local method like sequential quadratic programming (SQP) or interior point algorithm (IPA) to intensify and hasten the search in the narrow areas and reach almost close to optimal solution i.e. global minimum. We • For solving this system of nonlinear coupled differential equations representing between host HIV/AIDS model, we have used global search method as differential evolution (class of evolutionary algorithms) that is good for exploration over a wide space of local minimas and a global minimum. It has been hybridized with sequential quadratic programming which is a fast convergent local search method thus achieving the global minimum. Moreover, we have used Bernstein polynomials whose different linear combinations give us solutions for susceptible population, infected population, and population with AIDS.
• We have analyzed the stability of the undiseased equilibrium, as well as endemic equilibrium. We have found that with the realistic range values of the parameters, equilibrium without disease is stable. However, the endemic equilibrium is unstable by using the values of parameters within the realistic range. This result contradicts that in [19].

Between-host HIV epidemic models
Various epidemic models have been developed to describe the transmission dynamics at population level [27].
The epidemic models of HIV differ in many aspects; ranging from different classifications of the total population to various terms used to represent transmission rates. Total population is generally divided into susceptible individuals, infected individuals and individuals with AIDS [1,19,20]. The infected individuals may further be divided into the individuals who are aware of their infection and who are not aware of their infection. The transmission rate is considered constant in some models while other models may differ in terms of transmission rates depending upon stages of HIV infection [28][29][30].
The simplest form of a susceptible-infection model is given in Eqs. (1) and (2): We have considered the model that divides the infected population into two groups, I u , who do not know they are infected and I a , who know that they are infected. This model is given by Eqs. (3)(4)(5)(6) which includes more realistic infection rate term and the individuals with AIDS infection: In this model, the total sexually-active population at time t, is divided into four mutually exclusive groups to aids group, A, while δ 2 is the transfer rate of I a individuals to aids group, A. B is recruitment rate of sexually active adults to the population per year. All of the parameters are summarized in Table 2.

Parameter Description Values
Infection rate of I u group 0.28 [31] β 2 Infection rate of I a group 0.35 [32] ω Transfer rate of individuals of I u group into I a group 0.9 [32] δ 1 Transfer rate of individuals of I u group into A group 0.7 [19,20] δ 2 Transfer rate of individuals of I a group into A group 0.4 [19] µ Natural mortality rate 0.000027 [19] d Disease-induced death rate 0.00002 [20] B Susceptibles adding per year 30

Problem formulation
Since the exact solution to the set of coupled nonlinear differential equations, does not exist, the researchers have attempted the solution through numerical as well as heuristic techniques. We propose to employ the memetic computing approach which uses both global and local optimizers for getting the suboptimal solution.
Specifically, we have used differential evolution (DE) algorithm as global optimizer and as local optimizer, SQP is used. DE belongs to the class of evolutionary algorithms and was proposed by Rainer Storn & Kenneth Price in 1997 [33] for optimization problems in continuous domain. It has some distinguishing properties like simple structure, speed, ease of use, and robustness [34]. Differential evolution mainly relies on mutation while GAs rely on both crossover and mutation [35,36].
The proposed methodology comprises two major steps. The first step is the formulation of fitness (error) function. Bernstein polynomials are linearly combined to represent the solution and are substituted in the differential equations to give us fitness function. In the second step, the coefficients of the linear combination are optimized by employing the memetic computing so as to give us minimal possible error function.

Bernstein polynomials
The Bernstein polynomials [37] are generally used for the approximation of continuous functions, uniformly, on a closed interval of interest. These polynomials have some peculiar properties, such as positivity, continuity, excellent numerical stability, uniform convergence, partition of unity property, and expansion in terms of power basis, that differentiate them from other polynomials. Bernstein polynomials are defined as follows: with the properties that Bernstein polynomials provide greater flexibility to levy boundary conditions at the end points of the interval.
B i,n at t = 0 and t = 1 is defined as follows:  Bernstein polynomials have inherent and higher level of agility and provide greater flexibility to approximate any function due to its existence in the interval [0, 1] as shown in Figure 1. Bernstein polynomials have recursive formulae, whereby the higher order polynomials can be expressed in terms of lower-order polynomials. The same is true for the first-and higher-order derivatives.
The first step is to express S(t) , I u (t), I a (t) , and A(t) as linear combination of Bernstein polynomials.

Fitness or error function
The set of differential equations need to be satisfied at every point on the grid from 0 to 1. This initiates an error function for all differential equations given in Eqs. (3)(4)(5)(6). Similarly, We are interested in minimizing the sum of these four errors:

Sequential steps of the hybrid algorithm
For achieving optimal values of the coefficients, we propose the use of DE for global optimization and sequential quadratic programming as local optimization. DE like other heuristic computing algorithms (HCAs) are efficient to explore global search space because of their capability to quickly explore and discover desired results in the search space. The evolutionary algorithms, like DE, has the useful property of exploring the wide search space and not getting stuck in local minima [38]. However, it has the tendency of slow convergence to the optimal solution [39,40]. Therefore, it is hybridized with a local search method, like SQP or IPA, which has the probability of search intensification and quick convergence to near optimal solution in the global minimum. This creates a good balance between exploration and exploitation.
Step 1: Initialization: In general, the differential evolution algorithm starts with the initialization of the candidate solutions and it may be biased or random initialization. In our case, it is a random initialization with real numbers between 10 and -10. For our discussion, we assume that x i,g (j) is the j th element of the i th candidate solution in the generation g .
Step 2: Fitness evaluation: Before we go for applying DE operators, we have to evaluate fitness of each candidate solution. This is done by using Eq. (16) for each candidate within the initial generation. The lower the error E is, the more the fitness is.
Step 3: Mutation: A mutant vector is created by selecting three random numbers s.t. n 1 ̸ = n 2 ̸ = n 3 ̸ = i as follows: The mutant vector is given v g+1 where v g+1 i and n1 , n2, and n3 are random numbers within the range of 1 to size of population such that n1 ̸ = n2 ̸ = n3 . µ ∈ (0, 2] and it is important to mention that a large mutation factor µ expands the search space but results in slowing down the convergence.
Step 4: Recombination: The parent vector is combined with the mutated one to generate a trial vector.
Step 5: Selection: Step 6: Termination criterion: The algorithm is stopped if either error occurs or the number of iterations reached.
Step 7: Initialization of SQP: After the termination of DE, the candidate solution with the best fitness is selected and handed over to SQP which is one of the best local optimizer, in terms of fast convergence, for the numerical solution of nonlinear optimization problems.
Step 8: Fitness Evaluation: Fitness evaluation is performed by using eq. (16) as in DE.
Step 9: Termination: After the specified number of iterations the SQP algorithm is stopped and the final result (solution) is obtained.

Results and discussion
For the HIV interhost model described by coupled nonlinear differential equations, in Section 2, we have carried out simulations in MATLAB for various values of model parameters. These parameter values have been explained in Table 2, and their values have been taken from the contemporary research work. Figure 2 shows the impact of increasing the transfer rate between I u and I a . As we increase the transfer rate between the I u group and I a group, the number of susceptible individuals increases over the specified time interval, which is 100 days. For the highest transfer rate ω , we have the highest conversion rate from the I u group to the I a group.  Figure 3. Similarly if we keep δ 1 and δ 2 fixed and increased ω from 0.3 to 0.6, the number of susceptible individuals increases again. In summary, there is an increasing trend from I u to I a , I u to A, and I a to A as the parameter of transfer rates ( δ 1 , δ 2 , and ω ) are increased. Figure 4 shows the trend in the number of susceptible individuals in the typical simulation scenario which is exponential in nature, after the mid of the observation interval.     and δ 2 shows a considerable increase in the number of unaware infected individuals and this increase goes up to an order of magnitude. Similarly, keeping δ 1 and δ 2 fixed, and increasing ω results in an increase in the number of unaware infected individuals. In Figure 6, we have fixed the values of δ 1 and δ 2 and investigated the impact of increasing ω on the number of unaware infected individuals. The impact of changing ω starts becoming visible after the mid of the observation interval i.e. after 50 days and it becomes increasingly different by the end of the observation interval. In fact, the impact of different values of ω results in a nonlinear increase in the number of unaware infected individuals. Figure 7 shows the pattern of increase in the number of unaware infected individuals in the typical simulation scenario which is exponential in nature. Figure 8 shows the impact of changing ω by keeping δ 1 and δ 2 fixed. This is clearly a nonlinear relationship after the mid of the observation interval. Figure 9 shows the variation in I a , the number of people who know about their infection, with different values of the parameters. For fixed value of ω (transfer rate between I u and I a ), if we increase in δ 1 and δ 2 together, the graph shows a considerable increase in the number of aware infected individuals and this increase. Next we fixed δ 1 and δ 2 , and increased ω which resulted in an increase in the number of unaware infected individuals.   Figure 10 shows the pattern of increase in the number of aware infected individuals in the typical simulation scenario which is found to be exponential in nature.
In Figure 11 we can observe that the number of AIDS infected individuals increase as the parameter ω is increased which is quiet obvious that ω is the transfer rate between unaware infected and aware infected population.    Figure 13 shows the increase in the number of AIDS infected individuals in the typical simulation scenario and it is found to be exponentially rising. parameters ω and δ , used for getting the simulation results, shown above. In just 15 iterations, the mean square error becomes lower than 10 −7 , which is remarkable.
In Figure 15, ω is varied and other parameters have been kept constant. Similarly, for a complete set of parameters, Figure 16 shows the error curve which goes down from 2.5 to close to zero in 15 iterations.

Equilibrium points and stability analysis
For carrying out the stability analysis of equilibrium points, we have taken I(t) = I u (t) + I a (t) to make it simple. So the equations become as given below. For finding equilibrium points we put the derivatives equal to zero.
Here δ represents transfer rate between I and A .

Equilibrium points
Disease-free equilibrium: Hence, endemic equilibrium is X e = (S e , I e , A e ), where S e , I e , and A e are given in Eqs. (26)-(28), respectively.

Stability analysis of the two equilibria
First we carry out stability analysis of disease-free equilibrium X f = (S f , 0, 0). We find the Jacobian of this nonlinear system.
All the rates will be taken as per day. Hence, equilibrium point without disease is stable.
Let us look at λ 2 and λ 3 whether they are negative or not.
After inserting S e and I e from Eqs. (26) and (27) respectively, we get: Using the minimum values of µ and δ on RHS and maximum values of parameters on LHS, this inequality is still far from being satisfied.
Hence, λ 3 > 0 , which means that endemic equilibrium is not stable at all.

Conclusion
The simulation has been carried out for a range of infection and transfer rates and the results have shown that the infection may be controlled by governing the transfer rates from the infected class to the AIDS class. This may be achieved to some extent by applying appropriate drug therapy. A more effective way of controlling the infection is through public awareness. We also carried out stability analysis. The disease-free equilibrium was stable but endemic equilibrium was not stable in the whole practical range of parameters. Memetic computing has been employed to solve the interhost HIV infection model. The evolution of susceptible, infected and AIDS patients have been deduced from the mathematical model by using memetic computing and Bernstein polynomials. The mean square error has been brought down to 10 -7 in just 15 iterations, which shows excellent convergence of the memetic computing.
One of the future research direction may be to improve the HIV epidemiological mathematical model by incorporating the impact of movement of the infected individuals from one geographical location to another. Yet another research work may be carried out to assess the impact of awareness campaign on the spread of the infection.