Modelling Botswana’s HIV / AIDS response and treatment policy changes: Insights from a cascade of mathematical models

: The management of HIV / AIDS has evolved ever since advent of the disease in the past three decades. Many countries have had to revise their policies as new information on the virus, and its transmission dynamics emerged. In this paper, we track the changes in Botswana’s HIV / AIDS response and treatment policies using a piece-wise system of di ff erential equations. The policy changes are easily tracked in three epochs. Models for each era are formulated from a “grand model” that can be linked to all the epochs. The grand model’s steady states are determined and analysed in terms of the model reproduction number, R 0 . The model exhibits a backward bifurcation, where a stable disease-free equilibrium coexists with a stable endemic equilibrium when R 0 < 1 . The stability of the models for the other epochs can be derived from that of the grand model by setting some parameters to zero. The models are fitted to HIV / AIDS prevalence data from Botswana for the past three decades. The changes in the populations in each compartment are tracked as the response to the disease and treatment policy changed over time. Finally, projections are made to determine the possible trajectory of HIV / AIDS in Botswana. The implications of the policy changes are easily seen, and a discussion on how these changes impacted the epidemic are articulated. The results presented have crucial impact on how policy changes a ff ected and continue to influence the trajectory of the HIV / AIDS epidemic in


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 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.

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)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001). The second epoch captures the period from which the universal ART programme was rolled out (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(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".

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 Λ. Susceptible individuals are infected at a rate λ, 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 ϵ of those tested at a rate θ 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 γ. 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 ρ 1 , ρ 2 and ρ 3 respectively. Individuals in the AIDS class are assumed to transmit infection and they also die due to the disease at a rate δ. Individuals in each class are assumed to die naturally at a rate µ. The forces of infection are defined by λ = (λ 1 , λ 2 , λ 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. per time 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 λ 2 (t) is given by Here, β 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 ≤ p ≤ 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 β =βn(1 − p), whereβ 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.

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 θ, γ, ρ 2 , ρ 3 and setting them to zero results an S IA 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].

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, ϵ, ρ 2 and γ, 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, The force of infection λ 3 (t) is given by The parameters are as described in Table 1.

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.

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 λ = λ 2 , we have We thus have From the second equation of the model system (2.1) we obtain dI dt = λ 2 S − (µ + θ + ρ 1 )I ≥ −(µ + θ + ρ 1 )I, whose solution is given by In a similar manner we can demonstrate that 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 The region Ω is thus positively invariant and solutions are bounded. This means that every solution of (2.1) with initial conditions in Ω remain in Ω for all t ≥ 0. The model is thus epidemiologically and mathematically well-posed in the region
The steady state variables expressed in terms of I * are given by, where 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 * = Λ µ . This corresponds to the disease-free state. The model thus has a disease-free equilibrium (DFE) given by

The reproduction number
The reproduction number, R 0 , gives a threshold condition to the stability of the disease free diseasefree 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, θ, γ, ρ 2 , ρ 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 ϵ, γ and ρ 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, D 0 , as follows: Theorem 3.2.1. The DFE, D 0 , is locally asymptotically stable if R 0 < R c 0 < 1 and unstable otherwise. Note that the expression for R c 0 is defined in the next subsection.

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 ν 2 = mσ 1 ϕ 1 (σ 1 − µϕ 1 ), 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, ν 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.
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 C 0 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 c 0 ≈ 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]. The proof of the theorem in given in Appendix A.

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

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 (S IA) to the data using the inital conditions S (0) = 722000, I(0) = 35000 and A(0) = 5000, the SSE is equal to 8.1 × 10 3 . Figure 3(b) shows the fitting of the grand model (S II W I T A ) to the data with inital conditions S (0) = 722000, I(0) = 35000, I W = 0, I T = 0 and A(0) = 5000. The SSE equal to 3.5 × 10 3 . Figure 3(c) shows the fitting of epoch 3 model (S II T A ) to the data with inital conditions S (0) = 722000, I(0) = 35000, I T = 0 and A(0) = 5000 to data. The SSE in this case is equal to 7.01 × 10 3 . The best-fit curves for all the three models do not simulate the data ideally, as shown by Figure 3.

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 (ϵ, θ, γ, η 2 , η 3 , ρ 2 and ρ 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 (ϵ, γ, η 2 and ρ 2 ) to zero.

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 µ = 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 β, η 1 , η 2 , η 3 , ρ 1 and ϵ 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, β and the AIDS induced deaths, δ. 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.

Effects of varying m and θ 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 θ to evaluate their effects on the dynamics of HIV/AIDS for Botswana.  Table 4. 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 10(a) shows that increasing θ 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 θ 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 θ interventions. Therefore, this decrease can be interpreted as a cost-saving measure for the country.   Figure 10(b) shows that increasing both m and θ 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.

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.