STRATEGIES OF OPTIMAL CONTROL FOR HIV SPREADS PREVENTION WITH HEALTH CAMPAIGN

In the present paper, we discuss an HIV (Human Immunodeficiency Virus) transmission model with health information campaign about HIV control policies. The reason behind the conception of this model is the idea to divide the human population that gains awareness of HIV, due to this campaign. We assume that HIV will not infect people who are aware about the dangers of HIV. We analyze the existence and local stability of the equilibrium points. We found that the disease-free equilibrium point will be locally asymptotically stable (LAS) if the basic reproduction number (R0) is less than one, and unstable otherwise. A forward bifurcation of the system is shown numerically, depending on the intervention parameters. Some numerical simulations for the autonomous system are given to see the evolution of the system, with respect to some scenarios that might appear in the field. To accommodate the limitation of budget issue for implementation in the field, the model is reconstructed as an optimal control problem with two control variables. Numerical simulations for optimal control problems are presented for five different scenarios. Numerical simulation results suggest that controlling strategies by providing health campaigns is better, if they precede an endemic prevention strategy than endemic reduction, since the cost needed for endemic reduction is five times higher, compared with the endemic prevention. Optimal intervention should also note the value of R0. Larger control levels are needed when R0 > 1 if compared with when R0 < 1. The results of numerical simulation also show that the lower the cost of the health campaign, the more health campaigns can be provided. Based on the numerical simulation calculations, the optimal intervention of health campaign can raise public awareness about HIV, in order to reduce the number of HIV-infected individuals.


INTRODUCTION
Human Immunodeficiency Virus (HIV) is a virus that attacks the human immune system. As long as the virus impairs immune cell functions, especially CD4 cells, the body becomes immunodeficient. The HIV infection further develops into Acquired Immunodeficiency Syndrome (AIDS). Individuals with AIDS experience severe immune system damage and suffer from opportunistic infections. HIV can be transmitted through the exchange of body fluids, such as blood, breast milk, semen, and vaginal secretions. HIV is spread mainly through unprotected sexual contact with HIV-infected individuals and sharing HIV-infected syringes. HIV can also be inherited vertically from an HIV-infected mother to her child. This transmission can occur during pregnancy, the birth process or breastfeeding [1].
In Indonesia, some government efforts to control HIV/AIDS include expanding access to CD4 testing and viral load, increasing coverage of antiretroviral treatment and improving the quality of healthcare facilities [2]. The government also conducts health campaigns to provide education, information and communication to adolescents and adults. The campaigns aim to enhance their knowledge concerning HIV and adopt positive behavior to prevent HIV. One of the programs is "Aku Bangga, Aku Tahu" (I'm Proud, I Know)-an HIV prevention campaign aimed at young people aged 15-24 years [3]. There are many methods to understand how HIV spreads; one among them is through mathematical models. Some mathematical models have already been introduced to understand how HIV spreads in the human population, such as by considering ART as prevention and treatment, as proposed by Huo et al. [4], Aldila D. and Maimunah [5], and Silva and Tores [6]. On the other hand, Naresh, Tripathi and Sharma [7] modeled the spread of HIV in the population with immigration from HIV-infected individuals. Furthermore, Giamberardino et al. [8] discussed the HIV prevalence model with three controls, such as information campaign, test campaign, and an HIV/AIDS therapy action. Dubey Preeti et al. [9] introduced a model of HIV dynamics in the body.
In this manuscript, a mathematical model of HIV involving informative campaign about the danger of HIV will be constructed. The model is based on the deterministic model in [8] by adding a specific compartment for susceptible and infected humans, who may be aware or unaware of HIV. Let the human population be divided into six sub-populations; let them be called as unaware of HIV susceptible individuals (S 1 ), aware of HIV susceptible individuals (S 2 ), unaware of HIV-infected individuals in the acute stage (I 1 ), aware of HIV-infected individuals in the acute stage (I 2 ), infected individuals in the chronic stage (P) and individuals with AIDS (A). The total population at time t, denoted by N(t), is given by N(t) = S 1 (t) + S 2 (t) + I 1 (t) + I 2 (t) + P(t) + A(t).
Newborns will be entered into the unaware susceptible group (S 1 ). There are health campaigns (u 1 and u 2 ) conducted by the government to control the spread and increase public awareness of HIV/AIDS. Let us define that u 1 is the individuals transition rate from unaware to aware, while u 2 is the opposite. It is assumed that there is a consideration of individual consciousness due to health campaigns, for example, with electronic media campaigns. If the content of the campaign is focused on providing basic knowledge about HIV and its prevention, then the number of individuals aware of the dangers of HIV will increase. In our model, this kind of intervention denoted by the increase of the transition rate from the unaware to the aware compartment (u 1 ). Conversely, if the content of the campaign is focused on regular appeals to maintain healthy lifestyle habits as a preventive measure of HIV transmission, people are expected to remain aware of HIV. So it is expected that individuals transition rate from aware to unaware does not increase. Examples of implementations in the field include direct campaign in AIDS rehabilitation groups/organization, free distribution of condoms [11], etc. Increasing u 1 and reducing u 2 is the best way to control the HIV spread. Unfortunately, this means of intervention comes at a high cost. Therefore, in this article, both parameters will be treated as time-dependent variables, u 1 (t) and u 2 (t). Further explanation about the optimal control problem of these parameters will be discussed in section 3.
In this article, although there is a possibility of vertical transmission of HIV to newborn [12], it is assumed that HIV is only transmitted through sexual contact between unaware susceptible individuals and unaware infected individuals in the acute stage. The unaware susceptible individuals will be infected through sexual contact with unaware infected individuals in the acute stage at a rate of transmission β . Infected individuals in the chronic stage and individuals with AIDS are assumed to not transmit HIV because they are already aware of their status of the disease or because of their health conditions that do not permit to make sexual contact.  [10], it is assumed in this article that infected individuals in the acute stage will move to infected individuals in the chronic stage at a rate δ . In the other hand, infected individuals in the chronic stage will suffer AIDS at a rate γ and individuals with AIDS will die caused by AIDS at a rate of α. With these assumptions, the model is given as follows, together with the parameters described in Table 1 and the epidemiological scheme of the system (1) in Figure 1.

Parameters Description
Unit Transition rate from the unaware to the aware compartment caused by medical campaign  The paper is organized as follows: the state of the art of the model and model construction are carefully detailed in this section and followed by mathematical analysis about the existence and local stability of equilibrium points in section 2. The construction of basic reproduction number is also given in section 2. Mathematical results from the analysis of the basic reproduction number in section 2 are then used in section 3 to identify the optimal parameters as an optimal control problem. Numerical experiments on the autonomous model when control parameters are constant or varying with time are given in section 4. Finally the discussion and conclusion are presented in section 5.

ANALYSIS OF THE AUTONOMOUS SYSTEM
Existence of the equilibrium points. To analyze the behaviour of the model in system (1), let us consider the control variables are constant (u 1 (t) = u 1 , u 2 (t) = u 2 ), such that model (1) now reads as: Taking the right hand side of model (2), system (2) has two equilibrium points. First equilibrium is the HIV-free equilibrium point, which is given by: With the HIV-free equilibrium in hand, we are ready to calculate the basic reproduction number of model (2).
Basic reproduction number. The Basic reproduction number (R 0 ) is defined as the expected number of secondary infections caused by one primary infected individual during a period of infection in a population of all susceptible individuals [13,14]. Some methods can be used to construct R 0 , such as with the Next-Generation matrix approach [14] or with graph theory [15]. In this article, the Next-Generation Matrix is the approach we choose to find the R 0 of system (2). First, let's construct the infective sub-system of (2), which only involves the I 1 (t), I 2 (t), P(t) and A(t), i.e., Next, using the small domain of next-generation matrix, then the next-generation matrix of system (2) is given by : .
Therefore, the basic reproduction number R 0 of system (2) as the spectral radius of NGM is given by: .
Further discussion about R 0 will be given in the later part of this section.
The next equilibrium is the endemic HIV-equilibrium point, where all susceptible and infected populations coexist. This equilibrium is given by: , , From the expression above, it can be seen that all compartments in Ω 2 will be positive if R 0 > 1. These results are stated in the following theorem. Theorem 2.1. Model (2) has two equilibrium points, i.e., the HIV-free equilibrium (3), which exists without any constraint and the endemic-HIV-equilibrium point (6), which will exist if R 0 > 1, where R 0 is the basic reproduction number of model (2) given in equation (5).
Next, we will analyze how R 0 changes with respect to the change in other parameters. Since we know that R 0 will increase when β increases linearly since ∂ R 0 ∂ β does not depend on β . On the other hand, since we know that R 0 will be suppressed non linearly when u 1 increases. Please note that when u 1 is large enough, u 1 >> 0, the effect of u 1 on suppressing R 0 is no longer significant. Next, since then the increasing number of people who return to the unaware sub-population will increase the basic reproduction number.
In analysis we have previously conducted, R 0 determines the existence and local stability criteria of the equilibrium points Ω 1 and Ω 2 . Therefore, next we will show how the health campaign interventions u 1 and u 2 determine the magnitude of R 0 , as shown in Figure 2. The analysis of Figure 2 is given as follows. In region 1, u 1 ∈ [0, 0.033), the R 0 will always be greater than one, regardless of the magnitude of health campaign intensity. Therefore, the intervention of the health campaign cannot effectively reduce the spread of HIV to the HIVfree equilibrium point in region 1. On the other hand, in region 3, when u 1 ∈ [0.824, 1], the R 0 < 1 will always be smaller than one for all magnitude of u 1 and u 2 . This means that if u 1 > 0.0824, then the disease will always be extinct in the population. Interesting discussions appear in regions 2a and 2b, where the combination of u 1 and u 2 will determine the R 0 , i.e., whether it is greater or less than one. To achieve the persistence of HIV in the field (R 0 < 1), for specific values of u 1 , u 2 should fulfill the condition or in the numerical example using the same parameters for Figure 2 except u 1 and u 2 , previous expression can be written as u 2 > −0.22 + 0.538u 1 + 0.641 2.016u 2 1 + 0.097. Therefore, if the government can conduct a health campaign u 1 in the region of (0.033, 0.824), then they have to really consider the rate of u 2 to achieve the persistence of HIV.
Local stability of equilibrium points. Next, the local stability of HIV-free equilibrium point is investigated analytically, while the local stability for endemic-HIV-equilibrium is investigated numerically using specific values of R 0 . The local stability of the equilibrium point is determined by the eigenvalues of the Jacobian matrix at the corresponding equilibrium point of the system (2). The general form of the Jacobian matrix of model (2) can be written by: where N c = S 1 + S 2 + I 1 + I 2 and To analyze the local stability of the HIV-free equilibrium point, we evaluate J in HIV-free equilibrium point which yield: The HIV-free equilibrium point will be stable if all the resulting eigenvalues from the above matrix are negative. We have four explicit eigenvalues, all of which are negative, i.e., −µ, −(µ + u 1 + u 2 ), −(α + µ), and −(µ + γ). The other two eigenvalues come from the solution of the second order characteristic polynomial given by : It can be seen that the polynomial in (10) will only have negative eigenvalues, if . This result is given in the following theorem.
Due to the complexity of the model and the form of endemic-HIV-equilibrium (Ω 2 ) of system (2), the local stability criteria of this equilibrium will be investigated by choosing a specific set of parameters, such that R 0 > 1, i.e., Λ = 200, µ = 0.02, β = 0.6, δ = 0.1, γ = 0.5, α = 0.1, u 1 = 0.1 and u 2 = 0.05.With this data set, we have that R 0 > 1 which makes the HIV-free equilibrium point unstable (Theorem (2.1)), and the endemic equilibrium point exist (Theorem (2.2)). Using this set of parameters, we have the endemic equilibrium point given by Ω 2 = (1482, 2118, 672, 395, 205, 101), and when we evaluate it in the Jacobian matrix J , the result gives us six eigenvalues, all of which are negative. Therefore, we conclude that there is a set of parameters, such that R 0 > 1 that makes the endemic-HIV-equilibrium point Ω 2 locally stable. A forward bifurcation diagram of the equilibrium points depending on u 1 is given in Figure 3. It can be seen that the larger the magnitude of u 1 in the set of u 1 ∈ [0, 0.096], the endemic equilibrium point Ω 2 , is stable and the population size of I 1 , P and A is decreasing and finally reaches 0 when u 1 tends to 0.096. On the other hand, I 2 in the endemic equilibrium is initially increasing, but later decreases to 0, when u 1 → 0.096 is caused by the high intensity of government campaigns to educate people about HIV. When u > 0.096, the R 0 is decreasing to less than one. According to the Theorem (2.1) and Theorem (2.2), the HIV-endemic equilibrium no longer exist, while the HIV-free equilibrium becomes stable.

OPTIMAL CONTROL CHARACTERIZATION
To investigate the optimal application of the health campaign needed to eliminate the spread of HIV, we would like to minimize the number of infected population using as low as possible control interventions. Therefore, let us consider the following objective functional subject to model in system (1) where T is the final time of the simulation. Let ω i for i = 1, 2, 3, 4 represent the weight parameters for the infected population, while ϕ 1 and ϕ 2 represent the wight parameters for the control variables. Please note that ω 1 I 1 (t), ω 2 I 2 (t), ω 3 P(t) and ω 4 A(t), represent the related cost as a consequence of the high number of infected populations at time t, for example for hospitalization cost. On the other hand, ϕ 1 u 1 (t) and ϕ 2 u 2 (t) represent the cost related to the effort of the government to implement the health campaign in order to educate people to become aware about the spread of HIV. To guarantee the balance of the cost function, the condition of ω i X i ≈ ϕ j U j should be reached, where X i is the infected compartment, while U j is the control variable. In this article, the use of quadratic form in the cost function J is chosen, since there is no linear correlation between the cost and effect of intervention. Please see [19,21,22] for further example of optimal control in the epidemiological model, using a quadratic form on their cost function.
We aim to minimize J with the constraint given by with N(t) = S 1 (t) + S 2 (t) + I 1 (t) + I 2 (t) + P(t) + A(t), and the initial conditions given as follows while also minimizing the control variables u 1 (t) and u 2 (t). We seek an optimal control u * 1 , u * 2 , such that where A is the admissible control depending on the lower (u min i ) and upper (u max i ) bounds for each control variable. Please note that u min i ≥ 0 and u max i < ∞. To guarantee the existence of the optimal control, it is directly followed by standard results of optimal control theory [23].
To derive the solution of our optimal control problem, we use the Pontryangin's Maximum Principle and also to derive the necessary conditions [24]. The Hamiltonian H H H is defined as where λ λ λ k (t) for k = 1, 2, ..., 6 is the adjoint variables .
Proof. First, to get the adjoint system (15) with the terminal condition λ λ λ i (T ) = 0 for i = 1, 2, 3, 4, 5, 6. To obtain the optimality condition (14), we will also differentiate the Hamiltonian in equation (13) with respect to control variables u 1 and u 2 , which gives ∂ H ∂ u 1 (t) = −I 1 (t)λ 3 (t) + I 1 (t)λ 4 (t) − S 1 (t)λ 1 (t) + S 1 (t)λ 2 (t) + 2 u 1 (t)ϕ 1 = 0, and set these equations equal to zero. Solving them with respect to each control variable, we obtain u * 1 (t) = To determine an acceptable control variable value based on the needs and ability in field applications (lower and upper bounds), the optimal control variables now arê with u max(i) and u min(i) for i = 1, 2 being the upper and lower bounds for each control variable, respectively.
In the next section, numerical simulation for the autonomous model in system (2) and the optimal control problem explained in Theorem (3.1) will be conducted.     Table 2 are given in intervals. This approach raises the question of how sensitive those parameters are in determining the results of the system (2). This is important to give a better understanding of the results of the optimal control problem in the next subsection. Thus, we have a sensitivity contour of parameters with respect to R 0 , before we gave the simulation of the autonomous model. In Figure 4 we can see how different values of u 1 = 0.01 k, k = 1, 2, ..., 30 affect the change of the dynamic of model (2). On the other hand, we can see how different values of u 2 = 0.01 k, k = 1, 2, ..., 30 affect the change of the dynamic of model (2) in Figure 5. Increasing the value of u 1 for u 2 as a constant, will increase total susceptible population (S 1 + S 2 ) and suppress the total infected population (I 1 + I 2 + P + A).

NUMERICAL
On the other hand, for u 1 as a constant and u 2 increasing, it will increase the total infected population and decreasing the total susceptible population. The reason behind this is that larger values of u 2 make more people unaware of the risk of HIV, making them more vulnerable to infection by I 1 , I 2 , P or A individuals.

4.2.
Simulation of the optimal control problem. In this section, we examine the optimal control problem of system (1) with the cost function given in (11). The optimal control set is obtained by solving the optimality system with a forward-backward scheme [16]. Please see [17,18,19,20] for further examples of implementation of this method in epidemiological models. The rough algorithm is given in flowchart in Figure 6.
To investigate the optimal control simulation, the numerical experiments will be conducted into four different scenario, i.e : There are two different initial conditions used to conduct numerical simulations in this section.
The first initial condition is when the number of infected individuals already relatively high, which named as the "endemic reduction" scenario. The second initial condition is the "endemic prevention" scenario, where the number of infected individuals is still relatively small in t = 0. Therefore, we have (S 1 , S 2 , I 1 , I 2 , P, A) = (8000, 0, 2000, 0, 0, 0) for the endemic reduction scenario, and (S 1 , S 2 , I 1 , I 2 , P, A) = (9800, 0, 200, 0, 0, 0) for the endemic prevention scenario.
Recapitulation of the numerical experiments related to the cost functions and equilibrium points are given in Table 3, before the control variables are implemented in Table 4

4.2.1.
Scenario I: Different initial R 0 . The simulation in this chapter is provided to illustrate how environmental conditions (infectious rate, recovery rate and other factor that described in R 0 ) affect the dynamics of control variables needed to suppress the spread of HIV. The simulation is given for two different R 0 conditions at t = 0 when u 1 = u 2 = 0, but with the same initial value for each variable in system (1), i.e R 0 > 1 and R 0 < 1. According to the analytical results in Theorem (2.1) and (2.2), system (1) will tend to HIV-endemic equilibrium when R 0 > 1 and goes to HIV-free equilibrium when R 0 < 1. For this simulation, we use the set of parameters in Table 2 except β = 0.5 such that R 0 = 1.19 > 1 and β = 0.3 such that R 0 = 0.71 < 1. We use the "endemic reduction" scenario for the initial value in this section.
We also bound the value of control variables, i.e u 1 ∈ [0.01, 0.5] and u 2 ∈ [0.01, 0.5]. The result is given in Table 5, while the details for each compartment before and after implementation of controls can be seen in Table 3  Based on Table 5, it can be seen that the cost function, when R 0 > 1, is almost five times higher than when R 0 < 1. This result is not surprising, since the system will tend to endemic equilibrium when R 0 > 1, i.e., (S 1 , S 2 , I 2 , I 2 , P, A) = (2009, 2, 380, 0, 293, 143). Therefore, more intervention is needed to control the spread of HIV. After t = 40, the number of infected individuals could be suppressed to 1 when R 0 > 1 and to 2 when R 0 < 1. The dynamics of the total susceptible population, the total infected population, and control variables for each case of scenario I could be seen in Figure 7 .
In Figure 7, it can be seen that in both cases, the behavior of control variables is almost the same. u 1 (t) behaves monotonically decreasing for all t ∈ [0, 40] while u 2 (t) is increasing, in the beginning, to suppress the number of infected individuals and then decreasing when the dynamics of infected individuals show a decreasing trend. Even though the magnitude of controls when R 0 < 1 is slightly higher than when R 0 > 1, the cost function is higher when R 0 > 1. This is because the number of infected individuals when R 0 > 1 that should be reduced is much higher than when R 0 < 1, since the number of infected individuals determines the cost function. It can be seen from Table 6 that the cost function for endemic reduction is relatively large, compared to the endemic prevention case. This is simply because more intervention is needed to control the number of infected individuals in an endemic reduction case. The dynamics of the total susceptible population, total infected population, and control variables for scenario II are given in Figure 8. Although the dynamics of humans for both cases has no big difference, the extreme difference is shown in the dynamics of control variables. In the endemic reduction case, the dynamic of u 1 follows the dynamics of infected individuals. A huge amount of u 1 is needed when the number of infected individuals is increasing. Different from the endemic reduction case, the behavior of control variables in the endemic prevention case is almost always constant all the time. The value of u 1 ∈ [0, 40] is the lower bound of u 1 , which means that no major intervention is needed, which gave a small magnitude of the cost function. The dynamics of control and human population for this scenario can be seen in Figure 9,

Scenario
where the numerical result is given in Table 7. A large value of ω i will increase the cost function. It can also be seen in Figure 9(b) that since the cost function for hospitalization is four times larger than in Figure 9(a), then the solution will be emphasized in the control function.
Therefore, we can see that the dynamics of the control variable in Figure 9(b) is much larger than in Figure 9(a) to avoid the larger cost caused by hospitalization. Please note that the total number of infected humans in Figure 9(b) already tends to 0, when t > 20 years, while in  Table 8. The huge value of weighted parameters for control variables will intervene with controls not in high numbers (Figure 10(b)) if compared with a situation when the weight which is dominated by the cost for controls. On the other hand, the cost related to hospitalization is more dominant, when ϕ i is relatively large, i.e., 130.4865. This is because the preferred solution is to spend all the resources on hospitalization, rather than on the expense of controls.

CONCLUSION
Mathematical modelling is an important tool to understand how a disease spreads and how many factors effect it. In this article, we have developed a mathematical model about the spread of HIV in the population, where medical campaigns are required to encourage the human population to become more aware about HIV. The model was constructed as a system of six dimensional ordinary differential equations with two control variables. This paper aims to find an optimal dynamics of control variables (medical campaign) to reduce the number of infected individuals to as small as possible. This task is designed as an optimal control problem.
Mathematical analyses concerning the equilibrium points and their local stability, along with the related basic reproduction numbers have been conducted analytically and numerically. We find that the HIV-free equilibrium point is locally stable when the basic reproduction number is smaller than one, and unstable otherwise. The instability of HIV-free equilibrium points then make the existence of the HIV-endemic equilibrium exist and locally stable. From analysis of the basic reproduction number with respect to the medical campaign parameters, we find that these controls are likely success to control the spread of HIV significantly.
The optimal control problem is solved with the Pontryagin Maximum/Minimum Principle (PMP). Four different numerical simulation scenarios are conducted to describe a possible situation that might occur in the field. We find that conducting the intervention in the early stage of the endemic (endemic prevention scenario) can reduce the intervention cost significantly, rather than waiting for the endemic to already occur.
In this present study, we consider the medical campaign as the only intervention to control the spread of HIV. There are many more interventions that can be considered and modelled with mathematics, such as controlling the misuse of needles, mass campaign about digress sex behaviour, and many more. Also, many other factors should be considered to make the model become more realistic, such as the vertical transmission of HIV, impact of blood transfusion containing HIV, etc. Thus, we can consider more detailed models to accommodate these issues.