Optimal Control of an HIV/AIDS Epidemic Model with Infective Immigration and Behavioral Change

Mamatjan Mastahun, Xamxinur Abdurahman^*

College of Mathematics and System Sciences, Xinjiang University, Urumqi, China.

DOI: 10.4236/am.2017.81008                

Abstract

In order to find out the effect of human (sexual) behavior change and immigration in spreading the HIV/AIDS, a deterministic model of HIV/AIDS with infective immigration is formulated. First, basic properties of the model, including non-negativity and boundedness of the solutions, existence of the endemic equilibrium and the basic reproduction number, R₀ are analyzed. The geometrical approach is used to obtain the global asymptotic stability of endemic equilibrium. Then the basic model is extended to include several control efforts aimed at reducing infection and changing behavior. Pontryagin’s maximum principle is used to derive the optimality system and solve the system numerically. Our numerical findings are illustrated through simulations using MATLAB, which shows reliability of our model from the practical point of view.

Keywords

HIV/AIDS, Infective Immigration, Global Stability, Optimal Control

1. Introduction

Mathematical models used extensively to study the dynamics of epidemics both from the cellular level to the population level by many researchers [1] - [6] . Early models concerned immigration of infective individuals are studied by many researchers [7] - [12] . In terms of epidemiology, the wild prevalence of HIV/AIDS has often been associated with the the movement of people, especially, the movement of infected persons. For example, in China, internal migrants exceeds 260 million [13] and the size of this population is expected to continue to grow. The whole migration process offers ready environment to transport the HIV/AIDS disease from one place to another and vulnerability of migrants may end up with certain HIV/AIDS related risky behaviors. Unlike traditional infectious disease, the process of transmission of HIV/AIDS mostly requires the exchange of the body fluids, this often associated with the decision making process in which the participants can either choose to quit or to proceed. Spontaneously the human sexual behavior comes to the equation. All these aspects have featured a large number of new problems.

To cope with these problems, Megan Coffee et al. [10] studied the dynamics of HIV/AIDS with infective immigration both clinically and mathematically. A nonlinear mathematical model consisted of two stages of infection before full-blown AIDS with constant inflow of HIV infectives are developed by Ram Naresh et al. [8] . But they did not consider direct inflow of pre-AIDS individuals, also not consider HIV infectives who received treatment. Agraj Tripathi et al. [9] proposed a HIV/AIDS model with infective immigrants and time delay. But they have not take the protective measures and gradual behavioral change into consideration. In the absence of an effective and affordable vaccine and due to the non-curative nature of current antiretroviral therapy, behavioral and psychosocial prevention with the goal of limiting risky sexual behaviors remain central to the efforts to decrease the sexual transmission of HIV [4] . Swarnali Sharma and G. P. Samanta [2] also analyzed six compartmental HIV/AIDS model which is closely influential to our set of work.

The optimal control theory has been applied to quite a few HIV models (see [1] [4] [6] ). Baba Seidu and Oluwole D. Makinde [1] stress their emphases on the optimally reducing the spread of the disease and the increasing productivity in workplace and show the positive effect of applying a multifaceted approach in the fight against HIV/AIDS. Tunde T. Yusuf and Francis Benyah [4] formulated an optimality system for controlling HIV/AIDS. They consider the change in risky sexual habits and antiretroviral (ARV) therapy as control measures. Their results show that if more and more susceptible individuals practise safe sex, then we can ease the spread of the disease remarkably. In the light of literatures [2] [4] [7] and references therein we formulated an optimal control problem that considers behavioral change and screening as major control strategies. The reviews about the non-mathematical studies of transmission disease and the comments of recent mathematical works about HIV/AIDS have shaped the context of this paper.

The paper organized as follows: In section 2, we have developed our model and the non-negativity and the boundedness of the solutions are shown as a basic property of the system. Also, we have discussed existence of endemic equilibrium. We derive basic reproduction number $R_0$ , when there is no influxion of infective. The analysis of GAS of the endemic equilibrium is given in Section 3, while our modification of the basic model into an optimal control problem presented in Section 4. In the section 5, we give results of our numerical simulations.

2. The Mathematical Model

We divided the sexually active population $N\left(t\right)$ into six compartments, namely, susceptible individuals $S\left(t\right)$ , infected individuals who unaware of their HIV/AIDS infection $I_1\left(t\right)$ , infected individuals who aware of their HIV/AIDS infection $I_2\left(t\right)$ , infected individuals who receive treatment $T\left(t\right)$ , full-blown AIDS group $A\left(t\right)$ and removed class $R\left(t\right)$ at any time $t$ . The population dynamics is given by the following set of ordinary differential equations:

$(\begin{array}{c}\stackrel{\dot{}}{S}\left(t\right)=\left(1-{\epsilon }_1-{\epsilon }_2\right)Q_0-\left({\beta }_1{I}_1+{\beta }_2{I}_2\right)S-\left(d+\theta \right)S,\\ {\stackrel{\dot{}}{I}}_1\left(t\right)=\delta {\epsilon }_1{Q}_0+\left({\beta }_1{I}_1+{\beta }_2{I}_2\right)S-\left(k_1+k_2+d\right)I_1,\\ {\stackrel{\dot{}}{I}}_2\left(t\right)=\left(1-\delta \right){\epsilon }_1{Q}_0+k_1{I}_1-\left(k_3+{\mu }_1+d\right)I_2,\\ \stackrel{\dot{}}{T}\left(t\right)={\mu }_1{I}_2+{\mu }_{2}A-\left(d+{\sigma }_2\right)T,\\ \stackrel{\dot{}}{A}\left(t\right)=k_2{I}_1+k_3{I}_2-\left({\mu }_2+d+{\sigma }_1\right)A,\\ \stackrel{\dot{}}{R}\left(t\right)={\epsilon }_2{Q}_0+\theta S-dR.\end{array}$ (1)

with initial conditions

$S\left(0\right)>0,I_1\left(0\right)>0,I_2\left(0\right)>0,T\left(0\right)\ge 0,A\left(0\right)\ge 0,R\left(0\right)\ge 0.$ (2)

There are a brief description of the model parameters:

Q₀: Total number of newly recruited individuals by birth (who come of age) and by immigration;

ε₁: The proportion of infected individuals in the recruited population;

ε₂: The proportion of susceptible individuals who changed their sexual habits in the recruited population;

δ: The proportion of unaware infected individuals and $\left(1-\delta \right)$ the proportion of aware infected individuals in the recruited population;

β₁: The horizontal transmission rate for contact with the $I_1$ class;

β₂: The horizontal transmission rate for contact with the $I_2$ class;

d: The natural death rate of population;

θ: Proportion of susceptible individuals who changed their sexual habits;

κ₁: Proportion of unaware infected individuals who are screened;

κ₂: Progression rate of unaware infected individuals to the full-blown AIDS group;

κ₃: Progression rate of aware infected individuals to the full-blown AIDS group;

μ₁: Proportion of the $I_2$ class receiving treatment;

μ₂: Proportion of the $A$ class receiving treatment;

σ₁: Disease-induced death rate for full-blown AIDS individuals;

σ₂: Disease-induced death rate for the individuals who receive the treatment.

The model is formulated based on following assumptions:

・ Since our purpose in this model is to see what effect the human behavior (including movement, sexual habits) can play in the dynamics of HIV/AIDS disease, we avoid to consider detailed clinical stages of HIV/AIDS infection, instead we classed the population in two ways, uninfected and infected group. Uninfected group divided into two different compartments according to their behavior towards safe sex. Infected individuals divided four different compartments according to whether the infected individual aware of his/her HIV infection status, whether the infected individual received treatment and whether the infected individual has developed the last stage of the disease, the full blown AIDS.

・ Susceptible individuals are assumed to get infected by sexual contact with both aware and unaware infected individuals with different transmission rates. The assumption that aware infected individuals also take part in the transmission process is based on the fact that some aware infected individuals may practise low-efficiency safe sex measures (and a few of them may transmit the disease intentionally), and susceptible individuals may not be aware of the infected situation of his/her partner, which make them more vulnerable to the disease. So the new generated infected individuals by aware infective individuals are assumed to be not aware of his/her infection at first and go to the unaware infected individuals class.

・ The simplest conceptual framework based on homogeneous behavior gives us clear insights into how community based chemotherapy can influence epidemiological pattern and transmission success. Here the mixing of susceptibles with infectives is considered to be homogeneous and accordingly the incidence rate is assumed to be bilinear [2] [5] .

・ All new born are susceptible, i.e., in our model vertical transmission do not account for.

・ We assumed that individuals in the treatment class not only to receive the ART therapy, but also to be served with knowledge about the HIV/AIDS disease so that they were persuaded to avoid unsafe sexual behaviors. Full- blown AIDS individuals are assumed too ill to sexually active, So the suscep- tibles do not get infected through sexual contacts with individuals from these two groups.

・ Inclusion of compartment $R$ : It is true that an appreciable number of people are now changing their sexual habits sufficiently due to the awareness of the widespread nature of disease in society, the monumental deaths resulting from the disease, increasing knowledge of the agony and psychological trauma experienced by the infected individuals, and better enlightenment due to intense HIV/AIDS educational campaigns [4] .

2.1. Basic Properties of the Model

The model system (1) describes human population and therefore it is necessary to prove that all the variables $S\left(t\right)$ , $I_1\left(t\right)$ , $I_2\left(t\right)$ , $T\left(t\right)$ , $A\left(t\right)$ and $R\left(t\right)$ are non-negative for all time. Solutions of the model system (1) with positive initial conditions (2) remains positive for all time $t\ge 0$ and are bounded in $G$ , where

$\begin{array}{l}G=\{\left(S\left(t\right),I_1\left(t\right),I_2\left(t\right),T\left(t\right),A\left(t\right),R\left(t\right)\right)\subseteq {\mathcal{R}}_{+}^6,\frac{Q_0}{d+{\sigma }_1+{\sigma }_2}\\ <S\left(t\right)+I_1\left(t\right)+I_2\left(t\right)+T\left(t\right)+A\left(t\right)+R\left(t\right)\le \frac{Q_0}{d}\}\end{array}$

is defined based on biological considerations and positively invariant with respect to the model system (1).

Hence the theorem

Theorem 2.1 Every solution of the system (1) with initial conditions (2) exists in the interval $t\in \left[0,\infty \right)$ and $S\left(t\right)>0$ , $I_1\left(t\right)>0$ , $I_2\left(t\right)>0$ , $T\left(t\right)\ge 0$ ,

$A\left(t\right)\ge 0$ and $R\left(t\right)\ge 0$ , for all $t\ge 0$ . For the model system (1), the region $G$ is positively invariant and all solutions starting in $G$ approach, enter, or stay in $G$ .

Proof Since the right hand side of system (1) is completely continuous and locally Lipschitzian on C (space of continuous functions), the solution ( $S\left(t\right)$ , $I_1\left(t\right)$ , $I_2\left(t\right)$ , $T\left(t\right)$ , $A\left(t\right)$ , $R\left(t\right)$ ) of (1) with initial conditions (2) exists and is unique on [0; ξ), where $0<\xi \le +\infty$ . Under the given initial conditions (2), it is easy to prove that the components of solutions of the model system (1) are positive; if not, we assume a contradiction: that there exists a first time

$t_1=\inf \left\{t|S\left(t\right)=0or{I}_1\left(t\right)=0or{I}_2\left(t\right)=0orT\left(t\right)=0orA\left(t\right)=0orR\left(t\right)=0\right\}$

If $S\left(t_1\right)=0$ , $S\left(t\right)>0$ , $I_1\left(t\right)>0$ , $I_2\left(t\right)>0$ , $T\left(t\right)>0$ , $A\left(t\right)>0$ and $R\left(t\right)>0$ for $t\in \left(0;t_1\right)$ , then we have $S\text{'}\left(t_1\right)<0$ , but from the first equation of system (1) we have

$S\text{'}\left(t_1\right)=\left(1-{\epsilon }_1-{\epsilon }_2\right)Q_0>0$

which is a contradiction meaning that $S\left(t\right)$ remains positive; the other cases where $I_1\left(t_1\right)=0$ , $I_2\left(t_1\right)=0$ , $T\left(t_1\right)=0$ , $A\left(t_1\right)=0$ and $R\left(t_1\right)=0$ can also be discussed similarly as above. Thus in all cases $S\left(t\right)$ , $I_1\left(t\right)$ , $I_2\left(t\right)$ , $T\left(t\right)$ , $A\left(t\right)$ and $R\left(t\right)$ remain positive for all $t\ge 0$ . Since $N\left(t\right)\ge A\left(t\right)+T\left(t\right)$ , then

$Q_0-\left(d+{\sigma }_1+{\sigma }_2\right)N\le \stackrel{\dot{}}{N}\left(t\right)=Q_0-dN-{\sigma }_{1}A-{\sigma }_{2}T\le Q_0-dN$

which implies that $N\left(t\right)$ is bounded and all the solutions starting in G approach, enter or stay in G.

This completes the proof.

Since the variables $T$ , $A$ and $R$ of the system (1) do not appear in the first three equations of the system (1), in the subsequent analysis we only consider the following subsystem:

$(\begin{array}{l}\stackrel{\dot{}}{S}\left(t\right)=\left(1-{\epsilon }_1-{\epsilon }_2\right)Q_0-\left({\beta }_1{I}_1+{\beta }_2{I}_2\right)S-\left(d+\theta \right)S,\hfill \\ {\stackrel{\dot{}}{I}}_1\left(t\right)=\delta {\epsilon }_1{Q}_0+\left({\beta }_1{I}_1+{\beta }_2{I}_2\right)S-\left(k_1+k_2+d\right)I_1,\hfill \\ {\stackrel{\dot{}}{I}}_2\left(t\right)=\left(1-\delta \right){\epsilon }_1{Q}_0+k_1{I}_1-\left(k_3+{\mu }_1+d\right)I_2,\hfill \end{array}$ (3)

with initial conditions

$S\left(0\right)>0,I_1\left(0\right)>0,I_2\left(0\right)>0.$ (4)

2.2. Existence of Endemic Equilibrium

The system (3) does not exhibit a disease-free equilibrium due to direct inflow of population at a constant rate. However, there exists only one non-negative equilibrium point of the model (3), i.e., endemic equilibrium $E^{\ast }\left(S^{\ast },I_1^{\ast },I_2^{\ast }\right)$ , where $S^{\ast }$ , $I_1^{\ast }$ , $I_2^{\ast }$ are positive solutions of the following system of algebraic equations,

$\left(1-{\epsilon }_1-{\epsilon }_2\right)Q_0-\left({\beta }_1{I}_1^{\ast }+{\beta }_2{I}_2^{\ast }\right)S^{\ast }-\left(d+\theta \right)S^{\ast }=0,$ (5)

$\delta {\epsilon }_1{Q}_0+\left({\beta }_1{I}_1^{\ast }+{\beta }_2{I}_2^{\ast }\right)S^{\ast }-\left(k_1+k_2+d\right)I_1^{\ast }=0,$ (6)

$\left(1-\delta \right){\epsilon }_1{Q}_0+k_1{I}_1^{\ast }-\left(k_3+{\mu }_1+d\right)I_2^{\ast }=0.$ (7)

By solving Equations (5)-(7), we get

$S^{\ast }=\frac{\left(1-{\epsilon }_1-{\epsilon }_2\right)Q_0}{A{I}_1^{\ast }+B},$ (8)

$I_2^{\ast }=\frac{C}{{\beta }_2}+\frac{k_1{I}_1^{\ast }}{k_3+{\mu }_1+d},$ (9)

$I_1^{\ast }$ is a positive solution of following quadratic equation

$a{\left(I_1^{\ast }\right)}^2+b{I}_1^{\ast }+c=0,$ (10)

where

$A={\beta }_1+\frac{{\beta }_2{k}_1}{k_3+{\mu }_1+d},B=\frac{{\beta }_2\left(1-\delta \right){\epsilon }_1{Q}_0}{k_3+{\mu }_1+d}+d+\theta ,C=\frac{{\beta }_2\left(1-\delta \right){\epsilon }_1{Q}_0}{k_3+{\mu }_1+d},$

and

$\begin{array}{l}a=-\left(k_1+k_2+d\right)A,\\ b=\left[\delta {\epsilon }_1+\left(1-{\epsilon }_1-{\epsilon }_2\right)\right]Q_{0}A-\left(k_1+k_2+d\right)B,\\ c=\left[\delta {\epsilon }_{1}B+\left(1-{\epsilon }_1-{\epsilon }_2\right)C\right]Q_0.\end{array}$

Obviously $a$ is always negative and $c$ is always positive. By applying the Descartes’ rule of signs, one positive root of Equation (10) exists, whatever is the sign of $b$ , i.e., we always get a positive equilibrium of the system (3). Above discussions can be summarized as:

Theorem 2.2 The system (3) has a endemic equilibrium $E^{\ast }\left(S^{\ast },I_1^{\ast },I_2^{\ast }\right)$ , which exists for all parameter values.

Especially when ${\epsilon }_1=0$ , we can easily obtain basic reproduction number $R_0$ of the system (3) by using the next generation method as

$R_0=\frac{\left(1-{\epsilon }_2\right)Q_0}{d+\theta }\frac{{\beta }_1\left(k_3+{\mu }_1+d\right)+{\beta }_2{k}_1}{\left(k_1+k_2+d\right)\left(k_3+{\mu }_1+d\right)}$ .

3. Global Stability Analysis of the Endemic Equilibrium

In this section we shall discuss the global stability of the endemic equilibrium $E^{\ast }\left(S^{\ast },I_1^{\ast },I_2^{\ast }\right)$ .

Theorem 3.1 The endemic equilibrium $E^{\ast }\left(S^{\ast },I_1^{\ast },I_2^{\ast }\right)$ of the system (3) is

globally asymptotically stable, if $R_1\triangleq \frac{{\delta }_1{\epsilon }_{1}d+d}{2{\beta }_1\frac{Q_0}{d}+\theta }>1$ , where ${\delta }_1=\min \left\{\delta ,1-\delta \right\}$ .

Before start our proof, we first recall the following lemma by Li and Muldowney [14] .

Lemma If the system

$\frac{\text{d}x}{\text{d}t}=f\left(x\right),$

where $x\rightarrowtail f\left(x\right)\in {ℝ}^n$ , be a $C^1$ function for $x$ in an open set ${\Gamma }_1\subset {ℝ}^n$ such that

1) has a unique equilibrium $x^{\ast }$ in ${\Gamma }_1$ and

2) there exists a compact absorbing set $K\subset {\Gamma }_1$ ,

then the equilibrium $x^{\ast }$ is globally asymptotically stable provided that a

$\left(\begin{array}{c}2\\ n\end{array}\right)\times \left(\begin{array}{c}2\\ n\end{array}\right)$ matrix valued function $P\left(x\right)$ and a Lozinskii measure $\mathcal{L}$ of $B$ with respect to a vector norm $|\cdot |$ in ${ℝ}^N$ , $N=\left(\begin{array}{c}2\\ n\end{array}\right)$ exist such that

$\underset{t\to \infty }{{\displaystyle \lim \sup }}\underset{x\in K}{{\displaystyle \sup }}\left(\frac{1}{t}{\displaystyle {\int }_0^t}\mathcal{L}\left(B\left(x\left(s,x\right)\right)\right)\text{d}s\right)<0$

Here

$B=P_f{P}^{-1}+P{J}^{\left[2\right]}{P}^{-1},$

the matrix $P_f$ is obtained by replacing each entry $P_{ij}$ of $P$ by its derivative in the direction of $f$ and $J^{\left[2\right]}$ is the second additive compound matrix of the Jacobian matrix $J$ , i.e., $J\left(x\right)=Df\left(x\right)$ and

$\mathcal{L}\left(B\right)=\underset{h\to 0}{\lim }\frac{\left|I+hB\right|-1}{h}.$

Proof of Theorem 3.1 The system is uniformly persistent in $G$ (given in Theorem 2.1). The uniform persistence of system (3) in the bounded set $G$ is equivalent to the existence of a compact set $K\subset G$ that is absorbing for system (3). Therefore the system (3) satisfies the conditions $\left(i\right)$ and $\left(ii\right)$ of the pre- vious Lemma.

The Jacobian matrix of system (3) is given by

$J=\left(\begin{array}{ccc}-\left({\beta }_1{I}_1+{\beta }_2{I}_2+d+\theta \right)& -{\beta }_{1}S& -{\beta }_{2}S\\ {\beta }_1{I}_1+{\beta }_2{I}_2& {\beta }_{1}S-\left(k_1+k_2+d\right)& {\beta }_{2}S\\ 0& k_1& -\left(k_3+{\mu }_1+d\right)\end{array}\right)=\left(\begin{array}{ccc}{J}_{11}& J_{12}& J_{13}\\ J_{21}& J_{22}& J_{23}\\ J_{31}& J_{32}& J_{33}\end{array}\right)$

where $J_{ij}$ is the corresponding entry of the matrix $J$ .

The second additive compound matrix of $J$ is given by

$\begin{array}{c}{J}^{\left[2\right]}=\left(\begin{array}{ccc}{J}_{11}+J_{22}& J_{23}& -J_{13}\\ J_{32}& J_{11}+J_{22}& J_{12}\\ -J_{31}& J_{21}& J_{22}+J_{33}\end{array}\right)\\ =\left(\begin{array}{ccc}\begin{array}{l}{\beta }_{1}S-({\beta }_1{I}_1+{\beta }_2{I}_2+d+\theta \\ +k_1+k_2+d)\end{array}& {\beta }_{2}S& {\beta }_{2}S\\ k_1& \begin{array}{l}-({\beta }_1{I}_1+{\beta }_2{I}_2+d+\theta \\ +k_3+{\mu }_1+d)\end{array}& -{\beta }_{1}S\\ 0& {\beta }_1{I}_1+{\beta }_2{I}_2& \begin{array}{l}{\beta }_{1}S-(k_1+k_2+d\\ +k_3{\mu }_1+d)\end{array}\end{array}\right)\end{array}$

Consider the function

$P=P\left(S,I_1,I_2\right)=\text{diag}\left(1,\frac{I_1}{I_2},\frac{I_1}{I_2}\right)$

then

$\begin{array}{l}{P}^{-1}=\text{diag}\left(1,\frac{I_2}{I_1},\frac{I_2}{I_1}\right),\\ P_f=\text{diag}\left(0,\frac{{I^{\prime }}_1{I}_2-I_1{I^{\prime }}_2}{I_2^2},\frac{{I^{\prime }}_1{I}_2-I_1{I^{\prime }}_2}{I_2^2}\right)\end{array}$

Therefore

$P_f{P}^{-1}=\text{diag}\left(0,\frac{{I^{\prime }}_1}{I_1}-\frac{{I^{\prime }}_2}{I_2},\frac{{I^{\prime }}_1}{I_1}-\frac{{I^{\prime }}_2}{I_2}\right)$

and

$P{J}^{\left[2\right]}{P}^{-1}=\left(\begin{array}{ccc}\begin{array}{l}{\beta }_{1}S-({\beta }_1{I}_1+{\beta }_2{I}_2+d+\theta \\ +k_1+k_2+d)\end{array}& \frac{I_2}{I_1}{\beta }_{2}S& \frac{I_2}{I_1}{\beta }_{2}S\\ \frac{I_1}{I_2}{k}_1& \begin{array}{l}-({\beta }_1{I}_1+{\beta }_2{I}_2+d+\theta \\ +k_3+{\mu }_1+d)\end{array}& -{\beta }_{1}S\\ 0& {\beta }_1{I}_1+{\beta }_2{I}_2& \begin{array}{l}{\beta }_{1}S-(k_1+k_2+d\\ +k_3+{\mu }_1+d)\end{array}\end{array}\right).$

Therefore,

$B=P_f{P}^{-1}+P{J}^{\left[2\right]}{P}^{-1}=\left(\begin{array}{cc}{B}_{11}& B_{12}\\ B_{21}& B_{22}\end{array}\right),$

where

$\begin{array}{l}{B}_{11}={\beta }_{1}S-\left({\beta }_1{I}_1+{\beta }_2{I}_2+d+\theta +k_1+k_2+d\right),\\ B_{12}=\left(\begin{array}{cc}\frac{I_2}{I_1}{\beta }_{2}S& \frac{I_2}{I_1}{\beta }_{2}S\end{array}\right),\\ B_{21}=\left(\begin{array}{c}\frac{I_1}{I_2}{k}_1\\ 0\end{array}\right)\end{array}$

$B_{22}=\left(\begin{array}{cc}\begin{array}{l}\frac{{I^{\prime }}_1}{I_1}-\frac{{I^{\prime }}_2}{I_2}-({\beta }_1{I}_1+{\beta }_2{I}_2+d+\theta \\ +k_3+{\mu }_1+d)\end{array}& -{\beta }_{1}S\\ {\beta }_1{I}_1+{\beta }_2{I}_2& \begin{array}{l}\frac{{I^{\prime }}_1}{I_1}-\frac{{I^{\prime }}_2}{I_2}-(k_1+k_2+d\\ +k_3+{\mu }_1+d)\end{array}\end{array}\right)$

Let $\left(u,\upsilon ,\omega \right)$ be a vector in ${ℝ}^3$ and its norm $|\cdot |$ is defined as

$\left|\left(u,\upsilon ,\omega \right)\right|=\max \left\{\left|u\right|,\left|\upsilon \right|+\left|\omega \right|\right\}.$

Let $\mathcal{L}\left(\mathcal{B}\right)$ be the Lozinskii measure with respect to this norm. Then we choose

$\mathcal{L}\left(\mathcal{B}\right)\le \sup \left\{g_1,g_2\right\},$

where

$g_1={\mathcal{L}}_1\left(B_{11}\right)+\left|B_{12}\right|,g_2={\mathcal{L}}_1\left(B_{22}\right)+\left|B_{21}\right|,$

where $\left|B_{12}\right|$ , $\left|B_{21}\right|$ are matrix norms with respect to the $L^1$ vector norm and ${\mathcal{L}}_1$ denotes the the Lozinski measure with respect to the $L^1$ norm.

Therefore,

$\begin{array}{l}{\mathcal{L}}_1\left(B_{11}\right)={\beta }_{1}S-\left({\beta }_1{I}_1+{\beta }_2{I}_2+d+\theta +k_1+k_2+d\right),\\ \left|B_{12}\right|=\frac{I_2}{I_1}{\beta }_{2}S,\\ \left|B_{21}\right|=\frac{I_1}{I_2}{k}_1,\\ {\mathcal{L}}_1\left(B_{22}\right)=\max \left\{\frac{I_1^{\text{'}}}{I_1}-\frac{I_2^{\text{'}}}{I_2}-\left(d+\theta +k_3+{\mu }_1+d\right),\frac{I_1^{\text{'}}}{I_1}-\frac{I_2^{\text{'}}}{I_2}+2{\beta }_{1}S-\left(k_1+k_2+d+k_3+{\mu }_1+d\right)\right\}.\\ =\frac{I_1^{\text{'}}}{I_1}-\frac{I_2^{\text{'}}}{I_2}-\left(d+k_3+{\mu }_1+d\right)+\max \left\{-\theta ,2{\beta }_{1}S-\left(k_1+k_2\right)\right\}\end{array}$

So, we have

$\begin{array}{c}{g}_1={\mathcal{L}}_1\left(B_{11}\right)+\left|B_{12}\right|={\beta }_{1}S-\left({\beta }_1{I}_1+{\beta }_2{I}_2+d+\theta +k_1+k_2+d\right)+\frac{I_2}{I_1}{\beta }_{2}S\\ ={\beta }_{1}S+\frac{I_2}{I_1}{\beta }_{2}S+\left(k_1+k_2+d\right)-\left({\beta }_1{I}_1+{\beta }_2{I}_2+d+\theta \right)\\ =\frac{{I^{\prime }}_1}{I_1}-\frac{\delta {\epsilon }_1{Q}_0}{I_1}-\left({\beta }_1{I}_1+{\beta }_2{I}_2+d+\theta \right)\\ \le \frac{{I^{\prime }}_1}{I_1}-{\delta }_1{\epsilon }_{1}d-\left({\beta }_1+{\beta }_2\right)\frac{Q_0}{d+{\sigma }_1+{\sigma }_2}-d-\theta ,\\ g_2={\mathcal{L}}_1\left(B_{22}\right)+\left|B_{21}\right|=\frac{{I^{\prime }}_1}{I_1}-\frac{{I^{\prime }}_2}{I_2}-\left(d+k_3+{\mu }_1+d\right)+\max \left\{-\theta ,2{\beta }_{1}S-\left(k_1+k_2\right)\right\}+\frac{I_1}{I_2}{k}_1\\ =\frac{{I^{\prime }}_1}{I_1}-d+\max \left\{-\theta ,2{\beta }_{1}S-\left(k_1+k_2\right)\right\}-\frac{{I^{\prime }}_2}{I_2}-\left(k_3+{\mu }_1+d\right)+\frac{I_1}{I_2}{k}_1\\ =\frac{{I^{\prime }}_1}{I_1}-d+\max \left\{-\theta ,2{\beta }_{1}S-\left(k_1+k_2\right)\right\}-\frac{\left(1-\delta \right){\epsilon }_1{Q}_0}{I_2}\\ \le \frac{{I^{\prime }}_1}{I_1}-d+\max \left\{0,2{\beta }_{1}S+\theta -\left(k_1+k_2\right)\right\}-\left(1-\delta \right){\epsilon }_{1}d\\ \le \frac{{I^{\prime }}_1}{I_1}-d+2{\beta }_1\frac{Q_0}{d}+\theta -{\delta }_1{\epsilon }_{1}d\end{array}$

where ${\delta }_1=\min \left\{\delta ,1-\delta \right\}$ .

Furthermore, we obtain

$\begin{array}{c}\mathcal{L}\left(\mathcal{B}\right)\le \sup \left\{g_1,g_2\right\}\\ \le \sup \left\{\frac{{I^{\prime }}_1}{I_1}-{\delta }_1{\epsilon }_{1}d-\left({\beta }_1+{\beta }_2\right)\frac{Q_0}{d+{\sigma }_1+{\sigma }_2}-d-\theta ,\frac{{I^{\prime }}_1}{I_1}-d+2{\beta }_1\frac{Q_0}{d}+\theta -{\delta }_1{\epsilon }_{1}d\right\}\\ =\frac{{I^{\prime }}_1}{I_1}-{\delta }_1{\epsilon }_{1}d-d+\max \left\{-\left({\beta }_1+{\beta }_2\right)\frac{Q_0}{d+{\sigma }_1+{\sigma }_2}-d-\theta ,2{\beta }_1\frac{Q_0}{d}+\theta \right\}\\ =\frac{{I^{\prime }}_1}{I_1}-{\delta }_1{\epsilon }_{1}d-d+2{\beta }_1\frac{Q_0}{d}+\theta ,\end{array}$

then

$\begin{array}{c}\frac{1}{t}{\displaystyle {\int }_0^t}\mathcal{L}\left(B\right)\text{d}\zeta \le \frac{1}{t}{\displaystyle {\int }_0^t}\left(\frac{{I^{\prime }}_1}{I_1}-{\delta }_1{\epsilon }_{1}d-d+2{\beta }_1\frac{Q_0}{d}+\theta \right)\text{d}\zeta \\ =\frac{1}{t}\ln \frac{I_1}{I_1\left(0\right)}-{\delta }_1{\epsilon }_{1}d-d+2{\beta }_1\frac{Q_0}{d}+\theta .\end{array}$

Therefore

$\underset{t\to \infty }{{\displaystyle \lim \sup }}\underset{t\in K}{{\displaystyle \sup }}\left(\frac{1}{t}{\displaystyle {\int }_0^t}\mathcal{L}\left(B\left(x\left(s,x\right)\right)\right)\text{d}s\right)<\frac{1}{2{\beta }_1\frac{Q_0}{d}+\theta }\left(1-R_1\right)<0,\text{if}{R}_1>1$

Then based on Theorem 3.5 of [14] the endemic equilibrium $E^{\ast }$ of the system (3) is globally asymptotically stable. Hence the theorem.

4. Optimal Control Problem Formulation

4.1. Existence of an Optimal Control Pair

In this section we give following optimal control problem.

$(\begin{array}{l}\stackrel{\dot{}}{S}\left(t\right)=\left(1-{\epsilon }_1-{\epsilon }_2\right)Q_0-\left[{\beta }_1{I}_1+\left(1-u_1\right){\beta }_2{I}_2\right]S-\left(d+u_2\right)S,\hfill \\ {\stackrel{\dot{}}{I}}_1\left(t\right)=u_3{\epsilon }_1{Q}_0+\left[{\beta }_1{I}_1+\left(1-u_1\right){\beta }_2{I}_2\right]S-\left(u_4+k_2+d\right)I_1,\hfill \\ {\stackrel{\dot{}}{I}}_2\left(t\right)=\left(1-u_3\right){\epsilon }_1{Q}_0+u_4{I}_1-\left(k_3+{\mu }_1+d\right)I_2,\hfill \\ \stackrel{\dot{}}{T}\left(t\right)={\mu }_1{I}_2+{\mu }_{2}A-\left(d+{\sigma }_2\right)T,\hfill \\ \stackrel{\dot{}}{A}\left(t\right)=k_2{I}_1+k_3{I}_2-\left({\mu }_2+d+{\sigma }_1\right)A,\hfill \\ \stackrel{\dot{}}{R}\left(t\right)={\epsilon }_2{Q}_0+u_{2}S-dR.\hfill \end{array}$ (11)

satisfying initial conditions given in (2). We consider following optimal control parameters

1) $u_1$ is the control effort aimed at reducing the infection of susceptible individuals;

2) $u_2$ is the control effort aimed at persuading people to change or moderate their sexual behaviour;

3) $u_3$ is the control effort aimed at screening new arrivals;

4) $u_4$ is the control effort aimed at screening.

The objective functional [1] [4] is defined as

$J=\underset{u_1,u_2,u_3,u_4}{\min }\left({\displaystyle {\int }_0^{t_f}}\left[a_{1}S+a_2{I}_1+a_3{I}_2+a_{4}A+\frac{1}{2}\left({\omega }_1{u}_1^2+{\omega }_2{u}_2^2+{\omega }_3{u}_3^2+{\omega }_4{u}_4^2\right)\right]\text{d}t\right)$ (12)

the weight constants $a_1$ , $a_2$ , $a_3$ , $a_4$ , ${\omega }_1$ , ${\omega }_2$ , ${\omega }_3$ and ${\omega }_4$ are the relative weights and help to balance each term in the integrand so that any of the terms do not dominate. $t_f$ is the final time. Our aim is to minimize the objective functional or cost function $J\left(u_1,u_2,u_3,u_4\right)$ given in (15) so that the infective as well as the cost of implementing the control strategies can be minimized. So, we seek a set of optimal control $\left(u_1^{\ast },u_2^{\ast },u_3^{\ast },u_4^{\ast }\right)$ such that

$J\left(u_1^{\ast },u_2^{\ast },u_3^{\ast },u_4^{\ast }\right)=\min \left\{J\left(u_1,u_2,u_3,u_4\right):\left(u_1,u_2,u_3,u_4\right)\in U\right\},$

the admissible control set is given as

$U=\left\{\left(u_1,u_2,u_3,u_4\right)|u_i\text{isLebesguemeasurablewith}0\le u_i\le 1,\forall t\in \left[0,t_f\right]\right\}.$ (13)

Theorem 4.1 Given the objective function $J\left(u_1,u_2,u_3,u_4\right)$ as (15) with admissible control set $U$ , subject to the system (14) with initial conditions (2), then there exist an optimal control set $\left(u_1^{\ast },u_2^{\ast },u_3^{\ast },u_4^{\ast }\right)$ , such that (16), if the following conditions are satisfied:

1) The class of all initial conditions with the corresponding control functions in $U$ is non-empty;

2) The admissible control set $U$ is closed and convex;

3) Each RHS of the system (14) is continuous and is bounded above by sum of bounded control and state and can be written as a linear function of the control variables with coefficients dependent on time and state variables;

4) The integrand $L$ in equation $J$ is convex on $U$ and additionally satisfies $L\left(t,S,I_1,I_2,A,T,R\right)\ge c_1{\left|\left(u_1,u_2,u_3,u_4\right)\right|}^{\alpha }-c_2$ , where $c_1$ , $c_2>0$ and $\alpha >1$ .

Proof 1) We refer to Theorem 3.1 proposed by Picard-Lindelof in [15] . based on this theorem, if the solutions to the state equations are priori bounded and if the state equations are continuous and Lipschitz in the state variables, then there exists a unique solution, corresponding to the every admissible control set in $U$ . Using the fact that for all $L\left(t,S,I_1,I_2,A,T,R\right)\in G$ , all the model states are bounded below and above, the solutions to the state equations are bounded. In addition, the boundedness of the partial derivatives with respect to the state variables in the system can be directly shown, and this shows that the system is Lipschitz with respect to the state variables.

2) The control set $U$ is convex and closed by definition.

3) We observe that the integrand $L$ in our objective functional is convex, since it is quadratic in the controls.

4) $\begin{array}{c}L=a_{1}S+a_2{I}_1+a_3{I}_2+a_{4}A+\frac{1}{2}\left({\omega }_1{u}_1^2+{\omega }_2{u}_2^2+{\omega }_3{u}_3^2+{\omega }_4{u}_4^2\right)\\ \ge \frac{1}{2}\left({\omega }_1{u}_1^2+{\omega }_2{u}_2^2+{\omega }_3{u}_3^2+{\omega }_4{u}_4^2\right)\text{since}{\omega }_i,a_i>0,i=1,2,3,4\\ \ge \frac{1}{2}\left({\omega }_1{u}_1^2+{\omega }_2{u}_2^2+{\omega }_3{u}_3^2+{\omega }_4{u}_4^2\right)-{\omega }_1\text{since}{\omega }_1{u}_1^2-{\omega }_1\le 0\\ \ge \frac{1}{2}\min \left\{{\omega }_1,{\omega }_2,{\omega }_3,{\omega }_4\right\}\left(u_1^2+u_2^2+u_3^2+u_4^2\right)-{\omega }_1\\ \ge W{\Vert \left(u_1,u_2,u_3,u_4\right)\Vert }^2-{\omega }_1,\end{array}$

where $W=\frac{1}{2}\min \left\{{\omega }_1,{\omega }_2,{\omega }_3,{\omega }_4\right\}$ .

The above establishes a bound on $L$ .

Thus, we have a unique solution of the optimality system for small time intervals due to the opposite time orientations of the state equations and the adjoint equations [4] . Moreover, the uniqueness of the solution of the optimality system guarantees the uniqueness of the optimal control if it exists.

4.2. Characterization of the Optimal Control

Using Pontryagin’s Maximal principle [16] , we obtain

#Math_203# (14)

where Hamiltonian is defined as

$\begin{array}{c}H=L+{\lambda }_1\stackrel{\dot{}}{S}+{\lambda }_2{\stackrel{\dot{}}{I}}_1+{\lambda }_3{\stackrel{\dot{}}{I}}_2+{\lambda }_4\stackrel{\dot{}}{T}+{\lambda }_5\stackrel{\dot{}}{A}+{\lambda }_6\stackrel{\dot{}}{R}\\ =a_{1}S+a_2{I}_1+a_3{I}_2+a_{4}A+\frac{1}{2}\left({\omega }_1{u}_1^2+{\omega }_2{u}_2^2+{\omega }_3{u}_3^2+{\omega }_4{u}_4^2\right)\\ +{\lambda }_1\left\{\left(1-{\epsilon }_1-{\epsilon }_2\right)Q_0-\left[{\beta }_1{I}_1+\left(1-u_1\right){\beta }_2{I}_2\right]S-\left(d+u_2\right)S\right\}\\ +{\lambda }_2\left\{u_3{\epsilon }_1{Q}_0+\left[{\beta }_1{I}_1+\left(1-u_1\right){\beta }_2{I}_2\right]S-\left(u_4+k_2+d\right)I_1\right\}\\ +{\lambda }_3\left[\left(1-u_3\right){\epsilon }_1{Q}_0+u_4{I}_1-\left(k_3+{\mu }_1+d\right)I_2\right]\\ +{\lambda }_4\left[{\mu }_1{I}_2+{\mu }_{2}A-\left(d+{\sigma }_2\right)T\right]\\ +{\lambda }_5\left[k_2{I}_1+k_3{I}_2-\left({\mu }_2+d+{\sigma }_1\right)A\right]\\ +{\lambda }_6\left[{\epsilon }_2{Q}_0+u_{2}S-dR\right].\end{array}$

$\begin{array}{l}\frac{\partial H}{\partial u_1}={\omega }_1{u}_1+\left({\lambda }_1-{\lambda }_2\right){\beta }_2{I}_{2}S=0,\text{at}{u}_1=u_1^{\ast },\\ \frac{\partial H}{\partial u_2}={\omega }_2{u}_2+\left({\lambda }_6-{\lambda }_1\right)S=0,\text{at}{u}_2=u_2^{\ast },\\ \frac{\partial H}{\partial u_3}={\omega }_3{u}_3+\left({\lambda }_2-{\lambda }_3\right){\epsilon }_1{Q}_0=0,\text{at}{u}_3=u_3^{\ast },\\ \frac{\partial H}{\partial u_4}={\omega }_4{u}_4+\left({\lambda }_3-{\lambda }_2\right)I_1=0,\text{at}{u}_4=u_4^{\ast }.\end{array}$

Hence solving for $u_1^{\ast }$ , $u_2^{\ast }$ , $u_3^{\ast }$ and $u_4^{\ast }$ we get

$\begin{array}{l}{u}_1^{\ast }=\frac{{\beta }_2{I}_{2}S\left({\lambda }_2-{\lambda }_1\right)}{{\omega }_1},\\ u_2^{\ast }=\frac{S\left({\lambda }_1-{\lambda }_6\right)}{{\omega }_2},\\ u_3^{\ast }=\frac{{\epsilon }_1{Q}_0\left({\lambda }_3-{\lambda }_2\right)}{{\omega }_3},\\ u_4^{\ast }=\frac{I_1\left({\lambda }_2-{\lambda }_3\right)}{{\omega }_4},\end{array}$

We can now impose the bounds $0\le u_i\le u_{i\max },i=1,2,3,\mathrm{4,}$ on the controls to get

$\begin{array}{l}{u}_1^{\ast }=\min \left\{\max \left(0,\frac{{\beta }_2{I}_{2}S\left({\lambda }_2-{\lambda }_1\right)}{{\omega }_1}\right),u_{1\max }\right\},\\ u_2^{\ast }=\min \left\{\max \left(0,\frac{S\left({\lambda }_1-{\lambda }_6\right)}{{\omega }_2}\right),u_{2\max }\right\},\\ u_3^{\ast }=\min \left\{\max \left(0,\frac{{\epsilon }_1{Q}_0\left({\lambda }_3-{\lambda }_2\right)}{{\omega }_3}\right),u_{3\max }\right\},\\ u_4^{\ast }=\min \left\{\max \left(0,\frac{I_1\left({\lambda }_2-{\lambda }_3\right)}{{\omega }_4}\right),u_{4\max }\right\},\end{array}$

or

$u_1^{\ast }=(\begin{array}{cc}0& \text{when}{\lambda }_2-{\lambda }_1>0,\\ \min \left\{\frac{{\beta }_2{I}_{2}S\left({\lambda }_2-{\lambda }_1\right)}{{\omega }_1},u_{1\max }\right\}& \text{when}{\lambda }_2-{\lambda }_1<0,\end{array}$

$u_2^{\ast }=(\begin{array}{cc}0& \text{when}{\lambda }_1-{\lambda }_6>0,\\ \min \left\{\frac{S\left({\lambda }_1-{\lambda }_6\right)}{{\omega }_2},u_{2\max }\right\}& \text{when}{\lambda }_1-{\lambda }_6<0,\end{array}$

$u_3^{\ast }=(\begin{array}{cc}0& \text{when}{\lambda }_3-{\lambda }_2>0,\\ \min \left\{\frac{{\epsilon }_1{Q}_0\left({\lambda }_3-{\lambda }_2\right)}{{\omega }_3},u_{3\max }\right\}& \text{when}{\lambda }_3-{\lambda }_2<0,\end{array}$

$u_4^{\ast }=(\begin{array}{cc}0& \text{when}{\lambda }_2-{\lambda }_3>0,\\ \min \left\{\frac{I_1\left({\lambda }_2-{\lambda }_3\right)}{{\omega }_4},u_{4\max }\right\}& \text{when}{\lambda }_2-{\lambda }_3<0.\end{array}$

5. Numerical Simulations

In this section, first we present some numerical results of the system (1), when $R_1=1.622>1$ using initial conditions and parameter values given in Table 1(a)

(a) (b)

Table 1. (a) Initial states; (b) parameter values.

and Table 1(b). In our simulations, we consulted some stastical reports [17] [18] and similar works [2] [9] in the field for the parameter values Q₀, d, μ₁, ${\sigma }_1$ and ${\sigma }_2$ . All the other parameter values are fitted. The stability of the model (1) is presented in Figure 1(a) and Figure 1(b). In these figures we choose to separate the disease classes from the non-disease classes for their different population scales, we also provide the graphics for the different time scales to clarify the numerical findings. These figures shows that all model variables approaches to the endemic equilibrium of the model (1), namely, $E^{\ast }=\left(5.52\times {10}^6,416.919,144.81,104.69,21,2.8\times {10}^6\right)$ .

Now we solved the optimality system numerically by using fourth-order iterative Runge-Kutta scheme and presented some graphical results which mostly are comparative figures of susceptibles $\left(S\left(t\right)\right)$ , unaware infectives $\left(I_1\left(t\right)\right)$ , aware infectives $\left(I_2\left(t\right)\right)$ and infectives who receiving treatment $\left(T\left(t\right)\right)$ . We investigated the dynamics of compartments for varying combination of control measures $u_1$ , $u_2$ , $u_3$ and $u_4$ and present some graphics that show the most explicit difference between variables before and after applying the control measures.. The numerical results for optimality system obtained by using the parameter values given in Table 1(a) and Table 1(b), and for weight constants from the objective function, we set $a_1=35$ , $a_1=100$ , $a_1=50$ , $a_1=50$ , ${\omega }_1=70$ , ${\omega }_2=10$ , ${\omega }_3=10$ and ${\omega }_4=100$ .

In Figures 2(a)-(c), we present some graphics when we apply all combination of control measures to the system. By comparison its easy to see that we have some positive results immediately after the time break. The decrease in number of susceptible $\left(S\right)$ and unaware infectives $\left(I_1\right)$ can be interpreted as there are less people get infected and less people to spread the disease with higher transmission probability. In the contrary, population of recovered class $\left(R\right)$ increased and this is very anticipated, ideal state in the war against HIV/AIDS which has no cure. These effects of control measures to the above three classes may explicitly show the essence of this paper. Because persuading more people

(a) (b)

Figure 1. (a) The global asymptotically stability of E*. (b) The global asymptotically stability of E*.

(a) (b)(c) (d)(e) (f)(g) (h)(i)

Figure 2. (a) The control figure for susceptibles (S); (b) the control figure for unaware infectives (I₁); (c) the control figure for aware infectives (I₂); (d) the control figure for treatment class (T); (e) the control figure for recovered class (R); (f) the control u₁; (g) the control u₂; (h) the control u₃; (i) the control u₄.

to practise safe sex measures and efficient medical testing is far more cost- efficient. From Figure 2(c) and Figure 2(d), as an immediate effect of control measures, we depicted a tolerable increase in the number of aware infectives $\left(I_2\right)$ and infectives who receiving treatment $\left(T\right)$ .

In Figures 3(a)-(c), we consider implying two control measures $u_2$ and $u_3$ . In this case, number of susceptible significantly reduced and recovered class expand rapidly. The main characterization of this combination of control measures would be the static increase of aware infectives than the one there is no control.

6. Conclusion

In this paper, we formulated a deterministic model for controlling HIV/AIDS disease. we proved that our system only has one endemic equilibrium and it’s

(a) (b)(c) (d)(e) (f)

Figure 3. (a) The control figure for susceptibles (S); (b) the control figure for unaware infectives (I₁); (c) the control figure for aware infectives (I₂); (d) the control figure for recovered class (R); (e) the control u₂, (f) the control u₃.

globally asymptotically stable if threshold-like conditions satisfied. And then we presented an optimal control problem. Our aim is to investigate combined role of the human behavior change and medical screening in the transmission of HIV/AIDS. We proved the existence and uniqueness of the optimal control and characterized the controls using Pontryagon’s Maximal Principal. In the end we solved the optimality system numerically, and results once again shows us that the optimal way of preventing further prevalence of HIV/AIDS is practicing safe sexual behaviors and conducting efficient, affordable testing as soon as possible for all people. The initial conditions we used in our simulations come from the statistical data, we have fixed some parameters and introduced some of the parameter values from similar works in this field to show our analytical findings.

Funding

This work was supported by the National Natural Science Foundation of China (Grants Nos. 11261056, 11271312).

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Seidu, B. and Makinde, O.D. (2014) Optimal Control of HIV/AIDS in the Workplace in the Presence of Careless Individuals. Computational & Mathematical Methods in Medicine, 2014, 119-128.
[2]	Sharma, S. and Samanta, G.P. (2014) Dynamical Behaviour of an HIV/AIDS Epidemic Model. Differential Equations and Dynamical Systems, 22, 369-395. https://doi.org/10.1007/s12591-013-0173-7
[3]	Okosun, K.O., Makinde, O.D. and Takaidza, I. (2013) Impact of Optimal Control on the Treatment of HIV/AIDS and Screening of Unaware Infectives. Applied Mathematical Modelling, 37, 3802-3820. https://doi.org/10.1016/j.apm.2012.08.004
[4]	Yusuf, T.T. and Benyah, F. (1969) Optimal Strategy for Controlling the Spread of HIV/AIDS Disease: A Case Study of South Africa. Journal of Biological Dynamics, 6, 475-494. https://doi.org/10.1080/17513758.2011.628700
[5]	Garnett, G.P. and Anderson, R.M. (1996) Antiviral Therapy and the Transmission Dynamics of HIV-1. Journal of Antimicrobial Chemotherapy, 37, 135-50. https://doi.org/10.1093/jac/37.suppl_B.135
[6]	Castillo-Chavez, C. and Feng, Z. (1998) Global Stability of an Age-Structure Model for TB and Its Applications to Optimal Vaccination Strategies. Mathematical Biosciences, 151, 135-154. https://doi.org/10.1016/S0025-5564(98)10016-0
[7]	Li, G.H., Wang, W.D. and Jin, Z. (2006) Global Stability of an SEIR Epidemic Model with Constant Immigration. Chaos, Solitons and Fractals, 30, 1012-1019. https://doi.org/10.1016/j.chaos.2005.09.024
[8]	Naresh, R., Tripathi, A. and Sharma, D. (2009) Modelling and Analysis of the Spread of AIDS Epidemic with Immigration of HIV Infectives. Mathematical and Computer Modelling, 49, 880-892. https://doi.org/10.1016/j.mcm.2008.09.013
[9]	Tripathi, A., Naresh, R., Tchuenche, J.M. and Sharma, D. (2013) Modeling the Spread of HIV-AIDS with Infective Immigrants and Time Delay. International Journal of Nonlinear Science, 16, 313-322.
[10]	Coffeea, M., Lurieb, M.N. and Garnett, G.P. (2007) Modelling the Impact of Migration on the HIV Epidemic in South Africa. AIDS, 21, 343-350. https://doi.org/10.1097/QAD.0b013e328011dac9
[11]	Brauer, F. and van den Driessche, P. (2001) Models for Tranmission of Disease with Immigration of Infectives. Mathematical Biosciences, 171, 143-154. https://doi.org/10.1016/S0025-5564(01)00057-8
[12]	Tumwiine, J., Mugisha, J.Y.T. and Luboobi, L.S. (2010) A Host-Vector Model for Malaria with Infective Immigrants. Journal of Mathematical Analysis & Applications, 361, 139-149. https://doi.org/10.1016/j.jmaa.2009.09.005
[13]	Peng, X. (2011) China’s Demographic History and Future Challenges. Science, 333, 581-587. https://doi.org/10.1126/science.1209396
[14]	Li, M.Y. and Muldowney, J.S. (1996) A Geometric Approach to Global-Stability Problems. Siam Journal on Mathematical Analysis, 27, 1070-1083. https://doi.org/10.1137/S0036141094266449
[15]	Coddington, E.A., Levinson, N. and Teichmann, T. (1991) Theory of Ordinary Differential Equations. Mathematical Gazette, 75, 18.
[16]	Pontriagin, L.S., Boltyanskii, V.G, Gamkrelidze, R.V., et al. (1962) The Mathematical Theory of Optimal Processes. Interscience Publishers, New York, 493-494.
[17]	http://www.stats.gov.cn/tjsj/ndsj/2014/indexch.htm
[18]	http://www.xjtj.gov.cn/