Monotonicity properties arising in a simple model of Wolbachia invasion for wild mosquito populations

1.   Introduction

The first official case of HIV in Botswana was identified in 1985. The nation's response to the epidemic dates from the late 1980s with the Short Term Plan (STP) for the years 1987-1989, followed by the Medium Term Plan Ⅰ (MTP Ⅰ) for the years 1991–1996, which outlined a medical and health system response to HIV/AIDS ^[1]. The Medium Term Plan Ⅱ (MTP Ⅱ) for 1997–2002 developed the institutional structures necessary to organise and manage the national response. The national policy was revised in 1998 following new scientific, epidemiological and medical knowledge about the distribution and transmission of HIV and other sexually transmitted diseases (STDs) and proven effective interventions in prevention and care. The evolution of the national response in Botswana's continued with developing the National Strategic Framework (NSF) for HIV/AIDS 2003–2009. The review of MTP Ⅱ informed the NSF. The second national strategic framework (NSF Ⅱ) established in 2010 outlined the national priorities for the national response from 2010–2016. NSF Ⅱ focused on a costed multi-sectoral and rolling operational plan that guided the implementation of the national response. Due to increased global research on HIV/AIDS dynamics, new knowledge and information on the field rapidly kept changing. Therefore, a review of the HIV/AIDS policy was needed as new findings emerged. As a result, the Botswana national policy on HIV and AIDS was revised in 2012.

Botswana is now on its third national strategic framework (NSF Ⅲ) for 2019–2023. The NSF Ⅲ targets ending AIDS by 2030. This strategic framework has been informed by relevant global, regional and national policies and strategic frameworks. NSF Ⅲ goals include; prioritising an HIV combination prevention strategy, targeting and exceeding the global target (90-90-90) of HIV testing and treatment, reduce HIV related stigma, discrimination and gender vulnerability while promoting and protecting the human rights of people living with HIV.

According to the World Health Organization (WHO), health policy refers to decisions, plans, and actions undertaken to achieve specific health care goals within a society. An explicit health policy can achieve several objectives; it defines a vision for the future, outlines priorities and the expected roles of different population groups, and builds consensus and informs people ^[2]. Health policy categories may be global, national, local, public and workplace-based. National health policies, strategies and plans play an essential role in defining a country's vision, policy directions, and strategies for ensuring her population's health. For example, an HIV/AIDS policy is a health policy that provides the basic framework for reducing the spread of HIV/AIDS and managing its impacts. Countries have adopted targets for achieving universal access to HIV prevention, treatment, care and support by producing national HIV/AIDS policies.

Mathematical models have a vital role in assisting policymakers in developing HIV/AIDS prevention and treatment strategies. They help in bridging the gap between theory and policy. Some of the several ways models can assist HIV/AIDS policymakers are; comparing the impact and cost-effectiveness of different interventions, generalising the results of randomised controlled trials to different settings, identifying threats to programme success and opportunities for maximising intervention efficiency, evaluating and motoring the extent to which observed trends in HIV prevalence are attributable to HIV/AIDS programme success ^[3]. Mathematical models have also been used to obtain insights into other HIV/AIDS policy research. They include; evaluation of treatment as a prevention strategy ^[4,5], the scale-up of microbicides ^[6], the impact and cost-effectiveness of HIV Pre-exposure Prophylaxis (PrEP) targeting key populations ^[7,8], the efficacy of the Post-Exposure Prophylaxis (PEP) ^[9], projecting HIV diagnoses among children and adolescents ^[10,11], HIV vaccines in future with the presence of current interventions ^[12], the impact of voluntary medical male circumcision (VMMC) ^[13] and other recent recent developments in modelling techniques for HIV/AIDS models ^[14,15,16]. Many mathematical models of HIV/AIDS, with different structures, complexity and parameter choices have contributed to HIV/AIDS research for many years. However, none of the models including the recent work, has tracked how the change in HIV/AIDS policy over the years have impacted or influenced the trajectory of the epidemic. We thus present a unique model that traces the impact of policy changes over time.

In this study, we track the changes in policy for HIV/AIDS testing and treatment in Botswana. In the first phase of the national HIV/AIDS policy, voluntary counselling and testing are encouraged. However, the policy is silent on treatment but rather focusing on behavioural change. Botswana was the first African country to establish a national HIV/AIDS treatment programme, called "Masa", the Setswana word for "a new dawn", in 2002, leading to the second phase of national HIV/AIDS policies. Since the start of the Masa programme, the national guidelines changed to take into account the improved understanding of the biology of HIV, reduce adverse side effects associated with treatment, and accommodate the availability of improved drugs ^[17]. In 2008, the eligibility criteria threshold for Antiretroviral therapy (ART) changed from a CD4 cell count of 200 cells per cubic millimetre to 250 cells per cubic millimetre. In 2013, it changed again to 350 cells per cubic millimetre. In June 2016, Botswana launched the "Treat All" strategy, promoting universal health coverage and ensuring that all individuals who test positive for HIV get treatment immediately regardless of their viral load or CD4 count. Mathematical models have been developed and fitted to HIV/AIDS data in different countries, see for instance ^[18,19]. All these articles fit one model to data spanning over a long period. However, none of the models account for the change in dynamics due to the policy changes. We aim to develop deterministic models to integrate HIV/AIDS prevalence data amongst adults, 15 years and above, for Botswana and the change in policies. The objective of this study is to examine the impacts of change in HIV treatment policy over time.

This paper is arranged as follows: The models are developed in the next section and analysed in Section 3. Numerical simulations are presented in Section 4 and the paper is concluded in Section 5.

2.   Model development

The HIV/AIDS treatment policy for Botswana can thus be divided into three epochs with respect to the treatment policy. It is important to note that additional epochs can be considered depending on the nature of the response being considered. The first epoch captures the dynamics before introducing the universal ART programme (1990–2001). The second epoch captures the period from which the universal ART programme was rolled out (2002–2015) with criterion eligibility of CD4 count to enrol. Finally, the third epoch captures the era (2017–to date) where eligibility to enrol on ART is a positive test. Comparing the three models, it is notable that the first model changed by introducing treatment in 2002. As a result, two compartments were added in the second epoch to cater for the treatment policy. It is important to note that this modelling framework is similar to that used in ^[20,21], where individuals were classified according to their immunological staging when they test positive. In 2016, when the treat all policy was launched, the model changed by dropping out of the waiting compartment. We note that the second epoch model is the biggest, and the other two models can be deduced from it by setting some parameters to zero. Given that the first and third epochs can be obtained from the second epoch from a modelling perspective, we begin by looking at the model for the second epoch, which we dub here the "grand model".

2.1.   The grand model

We propose a model that captures the dynamics of HIV/AIDS with an ART programme that considers CD4 cell count as an eligibility criterion for treatment. The population of sexually active adults, $N(t),$ is divided into five classes: $S(t),$ for the adult population at risk of being infected by HIV and $I(t)$, for individuals who have contracted HIV and are unaware of their status. Due to unavailable data on the CD4 cell count of the infected individuals, screened individuals are divided into two groups, $I_W (t)$ and $I_T (t)$, where $I_W(t)$ represents those who tested positive for HIV but with a CD4 cell count above the threshold, therefore not legible to enrol in ART and $I_T(t)$ represents those individuals who tested positive for HIV and are legible to enrol into ART. $A(t)$, represent HIV infected individuals who will progress to the AIDS stage. Thus,

We assume that susceptible individuals are recruited, as they become sexually active, at a rate $\Lambda$. Susceptible individuals are infected at a rate $\lambda,$ whose full description is given with system (2.1). Infected individuals in class $I$ can either develop AIDS or get into treatment upon testing if their CD 4 cell count is below a set threshold or they wait till their CD4 count level decreases to the required threshold. A proportion $\epsilon$ of those tested at a rate $\theta$ will have a CD4 count level above the threshold while the remaining proportion will be treated upon testing. Individuals in $I_W$ will progress to the class $I_T$ after a decrease in their CD4 cell count meets the required CD4 cell count level at a rate $\gamma.$ All individuals with HIV, i.e those in $I$, $I_W (t)$ and $I_T (t)$ will progress to the AIDS class $A(t)$, at rates $\rho_1$, $\rho_2$ and $\rho_3$ respectively. Individuals in the AIDS class are assumed to transmit infection and they also die due to the disease at a rate $\delta.$ Individuals in each class are assumed to die naturally at a rate $\mu.$ The forces of infection are defined by ${\bf \lambda} = (\lambda_1, \lambda_2, \lambda_3)$ for the three epochs. The explicit expressions are given in the description of each epoch. The dynamics of the second epoch are presented in Figure 1.

A summary of the model parameters is given in Table 1.

The following system of ordinary differential equations thus governs the second epoch model,

where $X_2 = (S, I, I_W, I_T, A)^T$ and the force of infection $\lambda_2 (t)$ is given by

Here, $\beta$ reflects not only the infectiousness per contact and the rate of sexual contact per unit time but also incorporates the impact of condom use as a primary prevention tool. Assuming that the level of protection by condoms is given by $p$, where $0 \leq p \leq 1$. If $p = 0$ then condoms do not offer any protection while $p = 1$ implies perfect protection. Thus, $(1-p)$ measures condom failure in preventing HIV transmission, see ^[22] for a detailed explanation. The number of sexual contacts in which condoms fail is represented by $c = n (1-p)$ where $n$ is the number of sexual partners per unit time. We thus define $\beta = \hat{\beta}n(1-p)$, where $\hat{\beta}$ is the effective infectiousness per contact. This force of infection is reasonable for HIV since as the infection spreads, it is most likely that the remaining susceptible individuals become more cautious about their contacts and potential exposure to HIV, see also ^[23]. The force of infection will begin to decline as the number of HIV/AIDS cases and deaths due to the disease cases increases due to self-imposed fear ^[24]. Therefore $m$ represents the fear of infection that will result in behavioural change, driven by the prevalence of the infection in the population. The parameter can be viewed as the measure in the reduction of the infection rate induced by the fear generated per infected case.

2.2.   The first epoch

Given that there was no treatment in the period 1990 to 2001. We set the parameters related to treatment to zero in the grand model to obtain the first epoch model. The parameters related to treatment in the grand model are $\theta, \gamma, \rho_{2}, \rho_{3}$ and setting them to zero results an $SIA$ model for the first epoch.

Our model is thus governed by the following system of ordinary differential equations represented as follows,

where $X_1 = (S, I, A)^T$ and

Note that the mathematical analysis of system (2.2) is presented in ^[25].

2.3.   The third epoch

In June 2016, the Botswana government launched the "Treat All" strategy, promoting universal health coverage and ensuring that all individuals who test positive for HIV enrol to treatment immediately regardless of their viral load or CD4 cell count. An assumption is made that all individuals in $I_W$ will move to $I_T$ since they already know their status and have no reason to wait. Setting the parameters, $\epsilon$, $\rho_{2}$ and $\gamma$, related to the non-legible to treatment compartment $I_W$ to zero in the grand model, we obtain the third epoch model. As a result, we have a four-compartment deterministic model of susceptible individuals $(S)$, infected individuals who are unaware of their HIV infection status $(I)$, infected individuals who have tested and immediately enrolled into treatment $(I_{T})$ and infected who develop AIDS $(A)$.

The third epoch model is thus governed by the following system of ordinary differential equations,

where $X_3 = (S, I, T_T, A)^T.$ The force of infection $\lambda_3 (t)$ is given by

The parameters are as described in Table 1.

3.   Mathematical analysis of the models

Given that the first and the third epochs can be obtained from the grand model, we consider the mathematical analysis of system (2.1). The mathematical analysis of the other two models can be drawn from the big model by setting relevant parameters to zero. We present the mathematical analysis of model (2.1) subject to the initial conditions:

It is important to note that the analysis results of the grand model captures those of the models of first and third epochs given that they are "sub-models" of the grand model.

3.1.   Model properties

Lemma 3.1. The model (2.1) with the initial conditions (3.1), has non-negative solutions and the solution of the system will remain positive for all $t > 0$.

Proof. From the first equation of the model system (2.1), let $\lambda = \lambda_2$, we have

We thus have

From the second equation of the model system (2.1) we obtain

whose solution is given by

In a similar manner we can demonstrate that $I_{W}(t) > 0$, $I_{T}(t) > 0$ and $A(t) > 0$ for all $t > 0$.

Lemma 3.2. The model (2.1) solutions are uniformly bounded in the set

Proof. All parameters and initial conditions in the system (2.1) are assumed to be positive.

The equations of system (2.1) gives

The solution of the differential inequality is given by

If $N_0 > \frac{\Lambda}{\mu}$, the solutions of $N(t)$ decreases to $\frac{\Lambda}{\mu}$ and if $N_0 < \frac{\Lambda}{\mu}$ then $N(t) \rightarrow \frac{\Lambda}{\mu}$ as $t\rightarrow \infty$. The region $\Omega$ is thus positively invariant and solutions are bounded. This means that every solution of (2.1) with initial conditions in $\Omega$ remain in $\Omega$ for all $t\geq 0$. The model is thus epidemiologically and mathematically well-posed in the region

3.2.   Model equilibria

The steady states of the model are obtained by setting the equation of system (2.1) to zero such that

where $\sigma_1 = \theta + \rho_1 + \mu$, $\sigma_2 = \gamma +\rho_2 + \mu$, $\sigma_3 = \rho_3 + \mu$ and $\sigma_4 = \delta + \mu.$

The steady state variables expressed in terms of $I^*$ are given by,

where

We also have

where $\phi_{1} = 1+ \psi_1 + \psi_2+ \psi_3$ and $\phi_2 = \beta(1+\eta_{1} \psi_3+ \eta_{2} \psi_1 + \eta_{3} \psi_2).$

Substituting (3.4) into the second equation (3.2) and solving for $I^*$ we obtain

If $I^* = 0$ then we have $I_W^* = I_T^* = A_0^* = 0$ and $S^* = \frac{\Lambda}{\mu}$. This corresponds to the disease-free state. The model thus has a disease-free equilibrium (DFE) given by

3.2.1.   The reproduction number

The reproduction number, $R_0$, gives a threshold condition to the stability of the disease free disease-free equilibrium ^[23]. We compute the reproduction number of the grand model system 2.1 using the next generation method. Following ^[26], and adopting the matrix notation, the matrices for new infection terms, $F$, and the transfer terms, $V$, at the DFE are given by

$FV^{-1}$ is the next generation matrix. The spectral radius, of the matrix $FV^{-1}$ is the reproduction number of the model. Thus

where

$R_0$ is the sum of four terms representing the contribution of infection by classes, $I, I_W, I_T$ and $A$. We note that the $R_0$ for the first epoch model can be found by setting the parameters, $\theta, \; \gamma, \; \rho_{2}, \; \rho_{3},$ in Eq (3.6) to zero. The same applies to the reproduction number of the third epoch, which is obtained by setting the parameters $\epsilon, \gamma$ and $\rho_2$ in Eq (3.6) to zero. It is defined as the expected number of secondary cases produced by a typical infected individual during its entire period of infectiousness in a completely susceptible population and is mathematically defined as the dominant eigenvalue of a positive linear operator.

Following Theorem 2 in ^[26], we state the stability of the disease free equilibrium, $\mathcal{D}_{0},$ as follows:

Theorem 3.2.1. The DFE, $\mathcal{D}_{0},$ is locally asymptotically stable if $R_0 < R_0^c < 1$ and unstable otherwise.

Note that the expression for $R_0^c$ is defined in the next subsection.

3.2.2.   Existence of the endemic equilibrium

The total population at the steady state is given by

Substituting (3.7) into Eq (3.5) yields

Substituting Eqs (3.4), (3.7) and (3.8) into the first equation of the system (3.2) we obtain the quadratic polynomial in $I^*$ given by

with

Here, $R_{0}$ is the model reproduction number given by the simplified expression

where $R_1, \; R_2, \; R_3$ and $R_4$ are define in (3.6).

The existence of the endemic equilibrium (EE) is subjected to the roots of the quadratic equation (3.9) being positive, we thus have

The coefficient of $I^*$ in Eq (3.9) can either be positive or negative. Also, $\nu_0$ can be either positive or negative depending on whether $R_0$ is less or greater than unity. The results are presented in Table 2.

We note from Table 2 that for the existence and uniqueness of the endemic equilibrium, $\nu_1$ and $\nu_0$ must be negative as $\nu_2$ is always positive for $R_0 > 1$. Thus, we have the following result.

Theorem 3.2.2. The model has:

$i.$ a unique endemic equilibrium if $\nu_0 < 0, R_0 > 1,$

$ii.$ a unique endemic equilibrium if $\nu_1^2 - 4\nu_2\nu_0 = 0,$

$iii.$ two endemic equilibria if $\nu_1 < 0$ and $\nu_0 > 0, R_0 < 1$,

$iv.$ no endemic equilibrium otherwise.

From case $(i)$ of Theorem 3.2.2, there exist a unique endemic equilibrium whenever $R_0 > 1$. Furthermore, case $(iii)$ indicates the possibility of backward bifurcation. Backward bifurcation in epidemic models implies a co-existence of a stable DFE with a stable EE when the associated reproduction number is less than unity, see also ^[27]. The epidemiological implications of this phenomenon is that the requirement $R_0 < 1$ is not sufficient for the disease elimination effort rather the effort is described by the value of the critical parameter at the turning point ^[28]. We set the discriminant of Eq (3.9) to zero and solve for the critical value of $R_0$, denoted $R_0^C$ so that

The backward bifurcation is depicted in Figure 2 for the parameter values given in the caption. The backward bifurcation shows a locally stable DFE and an unstable and stable EE for $R_0 < 1$. To effectively control the disease, a globally stable DFE is desired for $R_0 < 1$. In this case, the DFE is globally stable below the threshold $R_0^c \approx 0.6$ making it challenging to eliminate the disease because it is far from unity for the chosen parameter values. Backward bifurcation has been observed in models for disease dynamics such as those for behavioural responses to perceived risk, treatment, in-host dynamics, the imperfect vaccine, public health education campaigns and reinfection in the transmission dynamics ^[27].

3.2.3.   The global stability of the endemic equilibrium

Theorem 3.2.3.

The endemic equilibrium point $\mathcal{D}_{1}$ of the model system equations is globally asymptotically stable if $R_0 > 1.$

The proof of the theorem in given in Appendix A.

4.   Numerical simulations

4.1.   Parameter estimation

This section estimates the model's unknown parameters by calibration or curve fitting. Our data is extracted from the UNAIDS website (https://aidsinfo.unaids.org/). We consider data on the number of adults (aged 15 year and above, living with HIV in Botswana from 1990 to 2019. We first attempt to fit each of the three models, described by systems (2.1), (2.2) and (2.3), to the prevalence data of Botswana for the entire modelling period. A Matlab built-in function ode45 is used to solve the systems and the fminsearch algorithm is used for the minimisation routine. The data points are to be compared with simulation results by minimizing the sum of square difference of the models' prevalence $P(t),$ (i.e the number of all HIV and AIDS-infected people over the modelling time) and the HIV/AIDS prevalence for the adults population data $(t_i, X_i)$ given as

4.1.1.   Fitting single models over the modelling time

For all the models, we use the pre-estimated and parameters ranges in Table 4. Figure 3 shows the model fits to data from UNAIDS website for Botswana from 1990–2019 for the three epoch models.

Figure 3(a) shows the fitting the epoch 1 model ($SIA$) to the data using the inital conditions $S(0) = 722 000$, $I(0) = 35 000$ and $A(0) = 5000$, the SSE is equal to $8.1 \times 10^3$. Figure 3(b) shows the fitting of the grand model ($S I I_W I_T A$) to the data with inital conditions $S(0) = 722 000$, $I(0) = 35 000$, $I_W = 0$, $I_T = 0$ and $A(0) = 5000$. The SSE equal to $3.5 \times 10^3$. Figure 3(c) shows the fitting of epoch 3 model ($S I I_T A$) to the data with inital conditions $S(0) = 722 000$, $I(0) = 35 000$, $I_T = 0$ and $A(0) = 5000$ to data. The SSE in this case is equal to $7.01 \times 10^3$. The best-fit curves for all the three models do not simulate the data ideally, as shown by Figure 3.

4.1.2.   Fitting using the piece-wise system

We note that UNAIDS collected the data for almost three decades for various responses to HIV/AIDS. In order to achieve a simulated best fit to a series of data points that captures the change in HIV treatment policy in Botswana for almost three decades, we propose to calibrate the second epoch model and allow some parameters to vary over the modelling period. We vary these parameters based solely on the implementation strategy of ART in Botswana.The second epoch model fit the data to Botswana's treatment policy for the period 2002–2016. The implementation of antiretroviral therapy in Botswana started in 2002. We, therefore, render all the parameters related to treatment ($\epsilon, \theta, \gamma, \eta_2, \eta_3, \rho_2$ and $\rho_3$) to zero for the period before 2002 because there was no HIV treatment. Botswana launched a "Treat All" strategy in 2016, that is anyone who tests positive for HIV enrols on treatment immediately. We, therefore, render all parameters related to waiting for eligibility to treatment based on CD4 count ($\epsilon, \gamma, \eta_2$ and $\rho_2$) to zero.

Below we show the implementation of shutting down parameters of the grandmodel to obtain the models of epoch 1 and 3 using Matlab;

We thus have a piecewise system,

Although many mathematical modelers have fitted models to data see for instance ^[29], they have ignored the fact that factors influencing disease transmission change during the modelling period. We aim to incorporate changes in policy in the fitting to data. We estimate the model parameters guided by the Botswana statistics, literature and curve fitting algorithm. The natural death rate, $\mu$, is taken to be the reciprocal of the average life expectancy of Botswana ^[30]. The recruitment rate, $\Lambda$, is estimated by considering that our initial modelling time is 1990. The World Population Prospects (www.population.un.org) estimates the total population of Botswana in the year 1990 to be 1 382 000 of which the adults (aged 15+) accounted for 55.2%. Therefore the total adult population, $N = 762 864$. The UNAIDS estimated the total number of adults living with HIV to be 40 000 in year 1990. Thus we have $I+A = 40000$ since the data does not distinguish between infectious and AIDS staged individuals. We thus estimate the initial conditions as follows;

$\bullet$ $N$ = 762 864, adult population found in the year 1990.

$\bullet$ $I(0)+ A(0)$ = 40 000, HIV adults prevalence in 1990.

$\bullet$ $I(0)$ = 35 000, assumed number of infected staged individuals.

$\bullet$ $A(0)$ = 5000, assumed number of AIDS patients.

$\bullet$ $S(0) = N-I(0)-I_T(0)-I_D(0)-A(0)$ = 722 864.

Figure 4 shows the fit of all three epoch models to the data of adults living with HIV/AIDS for Botswana. It shown from Figure 4(a) that the model fits very well to the data and the SSE = $9.07\times 10^2$. The estimated parameters from the fitting are given in Table 4. Figure 4(b) shows a steady increase in the number of infected adults by 2030, provided the HIV/AIDS policy for Botswana does not change.

Figure 5 shows the time series plots of all state variables of the model. Figure 5(a) shows infected increasing from 1990 and starts falling after 2002 before stabilizing. From 1990 to before 2002, the number of individuals not legible for treatment is zero. After 2002 the number of individuals increased starts falling after 2007 and stabilize from 2007 as shown in Figure 5(b). Individuals on treatment kick off from 2002, shown in Figure 5(c), and have been growing since. Figure 5(d) shows an increase of individuals on the AIDS stage, from 1990 to around 2002 before decreasing.

Figure 6 shows the projections of all state variables of the model provided the policy does not change. Figure 6(a), shows that susceptible individuals would continue to decrease steadily. Figure 5(b) shows that the infected population also continue to decrease, infected individuals on treatment would continue to increase every year as shown in Figure 6(c). The infected progressing to AIDS stage would steadily decrease as seen in Figure 6(d).

4.2.   Sensitivity analysis

The estimation of parameters through curve fitting is subjected to variation since the parameters are selected from a range. Sensitivity analysis aims to quantify the influence of parameters variation in the model output. We perform global sensitivity analysis to examine the sensitivity of $R_0$ to variation in parameters using the Latin Hypercube Sampling (LHS) and the partial rank correlation coefficients (PRCCs). LHS is a statistical sampling method that allows for an efficient analysis of parameter variation across simultaneous uncertainty ranges in each parameter ^[31]. PRCCs illustrate the degree of the effect that each parameter has on the outcome. Fixing $\mu = 1/55$ while taking ranges of the rest of the other parameters from Table 4, we compute the PRCCs of the parameters with respect to the $R_0$. The PRCCs for each parameter against $R_0$ obtained after 1000 LHS samples is presented in Figure 7. Parameters that have high absolute values of PRCCs are the most influential to the reproduction number. The sign of the PRCCS gives the qualitative relationship between the input parameters and the output variable, $R_0$. The positive PRCCs imply that as the parameter value increases, the reproduction number will also increase and vice versa. On the other hand, negative PRCCs imply that as the parameter value increases, the $R_0$ will decrease and vice versa. The parameters $\beta, \eta_1, \eta_2, \eta_3, \rho_1$ and $\epsilon$ have the potential to cause a burden to the disease if not mitigated. Whilst the rest of the parameters, when enhanced, can help ease the disease burden.

Figure 7(a) shows that the most influential parameters of the model are the transmission rate, $\beta$ and the AIDS induced deaths, $\delta$. The box plot Figure 7(b) shows the five number distribution of obtained values of $R_0$, the lower quartile, the median and the upper quartile being about 2.19, 4.29 and 6.53, respectively. The best fit parameter values yield an $R_0 = 6.8$. The median value of $R_0$ in African countries is 4.5, the 90th percentile value of $R_0$ is 6.3, and the maximum being 9.5 ^[32]. The best fit parameter values yield an $R_0 = 6.8$ for Botswana, which is within the estimated calculations given in ^[32]. We note that although some parameters may have small magnitudes of PRCCs, they may still make significant changes in the disease spread.

4.3.   Effects of varying $m$ and $\theta$ on the dynamics of HIV/AIDS

To investigate the effects of Botswana's HIV/AIDS policies, we vary the parameters linked to behavioural change and treatment. The first epoch policy was mainly focused on behavioural change. However, the second and third epochs policies included treatment at different rates. We, therefore, vary the parameters $m$ and $\theta$ to evaluate their effects on the dynamics of HIV/AIDS for Botswana.

Figure 8 shows the projected change in the number of HIV/AIDS infected adults. Increasing both the testing and treatment rate and behavioural change. A decrease in HIV/AIDS infected adults is observed as we keep doubling the testing and therapy efforts and behavioural change rate.

Figure 9(a) shows that increasing $\theta$ only, decreases the number of infected adults just from 2002 when treatment was introduced. Figure 9(b) shows that increasing both $m$ and $\theta$ would decrease the infected cases from the beginning of the modelling period. The shaded regions show the number of cases averted by the increase in treatment and behavioral change. More cases would be averted by increasing both $m$ and $\theta$.

Figure 10(a) shows that increasing $\theta$ only, increases the number of infected adults enrolled on treatment and the shaded area shows the number of individuals added into ART. Interestingly, Figure 10(b) shows that increasing both $m$ and $\theta$ would decrease the number of infected adults enrolled on treatment; this is as a result of more averted infected cases shown by Figure 9(b). The shaded regions show the number of treated cases averted by increasing both $m$ and $\theta$ interventions. Therefore, this decrease can be interpreted as a cost-saving measure for the country.

Figure 11(a) shows that increasing $\theta$ only decreases the number of infected adults progressing to the aids stage just from 2002 when treatment was introduced. The shaded regions show the number of averted AIDS cases. On the other hand, Figure 10(b) shows that increasing both $m$ and $\theta$ decrease the number of infected adults progressing to the aids stage from the beginning of the modelling period. This result shows that more AIDS cases are averted by increasing both behavioral change and treatment.

5.   Conclusions

In this paper, a piecewise system of three models is proposed and analysed to track Botswana's HIV/AIDS policy changes. Introducing policy change in mathematical models is critical for a holistic view of an epidemic influenced by policy dynamics. The structure of the models become dynamic to capture the changes in policy and response to the disease. In this study, we looked at how the treatment policies of Botswana evolved. As a result, this led to the development of models over three epochs. In the 1990s, there was no treatment, and many of the models only involved those at risk, the infected and those who would have developed AIDS. With the advent of treatment that depended on the CD4 cell count, the models had to adapt to capture the rollout of ART programmes. Many researchers developed several mathematical models, but the model presented in this paper represents the perfect summary of the general structure of models that capture how the response to the epidemic should be modelled. The treatment policies further developed into the "Treat All" strategy, which assumed that treatment had the potential to reduce the viral load to undetectable levels, resulting in reduced transmission. We also present the model that captures this policy.

Of the three models developed in this paper, the model of the second epoch captures the dynamics of the first and third epochs if some chosen parameters values are set to zero. The analysis of the grand model thus captures the study of the other two models. We, therefore, focused on the analysis of the grand model. The steady states were determined and analysed in terms of the model reproduction number $R_0$. The model is fitted to data from UNAIDS over three decades. It is clear that the use of a piece-wise models, captures the ideal scenario in which the disease transmission process is dynamic, changing as the policies change.

Many of the mathematical models presented with data fitting ignore aspects of changes in policies over time. They thus consider the parameter values that are averaged over the entire modelling period. It is however important to track the changes in the policies over time that will result in changes in the parameter values over time. We fitted the three models and showed that the fitting in which policy changes are not taken into consideration do not perfectly fit to the data. To address the changes in the parameters over time, we then fitted the grand model to prevalence data from UNAIDS on Botswana. Some of the parameters were considered as continuous piecewise functions to generate a perfect fit, as shown in Figure 4. Sensitivity analysis is performed to determine how each of the parameters impacts Botswana's HIV/AIDS epidemic.

The early policies were based on campaigns to target behavioural change, and the later policies were focused on the uptake of testing and treatment. However, the simulations results showed that a combination of policies leads to effective management of the disease. Therefore, the policymakers need to evaluate the policies and use the information to update the new policies adding to the latest research findings and recommendations. Also, it is imperative to ensure that when new policies are implemented, the earlier ones are not left behind but updated. Earlier policies may be updated to meet the latest trends. For example, the testing and treatment enrollment campaigns to reach the 95-95-95 target should also carry a message of prevention from contracting the disease targeting behavioural change as a package.

The models presented here capture the most straightforward scenarios in each epoch. It is thus essential to note that further refinements can be done by including additional compartments to capture the challenges of managing HIV/AIDS in Botswana. In particular, aspects of drop-out from treatment programmes, defaulting, delays in accessing health care services, and excluding the migrant population from accessing free HIV/AIDS treatment before September 2019 should be considered to give a more realistic picture of how policies have impacted the epidemics in Botswana. Our study, however, provides a significant impetus in the modelling of HIV/AIDS dynamics in the presence of policy changes.

A.   Proof of the global stability of the endemic equilibrium point

Proof. We propose a Lyapunov function of the form

where $a_2$, $b_2$, $c_2$ and $d_2$ are positive constants to be determined. Note that $V = 0$ when $(S_1, I_1, I_{W1}, I_{T1}, A_1) = (S^*, I^*, I_W^*, I_T^*, A^*)$ otherwise $V > 0$. We now show that the time derivative of $V$ is negative. Therefore, the time derivative of the Lyapunov function is given by

From the system (3.2), we also have

We express the derivative of the Lyapunov function as the sum of its components, $\dot{V} = \dot{V_{1}}+\dot{V_{2}}+\dot{V_{3}}+\dot{V_{4}}$, and evaluating the components separately, we obtain

Letting

the derivative of $V$ them reduces to

where

By setting

the coefficients of $F, G, H, J, EF, EG, EH$ and $EJ$ become zero and

We note that $\left(2-\frac{1}{E}-E \right)$ is less than or equal to zero by the arithmetic-mean-geometric-mean inequality, with equality if and only if $E = 1$. Also, the following expressions;

are less than or equal to zero by the arithmetic-mean-geometric-mean inequality, with equality if and only if $E = 1$ and $F = G = J = H$. After some tedious algebraic manipulations of replacing the constants $b_2, c_2$ and $d_2$ in Eq (6.4), we draw similar conclusions. We thus have

if and only if $S = S^*$ and $\frac{I}{I^*} = \frac{I_{W}}{I_{W}^*} = \frac{I_{T}}{I_{T}^*} = \frac{A}{A^*}$. Since $S$ must remain constant at $S^*$ then $\dot{S} = 0$. Therefore, $\frac{I}{I^*} = \frac{I_W}{I_W^*} = \frac{I_T}{I_T^*} = \frac{A}{A^*} = 1$. Thus, the invariant set in $\Omega$ is the singleton $(S^*, I^*, I_W^*, I_T^*, A^*)$. Thus the endemic equilibrium $\mathcal{D}_{1}$ is globally asymptotically stable if $R_0 > 1$ by the LaSalle invariance principle ^[33].

Conflict of interest

The authors declare there is no conflict of interest.