Modelling the transmission dynamics and optimal control strategies for HIV infection in China

ABSTRACT

In order to end the AIDS epidemic by 2030 that was put forward by the Joint United Nations Programme on HIV/AIDS in 2014, China needs to take more effective measures to achieve the three 90% goals (90-90-90). We establish a compartmental model to study the dynamics of HIV transmission with control strategies. The analytical results show the existence and stability of the disease-free equilibrium and endemic equilibrium. An optimal control model is constructed to evaluate the impacts of control measures. The simulation results show that the optimal control strategy proposed in this work can eradicate AIDS by 2030. The cost-effectiveness analysis indicates that the cost of the control strategy that combines screening for latent individuals and enhancing education for unaware infected individuals is the lowest. Our findings can provide guidance for public health authorities on effective mitigation strategies to achieve the goals proposed by the United Nations Program on HIV/AIDS.

KEYWORDS:

• HIV/AIDS
• equilibrium
• mitigation strategies
• cost-effectiveness

1. Introduction

Since the first case of Human Immunodeficiency Virus (HIV) infection was detected in the United States in 1981, Acquired Immune Deficiency Syndrome (AIDS) caused by HIV infection has widely spread in different countries and regions. AIDS is a serious life-threatening disease, and presents a great challenge to public health authorities [3,4,13]. In 2014, the United Nations Programme on AIDS (UNAIDS) proposed to achieve three 90%s by 2020, so as to end global AIDS epidemic by 2030. The three 90% goals include: 90% of people living with HIV need to be diagnosed, 90% of those diagnosed need to receive antiretroviral therapy (ART), and 90% of those under treatment need to achieve viral suppression [24]. At present, no vaccines are available to prevent HIV infection [2,6,23,34]. HIV/AIDS has long latent period, strong infectivity, and high death rate. Many countries are affected by the transmission of AIDS, such as China, Cape Verde, Kenya, Lesotho, Malawi, and Nigeria [11].

HIV can be transmitted through blood transmission, sexual transmission, and mother-to-child transmission. People infected with HIV eventually show immunodeficiency syndrome over time. HIV aims to invade the immune system, and the infection process includes four periods, namely, acute infectious period, latent period, pre AIDS period, and AIDS period. Latent individuals and unaware infected individuals can transmit AIDS to susceptible individuals, which makes it difficult to control the spread of AIDS [16]. In order to reduce the infection transmitted by latent individuals and unaware infected individuals, it is essential to improve the screening rate and education rate of infected individuals [23].

Mathematical modelling is a commonly used tool to study the spread of infectious diseases [9,18,32,33,38,42,43,46,48]. The factors that affect the transmission of infectious diseases can be identified by mathematical models [17]. Many mathematical models have been developed to study the transmission dynamics of HIV. Tripathi et al. analysed the impact of screening unconscious infected individuals [45]. Suryanto et al. proposed an ODE model to analyse the influence of awareness, education, screening, and treatment on the spread of AIDS, and proved that education, screening, and treatment can reduce the number of infected individuals [44]. An AIDS model with awareness and control strategy was proposed in [35], where susceptible individuals were mainly classified by gender, and infected individuals were classified according to the level of awareness. Makinde et al. analysed the effects of screening and treatment on the spread of HIV/AIDS infection in the population, and incorporated the awareness and treatment for HIV infected individuals into the model [39]. Nyabadza et al. proposed an HIV/AIDS model with screening and treatment, and showed that the number of aware infected individuals has a great impact on the spread of AIDS [30]. The above works all take into account the impacts of awareness on HIV transmission in their models. However, none of them explore how control strategies can be implemented to achieve the goal of three $90\mathrm{\%}$s and ending the AIDS epidemic by 2030 proposed by the UNAIDS in 2014. Therefore, the purpose of our work is to explore how to achieve these goals, and provide guidance to public health authorities on containing the spread of HIV.

Effective education and intervention for individuals at risk of HIV infection require a large amount of resources. Since HIV infection cannot be completely cured at present, the infected individuals who are receiving ART treatment need to take medicine for a long term. The resources for mitigating HIV infection become more scarce since the outbreaks of COVID-19. In this case, it is of great significance to study how to allocate limited healthcare resources effectively to improve the efficiency of prevention and control for AIDS.

We incorporate mitigation strategies, including providing treatment for aware infected individuals, screening latent individuals, and educating unaware infected individuals, into the model. We apply Pontryagin's Maximum Principle [31] to derive optimal control strategies for achieving the goals of end AIDS. The results show that the optimal control strategy is the combination of all three control measures. The control effect  is closely linked to the weight [15,19,21,25,26].

The rest of this work is organized as follows. In Section 2, we develop a mathematical model to study the transmission dynamics of HIV, derive the basic reproduction number, and analyse the stability of the disease free equilibrium and endemic equilibrium. We exhibit how transimation rate, screening rate, and education rate impact the spread of HIV through numerical simulations. In Section 3, we construct a model incorporating control strategies to find optimal control strategies for curbing the spread through numerical simulations. In Section 4, we summarize our work, and propose future work.

2. Mathematical formulation and dynamics analysis

In this section, we introduce an HIV model and verify the feasible region of the solution. We also calculate the basic reproduction number, and analyse the stability of the disease free equilibrium and the endemic equilibrium. In addition, we perform numerical simulations.

2.1. Model formulation

The infected individuals are divided into latent individuals, aware individuals, and unaware individuals. We assume that all aware infected individuals seek treatments and do not spread the infection. On the contrary, unaware infected individuals do not know that they are infected and are more likely to infect others. Therefore, unaware infected individuals are not treated unless they are screened. In addition, HIV testing has not been widely applied, which is one of the reasons why the unaware infected individuals are not treated.

In this work, we construct a mathematical model for HIV/AIDS transmission to analyse the spread of HIV among individuals. The population is divided into six compartments, namely, susceptible (S), latent (E), aware ($I_1$), unaware ($I_2$), treated (T), and AIDS (A) compartments. The total number of individuals is $N\left(t\right)$ at time t, where $N\left(t\right)=S\left(t\right)+E\left(t\right)+I_1\left(t\right)+I_2\left(t\right)+T\left(t\right)+A\left(t\right)$. The initial conditions are $S\left(0\right)$, $E\left(0\right)$, $I_1\left(0\right)$, $I_2\left(0\right)$, $T\left(0\right)$, and $A\left(0\right)$.

Typically, all infected individuals are likely to transmit the disease. The infectivities of infected individuals in different compartments are different. In this work, we assume that aware individuals do not spread the infection. Let $\lambda \left(t\right)$ be the force of infection for individuals who are infected by HIV as follows $\lambda \left(t\right)=\text{β}ϵE\left(t\right)+\text{β}{I}_2\left(t\right),$where β is the effective contact rate for HIV transmission and ε is a coefficient representing reduced infectivity of latent individuals compared with infectious individuals. Here, the parameter ε satisfies $0<ϵ\le 1$. The recruitment rate of susceptible individuals is Λ. Susceptible individuals may be infected by contacting infectious individuals at the rate $\lambda \left(t\right)$. A fraction of latent individuals, p, transfer to aware class, and 1−p of latent individuals proceed to unaware infected class. Latent individuals become symptomatic infected at the rate α. Unaware infected individuals become aware at the rate δ. Aware infected individuals receive ART treatment at the rate ξ. The infected individuals who receive ART treatment transfer to latent class due to treatment failure at the rate η. The infected individuals transfer to AIDS class at the rate γ. The natural death rate is μ and the mortality rate of infected individuals is ${\mu }_0$. The schematic diagram diagram is shown in Figure 1. The state variables and the parameters are listed in Tables 1 and 2, respectively. Based on the transmission mechanism, we propose the following ordinary differential equation model. (1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-\text{β}(ϵE+I_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=\text{β}(ϵE+I_2)S-(\alpha +\mu )E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E+\delta I_2-(\mu +\xi +\gamma )I_1,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-(\delta +\gamma +\mu )I_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-(\mu +\eta )T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma (I_1+I_2)+\eta T-{\mu }_{0}A-\mu A.\end{array}$(1) To facilitate calculation and analysis, let $k_1=\alpha +\mu ,k_2=\mu +\xi +\gamma ,k_3=\mu +\gamma +\delta$, and $k_4=\mu +\eta$, then Model (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-\text{β}(ϵE+I_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=\text{β}(ϵE+I_2)S-(\alpha +\mu )E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E+\delta I_2-(\mu +\xi +\gamma )I_1,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-(\delta +\gamma +\mu )I_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-(\mu +\eta )T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma (I_1+I_2)+\eta T-{\mu }_{0}A-\mu A.\end{array}$(1) ) is changed as follows (2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-k_{1}E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E-k_2{I}_1+\delta I_2,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-k_3{I}_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-k_{4}T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(2)

Figure 1. Schematic diagram of HIV transmission.

Table 1. The state variables for Model (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-\text{β}(ϵE+I_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=\text{β}(ϵE+I_2)S-(\alpha +\mu )E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E+\delta I_2-(\mu +\xi +\gamma )I_1,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-(\delta +\gamma +\mu )I_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-(\mu +\eta )T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma (I_1+I_2)+\eta T-{\mu }_{0}A-\mu A.\end{array}$(1) ).

Status variable	Description
S	The susceptible individuals
E	The latent individuals
$I_1$	The aware infected individuals
$I_2$	The unaware infected individuals
T	The HIV-infected individuals who receive ART treatment
A	The AIDS individuals

Table 2. Description of parameters in Model (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-\text{β}(ϵE+I_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=\text{β}(ϵE+I_2)S-(\alpha +\mu )E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E+\delta I_2-(\mu +\xi +\gamma )I_1,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-(\delta +\gamma +\mu )I_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-(\mu +\eta )T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma (I_1+I_2)+\eta T-{\mu }_{0}A-\mu A.\end{array}$(1) ).

Parameter	Description	Unit	Value	Source
Λ	Recruitment rate of the susceptible population	Year${}^{-1}$	1430	–
β	The effective contact rate for HIV transmission	Year${}^{-1}$	$1.7649e^{-6}$	MCMC
α	The rate at which latent individuals transfer to unaware or aware infected class	Year${}^{-1}$	0.333	[41]
ε	Reduced infectivity of latent individuals	Dimensionless	0.769231	[41]
$p_0$	The minimum fraction of aware infected individuals	Dimensionless	0.022591	MCMC
$p_{max}$	The maximum fraction of aware infected individuals	Dimensionless	0.45584	MCMC
a	The constant in the function $p\left(t\right)$	Dimensionless	0.055961	MCMC
δ	The rate at which unaware individuals become aware	Dimensionless	0.015697	MCMC
ξ	The rate at which aware infected individuals receive treatments	Dimensionless	0.34877	MCMC
$\gamma$	The rate at which unaware or aware infected individuals progress to AIDS stage	Dimensionless	0.056044	MCMC
${\mu }_0$	Mortality rate of individuals with AIDS	Year${}^{-1}$	0.2	[30]
μ	Natural death rate	Year${}^{-1}$	0.014286	[41]

2.2. Mathematical analysis

Before analysing Model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-k_{1}E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E-k_2{I}_1+\delta I_2,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-k_3{I}_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-k_{4}T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(2) ), we derive the feasible region, the basic reproduction number, and stability of equilibria.

2.2.1. The feasible region

Here, we analyse the non-negativity of solutions and the feasible region Ω.

Lemma 2.1

Let the initial values of state variables in Model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-k_{1}E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E-k_2{I}_1+\delta I_2,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-k_3{I}_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-k_{4}T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(2) ) be positive, i.e. $S\left(0\right)>0$, $E\left(0\right)>0$, $I_1\left(0\right)>0$, $I_2\left(0\right)>0$, $T\left(0\right)>0$, and $A\left(0\right)>0$. The solutions of $S\left(t\right)$, $E\left(t\right)$, $I_1\left(t\right)$, $I_2\left(t\right)$, $T\left(t\right)$, and $A\left(t\right)$ in Model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-k_{1}E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E-k_2{I}_1+\delta I_2,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-k_3{I}_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-k_{4}T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(2) ) are non-negative for all $t\ge 0$, and $\underset{t\to \mathrm{\infty }}{lim}supN\left(t\right)\le \frac{\mathrm{\Lambda }}{\mu }.$

Proof.

Let $H\left(t\right)$ be the minimum value of all state variables, i.e. $H\left(t\right)=min\left\{S\right(t),E(t),I_1(t),I_2(t),T(t),A(t\left)\right\},$ for t>0. Obviously, $H\left(t\right)>0$. We assume that there exists a minimum time $t_1$ satisfying that $H\left(t_1\right)=0$. Suppose $H\left(t_1\right)=S\left(t_1\right)$, then $E\left(t\right)>0,I_1\left(t\right)>0,I_2\left(t\right)>0,T\left(t\right)>0$, and $A\left(t\right)>0$ for $\mathrm{\forall }t\in [0,t_1)$. When $t\in [0,t_1]$, we obtain that $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda \left(t\right)S-\mu S\\ & \ge -\lambda \left(t\right)S-\mu S\\ & =-(\text{β}ϵE+\text{β}{I}_2)S-\mu S.\end{array}$By integrating both sides of the above inequation from 0 to t, we obtain that $S\left(t\right)\ge S\left(0\right)e^{-{\int }_0^t\left[\text{β}ϵE\right(\tau )+\text{β}{I}_2(\tau )+\mu ]\mathrm{d}\tau }>0,\mathrm{\forall }t\in [0,t_1],$which contradicts with $H\left(t_1\right)=S\left(t_1\right)=0$. Therefore, $S\left(t\right)>0,E\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$ when $t\ge 0$.

Next, we calculate the bounds of the state variables. Note that $N\left(t\right)=S\left(t\right)+E\left(t\right)+I_1\left(t\right)+I_2\left(t\right)+T\left(t\right)+A\left(t\right)$. Adding the equations of Model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-k_{1}E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E-k_2{I}_1+\delta I_2,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-k_3{I}_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-k_{4}T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(2) ), we obtain $\frac{\mathrm{d}N\left(t\right)}{\mathrm{d}t}=\mathrm{\Lambda }-\mu \left(S\right(t)+E(t)+I_1(t)+I_2(t)+T(t\left)\right)-{\mu }_{0}A\left(t\right)-\mu A\left(t\right).$Because ${\mu }_0\ge 0$ and $\begin{array}{rl}\frac{\mathrm{d}N\left(t\right)}{\mathrm{d}t}& =\mathrm{\Lambda }-\mu N\left(t\right)-{\mu }_{0}A\left(t\right),\\ & \le \mathrm{\Lambda }-\mu N\left(t\right),\end{array}$using the Comparison Principle Theorem [28], we obtain that $N\left(t\right)\le N\left(0\right)e^{-\mu t}+\frac{\mathrm{\Lambda }}{\mu }(1-e^{-\mu t}),$where $N\left(0\right)=S\left(0\right)+E\left(0\right)+I_1\left(0\right)+I_2\left(0\right)+T\left(0\right)+A\left(0\right)$. Hence, $\underset{t\to \mathrm{\infty }}{lim}supN\left(t\right)\le \frac{\mathrm{\Lambda }}{\mu }.$

Lemma 2.2

The feasible region Ω of Model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-k_{1}E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E-k_2{I}_1+\delta I_2,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-k_3{I}_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-k_{4}T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(2) ) can be defined by $\mathrm{\Omega }=\{\left(S\left(t\right),E\left(t\right),I_1\left(t\right),I_2\left(t\right),T\left(t\right),A\left(t\right)\right)\in {\mathbb{R}}_{+}^6:0\le S,E,I_1,I_2,T,A\le N,N(t)\le \frac{\mathrm{\Lambda }}{\mu }\}.$

2.2.2. The basic reproduction number

Next, we prove the existence of the disease-free equilibrium point and compute the basic reproduction number for Model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-k_{1}E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E-k_2{I}_1+\delta I_2,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-k_3{I}_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-k_{4}T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(2) ).

Let the right side of Model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-k_{1}E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E-k_2{I}_1+\delta I_2,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-k_3{I}_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-k_{4}T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(2) ) equal to zero, we obtain the following disease-free equilibrium $Q_0=(S_0,E_0,I_{10},I_{20},T_0,A_0)=(\frac{\mathrm{\Lambda }}{\mu },0,0,0,0,0).$The basic reproduction number of Model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-k_{1}E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E-k_2{I}_1+\delta I_2,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-k_3{I}_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-k_{4}T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(2) ) is obtained by the next generation matrix approach [10]. Let ${\mathcal{F}}_i$ be newly infected individuals in compartment i and ${\mathcal{V}}_i$ be the transfer of individuals in each compartment i, where i represents the compartment $S,E,I_1,I_2,T,$ and A. We define $x=(S,E,I_1,I_2,T,A)=(x_1,x_2,x_3,x_4,x_5,x_6)$. Model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-k_{1}E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E-k_2{I}_1+\delta I_2,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-k_3{I}_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-k_{4}T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(2) ) can be rewritten as follows $\frac{\mathrm{d}x}{\mathrm{d}t}=\mathcal{F}-\mathcal{V},$where $\mathcal{F}=\left(\begin{array}{c}0\\ \text{β}ϵES+\text{β}{I}_{2}S\\ 0\\ 0\\ 0\\ 0\end{array}\right),\mathcal{V}=\left(\begin{array}{c}-\mathrm{\Lambda }+(\text{β}ϵE+\text{β}{I}_2)S+\mu S\\ k_{1}E\\ -p\alpha E+k_2{I}_1-\delta I_2\\ -(1-p)\alpha E+k_3{I}_2\\ -\xi I_1+k_{4}T\\ -\gamma I_1-\gamma I_2-\eta T+{\mu }_{0}A+\mu A\end{array}\right).$According to the next generation matrix approach for calculating the basic reproduction number [10], we consider the infected compartments $x_i$, i = 2, 3, 4, 5. At the disease free equilibrium, $Q_0$, we have $F=\left[\frac{{\mathrm{\partial }\mathcal{F}}_i}{\mathrm{\partial }{x}_i}\right(Q_0\left)\right]=\left(\begin{array}{cccc}\text{β}ϵS_0& 0& \text{β}{S}_0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right),$and $V=\left[\frac{{\mathrm{\partial }\mathcal{V}}_i}{\mathrm{\partial }{x}_i}\right(Q_0\left)\right]=\left(\begin{array}{cccc}{k}_1& 0& 0& 0\\ -p\alpha & k_2& -\delta & 0\\ -(1-p)\alpha & 0& k_3& 0\\ 0& -\xi & 0& k_4\end{array}\right),$where $k_1=\alpha +\mu >0,k_2=\mu +\xi +\gamma >0,k_3=\mu +\gamma +\delta >0$, and $k_4=\mu +\eta >0$. In addition, matrices F and V satisfy the assumptions (A1)–(A5) in [10]. By calculating the spectral radius of the next generation matrix $F{V}^{-1}$, we derive the basic reproduction number, ${\mathcal{R}}_0$, as follows ${\mathcal{R}}_0=\rho \left(F{V}^{-1}\right)=\frac{\text{β}ϵS_0}{k_1}+\frac{\text{β}{S}_0\alpha (1-p)}{k_1{k}_3}.$Here, $\frac{\text{β}ϵS_0}{k_1}$ represents the number of humans infected by latent individuals, $\frac{\text{β}{S}_0\alpha (1-p)}{k_1{k}_3}$ indicates the number of humans infected by unaware infected individuals.

2.2.3. The stability of disease-free equilibrium

First, we study the local stability of the disease free equilibrium.

Theorem 2.1

For Model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-k_{1}E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E-k_2{I}_1+\delta I_2,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-k_3{I}_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-k_{4}T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(2) ), the disease free equilibrium, $Q_0$, is locally asymptotically stable when $R_0<1$ in the feasible region Ω.

Proof.

Based on Model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-k_{1}E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E-k_2{I}_1+\delta I_2,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-k_3{I}_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-k_{4}T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(2) ), we construct the Jacobian matrix as follows $\begin{array}{rl}J& =\left(\begin{array}{cccccc}-\left(\lambda \right(t)+\mu )& -\text{β}ϵS& 0& -\text{β}S& 0& 0\\ \lambda \left(t\right)& \text{β}ϵS-k_1& 0& \text{β}S& 0& 0\\ 0& p\alpha & -k_2& \delta & 0& 0\\ 0& (1-p)\alpha & 0& -k_3& 0& 0\\ 0& 0& \xi & 0& -k_4& 0\\ 0& 0& \gamma & \gamma & \eta & -{\mu }_0-\mu \end{array}\right)\\ & =\left(\begin{array}{cc}(J_1{)}_{4\times 4}& 0\\ (J_3{)}_{2\times 4}& (J_4{)}_{2\times 2}\end{array}\right),\end{array}$where $k_1=\alpha +\mu >0,k_2=\mu +\xi +\gamma >0,k_3=\mu +\gamma +\delta >0$, and $k_4=\mu +\eta >0$. By calculation, the roots of $|rI-J\left(Q_0\right)|=0$ are the roots of $|rI-J_1\left(Q_0\right)|=0$ and $|rI-J_4\left(Q_0\right)|=0$. Here, r represents the eigenvalue.

The Jacobian matrix $J_1\left(Q_0\right)$ is $J_1\left(Q_0\right)=\left(\begin{array}{cccc}-\mu & -\text{β}ϵS_0& 0& -\text{β}{S}_0\\ 0& \text{β}ϵS_0-k_1& 0& \text{β}{S}_0\\ 0& p\alpha & -k_2& \delta \\ 0& (1-p)\alpha & 0& -k_3\end{array}\right).$The characteristic equation is $\begin{array}{rl}\left|rI-J_1\left(Q_0\right)\right|& =\left|\begin{array}{cccc}r+\mu & \text{β}ϵS_0& 0& \text{β}{S}_0\\ 0& r-(\text{β}ϵS_0-k_1)& 0& -\text{β}{S}_0\\ 0& -p\alpha & r+k_2& -\delta \\ 0& -(1-p)\alpha & 0& r+k_3\end{array}\right|\\ & =(r+k_2)(r+\mu )\\ & \times [r^2+(k_1+k_3-\text{β}ϵS_0)r-k_3(\text{β}ϵS_0-k_1)-\text{β}{S}_0\alpha (1-p\left)\right]\\ & =0.\end{array}$Clearly, the first and second roots of $|rI-J_1\left(Q_0\right)|=0$ are $r_1=-k_2<0$ and $r_2=-\mu <0$. Next, we mainly analyse the roots of the following equation (3) $\begin{array}{rl}& r^2+(k_1+k_3-\text{β}ϵS_0)r-k_3(\text{β}ϵS_0-k_1)-\text{β}{S}_0\alpha (1-p)\\ & =r^2+m_{1}r+m_0\\ & =0.\end{array}$(3) Here, $\begin{array}{rl}{m}_1& =k_1+k_3-\text{β}ϵS_0\\ m_0& =k_3(k_1-\text{β}ϵS_0)-\text{β}{S}_0\alpha (1-p)\\ & =k_1{k}_3[1-(\frac{\text{β}ϵS_0}{k_1}+\frac{\text{β}\alpha S_0(1-p)}{k_1{k}_3})]\\ & =k_1{k}_3(1-{\mathcal{R}}_0).\end{array}$When $R_0<1$, it is easy to prove that $m_1>0$ and $m_0>0$. Here, $R_0=\frac{\text{β}ϵS_0}{k_1}+\frac{\text{β}{S}_0\alpha (1-p)}{k_1{k}_3}$. By the Routh–Herwitz criteria [8], the roots of Equation (3(3) $\begin{array}{rl}& r^2+(k_1+k_3-\text{β}ϵS_0)r-k_3(\text{β}ϵS_0-k_1)-\text{β}{S}_0\alpha (1-p)\\ & =r^2+m_{1}r+m_0\\ & =0.\end{array}$(3) ) have negative real parts.

For $J_4\left(Q_0\right)$, $\left|rI-J_4\left(Q_0\right)\right|=(r+k_4)(r+{\mu }_0+\mu )=0.$Obviously, the eigenvalues of $J_4\left(Q_0\right)$ are $r_5=-k_4<0$ and $r_6=-{\mu }_0-\mu <0$.

In summary, the disease free equilibrium, $Q_0$, is locally asymptotically stable in Ω when $R_0<1$.

Next, we explore the globally asymptotic stability of the disease free equilibrium.

Theorem 2.2

For Model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-k_{1}E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E-k_2{I}_1+\delta I_2,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-k_3{I}_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-k_{4}T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(2) ), when ${\mathcal{R}}_0<1$, the disease-free equilibrium point, $Q_0$ is globally asymptotically stable.

Proof.

By Theorem 2.1 of [40], we construct a Lyapunov function as follows $L=(ϵk_3+\alpha (1-p\left)\right)E+k_1{I}_2.$Obviously, $L\ge 0$.

Next, we take the derivative of L. $\begin{array}{rl}\frac{\mathrm{d}L}{\mathrm{d}t}& =(ϵk_3+\alpha (1-p\left)\right)\frac{\mathrm{d}E}{\mathrm{d}t}+k_1\frac{\mathrm{d}{I}_2}{\mathrm{d}t}\\ & =(ϵk_3+\alpha (1-p\left)\right)(\text{β}ϵSE+\text{β}S{I}_2-k_{1}E)+\alpha (1-p)k_{1}E-k_1{k}_3{I}_2\\ & =(ϵk_3+\alpha (1-p\left)\right)(\text{β}ϵSE+\text{β}S{I}_2)-k_1{k}_3ϵE-k_1{k}_3{I}_2\\ & \le (ϵk_3+\alpha (1-p\left)\right)(\frac{\mathrm{\Lambda }}{\mu }\text{β}ϵE+\frac{\mathrm{\Lambda }}{\mu }\text{β}{I}_2)-k_1{k}_3ϵE-k_1{k}_3{I}_2\\ & =(ϵE+I_2)k_1{k}_3({\mathcal{R}}_0-1).\end{array}$When ${\mathcal{R}}_0<1$, $\frac{\mathrm{d}L}{\mathrm{d}t}<0$. Therefore, the largest invariant set contained in $\left\{\right(S,E,I_1,I_2,T,A)\in \mathrm{\Omega }:\frac{\mathrm{d}L}{\mathrm{d}t}=0\}$is $\left\{Q_0\right\}$. According to the Lasalle's invariance principle [20], the disease-free equilibrium $Q_0$ is globally asymptotically stable in Ω when $R_0<1$.

2.2.4. Endemic equilibrium and stability of endemic equilibrium

In this section, we prove the existence of a unique endemic equilibrium point $Q^{\ast }$ in Model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-k_{1}E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E-k_2{I}_1+\delta I_2,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-k_3{I}_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-k_{4}T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(2) ). Let $Q^{\ast }=(S^{\ast },E^{\ast },I_1^{\ast },I_2^{\ast },T^{\ast },A^{\ast })$.

From the second to sixth equations in   Model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-k_{1}E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E-k_2{I}_1+\delta I_2,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-k_3{I}_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-k_{4}T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(2) ), we have $\begin{array}{rl}{E}^{\ast }& =\frac{k_3}{(1-p)\alpha }{I_2}^{\ast },{I_1}^{\ast }=\left[\frac{p{k}_3+(1-p)\delta }{(1-p)k_2}\right]{I_2}^{\ast },T^{\ast }=\frac{\xi (p{k}_3+(1-p\left)\delta \right)}{(1-p)k_2{k}_4}{I_2}^{\ast },\\ A^{\ast }& =\frac{\gamma I_2^{\ast }}{\mu +{\mu }_0}[\frac{p{k}_3+\delta (1-p)}{k_2(1-p)}+1]+\frac{\eta \xi {I_2}^{\ast }}{(\mu +{\mu }_0)k_4}\left[\frac{p{k}_3+\delta (1-p)}{k_2(1-p)}\right],\\ S^{\ast }& =\frac{\mathrm{\Lambda }}{\mu +\frac{k_1{k}_3{I_2}^{\ast }}{\alpha (1-p)}\frac{{\mathcal{R}}_0}{S_0}},I_2^{\ast }=\frac{({\mathcal{R}}_0-1)\mu \alpha (1-p)S_0}{k_1{k}_3{\mathcal{R}}_0}.\end{array}$Hence, we obtain that $I_2^{\ast }>0$ iff $R_0>1$. Based on the above analysis, Model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-k_{1}E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E-k_2{I}_1+\delta I_2,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-k_3{I}_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-k_{4}T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(2) ) has a unique endemic equilibrium, $Q^{\ast }$ when ${\mathcal{R}}_0>1$.

Next, we study the global stability of the endemic equilibrium, $Q^{\ast }$.

Theorem 2.3

For Model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-k_{1}E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E-k_2{I}_1+\delta I_2,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-k_3{I}_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-k_{4}T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(2) ), if ${\mathcal{R}}_0>1$, the endemic equilibrium $Q^{\ast }$ is globally asymptotically stable.

Proof.

We define a Lyapunov function [22] $\tilde{V}$ as follows $\begin{array}{rl}\tilde{V}& =(S-S^{\ast }-S^{\ast }ln\frac{S}{S^{\ast }})+B(E-E^{\ast }-E^{\ast }ln\frac{E}{E^{\ast }})+C(I_1-I_1^{\ast }-I_1^{\ast }ln\frac{I_1}{I_1^{\ast }})\\ & +D(I_2-I_2^{\ast }-I_2^{\ast }ln\frac{I_2}{I_2^{\ast }})+G(T-T^{\ast }-T^{\ast }ln\frac{T}{T^{\ast }}).\end{array}$where B, C, D, and G are nonnegative constants. The function $\tilde{V}$ is positive definite and has the first derivative. Next, we verify that $\tilde{V}$ is positive definite. All parameters of Model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-k_{1}E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E-k_2{I}_1+\delta I_2,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-k_3{I}_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-k_{4}T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(2) ) are bounded and positive. Obviously, $\tilde{V}\left(Q^{\ast }\right)=0$. When $S>0,E>0,I_1>0,I_2>0,T>0,A>0$ and $S\ne S^{\ast },E\ne E^{\ast },I_1\ne {I_1}^{\ast },I_2\ne {I_2}^{\ast },T\ne T^{\ast },A\ne A^{\ast }$, we define $\tilde{V}(S,E,I_1,I_2,T)=h\left(S\right)+Bh\left(E\right)+Ch\left(I_1\right)+Dh\left(I_2\right)+Gh\left(T\right),$where $h\left(S\right)=S-S^{\ast }-S^{\ast }ln\frac{S}{S^{\ast }}$, $h\left(E\right)=E-E^{\ast }-E^{\ast }ln\frac{E}{E^{\ast }}$, $h\left(I_1\right)=I_1-{I_1}^{\ast }-{I_1}^{\ast }ln\frac{I_1}{{I_1}^{\ast }}$, $h\left(I_2\right)=I_2-{I_2}^{\ast }-{I_2}^{\ast }ln\frac{I_2}{{I_2}^{\ast }}$, and $h\left(T\right)=T-T^{\ast }-T^{\ast }ln\frac{T}{T^{\ast }}$. Because $h\left(S^{\ast }\right)=0$ and $h^{\prime }\left(S\right)=1-\frac{S^{\ast }}{S}$, $h\left(S\right)\ge 0$. Similarly, $h\left(E\right)\ge 0$, $h\left(I_1\right)\ge 0$, $h\left(I_2\right)\ge 0$, and $h\left(T\right)\ge 0$.

Based on Model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-k_{1}E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E-k_2{I}_1+\delta I_2,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-k_3{I}_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-k_{4}T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(2) ), $\begin{array}{rl}{\tilde{V}}^{\prime }& =(1-\frac{S^{\ast }}{S})S^{\prime }+B(1-\frac{E^{\ast }}{E})E^{\prime }+C(1-\frac{{I_1}^{\ast }}{I_1}){I_1}^{\prime }\\ & +D(1-\frac{{I_2}^{\ast }}{I_2}){I_2}^{\prime }+G(1-\frac{T^{\ast }}{T})T^{\prime }\\ & =(1-\frac{S^{\ast }}{S})(\mathrm{\Lambda }-\text{β}(ϵE+I_2)S-\mu S)+B(1-\frac{E^{\ast }}{E})(\text{β}(ϵE+I_2)S-(\alpha +\mu )E)\\ & +C(1-\frac{{I_1}^{\ast }}{I_1})(p\alpha E+\delta I_2-(\mu +\xi +\gamma )I_1)\\ & +D(1-\frac{{I_2}^{\ast }}{I_2})((1-p)\alpha E-(\delta +\gamma +\mu )I_2)\\ & +G(1-\frac{T^{\ast }}{T})(\xi I_1-(\mu +\eta )T)\\ & =(1-\frac{S^{\ast }}{S})[\text{β}ϵE^{\ast }{S}^{\ast }+\text{β}{I_2}^{\ast }{S}^{\ast }+\mu S^{\ast }-(\text{β}ϵES+\text{β}{I}_{2}S+\mu S)]\\ & +B(1-\frac{E^{\ast }}{E})(\text{β}ϵES+\text{β}{I}_{2}S-\frac{\text{β}ϵE^{\ast }{S}^{\ast }+\text{β}{I_2}^{\ast }{S}^{\ast }}{E^{\ast }}E)\\ & +C(1-\frac{{I_1}^{\ast }}{I_1})(p\alpha E+\delta I_2-\frac{p\alpha E^{\ast }+\delta {I_2}^{\ast }}{I_1^{\ast }}{I}_1)\\ & +D(1-\frac{{I_2}^{\ast }}{I_2})((1-p)\alpha E-\frac{(1-p)\alpha E^{\ast }}{{I_2}^{\ast }}{I}_2)\\ & +G(1-\frac{T^{\ast }}{T})(\xi I_1-\frac{\xi {I_1}^{\ast }}{T^{\ast }}T).\end{array}$Let $\frac{S}{S^{\ast }}=x,\frac{E}{E^{\ast }}=y,\frac{I_1}{{I_1}^{\ast }}=z,\frac{I_2}{{I_2}^{\ast }}=m,\frac{T}{T^{\ast }}=n.$Then (4) $\begin{array}{rl}{\tilde{V}}^{\prime }& =-\mu S^{\ast }\frac{{(1-x)}^2}{x}+\text{β}ϵS^{\ast }{E}^{\ast }(1-\frac{1}{x})(1-xy)+\text{β}{I}_2^{\ast }{S}^{\ast }(1-\frac{1}{x})(1-xm)\\ & +B\text{β}ϵS^{\ast }{E}^{\ast }(1-\frac{1}{y})(xy-y)+B\text{β}{I}_2^{\ast }{S}^{\ast }(1-\frac{1}{y})(xm-y)\\ & +Cp\alpha E^{\ast }(1-\frac{1}{z})(y-z)+C\delta I_2^{\ast }(1-\frac{1}{z})(m-z)\\ & +D\alpha (1-p)E^{\ast }(1-\frac{1}{m})(y-m)+G\xi I_1^{\ast }(1-\frac{1}{n})(z-n)\\ & =-\mu S^{\ast }\frac{{(1-x)}^2}{x}+[\text{β}ϵS^{\ast }{E}^{\ast }+\text{β}{I}_2^{\ast }{S}^{\ast }+B\text{β}ϵS^{\ast }{E}^{\ast }+B\text{β}{I}_2^{\ast }{S}^{\ast }+Cpa{E}^{\ast }\\ & +C\delta I_2^{\ast }+D\alpha (1-p)E^{\ast }+G\xi I_1^{\ast }]+xy(-\text{β}ϵS^{\ast }{E}^{\ast }+B\text{β}ϵS^{\ast }{E}^{\ast })\\ & +\frac{1}{x}(-\text{β}ϵS^{\ast }{E}^{\ast }-\text{β}{I}_2^{\ast }{S}^{\ast })\\ & +y(\text{β}ϵS^{\ast }{E}^{\ast }-B\text{β}ϵS^{\ast }{E}^{\ast }-B\text{β}{I}_2^{\ast }{S}^{\ast }+Cp\alpha E^{\ast }+D\alpha (1-p)E^{\ast })\\ & +m(\text{β}{I_2}^{\ast }{S}^{\ast }+C\delta I_2^{\ast }-D\alpha (1-p)E^{\ast })+xm(-\text{β}{I}_2^{\ast }{S}^{\ast }+B\text{β}{I}_2^{\ast }{S}^{\ast })\\ & +z(-Cp\alpha E^{\ast }-C\delta {I_2}^{\ast }+G\xi I_1^{\ast })-xB\text{β}ϵS^{\ast }{E}^{\ast }-\frac{xm}{y}B\text{β}{I_2}^{\ast }{S}^{\ast }-\frac{y}{z}C\alpha p{E}^{\ast }-\frac{m}{z}C\delta {I_2}^{\ast }\\ & -\frac{y}{m}D\alpha (1-p)E^{\ast }-nG\xi {I_1}^{\ast }-\frac{z}{n}G\xi I_1^{\ast }.\end{array}$(4) Based on [22], we let the coefficients of xy, y, m, and z equal to zero and obtain that (5) $\{\begin{array}{l}-\text{β}ϵS^{\ast }{E}^{\ast }+B\text{β}ϵS^{\ast }{E}^{\ast }=0,\\ \text{β}ϵS^{\ast }{E}^{\ast }-B\text{β}ϵS^{\ast }{E}^{\ast }-B\text{β}{I_2}^{\ast }{S}^{\ast }+Cp\alpha E^{\ast }+D\alpha (1-p)E^{\ast }=0,\\ \text{β}{I_2}^{\ast }{S}^{\ast }+C\delta I_2^{\ast }-D\alpha (1-p)E^{\ast }=0,\\ -Cp\alpha E^{\ast }-C\delta {I_2}^{\ast }+G\xi {I_1}^{\ast }=0.\end{array}$(5) Solving Equation (5(5) $\{\begin{array}{l}-\text{β}ϵS^{\ast }{E}^{\ast }+B\text{β}ϵS^{\ast }{E}^{\ast }=0,\\ \text{β}ϵS^{\ast }{E}^{\ast }-B\text{β}ϵS^{\ast }{E}^{\ast }-B\text{β}{I_2}^{\ast }{S}^{\ast }+Cp\alpha E^{\ast }+D\alpha (1-p)E^{\ast }=0,\\ \text{β}{I_2}^{\ast }{S}^{\ast }+C\delta I_2^{\ast }-D\alpha (1-p)E^{\ast }=0,\\ -Cp\alpha E^{\ast }-C\delta {I_2}^{\ast }+G\xi {I_1}^{\ast }=0.\end{array}$(5) ), we get $B=1,C=0,D=\frac{\text{β}{I}_2^{\ast }{S}^{\ast }}{\alpha (1-p)E^{\ast }},G=0,$where $p\ne 1$.

With the above B, C, D, and G, the Lyapunov function $\tilde{V}$ is $\tilde{V}=(S-S^{\ast }-S^{\ast }ln\frac{S}{S^{\ast }})+(E-E^{\ast }-E^{\ast }ln\frac{E}{E^{\ast }})+\frac{\text{β}{I}_2^{\ast }{S}^{\ast }}{\alpha (1-p)E^{\ast }}(I_2-I_2^{\ast }-I_2^{\ast }ln\frac{I_2}{I_2^{\ast }}).$Thus, $\begin{array}{rl}{\tilde{V}}^{\prime }& =-\mu S^{\ast }\frac{{(1-x)}^2}{x}+[\text{β}ϵS^{\ast }{E}^{\ast }+\text{β}{I}_2^{\ast }{S}^{\ast }+B\text{β}ϵS^{\ast }{E}^{\ast }+B\text{β}{I}_2^{\ast }{S}^{\ast }+Cpa{E}^{\ast }+C\delta I_2^{\ast }\\ & +D\alpha (1-p)E^{\ast }+G\xi I_1^{\ast }]+\frac{1}{x}(-\text{β}ϵS^{\ast }{E}^{\ast }-\text{β}{I}_2^{\ast }{S}^{\ast })\\ & +Cp\alpha E^{\ast }+D\alpha (1-p)E^{\ast }-xB\text{β}ϵS^{\ast }{E}^{\ast }\\ & -\frac{xm}{y}B\text{β}{I_2}^{\ast }{S}^{\ast }-\frac{y}{z}C\alpha p{E}^{\ast }-\frac{m}{z}C\delta {I_2}^{\ast }-\frac{y}{m}D\alpha (1-p)E^{\ast }-nG\xi {I_1}^{\ast }-\frac{z}{n}G\xi I_1^{\ast }\\ & =-\mu S^{\ast }\frac{{(1-x)}^2}{x}+[2\text{β}ϵS^{\ast }{E}^{\ast }+2\text{β}{I}_2^{\ast }{S}^{\ast }+\frac{\text{β}{I}_2^{\ast }{S}^{\ast }}{\alpha (1-p)E^{\ast }}\alpha (1-p)E^{\ast }]\\ & +\frac{1}{x}(-\text{β}ϵS^{\ast }{E}^{\ast }-\text{β}{I}_2^{\ast }{S}^{\ast })\\ & -x\text{β}ϵS^{\ast }{E}^{\ast }-\frac{xm}{y}\text{β}{I_2}^{\ast }{S}^{\ast }-\frac{y}{m}\frac{\text{β}{I}_2^{\ast }{S}^{\ast }}{\alpha (1-p)E^{\ast }}\alpha (1-p)E^{\ast }\\ & =-\mu S^{\ast }\frac{{(1-x)}^2}{x}+\text{β}ϵS^{\ast }{E}^{\ast }(2-\frac{1}{x}-x)+\text{β}{I_2}^{\ast }{S}^{\ast }(3-\frac{1}{x}-\frac{xm}{y}-\frac{y}{m}).\end{array}$Because the arithmetic mean is greater than or equal to the geometrical mean, $2-\frac{1}{x}-x\le 0$ and $3-\frac{1}{x}-\frac{xm}{y}-\frac{y}{m}\le 0$ for $x,y,m>0$. Here, $2-\frac{1}{x}-x=0$ iff $x=1$, and $3-\frac{1}{x}-\frac{xm}{y}-\frac{y}{m}=0$ iff x = 1 and $y=m$. Therefore, ${\tilde{V}}^{\prime }\le 0$ for $x,y,m>0$, and ${\tilde{V}}^{\prime }=0$ iff $x=1$ and $y=m$. According to the LaSalle's Invariable Principle [20], it is easy to prove that the endemic equilibrium $Q^{\ast }$ is globally asymptotically stable in Ω. In summary, the endemic equilibrium $Q^{\ast }$ is globally asymptotically stable when $R_0>1$.

Since the endemic equilibrium $Q^{\ast }$ is globally asymptotically stable, $Q^{\ast }$ is locally asymptotically stable when $R_0>1$.

2.3. Numerical simulations

The HIV case data is obtained from the Chinese Center for Disease Control and Prevention [7]. We mainly consider the spread of HIV among people between 15 and 60 years old and fit the number of infected individuals to Model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-k_{1}E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E-k_2{I}_1+\delta I_2,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-k_3{I}_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-k_{4}T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(2) ). We apply the annual number of new HIV cases in China from 2002 to 2019 by the Markov chain Monte Carlo (MCMC) method to estimate the unknown parameters in Model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-k_{1}E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E-k_2{I}_1+\delta I_2,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-k_3{I}_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-k_{4}T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(2) ) as shown in Table 2. Moreover, we use the Partial Rank Correlation Coefficients (PRCC) to study the global sensitivity of the parameters of Model (2(2) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S,\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-k_{1}E,\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E-k_2{I}_1+\delta I_2,\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-k_3{I}_2,\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-k_{4}T,\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(2) ). The goal is to identify the most important parameter that affects HIV transmission. The PRCCs of the parameters with respect to $R_0$ are listed in Figure 2. Our results show that parameters ε and β are positively correlated with $R_0$, and parameters p, ξ, γ, and δ are negatively correlated. Moreover, we find that parameters β and p are the most sensitive to $R_0$. Hence, reducing the value of parameter β and increasing the value of parameter p will effectively mitigate HIV spread.

Figure 2. The global sensitivity analysis of parameters ε, p, ξ, γ, δ, and β.

We use the AIDS data in 2002 as the initial conditions and assume that $S_0=S\left(0\right)=924,860,000$, $E_0=E\left(0\right)=770,000$, ${I_1}_0=I_1\left(0\right)=40,000$, ${I_2}_0=I_2\left(0\right)=39,000$, $T_0=T\left(0\right)=560$, and $A_0=A\left(0\right)=0$. In the following numerical simulations, we choose the function $p=p\left(t\right)=(p_0-p_{max})e^{-at}+p_{max}$ to fit the fraction of latent individuals that become aware. Here, $p_0$ is the minimum fraction of aware infected individuals, $p_{max}$ is the maximum fraction of aware infected individuals, and a is a constant.

The stability of the model is shown in Figures 3 and 4. Figure 3 shows the stability of the endemic equilibrium point. When ${\mathcal{R}}_0>1$, the endemic equilibrium point tends to be stable. The number of susceptible individuals, $S\left(t\right)$, decreases first, then tends to be stable. The trends of E, $I_1$, $I_2$, and T are similar. They all increase obviously at the early stage, then tend to be stable. In Figure 4, the disease free equilibrium point tends to be stable when ${\mathcal{R}}_0<1$. The number of susceptible individual increases to a larger level, then tends to be stable, while the numbers of exposed, unaware infected, and aware infected individuals decrease first and tend to be stable afterwards.

Figure 3. The stability analysis when the parameters are $\text{β}=1.7649e^{-6}$, $\delta =0.015697$, $\gamma =0.056044$, $\xi =0.34877$, $p_0=0.022591$, a = 0.055961, $p_{max}=0.45584$, and ${\mathcal{R}}_0=1.8281>1$. (a) The number of susceptible individuals (S) varying with time. (b) The number of latent individuals (E) varying with time. (c) The number of aware infected individuals ($I_1$) varying with time. (d) The number of unaware infected individuals ($I_2$) varying with time. (e) The number of individuals who receive ART treatment individuals (T) varying with time.

Figure 4. The stability analysis when the parameters are $\text{β}=1.0e^{-06}$, $\delta =0.01$, $\gamma =0.15$, $\xi =0.34877$, $p_0=0.022591$, a = 0.055961, $p_{max}=0.45584$, and ${\mathcal{R}}_0=0.7699<1$. (a) The number of susceptible individuals (S) varying with time. (b) The number of latent individuals (E) varying with time. (c) The number of aware infected individuals ($I_1$) varying with time. (d) The number of unaware infected individuals ($I_2$) varying with time. (e) The number of individuals who receive ART treatment individuals (T) varying with time.

In the following numerical simulations, we use the parameters in Table 2. Figure 5 shows that the number of unaware infected and treated individuals increase from 2002 to 2019. The numbers of latent individuals decreases from 2002 to 2011, then slowly increases from 2012. In addition, we find that the number of aware infected individuals increases from 2002 to 2005, then this number remains almost the same from 2006 to 2019. Particularly, we find that the total number of $I_1$ and T tends to be 800,000 at the end of 2019, which is close to the data.

Figure 5. The number of exposed, unaware infected, aware infected, and treated individuals varying with time. The parameters in Table 2 are used in the simulation.

Next, we analyse the impact of a single parameter on the AIDS epidemic. We investigate the effects of the transmission rate, β, on the epidemic. In Figure 6, the number of infected individuals during the latent period of AIDS can be effectively controlled by reducing the transmission rate. Moreover, when the transmission rate is $\text{β}=0.6e^{-6}$, the epidemic dies out more quickly. The changes in the transmission rate indicate that the public health authorities need to continue strengthening control efforts to suppress the spread of HIV infection. Therefore, it is necessary to improve the awareness of latent individuals. Figure 7 shows the impacts of changing $p_0$ and $p_{max}$ on the number of latent individuals. We find that the number of latent individuals does not change much when the value of $p_0$ is between $0.022591$ and 0.4. However, the value of $p_{max}$ has a great influence on the number of latent individuals, that is, the larger the value of $p_{max}$, the smaller the number of latent individuals. Figure 8 explores the effect of the rate at which unaware individuals become aware, δ, on the number of infected individuals. We find that the numbers of latent individuals, unaware infected, aware infected individuals, and treated individuals decrease significantly when the value of δ is increased from 0.01 to 0.2. Moreover, the number of latent individuals increases from 2010 to 2030 when the value of δ is very small. The number of latent individuals decreases till the end of 2030 when δ is large enough.

Figure 6. The number of latent individuals, E, varying with β.

Figure 7. The number of latent individuals, E, varying with $p_0$ and $p_{max}$.

Figure 8. The effect of the rate at which unaware infected individuals become aware, δ, on the number of infected individuals. (a) δ= 0.01. (b) δ= 0.05. (c) δ= 0.1. (d) δ= 0.2.

To evaluate the combined effect of β, $p_{max}$, and δ on the spread of HIV, we simulate the number of infected individuals with different combinations of these three values as is shown in Figure 9. The numbers of latent individuals and unaware individuals are decreasing with the decrease of β and the increase of $p_{max}$ and δ as shown in Figure 9(a,c). However, the number of aware individuals increases from 2002 to 2006, then decreases from 2007 to 2019 when $\text{β}\le 1.5649e^{-6}$, $p_{max}\ge 0.5$, and $\delta \ge 0.03$ as shown in Figure 9(b). From Figure 9(d), we find that with the decrease of β and the increase of $p_{max}$ and δ, the number of infected individuals receiving ART treatment decreases from 2002 to 2015, while the situation is opposite after 2015.

Figure 9. The number of infected individuals $I_1$, $I_2$, E, and T when the values of parameters β, δ, and $p_{max}$ are varying. (a) The effect of changing β, δ, and $p_{max}$ on the number of latent individuals (E). (b) The effect of changing β, δ, and $p_{max}$ on the number of aware individuals ($I_1$). (c) The effect of changing β, δ, and $p_{max}$ on the number of unaware individuals ($I_2$). (d) The effect of changing β, δ, and $p_{max}$ on the number of infected individuals receiving ART treatment (T).

Figure 10 explores the effect of the transmission rate, β, the maximum fraction of aware infected individuals, $p_{max}$, and the rate at which unaware individuals become aware, δ, on the number of infected individuals. From Figure 10(a–d), we find that reducing transmission rate, increasing the maximum fraction of aware infected individuals and the rate at which unaware individuals become aware can effectively mitigate the spread of AIDS. Particularly, the numbers of latent, aware, and unaware individuals become zero by the end of 2030 when $\text{β}=1e^{-6}$, $p_{max}=0.7$, and $\delta =0.14$. Therefore, we can eradicate AIDS by reducing transmission rate and raising the awareness of infected individuals.

Figure 10. The variation trend of the number of the infected individuals under different parameter values of β, $p_{max}$, and δ. (a) $\text{β}=1.5649e^{-6}$, $p_{max}=0.5$, and $\delta =0.03$. (b) $\text{β}=1.4649e^{-6}$, $p_{max}=0.55$, and $\delta =0.06$. (c) $\text{β}=1.2e^{-6}$, $p_{max}=0.65$, and $\delta =0.1$. (d) $\text{β}=1e^{-6}$, $p_{max}=0.7$, and $\delta =0.14$.

From the above numerical simulations, the basic reproduction number $R_0=1.8281>1$ when the parameters in Table 2 are applied in the simulation. If the current control remains unchanged, HIV outbreaks will continue, and eventually approach the endemic equilibrium point (see Figure 3). In addition, we find that the fraction of latent individuals who transfer to aware individuals class, p, and the rate at which unaware infected individuals become aware, δ, have a great impact on reducing the number of infected individuals. When more HIV infected individuals are screened or become aware, the number of individuals receiving ART increases. This leads to a significant reduction in the incidence of HIV, which in turn reduces the burden on the health system [14]. Therefore, it is necessary to increase HIV screening rate and enhance media campaigns to increase the awareness about HIV infection.

3. Optimal control

We extend our model by including multiple control strategies. We prove the existence, boundedness and uniqueness of the optimal solution and apply the Pontryagin's Maximum Principle [31] to compute the optimal solution. Moreover, we analyse the efficiency and the costs of the control strategies and provide the most economical control strategy.

3.1. Control model formulation

Model (6(6) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-(\alpha +\mu )E-u_1\left(t\right)E\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E+\delta I_2-(\mu +\xi +\gamma )I_1+u_1\left(t\right)E+u_2\left(t\right)I_2-u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-(\delta +\gamma +\mu )I_2-u_2\left(t\right)I_2\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-(\eta +\mu )T+u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(6) ) has six state variables, that is, S, E, $I_1$, $I_2$, T, and A, and three control variables, namely, $u_1\left(t\right)$, $u_2\left(t\right)$, and $u_3\left(t\right)$. Here, $u_1\left(t\right)$ represents the intensity of screening latent individuals, $u_2\left(t\right)$ represents the intensity of education for unaware infected individuals, and $u_3\left(t\right)$ represents the treatments for aware infected individuals. The control set is $\mathrm{\Omega }=\left\{\right(u_1,u_2,u_3\left)|u_i\right(t)\in L^{\mathrm{\infty }}[0,t_f],0\le u_i(t)\le c_i,0<c_i\le 1,i=1,2,3\},$where $t_f$ is the end time of implementing controls.

Then, the optimal control model is described as follows (6) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-(\alpha +\mu )E-u_1\left(t\right)E\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E+\delta I_2-(\mu +\xi +\gamma )I_1+u_1\left(t\right)E+u_2\left(t\right)I_2-u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-(\delta +\gamma +\mu )I_2-u_2\left(t\right)I_2\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-(\eta +\mu )T+u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(6) The initial conditions satisfy (7) $S\left(0\right)\ge 0,E\left(0\right)\ge 0,I_1\left(0\right)\ge 0,I_2\left(0\right)\ge 0,T\left(0\right)\ge 0,A\left(0\right)\ge 0.$(7) We seek to minimize the number of infected individuals and the cost of applying screening, education and treatment controls. Thus, the objective function is given by (8) $J(u_1,u_2,u_3)={\int }_0^{t_f}\left[m_{1}E+m_2{I}_1+m_3{I}_2+\frac{n_1}{2}{u}_1^2\left(t\right)+\frac{n_2}{2}{u}_2^2\left(t\right)+\frac{n_3}{2}{u}_3^2\left(t\right)\right]\mathrm{d}t,$(8) where $m_1$, $m_2$, and $m_3$ represent the weights for the numbers of latent individuals, aware infected individuals and unaware infected individuals, respectively. The weights $n_1$, $n_2$, and $n_3$ measure the costs of control variables $u_1$, $u_2$, and $u_3$, respectively.

First, we apply the method in [27,47] to prove the existence of optimal solution.

Theorem 3.1

For the objective function $J\left(u\right)$, there exists an optimal solution $u^{\ast }=\{u_1^{\ast },u_2^{\ast },u_3^{\ast }\}\in \mathrm{\Omega }$ such that $J\left(u^{\ast }\right)=J(u_1^{\ast },u_2^{\ast },u_3^{\ast })=\underset{u\in \mathrm{\Omega }}{min}J(u_1,u_2,u_3).$

Proof.

We prove the existence of an optimal strategy $u^{\ast }$. By definition, the control set Ω is closed and convex. The integration of the function J is also concave on Ω. The control system is bounded, which implies the compactness of the optimal control. Furthermore, there exists a constant $\zeta >1$, and positive values $C_1$, and $C_2$, such that $J(u_1,u_2,u_3)\ge C_1{\left(|u_1\left(t\right){|}^2+|u_2\left(t\right){|}^2+|u_3\left(t\right){|}^2\right)}^{\frac{\zeta }{2}}-C_2,$which proves the existence of the optimal control.

In order to find the optimal solution, we construct the following Lagrange function L and Hamiltonian function H $L(t,\varphi ,u)=m_{1}E+m_2{I}_1+m_3{I}_2+\frac{n_1}{2}{u}_1^2\left(t\right)+\frac{n_2}{2}{u}_2^2\left(t\right)+\frac{n_3}{2}{u}_3^2\left(t\right),$ $\begin{array}{rl}H(t,\varphi ,u)& =m_{1}E+m_2{I}_1+m_3{I}_2+\frac{n_1}{2}{u}_1^2\left(t\right)+\frac{n_2}{2}{u}_2^2\left(t\right)+\frac{n_3}{2}{u}_3^2\left(t\right)\\ & +{\lambda }_1(\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S)+{\lambda }_2\left(\right(\text{β}ϵE+\text{β}{I}_2)S-(\alpha +\mu )E-u_1(t\left)E\right)\\ & +{\lambda }_3(p\alpha E+\delta I_2-(\mu +\xi +\gamma )I_1+u_1(t)E+u_2(t)I_2-u_3(t\left)I_1\right)\\ & +{\lambda }_4\left(\right(1-p)\alpha E-(\delta +\gamma +\mu )I_2-u_2(t\left)I_2\right)\\ & +{\lambda }_5(\xi I_1-(\eta +\mu )T+u_3(t\left)I_1\right)+{\lambda }_6(\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A),\end{array}$where $\varphi =(S,E,I_1,I_2,T,A{)}^T$, $u=(u_1,u_2,u_3{)}^T$, and ${\lambda }_i\left(t\right)$, $i=1,2,\dots ,6$, are adjoint variables.

Second, we apply the Pontryagin's Maximum Principle [31] to compute the optimal solution.

Theorem 3.2

For Model (6(6) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-(\alpha +\mu )E-u_1\left(t\right)E\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E+\delta I_2-(\mu +\xi +\gamma )I_1+u_1\left(t\right)E+u_2\left(t\right)I_2-u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-(\delta +\gamma +\mu )I_2-u_2\left(t\right)I_2\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-(\eta +\mu )T+u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(6) ), if the control variables $(u_1^{\ast },u_2^{\ast },u_3^{\ast })$ satisfy that $H(t,\varphi ,u^{\ast })<H(t,\varphi ,u)$, the states variables $S^{\ast \ast },E^{\ast \ast },I_1^{\ast \ast },I_2^{\ast \ast },T^{\ast \ast },$ and $A^{\ast \ast }$ are the solutions. Thus, there exist adjoint variables, ${\lambda }_i\left(t\right),i=1,2,\dots ,6$ satisfying (9) $\begin{array}{rl}{{\lambda }_1}^{\prime }& ={\lambda }_1(\text{β}ϵE+\text{β}{I}_2+\mu )-{\lambda }_2(\text{β}ϵE+\text{β}{I}_2),\\ {{\lambda }_2}^{\prime }& =-m_1+{\lambda }_1\text{β}ϵS-{\lambda }_2(\text{β}ϵS-(\alpha +\mu )-u_1(t\left)\right)-{\lambda }_3(p\alpha +u_1(t\left)\right)-{\lambda }_4(1-p)\alpha ,\\ {{\lambda }_3}^{\prime }& =-m_2+{\lambda }_3(\mu +\xi +\gamma +u_3(t\left)\right)-{\lambda }_5(\xi +u_3(t\left)\right)-{\lambda }_6\gamma ,\\ {{\lambda }_4}^{\prime }& =-m_3+{\lambda }_1\text{β}S-{\lambda }_2\text{β}S-{\lambda }_3(\delta +u_2(t\left)\right)+{\lambda }_4(\delta +\gamma +\mu +u_2(t\left)\right)-{\lambda }_6\gamma ,\\ {{\lambda }_5}^{\prime }& ={\lambda }_5(\eta +\mu )-{\lambda }_6\eta ,\\ {{\lambda }_6}^{\prime }& ={\lambda }_6\left({\mu }_0+\mu \right).\end{array}$(9) The boundary conditions are ${\lambda }_i\left(t_f\right)=0,i=1,2,\dots ,6.$The optimal controls $u_1^{\ast }$, $u_2^{\ast }$, and $u_3^{\ast }$ are $\{\begin{array}{l}{u}_1^{\ast }\left(t\right)=min\{max\{0,u_1^c\},1\},\\ u_2^{\ast }\left(t\right)=min\{max\{0,u_2^c\},1\},\\ u_3^{\ast }\left(t\right)=min\{max\{0,u_3^c\},1\},\end{array}$where $u_1^c={\displaystyle \frac{({\lambda }_2-{\lambda }_3)E^{\ast \ast }}{n_1}}$, $u_2^c={\displaystyle \frac{({\lambda }_4-{\lambda }_3)I_2^{\ast \ast }}{n_2}},$ and $u_3^c={\displaystyle \frac{({\lambda }_3-{\lambda }_5)I_1^{\ast \ast }}{n_3}}$. Hence, the optimal solution is $u_i^{\ast }=\{\begin{array}{lll}0,& if& u_i^c\le 0,\\ u_i^c,& if& 0<u_i^c<1,\\ 1,& if& u_i^c\ge 1,\end{array}$where i = 1, 2, 3.

Proof.

According to the existence of optimal control solutions based on the Pontryagin's Maximum Principle, we derive the differential equation system of the Hamiltonia function H as follows $\{\begin{array}{ll}{\displaystyle \frac{\mathrm{d}{\lambda }_1}{\mathrm{d}t}}& =-{\displaystyle \frac{\mathrm{\partial }H}{\mathrm{\partial }S}}={\lambda }_1(\text{β}ϵE+\text{β}{I}_2+\mu )-{\lambda }_2(\text{β}ϵE+\text{β}{I}_2),\\ {\displaystyle \frac{\mathrm{d}{\lambda }_2}{\mathrm{d}t}}& =-{\displaystyle \frac{\mathrm{\partial }H}{\mathrm{\partial }E}}=-m_1+{\lambda }_1\text{β}ϵS-{\lambda }_2(\text{β}ϵS-(\alpha +\mu )-u_1(t\left)\right)\\ & -{\lambda }_3(p\alpha +u_1(t\left)\right)-{\lambda }_4(1-p)\alpha ,\\ {\displaystyle \frac{\mathrm{d}{\lambda }_3}{\mathrm{d}t}}& =-{\displaystyle \frac{\mathrm{\partial }H}{\mathrm{\partial }{I}_1}}=-m_2+{\lambda }_3(\mu +\xi +\gamma +u_3(t\left)\right)-{\lambda }_5(\xi +u_3(t\left)\right)-{\lambda }_6\gamma ,\\ {\displaystyle \frac{\mathrm{d}{\lambda }_4}{\mathrm{d}t}}& =-{\displaystyle \frac{\mathrm{\partial }H}{\mathrm{\partial }{I}_2}}=-m_3+{\lambda }_1\text{β}S-{\lambda }_2\text{β}S-{\lambda }_3(\delta +u_2(t\left)\right)\\ & +{\lambda }_4(\delta +\gamma +\mu +u_2(t\left)\right)-{\lambda }_6\gamma ,\\ {\displaystyle \frac{\mathrm{d}{\lambda }_5}{\mathrm{d}t}}& =-{\displaystyle \frac{\mathrm{\partial }H}{\mathrm{\partial }T}}={\lambda }_5(\eta +\mu )-{\lambda }_6\eta ,\\ {\displaystyle \frac{\mathrm{d}{\lambda }_6}{\mathrm{d}t}}& =-{\displaystyle \frac{\mathrm{\partial }H}{\mathrm{\partial }A}}={\lambda }_6\left({\mu }_0+\mu \right),\end{array}$and ${\lambda }_i\left(t_f\right)=0$, $i=1,2,\dots ,6$.

Next, for the control variables $u_1$, $u_2$, and $u_3$ in the control set Ω, we derive the partial derivatives of the Hamiltonian function with respect to $u_1$, $u_2$, and $u_3$ as follows $\begin{array}{rl}\frac{\mathrm{\partial }H}{\mathrm{\partial }{u}_1}& =n_1{u}_1\left(t\right)+({\lambda }_3-{\lambda }_2)E^{\ast \ast },\\ \frac{\mathrm{\partial }H}{\mathrm{\partial }{u}_2}& =n_2{u}_2\left(t\right)+({\lambda }_3-{\lambda }_4)I_2^{\ast \ast },\\ \frac{\mathrm{\partial }H}{\mathrm{\partial }{u}_3}& =n_3{u}_3\left(t\right)+({\lambda }_5-{\lambda }_3)I_1^{\ast \ast }.\end{array}$Let $\frac{\mathrm{\partial }H}{\mathrm{\partial }{u}_1}=0$, $\frac{\mathrm{\partial }H}{\mathrm{\partial }{u}_2}=0$, and $\frac{\mathrm{\partial }H}{\mathrm{\partial }{u}_3}=0$, we obtain $\{\begin{array}{ll}{u}_1^{\ast }\left(t\right)& ={\displaystyle \frac{({\lambda }_2-{\lambda }_3)E^{\ast \ast }}{n_1}},\\ u_2^{\ast }\left(t\right)& ={\displaystyle \frac{({\lambda }_4-{\lambda }_3)I_2^{\ast \ast }}{n_2}},\\ u_3^{\ast }\left(t\right)& ={\displaystyle \frac{({\lambda }_3-{\lambda }_5)I_1^{\ast \ast }}{n_3}}.\end{array}$Since the control variable $0\le u_i\le 1,i=1,2,3$, the optimal controls are $\{\begin{array}{l}{u}_1^{\ast }\left(t\right)=min\left\{max\left\{0,{\displaystyle \frac{({\lambda }_2-{\lambda }_3)E^{\ast \ast }}{n_1}}\right\},1\right\}\\ u_2^{\ast }\left(t\right)=min\left\{max\left\{0,{\displaystyle \frac{({\lambda }_4-{\lambda }_3)I_2^{\ast \ast }}{n_2}}\right\},1\right\}\\ u_3^{\ast }\left(t\right)=min\left\{max\left\{0,{\displaystyle \frac{({\lambda }_3-{\lambda }_5)I_1^{\ast \ast }}{n_3}}\right\},1\right\}.\end{array}$

Next, we show that the solutions of Model (6(6) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-(\alpha +\mu )E-u_1\left(t\right)E\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E+\delta I_2-(\mu +\xi +\gamma )I_1+u_1\left(t\right)E+u_2\left(t\right)I_2-u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-(\delta +\gamma +\mu )I_2-u_2\left(t\right)I_2\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-(\eta +\mu )T+u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(6) ) satisfying initial condition (7(7) $S\left(0\right)\ge 0,E\left(0\right)\ge 0,I_1\left(0\right)\ge 0,I_2\left(0\right)\ge 0,T\left(0\right)\ge 0,A\left(0\right)\ge 0.$(7) ) are bounded by applying the method in [12,36].

Theorem 3.3

The corresponding absolutely continuous solution $\left(S\right(t),E(t),I_1(t),I_2(t),T(t),A(t\left)\right)$ to Model (6(6) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-(\alpha +\mu )E-u_1\left(t\right)E\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E+\delta I_2-(\mu +\xi +\gamma )I_1+u_1\left(t\right)E+u_2\left(t\right)I_2-u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-(\delta +\gamma +\mu )I_2-u_2\left(t\right)I_2\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-(\eta +\mu )T+u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(6) ) is defined on the entire interval $[0,t_f]$ for any admissible controls $\left(u_1\right(t),u_2(t),u_3(t\left)\right)$. The state variables $S\left(t\right),E\left(t\right),I_1\left(t\right),I_2\left(t\right),T\left(t\right),A\left(t\right)$ satisfy the following inequalities $\begin{array}{rl}& 0\le S\left(t\right)\le \frac{\mathrm{\Lambda }}{\mu },0\le E\left(t\right)\le \frac{\mathrm{\Lambda }}{\mu },0\le I_1\left(t\right)\le \frac{\mathrm{\Lambda }}{\mu },0\le I_2\left(t\right)\le \frac{\mathrm{\Lambda }}{\mu },\\ & 0\le T\left(t\right)\le \frac{\mathrm{\Lambda }}{\mu },0\le A\left(t\right)\le \frac{\mathrm{\Lambda }}{\mu }.\end{array}$

Proof.

We assume that the solution $\left(S\right(t),E(t),I_1(t),I_2(t),T(t),A(t\left)\right)$ of Model (6(6) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-(\alpha +\mu )E-u_1\left(t\right)E\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E+\delta I_2-(\mu +\xi +\gamma )I_1+u_1\left(t\right)E+u_2\left(t\right)I_2-u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-(\delta +\gamma +\mu )I_2-u_2\left(t\right)I_2\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-(\eta +\mu )T+u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(6) ) is defined in the interval $[0,t^{\ast })$, which is the maximum interval of existence. Without loss of generality, we assume that $t^{\ast }\le t_f$. According to Model (6(6) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-(\alpha +\mu )E-u_1\left(t\right)E\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E+\delta I_2-(\mu +\xi +\gamma )I_1+u_1\left(t\right)E+u_2\left(t\right)I_2-u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-(\delta +\gamma +\mu )I_2-u_2\left(t\right)I_2\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-(\eta +\mu )T+u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(6) ) and the conditions $0\le u_i\le 1$, i = 1, 2, 3, the nonnegativity of the solutions $\left(S\right(t),E(t),I_1(t),I_2(t),T(t),A(t\left)\right)$ is verified in Lemma 2.2. Simultaneously, we obtain the upper bound of the solution of Model (6(6) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-(\alpha +\mu )E-u_1\left(t\right)E\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E+\delta I_2-(\mu +\xi +\gamma )I_1+u_1\left(t\right)E+u_2\left(t\right)I_2-u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-(\delta +\gamma +\mu )I_2-u_2\left(t\right)I_2\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-(\eta +\mu )T+u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(6) ) by Lemma 2.2, that is, $0\le S\left(t\right),E\left(t\right),I_1\left(t\right),I_2\left(t\right),T\left(t\right),A\left(t\right)\le \frac{\mathrm{\Lambda }}{\mu }.$Therefore, the solutions of Model (6(6) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-(\alpha +\mu )E-u_1\left(t\right)E\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E+\delta I_2-(\mu +\xi +\gamma )I_1+u_1\left(t\right)E+u_2\left(t\right)I_2-u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-(\delta +\gamma +\mu )I_2-u_2\left(t\right)I_2\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-(\eta +\mu )T+u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(6) ) satisfying initial condition (7(7) $S\left(0\right)\ge 0,E\left(0\right)\ge 0,I_1\left(0\right)\ge 0,I_2\left(0\right)\ge 0,T\left(0\right)\ge 0,A\left(0\right)\ge 0.$(7) ) are bounded.

Next, we apply the method in [35] to prove the uniqueness of optimal control.

Theorem 3.4

The solution of the optimal system (6(6) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-(\alpha +\mu )E-u_1\left(t\right)E\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E+\delta I_2-(\mu +\xi +\gamma )I_1+u_1\left(t\right)E+u_2\left(t\right)I_2-u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-(\delta +\gamma +\mu )I_2-u_2\left(t\right)I_2\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-(\eta +\mu )T+u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(6) ) is unique for sufficiently small interval $[t_0,t_f]$.

Proof.

We suppose that $(S,E,I_1,I_2,T,A,{\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4,{\lambda }_5,{\lambda }_6)$ and $(\overline{S},\overline{E},\overline{I_1},\overline{I_2},\overline{T},\overline{A},{\overline{\lambda }}_1,{\overline{\lambda }}_2,{\overline{\lambda }}_3,{\overline{\lambda }}_4,{\overline{\lambda }}_5,{\overline{\lambda }}_6)$ are two solutions of models (6(6) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-(\alpha +\mu )E-u_1\left(t\right)E\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E+\delta I_2-(\mu +\xi +\gamma )I_1+u_1\left(t\right)E+u_2\left(t\right)I_2-u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-(\delta +\gamma +\mu )I_2-u_2\left(t\right)I_2\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-(\eta +\mu )T+u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(6) ) and (9(9) $\begin{array}{rl}{{\lambda }_1}^{\prime }& ={\lambda }_1(\text{β}ϵE+\text{β}{I}_2+\mu )-{\lambda }_2(\text{β}ϵE+\text{β}{I}_2),\\ {{\lambda }_2}^{\prime }& =-m_1+{\lambda }_1\text{β}ϵS-{\lambda }_2(\text{β}ϵS-(\alpha +\mu )-u_1(t\left)\right)-{\lambda }_3(p\alpha +u_1(t\left)\right)-{\lambda }_4(1-p)\alpha ,\\ {{\lambda }_3}^{\prime }& =-m_2+{\lambda }_3(\mu +\xi +\gamma +u_3(t\left)\right)-{\lambda }_5(\xi +u_3(t\left)\right)-{\lambda }_6\gamma ,\\ {{\lambda }_4}^{\prime }& =-m_3+{\lambda }_1\text{β}S-{\lambda }_2\text{β}S-{\lambda }_3(\delta +u_2(t\left)\right)+{\lambda }_4(\delta +\gamma +\mu +u_2(t\left)\right)-{\lambda }_6\gamma ,\\ {{\lambda }_5}^{\prime }& ={\lambda }_5(\eta +\mu )-{\lambda }_6\eta ,\\ {{\lambda }_6}^{\prime }& ={\lambda }_6\left({\mu }_0+\mu \right).\end{array}$(9) ). Let $S=e^{\theta t}{p}_1$, $E=e^{\theta t}{p}_2$, $I_1=e^{\theta t}{p}_3$, $I_2=e^{\theta t}{p}_4$, $T=e^{\theta t}{p}_5$, $A=e^{\theta t}{p}_6$, ${\lambda }_1=e^{-\theta t}{q}_1$, ${\lambda }_2=e^{-\theta t}{q}_2$, ${\lambda }_3=e^{-\theta t}{q}_3$, ${\lambda }_4=e^{-\theta t}{q}_4$, ${\lambda }_5=e^{-\theta t}{q}_5$, ${\lambda }_6=e^{-\theta t}{q}_6$, $\overline{S}=e^{\theta t}{\overline{p}}_1$, $\overline{E}=e^{\theta t}{\overline{p}}_2$, $\overline{I_1}=e^{\theta t}{\overline{p}}_3$, $\overline{I_2}=e^{\theta t}{\overline{p}}_4$, $\overline{T}=e^{\theta t}{\overline{p}}_5$, $\overline{A}=e^{\theta t}{\overline{p}}_6$, ${\overline{\lambda }}_1=e^{-\theta t}{\overline{q}}_1$, ${\overline{\lambda }}_2=e^{-\theta t}{\overline{q}}_2$, ${\overline{\lambda }}_3=e^{-\theta t}{\overline{q}}_3$, ${\overline{\lambda }}_4=e^{-\theta t}{\overline{q}}_4$, ${\overline{\lambda }}_5=e^{-\theta t}{\overline{q}}_5$, and ${\overline{\lambda }}_6=e^{-\theta t}{\overline{q}}_6$, where $\theta >0$ is to be determined. $p_i$, $q_i$, ${\overline{p}}_i$, and ${\overline{q}}_i$ $(i=1,2,\dots ,6)$ are variables with respect to t.

Moreover, we have $\begin{array}{rl}{u}_1^{\ast }\left(t\right)& =min\left\{max\left\{0,\frac{(q_2-q_3)p_2}{n_1}\right\},1\right\}\\ u_2^{\ast }\left(t\right)& =min\left\{max\left\{0,\frac{(q_4-q_3)p_4}{n_2}\right\},1\right\}\\ u_3^{\ast }\left(t\right)& =min\left\{max\left\{0,\frac{(q_3-q_5)p_3}{n_3}\right\},1\right\}\end{array}$and $\begin{array}{rl}\overline{u_1^{\ast }}\left(t\right)& =min\left\{max\left\{0,\frac{(\overline{q_2}-{\overline{q}}_3){\overline{p}}_2}{n_1}\right\},1\right\}\\ \overline{u_2^{\ast }}\left(t\right)& =min\left\{max\left\{0,\frac{({\overline{q}}_4-{\overline{q}}_3){\overline{p}}_4}{n_2}\right\},1\right\}\\ \overline{u_3^{\ast }}\left(t\right)& =min\left\{max\left\{0,\frac{({\overline{q}}_3-{\overline{q}}_5){\overline{p}}_3}{n_3}\right\},1\right\}.\end{array}$Next, we substitute the values of $S,E,I_1,I_2,T,A,{\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4,{\lambda }_5,{\lambda }_6$ into Equations (6(6) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-(\alpha +\mu )E-u_1\left(t\right)E\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E+\delta I_2-(\mu +\xi +\gamma )I_1+u_1\left(t\right)E+u_2\left(t\right)I_2-u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-(\delta +\gamma +\mu )I_2-u_2\left(t\right)I_2\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-(\eta +\mu )T+u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(6) ) and (9(9) $\begin{array}{rl}{{\lambda }_1}^{\prime }& ={\lambda }_1(\text{β}ϵE+\text{β}{I}_2+\mu )-{\lambda }_2(\text{β}ϵE+\text{β}{I}_2),\\ {{\lambda }_2}^{\prime }& =-m_1+{\lambda }_1\text{β}ϵS-{\lambda }_2(\text{β}ϵS-(\alpha +\mu )-u_1(t\left)\right)-{\lambda }_3(p\alpha +u_1(t\left)\right)-{\lambda }_4(1-p)\alpha ,\\ {{\lambda }_3}^{\prime }& =-m_2+{\lambda }_3(\mu +\xi +\gamma +u_3(t\left)\right)-{\lambda }_5(\xi +u_3(t\left)\right)-{\lambda }_6\gamma ,\\ {{\lambda }_4}^{\prime }& =-m_3+{\lambda }_1\text{β}S-{\lambda }_2\text{β}S-{\lambda }_3(\delta +u_2(t\left)\right)+{\lambda }_4(\delta +\gamma +\mu +u_2(t\left)\right)-{\lambda }_6\gamma ,\\ {{\lambda }_5}^{\prime }& ={\lambda }_5(\eta +\mu )-{\lambda }_6\eta ,\\ {{\lambda }_6}^{\prime }& ={\lambda }_6\left({\mu }_0+\mu \right).\end{array}$(9) ). Then, we have $\begin{array}{rl}{p}_1^{\prime }+\theta p_1& =\mathrm{\Lambda }{e}^{-\theta t}-(\text{β}ϵp_2+\text{β}{p}_4)p_1{e}^{\theta t}-\mu p_1,\\ p_2^{\prime }+\theta p_2& =(\text{β}ϵp_2+\text{β}{p}_4)p_1{e}^{\theta t}-(\alpha +\mu )p_2-\frac{(q_2-q_3)p_2}{n_1}{p}_2,\\ p_3^{\prime }+\theta p_3& =p\alpha p_2+\delta p_4-(\mu +\xi +\gamma )p_3+\frac{(q_2-q_3)p_2}{n_1}{p}_2\\ & +\frac{(q_4-q_3)p_4}{n_2}{p}_4-\frac{(q_3-q_5)p_3}{n_1}{p}_3,\\ p_4^{\prime }+\theta p_4& =(1-p)\alpha p_2-(\delta +\gamma +\mu )p_4-\frac{(q_4-q_3)p_4}{n_2}{p}_4,\\ p_5^{\prime }+\theta p_5& =\xi p_3-(\eta +\mu )p_5+\frac{(q_3-q_5)p_3}{n_1}{p}_3,\\ p_6^{\prime }+\theta p_6& =\gamma p_3+\gamma p_4+\eta p_5-{\mu }_0{p}_6-\mu p_6,\\ q_1^{\prime }-\theta q_1& =q_1(\text{β}ϵp_2{e}^{\theta t}+\text{β}{p}_4{e}^{\theta t}+\mu )-q_2{e}^{\theta t}(\text{β}ϵp_2+\text{β}{p}_4),\\ q_2^{\prime }-\theta q_2& =-m_1{e}^{\theta t}+q_1{p}_1\text{β}ϵe^{\theta t}-q_2(\text{β}ϵp_1{e}^{\theta t}-(\alpha +\mu )-\frac{(q_2-q_3)p_2}{n_1}{p}_2)\\ & -q_3(p\alpha +\frac{(q_2-q_3)p_2}{n_1}{p}_2)-q_4(1-p)\alpha ,\\ q_3^{\prime }-\theta q_3& =-m_2{e}^{\theta t}+q_3(\mu +\xi +\gamma +\frac{(q_3-q_5)p_3}{n_1}{p}_3)-q_5(\xi +\frac{(q_3-q_5)p_3}{n_1}{p}_3)-q_6\gamma ,\\ q_4^{\prime }-\theta q_4& =-m_3{e}^{\theta t}-q_3(\delta +\frac{(q_4-q_3)p_4}{n_2}{p}_4)+q_4(\delta +\gamma +\mu ++\frac{(q_4-q_3)p_4}{n_2}{p}_4)-q_6\gamma ,\\ q_5^{\prime }-\theta q_5& =q_5(\eta +\mu )-q_6\eta ,\\ q_6^{\prime }-\theta q_6& =q_6({\mu }_0+\mu ).\end{array}$The equations for S and $\overline{S}$, E and $\overline{E}$, $I_1$ and $\overline{I_1}$, $I_2$ and $\overline{I_2}$, T and $\overline{T}$, A and $\overline{A}$ are subtracted. Each equation is multiplied by an appropriate function and integrated from $t_0$ to $t_f$. Thus, $\begin{array}{rl}& {\int }_{t_0}^{t_f}(u_1^{\ast }-\overline{u_1^{\ast }}{)}^2\mathrm{d}t\le B_1{e}^{2\theta t_f}{\int }_{t_0}^{t_f}[|p_2-{\overline{p}}_2{|}^2+|q_2-{\overline{q}}_2{|}^2+|q_3-{\overline{q}}_3{|}^2]\mathrm{d}t,\\ & {\int }_{t_0}^{t_f}(u_2^{\ast }-\overline{u_2^{\ast }}{)}^2\mathrm{d}t\le B_2{e}^{2\theta t_f}{\int }_{t_0}^{t_f}[|p_4-{\overline{p}}_4{|}^2+|q_3-{\overline{q}}_3{|}^2+|q_4-{\overline{q}}_4{|}^2]\mathrm{d}t,\\ & {\int }_{t_0}^{t_f}(u_3^{\ast }-\overline{u_3^{\ast }}{)}^2\mathrm{d}t\le B_3{e}^{2\theta t_f}{\int }_{t_0}^{t_f}[|p_3-{\overline{p}}_3{|}^2+|q_3-{\overline{q}}_3{|}^2+|q_5-{\overline{q}}_5{|}^2]\mathrm{d}t,\end{array}$where $B_1$, $B_2$, and $B_3$ are constants.

Therefore, we obtain $\begin{array}{rl}& \frac{1}{2}(p_1-{\overline{p}}_1{)}^2(t_f)+\theta {\int }_{t_0}^{t_f}|p_1-{\overline{p}}_1{|}^2\mathrm{d}t\\ & \le \text{β}ϵe^{-\theta t_f}{\int }_{t_0}^{t_f}[|p_1-{\overline{p}}_1{|}^2+|p_2-{\overline{p}}_2{|}^2]\mathrm{d}t+\text{β}{e}^{-\theta t_f}{\int }_{t_0}^{t_f}[|p_1-{\overline{p}}_1{|}^2\\ & +|p_4-{\overline{p}}_4{|}^2]\mathrm{d}t+\mu {\int }_{t_0}^{t_f}|p_1-{\overline{p}}_1{|}^2\mathrm{d}t,\end{array}$and $\begin{array}{rl}& \frac{1}{2}(q_1-{\overline{q}}_1{)}^2(t_f)+\theta {\int }_{t_0}^{t_f}|q_1-{\overline{q}}_1{|}^2\mathrm{d}t\le \text{β}ϵe^{2\theta t_f}{\int }_{t_0}^{t_f}[|p_2-{\overline{p}}_2{|}^2+|q_1-{\overline{q}}_1{|}^2]\mathrm{d}t\\ & +\text{β}{e}^{2\theta t_f}{\int }_{t_0}^{t_f}[|p_4-{\overline{p}}_4{|}^2+|q_1-{\overline{q}}_1{|}^2]\mathrm{d}t+\text{β}ϵe^{2\theta t_f}{\int }_{t_0}^{t_f}[|p_2-{\overline{p}}_2{|}^2+|q_2-{\overline{q}}_2{|}^2]\mathrm{d}t\\ & +\text{β}{e}^{2\theta t_f}{\int }_{t_0}^{t_f}[|p_4-{\overline{p}}_4{|}^2+|q_2-{\overline{q}}_2{|}^2]\mathrm{d}t+\mu e^{\theta t_f}{\int }_{t_0}^{t_f}[|p_1-{\overline{p}}_1{|}^2+|q_1-{\overline{q}}_1{|}^2]\mathrm{d}t.\end{array}$Similarly, we derive the inequalities for $p_2$ and ${\overline{p}}_2$, $p_3$ and ${\overline{p}}_3$, $p_4$ and ${\overline{p}}_4$, $p_5$ and ${\overline{p}}_5$, $p_6$ and ${\overline{p}}_6$, $q_2$ and ${\overline{q}}_2$, $q_3$ and ${\overline{q}}_3$, $q_4$ and ${\overline{q}}_4$, $q_5$ and ${\overline{q}}_5$, $q_6$ and ${\overline{q}}_6$. Adding all the inequalities, we have $\begin{array}{rl}& \frac{1}{2}\left[\right(p_1-{\overline{p}}_1{)}^2\left(t_f\right)+(p_2-{\overline{p}}_2{)}^2(t_f)+(p_3-{\overline{p}}_3{)}^2\left(t_f\right)\\ & +(p_4-{\overline{p}}_4{)}^2(t_f)+(p_5-{\overline{p}}_5{)}^2\left(t_f\right)+(p_6-{\overline{p}}_6{)}^2(t_f)\\ & +(q_1-{\overline{q}}_1{)}^2(t_0)+(q_2-{\overline{q}}_2{)}^2\left(t_0\right)+(q_3-{\overline{q}}_3{)}^2(t_0)\\ & +(q_4-{\overline{q}}_4{)}^2(t_0)+(q_5-{\overline{q}}_5{)}^2\left(t_0\right)+(q_6-{\overline{q}}_6{)}^2(t_0\left)\right]\\ & +\theta {\int }_{t_0}^{t_f}[|p_1-{\overline{p}}_1{|}^2+|p_2-{\overline{p}}_2{|}^2+|p_3-{\overline{p}}_3{|}^2\\ & +|p_4-{\overline{p}}_4{|}^2+|p_5-{\overline{p}}_5{|}^2+|p_6-{\overline{p}}_6{|}^2+|q_1-{\overline{q}}_1{|}^2\\ & +|q_2-{\overline{q}}_2{|}^2+|q_3-{\overline{q}}_3{|}^2+|q_4-{\overline{q}}_4{|}^2+|q_5-{\overline{q}}_5{|}^2+|q_6-{\overline{q}}_6{|}^2]\mathrm{d}t\\ & \le ({\tilde{K}}_1+{\tilde{K}}_2{e}^{3\theta t_f}){\int }_{t_0}^{t_f}[|p_1-{\overline{p}}_1{|}^2+|p_2-{\overline{p}}_2{|}^2\\ & +|p_3-{\overline{p}}_3{|}^2+|p_4-{\overline{p}}_4{|}^2+|p_5-{\overline{p}}_5{|}^2+|p_6-{\overline{p}}_6{|}^2\\ & +|q_1-{\overline{q}}_1{|}^2+|q_2-{\overline{q}}_2{|}^2+|q_3-{\overline{q}}_3{|}^2+|q_4-{\overline{q}}_4{|}^2+|q_5-{\overline{q}}_5{|}^2+|q_6-{\overline{q}}_6{|}^2]\mathrm{d}t.\end{array}$From the above equation, we have $\begin{array}{rl}& (\theta -{\tilde{K}}_1-{\tilde{K}}_2{e}^{3\theta t_f}){\int }_{t_0}^{t_f}[|p_1-{\overline{p}}_1{|}^2+|p_2-{\overline{p}}_2{|}^2+|p_3-{\overline{p}}_3{|}^2+|p_4-{\overline{p}}_4{|}^2+|p_5-{\overline{p}}_5{|}^2\\ & +|p_6-{\overline{p}}_6{|}^2+|q_1-{\overline{q}}_1{|}^2+|q_2-{\overline{q}}_2{|}^2+|q_3-{\overline{q}}_3{|}^2+|q_4-{\overline{q}}_4{|}^2\\ & +|q_5-{\overline{q}}_5{|}^2+|q_6-{\overline{q}}_6{|}^2]\mathrm{d}t\le 0,\end{array}$where ${\tilde{K}}_1$ and ${\tilde{K}}_2$ depend on the coefficients and the bounds on $p_i$ and $q_i$, $i=1,2\dots ,6$.

We choose θ such that $\theta >{\tilde{K}}_1+{\tilde{K}}_2$ and $t_f<\frac{1}{3\theta }ln\left(\frac{\theta -{\tilde{K}}_1}{{\tilde{K}}_2}\right)$, then $p_1={\overline{p}}_1$, $p_2={\overline{p}}_2$, $p_3={\overline{p}}_3$, $p_4={\overline{p}}_4$, $p_5={\overline{p}}_5$, $p_6={\overline{p}}_6$, $q_1={\overline{q}}_1$, $q_2={\overline{q}}_2$, $q_3={\overline{q}}_3$, $q_4={\overline{q}}_4$, $q_5={\overline{q}}_5$, $q_6={\overline{q}}_6$. Therefore, Model (6(6) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=\mathrm{\Lambda }-(\text{β}ϵE+\text{β}{I}_2)S-\mu S\\ {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=(\text{β}ϵE+\text{β}{I}_2)S-(\alpha +\mu )E-u_1\left(t\right)E\\ {\displaystyle \frac{\mathrm{d}{I}_1}{\mathrm{d}t}}=p\alpha E+\delta I_2-(\mu +\xi +\gamma )I_1+u_1\left(t\right)E+u_2\left(t\right)I_2-u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}{I}_2}{\mathrm{d}t}}=(1-p)\alpha E-(\delta +\gamma +\mu )I_2-u_2\left(t\right)I_2\\ {\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}=\xi I_1-(\eta +\mu )T+u_3\left(t\right)I_1\\ {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=\gamma I_1+\gamma I_2+\eta T-{\mu }_{0}A-\mu A.\end{array}$(6) ) has a unique optimal solution within a small time interval.

3.2. Numerical simulations

In order to illustrate the feasibility of the theoretical results and the control strategies, numerical simulations are carried out to study the influence of the optimal control strategy on the transmission of HIV.

The weights of state variables and control variables in the Hamiltonian function are set as $m_1=20$, $m_2=30$, $m_3=30$, $n_1=45$, $n_2=25$, and $n_3=10$. Based on the three control variables, $u_1\left(t\right)$, $u_2\left(t\right)$, and $u_3\left(t\right)$, we propose the following four control strategies.

To compare the control effect after implementing the control strategies, we first simulate the number of infected individuals without control as shown in Figure 11. The infected individuals include latent, aware, and unaware infected individuals. Figure 11 shows that the number of infected individuals increases from 2002 to 2014, then decreases from 2015 to 2030. By the end of 2030, the number of infected individuals is 45905000. Next, we explore the impact of control strategies on HIV transmission.

Figure 11. The number of infected individuals when $u_1=0$, $u_2=0$, and $u_3=0$.

Strategy 1: A combination of screening for latent individuals $u_1\left(t\right)$ and education for unaware infected individuals $u_2\left(t\right)$ is applied. In other words, $u_1\ne 0$, $u_2\ne 0$, and $u_3=0$. Then, we optimize the function J in Equation (8(8) $J(u_1,u_2,u_3)={\int }_0^{t_f}\left[m_{1}E+m_2{I}_1+m_3{I}_2+\frac{n_1}{2}{u}_1^2\left(t\right)+\frac{n_2}{2}{u}_2^2\left(t\right)+\frac{n_3}{2}{u}_3^2\left(t\right)\right]\mathrm{d}t,$(8) ).

Figure 12(a) shows the effects of control strategies $u_1$ and $u_2$ over time. Controls $u_1$ and $u_2$ maintain their maximum values for the first four months and five months, respectively. Then, $u_1$ remains constant at the value of 0.95 until 2009. From 2010, control $u_1$ begins to decrease slowly. From May 2002 to January 2030, the control $u_2$ remains constant at the value of 0.95, then slowly decreases to zero. Figure 12(b) shows the effect of controls $u_1$ and $u_2$ on the number of infected individuals. Comparing with the case without control strategies (see Figure 11), implementing controls $u_1$ and $u_2$ can significantly reduce the number of infected individuals. In particular, HIV can be eradicated in 2003.

Figure 12. The impacts of implementing control Strategies 1. (a) The optimal control solutions $u_1$ and $u_2$. (b) The size of $E,I_1,$ and $I_2$ under the controls of $u_1$ and $u_2$.

Strategy 2: A combination of screening for latent individuals $u_1\left(t\right)$ and treatments for aware infected individuals $u_3\left(t\right)$ is considered. In other words, $u_1\ne 0$, $u_3\ne 0$, and $u_2=0$. Then, we optimize the cost function J in Equation (8(8) $J(u_1,u_2,u_3)={\int }_0^{t_f}\left[m_{1}E+m_2{I}_1+m_3{I}_2+\frac{n_1}{2}{u}_1^2\left(t\right)+\frac{n_2}{2}{u}_2^2\left(t\right)+\frac{n_3}{2}{u}_3^2\left(t\right)\right]\mathrm{d}t,$(8) ) to analyse the control effects.

Figure 13(a) shows the effects of controls $u_1$ and $u_3$ over time. As shown in the figure, controls $u_1$ and $u_3$ are at the maximum control level of $100\mathrm{\%}$ for 37 months, then these control measures gradually decrease until reaching $95\mathrm{\%}$. Eventually, the controls $u_1$ and $u_3$ decrease to zero. Figure 13(b) illustrates the effect of controls $u_1$ and $u_3$ on the number of infected individuals. Compared with the case without control strategies, applying controls $u_1$ and $u_3$ could reduce the total number of latent, aware, and unaware infected individuals. In addition, we find that the combined use of controls $u_1$ and $u_3$ could eradicate HIV in 2020. Therefore, optimal control effects can be achieved by screening latent infection individuals and treating the infected individuals during a long period. Then, we can control the spread of HIV by improving the awareness and increasing the opportunities of being treated.

Figure 13. The impacts of implementing control Strategies 2. (a) The optimal control solutions $u_1$ and $u_3$. (b) The total size of $E,I_1,$ and $I_2$ under the controls of $u_1$ and $u_3$.

Strategy 3: A combination of education for infected individuals $u_2\left(t\right)$ and treatment for aware infected individuals $u_3\left(t\right)$ is considered. In other words, $u_1=0$, $u_2\ne 0$, and $u_3\ne 0$. Then, we optimize the cost function J in Equation (8(8) $J(u_1,u_2,u_3)={\int }_0^{t_f}\left[m_{1}E+m_2{I}_1+m_3{I}_2+\frac{n_1}{2}{u}_1^2\left(t\right)+\frac{n_2}{2}{u}_2^2\left(t\right)+\frac{n_3}{2}{u}_3^2\left(t\right)\right]\mathrm{d}t,$(8) ) to analyse the control effect as follows.

Figure 14(a) shows the effects of control strategies $u_2$ and $u_3$. At the early stage, the controls $u_2$ and $u_3$ are at the maximum control level of $100\mathrm{\%}$ for 19 months and 23 months, respectively. Then, the two control strategies keep at $95\mathrm{\%}$, and eventually decrease to zero. Figure 14(b) indicates the effect of controls $u_2$ and $u_3$ on the number of infected individuals. The results show that under the controls $u_2$ and $u_3$, the total numbers of latent, aware, and unaware infected individuals greatly decrease compared with the case without control (see Figure 11). Moreover, we find that HIV is eradicated in 2005 when the controls $u_2$ and $u_3$ are implemented. This indicates that these controls are very effective to mitigate HIV transmission.

Figure 14. The impacts of implementing control Strategies 3. (a) The optimal control solutions $u_2$ and $u_3$. (b) The total size of $E,I_1,$ and $I_2$ under the controls of $u_2$ and $u_3$.

Strategy 4: A combination of screening $u_1\left(t\right)$, education for aware infected individuals $u_2\left(t\right)$ and treatment for aware infected individuals $u_3\left(t\right)$ is considered. In other words, $u_1\ne 0$, $u_2\ne 0$, $u_3\ne 0$. Then, we optimize the cost function J in Equation (8(8) $J(u_1,u_2,u_3)={\int }_0^{t_f}\left[m_{1}E+m_2{I}_1+m_3{I}_2+\frac{n_1}{2}{u}_1^2\left(t\right)+\frac{n_2}{2}{u}_2^2\left(t\right)+\frac{n_3}{2}{u}_3^2\left(t\right)\right]\mathrm{d}t,$(8) ) to analyse the control effect.

Figure 15 indicates the change of the number of infected individuals when three controls are implemented together. In Figure 15(a), the controls $u_1$, $u_2$, and $u_3$ achieve the maximum control level of $100\mathrm{\%}$ and last for about four, five, and seven months, respectively. Then, control $u_1$ decreases to $95\mathrm{\%}$ level until 2011 and slowly decreases to zero in 2030. For control $u_2$, it has maintained a $95\mathrm{\%}$ control level from August 2002 to February 2029, then slowly reduces to zero. However, the control strength of control $u_3$ first decreases and increases from August 2002 to May 2005. Then, control $u_3$ slowly reduces to zero in 2030. Figure 15(b) shows the effect of controls $u_1$, $u_2$, and $u_3$ on the number of infected individuals. The results show that compared with the case without control strategies (see Figure 11), controls $u_1$, $u_2$, and $u_3$ can reduce the total numbers of latent, aware, and unaware infected individuals. It means that under the control strategies, more infected individuals receive treatments, and more latent and unaware infected individuals become aware.

Figure 15. The impacts of implementing control Strategies 4. (a) The optimal control solutions $u_1$, $u_2$, and $u_3$. (b) The total size of $E,I_1,$ and $I_2$ under the control of $u_1$, $u_2$, and $u_3$.

Further, we compare the impact of four control strategies on the number of infected individuals, that is, latent, aware, and unaware infected individuals as shown in Figure 16. We find that Strategy 4: $u_1\ne 0$, $u_2\ne 0$, $u_3\ne 0$ is the most efficient. This is because Strategy 4 can eradicate HIV in the shortest time (about 1.5 years). Other strategies need take more time to achieve the same control effect. To judge whether the four control strategies can achieve the three goals of $90\mathrm{\%}$ reduction and eradicate AIDS by 2030 proposed by the UNAIDS, we summarize the control results in Table 4. Table 4 indicates that four control strategies can achieve the control objectives proposed by the UNAIDS.

Figure 16. The impact of different control strategies on the number of infected individuals.

Next, we study the influence of control weights on Strategy 4. We compare the optimal control curves and the numbers of infected individuals with $n_1=50$, $n_1=100$, and $n_1=150$ when $n_2=25$ and $n_3=10$. In Figure 17, with the increase of $n_1$, the duration of the maximum control level of $u_1$ becomes shorter, while the duration of the maximum control levels of $u_2$ and $u_3$ becomes longer. When $n_1=50$, the controls $u_1$, $u_2$, and $u_3$ reach the $100\mathrm{\%}$ control level at the early stage and last for about four, five, and seven months, respectively, then the control levels reduce to $95\mathrm{\%}$ as shown in Figure 17(a). The control $u_1$ starts to slowly decrease to zero from 2011. The control level of control $u_2$ remains at $95\mathrm{\%}$ level from June 2002 to March 2030, then the control level gradually decreases to zero. The control $u_3$ first decreases, then increases from August 2002 to May 2005. Then, control $u_3$ slowly reduces to zero in 2030. In Figure 17(b), when $n_1=100$, the controls $u_1$, $u_2$, and $u_3$ reach the $100\mathrm{\%}$ control level at the early stage and last for about three, five, and seven months, respectively. Then, the controls $u_1$, $u_2$, and $u_3$ maintain $95\mathrm{\%}$ control level until 2008, 2030, and 2009, respectively, after which they slowly decrease to zero. In Figure 17(c), when $n_1=150$, the controls $u_1$, $u_2$, and $u_3$ reach the $100\mathrm{\%}$ control level and last for about three, five, and seven months, respectively. Then, the controls $u_1$, $u_2$, and $u_3$ maintain $95\mathrm{\%}$ control level until 2008, 2030, and 2009, respectively, after which they slowly decrease to zero. In addition, we explore the effect of different $n_1$ values on the number of infected individuals. Figure 17(d) shows that the smaller the weight $n_1$, the earlier the HIV eradication can be achieved.

Figure 17. Impact of changes on the cost weight of control $u_1\left(t\right)$, $n_1$, on control strategies and the number of infected individuals. (a) Weight: $n_1=50,n_2=25,n_3=75$. (b) Weight: $n_1=100,n_2=25,n_3=10$. (c) Weight: $n_1=150,n_2=25,n_3=10$. (d) The number of infected individuals under different cost weights of $n_1$.

The effect of weight $n_2$ on the control curves and the number of infected individuals with $n_2=20$, $n_2=100$, and $n_2=180$ is shown in Figure 18. With the increase of the weight $n_2$, the maximum control level of control $u_2$ decreases gradually. In Figure 18(a), when $n_2=20$, the controls $u_1$, $u_2$, and $u_3$ reach the $100\mathrm{\%}$ control level and last for about four, five, and seven months, respectively. After the maximum control level, the control $u_1$ starts to decrease slowly after eight years of implementing the control level of $95\mathrm{\%}$, the control $u_2$ maintains a control level of $95\mathrm{\%}$ all the time, while the control $u_3$ first decreases in July 2002 and increases in August 2002, then gradually decrease after June 2004. In Figure 18(b), when $n_2=100$, the controls $u_1$, $u_2$, and $u_3$ reach the $100\mathrm{\%}$ control level at the early stage and last for about four, three, and seven months, respectively. Then, the controls $u_1$, $u_2$, and $u_3$ maintain $95\mathrm{\%}$ control level until 2024, 2016, and 2013, respectively, after which they slowly decrease to zero. In Figure 18(c), when $n_2=180$, the controls $u_1$, $u_2$, and $u_3$ reach the $100\mathrm{\%}$ control level and last for about four, two, and seven months, respectively. Then, the controls $u_1$, $u_2$, and $u_3$ maintain $95\mathrm{\%}$ control level until 2029, 2015, and 2022, respectively, after which they slowly decrease to zero. Moreover, we explore the effect of different values of $n_2$ on the number of infected individuals. Figure 18(d) illustrates that the smaller the weight $n_2$, the better control effect.

Figure 18. Impact of changes on the cost weight of control $u_2\left(t\right)$, $n_2$, on control strategies and the number of infected individuals. (a) Weight: $n_1=45,n_2=20,n_3=10$. (b) Weight: $n_1=45,n_2=100,n_3=10$. (c) Weight: $n_1=45,n_2=180,n_3=10$. (d) The number of infected individuals under different cost weights of $n_2$.

In addition, we explored the effect of different values of $n_3$ on the control strategies and the number of infected individuals. Here, we explore the cases of $n_3=4$, $n_3=8$, and $n_3=12$, respectively. Figure 19 shows that with the increase of the weight $n_3$, the maximum control level of control strategy $u_3$ decreases gradually. In Figure 19(a), when $n_3=4$, the controls $u_1$, $u_2$, and $u_3$ reach the $100\mathrm{\%}$ control level at the early stage and last for about four, five, and eight months, respectively. Then, the controls $u_1$, $u_2$, and $u_3$ maintain $95\mathrm{\%}$ control level until 2010, 2029, and 2014, respectively, after which they slowly decrease to zero. Figure 19(b) shows that when $n_3=8$, the controls $u_1$, $u_2$, and $u_3$ reach the $100\mathrm{\%}$ control level at the early stage and last for about four, five, and seven months, respectively. Then, the controls $u_1$ and $u_2$ maintain $95\mathrm{\%}$ control level until 2010 and 2029, respectively, after which they slowly decrease to zero. For control $u_3$, it decreases from the maximum control level to $70.7\mathrm{\%}$, then starts to gradually increase to $95\mathrm{\%}$ in October 2002. After 43 months, control $u_3$ begins to decrease. Figure 19(c) shows that when $n_3=12$, the controls $u_1$, $u_2$, and $u_3$ reach the $100\mathrm{\%}$ control level and last for about four, five, and seven months, respectively. Then, the controls $u_1$ and $u_2$ maintain $95\mathrm{\%}$ control level until 2011 and 2029, respectively, after which they slowly decrease to zero. After the maximum control level, the control $u_1$ starts to decrease slowly after nine years of implementing the control level of $95\mathrm{\%}$, the control $u_2$ maintains the control level of $95\mathrm{\%}$ all the time, while the control $u_3$ first decreases in July 2002 and increases in September 2003, then gradually decreases after November 2005. At last, we explore the effect of different $n_3$ values on the number of infected individuals as shown in Figure 19(d). We find that the smaller the weight $n_3$, the better control effect.

Figure 19. Impact of changes on the cost weight of control $u_3\left(t\right)$, $n_3$, on control strategies and the number of infected individuals. (a) Weight: $n_1=45,n_2=20,n_3=4$. (b) Weight: $n_1=45,n_2=20,n_3=8$. (c) Weight: $n_1=45,n_2=20,n_3=12$. (d) The number of infected individuals under different cost weights of $n_3$.

In summary, as the weight $n_1$ increases, the control level of $u_1$ gradually decreases. At the same time, the control levels of $u_2$ and $u_3$ also decrease. As the weight $n_2$ increases, the control level of $u_2$ gradually decreases, and the control levels of $u_1$ and $u_3$ decrease. As the weight $n_3$ increases, the control level of $u_3$ gradually decreases, and the control levels of $u_1$ and $u_2$ decrease.

3.3. Efficiency analysis and cost-effectiveness analysis

Although all four control strategies can achieve the control objectives proposed by the UNAIDS, we also need to consider the corresponding efficiency and costs. Therefore, we explore the efficiency and cost-effectiveness of four control strategies in this section.

First, we perform the efficacy analysis using the method in references [12,27,47]. The efficiency index is defined as $\mathrm{\Sigma }=(1-\frac{{\mathrm{\Omega }}_c}{{\mathrm{\Omega }}_s})\times 100\mathrm{\%},$where ${\mathrm{\Omega }}_c$ is the cumulative number of infected individuals when different control strategies are implemented and ${\mathrm{\Omega }}_s$ is the cumulative number of infected individuals in absence of any control strategies. The cumulative number of infected individuals during the time interval $[0,t_f]$ is denoted by $\mathrm{\Omega }={\int }_0^{t_f}\left[E\right(t)+I_1(t)+I_2(t\left)\right]\mathrm{d}t,$where $t_f$ indicates the end time of the implementation controls. In this subsection, $t_f$ equals to 29 years.

The best strategy has the biggest efficiency index [1,5]. It can be seen from Table 3 that the effective indexes of four control strategies. The effective indexes show that Strategies 1, 2, 3, and 4 are effective to mitigate HIV transmission, but the most effective control strategy is Strategy 4. Next, we explore the costs of four control strategies.

Table 3. Indices for the effectivenesses of control strategies.

Strategies	$\mathrm{\Omega }={\int }_0^{t_f}E\left(t\right)+I_1\left(t\right)+I_2\left(t\right)\mathrm{d}t$	$\mathrm{\Sigma }=(1-\frac{{\mathrm{\Omega }}_c}{{\mathrm{\Omega }}_s})\times 100\mathrm{\%}$
No control	45,905,000	0
Strategy 1	3,553,479	$92.26\mathrm{\%}$
Strategy 2	7,167,623	$84.39\mathrm{\%}$
Strategy 3	7,041,835	$84.66\mathrm{\%}$
Strategy 4	1,902,933	$95.85\mathrm{\%}$

Table 4. Comparison of impacts of four strategies by 2020.

Strategies	Total infectious individuals	Total infectious averted	Total costs	Incidence reduction
No control	45,905,000	–	–	–
Strategy 1	3,553,479	42,351,521	8408	92.26%
Strategy 2	7,167,623	38,737,377	14,351	84.39%
Strategy 3	7,041,835	38,863,165	36,270	84.66%
Strategy 4	1,902,933	44,002,067	19,735	95.85%

The total costs is $\mathrm{C}\mathrm{o}\mathrm{s}\mathrm{t}\mathrm{s}={\int }_0^{t_f}\left[\frac{n_1}{2}{u}_1^2\right(t)+\frac{n_2}{2}{u}_2^2(t)+\frac{n_3}{2}{u}_3^2(t\left)\right]\mathrm{d}t$, where $n_1=45$, $n_2=25$, $n_3=10$, and $t_f=29$. The total costs of different control strategies as shown in Table 4. The result shows that Strategy 1 is the most cost-effective. However, the most expensive strategy is Strategy 3, since its cost is 4.3 times as that of Strategy 1. This may be because the control of screening latent individuals is not included in Strategy 3, it takes a long time to implement the control to achieve the control goal. Since the strategies that include control of screening latent individuals can end the epidemic in a relatively short time, the cost is relatively low.

4. Discussion and conclusion

This work analyses the spread of HIV via a mathematical model incorporating susceptible individuals, latent individuals, unaware infected individuals, aware infected individuals, infected individuals receiving ART treatment, and AIDS individuals. Based on the model, we carry out theoretical analysis and numerical simulations. The theoretical analysis shows that when $R_0<1$, the disease free equilibrium is globally asymptotically stable, and the disease gradually dies out. When $R_0>1$, the endemic equilibrium is globally asymptotically stable, and the disease spreads. We obtain the parameter values by the MCMC method. Through numerical simulations, we verify the stability of the disease free equilibrium and the endemic equilibrium. We further analyse the effects of a single measure and multiple measures on disease transmission. Compared with a single measure, multiple measures can control the transmission of disease more effectively.

We incorporate mitigation strategies, including treating aware infected individuals, screening latent individuals, and educating unaware individuals into the model. We apply Pontryagin's Maximum Principle to derive optimal control strategies. The results show that the optimal control strategy is the combination of all three control measures. Since people who receive ART treatments and aware infected individuals are not likely to transmit the infection, the epidemic can be brought under control when more infected individuals are aware of their infection and more people receive ART treatments. When all three control measures are implemented, AIDS can be end in the shortest time. In our study, we find that each of the four different control strategies in this work can achieve the three $90\mathrm{\%}$ reduction goals and end the AIDS epidemic by 2030 proposed by the UNAIDS in 2014. Simultaneously, the earlier the control strategies are implemented, the sooner AIDS can be eradicated. Particularly, we find that the proportion of aware infected individuals has an important influence on the transmission of HIV. The greater the proportion of humans who are aware infected, the easier it is to control the spread of HIV. This is because ART therapy can achieve a great reduction in HIV incidence, while more efficient HIV testing is needed [37].

Finally, we analyse the impact of cost changes according to the optimal control strategy. The results show that the number of infected individuals can be greatly reduced when the cost of implementing all control strategies is lower. When more medicines for HIV are covered by medical insurances, more infected individuals can receive treatments [29]. Moreover, we analyse the total costs of control strategies and find that the most cost-effective strategy is a combination of enhanced screening for latent individuals and education for unaware infected. The reason is that the combined implementation of these two control strategies can avoid more infected individuals, which in turn can save the costs for testing and treatment [37]. Therefore, we can mitigate the spread of AIDS by improving medical technology, increasing testing rates, and intensify education for unaware infected individuals.

The duration of AIDS is long, and it is difficult to obtain the number of infected individuals at each stage. Once such data are available, we can evaluate whether the goal of eradicating AIDS can be achieved by 2030 more accurately. Moreover, the model can be applied to study many other diseases with similar transmission mechanisms, such as COVID-19, Syphilis, and hepatitis B. We did not take into account detection rate and the effect of HIV co-infection with other diseases on HIV transmission in this work. In future work, we will study the effects of the detection rate and the co-infection of HIV with other diseases, and discuss whether the eradication of AIDS can be achieved earlier once an effective vaccine is available.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

LX is funded by the National Natural Science Foundation of China [grant number 12171116] and Fundamental Research Funds for the Central Universities of China [grant numbers 3072020CFT2402 and 3072022TS2404]. WS is funded by the Fundamental Research Funds for the Central Universities of China [grant number 3072021CFP2401].