Mathematical modelling of the impact of testing, treatment and control of HIV transmission in Kenya

Abstract

The control of spread of HIV to reduce its effects on a population is an important role of public health. HIV testing and counselling (HTC) and eventual enrolment of infected individuals on anti retro-viral treatment (ART) as soon as possible to reduce the risk of dying is currently the main intervention against HIV. Mathematical models can be used to study the effects of HIV prevention, testing and treatment with ART on HIV patients. In this study, we employ a deterministic model to provide a quantification of HIV prevention, testing and treatment with ART as public health measurements in the fight against HIV infection. Lyapunov function has been used to derive a condition that ensures that the model system is globally asymptotically stable when R₀ is less than unity. Through sensitivity analysis, we determine the relative importance of model parameters for disease transmission. The sensitivity analysis results suggest that the effective contact rates are mechanisms fuelling HIV epidemic proliferation while ART efficacy reduces the incidence. The model is fitted to HIV surveillance data obtained from world in data website. Although the results show a high proportion of individuals with HIV in Kenya, the incidence curve is indicative of a declining HIV infection and settling at an endemic steady state. The results are suggestive of the need to promote preventive mechanism against the occurrence of new infections. Moreover, the results show that the combination of several control mechanisms would significantly reduce the spread of the disease, if we maintain the level of each control high.

Keywords:

• Human Immune-deficiency Virus (HIV)
• Anti Retro-viral Therapy (ART)
• basic reproduction number
• global stability
• sensitivity analysis
• simulation

PUBLIC INTEREST STATEMENT

This paper provides a deterministic model usable to study the dynamics of HIV and its treatment geared towards attaining a HIV free world by year 2030. The dynamics are clearly unique and novel. It provides knowledge on how effective a combination of control strategies is in limiting cases of HIV infections. Furthermore, sensitivity analysis shows that the epidemic is driven by one on one contact between the infected population and the population at high risk of infection. It can be argued that effective PrEP up-take and high ART efficacy has the potential of reducing the HIV infection in the population.

1. Introduction

Diagnosis of HIV infection is a prerequisite for its treatment. This begins from voluntary counselling and testing which refers to the process initiated by an individual who wishes to know his/her HIV status. According to the WHO regulations, all forms of HIV testing and counselling should be voluntary and adhere to the five Cs, namely consent, confidentiality, counselling, correct test results and connection to care, treatment and prevention services (WHO, 2014). Identification of recent HIV infections and their role in driving the HIV epidemics is critical in understanding the spread of HIV and the deployment of relevant control measures (Quinn et al., 2000; Yerly et al., 2001). According to NACC (2014a), an estimated 53% of the 1.6 million people living with HIV in Kenya are not aware of their HIV status. In addition, approximately 260,000 couples in HIV sero-discordant couples significantly contribute to new infections.

HTC has been adopted through targeted community-based testing and door-to-door testing initiatives by the Kenya government in order to reduce the spread of HIV. In order to boost these efforts, the government of Kenya announced plans to introduce self-test kits in 2015 (UNAIDS, 2016). In 2008, approximately 860,000 people were being tested annually for HIV totalling to about 6.4 million people by the end of 2013 (NACC, 2014b). Furthermore, the report from NACC (2014b) indicates that HTC coverage among the pregnant women has risen substantially. Between 2009 and 2013, the number of pregnant women tested for HIV increased from 68% to 92%. Similarly, there has been a concerted effort to increase testing rates among Kenyan men, with community-based testing programmes proving particularly successful.

In order to aid the reduction of new HIV infections in Kenya, in 2016, the government issued full regulatory approval of pre-exposure prophylaxis (PrEP) thus becoming the second country in Sub-Saharan Africa to make such an approval after South Africa (UNAIDS, 2016). Pre-Exposure Prophylaxis (PrEP) is the prophylactic use of anti retro-viral drugs (ARVS) by people who do not have HIV but are at high risk of acquiring it to prevent HIV infection (AIDS, 2016; HIV/AIDS, 2016). However, successes are yet to be reported since research into the uptake and impact of PrEP, specifically with young women and girls in high-incidence areas, is still on-going (UNAIDS, 2016).

The major intervention for the HIV infection is that when people are infected with HIV they should start anti retro-viral treatment (ART) as soon as possible to reduce the risk of dying (Williams, 2014; Williams, Hargrove, & Humphrey, 2010). The fact that ART has the potential to stop transmission of HIV, United Nations Programme on HIV and AIDS (UNAIDS) called for an end to HIV/AIDS centred on 90–90–90 target by 2020 (United Nations Programme on HIV/AIDS, 2014). The strategy implies that 90% of all those living with HIV should have been tested within the last one year, unless they already know that they are infected with HIV, 90% of these should be on treatment and 90% of these should have plasma viral loads below 1,000 copies/mL. In line with this strategy, Kenya began to adopt 2015 World Health Organization (WHO) recommendations to immediately offer treatment to people diagnosed with HIV in order to increase the ART access. It is estimated that about 826,000 adults and 71,500 children were receiving ART treatment in 2015 (UNAIDS, 2016). This paper presents a deterministic model for predicting the epidemiological trends of HIV that exploits HIV surveillance data to model the disease evolution in Kenya. Furthermore, we present results showing the impact of combination of HTC, PrEP and ART treatment as HIV intervention measures.

The rest of the paper is organised as follows. In Section 2 a brief literature related to this work is reviewed. In Section 3 a deterministic model is presented while in Section 4 detailed mathematical analysis is provided. Numerical simulations are carried out in Section 5. Finally, Section 6 gives the conclusion and suggestions for future work.

2. Related literature

Mathematical modelling is an important tool in understanding the spread and control of HIV. Early diagnosis and immediate treatment of HIV infection as well as suppression of the viral load remain the key interventions in reducing the HIV incidence. To understand the severity of the HIV epidemics as well as establishing a baseline, Su et al. (2016) constructed a deterministic transmission model of HIV. The model was calibrated using surveillance and treatment data for the period 2005–2008. The authors then validated the model by comparing its predicted value of HIV prevalence in 2010 to the prevalence data of 2010. The validated model produced estimations on the new infections, people living with HIV (PLWHIV) and HIV-AIDS related deaths. They established that HIV control programs have drastically scaled up the HIV testing and treatment.

Kok et al. (2015) constructed a dynamic model of the continuum of HIV care in Vancouver, Canada. This model incorporated the main activities and decisions in the delivery of anti retro-viral therapy that included HIV testing, linkage to care and long-term retention in care and treatment. In the development of the model, the authors considered four sub-populations namely men who have sex with men (MSM), injection drug users, female sex workers and the general population. The model was fitted to data obtained from public health sources that included British Columbia Centre for Excellence in HIV/AIDS. The findings predicted a lower than expected number of new diagnoses among MSM. Furthermore, their findings suggested that optimal resource allocation favours routine testing in high prevalence settings over targeted testing and that a greater impact would be achieved by allocating more resources to routine testing in high prevalence settings for MSM. Okosun, Makinde, and Takaidza (2013) used a deterministic model based on ordinary differential equations to investigate the impact of combining prevention intervention strategies on the spread of HIV in the presence of unaware infectives. Thiébaut and May (2013) highlight the use of mathematical models to predict the impact of early antiretroviral treatment (ART) of HIV-infected individuals and pre-exposure prophylaxis (PrEP) of uninfected individuals.

Pedro and Tchuenche (2010) developed a mathematical model to study the effect of socio-economic classes in the transmission and spread of HIV. The results from this model established that men with better livelihood tend to have higher female sexual partners from the poor settings leading to high HIV incidence as opposed to the low income earners. He further argued that targeted HIV information campaign was needed to encourage multiple partners from the well-off settings to regularly or occasionally observe sexual hygiene. However, this study did not include nor investigate the potential impact of control strategies that are targeting the reduction of HIV incidence. On the other hand, Yang, Shu, Lam, and Sun (2017) performed global analysis on the dynamics of HIV in female sex workers by factoring the effect of senior male clients. Graphical representation of their results showed that an increase in the number of male clients to the female sex workers results in an increase in the generation of new cases of infections. The major weakness of this model was failure to incorporate health education campaign and other treatment control strategies. The assessment of these strategies is essential in understanding the HIV dynamics and establishing how effective they can limit its spread to realise a HIV free population. Other modelling work that are pertinent to this study include Birger, Hallett, Sinha, Grenfell, and Hodder (2013); Lou et al. (2014); Luo et al. (2015); Sun, Xiao, Peng, and Wang (2013); Zhang, Gray, and Wilson (2012).

None of the above mentioned models in the literature take into account the combination of HIV intervention strategies such as HTC, PrEP and ART treatment as HIV. Stopping the transmission chain of HIV remains the only viable form of control against the spread of infection. The potential effect of testing and immediate treatment with ART is thus of great interest. In the current paper, we use a mathematical model to describe the transmission dynamics of HIV, the impact of testing and treatment on the disease transmission. As mentioned previously, in our model we consider surveillance data to describe the spread of HIV. We assume a scenario in which the anti retro-viral therapy (ART) is initiated using the World Health Organization (WHO) viral load suppression strategy centred on 90-90-90 target by 2020. This is to imply that all new HIV infected individuals are immediately linked to treatment and retained in the system through regular testing.

3. Mathematical model

3.1. Model formulation

According to WHO (2015) ART guidelines, the strategy of HIV treatment changed from the CD4 cell count $\le$ 350 cells/μL for eligibility to enrol for ART treatment to achieve viral load suppression through targeted testing and immediate enrolment into ART treatment to achieve the 90-90-90 plan. A deterministic model is proposed to represent the effect of HIV testing, uptake of PrEP and ART as HIV intervention method. The model comprises susceptible class (S) and infection class. The infection class is further divided into four compartments. The undiagnosed class (I) comprising individuals who are unaware of their HIV status. These individuals are tested and then grouped into two classes based on the CD4 cell counts. These include individuals with CD4 cell count $\ge 350$ cells/μL considered as HIV infected individuals with lower plasma viral load $(I_1)$ and individuals with CD4 cell count $<350$ cells/μL considered as HIV infected patients with high plasma viral load $(I_2)$. These diagnosed patients move into treatment classes $I_1^T$ and $I_2^T$, respectively. We consider a population of size $N(t)$ at time $t$ given in (1)

(1) $N(t)=S(t)+I(t)+I_1(t)+I_2(t)+I_1^T(t)+I_2^T(t),$(1)

The recruitment of susceptible is assumed to occur at a constant rate $\mathrm{\Pi }$. During this modelling time, we only consider individuals aged 15 years and above. This group is assumed to be adult population and hence the recruitment rate defines the demographic process of individuals attaining the age 15. The rate at which susceptible individuals acquire infection is given by

(2) $\begin{array}{c}\lambda =\left(\frac{{\beta }_{1}I+{\beta }_2{I}_1+{\beta }_3{I}_2+{\beta }_4{I}_1^T+{\beta }_5{I}_2^T}{N}\right).\end{array}$(2)

Here ${\beta }_i(i=1,2,3,4,5)$ are the transmission contact rates for susceptible individuals with $I,I_1,I_2$, $I_1^T\mathrm{a}\mathrm{n}\mathrm{d}{I}_2^T$ respectively. The natural removal rate $\mu$ represent the death-rate while ${\delta }_1\mathrm{a}\mathrm{n}\mathrm{d}{\delta }_2$ represent AIDS induced death-rate in the classes $I_2$ and $I_2^T$ respectively. The progression rates from $I$ to $I_1$ and from $I$ to $I_2$ are respectively given by ${\alpha }_1$ and ${\alpha }_2$ upon diagnosis. The progression rate from $I_1$ to $I_2$ is given by ${\eta }_1$. The treatment rates from $I_1$ to $I_1^T$ and from $I_2$ to $I_2^T$ are given by ${\gamma }_1$ and ${\gamma }_2$ respectively. These rates can be viewed as the up-take of ART. The proportion of the infected individuals on treatment is bi-directional due to attrition or adherence to ART and decline or improvement of immunological status. The progression between these classes representing improvement or decline in immunological status is given by ${\eta }_2\mathrm{a}\mathrm{n}\mathrm{d}{\eta }_3$ respectively. Furthermore, it is assumed that the transmission of HIV is assumed to be mainly through heterosexual means. In addition, an exit due to death as a consequence of development of AIDS has been included, hence AIDS class is considered redundant and thus left out. The following interventions have been added to the model:

1. The first intervention $0\le {\varphi }_1\le 1$ measures the effectiveness of the efforts made to protect susceptible individuals from infections. It mainly consists of the use of PrEP and or educational campaigns. Thus the infection term given in (2) is modified as follows

   (3) ${\lambda }_c=(1-{\varphi }_1)\lambda ,$(3)

2. The second intervention $0\le {\varphi }_2\le 1$ represents efforts made for counselling and testing. Thus we modify ${\alpha }_i$ to ${\alpha }_i^c={\varphi }_2{\alpha }_i$, for $i=1,2$.

3. The third intervention $0\le {\varphi }_3\le 1$ represents the efficacy of the ART treatment. It mainly consists of up-take of ARV drugs to reduce the risk of dying (Williams, 2014; Williams et al., 2010). Thus we modify ${\gamma }_i$ to ${\gamma }_i^c={\varphi }_3{\gamma }_i$, for $i=1,2$.

Figure 1 describes the five compartments and model variables.

Figure 1. Graphical representation of HIV transmission. Susceptible individuals flow into the undiagnosed (I) via infection rate. Unaware individuals recieve diagnosis and flow into ${\mathit{I}}_{\mathbf{1}}$ and ${\mathit{I}}_{\mathbf{2}}$ based on the CD4 cell counts. The tested individuals are then linked to ART treatment in classes ${\mathit{I}}_{\mathbf{1}}^{\mathit{T}}$ and ${\mathit{I}}_{\mathbf{2}}^{\mathit{T}}$.

3.1.1. Model equations

Given the dynamics described in Figure 1, the following system of non-linear ordinary differential equations, with non-negative initial conditions, are derived.

(4) $\left.\begin{array}{cc}\frac{dS}{dt}& =\mathrm{\Pi }-{\lambda }_{c}S-\mu S,\\ \frac{dI}{dt}& ={\lambda }_{c}S-(({\alpha }_1+{\alpha }_2){\varphi }_2+\mu )I,\\ \frac{d{I}_1}{dt}& ={\varphi }_2{\alpha }_{1}I-({\eta }_1+{\varphi }_3{\gamma }_1+\mu )I_1,\\ \frac{d{I}_2}{dt}& ={\varphi }_2{\alpha }_{2}I+{\eta }_1{I}_1-({\varphi }_3{\gamma }_2+\mu +{\delta }_1)I_2,\\ \frac{d{I}_1^T}{dt}& ={\varphi }_3{\gamma }_1{I}_1+{\eta }_2{I}_2^T-({\eta }_3+\mu )I_1^T,\\ \frac{d{I}_2^T}{dt}& ={\varphi }_3{\gamma }_2{I}_2+{\eta }_3{I}_1^T-({\eta }_2+\mu +{\delta }_2)I_2^T,\end{array}\right\}$(4)

subject to the following initial conditions

(5) $\begin{array}{c}S(0)\ge 0,I(0)\ge 0,I_1(0)\ge 0,I_2(0)\ge 0,I_1^T(0)\ge 0,I_2^T(0)\ge 0.\end{array}$(5)

4. Mathematical analysis

4.1. Well–posedness of the model

In this sub-section, we prove that the system (4) is mathematically well defined and biologically feasible.

Lemma 4.1. Given that the initial conditions of the system (4) are as given in (5), the solutions $S(t),I(t),I_1(t),I_2(t),I_1^T(t)\mathrm{a}\mathrm{n}\mathrm{d}{I}_2^T(t)$ are non-negative for all $t>0$.

Proof. Assume that

$\hat{t}=sup\left\{t>0:S(t)>0,I(t)>0,I_1(t)>0,I_2(t)>0,I_1^T(t)>0,I_2^T(t)>0\right\}\in [0,t].$

Thus $\hat{t}>0$, and it follows directly from the first equation of the system (4) that

$\frac{dS}{dt}=\mathrm{\Pi }-({\lambda }_c+\mu )S.$

Therefore we have

$\begin{array}{c}\frac{d}{dt}\left[S(t)exp\left\{\mu t+{\int }_0^t{\lambda }_c(s)ds\right\}\right]=\mathrm{\Pi }exp\left[\mu t+{\int }_0^t{\lambda }_c(s)ds\right].\end{array}$

Thus

$\begin{array}{c}S(t)exp\left[\mu t+{\int }_0^t{\lambda }_c(s)ds\right]-S(0)={\int }_0^t\mathrm{\Pi }exp\left[\mu t+{\int }_0^t{\lambda }_c(f)df\right]dt,\end{array}$

so that

$\begin{array}{c}S(t)=exp\left[-\left(\mu t+{\int }_0^t{\lambda }_c(s)ds\right)\right]\left\{S(0)+{\int }_0^t\mathrm{\Pi }exp\left[\mu t+{\int }_0^t{\lambda }_c(f)df\right]dt\right\}>0.\end{array}$

Similarly, it can be shown that $I(t)>0,I_1(t)>0,I_2(t)>0,I_1^T(t)>0\mathrm{a}\mathrm{n}\mathrm{d}{I}_2^T(t)>0$ for $t>0$.

Lemma 4.2. Let the feasible region $\mathrm{\Gamma }$ be defined by

$\mathrm{\Gamma }=\left\{(S(t),I(t),I_1(t),I_2(t),I_1^T(t),I_2^T(t))\in {\mathbb{R}}_{+}^6|0\le N\le \frac{\mathrm{\Pi }}{\mu }\right\},$

with the initial conditions as given in (5). The region $\mathrm{\Gamma }$ is positively invariant and attracting with respect to the system (4) for all $t>0$.

Proof. The evolution equation change in the total population is given by

$\frac{dN}{dt}=\mathrm{\Pi }-\mu N.$

So that

(6) $N(t)=\frac{\mathrm{\Pi }}{\mu }+\left(N(0)-\frac{\mathrm{\Pi }}{\mu }\right){\mathrm{e}}^{-\mu t}.$(6)

From (6) we can deduce that $N(t)\le max\left\{N(0),\frac{\mathrm{\Pi }}{\mu }\right\}$. Thus, $\mathrm{\Gamma }$ is positively invariant under the flow induced by the system (4). Therefore in $\mathrm{\Gamma }$, the system (4) is well-posed epidemiologically and hence sufficient to study the system (4) in $\mathrm{\Gamma }$.

4.2. The reproduction number

Supposing that the community has not experienced HIV infection, then it is expected to remain free of the disease. This implies that the entire community is rendered susceptible to the infection. In this respect, there are no infectious individuals. Therefore the system (4) has a disease-free equilibrium given by

${\mathcal{E}}_0=\left[\frac{\mathrm{\Pi }}{\mu },0,0,0,0,0\right].$

The basic reproduction number, ${\mathcal{R}}_0$, is defined as the average number of new infections generated by a single infected individual in a completely susceptible population (Van Den Driessche & Watmough, 2002). We also emphasize that the reproduction number computed should reflect all the routes contributing to new infections. We use the next-generation matrix method as defined in Van Den Driessche and Watmough (2002), such that the matrices of new infections evaluated at the disease-free equilibrium are given by

$\mathbf{F}=(1-{\varphi }_1)\left[\begin{array}{ccccc}{\beta }_1& {\beta }_2& {\beta }_3& {\beta }_4& {\beta }_5\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right],\mathbf{V}=\left[\begin{array}{ccccc}{Q}_1& 0& 0& 0& 0\\ -Q_2& Q_3& 0& 0& 0\\ -Q_4& -{\eta }_1& Q_5& 0& 0\\ 0& -Q_6& 0& Q_7& -{\eta }_2\\ 0& 0& -Q_8& -{\eta }_3& Q_9\end{array}\right].$

Here,

$Q_1=({\alpha }_1+{\alpha }_2){\varphi }_2+\mu ,Q_2={\varphi }_2{\alpha }_1,Q_3={\eta }_1+{\varphi }_3{\gamma }_1+\mu ,Q_4={\varphi }_2{\alpha }_2,$

$Q_5={\varphi }_3{\gamma }_2+\mu +{\delta }_1,Q_6={\varphi }_3{\gamma }_1,Q_7={\eta }_3+\mu ,Q_8={\varphi }_3{\gamma }_2,Q_9={\eta }_2+\mu +{\delta }_2.$

Thus, the basic reproduction number denoted by ${\mathcal{R}}_0$, is given by the spectral radius of $\mathbf{F}{\mathbf{V}}^{-1}$ so that we have

${\mathcal{R}}_0=(1-{\varphi }_1)[R_1+R_2+R_3+R_4+R_5],$

where

$R_1=\frac{{\beta }_1}{Q_1},R_2=\frac{{\beta }_2{Q}_2}{Q_1{Q}_3},R_3={\beta }_3\left[\frac{{\eta }_1{Q}_2+Q_3{Q}_4}{Q_1{Q}_3{Q}_5}\right],$

$R_4={\beta }_4\left[\frac{{\eta }_2{Q}_8\left({\eta }_1{Q}_2+Q_3{Q}_4\right)+Q_2{Q}_5{Q}_6{Q}_9}{Q_1{Q}_3{Q}_5{Q}_7{Q}_9(1-{\mathrm{\Phi }}_1)}\right],$

$R_5={\beta }_5\left[\frac{Q_7{Q}_8\left({\eta }_1{Q}_2+Q_3{Q}_4\right)+{\eta }_3{Q}_2{Q}_5{Q}_6}{Q_1{Q}_3{Q}_5{Q}_7{Q}_9(1-{\mathrm{\Phi }}_1)}\right].$

with ${\mathrm{\Phi }}_1=\frac{{\eta }_2{\eta }_3}{Q_7{Q}_9}$. It is important to note that the reproduction number computed is the sum of sub-population reproduction numbers representing the contributions of individuals in compartments $I,I_1,I_2,I_1^T$ and $I_2^T$. From Theorem 2 in Van Den Driessche and Watmough (2002), the following result is established.

Theorem 4.3. The disease-free equilibrium ${\mathcal{E}}_0$ of the system (4) is locally asymptotically stable whenever ${\mathcal{R}}_0<1$ and unstable when ${\mathcal{R}}_0>1$.

Remark 1. The results in Theorem 4.3 implies that HIV can be eliminated from the community when ${\mathcal{R}}_0<1$ if the initial sizes of the sub-populations of the model are in the basin of attraction of the disease-free equilibrium ${\mathcal{E}}_0$. To ensure that elimination of HIV is independent of the initial sizes of the sub-populations, it is necessary to show that the disease-free equilibrium, ${\mathcal{E}}_0$, is globally asymptotically stable. Thus, we have the following Theorem.

Theorem 4.4. The disease-free equilibrium ${\mathcal{E}}_0$ of the system (4) is globally asymptotically stable if ${\mathcal{R}}_0<1$ and unstable otherwise. The disease-free equilibrium ${\mathcal{E}}_0$ is the only equilibrium when ${\mathcal{R}}_0\le 1$.

Proof. Let

(7) $W={\mathrm{\Psi }}_{1}I+{\mathrm{\Psi }}_2{I}_1+{\mathrm{\Psi }}_3{I}_2+{\mathrm{\Psi }}_4{I}_1^T+{\mathrm{\Psi }}_5{I}_2^T,$(7)

be the Lyapunov function which involves individuals who contribute to infection in the population. The constants ${\mathrm{\Psi }}_1,{\mathrm{\Psi }}_2,{\mathrm{\Psi }}_3,{\mathrm{\Psi }}_4\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{\Psi }}_5$ are all non negative and we ought to find them. The time derivative of the Lyapunov function (7) is given by

(8) $\left.\begin{array}{cc}\frac{dW}{dt}& ={\mathrm{\Psi }}_1\frac{dI}{dt}+{\mathrm{\Psi }}_2\frac{d{I}_1}{dt}+{\mathrm{\Psi }}_3\frac{d{I}_2}{dt}+{\mathrm{\Psi }}_4\frac{d{I}_1^T}{dt}+{\mathrm{\Psi }}_5\frac{d{I}_2^T}{dt},\\ \hspace{0.17em}& \le {\mathrm{\Psi }}_1\left[(1-{\varphi }_1)\left({\beta }_{1}I+{\beta }_2{I}_1+{\beta }_3{I}_2+{\beta }_4{I}_1^T+{\beta }_5{I}_2^T\right)-Q_{1}I\right]+{\mathrm{\Psi }}_2\left[Q_{2}I-Q_3{I}_1\right]\\ \hspace{0.17em}& +{\mathrm{\Psi }}_3\left[Q_{4}I+{\eta }_1{I}_1-Q_5{I}_2\right]+{\mathrm{\Psi }}_4\left[Q_6{I}_1+{\eta }_2{I}_2^T-Q_7{I}_1^T\right]+{\mathrm{\Psi }}_5\left[Q_8{I}_2+{\eta }_3{I}_1^T-Q_9{I}_2^T\right],\\ \hspace{0.17em}& =\left[{\mathrm{\Psi }}_1(1-{\varphi }_1){\beta }_1+{\mathrm{\Psi }}_2{Q}_2+{\mathrm{\Psi }}_3{Q}_4-{\mathrm{\Psi }}_1{Q}_1\right]I+\left[{\mathrm{\Psi }}_1(1-{\varphi }_1){\beta }_2+{\mathrm{\Psi }}_3{\eta }_1+{\mathrm{\Psi }}_4{Q}_6-{\mathrm{\Psi }}_2{Q}_3\right]I_1\\ \hspace{0.17em}& +\left[{\mathrm{\Psi }}_1(1-{\varphi }_1){\beta }_3+{\mathrm{\Psi }}_5{Q}_8-{\mathrm{\Psi }}_3{Q}_5\right]I_2+\left[{\mathrm{\Psi }}_1(1-{\varphi }_1){\beta }_4+{\mathrm{\Psi }}_5{\eta }_3-{\mathrm{\Psi }}_4{Q}_7\right]I_1^T\\ \hspace{0.17em}& +\left[{\mathrm{\Psi }}_1(1-{\varphi }_1){\beta }_5+{\mathrm{\Psi }}_4{\eta }_2-{\mathrm{\Psi }}_5{Q}_9\right]I_2^T.\end{array}\right\}$(8)

We equate the coefficients of the components $I_1,I_2,I_1^T\mathrm{a}\mathrm{n}\mathrm{d}{I}_2^T$ to zero and solve for the coefficients of the Lyapunov function obtaining

$\begin{array}{cc}{\mathrm{\Psi }}_1& =\frac{Q_5{Q}_7{Q}_9(1-{\mathrm{\Phi }}_4)}{(1-{\varphi }_1)},\\ {\mathrm{\Psi }}_2& =\frac{Q_7{Q}_9(1-{\mathrm{\Phi }}_1)[Q_5{\beta }_2+{\eta }_1{\beta }_3]+(Q_5{Q}_6{Q}_9+Q_8{\eta }_1{\eta }_2){\beta }_4+(Q_7{Q}_8{\eta }_1+Q_5{Q}_6{\eta }_3){\beta }_5}{Q_3},\\ {\mathrm{\Psi }}_3& =Q_7{Q}_9{\beta }_3(1-{\mathrm{\Phi }}_1)+Q_8({\eta }_2{\beta }_4+Q_7{\beta }_5),{\mathrm{\Psi }}_4=Q_5({\beta }_5{\eta }_3+{\beta }_4{Q}_9),\\ {\mathrm{\Psi }}_5& =Q_5({\beta }_4{\eta }_2+{\beta }_5{Q}_7).\end{array}$

Substituting the constants ${\mathrm{\Psi }}_1,{\mathrm{\Psi }}_2,{\mathrm{\Psi }}_3,{\mathrm{\Psi }}_4\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{\Psi }}_5$ into (8), we obtain

$\begin{array}{c}\frac{dW}{dt}\le \left[\frac{Q_1{Q}_5{Q}_7{Q}_9(1-{\mathrm{\Phi }}_1)}{(1-{\varphi }_1)}\right][{\mathcal{R}}_0-1].\end{array}$

Thus, $\frac{dW}{dt}\le 0$ when ${\mathcal{R}}_0\le 1$. Furthermore, $\frac{dW}{dt}=0$ if and only if either $I=I_1=I_2=I_1^T=I_2^T=0$ or ${\mathcal{R}}_0=1$ are satisfied. In either case, the largest compact invariant subset of the set

$\mathcal{E}=\left\{(S(t),I(t),I_1(t),I_2(t),I_1^T(t),I_2^T(t))\in \mathrm{\Gamma }:\frac{dW}{dt}=0\right\},$

is the singleton ${\mathcal{E}}_0$. By LaSalle’s Invariance Principle given in La Salle (1976), this implies that ${\mathcal{E}}_0$ is globally stable in $\mathrm{\Gamma }$ provided ${\mathcal{R}}_0\le 1$.

4.3. Existence of the disease-persistent equilibrium

Theorem 4.5. If ${\mathcal{R}}_0>1$, the system (4) has a unique disease-persistent equilibrium given

${\mathcal{E}}_1=\left[S^{\ast },I^{\ast },I_1^{\ast },I_2^{\ast },I_1^{T\ast },I_2^{T\ast }\right].$

Proof. To compute the endemic equilibrium points of the system (4), we equate the right-hand side of the system to zero and solve in terms of force of infection described in (3). Thus, we have

$\begin{array}{rl}& S^{\ast }=\frac{\mathrm{\Pi }}{\mu +{\lambda }_c^{\ast }},I^{\ast }=\frac{\mathrm{\Pi }{\lambda }_c^{\ast }}{Q_1(\mu +{\lambda }_c^{\ast })},I_1^{\ast }=\frac{{\lambda }_c^{\ast }\mathrm{\Pi }{Q}_2}{Q_1{Q}_3({\lambda }_c^{\ast }+\mu )},I_2^{\ast }=\frac{{\lambda }_c^{\ast }\mathrm{\Pi }\left({\eta }_1{Q}_2+Q_3{Q}_4\right)}{Q_1{Q}_3{Q}_5({\lambda }_c^{\ast }+\mu )},\\ & I_1^{T\ast }=\frac{{\lambda }_c^c\mathrm{\Pi }\left({\eta }_2{Q}_8\left({\eta }_1{Q}_2+Q_3{Q}_4\right)+Q_2{Q}_5{Q}_6{Q}_9\right)}{Q_1{Q}_3{Q}_5{Q}_7{Q}_9({\lambda }_c^{\hspace{0.17em}}j+\mu )\left(1-{\mathrm{\Phi }}_1\right)},\\ & I_2^{T\ast }=\frac{{\lambda }_c^{\ast }\mathrm{\Pi }\left(Q_7{Q}_8\left({\eta }_1{Q}_2+Q_3{Q}_4\right)+{\eta }_3{Q}_2{Q}_5{Q}_6\right)}{Q_1{Q}_3{Q}_5{Q}_7{Q}_9({\lambda }_c^{\ast }+\mu )\left(1-{\mathrm{\Phi }}_1\right)}.\end{array}$

Substituting the expressions for the state variables $I^{\ast },I_1^{\ast },I_2^{\ast },I_1^{T\ast }\mathrm{a}\mathrm{n}\mathrm{d}{I}_2^{T\ast }$ into (3) and simplifying, we obtain the following polynomial

(9) $B_1{\lambda }_c^{\ast }+B_0=0,$(9)

where

$B_1=Q_1{Q}_3{Q}_5{Q}_7{Q}_9(1-{\mathrm{\Phi }}_1),B_0=\mu (1-{\mathcal{R}}_0).$

We note that $B_0>0$ if ${\mathcal{R}}_0<1$ and expression (9) does not have a positive solution. However, $B_0<0$ if ${\mathcal{R}}_0>1$ and expression (9) does have a positive solution implying that there exists a unique disease-persistent equilibrium if and only if ${\mathcal{R}}_0>1$.

4.4. Global stability of the disease-persistent equilibrium

Theorem 4.6. The disease-persistent equilibrium ${\mathcal{E}}_1$ is globally asymptotically stable whenever ${\mathcal{R}}_0$ is greater than unity

Proof. We know that $N=\frac{\mathrm{\Pi }}{\mu }\mathrm{a}\mathrm{s}t\to \mathrm{\infty }$. Therefore, using the relation $S=\frac{\mathrm{\Pi }}{\mu }-I-I_1-I_2-I_1^T-I_2^T$ and substituting in (4) gives the following limiting system

(10) $\left.\begin{array}{cc}\frac{dI}{dt}& ={\lambda }_c\left[\frac{\mathrm{\Pi }}{\mu }-I-I_1-I_2-I_1^T-I_2^T\right]-Q_{1}I,\\ \frac{d{I}_1}{dt}& =Q_{2}I-Q_3{I}_1,\\ \frac{d{I}_2}{dt}& =Q_{4}I+{\eta }_1{I}_1-Q_5{I}_2,\\ \frac{d{I}_1^T}{dt}& =Q_6{I}_1+{\eta }_2{I}_2^T-Q_7{I}_1^T,\\ \frac{d{I}_2^T}{dt}& =Q_8{I}_2+{\eta }_3{I}_1^T-Q_9{I}_2^T.\end{array}\right\}$(10)

Using the Dulacs multiplier ${\mathrm{\Phi }}_2(I,I_2)=\frac{1}{I{I}_2}$, (see Mukandavire, Gumel, Garira, and Tchuenche (2009)), it follows that

$\begin{array}{rl}& \mathrm{\partial }{\mathrm{\Phi }}_2=\frac{\mathrm{\partial }}{\mathrm{\partial }I}\left[\left(\frac{\mu (1-{\varphi }_1)({\beta }_{1}I+{\beta }_2{I}_1+{\beta }_3{I}_2+{\beta }_4{I}_1^T+{\beta }_5{I}_2^T)}{I{I}_1{I}_2{I}_1^T{I}_2^T\mathrm{\Pi }}\right)\right.\\ & \left(\frac{\mathrm{\Pi }}{\mu }-I-I_1-I_2-I_1^T-I_2^T\right)-\frac{Q_1}{I_1{I}_2{I}_1^T{I}_2^T}\\ & +\frac{\mathrm{\partial }}{\mathrm{\partial }{I}_2}\left[\frac{Q_4}{I_1{I}_2{I}_1^T{I}_2^T}+\frac{{\eta }_1}{I{I}_2{I}_1^T{I}_2^T}-\frac{Q_5}{I{I}_1{I}_1^T{I}_2^T}\right],\\ & =-\left[\frac{\mathrm{\Pi }{Q}_4{I}^2+{\eta }_1\mathrm{\Pi }I{I}_1+{\mathrm{\Phi }}_3}{\mathrm{\Pi }{I}^2{I}_1{I}_2^2{I}_1^T{1}_2^T}\right]<0,\end{array}$

where

${\mathrm{\Phi }}_3=I^2{I}_2{\beta }_1\mu (1-{\varphi }_1)+I_2(1-{\varphi }_1)(\mathrm{\Pi }-\mu (I_1+I_2+I_1^T+I_2^T))[{\beta }_2{I}_1+{\beta }_3{I}_2+{\beta }_4{I}_1^T+{\beta }_5{I}_2^T].$

Thus, by Dulacs criterion, there are no periodic orbits in $\mathrm{\Gamma }$. Since $\mathrm{\Gamma }$ is positively invariant, and the disease-persistent equilibrium exists whenever ${\mathcal{R}}_0>1$, then it follows from the Poincare-Bendixson Theorem (Perko, 2013) that all solutions of the limiting system originating in $\mathrm{\Gamma }$ remain in $\mathrm{\Gamma }$ for all $t$. Furthermore, the absence of periodic orbits in $\mathrm{\Gamma }$ implies that the unique disease-persistent equilibrium of the system (4) is globally asymptotically stable whenever ${\mathcal{R}}_0>1$.

5. Numerical simulations

5.1. Model validation

To estimate the parameter values and validate the model, we fit the system (4) to the data for the estimated new HIV infection in Kenya. The data were obtained from our world in data website (Roser & Ortiz-Ospina, 2017). These data were reported on an yearly estimate from 1990 to 2017. The least-squares curve fit is used to fit the model to data. Many parameters are known to lie within limits. Only a few parameters are known exactly and it is thus important to estimate the others. The estimation process attempts to find the best concordance between computed and observed data. It can be carried out by trial and error or by the use of software programs designed to find parameters that give the best fit. Here, the fitting process involves the use of the least squares-curve fitting method. A Matlab code is used where unknown parameter values are given a lower and upper bound from which the set of parameter values that produce the best fit are obtained. The intervals used and a few parameters obtained from literature are given in Table 1. The following initial values were used.

Table 1. Parameter values and ranges per year obtained from data fitting

Par	Min	Max	Value	Source
$\mathrm{\Pi }$	0.0308N	0.0409N	0.0345N	WB (2017)
${\beta }_1$	0.0	2.0	1.9652	Okosun et al. (2013); Tripathi, Naresh, and Sharma (2007)
${\beta }_2$	0.0	2.0	2.000	Okosun et al. (2013); Tripathi et al. (2007)
${\beta }_3$	0.0	2.0	1.4153	Okosun et al. (2013); Tripathi et al. (2007)
${\beta }_4$	0.0	1.5	0.0160	Okosun et al. (2013)
${\beta }_5$	0.0	1.5	0.0376	Okosun et al. (2013)
${\alpha }_1$	0.1	1.5	1.4143	Okosun et al. (2013)
${\alpha }_2$	0.1	1.5	1.4363	Okosun et al. (2013)
${\eta }_1$	0.1	1.5	1.4732	Okosun et al. (2013)
${\gamma }_1$	0.0	1.2	0.5186	Assumed
${\gamma }_2$	0.0	1.5	1.4360	Assumed
${\eta }_2$	0.01	2.0	1.9960	Assumed
${\eta }_3$	0.0	1.0	0.3785	Assumed
$\mu$	0.014	0.02	0.0158	WB (2017)
${\delta }_1$	0.083	0.125	0.0845	Okosun et al. (2013)
${\delta }_2$	0.083	0.125	0.0835	Okosun et al. (2013)
ϕ	0.001	0.9	0.0511	Assumed
ϕ	0.001	0.9	0.7916	Assumed
ϕ	0.01	0.9	0.7990	Assumed

(i) $N=23,000,000$, population that were found in 1990 (KNBS, 2017; KNBS2, 2017).

(ii) $S(0)=N-I(0)-I_1(0)-I_2(0)-I_1^T(0)-I_2^T(0)=22,440,000$.

(iii) $I(0)=330,000$ is assumed number of asymptomatic HIV individuals (undiagnosed).

(iv) $I_1(0)=200,000$, reported number of new HIV infections at the end of 1990.

(v) $I_2(0)=30,000$ (assumed), $I_1^T(0)=0,I_2^T(0)$ (all assumed).

We observe in Figure 2(a) that the model fits well with the data. It is important to observe that the cases of new HIV infection peaked in the year 1994. The results show that there was a rise in HIV infection between 1990 and 1995, followed by a significant slow down in the occurrences ofnew cases of infection. The data show an epidemic stabilizing at about 80,000 cases of HIV infection. The model also shows a steady state solution close to this value. The projected cases of new infection in Kenya till 2030 are given in Figure 2(b). We observe that the cases of new infection remain around the steady state registering about 80,000 cases of infection per year towards 2030. This is supported by the data that seem to fluctuate around the steady state. Our results are suggestive of the fact that interventions aimed at averting new infections will be most appropriate in achieving Kenya’s vision 2030 of HIV eradication.

Figure 2. (a) shows the trend of new HIV infections while (b) shows the projected new HIV infections in Kenya by 2030.

5.2. Sensitivity analysis

In this sub-section, we investigate the effect of parameter values on the basic reproduction number, ${\mathcal{R}}_0$. The uncertainties involved in data collection and the estimation parameter values lead us to evaluating the robustness of the model predictions with the parameter values as well as the effect of varying a parameter value on ${\mathcal{R}}_0$. To examine the sensitivity of ${\mathcal{R}}_0$ to variation in parameters, we used the normalised forward sensitivity index also known as elasticity (Helton, Iman, & Brown, 1985). If we let ${\mathcal{R}}_0$ to be a differentiable function with respect to each of the parameters, then the normalised forward sensitivity index of ${\mathcal{R}}_0$ with respect to the parameter $\theta$ is given by

(11) $S_{\theta }^{{\mathcal{R}}_0}=\frac{\mathrm{\partial }{\mathcal{R}}_0}{\mathrm{\partial }\theta }\times \frac{\theta }{{\mathcal{R}}_0}.$(11)

It is important to note that the quotient $\frac{\theta }{{\mathcal{R}}_0}$ is introduced to normalise the coefficient by removing the effect of units. Furthermore, this computation is based on the assumption that higher order partial derivatives are negligible and that there is no existence of correlation between input parameters (Hamby, 1994; Omondi, Orwa, & Nyabadza, 2017). Using the elasticity formula given in (11) with the parameter values in Table 1, we now obtain the numerical values for the elasticities and display the results in Table 2.

Table 2. Parameter values and ranges used in numerical simulations as well as determination of sensitivity indices of ${\mathcal{R}}_0$. The parameter rates are given per annum

Par	Index	Par	Index	Par	Index	Par	Index
${\beta }_1$	+0.4331	${\alpha }_1$	−0.0360	${\eta }_2$	−0.0787	ϕ	−0.0539
${\beta }_2$	+0.3864	${\alpha }_2$	−0.1117	${\eta }_3$	+0.0772	ϕ	−0.0817
${\beta }_3$	+0.2433	${\eta }_1$	−0.0606	$\mu$	−0.1933	ϕ	−0.1542
${\beta }_4$	+0.2319	${\gamma }_1$	−0.0486	${\delta }_1$	−0.0180
${\beta }_5$	+0.1053	${\gamma }_2$	−0.1426	${\delta }_2$	−0.1547

The sensitivity index values of ${\mathcal{R}}_0$ in Table 2 show that an increase/decrease of 1% in any of the parameter values results in a percentage increase/decrease in ${\mathcal{R}}_0$. Model system (4) is most sensitive to ${\beta }_1,{\beta }_2,{\beta }_3$ and ${\beta }_4$. These are the parameters with the greatest potential to increase the HIV infection. Of the three controls added, ${\varphi }_3$ has the greatest effect on the basic reproduction number. This implies that ART efficacy is essential in reducing the HIV infection in the population.

5.3. Numerical results

In this section, we carry out numerical simulation by fourth order Runge-Kutta method in Matlab 2017a to see the dynamical behaviour of the system (4) in the presence of controls using parameters: $\mathrm{\Pi }=0.02\mathrm{N}$, where $N$ is the sum of the initial conditions of the state variables presented in (12), $\begin{array}{c}{\beta }_1=0.981,{\beta }_2=0.725,{\beta }_3=0.6670,{\beta }_4=0.4182,{\beta }_5=0.351,{\alpha }_1=0.8749,{\alpha }_2=0.2794,{\eta }_1=0.900,\\ {\gamma }_1=0.6000,{\gamma }_2=0.8758,{\eta }_2=0.5998,{\eta }_3=0.0109,\mu =0.02,{\delta }_1=0.0845,{\delta }_2=\mathrm{0.0835.}\end{array}$ The computer simulations are performed using the following initial conditions.

(12) $\begin{array}{cc}S(0)& =800,000,I(0)=50,000,I_1(0)=10,000,I_2(0)=25,000,I_1^T(0)=5000,\\ I_2^T(0)& =2000.\end{array}$(12)

5.3.1. Control of new infection combined with ART treatment

In this strategy, we set ${\varphi }_2={\varphi }_3=1$ and vary the parameter related to protection of susceptible from infection, namely ${\varphi }_1$, between 0.2 and 0.9. Figure 3 shows that the increase of the individual protection level, reduces the disease prevalence. From this result, it is clear that this reduction is significant if the level of protection is kept high over a long period. Thus, continuous education campaigns of local populations, on how to protect themselves individually using the available protection measures such as PrEP, are important in the fight against the spread of HIV.

Figure 3. Simulation result showing the effect of efforts made to protect susceptible individuals from infections on the disease prevalence.

5.3.2. Counselling and testing combined with ART treatment

Since the enrolment of an individual into ART treatment programme depends on the diagnosis of HIV, we explore the benefits of enhancing the counselling and testing of individuals in the presence of ART treatment. We fix the efficacy of ART treatment (${\varphi }_3)$ to 1 and the control of new infections (${\varphi }_1)$ to zero. Simulation results in Figure 4 show that an increase in counselling and testing of individuals, ${\varphi }_2$, has a significant impact on the decrease of the disease prevalence. However, this decrease settles asymptotically at prevalence slightly above 1% implying that further increase of HIV counselling and testing without controlling new infections may not halt the disease spread.

Figure 4. Simulation result showing the testing and counselling on the disease prevalence.

5.3.3. Control of new infections, counselling and testing combined with ART treatment

The results in Figure 5 show that combining all three control strategies yield a positive impact in the control of HIV. The results indicate that if the controls are maintained at over 90%, then it will be very easy to halt the further spread of HIV.

Figure 5. Simulation results showing the advantage of combining personal protection, testing and counselling and ART treatment on the disease prevalence.

6. Conclusion

In this paper, an HIV transmission model is derived and analyzed. The basic reproduction number ${\mathcal{R}}_0$ is computed and stability analysis of the model equilibria carried out. By constructing a suitable Lyapunov function, the analysis shows that the disease-free equilibrium of the model is globally asymptotically stable whenever ${\mathcal{R}}_0$ is less than unity. The health implications of this observation are that keeping the reproduction number below unity may be necessary in the control of the spread of HIV epidemic. However, when ${\mathcal{R}}_0>1$, the disease-free equilibrium loses its stability and there exists endemic equilibrium which is globally asymptotically stable.

The least squares curve fit routine in Matlab with optimisation has been used to fit the model to data on new HIV infections with the objective of using the model parameters that give the best fit to predict the cases of new infections by 2030. Using parameter values obtained from the curve fitting, we have calculated the sensitivity indices of the basic reproduction number, ${\mathcal{R}}_0$, to the parameters in the model using normalised sensitivity index method. The results of sensitivity analysis showed that the model system is most sensitive to infection contact rates, the testing and counselling rates as well as treatment rates.

To assess the impact of combination of different controls, simulations have been conducted. According to the numerical results, we conclude that both high individual protection and counselling and testing combined with ART treatment, can effectively reduce the transmission and the proliferation of HIV in the population. Thus, more resources should focus on the education of the local populations so as to involve more people in the fight against HIV.

The model is not without limitations. The model did not take into account immigration of infected individuals into the system as well as other transmission routes such as the one in the people who inject drugs. Incorporating these processes will undoubtedly facilitate in the understanding of HIV transmission and control in Kenya. The models greatest value is not as a predictive tool but rather a tool for organizing the thinking around HIV infections and the generation of useful insights into the dynamics of HIV transmission. The dynamics are much more complex than described here. We argue that the presented model provides the required insights to understanding HIV spread. The results are relevant in designing time-dependent intervention strategies instead of the constant controls using optimal control theory to effectively examine the cost of each control aimed at combating spread of HIV. Thus, the next stage of this analysis will extend the dynamical model of the epidemic in order to further test the ideas and suggestions discussed here and at the same time project future time trends in the epidemic and to carry out detailed cost benefit analyses.

Acknowledgements

The authors acknowledge, with thanks, the support of Strathmore University, Institute of Mathematical Sciences for the production of this manuscript. The authors are also grateful to the reviewer for the insightful comments.

Additional information

Funding

The authors received no direct funding for this research.

Notes on contributors

E. O. Omondi

E. O. Omondi is a doctoral fellow at Strathmore University in Kenya. His research interest is in the field of mathematical biology. He is also a teaching fellow at the Institute of Mathematical Sciences of Strathmore University which delivers both undergraduate and postgraduate education and training.

R. W. Mbogo

R. W. Mbogo is the Director in charge of Graduate Research and Training at the Institute of Mathematical Sciences of Strathmore University, Kenya. Her area of interest includes Mathematical Epidemiology and specifically the application of mathematics to disease dynamics.

L. S. Luboobi

L. S. Luboobi is a visiting professor in the Institute of Mathematical Sciences of Strathmore University, Kenya. He has published extensively on models of tropical infections. Our joint research contribution is towards modelling the effectiveness of a combination of HIV control strategies.