ABSTRACT

The aim of this study is to develop two-dimensional cellular automata model of HIV infection that depicts the dynamics involved in the interactions between acquired immune system and HIV infection in the peripheral blood stream. The appropriate biological rules of cellular automata model have been extracted from expert knowledge and the model has been simulated with determined initial conditions. Obtained results have been validated through comparing with the accepted AIDS reference curve. The new rules and states were added to the proposed model to show the effects of applying combined antiretroviral therapy. Our results showed that by applying RTI and PI drugs with maximum drug effectiveness, comparing with cases in which no treatment was applied, the steady state concentrations of healthy (infected) CD₄⁺T cells were increased (decreased) 53% (41%). Also, the use of cART with maximum drug effectiveness led to a 69% reduction in the steady state level of viral load. At this time, obtained results have been validated through comparing with available clinical data. Our results showed good agreement with both reference curve and the clinical data. In the second phase of this study, by applying genetic algorithms, a therapeutic schedule has been provided that its use, while maintaining the quality of the treatment, leads to a 47% reduction in both drug dosage and the side effects of antiretroviral drugs.

Keywords:

HIV Infection, Cellular Automata Model, Combined Antiretroviral Therapy, Genetic Algorithms

1. Introduction

Human Immunodeficiency Virus (HIV) is a subgroup of retrovirus that causes the Acquired Immune Deficiency Syndrome (AIDS), a condition in humans in which progressive failure of the immune system allows life-threatening opportunistic infections and cancers to thrive. HIV infects vital cells in the human immune system such as helper T cells (specifically CD₄⁺T cells), macrophages and dendritic cells [1] . Figure 1 shows HIV life cycle. According to the World Health Organization, millions of people die each year from Acquired Immune Deficiency Syndrome (AIDS) [2] . That’s why most of research teams have focused their studies on AIDS and their knowledge about Human Immunodeficiency Virus (HIV) has been improved significantly. Although AIDS is not treatable definitely, it is controllable [3] . Antiretroviral drugs have been effectively applied to control the progression of disease in the HIV infected individuals because they increase the time interval between entry of HIV and onset of AIDS. Two of the most important categories of antiretroviral drugs for AIDS control are Reverse Transcriptase Inhibitors (RTI) and Protease Inhibitors (PI). The RTI antiretroviral drugs interfere with ability of HIV enzyme reverse transcriptase to convert HIV RNA into HIV DNA. Thus they can stop the HIV replication process. The PI drugs prevent the assembly of particles of new virus by interfering with HIV enzyme protease. Since by transforming the HIV as a result of mutation, some strains of HIV will be resistant to the type of drug used, therefore, the chance to control the HIV will be increased if the patient uses multiple types of drugs simultaneously. One of the conventional methods in HIV control is Combined Anti Retroviral Therapy (cART) in which the patient uses combination of RTI and PI antiretroviral drugs [4] .

By finding an appropriate model for a biological system, experts will be needless to do many relevant clinical trials. If we also consider the effects of drug therapy in such models, we will be able to introduce the best treatment by simulating the effects of different therapeutic methods. Although antiretroviral therapy has significant effects on prolonging survival of patients, the side effects of long-term applying of these drugs may put the patient at risk [5] . Therefore, in this study, after development of an appropriate model of HIV infection in the peripheral blood stream of patients who are submitted to antiretroviral therapy, a therapeutic schedule is provided that its use causes to a significant decrease in drug dosage. As a result, by applying obtained therapeutic schedule, while maintaining the quality of treatment, the side effects of antiretroviral drugs will be decreased significantly.

Most of researchers have used models with Ordinary Differential Equations (ODEs) or Partial Differential Equations (PDEs) to simulate the dynamics of HIV infection [6] [7] [8] . The biological systems such as human immune system are so complicated and their components interact with each other through specific rules that are almost impossible to depict such rules using ODEs approach. On the other hand, ODEs modeling approaches are unable to depict the specific properties of immune systems such as memory and emerging [9] . Also, especially when the drug is administered to control the disease, behavioral variations, and special states arise which ODEs modeling approaches are unable to depict all of these variations and they can’t consider appropriate parameters for showing all of these states. Because of the reasons mentioned, in this study, instead of using ODEs approaches we tried to use a qualitative modeling method based on rules to simulate the complicated and nonlinear dynamics of a human immune system and its response to the virus entered the human body.

The Cellular Automata (CA) modeling approach is a perfect choice for this type of modeling techniques. Cellular Automata (CA) is discrete dynamical systems for which its behavior is based entirely on local communications [10] . The main idea of CA is that each region of grid is considered as a cell, and specific status is assigned to that cell. In each iteration, the state of each cell will be updated, according to own state and the neighboring cells state. The neighboring cells of one cell include the cell itself and its nearest adjacent cells. These states are correlated to each other through certain rules. These rules are obtained from the biological and physical rules which are responsible for guiding the system behavior. Since in most biological systems a large degree of uncertainty is observed, therefore, for such systems, using the deterministic rules is not reasonable. Thus in modeling the biological systems with CA method, sometimes it will be necessary to consider probable rules [11] . CA modeling power is from its simplicity and at the same time its ability to model complex systems. Also, due to the direct use of real biological rules governing the issue, there will be the possibility of full compliance of model to the reality of issue [12] . Some researchers have developed different CA models to explain the dynamics of HIV infection and most of them have focused on lymphoid tissue. Dos Santos used two-di- mensional CA modeling approach to model the dynamics of HIV proliferation and infection in the lymphoid tissue [12] . They dispersed a fraction of infected CD₄⁺T lymphocytes as the initial infection among of the healthy CD₄⁺T lymphocytes and depicted the interactions between these healthy and infected lymphocytes through a set of simple rules. Then they validated their model through comparing the behavioral quality of their obtained results with a reference curve which shows the natural history of HIV infection dynamics currently accepted. This common pattern is depicted in Figure 2. The study was conducted by dos Santos was the basis of many articles published in the field of modeling the dynamics of HIV infection in lymphoid tissue using CA approach [13] [14] [15] [16] .

In the previous mentioned literature, the dynamics of HIV infection have been studied in the lymph node. However, whenever most physicians try to assess the progress of disease, they would prefer to use patient’s blood data because such data can be accessed more easily. Moreover, they can find all types of cells in the peripheral bloodstream. Because of these reasons, some of the authors selected peripheral bloodstream to model the dynamics of HIV infection using CA method [17] [18] . Jafelice et al. developed a CA model to simulate the evolution of HIV in the bloodstream of positive individuals [17] . They studied the effects of drug therapy in their model using a fuzzy rule based system with two inputs, the medication potency and patient’s treatment adhesion. They considered four types of cells namely uninfected cells, infected cells of lymphocytes T of CD₄⁺, free virus particles and specific antibodies CTL. These cells were placed at random positions initially in a two dimensional cell grid and the states of each cell in the grid were updated according to the defined local dynamic rules of each cell. Obtained results were validated through the reference HIV history pattern. It should be noted that their work was invaluable but they just had considered cellular response of immune system and they ignored the role of humoral response in the adaptive system. Moreover, the effects of intracellular biochemical factors that influence the effectiveness, such as the susceptibility of inhibiting

drugs was not considered. Khabouze et al., developed the study done by Jafelice by taking into account the role of humoral response as one of the basic part of adaptive immune system by considering a new type of cell called antibody but the effects of drug treatment in their model was not studied [18] .

In this study we model the interactions between HIV enters the human’s body and the immune system response, in the peripheral bloodstream, using two-dimensional CA technique. We also consider the effects of antiretroviral therapy in our model rules. Actually our CA model is an extended version of model proposed by Khabouze [18] by taking into account the following important points: a) the effects of either mono or combined antiretroviral therapy. For this purpose we considered new cell states associated with reverse transcriptase inhibitors and protease inhibitors drug absorption by the healthy cells of CD₄⁺T lymphocytes and we defined new rules to show how these inhibitor drugs can prevent the healthy cells from infection. We also investigated the effects of drug efficiency and initial time of onset of treatment in improving the quality of treatment. b) The phenomenon of drug resistance, c) The effect of latently infected cells in clinical latency phase of HIV infection and their important role in taking the patient to the onset of AIDS phase, d) The phenomenon of virus mutation, e) Infecting the healthy cells through contact with either free viruses or infected cells, while in the previous studies [17] [18] , authors considered that infection of healthy cells can only occur through free viruses. The simulation results obtained from our CA model reproduces three phase pattern of HIV infection and shows how the use of antiretroviral therapy is able to delay the arrival of the patient to the onset of AIDS phase. We validated our results through comparison with available clinical data.

In the next step, we used an intelligent optimization technique to provide an optimal therapeutic schedule that simultaneously keeps disease progression under control and allows the restoration of immune system. So far, in often models that have been used to depict the dynamic behavior of HIV, the time-continuous equations have been applied, and in this case, using optimal control methods (despite the high computational volume) were the best treatment strategy. However, since the CA model that used in this study is described by biological rules and there is no analytical equation, applying the classic optimal control methods is impossible. Therefore, in order to optimize treatment protocol and to provide optimal therapeutic schedule in terms of number and order of days of applying antiretroviral drugs and drug dosage, utilizing of algorithms which are search-based in discrete-time space, is recommended. Since the genetic algorithm is a good candidate of such algorithms, in this study we used genetic algorithms to provide the optimal therapeutic schedule.

After this introduction, this paper is organized as follow: in the next section, we introduce the proposed CA model by incorporating the antiretroviral therapy process into the rules of this model. Then we present the simulation results and we compare our obtained results with available clinical. After ensuring the validity of proposed model, in the next step, by utilizing the genetic algorithms, the optimal therapeutic schedule is obtained. In order to prove the effectiveness of the proposed treatment program, the results obtained from applying this optimum schedule is compared with the results obtained from two other therapeutic programs named “Full therapy” and “Random therapy”, respectively. In the last section we have concluded our work and suggestions for future works are expressed.

2. Material and Methods

2.1. Proposed Cellular Automata Model

The structure of proposed model in this study is two-dimensional CA model with Moor neighboring. For developing such a structure, we used two-dimen- sional cell grids with size of 250 × 250 to depict the bloodstream toroidal system.

Since CA models belong to subset of agent based modeling methods, defining the agents of proposed CA model is the first step of organizing this model. In this study, following agents have been defined:

1) Healthy cells (H): this is the state of a normal healthy cell which has not absorbed any type of antiretroviral drugs. The cell in this state is capable of becoming infected and this happens, it can infect other healthy cells in own neighboring.

2) Healthy cells resistant to infection (H_R): this is the state of healthy cells in which due to the absorption of both RTI and PI antiretroviral drugs, they will not become infected even if they are exposed to infectious factors.

3) Active infected cells (A₁): this state refers to the active infected cells which are able to infect the normal healthy cells in their own neighboring [12] .

4) Inactive infected cells (A₂): this state corresponds to the infected cells that after Ƭ₁ time delay, the immune system has developed own cellular response against it. With considering such time delay we can depict the mutation rate of HIV. In fact we consider that since each infected cell contains a new strain of HIV, thus the immune system needs time to identify each of them. On the other hand, as a result of cellular response of immune system, such infected cells alone are not able to infect the normal healthy cells in their own neighboring but local aggregation of them will cause the neighbor healthy cells to become infected. It means that if, at least, R numbers of A₂ cells have been located in the neighboring of a normal healthy cell, the healthy cell will become an active infected cell [12] .

5) Latently infected cells (A₀): this is the state of infected cells which have kept the latent infection in to themselves for a time and they are not capable of spreading the infection during this period of time. But they may be re-activated after a time delay Ƭ₂ and if it happens, they will be able to infect their own adjacent healthy cells [14] .

6) Inactivated infected cells in effect of PI drug absorption (A_PI): this state corresponds to the infected cells which have arisen due to infecting healthy cells that have absorbed PI antiretroviral drugs. The cells in this state are not able to spread the infection and they will not produce any new HIV strain.

7) HIV _Infectious: this is the state of active and infectious HIV particles which are capable of infecting their own adjacent normal healthy cells [16] .

8) HIV _{non-Infectious}: when a healthy cell that has absorbed PI antiretroviral drugs becomes infected, this type of infected cell produces free viruses that are often defective and are not capable of infecting other healthy cells. This state is corresponds to such detective and non-infectious viruses [16] .

9) HIV _Inactive: this state corresponds to viruses that the immune system has developed own humoral response against them. Therefore, due to the performance of immune system, the ability of such viruses is decreased for infecting other healthy cells [16] .

10) Dead cells (D): this state corresponds to both of the healthy cells and free viruses that are dead due to natural death. This state also contains the infected cells that will die in next step due to the performance of immune system against them [12] .

11) CTLs: the state of immune system cells that is responsible for identifying killing the infected cells [17] .

12) Anti bodies (AB): the state of anti bodies that are responsible for killing the viruses [18] .

In second step, we considered following percentages as an initial number of each of cell types in our model: 25% ± 10% of the entire cell grid was initially considered as normal healthy CD4⁺T lymphocytes [17] [18] . To depict the initial infection, we considered 0.05% of the entire cell grid as infected cells [12] . 5% ± 1% of the entire cell grid was considered as free infectious viruses [17] [18] . CTLs and anti bodies, each of them initially constituted 3% ± 1% of the entire cell grid [18] . Mentioned cells, initially were randomly distributed in the cell grid space.

2.1.1. The Rules of Proposed CA Model

As previously mentioned, in this study we have developed a CA model for modeling the HIV infection in the peripheral bloodstream with taking into account the effects of applying antiretroviral therapy. For this purpose, the effectiveness of RTI and PI drugs was considered by P_RTI and P_PI, respectively, and the value of them was selected between 0 to 0.9 (the value of 0 indicates no use of drug of interest and the value 0.9 indicates maximum effectiveness of the drug). We defined following rules for our proposed CA model based on neighboring. It is worth noting that in the defined rules where ever the phrase “in the presence of infecting factors” is used, it means that, at least, one active infected cell (A₁) or one infectious HIV or R numbers of inactivated infected cells (A₂) exist in the neighborhood of target cell.

Rule 1) The rules for updating the state of a normal healthy cell (H):

a) We have defined a lifespan limit for normal healthy cells and a random initial age (between 0 and the lifespan limit) has been assigned to each of H cells. At the beginning of each of iterations, if an H cell reaches its lifespan then it will die, otherwise the following updating rules will be used [17] [18] .

b) The H cell will choose randomly an empty place in the neighborhood and moves there for the next iteration. If there is no empty place then it remains at own place [17] .

c) For the next iteration, each of H cells, becomes an H_R cell with probability of P_RTI×P_PI in the presence or absence of infecting factors―it shows that when a normal healthy cell absorbs both of RTI and PI drugs simultaneously, it becomes a resistant healthy cell which none of the infecting factors are not able to infect that and it is protected from infection for longer period of time than the normal healthy cells which have absorbed only RTI drug;

Or remains an H cell with probability of P_RTI × (1 − P_PI) in the presence or absence of infecting factors and its age will be increased one unit?this shows that if a normal healthy cell only absorbs RTI drug, in the effect of this type of drug the virus prevents from entering to healthy cell. So the healthy cell will be protected from infection, but since it has not absorbed both RTI and PI drugs, this protection lasts only until the next iteration-;

Or with probability of (1 − P_RTI) × P_PI, two modes may occur: 1) in the absence of infecting factors, it remains an H cell and its age will be increased one unit, 2) in the presence of infecting factors, it becomes an A_P cell with age zero?it shows that since the healthy cell has not absorbed RTI drug, it will not protected from entering virus and, therefore, it will be infected in presence of infecting factors. But since it has absorbed PI drug, even if it becomes infected by infecting factors, the risen infected cell (A_PI) will not be able to infect other healthy cells-;

Or with probability of (1 − P_RTI) × (1 − P_PI) two modes may occur: 1) in the absence of infecting factors it remains at H state and its age will be increased one unit for the next iteration, 2) in the presence of infecting factors it becomes an active infected cell (A₁) with age of zero with probability of P_inf, or it becomes a latently infected cell (A₀) with age of zero with probability of (1 − P_inf) [14] [15] .

d) The normal healthy cells are proliferating with a constant rate during the simulation [17] [18] .

Rule 2) The rules for updating the state of a resistant healthy cell (H_R):

Since a resistant healthy cell has absorbed both of RTI and PI drugs, therefore, even if the infecting factors exist in its neighborhood, it will not become infected. This type of healthy cell becomes a normal healthy cell for the next iteration [15] [16] .

Rule 3) The rules for updating the state of active infected cells (A₁):

If an A1 cell is not identified by the immune system cells, it releases a new infectious HIV from within itself with probability of P_HIV-PR [16] [18] [19] . According to empirical findings, each of active infected cells may be identified by effective immune cells about 3 to 5 weeks after infection and loses its ability to spread the infection and becomes an inactivated infected cell (A₂) [12] [19] . In the previous studies, the time required to identify A₁ cells by the immune system cells was considered 4 weeks (the average of 3 and 5 weeks) [12] . But in this study, instead of averaging, we considered that there is the probability of identifying an A₁ cell by the immune system cells during the third until fifth weeks after infection of healthy cell. From the third week until fifth week, this probability is increasing. But due to the occurrence of mutation, an infected cell may contain a new strain of HIV into itself which has not previously been introduced to the immune system. For this reason, an infected cell may not be identified by the immune cells even five weeks after occurrence of infection. In this case, due to release of several HIVs from infected cell, its membrane will tear and thus the infected cell becomes dead cell for the next iteration [20] .

The state of A₁ cell is updated according to the following rules:

a) Until three weeks after infection, an A₁ cell produces HIV_Infectious with probability of P_HIV-PR or produces HIV_Inactive with probability of (1-P_HIV-PR) for the next iteration. The new released HIV is placed at one of empty sites in A1’s neighborhood. Also the age of an A₁ cell is incremented one unit for the next iteration [16] [18] .

b) Between three until five weeks after infection, following modes may occur:

1) If there is, at least, one CTL in the neighborhood, for the next iteration, an A₁ cell becomes A₂ with probability of P_{identification}; or produces HIV_Infectious with probability of (1 − P_{identificatin}) × P_HIV-PR; or produces HIV_Inactive with probability of (1 − P_{identificatin}) × (1 − P_HIV-PR). Also the age of an A₁ cell is incremented one unit for the next iteration [12] [16] [17] [18] [19] .

2) If there is no CTL in the neighborhood, for the next iteration, an A₁ cell produces HIV_Infectious with probability of P_HIV-PR; or produces HIV_Inactive with probability of (1-P_HIV-PR). Also the age of an A₁ cell is incremented one unit for the next iteration [16] [18] .

c) If due to mutations, an A₁ cell is not detected by the CTLs up to five weeks after infection, due to releasing new viruses from within it, finally the cell membrane of this A₁ cell will tear and it will die for the next iteration [20] .

d) Before an infected cell dies due to membrane rupture, if it is not detected at each of iterations by CTLs then for the next iteration it will select randomly an empty place in its neighborhood to move there and the age of this A₁ cell is incremented one unit. But if there is no empty place in its neighborhood then it remains at an own place and its age is incremented one unit for the next iteration [17] [18] .

Rule 4) The rules for updating the state of an inactivated infected cell (A₂):

An A₂ cell becomes a dead cell for the next iteration [12] .

Rule 5) The rules for updating the state of a latently infected cell (A₀):

a) Each of A₀ cells after a time delay T_{Re-activation} will be reactivated with probability of P_{Re-activation} and becomes a new A₁ cell with the age of zero for the next iteration; or remains an A₀ with probability of (1 − P_Reactivation) and dies for the next iteration [14] .

b) Before reaching the time of T_{Re-activation}, for the next iteration each of A₀ cells will randomly select an empty place in its neighborhood to move there and if there is no empty place, it remains at own place. Also the age of A₀ cell is incremented one unit for the next iteration.

Rule 6) The rules for updating the state of an A_PI cell:

As we mentioned before, whenever a healthy that has absorbed PI drugs becomes infected by infecting factors in its neighborhood, it becomes an A_PI cell with the age of zero for the next iteration. We also considered a lifespan limit for A_PI cells similar to what we had previously defined for A₁ cells (5 iterations). We defined following rules for updating the state of A_PI cells:

a) If an A_PI cell reaches its lifespan limit, it becomes a dead cell for the next iteration.

b) Before reaching the lifespan limit, following modes may occur for the next iteration:

1) It remains an A_PI cell and produces HIV_{noninfectious} with probability of P_PI × P_HIV-PR; or it remains an A_PI cell and produces HIV_Inactive with probability of P_PI × (1 − P_HIV-PR). The new released HIV is placed at one of the empty sites in A_PI’s neighborhood. In this case, the A_PI cell will randomly choose an empty place in the own neighborhood to move there and if there is no empty place, it remains at own place. Also the age of A_PI cell is incremented one unit for the next iteration.

2) It becomes an A₁ cell with the last age of A_PI cell and produces HIV_Infectious with probability of (1 − P_PI) × P_HIV-PR; or it becomes an A₁ cell with the last age of A_PI and produces HIV_Inactive with probability of (1 − P_PI) × (1 − P_HIV-PR). In this case, as we mentioned before, the arisen A₁ cell continues its life with the last age of corresponding A_PI cell and from now on, the rules that were defined for updating the state of A₁ cell will be used for it.

The rules defined for updating the state of A_PI cell is inspired from [4] [21] .

It is worth mentioning that the viruses which are produced by A_PI cells or by A₁ cells will be placed randomly in one of the neighboring places of infected cell producing their own. Depending on the type of the selected place, following modes may occur:

I. If this place is an empty place then the produced virus will be located there for the next iteration.

II. Since only the normal healthy lymphocytes are hosts for viruses, thus if this place be relevant to any state except H, the virus will not have no effect on the state of the place.

III. If this place is relevant to a normal healthy cell, since both of HIV_{non-Infectious} and HIV_Inactive are not capable enough to infect other cells, thus only HIV_Infectious particles can infect the H cell.

Rule 7) The rules for updating the state of dead cells:

A dead cell can be replaced by a new healthy cell for the next iteration. The probability of this replacement is shown by P_replication. In fact, P_replication depicts the ability of the immune system for reconstructing the destroyed cells. Among new produced healthy cells, a very small fraction of them contain the infection. The new produced cell may be infected with probability of P_{new-infection} [12] . Therefore, following rules have been defined for updating the sate of dead cells:

a) Each of dead cells becomes a new normal healthy cell (H) with probability of P_replication × (1 − P_{new-infection}) or becomes a new infected cell (A₁) with probability of P_replication × P_{new-infection} or remains a dead cell with probability of (1 − P_replication) for the next iteration [12] .

Rule 8) The rules of HIV proliferation and updating its state:

As mentioned before, HIV_Infectious is a virus particle that is able to infect own adjacent healthy cells. It means if there is, at least, one HIV_Infectious in the neighborhood of a normal healthy cell, this virus transcribes own RNA in to the DNA of host cell through reverse transcriptase enzyme and infect the healthy cell in this way. On the other hand, the rules numbers 3a, 3b and 6b show how one A₁ or A_PI cell causes to the proliferation of this type of HIV. We also defined following rules for updating the state of HIV_Infectious:

a) We defined a lifespan limit for HIV_Infectious particles that is shown by the parameter “V_LSL”. If one HIV_Infectious particle reaches the V_LSL, this virus particle will die due to the natural death for the next iteration [16] [17] [18] .

b) Before the HIV_Infectious reaches V_LSL, following modes may occur:

1) At each iteration, if there is, at least, one anti body cell (AB) in the neighborhood of an HIV_Infectious and detects this HIV particle with probability of P_{identification} then the HIV_Infectious will lose its ability for spreading the infection due to the performance of humoral response of immune system. Therefore, the HIV_Infectious becomes HIV_inactive for the next iteration [16] [18] [20] .

2) If the HIV_Infectious is not identified by the AB cell in its adjacent with probability of (1 − P_{identification}), or if there is not any AB cell in the neighborhood of this HIV_Infectious, then for the next iteration the HIV_Infectious will randomly select one of the empty places in own neighborhood to move there and if there is no empty place, it will remain at own place. Also, the age of this HIV_Infectious is incremented one unit for the next iteration [18] [20] .

As mentioned before, whenever a normal healthy cell that has absorbed PI drugs, becomes infected by infecting factors, the arisen infected cell (A_PI) produces the HIV particles which are often faulty and are not able to infect other healthy cells. We called this type of HIV particles HIV_{non-infectious}. The rule 7b shows how A_PI cells produce this type of HIV. We also defined following rule for HIV_{non-infectious}:

c) Each of HIV_{non-infectious} becomes HIV_inactive for the next iteration [16] .

Also we defined following rule for HIV_inactive:

d) Each of HIV_inactive becomes a dead cell for the next iteration [16] .

Rule 9) The rules for updating the state of CTLs:

The immune system develops own cellular response against infection through CTLs. This type of immune system cell is responsible for killing the infected cells. The rule number 3 shows how CTLs detect the infected cells and kill them. Following rules has been defined for updating the state of CTLs:

a) A lifespan limit for CTLs (CTL_LSL) has been defined. If a CTL reaches the lifespan limit then it becomes a dead cell for the next iteration. Otherwise, for the next iteration it will randomly select an empty place in the own neighborhood to move there and if there is no empty place then it remains at own place. Also, the age of CTL is incremented one unit for the next iteration [17] [18] .

b) At each of iterations, CTLs are being reproduced by a constant rate (CTL_PR).

Rule 10) The rules for updating the state of anti bodies:

The main role of antibodies (AB) is identifying virus particles and killing them. In rule 8b, we described how AB cells do this task. Following rules has been defined for updating the state of AB:

a) A lifespan limit for AB cells (AB_LSL) has been defined. If an AB reaches the lifespan limit then it becomes a dead cell for the next iteration. Otherwise, for the next iteration it will randomly select an empty place in the own neighborhood to move there and if there is no empty place then it remains at own place. Also, the age of AB is incremented one unit for the next iteration [18] .

b) At each of iterations, AB cells are being reproduced by a constant rate (AB_PR).

Figure 3 shows the block diagram of the interaction between the states of all agents.

2.1.2. The Parameters of Proposed Model

In order to simulate the proposed model, we used a set of parameters some of which have been extracted from reference articles and the rest of them are derived from expert knowledge and experimental findings. The expert knowledge includes medical reference books on AIDS [22] [23] , medical authoritative articles on AIDS [19] [20] [24] [25] [26] [27] and knowledge of infectious disease specialist. Table 1 shows the value of parameters used for simulating the proposed model.

2.1.3. Consideration of Drug Resistance in the Model

By using the antiretroviral drugs, the concentration of healthy CD⁴⁺T cells is increased and both concentration of infected CD⁴⁺T cells and viral load are decreased. But at the same time, due to occurrence of virus’s mutation and due to appearance of new strains of HIV, after a period of time from start of treatment, the patient becomes drug-resistant and the effectiveness of these drugs will be reduced. Due to the reduction of drug effectiveness, trend of increasing (reducing) of concentration of healthy cells (concentration of infected cells and viral load) stops and finally the concentration of these cells reaches to steady state level. The saturation of drug activity has been modeled through a negative exponential function by Caetano Marco [28] . Therefore, to describe drug-resistant with considering the interactions between drug effectiveness and concentration of healthy cells, infected cells and viral load, following exponential equations has been defined. Since the PI (RTI) drugs has more effects on the concentration of viral load (healthy and infected cells), we defined two equations separately for describing RTI and PI drug resistance.

(1)

(2)

In the above equation, P_PI (t) and P_RTI (t) denote PI and RTI drug effectiveness at iteration number t, respectively. P_0PI and P_0RTI denote the initial value of PI and RTI drug effectiveness, respectively. V(t − 1), I(t − 1) and H(t − 1) denote the viral load, concentration on infected cells and the concentration of healthy cells at iteration number t − 1, respectively. V(t_s), I(t_s) and H(t_s) denote viral load, concentration of infected cells and concentration of healthy cells at time of onset of treatment, respectively.

2.1.4. The Inputs and Outputs of Proposed Model

The RTI and PI antiretroviral drugs were considered as the inputs of our model. To consider these inputs, we used a parameter named effectiveness as P_RTI and P_PI for RTI and PI antiretroviral drugs, respectively [15] [29] . The time evolution of concentrations of normal healthy CD⁴⁺T cells, infected CD⁴⁺T cells (A₁ + A₂), CTLs, antibodies, dead cells and viral load, were considered as the outputs of our model. The curve of these variations has been drawn in the results section and each of them has been explained in detail.

2.2. Optimizing the Therapeutic Schedule by Applying Genetic Algorithms

After ensuring the validity of proposed CA model (according to the results obtained in Section 3.2), an optimal therapeutic schedule was presented by using Genetic Algorithms (GA). An introduction of GA can be found in [30] [31] . Since for HIV infected patients, long-term medication is prescribed, these patients suffer from side effects of antiretroviral drugs. Therefore, the most important advantage of using the obtained therapeutic schedule by GA is while maintaining the quality of the treatment, leading to a reduction in both drug dosage and the side effects of antiretroviral drugs. As we mentioned before, in order to prevent the patients from occurrence of drug-resistant, most of physicians prefer to apply cART method for HIV infected patients [32] . Also, as we will prove in section3.1, the use of cART method with maximal drug effectiveness, leading to better results (i.e. higher level of the steady state concentrations of healthy CD⁴⁺T cells, lower level of the steady state concentrations of infected CD⁴⁺T cells, as well as, lower level of viral load at steady state). For these reasons, the cART method with maximum drug effectiveness was used in our proposed therapeutic schedule. In presenting the optimal therapeutic schedule, in order to simply the simulations, we assumed that antiretroviral drugs (both RTI & PI drugs) are administered or interrupted weekly.

In GA, each of the therapeutic schedules was defined by a chromosome with m bits, that m defines the number of weeks of therapeutic schedule [33] . In this study, we assumed that each of therapeutic schedules contains a period of 32 weeks (i.e. weeks in 8 months) and the treatment was started at 160^th week. A chromosome is a 2-state variable. Therefore, for each of weeks (bits) of therapeutic schedules (chromosomes), following modes may occur:

a) If the ith bit of a chromosome is set equal to zero, then it means that during the ith week of therapeutic schedule, no treatment is applied.

b) If the ith bit of a chromosome is set equal to one, then it means that during the ith week of therapeutic schedule, cART method with maximum drug effectiveness is applied.

Therefore, one of the values 0 or 1, was assigned randomly to each of bits of chromosome. Then, corresponding to the assigned value, the values of drug effectiveness were adjusted in CA model. Accordingly, if the assigned value was equal to zero, then in CA model, the values of both P_0RTI and P_0PI was set equal to zero (i.e. P_0RTI = P_0PI = 0). And if the assigned value was equal to one, then in CA model, the values of both P_0RTI and P_0PI was set equal to 0.9 (i.e. P_0RTI = P_0PI = 0.9).

In the GA optimization, the population size was set equal to 32 × 10 = 320. The reason of this choice was that each of chromosomes has 32 bits, and 10 different states can be assigned to each of bits. Therefore, with setting the population size equal to 320, many different therapeutic schedules will be evaluated by GA and obtained optimal therapeutic schedule will be more reliable. The GA optimization procedure ran on 50 generations. The operations and parameters of genetic algorithms that were applied in our simulations are shown in Table 2.

The drug dosage and concentrations of healthy CD⁴⁺T cells and viral load at steady state due to applying a therapeutic schedule are the most important factors in evaluating the quality of that therapeutic schedule. Therefore, in order to evaluate the quality of each of chromosomes (therapeutic schedules), the following functions were calculated for each of chromosomes [33] .

(3)

(4)

(5)

In the above equations, t_s indicates the time of onset of treatment; t_e indicates the time of end of treatment., and indicate the concentrations of healthy cells, viral load and drug dosage, respectively, that obtained from applying the therapeutic schedule related to ith chromosome of jth generation. Then for ith chromosome of jth generation, the following cost function was defined:

(6)

A desirable therapeutic schedule is a schedule that by applying it, with minimum drug dosage, the concentration of healthy CD⁴⁺T cells and viral load at steady state reaches to maximum and minimum level, respectively. Considering this description and according to functions described above, that is more worthy chromosome for which the calculated value of the cost function is smaller (i.e.

. Then the fitness function for ith chromosome of jth generation was defined as inverse of cost function:

(7)

Since the fitness function is defined as the inverse of cost function, therefore, the highest calculated value of cost function is related to the fittest chromosome. According to GA optimization procedure, all of the chromosomes of first generation are sorted according to the value of cost function that was calculated for each of them. Chromosomes with smaller calculated value of cost functions were transferred to the next generation as the fittest chromosomes of first generation.

3. Simulation Results and Discussion

3.1. Results Obtained from Simulating the Proposed CA Model

In order to simulate the proposed CA model, a cell grid with size of 250 × 250 was considered. According to the defined initial numbers for each type of cells, an initial numbers of normal healthy cells, infected cells (A₁), infectious viruses, CTLs and antibodies were distributed randomly in the grid space. According to the defined rules in Section 2.1.1 and with considering the values of Table 1 for parameters, we simulated the model for 600 iterations (each of iterations was considered as one week). The simulation was run 20 times with 20 different random initial configurations. Then the average of results obtained from 20 times run the simulation was shown for each of iterations. First we assumed that any types of antiretroviral drugs are not used. For this purpose, the values of parameters P_RTI and P_PI were considered as zero. Figure 4 shows the results obtained from this simulation.

For validating the results obtained from proposed model, we compared our results with the standard reference AIDS curve which has been shown in Figure 2. As this Figure 4 shows, the results obtained from our model are qualitatively in a good agreement with the reference curve. In the progress of AIDS, the concentrations of healthy and infected cells are of special importance. As long as the concentrations of healthy cells are higher than the threshold of AIDS (i.e. 20% of maximum number of healthy cells before infection occurs), the risk of death not threaten the patient. Therefore, from now on only the concentrations of healthy and infected cells will be investigated.

In the next step, we studied the effects of applying medication on the concentrations of healthy and infected cells. As we mentioned before, due to occurrence of drug-resistance and as a result of HIV mutation phenomenon, combined antiretroviral therapy (cART) method (the therapy method in which a combination of both of RTI and PI drugs is used) will be more effective in control of progression of AIDS. Therefore, in this study, the effects of applying cART method, with consideration of different values for drug effectiveness, are studied. The concentrations of healthy CD⁴⁺T cells and infected CD⁴⁺T cells (A₁ + A₂)

obtained from simulating our CA model have been shown in Figure 5 and Figure 6, respectively. In our simulations, following conditions have been consi- dered: at first, we did not apply any medication and we called this situation as “without treatment”, then we considered applying cART method with minimum drug effectiveness (P_0RTI = P_0PI = 0.5), medium drug effectiveness (P_0RTI = P_0PI = 0.7) and maximum drug effectiveness (P_0RTI = P_0PI = 0.9). It should be noted that in all these conditions the treatment has been started from 160^th weeks. In order to normalize the obtained results, the following measures were taken: a) in Figure 5, the average number of healthy cells at each of iterations, divided into the initial number of healthy cells. b) in Figure 6, the average number of infected cells at each of iterations, divided into the initial number of healthy cells to determine what fraction of normal healthy cells have became infected cells.

As Figure 5 shows, in terms of without treatment, the concentration of healthy cells will decline below the threshold of AIDS after about 250 weeks. While, in terms of applying cART, the concentration of healthy cells at steady state, remain in higher level than the threshold of AIDS. Also, we can see from this figure that the greater drug effectiveness causes more increase in the level of concentration of healthy cells at steady state. When we considered minimum, medium and maximum drug effectiveness, in compare with the case in which no treatment method was considered, the concentration of healthy cells at steady state has been increased 28%, 46% and 53%, respectively.

As Figure 6 shows, applying cART method causes a significant reduction in the level of infected cells at steady state. When we applied the cART method with minimum, medium and maximum drug effectiveness, in compare with condition in which no therapy method has been applied, the level of infected cells have been decreased 24%, 37% and 41%, respectively.

Since the viral load plays an important role in spread of infection, the level of infectious HIV particles was studied under following conditions: a) without treatment, b) with applying cART method with minimum, medium and maximum drug effectiveness. Figure 7 shows the level of load under conditions mentioned above. As this figure shows, by applying cART method with minimum, medium and maximum drug effectiveness in compare with the case of “no treatment”, the level of viral at steady state was reduced 55%, 65% and 69%, respectively. Thus, the best result (i.e. the lowest level of infectious viral load) obtained in presence of combined antiretroviral therapy with maximum drug effectiveness. Therefore, Comparing the results obtained in case of “without treatment” with other conditions in which a therapy method has been used shows that antiretroviral drugs, in addition to increasing the concentration of healthy cells, is leading to a reduction in the level of infectious viral load in the patient’s peripheral blood.

3.2. Validating the Obtained Results from CA Model through Comparing with Clinical Data

Clinical data collection in the field of AIDS is not a simple task. Gonzalez et al.,

compared the results obtained from their CA model with clinical data which had been published in some of the authentic articles [15] . So in present work in order to validation of the results obtained from proposed CA model, we also apply the clinical data that had been used by Gonzales et al. In [35] nine patients were selected and were treated by applying cART. Treatment was started for these patients when the concentration of healthy CD⁴⁺T cells in their peripheral blood reached to 25% of concentration of this kind of cells in uninfected individual’s blood. These patients were followed over a period of eight months. The concentration of healthy CD⁴⁺T cells per microliter of patient’s peripheral blood was measured three times: once, right at the time of onset of treatment, second time, three weeks after onset of treatment and third time, eight months (24 weeks) after onset of treatment. In this reference, the average number of healthy CD⁴⁺T cells per microliter of an uninfected individual’s peripheral blood has been reported 970 ± 250 cells. In order to unify the scale of clinical data and results obtained from CA model, the clinical data were normalized by dividing them to average concentration of healthy CD⁴⁺T cells per microliter of an uninfected individual’s blood that has been reported in this reference. Figure 8(a) shows data relating to each of these patients. In order to comply with the conditions under which clinical data were recorded, we considered following assumptions: 1) In our model also, the treatment was started when the concentration of healthy CD⁴⁺T cells reached to 25% of initial concentration of this kind of cells. This

happened at 220^th week. 2) The simulations were ran nine times (to show nine patients) and at each time, the cART was started at 220^th week. During a period of 24 weeks, every three weeks, the concentration of healthy CD⁴⁺T cells obtained from our model was recorded. 3) This procedure was done once with consideration of minimum and once again with maximum drug effectiveness. Figure 8(b) shows the comparison between the average of results obtained from our model and average of values obtained from clinical data. As can be observed in Figure 8(b), the highest compliance of clinical data with results obtained from proposed CA model arises when the drug effectiveness (P_0RTI & P_0PI) are considered minimal.

In [36] , ten patients were selected and were treated by cART method. These patents were followed over a period of 168 day (24 weeks). The viral load per milliliter of patient’s peripheral blood was measured at following times: right at time of onset of treatment and respectively, 21 days (3 weeks) and 168 days (24 weeks) after onset of treatment. Figure 9(a) shows the recorded clinical data for each of 10 patients with an average of them at each measured time interval. In order to compare results obtained from our model with clinical data, following factors were considered in our model: 1) the simulation was ran 10 times (to show 10 patients) and at each time the cART therapy was started at 160^th week. 2) The concentrations of viral load obtained from CA model were recorded at same time intervals in which clinical data were measured. 3) This procedure was repeated once with minimum and once again with maximum drug effectiveness. Figure 9(b) shows the average of results obtained from our model in compare with real data. In order to unify the scale of clinical data and results obtained from CA model, we normalized both results and clinical data by dividing them to their values at the beginning of treatment. Therefore, we can show that the viral load at every week is what fraction of its value at the beginning of treatment. As both Figure 9(a) and Figure 9(b) show, clinical data, as well as results obtained from our CA model, reflect this fact that due to use of cART method the viral load is declined sharply about 3 weeks after onset of treatment. This proves that cART method can inhibit the viral load in patient’s peripheral blood quickly and successfully. As Figure 9(b) shows, the highest compliance of clinical data with our results arises when the drug effectiveness (P_0RTI & P_0PI) are considered maximal.

3.3. Results Obtained from GA Optimization

All steps of GA optimization were repeated for 50 generations. Figure 10(a) and Figure 10(b) show the average cost function and fitness function over generations, respectively. As Figure 10(a) shows, during the GA optimization over generations, the average cost function is decreasing and according to the inverse relation between cost function and fitness function, as Figure 10(b) shows, the average fitness function is increasing. Reduction in cost function and increase in fitness function during generation, are related to optimization of therapeutic schedule and leading to obtain the optimal schedule. In the last generation, the chromosome with smallest (biggest) value of cost function (fitness function) was presented as the best chromosome. Therefore, by applying the therapeutic schedule related to fittest chromosome, the highest concentration of CD⁴⁺T cells, as well as, the lowest viral load at steady state and simultaneously the lowest drug dosage will be obtained. Figure 11 shows the optimal therapeutic schedule obtained from GA optimization.

3.4. Validating the Results Obtained from GA Optimization

In order to prove desirable performance of optimal therapeutic schedule obtained from GA optimization, three different therapeutic schedules were defined and results obtained from applying these schedules were compared with each other. These therapeutic schedules are as below. In each of the therapeutic schedules the treatment was started from 160^th week.

a) Full therapy with maximum drug dosage: in this case the cART with maximum drug dosage (i.e. P_PI = P_RTI = 0.9) was applied at every week.

b) Best therapy with optimal drug dosage: in this case the optimal therapeutic schedule obtained from GA optimization was applied.

c) Random therapy: in this case, the number of weeks in which cART with minimum drug dosage is applied, was set equal to the number of therapy weeks in the case of best therapy, but a random therapeutic schedule was applied during these weeks.

Figure 12 shows a comparison between levels of drug dosage, as well as, the concentrations of healthy and infected CD⁴⁺T cells at steady state resulting from applying these three therapeutic schedules.

As Figure 12 shows, the steady state concentration of healthy CD⁴⁺T cells (infected CD⁴⁺T cells) obtained from “Best therapy” is just 3.9% (2%) lower (higher) than the steady state concentrations of these cells obtained from “Full therapy”. Also, the level of drug dosage in “Best therapy” is 47% lower than the level of drug dosage in “Full therapy”. Therefore, the steady state concentration of healthy and infected CD⁴⁺T cells obtained from “Best therapy” is very close to the steady state concentrations of these cells obtained from “Full therapy”, while, the drug dosage obtained from applying “Best therapy”, in compare with applying “Full therapy”, is decreased significantly. By comparing results obtained from optimum therapy schedule and results obtained from random therapy schedule we find that by applying random therapeutic schedule, although the level of drug dosage is identical with the case in which optimum therapy schedule has been applied, the steady state concentration of healthy (infected) CD⁴⁺T cells is 12.9% (9.4%) lower (higher) than the steady state concentrations of these cells obtained from optimum therapy.

Therefore, by concluding our results we find that by applying full therapy schedule, although the quality of treatment is maintained, due to the high level of drug dosage, the patient will suffer from side effects of antiretroviral drugs. On the other hand, by applying random therapy schedule, although the level of drug dosage has been decreased, the quality of treatment also decreased. But by

applying optimal therapy schedule that obtained from GA optimization, while maintaining the quality of treatment, leading to reduction in both, level of drug dosage and the side effects of antiretroviral drugs.

4. Conclusions

In this study, by using genetic algorithms (GA), we found an optimum therapeutic schedule for HIV infected patients based on combined antiretroviral therapy (cART) method. In order to simulate the effects of therapy, two-dimen- sional cellular automata (CA) model was used. This CA model depicts the dynamics involved in the interactions between acquired immune system and HIV, in the peripheral blood of HIV infected individuals. This model was developed taking into account the following points: a) occurrence of mutation phenomena, b) existing time delay in the immune system, c) the role of latently infected cells in spread of infection, d) the effects of antiretroviral drugs such as reverse transcriptase inhibitors (RTI) and protease inhibitors (PI), e) occurrence of drug resistance in patients who are under antiretroviral therapy.

For this purpose, first, we collected appropriate rules related to HIV infection in patient’s peripheral blood according to reference articles and base on expert knowledge including books and medical journals in the field of AIDS and we developed our CA model using these rules. At this step, by comparing the results obtained from our CA model with reference HIV curve, we showed that our results are in a good agreement with the existing reality in the interactions between acquired immune system and HIV. At next step, we studied the effects of cART by adding some new states and parameters to the proposed CA model. For this purpose, we considered a parameter named drug effectiveness for each of RTI and PI drugs that were shown by P_RTI and P_PI, respectively. Occurrence of drug resistance was considered by a negative exponential equation. The concentrations of healthy and infected CD4+T cells over the weeks and the level of viral load were obtained by consideration of different values for P_RTI and P_PI. Our results showed that by applying maximum effectiveness of RTI and PI drugs, compared with cases in which no treatment was applied, the steady state concentrations of healthy (infected) CD⁴⁺T cells were increased (decreased) 53% (41%). Also, the use of cART with maximum drug effectiveness led to a 69% reduction in the steady state level of viral load. In order to validate the results obtained from our proposed CA model in the cases in which cART was applied, we used clinical data which are provided in authoritative articles. By comparing results obtained from our model with clinical data, we prove that the best compliance between our results and clinical data arises when we consider minimum effectiveness for antiretroviral drugs.

In the second phase of this study, we used GA optimization to provide an optimal therapeutic schedule of cART with maximum drug effectiveness based on structured interruptions. To prove the effectiveness of obtained therapeutic schedule from GA optimization, three different therapy schedules were defined as following: a) Full therapy, b) Best therapy with optimum drug dosage, c) Random therapy. According to each of these therapy schedules, cART with minimum drug effectiveness was applied or interrupted at each of weeks in the proposed CA model. Then the steady state concentrations of healthy and infected CD⁴⁺T cells, and the amount of drug dosage obtained from applying each of these schedules were presented as our results. Our results showed that by applying the optimum therapeutic schedule obtained from GA optimization, the quality of treatment (i.e. high steady state concentration of healthy cells and low steady state concentration on infected cells) is very close when the full therapy schedule was applied. But the main point is that applying the optimum therapeutic schedule, compared with full therapy schedule, led to 47% reduction in the amount of drug dosage. Therefore, the optimum therapeutic schedule provides the high quality treatment, low amount of drug dosage and low side effects of antiretroviral drugs, simultaneously.

The CA model that was developed for AIDS in this study, involves the biological facts of this disease and it could depict all three patterns of progression of HIV infection properly. Most importantly, the effects of applying different antiretroviral drugs were also considered in this model. Therefore, as an application of the model developed in this study can be said, by using this model, physicians can evaluate the effects of applying different therapeutic methods on virtual patients through model, without requiring costly clinical trials on real patients. Also, by using Genetic Algorithms (GA), the best therapeutic schedule that provides the high quality treatment with low amount of drug dosage, can be obtained through this model. Therefore, results obtained from this model can be useful for physicians. However, this model is still developing. Finding newer and more accurate rules, more accurate estimating of model parameters, adding new states to show different aspects of disease, investigating more varied therapeutic methods and adjusting the drug dosage in the therapeutic schedule can be considered as future works.

Cite this paper

Golpayegani, G.N., Jafari, A.H. and Dabanloo, N.J. (2017) Providing a Therapeutic Scheduling for HIV Infected Individuals with Genetic Algorithms Using a Cellular Automata Model of HIV Infection in the Peripheral Blood Stream. J. Biomedical Science and Engineering, 10, 77-106. https://doi.org/10.4236/jbise.2017.103008

References

Abbreviation Note List

HIV: Human Immunodeficiency Virus

AIDS: Acquired Immunodeficiency Syndrome

RTI: Reverse Transcriptase Inhibitors

PI: Protease Inhibitors

cART: Combined Antiretroviral Therapy

CA: Cellular Automata

ODE: Ordinary Differential Equations

PDE: Partial Differential Equations

AB: Anti Bodies

GA: Genetic Algorithms