∎

11institutetext: J. P. Romero-Leiton & Julien Arino 22institutetext:
Department of Mathematics, University of Manitoba,
Winnipeg MB R3T 1E9 Canada.
22email<EMAIL_ADDRESS>&<EMAIL_ADDRESS>33institutetext: I. Sekkak & B. Nasri 44institutetext: Département de Médecine
Sociale et Préventive,
École de Santé Publique de l’Université de Montréal, Montreal QC H3N 1X9
Canada.
44email<EMAIL_ADDRESS>&<EMAIL_ADDRESS>
# Mathematical modelling of the first HIV/ZIKV co-infection cases in Colombia
and Brazil

Jhoana P. Romero-Leiton  Idriss Sekkak  Julien Arino  Bouchra Nasri

(Received: date / Accepted: date)

###### Abstract

This paper presents a mathematical model to investigate co-infection with
HIV/AIDS and zika virus (ZIKV) in Colombia and Brazil, where the first cases
were reported in 2015-2016. The model considers the sexual transmission
dynamics of both viruses and vector-host interactions. We begin by exploring
the qualitative behaviour of each model separately. Then, we analyze the
dynamics of the co-infection model using the thresholds and results defined
separately for each model. The model also considers the impact of intervention
strategies, such as, personal protection, antiretroviral therapy (ART), and
sexual protection (condoms use). Using available parameter values for Colombia
and Brazil, the model is calibrated to predict the potential effect of
implementing combinations of those intervention strategies on the co-infection
spread. According to these findings, transmission through sexual contact is a
determining factor in the long-term behaviour of these two diseases.
Furthermore, it is important to note that co-infection with HIV and ZIKV may
result in higher rates of HIV transmission and an increased risk of severe
congenital disabilities linked to ZIKV infection. As a result, control
measures have been implemented to limit the number of infected individuals and
mosquitoes, with the aim of halting disease transmission. This study provides
novel insights into the dynamics of HIV/ZIKV co-infection and highlights the
importance of integrated intervention strategies in controlling the spread of
these viruses, which may impact public health.

###### Keywords:

Stability Equilibrium points Optimal control Personal protection Sexual
protection Antiretroviral therapyModel calibration.

##  1 Introduction

Human immunodeficiency virus (HIV) and zika virus (ZIKV) are two major public
health concerns worldwide, particularly in Latin America and Caribbean
countries villamil2018 . While HIV is a chronic infection that attacks the
immune system, ZIKV is transmitted by mosquitoes and can even cause congenital
malformations, such as Guillain-Barré syndrome cao2016 ; oehler2014 ;
smith2016  and microcephaly mlakar2016 ; calvet2016 . On the one hand, if HIV
is not treated, it can cause Acquired Immunodeficiency Syndrome (AIDS). This
virus can be transmitted through sexual contact, syringe misuse, and
vertically (from mother to child) kabapy2020attributes . HIV/AIDS still has no
cure, so treatments seek to lower or reduce the level of virus replication in
the body of infected people, which consists of a mixture of different drugs,
commonly called antiretroviral (or combined) therapy (ART) tilton2010entry .
On the other hand, ZIKV, unlike other arboviruses, also presents transmission
through sexual contact. Some studies have demonstrated its detection and
transmission through semen, urine, and saliva atkinson2016 ; gourinat2015 ;
musso2015 .

Thus, although HIV is not a zoonotic disease, it shares specific transmission
characteristics with ZIKV. Both viruses can be transmitted through sexual
contact and vertically from mother to fetus. These additional modes of
transmission can significantly increase the risk of other sexually transmitted
infections in endemic regions, particularly HIV, for both the mother and the
fetus rothan2018 . There are few case reports of ZIKV infection in HIV-
infected individuals sherman2018zika . The first recorded case of HIV/ZIKV co-
infection was confirmed in a 38-year-old patient in a Rio de Janeiro (Brazil)
laboratory in 2016 calvet2016 . In the same region, a zika case was reported
in an HIV-infected pregnant woman joao2018 ; zorrilla2016 . The fetus
exhibited significant abnormalities, consistent with research on pregnant
women in Brazil who were infected with ZIKV brasil2016 , and the fetus
ultimately died. In 2018, five individuals from the departments of Risaralda
and Sucre were informed of HIV/ZIKV co-infection in Colombia villamil2018 .
Fortunately, these patients had reasonable immune and virologic controls when
infected with ZIKV, with no relevant differences from patients infected with
ZIKV alone.

However, additional research is necessary to understand better the
interactions between HIV/ZIKV co-infection and their impact on immune
response, disease severity, and control villamil2018 . The question remains
unanswered as to whether HIV infection increases the likelihood of contracting
ZIKV and whether ZIKV infection could worsen HIV infection, particularly
during pregnancy. Nevertheless, laboratory studies have revealed that
placental tissues are vulnerable to ZIKV infection rothan2018 . In addition to
causing placental dysfunction in co-infections, ZIKV infection significantly
affects the host’s immune response. A recent study indicated that ZIKV
primarily targets CD14+ monocytes, particularly during pregnancy, resulting in
inflammatory reactions and immune tolerance rothan2018 . However,ZIKV has the
potential to facilitate HIV infection by promoting HIV replication through the
release of cytokines that can activate CD4+ T cells or by binding directly to
HIV proteins that encourage HIV replication. Unlike HIV, the relationship
between mother-to-child transmission of ZIKV and fetal disease caused by ZIKV
infection has not yet been established. However, even in the presence of ART,
acute viral infection is likely to intensify immune system dysfunction in HIV-
infected pregnant women, increasing the risk and consequences of mother-to-
child transmission of HIV and ZIKV aschengrau2021international . Therefore,
the potential interaction between HIV and ZIKV has recently garnered
significant attention. These interactions can modify the epidemiology,
pathogenesis, immune response, and therapy of both infections. For instance,
co-infection can expedite HIV pathogenesis and enhance transmission by
boosting viral replication efficiency. Additionally, the spread of ZIKV
infection may reveal new and more severe clinical manifestations in
immunocompromised individuals exposed to the virus.

Given the potential impact of HIV/ZIKV co-infection on public health, it is
crucial to understand the transmission dynamics of these viruses and evaluate
the effectiveness of intervention strategies. Mathematical models are crucial
for understanding and providing valuable insights into public health policy
decisions. To our knowledge, there is no evidence of mathematical models
studying this phenomenon in the literature. Therefore, this study aimed to
formulate and analyze an HIV/ZIKV co-infection model, assuming that both
diseases are sexually transmitted, and ZIKV is also mosquito-transmitted. The
analysis of this model is expected to identify important transmission factors
that will help design and evaluate different control and prevention strategies
to minimize their impact on public health.

The organization of this study is as follows: In Section 2, the co-infection
model is introduced, providing a comprehensive understanding of its structure
and dynamics. Subsequently, Sections 2.1 and 2.2 present the individual
dynamics of the HIV-only and ZIKV-only models, respectively, highlighting the
unique characteristics of each infection. Section 3 focuses on the analysis of
the co-infection model, delving into the interplay and synergistic effects
between HIV and ZIKV. The optimal control problem is addressed analytically in
Section 4, emphasizing the exploration of strategies to mitigate the spread
and impact of the co-infection. Furthermore, in Section 5, a case study
centred in Colombia and Brazil is presented, whereby the uncontrolled and
controlled models are numerically analyzed using data derived from the
comprehensive literature review. Finally, in Section 6, a discussion and
concluding remarks are provided, encompassing the modelling approach, its
results, along with prospects for future research and unresolved questions.

##  2 The HIV/ZIKV mathematical model formulation

This model examines two distinct groups: the human (host) and mosquito
(vector) populations. We assumed that an individual who is susceptible to
either disease cannot be infected with both diseases simultaneously at the
same time. The total number $N(t)$ of people in the human population at a
given time is divided into six categories: those who are susceptible to both
viruses, denoted as $S(t)$; those who are infected with only ZIKV but still
susceptible to HIV, denoted as $I_{z}(t)$; those who are infected with only
HIV but still susceptible to ZIKV, denoted as $I_{h}(t)$; those who are
infected with both HIV and ZIKV simultaneously, denoted as $I_{hz}(t)$; those
who are infected with AIDS, denoted as $A(t)$; and those who have recovered
from ZIKV, denoted as $R(t)$. Thus, the total human population can be
represented as the sum of all these categories
($N(t)=S(t)+I_{z}(t)+I_{h}(t)+I_{hz}(t)+A(t)+R(t)$), as shown in Tables
2.1-2.2.

The total mosquito population $N_{m}(t)$ at time $t$ can be classified into
two compartments: the susceptible mosquito population $S_{m}(t)$ and the ZIKV-
carrying mosquitoes $I_{m}(t)$. Thus, $N_{m}(t)=S_{m}(t)+I_{m}(t)$ (see Tables
2.1 and 2.2).

Variable | Description
---|---
$N(t)$ |  The total human population at time $t$
$S(t)$ |  Susceptible human population at time $t$
$I_{z}(t)$ |  Infected human population with only ZIKV at time $t$
$I_{h}(t)$ |  Infected human population with only HIV at time $t$
$I_{hz}(t)$ |  Infected human population with ZIK/HIV at time $t$
$A(t)$ |  Infected human population with AIDS at time $t$
$R(t)$ |  Recovered human population of ZIKV at time $t$
$N_{m}(t)$ |  The total mosquito population at time $t$
$S_{m}(t)$ |  Susceptible mosquito population at time $t$
$I_{m}(t)$ |  ZIKV-carrying mosquito population at time $t$
Table 2.1: Description of the state variables involved in Model (2.1).

Parameter | Description | Dimension
---|---|---
$\Lambda$ | Recruitment rate of humans | $pop\times time^{-1}$
$\beta_{m}$ | Infection rate of humans by contact with infected mosquitoes with ZIKV | $(pop\times time)^{-1}$
$\beta_{z}$ | Infection rate of humans by contact with humans infected with ZIKV through sexual contact | $(pop\times time)^{-1}$
$\beta_{h}$ | Infection rate of humans by contact with humans infected with HIV through sexual contact | $(pop\times time)^{-1}$
$1/\sigma_{1}$ | Mean duration of the immunodeficiency period | time
$1/\sigma_{2}$ | Mean duration of the immunodeficiency period in co-infected individuals |  time-1
$1/\mu$ | Human mean lifespan | time
$\mu_{z}$ | Mortality rate by zika |  time-1
$\mu_{h}$ | Mortality rate by AIDS |  time-1
$\mu_{hz}$ | Mortality rate by ZIKV/AIDS |  time-1
$1/\delta_{z}$ | Mean duration of the zika infection | time
$\omega_{1}$ | Transition probability from zika to HIV/ZIKV co-infection | Dimensionless
$\omega_{2}$ | Transition probability from HIV to HIV/ZIKV co-infection | Dimensionless
$\epsilon$ | Modification parameter or decrease recovery rate factor | Dimensionless
$\Lambda_{m}$ | Recruitment rate of mosquitoes | $pop\times time^{-1}$
$\alpha_{m}$ | Infection rate of mosquitoes by contact with infected humans with ZIKV | $(pop\times time)^{-1}$
$1/\mu_{m}$ | Mosquito mean lifespan | time

Table 2.2: Description and dimension of the parameters involved in Model (2.1.

$\left\\{\begin{array}[]{ll}\dfrac{dS}{dt}&=\Lambda-\left(\beta_{m}I_{m}+\beta_{z}I_{z}+\beta_{h}I_{h}\right)S-\mu
S\\\ \\\
\dfrac{dI_{z}}{dt}&=\left(\beta_{m}I_{m}+\beta_{z}I_{z}\right)S-\omega_{2}\beta_{h}I_{z}I_{h}-(\mu_{z}+\delta_{z}+\mu)I_{z}\\\
\\\
\dfrac{dI_{h}}{dt}&=\beta_{h}I_{h}S-\omega_{1}(\beta_{m}I_{m}+\beta_{z}I_{z})I_{h}-(\sigma_{1}+\mu)I_{h}\\\
\\\
\dfrac{dI_{hz}}{dt}&=\omega_{1}(\beta_{m}I_{m}+\beta_{z}I_{z})I_{h}+\omega_{2}\beta_{h}I_{h}I_{z}-\epsilon\delta_{z}I_{hz}-(\sigma_{2}+\mu_{hz}+\mu)I_{hz}\\\
\\\ \dfrac{dA}{dt}&=\sigma_{1}I_{h}+\sigma_{2}I_{hz}-(\mu_{h}+\mu)A\\\ \\\
\dfrac{dR}{dt}&=\delta_{z}I_{z}+\epsilon\delta_{z}I_{hz}-\mu R\\\ \\\
\dfrac{dS_{m}}{dt}&=\Lambda_{m}-\alpha_{m}(I_{z}+I_{hz})S_{m}-\mu_{m}S_{m}\\\
\\\ \dfrac{dI_{m}}{dt}&=\alpha_{m}(I_{z}+I_{hz})S_{m}-\mu_{m}I_{m}\\\ \\\
\end{array}\right.$ (2.1)

Figure 2.1 illustrates the dynamics depicted by the equations involved in
Model (2.1)..

Figure 2.1: HIV/ZIKV co-infection model represented in Model (2.1).

In the following two sections, we qualitatively analyse the properties of
System (2.1). We will start by analysing the dynamics of the two component
models: the HIV/AIDS model and the ZIKV model.

###  2.1 Qualitative behaviour of the HIV/AIDS model

The HIV/AIDS model, is obtained by setting $I_{z}=I_{hz}=R=S_{m}=I_{m}=0$ in
System (2.1). Thus, the ODEs described in (2.1) can be rewritten as:

$\left\\{\begin{array}[]{ll}\dfrac{dS}{dt}&=\Lambda-\beta_{h}I_{h}S-\mu S\\\
\\\ \dfrac{dI_{h}}{dt}&=\beta_{h}I_{h}S-(\sigma_{1}+\mu)I_{h}\\\ \\\
\dfrac{dA}{dt}&=\sigma_{1}I_{h}-(\mu_{h}+\mu)A,\\\ \\\ \end{array}\right.$
(2.2)

where the total human population is $N_{h}(t)=S(t)+I_{h}(t)+A(t)$. For this
model, the region of biological interest is

$\Omega_{h}=\left\\{(S,I_{h},A)\in\mathbb{R}_{+}^{3}:0\leq
N_{h}\leq\frac{\Lambda}{\mu}\right\\}.$ (2.3)

It can be proved that $\Omega_{h}$ is positively-invariant under the flow of
(2.2) (see e.g., romero2019optimal ), that is, all solutions of System (2.2)
starting in $\Omega_{h}$ remain in $\Omega_{h}$ for all $t\geq 0$. Therefore,
it is enough to consider the dynamics of (2.2) in $\Omega_{h}$.

The disease-free equilibrium (DFE) of Model (2.2) is given by

$\textbf{E}_{h_{0}}=\left(\frac{\Lambda}{\mu},0,0\right),$ (2.4)

which is obtained when $I_{h}=A=0$. The stability of this equilibrium point
can be analyzed in terms of the basic reproduction number for the HIV/AIDS
model ($\mathcal{R}_{h}$), which can be computed using the next-generation
operator van2002reproduction . Using the notation of romero2019optimal  in
Model (2.2) the matrices F and V are given by

$\textbf{F}=\begin{bmatrix}\beta_{h}\frac{\Lambda}{\mu}&0\\\
0&0\end{bmatrix},\quad\text{and}\quad\textbf{V}=\begin{bmatrix}\sigma_{1}+\mu&0\\\
-\sigma_{1}&\mu_{h}+\mu\end{bmatrix}.$

Therefore, it follows that $\mathcal{R}_{h}$ associated to Model (2.2) is
given by

$\mathcal{R}_{h}=\rho(\textbf{FV}^{-1})=\dfrac{\beta_{h}}{\sigma_{1}+\mu}\dfrac{\Lambda}{\mu},$
(2.5)

where $\rho$ represents the spectral radius of the matrix $\textbf{FV}^{-1}$.

To determine the endemic equilibrium points of Model(2.2), we must solve the
system of algebraic equations

$\left\\{\begin{array}[]{ll}&0=\Lambda-\beta_{h}I_{h}S-\mu S\\\ \\\
&0=\beta_{h}I_{h}S-(\sigma_{1}+\mu)I_{h}\\\ \\\
&0=\sigma_{1}I_{h}-(\mu_{h}+\mu)A.\end{array}\right.$ (2.6)

After some algebraic manipulations for $I_{h},A\neq 0$, we find that the
solutions of System (2.6) are

$S^{*}=\frac{\sigma_{1}+\mu}{\beta_{h}},\quad
I_{h}^{*}=\frac{\mu}{\beta_{h}}(\mathcal{R}_{h}-1)\quad\text{and}\quad
A^{*}=\frac{\mu\sigma_{1}}{\beta_{h}(\mu_{h}+\mu)}(\mathcal{R}_{h}-1).$

Thus, as long as $\mathcal{R}_{h}>1$, System (2.2) has an endemic equilibrium
point given by

$\textbf{E}^{*}_{h}=\left(\frac{\sigma_{1}+\mu}{\beta_{h}},\frac{\mu}{\beta_{h}}(\mathcal{R}_{h}-1),\frac{\mu\sigma_{1}}{\beta_{h}(\mu_{h}+\mu)}(\mathcal{R}_{h}-1)\right).$
(2.7)

The local stability of $\textbf{E}_{h_{0}}$ defined on (2.4) and
$\textbf{E}^{*}_{h}$ defined on (2.7) is determined by the sign of the
eigenvalues of the linearisation matrix (Jacobian matrix) of System (2.2)
around them. The Jacobian matrix of System (2.2) at an arbitrary point
$\textbf{E}=(S,I_{h},A)$ is

$\textbf{J}(\textbf{E})=\begin{bmatrix}-(\beta_{h}I_{h}+\mu)&-\beta_{h}S&0\\\
\beta_{h}I_{h}&\beta_{h}S-(\sigma_{1}+\mu)&0\\\
0&\sigma_{1}&-(\mu_{h}+\mu)\end{bmatrix}=\begin{bmatrix}\mathbf{J}_{11}(\mathbf{E})&\mathbf{0}\\\
\mathbf{\star}&-(\mu_{h}+\mu)\end{bmatrix},$ (2.8)

so eigenvalues of $\mathbf{J}(\mathbf{E})$ are $-(\mu_{h}+\mu)<0$ and those of
$\mathbf{J}_{11}(\mathbf{E})$. At $\mathbf{E}_{h_{0}}$,

$\mathbf{J}_{11}(\mathbf{E}_{h_{0}})=\begin{bmatrix}-\mu&-\beta_{h}\frac{\Lambda}{\mu}\\\
0&(\sigma_{1}+\mu)(\mathcal{R}_{h}-1)\end{bmatrix}$

is upper triangular with eigenvalues $-\mu<0$ and
$(\sigma_{1}+\mu)(\mathcal{R}_{h}-1)$. At $\textbf{E}^{*}_{h}$,

$\mathbf{J}_{11}(\mathbf{E}_{h}^{*})=\begin{bmatrix}-\mu\mathcal{R}_{h}&-(\sigma_{1}+\mu)\\\
\mu\mathcal{R}_{h}&0\end{bmatrix}.$

Thus $\mathbf{J}_{11}(\mathbf{E}_{h}^{*})$ has positive determinant and
negative trace, implying eigenvalues with negative real parts when
$\mathbf{E}_{h}^{*}$ is biologically relevant, i.e., when $\mathcal{R}_{h}>1$.
These results are summarized in the following proposition.

######  Proposition 1.

System (2.2) always has a DFE $\textbf{E}_{h_{0}}$ given in (2.4), and for
$\mathcal{R}_{h}>1$ defined in (2.5), also exists an endemic equilibrium point
$\textbf{E}^{*}_{h}$ given in (2.7). Additionally, if $\mathcal{R}_{h}<1$ the
DFE is locally asymptotically stable (LAS) in $\Omega_{h}$ defined in (2.3),
whereas $\textbf{E}^{*}_{h}$ is unestable. If $\mathcal{R}_{h}>1$, then
$\textbf{E}^{*}_{h}$ is LAS in $\Omega_{h}$ and the DFE becomes an unstable
hyperbolic point.

As it is well known, the basic reproduction number of an infection
$\mathcal{R}_{h}$ is the average number of new cases generated by a given case
throughout an infectious period. Thus, Proposition 1 tells us that HIV
infection can be controlled in the community if the initial values of the
subpopulation of the model are in the region of attraction of
$\textbf{E}_{h_{0}}$. Therefore, to ensure that the control of the virus is
independent of initial conditions, we must prove the global asymptotic
stability (GAS) of the equilibrium points. We obtain the following result.

######  Proposition 2.

If $\mathcal{R}h<1$, the DFE $\textbf{E}{h_{0}}$ of Model (2.2), is GAS. If
$\mathcal{R}_{h}>1$, the endemic equilibrium point $\textbf{E}^{*}_{h}$, is
GAS.

The proof of global asymptotic stability of the DFE when $\mathcal{R}_{h}<1$
is found in Appendix A.1. The global asymptotic stability of
$\textbf{E}^{*}_{h}$ when $\mathcal{R}_{h}>1$ is found in Appendix A.2.

###  2.2 Qualitative behaviour of the ZIKV model

By setting $I_{h}=I_{hz}=A=0$ in System (2.1), we obtain the ZIKV model as
follows

$\left\\{\begin{array}[]{ll}\dfrac{dS}{dt}&=\Lambda-\left(\beta_{m}I_{m}+\beta_{z}I_{z}\right)S-\mu
S\\\ \\\
\dfrac{dI_{z}}{dt}&=\left(\beta_{m}I_{m}+\beta_{z}I_{z}\right)S-(\mu_{z}+\delta_{z}+\mu)I_{z}\\\
\\\ \dfrac{dR}{dt}&=\delta_{z}I_{z}-\mu R\\\ \\\
\dfrac{dS_{m}}{dt}&=\Lambda_{m}-\alpha_{m}I_{z}S_{m}-\mu_{m}S_{m}\\\ \\\
\dfrac{dI_{m}}{dt}&=\alpha_{m}I_{z}S_{m}-\mu_{m}I_{m}\end{array}\right.$ (2.9)

For this model, the total human population is $N_{z}(t)=S(t)+I_{z}(t)+R(t)$,
and the total mosquito population is $N_{m}(t)=S_{m}(t)+I_{m}(t)$. We are
interested in analyzing the solutions within the biological region of interest

$\Omega_{z}=\left\\{(S,I_{z},R,S_{m},I_{m})\in\mathbb{R}_{+}^{5}:0\leq
N_{h}\leq\frac{\Lambda}{\mu};0\leq
N_{m}\leq\frac{\Lambda_{m}}{\mu_{m}}\right\\}.$ (2.10)

It can also be shown that $\Omega_{z}$ is positively invariant under the flow
of (2.9).

The DFE for Model (2.9) is given by

$\textbf{E}_{z_{0}}=\left(\frac{\Lambda}{\mu},0,0,\frac{\Lambda_{m}}{\mu_{m}},0\right).$
(2.11)

Now, we will determine the basic reproduction number $\mathcal{R}_{z}$ for
Model (2.9). Similarly as in Section 2.1, the matrices F and V are given by

$\textbf{F}=\begin{bmatrix}\beta_{z}\frac{\Lambda}{\mu}&\beta_{m}\frac{\Lambda}{\mu}\\\
\alpha_{m}\frac{\Lambda_{m}}{\mu_{m}}&0\\\
\end{bmatrix},\quad\text{and}\quad\textbf{V}=\begin{bmatrix}\delta_{z}+\mu+\mu_{z}&0\\\
0&\mu_{m}\end{bmatrix}.$

Thus, the basic reproduction number for Model (2.9) is given by

$\mathcal{R}_{z}=\rho(\textbf{FV}^{-}1)=\mathcal{R}_{z_{1}}+\sqrt{\left(\mathcal{R}_{z_{1}}\right)^{2}+\bar{\mathcal{R}}_{z_{2}}}=:\mathcal{R}_{z_{1}}+\mathcal{R}_{z_{2}},$
(2.12)

where

$\mathcal{R}_{z_{1}}=\dfrac{\beta_{z}\Lambda}{2\mu(\mu_{z}+\delta_{z}+\mu)},\quad\bar{\mathcal{R}}_{z_{2}}=\dfrac{\beta_{m}\alpha_{m}\Lambda_{m}\Lambda}{\mu\mu_{m}^{2}(\mu_{z}+\delta_{z}+\mu)}\quad\text{and}\quad\mathcal{R}_{z_{2}}=\sqrt{\left(\mathcal{R}_{z_{1}}\right)^{2}+\bar{\mathcal{R}}_{z_{2}}}.$
(2.13)

######  Remark 1.

Note that if the transmission of ZIKV by sexual contact is not considered
($\beta_{z}=0$), $\mathcal{R}_{z_{1}}=0$ and $\mathcal{R}_{z}$ reduces to

$\left.\mathcal{R}_{z}\right|_{(\beta_{z}=0)}=\sqrt{\bar{\mathcal{R}}_{z_{2}}}=\sqrt{\dfrac{\beta_{m}\alpha_{m}\Lambda_{m}\Lambda}{\mu\mu_{m}^{2}(\mu_{z}+\delta_{z}+\mu)}},$

indicating that sexual contact transmission of ZIKV has an impact on
$\mathcal{R}_{z}$.

The following lemma makes it easier to determine the sign of
$\mathcal{R}_{z}$.

######  Lemma 1.

Let us define

$\begin{array}[]{ll}\mathcal{R}_{z}^{*}&=2\mathcal{R}_{z_{1}}+\bar{\mathcal{R}}_{z_{2}}\\\
\\\
&=:\mathcal{R}_{z}^{2}+2\mathcal{R}_{z_{1}}(1-\mathcal{R}_{z}).\end{array}$
(2.14)

  1. i.

If $\mathcal{R}_{z}^{*}<1$, then $2\mathcal{R}_{z_{1}}<1$ and
$\mathcal{R}_{z}<1$.

  2. ii.

If $\mathcal{R}_{z}^{*}>1$ and $2\mathcal{R}_{z_{1}}<1$, then
$\mathcal{R}_{z}>1$.

  3. iii.

If $\mathcal{R}_{z}^{*}>1$ and $2\mathcal{R}_{z_{1}}>1$, then
$\mathcal{R}_{z}<1$.

  4. iv.

If $\mathcal{R}_{z}^{*}=1$, then $\mathcal{R}_{z}=1$.

The proof can be found in Appendix A.3.

From the above lemma we can see that the sign of $\mathcal{R}_{z}$ is
determined by the signs of $\mathcal{R}^{*}_{z}$ and $2\mathcal{R}_{z_{1}}$.

We can now determine the local asymptotic stability of the DFE
$\textbf{E}_{z_{0}}$. To this end, we order equations and variables as
$S,S_{m},I_{z},I_{m},R$ and compute the Jacobian matrix of System (2.9) at an
arbitrary point $\textbf{E}=(S,S_{m},I_{z},I_{m},R)$, which is given by

$\textbf{J}(\textbf{E})=\begin{bmatrix}\textbf{J}_{11}(\textbf{E})&\mathbf{0}\\\
\mathbf{\star}&-\mu,\end{bmatrix}.$ (2.15)

Thus, the eigenvalues of J(E) are $-\mu<0$ and those of
$\textbf{J}_{11}(\textbf{E})$. At $\textbf{E}_{z_{0}}$

$\textbf{J}_{11}(\textbf{E}_{z_{0}})=\begin{bmatrix}-\mu&0&-\beta_{z}\Lambda/\mu&-\beta_{m}\Lambda/\mu\\\
0&-\mu_{m}&-\alpha_{m}\Lambda_{m}/\mu_{m}&0\\\
0&0&\beta_{z}\Lambda/\mu-(\mu_{z}+\delta_{z}+\mu)&\beta_{m}\Lambda_{m}/\mu_{m}\\\
0&0&\alpha_{m}\Lambda_{m}\mu_{m}&-\mu_{m}\end{bmatrix}=\begin{bmatrix}\mathbf{D}_{11}&\mathbf{\star}\\\
0&\mathbf{D}_{22}\end{bmatrix}.$

Again, the above matrix is $2\times 2$ block upper triangular with (1,1) block
being a diagonal matrix with diagonal entries (eigenvalues) $-\mu<0$ and
$-\mu_{m}<0$. So all that remains to determine the eigenvalues of the $2\times
2$ matrix

$\mathbf{D}_{22}=\begin{bmatrix}\beta_{z}\dfrac{\Lambda}{\mu}-(\mu_{z}+\delta_{z}+\mu)&\beta_{m}\dfrac{\Lambda_{m}}{\mu_{m}}\\\
\\\ \alpha_{m}\dfrac{\Lambda_{m}}{\mu_{m}}&-\mu_{m}\end{bmatrix}.$

After some algebraic manipulations, we found that the above matrix has
associated the following characteristic polynomial.

$q(x)=x^{2}+a_{1}x+a_{2},\quad\text{where}$ (2.16)

$\begin{array}[]{ll}&a_{1}=\mu_{m}-(\mu_{z}+\delta_{z}+\mu)(2\mathcal{R}_{z_{1}}-1)\\\
\\\ &a_{2}=\mu_{m}(\mu_{z}+\delta_{z}+\mu)(1-R_{z}^{*}).\end{array}$

Note that all roots of the characteristic polynomial $q(x)$ have a negative
real part if $\mathcal{R}^{*}_{z}$ defined in (2.14) satisfies
$\mathcal{R}^{*}_{z}<1$ (see Lemma 1). We have the following result.

######  Proposition 3.

If $\mathcal{R}^{*}_{z}<1$, then $\textbf{E}_{z_{0}}$ defined in (2.11) is LAS
in $\Omega_{z}$ defined in (2.10).

Now, in order to determine the endemic equilibrium points of System (2.9), we
have to solve the following system of algebraic equations

$\left\\{\begin{array}[]{ll}0=\Lambda-\left(\beta_{m}I_{m}+\beta_{z}I_{z}\right)S-\mu
S\\\ \\\
0=\left(\beta_{m}I_{m}+\beta_{z}I_{z}\right)S-(\mu_{z}+\delta_{z}+\mu)I_{z}\\\
\\\ 0=\delta_{z}I_{z}-\mu R\\\ \\\
0=\Lambda_{m}-\alpha_{m}I_{z}S_{m}-\mu_{m}S_{m}\\\ \\\
0=\alpha_{m}I_{z}S_{m}-\mu_{m}I_{m}.\end{array}\right.$ (2.17)

To this end, we first define the following threshold

$\begin{array}[]{ll}&R_{max}=\dfrac{\Lambda}{\mu}\dfrac{\delta_{z}}{\mu_{z}+\delta_{z}+\mu}.\end{array}$
(2.18)

Therefore, after some algebraic manipulations in (2.17), we can write the
variables $S$, $I_{z}$, $S_{m}$ and $I_{m}$ in terms of the variable $R$ and
the thresholds defined in (2.18) as follows.

$\begin{array}[]{llll}&S=\dfrac{\Lambda}{\mu}\left(1-\dfrac{R}{R_{max}}\right),&I_{z}=\dfrac{\mu}{\delta_{z}}R,\\\
\\\ &S_{m}=\dfrac{\Lambda_{m}}{\mu_{m}}\left(1-\dfrac{\alpha_{m}\mu
R}{\mu_{m}\delta_{z}+\alpha_{m}\mu
R}\right),&I_{m}=\dfrac{\alpha_{m}\Lambda_{m}\mu}{\mu_{m}\delta_{z}+\alpha_{m}\mu
R}R.\end{array}$ (2.19)

From the above expressions, it can be observed that $S_{m}$ is always positive
(because $\mu_{m}\delta_{z}+\alpha_{m}\mu R>\alpha_{m}\mu R$ and, if
$R<R_{max}$ then $S>0$). Replacing the above values into the second equation
of (2.17), we get the following quadratic equations in the variable $R$:

$aR^{2}+bR+c=0,\quad\text{where}$ (2.20)

$\begin{array}[]{ll}a&=\dfrac{\beta_{z}\Lambda\alpha_{m}\mu}{\delta_{z}R_{max}}\\\
\\\
b&=\dfrac{\mu(\mu_{z}+\delta_{z}+\mu)}{\delta_{z}}\left[\dfrac{\beta_{m}\alpha_{m}\Lambda_{m}}{\mu_{m}}+\beta_{z}\mu_{m}+\alpha_{m}\mu(1-2\mathcal{R}_{z_{1}})\right]\\\
\\\
c&=(\mu_{z}+\delta_{z}+\mu)\mu\mu_{m}\left(1-\mathcal{R}_{z}^{*}\right).\end{array}$
(2.21)

Note that $a>0$, and the signs of $b$ and $c$ depends on the sign of
$\mathcal{R}^{*}_{z}$ and $2\mathcal{R}_{z_{1}}$. We have the following
possibilities.

  * P1)

If $\mathcal{R}^{*}_{z}<1$, then $2R_{z_{1}}<1$ (Lemma 2.14 item $i$).
Therefore $b>0$ y $c>0$ and thus, the quadratic equation (2.20) has not any
positive root.

  * P2)

If $\mathcal{R}^{*}_{z}>1$, then $c<0$ and regardless of the sign of $b$, the
quadratic equation (2.20) has only one positive root given by

$R^{*}=\frac{-b+\sqrt{b^{2}-4ac}}{2a}.$ (2.22)

  * P3)

If $\mathcal{R}^{*}_{z}=1$ and $2R_{z_{1}}<1$, then $c=0$ and $b>0$.
Therefore, the quadratic equation (2.20) has as solution $R=-b/a$. Thus, there
are no positive roots.

  * P4)

If $\mathcal{R}^{*}_{z}=1$ and $b<0$, the quadratic equation (2.20) has a
positive root given by $R=-b/a$.

######  Remark 2.

The case $\mathcal{R}^{*}_{z}<1$ and $2R_{z_{1}}>1$ ($c>0$ and $b<0$), which
gives the possibility of the existence of two positive roots for the quadratic
equation (2.20) whenever $b^{2}-4ac>0$, is not considered since by Lemma 2.14
item $i$, when $\mathcal{R}^{*}_{z}<1$ leads to $2R_{z_{1}}<1$.

Based on the information presented earlier, the quadratic equation (2.20) has
only one positive root defined in (2.22) whenever $\mathcal{R}^{*}_{z}>1$. The
proof that the condition $R^{*}<R_{max}$ is satisfied can be found in Appendix
A.4. Thus, we can ensure the conditions for the existence of the endemic
equilibrium point.

######  Proposition 4.

If $\mathcal{R}_{z}^{*}>1$ System (2.9) has an endemic equilibrium point given
by

$\textbf{E}^{*}_{z}=\left(\dfrac{\Lambda}{\mu}\left(1-\dfrac{R^{*}}{R_{max}}\right),\dfrac{\mu}{\delta_{z}}R^{*},R^{*},\dfrac{\Lambda_{m}}{\mu_{m}}\left(1-\dfrac{\alpha_{m}\mu
R^{*}}{\mu_{m}\delta_{z}+\alpha_{m}\mu
R^{*}}\right),\dfrac{\alpha_{m}\Lambda_{m}\mu}{\mu_{m}\delta_{z}+\alpha_{m}\mu
R^{*}}R^{*}\right),$ (2.23)

with $R_{max}$ given in (2.18) and $R^{*}$ in (2.22).

We have already stated in Proposition 3 the local asymptotic stability of the
DFE. In a similar way, we can prove the local asymptotic stability of the
endemic equilibrium $\textbf{E}^{*}_{z}$. To this end, we determine the sign
of the eigenvalues of the matrix defined in (2.15) evaluated in
$\textbf{E}^{*}_{z}$. In order to simplify algebraic manipulations we make the
following change of variables

$\begin{array}[]{lll}&\kappa=\mu_{z}+\delta_{z}+\mu,&t_{1}=\beta_{m}I_{m}^{*}+\beta_{z}I_{z}^{*},\\\
\\\ &t_{2}=\beta_{z}S^{*},&t_{3}=\beta_{m}S^{*},\\\ \\\
&t_{4}=\alpha_{m}S_{m}^{*},&t_{5}=\alpha_{m}I_{z}^{*},\end{array}$

thus, the matrix $\textbf{J}(\textbf{E}^{*}_{z})$ can be rewritten as

$\textbf{J}(\textbf{E}^{*}_{z})=\begin{bmatrix}-(t_{1}+\mu)&-t_{2}&0&0&-t_{3}\\\
t_{1}&t_{2}-\kappa&0&0&t_{3}\\\ 0&\delta_{z}&-\mu&0&0\\\
0&-t_{4}&0&-(t_{5}+\mu_{m})&0\\\ 0&t_{4}&0&t_{5}&-\mu_{m}\end{bmatrix}.$

After some algebraic manipulations, we find that the characteristic polynomial
associated to the matrix $\textbf{J}(\textbf{E}^{*}_{z})$ is

$r(y)=(y+\mu)(y+\mu_{m})(y^{3}+b_{1}y^{2}+b_{2}y+b_{3}),\quad\text{where}$
(2.24)

$\begin{array}[]{ll}&b_{1}=t_{6}+\kappa(1-2\mathcal{R}_{z_{1}})\\\ \\\
&b_{2}=t_{7}+\kappa(\mathcal{R}_{z}^{*}-1)\\\ \\\
&b_{3}=t_{8}(\mathcal{R}_{z}^{*}-1),\end{array}$

and $t_{6}$, $t_{7}$ and $t_{8}$ defined as a positive linear combination of
$t_{i}$, $i=1,...,5$. The characteristic polynomial $r(y)$ gives five roots,
two of them are $y=-\mu$ and $y=-\mu_{m}$, whereas the Routh-Hurtwiz criteria
assures that the other three roots have negative real part if $b_{i}>0$ for
$i=1,2,3$ and $b_{1}b_{2}-b_{3}>0$. Clearly the coefficients $b_{i}>0$,
$i=1,2,3$ are all positive if $\mathcal{R}_{z}^{*}>1$ and
$1-2\mathcal{R}_{z_{1}}>0$. Thus, it follows that the endemic equilibrium
$\textbf{E}^{*}_{z}$ of System (2.9) is LAS if $\mathcal{R}_{z}^{*}>1$. To
conclude the qualitative analysis of Model (2.9), we prove the global
stability of the DFE ($\textbf{E}_{z_{0}}$) using similar techniques to those
used in Section 2.1.

######  Proposition 5.

If $\mathcal{R}^{*}_{z}>1$, the endemic equilibrium point $\textbf{E}^{*}_{z}$
defined in (2.23), is LAS in $\Omega_{z}$ defined in (2.10). Conversely, if
$\mathcal{R}_{z}^{*}<1$, the DFE $\textbf{E}_{z_{0}}$ of Model (2.2) defined
in (2.11), is GAS.

The proof of the global stability of $\textbf{E}_{z_{0}}$ when
$\mathcal{R}_{z}^{*}<1$ can be found in Appendix A.5.

##  3 Qualitative behaviour of the HIV/ZIKV model

In this section, we discuss the qualitative properties of the HIV/ZIKV co-
infection model (2.1). To achieve this purpose, we use the existence and
stability results as well as the definition of the basic reproduction number
for the HIV model $\mathcal{R}_{h}$ in (2.5) and $\mathcal{R}_{z}$ in (2.12)
obtained in Sections 2.1 and 2.2.

In Model (2.1), the total human population is denoted by
$N(t)=N_{h}(t)+N_{z}(t)$, and the total mosquito population is
$N_{m}(t)=S_{m}(t)+I_{m}(t)$. Additionally, to simplify algebraic calculations
we rename parameters:

$\begin{array}[]{ll}\kappa_{1}=\mu_{z}+\delta_{z}+\mu,&\kappa_{2}=\sigma_{1}+\mu,\\\
\\\ \kappa_{3}=\sigma_{2}+\mu_{hz}+\mu,&\kappa_{4}=\mu_{h}+\mu.\end{array}$
(3.1)

The interest region set is given by

$\Omega=\left\\{(S,I_{z},I_{h},I_{hz},A,R,S_{m},I_{m})\in\mathbb{R}_{+}^{8}:0\leq
N\leq\frac{\Lambda}{\mu};0\leq N_{m}\leq\frac{\Lambda_{m}}{\mu_{m}}\right\\}.$
(3.2)

As in the previous sections, it can be proved that $\Omega$ is positively
invariant under the flow of (2.1).

###  3.1 Computation of the basic reproduction number

The DFE for Model (2.1) is given by

$\textbf{E}_{0}=\left(\frac{\Lambda}{\mu},0,0,0,0,0,\frac{\Lambda_{m}}{\mu_{m}},0\right).$
(3.3)

Similarly to Sections 2.1 and 2.2, the basic reproduction number associated to
Model (2.1) can be determined through the matrices F and V given by

$\textbf{F}=\begin{bmatrix}\beta_{z}\frac{\Lambda}{\mu}&0&0&0&\beta_{m}\frac{\Lambda}{\mu}\\\
0&\beta_{m}\frac{\Lambda}{\mu}&0&0&0\\\ 0&0&0&0&0\\\
0&\sigma_{1}&\sigma_{2}&0&0\\\
\alpha_{m}\frac{\Lambda_{m}}{\mu_{m}}&0&\alpha_{m}\frac{\Lambda_{m}}{\mu_{m}}&0&0\end{bmatrix},\quad\text{and}\quad\textbf{V}=\begin{bmatrix}\kappa_{1}&0&0&0&0\\\
0&\kappa_{2}&0&0&0\\\ 0&0&\epsilon\delta_{z}+\kappa_{3}&0&0\\\
0&0&0&\kappa_{4}&0\\\ 0&0&0&0&\mu_{m}\end{bmatrix}.$

Here, the matrix $\textbf{FV}^{-1}$ is given by

$\textbf{FV}^{-}1=\begin{bmatrix}\beta_{z}\frac{\Lambda}{\mu\kappa_{1}}&0&0&0&\beta_{m}\frac{\Lambda}{\mu\mu_{m}}\\\
0&\beta_{m}\frac{\Lambda}{\mu\kappa_{2}}&0&0&0\\\ 0&0&0&0&0\\\
0&\frac{\sigma_{1}}{\kappa_{2}}&\frac{\sigma_{2}}{\epsilon\delta_{z}+\kappa_{3}}&0&0\\\
\alpha_{m}\frac{\Lambda_{m}}{\mu_{m}\kappa_{1}}&0&\alpha_{m}\frac{\Lambda_{m}}{\mu_{m}(\epsilon\delta_{z}+\kappa_{3})}&0&0\end{bmatrix}.$

The above matrix has as eigenvalues $\lambda_{1,2}=0$ (twice) and
$\lambda_{3}=\dfrac{\beta_{h}\Lambda}{\mu\kappa_{2}}=\mathcal{R}_{h}$, whereas
the other two eigenvalues are

$\lambda_{4,5}=\dfrac{\beta_{z}\Lambda}{2\mu\kappa_{1}}\pm\sqrt{\left(\dfrac{\beta_{z}\Lambda}{2\mu\kappa_{1}}\right)^{2}+\dfrac{\alpha_{m}\beta_{m}\Lambda_{m}\Lambda}{\mu\mu_{m}^{2}\kappa_{1}}}=\mathcal{R}_{z_{1}}\pm\mathcal{R}_{z_{2}},$

with the positive eigenvalues being

$\lambda_{4}=\dfrac{\beta_{z}\Lambda}{2\mu\kappa_{1}}+\sqrt{\left(\dfrac{\beta_{z}\Lambda}{2\mu\kappa_{1}}\right)^{2}+\dfrac{\alpha_{m}\beta_{m}\Lambda_{m}\Lambda}{\mu\mu_{m}^{2}\kappa_{1}}}=\mathcal{R}_{z_{1}}+\mathcal{R}_{z_{2}}=\mathcal{R}_{z}.$

Thus, the basic reproduction number for Model (2.1) is given by

$\mathcal{R}_{0}=\rho(\textbf{FV}^{-}1)=\max\\{\mathcal{R}_{h},\mathcal{R}_{z}\\},$
(3.4)

where $\mathcal{R}_{h}$ is the basic reproduction number for the HIV model
defined in the equation (2.5) and
$\mathcal{R}_{z}=\mathcal{R}_{z_{1}}+\mathcal{R}_{z_{2}}$ is the basic
reproduction number of the ZIKV model, which is defined in the equation
(2.12).

###  3.2 Local sensitivity analysis of the parameters

A local sensitivity analysis of $\mathcal{R}_{0}$ with respect to model
parameters allows quantifying parameter variations’ effect on the value of
$\mathcal{R}_{0}$. The sign of the sensitivity index denotes the direction of
the change, where a positive index for a particular parameter indicates that
increasing that parameter will increase $\mathcal{R}_{0}$ and vice versa. In
addition, the magnitude of the sensitivity index provides insight into the
relative significance of each parameter in the model predictions
banuelos2019using .

The normalized sensitivity index of a variable concerning a parameter is a
measure of how much the variable relatively changes to the change in the
parameter chitnis2008determining . In other words, the normalized sensitivity
index of the variable $X$, which is differentially dependent on parameter $p$,
is defined as:

$\Gamma_{p}^{X}=\dfrac{\partial X}{\partial p}\dfrac{p}{X}.$ (3.5)

Because $\mathcal{R}_{0}$ is defined as
$\max\\{\mathcal{R}_{h},\mathcal{R}_{z}\\}$, the sensitivity indices of
$\mathcal{R}_{0}$ with respect to the eleven parameters
$\\{\beta_{h},\Lambda,\mu,\sigma_{1},\beta_{z},\mu_{z},\delta_{z},\beta_{m},\alpha_{m},\Lambda_{m},\mu_{m}\\}$
in the expression of $\mathcal{R}_{0}$ in (3.4), can be determined for the
sensitivity indices of $\mathcal{R}_{h}$ and $\mathcal{R}_{z}$, respectively.
A calculation example of the sensitivity index of $\mathcal{R}_{z}$ with
respect to the parameter $\beta_{z}$ can be found in Appendix A.6.

###  3.3 Stability analysis

Similarly to in Sections 2.1 and 2.2, the local stability of the DFE
$\textbf{E}_{0}$ is determined by the linearization of (2.1) around an
arbitrary equilibrium. By ordering the equations and variables as
$S,S_{m},I_{z},I_{m},I_{h},I_{hz},A,R$ and computing the Jacobian matrix of
System (2.9) at an arbitrary equilibrium point
$\textbf{E}=(S,S_{m},I_{z},I_{m},I_{h},I_{hz},A,R)$, we obtain

$\textbf{J}(\textbf{E})=\begin{bmatrix}\textbf{J}_{11}(\textbf{E})&\textbf{0}&\textbf{0}\\\
\textbf{0}&-\kappa_{4}&\textbf{0}\\\ \textbf{0}&\textbf{0}&-\mu\end{bmatrix}.$
(3.6)

The matrix shown above has eigenvalues $-\mu<0$,
$-\kappa_{4}=-(\mu_{h}+\mu)<0$ and the eigenvalues of
$\textbf{J}_{11}(\textbf{E})$. By evaluating the matrix $\textbf{J}_{11}$ at
the DFE $\textbf{E}_{0}$ we obtain

$\textbf{J}_{11}(\textbf{E}_{0})=\begin{bmatrix}-\mu&\textbf{0}&\mathbf{\star}&\mathbf{\star}\\\
\textbf{0}&-\mu_{m}&\mathbf{\star}&\mathbf{\star}\\\
\textbf{0}&\textbf{0}&\textbf{M}_{33}&\mathbf{\star}\\\
\textbf{0}&\textbf{0}&\textbf{0}&\textbf{D}_{44}\end{bmatrix},\,\text{where}\,\textbf{M}_{33}=\begin{bmatrix}-\kappa_{1}+\dfrac{\Lambda\beta_{z}}{\mu}&\dfrac{\Lambda\beta_{m}}{\mu}\\\
\\\ \dfrac{\Lambda_{m}\alpha_{m}}{\mu_{m}}&-\mu_{m}\end{bmatrix}$

and $\textbf{D}_{44}$ is a diagonal matrix with entries
$-(\epsilon\delta_{z}+\kappa_{3})<0$ and
$\kappa_{2}-\frac{\beta_{h}\Lambda}{\mu}=\sigma_{1}+\mu-\frac{\beta_{h}\Lambda}{\mu}=(\sigma_{1}+\mu)(1-\mathcal{R}_{h})<0$
if and only if $\mathcal{R}_{h}<1$. Therefore, the remaining eigenvalues of
$\textbf{J}_{11}(\textbf{E}_{0})$ are those of the matrix $\textbf{M}_{33}$,
whose characteristic equation is

$\lambda^{2}+c_{1}\lambda+c_{2}=0,$ (3.7)

where

$\begin{array}[]{ll}&c_{1}=\kappa_{1}+\mu_{m}-\frac{\beta_{z}\Lambda}{\mu}=\mu_{m}+\kappa_{1}(1-2\mathcal{R}_{z_{1}})\\\
\\\
&c_{2}=\kappa_{1}\mu_{m}-\frac{\alpha_{m}\beta_{m}\Lambda_{m}\Lambda}{\mu\mu_{m}}-\frac{\beta_{z}\Lambda\mu_{m}}{\mu}=k_{1}\mu_{m}(1-\mathcal{R}_{z}^{*}).\end{array}$

Thus, the roots of the equation (3.7) have negative real part if and only if
$\mathcal{R}_{z}^{*}<1$ and $2\mathcal{R}_{z_{1}}<1$, where
$\mathcal{R}_{z}^{*}$ and $\mathcal{R}_{z_{1}}$ are defined in (2.14) and
(2.13), respectively. We have the following result:

######  Proposition 6.

If $\mathcal{R}_{0}=\max\\{\mathcal{R}_{h},\mathcal{R}_{z}\\}<1$, then
$\textbf{E}_{0}$ defined in (3.3) is LAS in $\Omega$ defined in (3.2).

Techniques similar to those used in Sections 2.1 and 2.2 can be applied to
confirm the presence of endemic solutions and assess their local and global
stability.

##  4 The control problem analysis

In this section, we formulate an optimal control problem (OCP) by adding three
intervention strategies to control the spread of HIV/ZIKV co-infection to
Model (2.1). The proposed approach mitigates HIV and zika infections through
the implementation of personal protection measures (such as the use of
repellents) using control $\eta_{1}$, the use of ART with control $\eta_{2}$
and preventive sexual contact (such as condom use) with control $\eta_{3}$.
The control functions $\eta_{1}$, $\eta_{2}$ and $\eta_{3}$ are defined in the
interval $[0,T$], where $T$ denotes the final time of the controls,
$0\leq\eta_{i}(t)\leq 1$ and $t\in[0,T]$ for $i=1,2,3$. Based on the
aforementioned considerations, the following OCP is formulated, where the
controls are shown in red for emphasis.

$\left\\{\begin{array}[]{ll}&\min\mathcal{J}(\eta)=\int_{0}^{T}\left(c_{1}I_{z}+c_{2}I_{h}+c_{3}I_{hz}+c_{4}I_{m}+c_{5}A+d_{1}\dfrac{\eta_{1}^{2}}{2}+d_{2}\dfrac{\eta_{2}^{2}}{2}+d_{3}\dfrac{\eta_{3}^{2}}{2}\right)dt\\\
\\\
&\dfrac{dS}{dt}=\Lambda-\left[{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}(1-\eta_{1})}\beta_{m}I_{m}+{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}(1-\eta_{3})}(\beta_{z}I_{z}+\beta_{h}I_{h})\right]S-\mu
S\\\ \\\
&\dfrac{dI_{z}}{dt}=\left[{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}(1-\eta_{1})}\beta_{m}I_{m}+{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}(1-\eta_{3})}\beta_{z}I_{z}\right]S-\omega_{2}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}(1-\eta_{3})}\beta_{h}I_{z}I_{h}-(\mu_{z}+\delta_{z}+\mu)I_{z}\\\
\\\
&\dfrac{dI_{h}}{dt}={\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}(1-\eta_{3})}\beta_{h}I_{h}S-\omega_{1}[{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}(1-\eta_{1})}\beta_{m}I_{m}+{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}(1-\eta_{3})}\beta_{z}I_{z}]I_{h}-[{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}(1-\eta_{2})}\sigma_{1}+\mu]I_{h}\\\
\\\
&\dfrac{dI_{hz}}{dt}=\omega_{1}[{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}(1-\eta_{1})}\beta_{m}I_{m}+{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}(1-\eta_{3})}\beta_{z}I_{z}]I_{h}+\omega_{2}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}(1-\eta_{3})}\beta_{h}I_{h}I_{z}-\epsilon\delta_{z}I_{hz}-[{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}(1-\eta_{2})}\sigma_{2}+\mu_{hz}+\mu]I_{hz}\\\
\\\
&\dfrac{dA}{dt}={\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}(1-\eta_{2})}\sigma_{1}I_{h}+{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}(1-\eta_{2})}\sigma_{2}I_{hz}-(\mu_{h}+\mu)A\\\
\\\ &\dfrac{dR}{dt}=\delta_{z}I_{z}+\epsilon\delta_{z}I_{hz}-\mu R\\\ \\\
&\dfrac{dS_{m}}{dt}=\Lambda_{m}-{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}(1-\eta_{1})}\alpha_{m}(I_{z}+I_{hz})S_{m}-\mu_{m}S_{m}\\\
\\\
&\dfrac{dI_{m}}{dt}={\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}(1-\eta_{1})}\alpha_{m}(I_{z}+I_{hz})S_{m}-\mu_{m}I_{m}\\\
\\\ &{\bf X}(0)=(S_{0},I_{z0},I_{h0},I_{hz0},A_{0},R_{0},S_{m0},I_{m0})={\bf
X}_{0}\\\ \\\ &{\bf
X}(T)=(S_{f},I_{zf},I_{hf},I_{hzf},A_{f},R_{f},S_{mf},I_{mf})={\bf
X}_{f}.\end{array}\right.$ (4.1)

In the above formulation $\eta=(\eta_{1}(t),\eta_{2}(t),\eta_{3}(t))$, and
$c_{1},c_{2},c_{3},c_{4},c_{5},d_{1},d_{2},$ and $d_{3}$ are positive weights.
Therefore, we seek an optimal control $\eta^{*}(t)$ determined as

$\mathcal{J}(\eta^{*}(t))=\min\left\\{\mathcal{J}(\eta(t)|\eta\in\mathcal{A})\right\\},$
(4.2)

with a set $\mathcal{A}$ of controls defined as

$\mathcal{A}=\left\\{\eta(t)=(\eta_{1}(t),\eta_{2}(t),\eta_{3}(t))|0\leq\eta_{1}(t)\leq\eta_{1}^{\max},0\leq\eta_{2}(t)\leq\eta_{2}^{\max},0\leq\eta_{3}(t)\leq\eta_{3}^{\max}\right\\},$

where $\eta_{i}^{\max}\leq 1$, $i=\\{1,2,3\\}$ and $\eta$ is Lebesgue
measurable. In order to define the formulation of our OCP using Pontryagin’s
Maximum Principle (PMP) pontryagin2018ls , we have the Lagrangian as

$L=c_{1}I_{z}(t)+c_{2}I_{h}(t)+c_{3}I_{hz}(t)+c_{4}I_{m}(t)+c_{5}A(t)+d_{1}\dfrac{\eta_{1}^{2}(t)}{2}+d_{2}\dfrac{\eta_{2}^{2}(t)}{2}+d_{3}\dfrac{\eta_{3}^{2}(t)}{2},$
(4.3)

and we determine the Hamiltonian function as

$H=L(I_{z},I_{h},I_{hz},I_{m},A,\eta)+p_{1}\dfrac{dS}{dt}+p_{2}\dfrac{dI_{z}}{dt}+p_{3}\dfrac{dI_{h}}{dt}+p_{4}\dfrac{dI_{hz}}{dt}+p_{5}\dfrac{dA}{dt}+p_{6}\dfrac{dR}{dt}+p_{7}\dfrac{dS_{m}}{dt}+p_{8}\dfrac{dI_{m}}{dt}.$
(4.4)

In the remainder, we investigate the minimum value of Lagrangian (4.3).
Firstly, we must prove existence of the optimal control $\eta^{*}$ according
to the controlled system (4.1).

######  Proposition 7.

There exists an optimal control $\eta^{*}$ such that

$\mathcal{J}\left(\eta^{*}(t)\right)=\min_{\eta\in\mathcal{A}}(\mathcal{J}(\eta(t))),$

subject to Control System (4.1) with initial conditions as $\bf{X}_{0}$.

The proof of the above proposition can be found in Appendix A.7.

In the following, we apply PMP pontryagin2018ls  to provide a characterization
of an optimal control solution to the Hamiltonian (4.4) subject to the OCP
(4.1). If $(X^{*},\eta^{*})$ is an optimal solution for the controlled system
(4.1), then there exists a non trivial vector function
$p=(p_{1},p_{2},p_{3},p_{4},p_{5},p_{6},p_{7},p_{8})$, such that

$\dfrac{\partial H}{\partial\eta_{i}}=0,\,i=1,2,3\quad\text{ and
}\dot{p}_{i}=\dfrac{dp_{i}}{dt}=-\dfrac{\partial H}{\partial
X_{i}},\,i=1,\ldots,8.$ (4.5)

######  Proposition 8.

Let $(S^{*},I_{z}^{*},I_{h}^{*},I_{hz}^{*},A^{*},R^{*},S_{m}^{*},I_{m}^{*})$
be the optimal state variables solution associated to the optimal control
variable $\eta^{*}$ subject to the control problem (4.2). Then, there exists
an adjoint vector $p$ that satisfies the controlled system (4.1), with
transversality conditions $p_{i}(T)=0$, for $i=1,\ldots,8$, where the optimal
controls are

$\begin{array}[]{ll}\eta_{1}^{*}=&\dfrac{\left((p_{2}-p_{1})S+(p_{4}-p_{3})I_{h}\right)\beta_{m}I_{m}+(p_{8}-p_{7})\alpha_{m}(I_{z}+I_{hz})S_{m}}{d_{1}}\\\
\\\
\eta_{2}^{*}=&\dfrac{\left(p_{5}-p_{3}\right)\sigma_{1}I_{h}+\left(p_{5}-p_{4}\right)\sigma_{2}I_{hz}}{d_{2}}\\\
\\\
\eta_{3}^{*}=&\dfrac{(p_{2}-p_{1})\beta_{z}I_{z}S+(p_{3}-p_{1})\beta_{h}I_{h}S+(p_{4}-p_{3})\omega_{1}\beta_{z}I_{z}I_{h}+(p_{4}-p_{2})\omega_{2}\beta_{h}I_{z}I_{h}}{d_{3}}.\end{array}$
(4.6)

The proof is in Appendix A.8.

##  5 HIV/ZIKV co-infection in Brazil and Colombia

As mentioned previously, during 2015 and 2016, Brazil and Colombia experienced
an alarming occurrence involving the co-infection of two significant viral
diseases: HIV/AIDS and zika. This phenomenon emerged as a complex public
health issue challenging the healthcare system in these regions. In fact, the
co-infection highlighted the need for further research, comprehensive
surveillance, and implementation of prevention and early reaction strategies
to mitigate the impact of these dual infections.

Owing to the lack of temporal records on HIV/ZIKV co-infection to date, it was
no possible to estimate the parameters of Model (2.1) in this section.
However, specific parameter values were derived from available demographic
information and previous research conducted on zika and HIV/AIDS in Colombia
and Brazil. In cases where there was a lack of available data, these values
were either estimated based on specific assumptions or adapted from research
conducted on different regions or diseases. The following outlines the main
assumptions for extracting these parameter values.

Firstly, the precise relationship between humans and population densities of
Aedes aegypti remains uncertain and has been a subject of considerable
interest and investigation within the field of entomology. Examining the
existing literature, it can be inferred that during 2015-2016, there existed
an approximate ratio of one human to three female Aedes aegypti mosquitoes
(1:3) in Colombia and Brazil. Secondly, the modification parameter related to
the recovery rate of humans $\epsilon$ was assumed to be much smaller than one
because, according to research in this area, an individual with HIV/ZIKV co-
infection recovers more slowly than an individual with only zika
joob2017concurrent ; sherman2018zika ; rothan2018 ; bidokhti2018siv ;
aschengrau2021international . Thirdly, the parameters associated with the co-
infection probability $w_{1}$ and $w_{2}$ were assumed to be small because by
2016, only five cases had been reported in Colombia and one in Brazil
calvet2016 ; villamil2018 . However, a hypothetical scenario in which these
parameters are escalated is depicted to demonstrate the correlation between
rising probabilities of co-infection and the subsequent increase in the number
of co-infected individuals. Finally, certain parameter values pertaining to
HIV/AIDS were predominantly adjusted at a population level, drawing upon
estimations derived from Luxembourg, the Czech Republic, Japan, Croatia, the
United Kingdom, and Mexico prieto2021current .

As in gao2016 , we assume that Colombia and Brazil share some common parameter
values (see Table 5.1), except for those related to country population size
and initial conditions (see Table 5.2).

Due to the dispersion of the units of measurement of the parameters in the
different sources consulted, all parameter values were adapted to day as the
unit of measurement.

Parameter | Range | Reference
---|---|---
$\delta_{z}$ | [1.11e-2,7.14e-2] | gao2016
$\sigma_{1}$ | [2.1e-4,5.46e-4] | prieto2021current
$\mu_{h}$ | 5.46e-5 | prieto2021current
$\mu_{m}$ | [5e-2,1.2e-1] | gao2016
$\mu_{z}$ | [3.42e-5,5.46e-5] | Assumed
$\mu_{hz}$ | 0 | calvet2016 ; villamil2018
$\sigma_{2}$ | 4.55e-5 | calvet2016 ; villamil2018
$\omega_{1}$ | [1e-3,1e-2] | calvet2016 ; villamil2018
$\omega_{2}$ | [1e-3,1e-2] | calvet2016 ; villamil2018
$\epsilon$ | 1e-10 | calvet2016 ; villamil2018
Table 5.1: Shared (Colombia-Brazil) range of parameter values (minimum and maximum) involved in Model (2.1). Time in days. Parameter | Colombia | Brazil | Reference
---|---|---|---
$\beta_{m}$ | [6.25e-10,1.25e-9] | [2.42e-10,2.9e-10] | gao2016
$\beta_{z}$ | [2.08e-12,1.04e-11] | [4.83e-13,2.42e-12] | gao2016
$\beta_{h}$ | [6.25e-12,1.46e-11] | [1.45e-12,3.38e-12] | prieto2021current
$\alpha_{m}$ | [2.08e-10,4.17e-10] | [4.83e-11,9.66e-11] | gao2016
$\Lambda$ | [1.73e+03 1.73e+03] | [7.44e+03,7.44e+03] | biswas2020mathematical
$\mu$ | 3.21e-5 | 3.21e-5 | biswas2020mathematical
$\Lambda_{m}$ | [1.03e+7,1.03e+7] | [4.43e+7,4.43e+7] | biswas2020mathematical
Table 5.2: Range of parameter values (minimum and maximum) involved in Model (2.1) that differ in Colombia and Brazil. Time in days.  Parameter | Colombia | Brazil | Reference
---|---|---|---
$S(0)$ | 28,800,000 | 124,200,000 | gao2016
$I_{z}(0)$ | 9,600,000 | 41,400,000 | gao2016
$I_{h}(0)$ | 4,320,000 | 18,630,000 | prieto2021current
$I_{hz}(0)$ | 4.8 | 20.7 | calvet2016 ; villamil2018
$A(0)$ | 2,832,000 | 12,213,000 | calvet2016 ; villamil2018
$R(0)$ | 2,448,000 | 10,557,000 | calvet2016 ; villamil2018
$S_{m}(0)$ | 33,600,000 | 144,900,000 | gao2016
$I_{m}(0)$ | 14,400,000 | 62,100,000 | gao2016
Table 5.3: Initial conditions involved in Model (2.1).

Once parameter values have been set, the case study is organized into three
separate stages. The initial stage involves numerical determination of the
basic reproduction number defined in Section 3.1. Furthermore, numerical
values of the sensitivity indices determined analytically in Section 3.2 are
also presented. In the second stage, the uncontrolled model defined by (2.1)
is simulated numerically. Finally, the third stage focuses on numerically
simulating the control problem described in (4.1). This sequential framework
allows for a thorough exploration of the various aspects under investigation,
providing a comprehensive analysis of the case study.

###  5.1 Value of the basic reproduction number and its sensitivity indices

In Section 3.1, we calculated the basic reproduction number (3.4) for Model
(2.1). Figure 5.1 shows the possible values for the basic reproduction number
using the extreme parameter values given in Tables 5.1 and 5.2 for Colombia
and Brazil, that result in the minimum or maximum value of $\mathcal{R}_{0}$,
respectively.

(a) Colombia (b) Brazil

Figure 5.1: Possible values for the thresholds using the extreme parameter
values given in Tables 5.1 and 5.2 for Colombia and Brazil, that result in the
minimum or maximum value of $\mathcal{R}_{0}$, respectively. In each case,
$\mathcal{R}_{0}=\max\\{\mathcal{R}_{h},\mathcal{R}_{z}\\}$. The vertical red
line represents $\mathcal{R}_{h,z}=1$.

The normalized sensitivity indices of $\mathcal{R}_{h}$ and $\mathcal{R}_{z}$
summarized in Figure 5.2, are obtained using the values of the parameters in
Tables 5.1 and 5.2.

Figure 5.2: Normalized sensitivity index of $\mathcal{R}_{h}$ to the
parameters $\\{\beta_{h},\Lambda,\mu,\sigma_{1}\\}$ and $\mathcal{R}_{z}$ to
the parameters
$\\{\Lambda,\mu,\beta_{z},\mu_{z},\delta_{z},\beta_{m},\alpha_{m},\Lambda_{m},\mu_{m}\\}$.
Thus, the red dot line represents sensitivity of zero. A calculation example
using Equation (3.5) can be found in Appendix A.6.

.

Regarding $\mathcal{R}_{h}$, the primary driver is the human death rate $\mu$,
followed by the infection rate $\beta_{h}$ of humans through sexual contact
with HIV-infected individuals and the human recruitment rate $\Lambda$. This
implies that a reduction in the the human death rate $\mu$ would result in the
most significant increase in the likelihood of HIV infection among humans.
Regarding $\mathcal{R}_{z}$, the outcomes reveal that the parameter with the
highest impact is the mosquitoes death rate $\mu_{m}$: longer surviving
mosquitoes lead to an increase in $\mathcal{R}_{z}$. Additionally, the human
recruitment rate $\Lambda$ and death rate $\mu$ are also essential factors
influencing $\mathcal{R}_{z}$.

###  5.2 Evaluation of uncontrolled population behaviour over time

We then performed numerical simulations of the uncontrolled Model (2.1). We
begin by representing the DFE for Colombia and Brazil. The DFE is obtained
when the value of $\mathcal{R}_{0}$ is less than one. As shown in Figure 5.1,
for Colombia and Brazil, $\mathcal{R}_{0}$=0.47192. Figure 5.3 shows the
behaviour of the solutions for both human and mosquito populations. In both
cases, on day 50 after the first observation, infected individuals (with
either disease) tended to zero, whereas susceptible humans are stabilized. A
similar scenario was observed in the mosquito population. On day 25 of the
first observation, infected mosquitoes decreased to zero, and susceptible
mosquitoes stabilized. In this scenario, HIV/AIDS cases outnumbered ZIKV
cases. No outbreaks of infected individuals were observed under the initial
conditions used. Similarly, the endemic equilibrium is obtained when the value
of $\mathcal{R}_{0}$ is greater than one (see Figure 5.1). This scenario is
simulated and depicted in Figure 5.4. For Colombia $\mathcal{R}_{0}$=3.2338
and for Brazil $\mathcal{R}_{0}$=3.2326. We can observe a significant increase
in the number of individuals infected with ZIKV during the first 70 days of
the first observation, with the number of HIV/AIDS-infected individuals being
lower. Individuals infected with ZIKV reached maximum values on day 70 of the
first observation. In contrast, individuals infected with HIV/AIDS did not
generate peaks and appear to have a constant behaviour throughout the observed
period (approximately a year). During the first 250 days from the first
observation, ZIKV-infected individuals outnumber HIV/AIDS-infected
individuals. However, after day 250, HIV/AIDS-infected individuals constantly
outnumber ZIKV-infected individuals.

For the mosquito population, there is always a higher population density of
susceptible than infected mosquitoes, but there is stable coexistence between
them. Both countries differ mainly in population density (Brazil has a larger
population than Colombia).The simulations indicated that the behaviour of both
populations in Colombia and Brazil was highly similar, with a notable
distinction in terms of population densities. This suggests that the dynamics
and patterns of infection transmission were comparable between the two
countries, with variations primarily related to the size of the populations
being studied. An interesting observation from the simulations was that the
behaviour of individuals infected with ZIKV displayed more significant
variability than those infected with HIV/AIDS. This implies that the
progression and manifestations of the ZIKV infection were more diverse and
fluctuating, potentially influenced by various factors such as environmental
conditions, viral load, or other epidemiological characteristics. In contrast,
the behaviour of HIV/AIDS-infected individuals appeared to be more consistent
and stable over time. These findings highlight the distinctive characteristics
of ZIKV and HIV/AIDS infections, particularly in terms of the dynamic nature
of zika infection compared to the relatively more constant nature of HIV/AIDS
infection.

(a) Colombia (b) Colombia (c) Brazil (d) Brazil

Figure 5.3: Simulations of the DFE for Colombia and Brazil. For Colombia,
$\mathcal{R}_{h}$= 0.4719, $\mathcal{R}_{z}$= 0.2306. For Brazil,
$\mathcal{R}_{h}$= 0.4719, $\mathcal{R}_{z}$= 0.2977. In both countries,
$\mathcal{R}_{0}$= 0.4719. The solutions tend to the DFE.

A hypothetical scenario for HIV/ZIKV co-infected individuals is shown in
Figure 5.5. We contrasted three different possibilities for the probability of
co-infection: $\omega_{1,2}$=0.001, 0.0055, 0.01. The first value (0.001) was
named low probability, the second value (0.0055) medium probability and the
last value (0.01) high probability. It can be seen that for higher values of
this pair of parameters, a higher population density of individuals co-
infected with ZIKV and HIV follows. Evidently, an increase in these two
probabilities increased the density of co-infected individuals.

(a) Colombia (b) Colombia (c) Brazil (d) Brazil

Figure 5.4: Simulations of the endemic equilibrium. For Colombia,
$\mathcal{R}_{h}$= 3.2308, $\mathcal{R}_{z}$= 3.2338. For Brazil,
$\mathcal{R}_{h}$= 3.2308, $\mathcal{R}_{z}$= 3.2326. The solutions tend to
the endemic equilibrium.

(a) Colombia (b) Brazil

Figure 5.5: HIV/ZIKV co-infected human population for different values of the
co-infection probabilities $\omega_{1}$ and $\omega_{2}$.

These numerical experiments showed that for values of $\mathcal{R}_{0}<1$, the
HIV/ZIKV co-infection epidemics tend to fade away after some time, with
infected individuals tending to zero and the susceptible population reaching a
stable level. The same trend was observed for the mosquito population, with
infected mosquitoes decreasing to zero and susceptible mosquitoes stabilizing.
A noticeable increase in the number of individuals infected with ZIKV, with a
lower number of individuals infected with HIV/AIDS, was observed during the
first 70 days of the first observation for the values of $\mathcal{R}_{0}>1$.
Infected individuals with ZIKV peaked on day 70, whereas HIV-infected
individuals did not generate any peaks during the 300 simulated days.
Regarding the mosquito population, susceptible mosquitoes always outnumber
infected mosquitoes; however, both coexist.

In addition, Figure 5.5 provides insights into the hypothetical scenario of
ZIKV/HIV co-infection. Comparing three different scenarios for co-infection
probabilities (high, medium and low), an increase in these probabilities above
0.01, leads to an increase in the density of co-infected individuals. These
findings suggest that the values of these parameters can significantly impact
the spread and severity of HIV/ZIKV co-infection, highlighting the importance
of considering multiple factors when designing effective strategies to prevent
and control the spread of these viruses, particularly in pregnant women to
avoid congenital malformations in children.

These findings suggest that the behaviour of populations in Colombia and
Brazil was similar but different in population densities. This suggests that
the patterns of infection transmission were comparable, primarily influenced
by population size. Notably, individuals infected with ZIKV exhibited more
significant behaviour variability than those infected with HIV/AIDS. This
indicates that ZIKV infection showed diverse and fluctuating progression and
manifestations, potentially influenced by various factors. In contrast,
HIV/AIDS-infected individuals displayed more consistent and stable behaviour
over time.

###  5.3 Evaluation of controlled population behaviour over time

This stage focuses on numerically simulating the control problem described by
(4.1). To this end, we investigate the effects of mixing strategies to control
HIV/ZIKV co-infection in Colombia and Brazil. Table 5.4 shows the balancing
and weighting constants values associated with the OPC (4.1).

Parameter | $c_{1}$ | $c_{1}$ | $c_{3}$ | $c_{4}$ | $d_{1}$ | $d_{2}$ | $d_{3}$
---|---|---|---|---|---|---|---
Value | 0.5 | 0.5 | 0.5 | 0.5 | $10^{5}$ | $10^{7}$ | $7\times 10^{6}$
Table 5.4: Balancing and weighting constants values associated to the OPC
(4.1).

On the one hand, controlling the spread of ZIKV and HIV in Colombia requires
specialized actions to address the particular problems that the population has
to overcome. To combat ZIKV, efforts must be directed toward preventing
mosquito bites using insect repellents, protective clothing, and mosquito
control measures. Therefore, the zika Fever Response Plan world2016zika
enacted in 2016 was a priority for the Colombian government and focused on
improving public health surveillance systems for vector-transmitted diseases
and raising awareness of health issues forero2020zika . To prevent sexually
transmitted infections from spreading, safe sex practices must be emphasized
in public health campaigns among vulnerable populations prachniak2016hiv .
Therefore, ART and condoms are critical for controlling the spread of the
human immunodeficiency virus. To ensure that people living with HIV in
Colombia receive the care they require to be healthy and to lower their risk
of spreading the virus to others, it is essential to increase access to HIV
testing and treatment. Condom use and other harm reduction techniques must
also be promoted in public health campaigns, especially among high-risk
groups, such as drug injectors. By addressing these challenges through a
comprehensive approach, including education, prevention, and access to care,
Colombia can make significant progress in controlling the spread of ZIKV and
HIV.

On the other hand, Brazil has implemented a countrywide strategy to manage
Aedes mosquito populations and to stop the spread of ZIKVbancroft2022vector .
This effort has incorporated activities, including insecticide application,
the elimination of standing water, and public awareness campaigns to motivate
people to take action to lessen mosquito breeding grounds, which are
incorporated as $\eta_{1}$ in the epidemic model (4.1). Additionally, it has
been noted that zika can be transmitted sexually; hence, condom use is advised
as a preventative measure. Brazil has implemented a comprehensive strategy for
HIV prevention and care benzaken2019antiretroviral . This has involved
administering ART to all HIV-positive individuals since 2013. Brazil started a
campaign to distribute condoms, encouraging condom use among vulnerable
populations. As a result, the number of AIDS-related deaths in Brazil has
significantly decreased.

Hence, using the controls $\eta_{1}^{*},\eta_{2}^{*},$ and $\eta_{3}^{*}$
determined in Proposition 8 to optimize the controlled dynamical system (4.1),
and by simulating the endemic scenario for Colombia and Brazil (Figure 5.4),
we obtain in Figures 5.6–5.8 a timely controlled pattern. The sizes of humans
infected with ZIKV, HIV/ZIKV, and HIV/AIDS are reduced due to mixed strategies
to maintain the level of infection. However, the size of the HIV-infected
population is increasing because of ART, resulting in better healthcare for
the infected population, which helped to reduce the number of individuals with
AIDS. In contrast, in Figure 5.6, ZIKV prevention measures helped reduce the
number of infected mosquitoes and their contact with susceptible individuals.

In Figure 5.7, the control $\eta_{1}$ is required at its full maximum at the
beginning by establishing measures such as mosquito repellents and nets,
elimination of mosquito breeding sites, and protective clothing. $\eta_{2}$
goes in three months to its maximum capacity to contain the number of
individuals who attended the AIDS stage, which slows the progression of the
disease, prevents further damage to the immune system, and reduces the risk of
HIV transmission. Additionally, $\eta_{3}$ describes a gradual campaign to
promote condom use can help reduce the transmission of HIV and increase
awareness about the importance of protecting oneself and one’s partners from
sexually transmitted infections, such as educating people about the benefits
of using condoms to prevent the spread of HIV, increasing easy access to
condoms, and recruiting community leaders, influential individuals, or social
media influencers to encourage the use of condoms.

(a) Colombia (b) Brazil

Figure 5.6: Simulations of the mosquito population in Colombia and Brazil with
the activation of different controls.

(a) Colombia (b) Brazil

Figure 5.7: Simulations of the behaviour of the controls in Colombia and
Brazil.

Compartments $I_{z}$ and $A$ are controlled to a smaller size than the
$\eta_{1}=\eta_{2}=\eta_{3}=0$ simulation, while compartments $I_{h}$ and
$I_{hz}$ have a larger size owing to the ART that prevents infected people
from reaching the AIDS stage. However, the mixed strategies were maintained at
their approximate total capacities, as illustrated in Figure 5.8.

In addition, the size of the infected mosquito population remained below the
initial condition which prevented the spread of ZIKV due to the effect of the
$\eta_{1}$ control.

(a) Colombia (b) Brazil

Figure 5.8: Simulations of controlling infected human sizes in Colombia and
Brazil.

Finally, to corroborate the behaviour of the co-infected human and the
mosquitoes populations with respect to the variation of the combined controls
$\eta_{1}$, $\eta_{2}$ and $\eta_{3}$, (see Figures 5.9 and 5.10), we have
performed some additional numerical experiments. Figure 5.9 shows that the
concurrent activation of $\eta_{1}$ (personal protection) and $\eta_{2}$ (ART)
generate comparable trends in the co-infected individuals with zika and HIV
within the contexts of Colombia and Brazil. Hence, in Figures 5.9 (a)-(b), we
investigated the co-infection patterns by fixing the control $\eta_{1}$ for
personal protection and increasing gradually the ART control from 0 to its
maximum value 1 with an increment value of
$\Delta\eta_{1}=\Delta\eta_{2}=0.003$. In both cases for Brazil and Columbia,
the increase of the personal protection did not play a significant role as the
gradual increase of ART for the co-infection. However, the personal protection
plays an important role to reduce the overall ZIKV infected population.
Analogously, the simultaneous implementation of ART and condoms use
($\eta_{3}$) in Colombian and Brazil reduce the number of co-infected
individuals burden. Hence, in Figures 5.9 (c-d), the condom use plays a more
significant role in reducing the size of co-infected size.

(a) Colombia ($\eta_{1}$-$\eta_{2}$)  (b) Brazil ($\eta_{1}$-$\eta_{2}$)  (c)
Colombia ($\eta_{2}$-$\eta_{3}$)  (d) Brazil ($\eta_{2}$-$\eta_{3}$)

Figure 5.9: Simulation of the effects of combined $\eta_{1}$-$\eta_{2}$ and
$\eta_{2}$-$\eta_{3}$ on the co-infected human population in Colombia and
Brazil.

This supports the outcomes derived from the simulations presented in Figure
5.8, wherein the control strategy involving the simultaneous application of
all controls above was previously established. The dynamics exhibited by the
population of zika virus-carrying mosquitoes, when subjected to concurrent
interventions involving personal protection, ART, and condom use, diverge from
those observed within the co-infected human population. This contrast is
evident upon examining Figure 5.10, wherein distinct patterns manifest for
Colombia and Brazil. Notably, the impact of personal protection emerges as
more pronounced than that of ART, as underscored by the observation that
mosquito density exhibits heightened levels when the efficacy of personal
protection approaches to zero. The limited influence of sexual transmission
dynamics of zika and HIV on mosquito density substantiates this outcome.
Conversely, upon simultaneous implementation of ART and condom use in Colombia
and Brazil, the peaks in mosquito densities materialize when both controls are
administered at 50%. In contrast, the lowest densities align with situations
with no controls or controls implemented at 100%. These observations
significantly explain the complex interplay between control strategies and
their varying effects on human and mosquito populations in the HIV/ZIKV co-
infection context. This information could guide policymakers and public health
authorities in making informed decisions regarding the most effective
strategies for managing and controlling these diseases.

(a) Colombia ($\eta_{1}$-$\eta_{2}$) (b) Brazil($\eta_{1}$-$\eta_{2}$) (c)
Colombia ($\eta_{2}$-$\eta_{3}$) (d) Brazil ($\eta_{2}$-$\eta_{3}$)

Figure 5.10: Simulation of the effects of combined $\eta_{1}$-$\eta_{2}$ and
$\eta_{2}$-$\eta_{3}$ on the infected mosquito population in Colombia and
Brazil.

##  6 Discussion and concluding remarks

The first reported cases of HIV/ZIKV co-infection in Colombia and Brazil
highlighted the potential interactions between these two viruses. The co-
circulation of these infections in South America presented a unique challenge
for public health officials, as both viruses were known to be transmitted
through sexual contact and had overlapping clinical symptoms.

The mathematical modelling conducted in this study was used to consider the
transmission dynamics of both viruses in co-infected individuals, allowing for
the evaluation of different intervention scenarios. The findings of this study
were particularly relevant during the 2015-2016 zika outbreak in South
America, which overlapped with the ongoing HIV epidemic in the region. By
identifying the patterns of co-infection and their potential impact on
transmission, the results of this study could inform public health strategies
to control and prevent the spread of both viruses.

On the one hand, numerical simulations conducted with parameter values for
Colombia provided valuable insights into the potential impact of HIV/ZIKV co-
infection. By comparing three distinct scenarios about co-infection
probabilities (high, medium, and low), the study demonstrated that an increase
in these parameters led to a higher density of co-infected individuals. By
comparing the qualitative behaviour of zika and HIV in Brazil and Colombia, we
found the following results: First, both countries exhibited similar trends in
spreading and stabilizing these viruses, with differences primarily
attributable to population density. Second, individuals infected with ZIKV
demonstrated greater variability in their behaviour than those infected with
HIV/AIDS. This suggests that ZIKV infection displayed a diverse and
fluctuating progression with various manifestations, potentially influenced by
multiple factors. In contrast, individuals infected with HIV/AIDS exhibited
more consistent and stable behaviour over time. Third, the simulation results
demonstrated that increasing the values of the co-infection probabilities
$\omega_{1}$ and $\omega_{2}$ could lead to a higher population density of co-
infected individuals. This finding underscores the importance of implementing
effective strategies to prevent and control the spread of these viruses,
particularly in pregnant women, to prevent congenital complications.

The primary preventative methods for zika and HIV infection included avoiding
mosquito bites, using condoms during sexual activity, and a treatment strategy
such as antiretroviral therapy (ART). Given that zika is mainly spread by
mosquitoes, we established a timely control to describe measures such as using
insect repellent, donning long sleeves, and staying indoors during prime
mosquito-biting hours are all personal protection measures. In addition, we
illustrated a control for condoms use to lower the risk of transmission of
zika, which can also spread sexually. Similarly, since sexual contact is the
primary way HIV spreads, the most effective strategy to prevent HIV
transmission during sex is to use condoms appropriately and regularly.
Additionally, we incorporated a control variable for ART, which included using
a variety of HIV medications to successfully suppress the virus, maintain or
enhance immune function, and lower the risk of transmission.

In summary, this study highlighted the need for continued research on the
transmission dynamics of zika and HIV/AIDS and developing effective
intervention strategies to control and prevent their spread. Future work in
this field plans to incorporate compartments of women giving birth to babies
with and without congenital malformations to understand the impact of co-
infection in children better. Including these compartments would allow for a
more detailed assessment of the long-term effects of co-infection on child
health outcomes, including the potential for developmental delays,
neurological deficits, and other complications. This approach could also
facilitate the development of more targeted prevention and treatment
strategies for the affected children and their families. Ultimately, this
research is critical for improving our understanding of the complex
interactions between HIV and zika and developing effective public health
interventions to mitigate their impact on affected individuals and
communities.

## Acknowledgements

The authors appreciate the support provided by the One Health Modelling
Network for Emerging Infections (OMNI-RÉUNIS) and Mathematics for Public
Health (MfPH), which are financially supported by the Natural Sciences and
Engineering Research Council of Canada (NSERC) and the Public Health Agency of
Canada (PHAC).

## References

  * (1) W. E. Villamil-Gómez, Álvaro Rosendo Sánchez-Herrera, H. Hernández-Prado, J. Hernández-Iriarte, K. Díaz-Ricardo, O. Vergara-Serpa, J. Castellanos, J. E. Portela-Gaviria, S. Patiño-Valencia, D. C. Vargas-Bedoya, A. F. Orozco-Villa, S. S. Reyes-Guerrero, J. D. Castrillón-Spitia, D. R. Murillo-García, V. Henao-SanMartin, J. J. Londoño, H. Bedoya-Rendón, J. de Jesús Cárdenas-Pérez, J. A. Cardona-Ospina, G. J. Lagos-Grisales, and A. J. Rodríguez-Morales, “Zika virus and HIV co-infection in five patients from two areas of Colombia,” 2018.
  * (2) V. M. Cao-Lormeau, A. Blake, S. Mons, S. Lastère, C. Roche, J. Vanhomwegen, T. Dub, L. Baudouin, A. Teissier, P. Larre, A. L. Vial, C. Decam, V. Choumet, S. K. Halstead, H. J. Willison, L. Musset, J. C. Manuguerra, P. Despres, E. Fournier, H. P. Mallet, D. Musso, A. Fontanet, J. Neil, and F. Ghawché, “Guillain-Barré syndrome outbreak associated with Zika virus infection in French Polynesia: a case-control study,” The Lancet, vol. 387, 2016.
  * (3) E. Oehler, L. Watrin, P. Larre, I. Leparc-Goffart, S. Lastãre, F. Valour, L. Baudouin, H. P. Mallet, D. Musso, and F. Ghawche, “Zika virus infection complicated by Guillain-Barré syndrome–case report, French Polynesia, December 2013,” Eurosurveillance, vol. 19, 2014.
  * (4) D. W. Smith and J. Mackenzie, “Zika virus and Guillain-Barré syndrome: another viral cause to add to the list,” 2016.
  * (5) J. Mlakar, M. Korva, N. Tul, M. Popović, M. Poljšak-Prijatelj, J. Mraz, M. Kolenc, K. R. Rus, T. V. Vipotnik, V. F. Vodušek, A. Vizjak, J. Pižem, M. Petrovec, and T. A. Županc, “Zika virus associated with microcephaly,” New England Journal of Medicine, vol. 374, 2016.
  * (6) G. Calvet, R. S. Aguiar, A. S. Melo, S. A. Sampaio, I. De Filippis, A. Fabri, E. S. Araujo, P. C. de Sequeira, M. C. de Mendonça, L. de Oliveira, et al., “Detection and sequencing of Zika virus from amniotic fluid of fetuses with microcephaly in Brazil: a case study,” The Lancet Infectious Diseases, vol. 16, no. 6, pp. 653–660, 2016.
  * (7) A. F. Kabapy, H. Z. Shatat, and E. W. A. El-Wahab, “Attributes of HIV infection over decades (1982–2018): A systematic review and meta-analysis,” 2020.
  * (8) J. C. Tilton and R. W. Doms, “Entry inhibitors in the treatment of HIV-1 infection,” 2010.
  * (9) B. Atkinson, P. Hearn, B. Afrough, S. Lumley, D. Carter, E. J. Aarons, A. J. Simpson, T. J. Brooks, and R. Hewson, “Detection of Zika virus in semen,” 2016\.
  * (10) A. C. Gourinat, O. O’Connor, E. Calvez, C. Goarant, and M. Dupont-Rouzeyrol, “Detection of Zika virus in urine,” Emerging Infectious Diseases, vol. 21, 2015.
  * (11) D. Musso, C. Roche, E. Robin, T. Nhan, A. Teissier, and V.-M. Cao-Lormeau, “Potential sexual transmission of Zika virus,” Emerging infectious diseases, vol. 21, no. 2, p. 359, 2015.
  * (12) H. A. Rothan, M. R. Bidokhti, and S. N. Byrareddy, “Current concerns and perspectives on zika virus co-infection with arboviruses and HIV,” Journal of autoimmunity, vol. 89, pp. 11–20, 2018.
  * (13) K. E. Sherman, S. D. Rouster, L. X. Kong, T. M. Shata, T. Archampong, A. Kwara, M. T. Aliota, and J. T. Blackard, “Zika virus exposure in an HIV-infected cohort in Ghana,” 2018.
  * (14) E. C. João, O. da C. Ferreira, M. I. Gouvêa, M. de Lourdes B. Teixeira, A. Tanuri, L. M. Higa, D. A. Costa, R. Mohana-Borges, M. B. Arruda, H. J. Matos, M. L. Cruz, W. Mendes-Silva, and J. S. Read, “Pregnant women co-infected with HIV and Zika: Outcomes and birth defects in infants according to maternal symptomatology,” PLoS ONE, vol. 13, 2018.
  * (15) C. D. Zorrilla, A. M. Mosquera, S. Rabionet, and J. Rivera-Viñas, “HIV and Zika in pregnancy: parallel stories and new challenges,” Obstetrics & Gynecology International Journal, vol. 5, no. 6, 2016.
  * (16) P. Brasil, J. P. Pereira, M. E. Moreira, R. M. R. Nogueira, L. Damasceno, M. Wakimoto, R. S. Rabello, S. G. Valderramos, U.-A. Halai, T. S. Salles, A. A. Zin, D. Horovitz, P. Daltro, M. Boechat, C. R. Gabaglia, P. C. de Sequeira, J. H. Pilotto, R. Medialdea-Carrera, D. C. da Cunha, L. M. A. de Carvalho, M. Pone, A. M. Siqueira, G. A. Calvet, A. E. R. Baião, E. S. Neves, P. R. N. de Carvalho, R. H. Hasue, P. B. Marschik, C. Einspieler, C. Janzen, J. D. Cherry, A. M. B. de Filippis, and K. Nielsen-Saines, “Zika virus infection in pregnant women in Rio de Janeiro,” New England Journal of Medicine, vol. 375, 2016.
  * (17) A. Aschengrau, M. M. Mussi-Pinhata, J. Moye, N. Chakhtoura, K. Patel, P. L. Williams, B. Karalius, P. A. Garvie, D. Monte, F. Whalen, J. Lebov, and G. R. Seage, “An international prospective cohort study of HIV and Zika in infants and pregnancy (HIV-ZIP): Study protocol,” Frontiers in Global Women’s Health, vol. 2, 2021.
  * (18) J. P. Romero-Leiton, J. E. Castellanos, and E. Ibargüen-Mondragón, “An optimal control problem and cost-effectiveness analysis of malaria disease with vertical transmission applied to San andrés de Tumaco (Colombia),” Computational and Applied Mathematics, vol. 38, 2019.
  * (19) P. V. D. Driessche and J. Watmough, “Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission,” Mathematical Biosciences, vol. 180, 2002.
  * (20) S. Bañuelos, M. Martinez, C. Mitchell, and A. Prieto-Langarica, “Using mathematical modelling to investigate the effect of the sexual behaviour of individuals and vector control measures on Zika,” Letters in Biomathematics, vol. 6, no. 1, pp. 1–19, 2019.
  * (21) N. Chitnis, J. M. Hyman, and J. M. Cushing, “Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model,” Bulletin of Mathematical Biology, vol. 70, 2008.
  * (22) L. S. Pontryagin, V. Boltyanskii, R. Gamkrelidze, E. Mishchenko, K. Trirogoff, and L. Neustadt, LS Pontryagin Selected Works: The Mathematical Theory of Optimal Processes.  Routledge, 2018.
  * (23) B. Joob and V. Wiwanitkit, “Concurrent HIV and Zika virus infection, a problem to be address in HIV medicine,” Int J Virol AIDS, vol. 3, p. 035, 2017.
  * (24) M. R. Bidokhti, D. Dutta, L. S. Madduri, S. M. Woollard, R. Norgren, L. Giavedoni, and S. N. Byrareddy, “SIV/SHIV-Zika co-infection does not alter disease pathogenesis in adult non-pregnant rhesus macaque model,” PLoS Neglected Tropical Diseases, vol. 12, 2018.
  * (25) K. Prieto and J. P. Romero–Leiton, “Current forecast of HIV/AIDS using Bayesian inference,” Natural Resource Modeling, vol. 34, 2021.
  * (26) D. Gao, Y. Lou, D. He, T. C. Porco, Y. Kuang, G. Chowell, and S. Ruan, “Prevention and control of zika as a mosquito-borne and sexually transmitted disease: a mathematical modeling analysis,” Springer, 2016.
  * (27) S. K. Biswas, U. Ghosh, and S. Sarkar, “Mathematical model of Zika virus dynamics with vector control and sensitivity analysis,” Infectious Disease Modelling, vol. 5, pp. 23–41, 2020.
  * (28) W. H. Organization et al., “Zika strategic response plan quarterly update,” tech. rep., World Health Organization, 2016.
  * (29) L. J. Forero-Martínez, R. Murad, M. Calderón-Jaramillo, and J. C. Rivillas-García, “Zika and women’s sexual and reproductive health: Critical first steps to understand the role of gender in the colombian epidemic,” International Journal of Gynecology and Obstetrics, vol. 148, 2020.
  * (30) C. Prachniak-Rincón and J. V. D. Onís, “HIV and the right to health in Colombia,” Health and Human Rights, vol. 18, 2016.
  * (31) D. Bancroft, G. M. Power, R. T. Jones, E. Massad, J. B. Iriat, R. Preet, J. Kinsman, and J. G. Logan, “Vector control strategies in Brazil: A qualitative investigation into community knowledge, attitudes and perceptions following the 2015-2016 Zika virus epidemic,” BMJ Open, vol. 12, 2022\.
  * (32) A. S. Benzaken, G. F. Pereira, L. Costa, A. Tanuri, A. F. Santos, and M. A. Soares, “Antiretroviral treatment, government policy and economy of HIV/AIDS in Brazil: Is it time for HIV cure in the country?,” 2019.
  * (33) J. Hainzl, M. D. Leko, A. Pignedoli, and B. Vrdoljak, 12: Ordinary Differential Equations and Dynamical Systems.  2022\.
  * (34) C. Castillo-Chavez and H. R. Thieme, “Asymptotically autonomous epidemic models,” Science, vol. 94, 1994.
  * (35) C. C. McCluskey and J. S. Muldowney, “Bendixson-Dulac criteria for difference equations,” Journal of Dynamics and Differential Equations, vol. 10, 1998.
  * (36) E. Roxin and D. L. Lukes, “Differential equations: Classical to controlled.,” The American Mathematical Monthly, vol. 92, 1985.

##  Appendix A Appendices

###  A.1 Proof that $\mathbf{E}_{h_{0}}$ is GAS when $\mathcal{R}_{h}<1$

Let X be the vector field given by the right side hand of System (A.2). It is
enough to prove the existence of a Lyapunov function for the translated system

$\dot{\textbf{X}}=F(\textbf{X}+\textbf{E}_{h_{0}})-F(\textbf{E}_{h_{0}})=f(\textbf{X}),$

where the system $\dot{\textbf{Y}}=F(\textbf{Y})$ has $\textbf{Y}=0$ as
equilibrium solution. Let us consider the function

$V^{*}(S,I_{h},A)=\dfrac{1}{\sigma_{1}+\mu}I_{h},$

and define

$V(\bar{S},\bar{I}_{h},\bar{A})=V^{*}\left(S-\frac{\Lambda}{\mu},I_{h},A\right).$
(A.1)

We verify that the function $V$ defined on (A.1) is a Lyapunov function.
Indeed, $V$ is positive definite, that is,
$V(\textbf{E}_{h_{0}})=V^{*}(\textbf{0})=0$ and $V>0$ for all
$(\bar{S},\bar{I}_{h},\bar{A})\neq\textbf{E}_{h_{0}}in\,\Omega_{h}$.
Additionally, the orbital derivative of $V$ along the trajectories of system
(2.2) is

$\begin{array}[]{ll}\dot{V}&=\dfrac{\partial
V^{*}}{\partial\left(S-\frac{\Lambda}{\mu}\right)}\dot{S}+\dfrac{\partial
V^{*}}{\partial I_{h}}\dot{I}_{h}+\dfrac{\partial V^{*}}{\partial A}\dot{A}\\\
\\\
&=\dfrac{1}{\sigma_{1}+\mu}\left(\beta_{h}I_{h}S-(\sigma_{1}+\mu)I_{h}\right)\\\
\\\
&\leq\dfrac{1}{\sigma_{1}+\mu}(\sigma_{1}+\mu)\left(\dfrac{\beta_{h}\frac{\Lambda}{\mu}}{\sigma_{1}+\mu}-1\right)I_{h}\\\
\\\ &=(\mathcal{R}_{h}-1)I_{h}\leq 0,\quad\forall I_{h}\geq 0.\end{array}.$

Let
$\triangle=\\{(S,I_{h},A)\in\mathbb{R}_{+}^{3}:\dot{V}=0\\}\subset\\{(S,I_{h},A)\in\mathbb{R}_{+}^{3}:I_{h}=0\\}$
and $\triangle^{\prime}\subset\triangle$ the largest invariant set respect to
(2.2). It can be easily showed that
$\triangle^{\prime}=\\{\textbf{E}_{h_{0}}\\}$. Therefore by the LaSalle’s
invariance principle teschl2012ordinary , $\textbf{E}_{h_{0}}$ is a global
attractor whenever $\mathcal{R}_{h}<1$.

###  A.2 Proof that $\mathbf{E}_{h}^{*}$ is GAS when $\mathcal{R}_{h}>1$

Note that the third equation of System (2.2) is uncoupled in the variable $S$
and its equilibrium is $I_{h}=\frac{\mu_{h}+\mu}{\sigma_{1}}A$. Replacing this
value in the two first equations, we obtain the planar system

$\left\\{\begin{array}[]{ll}\dfrac{dS}{dt}&=\Lambda-\beta_{h}\dfrac{\mu_{h}+\mu}{\sigma_{1}}SA-\mu
S\\\ \\\ \dfrac{dA}{dt}&=\beta_{h}AS-(\sigma_{1}+\mu)A.\end{array}\right.$
(A.2)

Therefore, solutions to system (2.2) tend asymptotically to those of the
planar system (A.2) (see e.g. castillo1994asymptotically ). The Dulac
criterion mccluskey1998bendixson  claims that if there exists a real
continuously differentiable function $\phi(S,A)$ such that
$\nabla\cdot[\phi(S,A)\textbf{X}(S,A)]\neq 0$, where
$\textbf{X}(S,A)=(F(S,A),G(S,A)$ is the right side hand of system (A.2), then
there are no periodic orbits contained entirely inside $\Omega_{h}$. Let

$\phi(S,A)=\dfrac{1}{SA}\quad\text{for}\quad S>0,\;A>0,$

then

$\begin{array}[]{ll}\nabla\cdot[\phi(S,A)\textbf{X}(S,A)]&=\dfrac{\partial(F\phi)}{\partial
S}+\dfrac{\partial(G\phi)}{\partial A}\\\ \\\ &=\dfrac{\partial}{\partial
S}\left(\frac{\Lambda}{SA}-\frac{\beta_{h}(\mu_{h}+\mu)}{\sigma_{1}SA}SA-\frac{\mu
S}{SA}\right)+\dfrac{\partial}{\partial
A}\left(\frac{\beta_{h}AS}{SA}-\dfrac{\sigma_{1}+\mu}{SA}A\right)\\\ \\\
&=\dfrac{\partial}{\partial
S}\left(\frac{\Lambda}{SA}-\frac{\beta_{h}(\mu_{h}+\mu)}{\sigma_{1}}-\frac{\mu}{A}\right)+\dfrac{\partial}{\partial
A}\left(\beta_{h}-\frac{\sigma_{1}+\mu}{S}\right)\\\ \\\
&=-\dfrac{\Lambda}{AS^{2}}<0,\,\text{for}\,S,A>0.\end{array}.$

Thus, there are no periodic orbits in $\Omega_{h}$. Given that $\Omega_{h}$ is
positively invariant, and the endemic equilibrium exists if
$\mathcal{R}_{h}>1$, it follows from the Poincaré-Bendixson Theorem
mccluskey1998bendixson  that all solutions of the system starting in
$\Omega_{h}$ remain in $\Omega_{h}$ for all $t$. Thus, because of the absence
of periodic orbits in $\Omega_{h}$, this implies that the unique endemic
equilibrium of System (2.2) is GAS when $\mathcal{R}_{h}>1$.

###  A.3 Proof of Lemma 1

We prove item $i)$; all others are proved in a similar way. On the one hand,
from the definition of $\mathcal{R}_{z}^{*}$, we immediately have that if
$\mathcal{R}_{z}^{*}=2\mathcal{R}_{z_{1}}+\bar{\mathcal{R}}_{z_{2}}<1$, then
$2\mathcal{R}_{z_{1}}<1.$ On the other hand,

$\begin{array}[]{ll}\mathcal{R}_{z}^{*}<1&\Rightarrow\mathcal{R}_{z}^{2}+2\mathcal{R}_{z_{1}}(1-\mathcal{R}_{z})<1\\\
\\\ &\Rightarrow\mathcal{R}_{z}(1-2\mathcal{R}_{z_{1}})+2\mathcal{R}_{z}\\\
\\\
&\Rightarrow\mathcal{R}_{z}(1-2\mathcal{R}_{z_{1}})<1-2\mathcal{R}_{z},\quad(2\mathcal{R}_{z_{1}}<1)\\\
\\\ &\Rightarrow\mathcal{R}_{z}<1.\end{array}$

###  A.4 Proof of the condition $R^{*}<R_{max}$

We have to prove that $R^{*}<R_{max}$ is satisfied. To this end, we consider
the polynomial $p(R)$ associated to system (2.20), whose coefficients are in
equation (2.21). The graph of $p(R)$ is a parabola that opens upward (see
Figure A.1). Since $R^{*}>0$ and $p(0)=c<0$ as long as
$\mathcal{R}_{z_{1}}+\mathcal{R}_{z_{2}}>1$, then to ensure that
$R^{*}<R_{max}$, it is enough to prove that $p(R_{max})>0$. In fact

$\begin{array}[]{ll}p(R_{max})&=aR_{max}^{2}+bR_{max}+c\\\ \\\
&=\dfrac{\beta_{z}\Lambda\alpha_{m}\mu}{\delta_{z}R_{max}}R_{max}^{2}+\left[\left(\dfrac{\beta_{m}\alpha_{m}\Lambda_{m}}{\mu_{m}}+\beta_{z}\mu_{m}\right)\dfrac{\Lambda}{R_{max}}+\left[\alpha_{m}\mu(\mu_{z}+\delta_{z}+\mu)-\beta_{z}\Lambda\alpha_{m}\right]\dfrac{\mu}{\delta_{z}}\right]R_{max}\\\
\\\ &+(\mu_{z}+\delta_{z}+\mu)\mu\mu_{m}\left(1-\mathcal{R}_{z}^{*}\right)\\\
\\\
&=\dfrac{\beta_{z}\Lambda\alpha_{m}\mu}{\delta_{z}}R_{max}+\left(\dfrac{\beta_{m}\alpha_{m}\Lambda_{m}}{\mu_{m}}+\beta_{z}\mu_{m}\right)\Lambda+\left[\alpha_{m}\mu(\mu_{z}+\delta_{z}+\mu)-\beta_{z}\Lambda\alpha_{m}\right]\dfrac{\mu}{\delta_{z}}R_{max}\\\
\\\ &+(\mu_{z}+\delta_{z}+\mu)\mu\mu_{m}\left(1-\mathcal{R}_{z}^{*}\right)\\\
\\\
&=\left(\dfrac{\beta_{m}\alpha_{m}\Lambda_{m}}{\mu_{m}}+\beta_{z}\mu_{m}\right)\Lambda+\dfrac{\alpha_{m}\mu^{2}(\mu_{z}+\delta_{z}+\mu)R_{max}}{\delta_{z}}+(\mu_{z}+\delta_{z}+\mu)\mu\mu_{m}\left(1-\mathcal{R}_{z}^{*}\right)\\\
\\\
&=(\mu_{z}+\delta_{z}+\mu)\mu\mu_{m}(\mathcal{R}_{z}^{*})(\mathcal{R}_{z}^{*}-\mathcal{R}_{z}^{*})+\alpha_{m}\Lambda\mu^{2}+(\mu_{z}+\delta_{z}+\mu)\mu\mu_{m}\\\
\\\
&=\alpha_{m}\Lambda\mu^{2}+(\mu_{z}+\delta_{z}+\mu)\mu\mu_{m}>0.\end{array}$

Figure A.1: Polynomial $p(R)$ associated to system (2.20).

###  A.5 Proof that $\textbf{E}_{z_{0}}$ is GAS when $\mathcal{R}_{z}^{*}<1$

We will use $\kappa=\mu_{z}+\delta_{z}+\mu$ and the notation and reasoning
analogous to that used in the proof of Proposition LABEL:global_DFE-hiv, with
functions

$V^{*}(S,I_{z},R,S_{m},I_{m})=\left(S-\dfrac{\Lambda}{\mu}\log{\dfrac{S\mu}{\Lambda}}\right)+I_{z}+\dfrac{\beta_{m}\Lambda}{\mu\mu_{m}}\left(S_{m}-\dfrac{\Lambda_{m}}{\mu_{m}}\log{\dfrac{S_{m}\mu_{m}}{\Lambda_{m}}}\right)+\dfrac{\beta_{m}\Lambda}{\mu\mu_{m}}I_{m}$

and

$V(\bar{S},\bar{I}_{z},\bar{R},\bar{S}_{m},\bar{I}_{m})=V^{*}\left(S-\frac{\Lambda}{\mu},I_{z},R,S_{m},I_{m}-\frac{\Lambda_{m}}{\mu_{m}}\right).$
(A.3)

We must prove that $V$ defined in (A.3) is a Lyapunov function. $V$ is
positive definite and the orbital derivative of $V$ along the trajectories of
(2.9) is

$\begin{array}[]{ll}\dot{V}&=\left(1-\frac{\Lambda}{\mu
S}\right)\dot{S}+\dot{I}_{z}+\frac{\beta_{m}\Lambda}{\mu\mu_{m}}\left(1-\frac{\Lambda_{m}}{\mu_{m}S_{m}}\right)\dot{S}_{m}+\frac{\beta_{m}\Lambda}{\mu\mu_{m}}\dot{I}_{m}\\\
\\\ &=\left(1-\frac{\Lambda}{\mu
S}\right)(\Lambda-\beta_{m}I_{m}S-\beta_{z}I_{z}S-\mu
S)+(\beta_{m}I_{m}S+\beta_{z}I_{z}S-\kappa I_{z})\\\ \\\
&+\frac{\beta_{m}\Lambda}{\mu\mu_{m}}\left(1-\frac{\Lambda_{m}}{S_{m}\mu_{m}}\right)(\Lambda_{m}-\alpha_{m}I_{z}S_{m}-\mu_{m}S_{m})+\frac{\beta_{m}\Lambda}{\mu\mu_{m}}(\alpha_{m}I_{z}S_{m}-\mu_{m}I_{m})\\\
\\\ &=(\Lambda-\mu S)-\frac{\Lambda}{\mu S}(\Lambda-\mu S)+\frac{\Lambda}{\mu
S}(\beta_{m}I_{m}+\beta_{z}I_{z})S-\kappa
I_{z}+\frac{\beta_{m}\Lambda}{\mu\mu_{m}}(\Lambda_{m}-\mu_{m}S_{m})\\\ \\\
&-\frac{\beta_{m}\Lambda\Lambda_{m}}{\mu\mu_{m}S_{m}\mu_{m}}(\Lambda_{m}-\mu_{m}S_{m})+\frac{\beta_{m}\Lambda\Lambda_{m}\alpha_{m}S_{m}I_{z}}{\mu\mu_{m}S_{m}\mu_{m}}-\frac{\beta_{m}\Lambda\mu_{m}I_{m}}{\mu\mu_{m}}\\\
\\\ &=-\frac{(\Lambda-\mu S)^{2}}{\mu
S}-\frac{\beta_{m}\Lambda}{\mu\mu_{m}^{2}}\frac{(\Lambda_{m}-\mu_{m}S_{m})^{2}}{S_{m}}+\left(\frac{\Lambda}{\mu}\beta_{z}+\frac{\beta_{m}\alpha_{m}\Lambda_{m}\Lambda}{\mu\mu_{m}^{2}}-\kappa\right)I_{z}\\\
\\\ &=-\frac{(\Lambda-\mu S)^{2}}{\mu
S}-\frac{\beta_{m}\Lambda}{\mu\mu_{m}^{2}}\frac{(\Lambda_{m}-\mu_{m}S_{m})^{2}}{S_{m}}+\kappa\left(\mathcal{R}_{z}^{*}-1\right)I_{z}.\end{array}$

Note that the last expression in the above inequality is negative if
$\mathcal{R}_{z}^{*}<1$ and for $S=\frac{\Lambda}{\mu}$,
$S_{m}=\frac{\Lambda_{m}}{\mu_{m}}$ and $I_{z}=I_{m}=0$. Finally, using the
LaSalle’s invariance principle teschl2012ordinary , $\textbf{E}_{z_{0}}$
defined in (2.11) is a global attractor whenever $\mathcal{R}_{z}^{*}<1$.

###  A.6 Sensitivity index of $\mathcal{R}_{z}$ with respect to $\beta_{z}$

We have

$\begin{array}[]{ll}\dfrac{\partial\mathcal{R}_{z}}{\partial\beta_{z}}&=\dfrac{\partial}{\partial\beta_{z}}\left(\dfrac{\beta_{z}\Lambda}{2\mu\kappa_{1}}+\sqrt{\left(\dfrac{\beta_{z}\Lambda}{2\mu\kappa_{1}}\right)^{2}+\dfrac{\alpha_{m}\beta_{m}\Lambda_{m}\Lambda}{\mu\mu_{m}^{2}\kappa_{1}}}\right)\\\
\\\
&=\dfrac{\Lambda}{2\mu(\delta_{z}+\mu+\mu_{z})}+\dfrac{\beta_{z}\Lambda^{2}}{4\mu^{2}(\delta_{z}+\mu+\mu_{z})^{2}\sqrt{\dfrac{\beta_{z}^{2}\Lambda^{2}}{4\mu^{2}(\delta_{z}+\mu+\mu_{z})^{2}}+\dfrac{\alpha_{m}\beta_{m}\Lambda\Lambda_{m}}{\mu\mu_{m}^{2}(\delta_{z}+\mu+\mu_{z})}}}.\end{array}$

Then, it is enough to compute the expression for

$\dfrac{\partial\mathcal{R}_{z}}{\partial\beta_{z}}\dfrac{\beta_{z}}{\mathcal{R}_{z}}=\dfrac{\partial\mathcal{R}_{z}}{\partial\beta_{z}}\dfrac{\beta_{z}}{\dfrac{\beta_{z}\Lambda}{2\mu\kappa_{1}}+\sqrt{\left(\dfrac{\beta_{z}\Lambda}{2\mu\kappa_{1}}\right)^{2}+\dfrac{\alpha_{m}\beta_{m}\Lambda_{m}\Lambda}{\mu\mu_{m}^{2}\kappa_{1}}}}.$

###  A.7 Proof of Proposition 7

All state variables and controls are non-negative and, for $i=\\{1,2,3\\}$,
the set of control variables $\eta_{i}\in\mathcal{A}$ is also convex and
closed. We note that the boundedness of the optimal system (4.1) determines
the compactness for the existence of the optimal control. Moreover, there
exists a constant $\nu>1$, $\omega_{1}=\min(d_{1},d_{2},d_{3})$, and
$\omega_{2}>0$ such that

$\mathcal{J}(\eta)\geq\omega_{1}||\eta||^{\nu}-\omega_{2}.$ (A.4)

Therefore, according to lukes1982differential , the controlled system (4.1)
admits an optimal control solution $\eta^{*}$.

###  A.8 Proof of Proposition 8

We have

$\displaystyle\dfrac{p_{1}}{dt}$ $\displaystyle=-\dfrac{\partial H}{\partial
S}=$ $\displaystyle
p_{1}\left[(1-\eta_{1})\beta_{m}I_{m}+(1-\eta_{3})(\beta_{z}I_{z}+\beta_{h}I_{h})+\mu\right]-p_{2}\left[(1-\eta_{1})\beta_{m}I_{m}+(1-\eta_{3})\beta_{z}I_{z}\right]$
$\displaystyle-p_{3}\left[(1-\eta_{3})\beta_{h}I_{h}\right],$
$\displaystyle\dfrac{p_{2}}{dt}$ $\displaystyle=-\dfrac{\partial H}{\partial
I_{z}}=$ $\displaystyle-
c_{1}+p_{1}(1-\eta_{3})\beta_{z}S-p_{2}\left[(1-\eta_{3})\beta_{z}S-\omega_{2}(1-\eta_{3})\beta_{h}I_{h}-(\mu_{z}+\delta_{z}+\mu)\right]-p_{3}(1-\eta_{3})\beta_{z}I_{h}$
$\displaystyle-
p_{6}\delta_{z}+p_{7}(1-\eta_{1})\alpha_{m}S_{m}-p_{8}(1-\eta_{1})\alpha_{m}S_{m},$
$\displaystyle\dfrac{p_{3}}{dt}$ $\displaystyle=-\dfrac{\partial H}{\partial
I_{h}}=$ $\displaystyle-
c_{2}+p_{1}(1-\eta_{3})\beta_{h}S+p_{2}\omega_{2}(1-\eta_{3})\beta_{h}I_{z}-p_{3}\left[\omega_{1}[(1-\eta_{1})\beta_{m}I_{m}+(1-\eta_{3})\beta_{z}I_{z}]\right.$
$\displaystyle\left.-(1-\eta_{2})\sigma_{1}-\mu\right]-p_{4}\left[\omega_{1}[(1-\eta_{1})\beta_{m}I_{m}+(1-\eta_{3})\beta_{z}I_{z}]+\omega_{2}(1-\eta_{3})\beta_{h}I_{z}\right]$
$\displaystyle-p_{5}(1-\eta_{2})\sigma_{1},$ $\displaystyle\dfrac{p_{4}}{dt}$
$\displaystyle=-\dfrac{\partial H}{\partial I_{hz}}=$ $\displaystyle-
c_{3}+p_{4}\left(\epsilon\delta_{z}+(1-\eta_{2})\sigma_{2}+\mu_{hz}+\mu\right)-p_{5}(1-\eta_{2})\sigma_{2}-p_{6}\epsilon\delta_{z}+p_{7}(1-\eta_{1})\alpha_{m}S_{m}$
$\displaystyle-p_{8}(1-\eta_{1})\alpha_{m}S_{m},$
$\displaystyle\dfrac{p_{5}}{dt}$ $\displaystyle=-\dfrac{\partial H}{\partial
A}=$ $\displaystyle-c_{5}+p_{5}(\mu_{h}+\mu),$
$\displaystyle\dfrac{p_{6}}{dt}$ $\displaystyle=-\dfrac{\partial H}{\partial
R}=$ $\displaystyle p_{6}\mu,$ $\displaystyle\dfrac{p_{7}}{dt}$
$\displaystyle=-\dfrac{\partial H}{\partial S_{m}}=$ $\displaystyle
p_{7}\left((1-\eta_{1})\alpha_{m}(I_{z}+I_{hz})+\mu_{m}\right)-p_{8}(1-\eta_{1})\alpha_{m}(I_{z}+I_{hz}),$
$\displaystyle\dfrac{p_{8}}{dt}$ $\displaystyle=-\dfrac{\partial H}{\partial
I_{m}}=$ $\displaystyle-
c_{4}+p_{1}(1-\eta_{1})\beta_{m}S-p_{2}(1-\eta_{1})\beta_{m}S+p_{3}(1-\eta_{1})\beta_{m}I_{h}-p_{4}\omega_{1}(1-\eta_{1})\beta_{m}I_{h}+p_{8}\mu_{m},$

with transversality conditions $p_{i}(T)=0$, for $i=\\{1,2,3,4,5,6,7,8\\}$.
According to PMP, the optimal conditions are

$\displaystyle\dfrac{\partial H}{\partial\eta_{1}}$ $\displaystyle=$
$\displaystyle
d_{1}\eta_{1}-\left((p_{2}-p_{1})S+(p_{4}-p_{3})I_{h}\right)\beta_{m}I_{m}-(p_{8}-p_{7})\alpha_{m}(I_{z}+I_{hz})S_{m}=0,$
$\displaystyle\dfrac{\partial H}{\partial\eta_{2}}$ $\displaystyle=$
$\displaystyle
d_{2}\eta_{2}-\left(p_{5}-p_{3}\right)\sigma_{1}I_{h}-\left(p_{5}-p_{4}\right)\sigma_{2}I_{hz}=0,$
$\displaystyle\dfrac{\partial H}{\partial\eta_{3}}$ $\displaystyle=$
$\displaystyle
d_{3}\eta_{3}-(p_{2}-p_{1})\beta_{z}I_{z}S-(p_{3}-p_{1})\beta_{h}I_{h}S-(p_{4}-p_{3})\omega_{1}\beta_{z}I_{z}I_{h}-(p_{4}-p_{2})\omega_{2}\beta_{h}I_{z}I_{h}=0.$

Hence, we get assertions (4.6). This completes the proof.