Mathematical analysis of sex-structured population model of HIV infection in Kenya

ABSTRACT In this paper, we develop a mathematical model describing the dynamics of HIV transmission by incorporating sexual orientation of individuals. Equilibrium points and the basic reproduction number are derived. The basic reproduction number provides a threshold that determines whether or not the disease fades away. The model, described by non-linear ODEs, shows existence of unique disease-free and disease-persistent equilibria. Least squares curve fitting is presented to quantitatively investigate the trend of infection within each gender. The results are indicative of a higher infectivity in the female population. We further investigated the effect of the introduction of pre-exposure prophylaxis (PrEP) on the dynamics of the HIV. Our results show that the introduction of PrEP has had a positive effect on the limitation of spread of HIV. Sensitivity analysis results show that control of effective contacts can result in control of the disease across gender divide. The model provides a unique opportunity to influence policy on HIV treatment and management.


Introduction
HIV remains a major global health problem affecting approximately 70 million people worldwide causing significant morbidity and mortality (WHO, 2018). In Kenya, the number of persons living with HIV was approximated to be 1.6 million in the year 2016, with women accounting for 910,000 (OPTIONS, 2016;UNAIDS, 2016a). The main transmission route is through heterosexual means with estimated 30% of the new infections deeply rooted among sex workers (National Aids Control Council of Kenya [NACC], 2014a). High HIV infections among sex workers in Kenya is attributed to unprotected sex. According to Shields (2012), most sex workers are constantly harrased by the police in which they are physically and sexually abused for carrying condoms. In addition, the sex workers have no power to negotiate for safe sex. This is due to the fact that clients may decline to pay if they have to use condom and hence use intimidation or violence to force unsafe sex (Ghimire, Smith, van Teijlingen, Dahal & Luitel, 2011). Furthermore, the clients may offer more money for unprotected sex, a proposal that is unlikely to be rejected by the sex workers. According to National Aids Control Council of Kenya (NACC, 2014b) an estimated 29.3% of female sex workers were living with HIV in 2011. In addition, the findings from the Sex Workers Outreach Project reported a HIV prevalence of 30% among female sex workers and 40% among male sex workers in 2011 (UNAIDS, 2015).
Young women are more than three times more likely to be exposed to sexual violence than young man (Kenya National Bureau of Statistics [KNBS], 2015). Among women aged 15-19 years, HIV prevalence was 23.0%, compared to 3.5% among men of the same age (KNBS, 2015;Voeten, Egesah, Varkevisser, & Habbema, 2007). Approximately 33% of the girls in Kenya have been raped by the time they reach 18 years with about 22% aged 15-19 years reporting that their first sexual encounter was forced (NACC, 2014a). In addition, women continue to face discrimination in terms of access to education, employment and healthcare (UNAIDS, 2018). Consequently, men often dominate sexual relationships, with women not always able to practice safer sex even when they know the risks involved (UNAIDS, 2015(UNAIDS, , 2018. For instance, in 2014, 35% of adult women (aged 15-49) who were or had been married had experienced spousal violence and 14% had experienced sexual violence (UNAIDS, 2018). Although progress is being made to control new HIV infections, Kenya still remains one of the six high HIV burdened countries in Africa (UNAIDS, 2014). The main treatment received by HIV patients is the uptake of ART which refers to the use of a combination of three or more ARVs to achieve the suppression of the viral load (WHO, 2014). The drugs do not kill or cure the virus and is a lifelong treatment for HIV patients to reduce the risk of dying (Williams, 2014;Williams, Hargrove, & Humphrey, 2010). In its efforts to reduce HIV infections, in 2015, Kenya began to adopt 2015 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 by end of 2015 (UNAIDS, 2016b).
In the recent years, quite a number of mathematical models have been developed to investigate the dynamics of HIV transmission and spread in the population (Athithan & Ghosh, 2014;Baryarama, Mugisha, & Luboobi, 2006;Doyle, Greenhalgh, & Blythe, 1998;Kaur, Ghosh, & Bhatia, 2014;Okango, Mwambi, & Ngesa, 2016;Omondi, Mbogo, & Luboobi, 2018b;Van Sighem et al., 2012;Wodarz & Nowak, 2002). However, to the authors knowledge, none of these existing mathematical models has explicitly considered incorporating sexual orientation of individuals using real-time series of HIV infection trend data. For example, the mathematical models in Garira (2007a, 2007b) describe the heterosexual interactions of males and females using integro-differential equations with a time delay due to incubation period. While these two models incorporated the effects of male and female condom use as the main mode of preventing HIV infection, they did not consider applying real-time surveillance data to establish the trend of infection. Additionally, in our modelling attempt, we consider a scenario in which individuals with HIV are linked to ART treatment. A study by Mukandavire, Chiyaka, Garira, and Musuka (2009) modelled a sex-structured model for heterosexual transmission of HIV/AIDS with explicit incubation period and provided an in-depth and complete qualitative mathematical analysis. However, there was no numerical results to show the effect HIV prevention intervention measures as well as trend within these two sexes. Bhunu, Mushayabasa, Kojouharov, and Tchuenche (2011) derived a HIV model and examined the effect of counselling and testing coupled with decrease in sexual activity on HIV spread in resource-limited communities. The results from this model suggested that effective counselling and testing have a great potential to partially control the epidemic especially when HIV positive individuals willingly withdraw from sexual activities.
While these models elucidated the importance of prevention to reduce HIV burden, they did not take into account the fact that the immediate linkage of HIV patients to ART treatment can reduce the risk of dying as well as preventing further occurrence of new infections. Our proposed model aims to integrate data on HIV within males and females in Kenya so as to develop a better understanding of the epidemic. In this study, we analyse the model and apply it to data for male and female. Our main goals are: first, to develop a mathematical model that takes into account the treatment of HIV patients with ART and carry out qualitative mathematical analysis of the model and find the necessary threshold for controlling the spread of the disease. Second, to fit the model to observed data of new infections to show the trend between males and females using the parameter values that produce the best fit to the data. The model can help in planning and designing effective control measures to reduce the risk of new infections. In addition, the model and its predictions provides a framework for evaluating the success or failure of the current HIV intervention measures.

Model formulation
The model proposed in this work is an extension of the previous models studied in Mukandavire et al. (2009);Garira (2007a, 2007b). We keep the model as simple as possible consistent with the available data. The model proposed here represents the sexually mature age group in Kenya (age 15 years and over). It is believed that this is the age group that is responsible for the spread of HIV (UNAIDS, 2015). It is assumed that the male and female populations are divided into compartments described by time-dependent state variables. These compartments are: Susceptible (S i ), infected individuals with CD4 cell counts ≥ 350/μL who are considered to be new infections (I i1 ) and infected individuals with CD4 cell counts < 350/μL (I i2 ). The infected individuals are then enrolled on ART which is divided into two depending the CD4 cell counts, that is, T i1 and T i2 . This is consistent with the WHO recommendation of immediate treatment with ART of HIV to attain viral load suppression (AIDS, 2015;Williams, 2014). Here, i = m, f denoting male and female, respectively. The total variable population at time t , denoted by N(t) is described by (1). where The recruitment in the population is at the constant rate of which a proportion α are assumed to be males and (1 − α) are assumed to be females. The newly recruited individuals enter the susceptible compartment S i . Each individual compartment goes out from the dynamics at natural mortality rate μ. The susceptible in either male or female populations is decreased following infection, which can be acquired via effective contact with an infectious person at the respective rates given by λ m and λ f . These rates are obtained as follows the probability that a male or female person chooses a particular partner can be assumed as 1/N m and 1/N f , respectively. Thus, a male or female receives in average c mj N f /N m and c fj N m /N f partners per unit of times, respectively. Then, the infection rate per susceptible male or female are, respectively, given by Here, k = 1,2 representing the two stages of infectious population, β ij , where i = m, f and j = 1, . . . , 4, are the probabilities of HIV transmission per partnership and c ij are the rates at which an individual acquires sexual partners associated with the four infectious compartments respectively. Thus, the expressions in (2) can be simplified to the following Infectious individuals in class (I i1 ) progress to class (I i2 ) at rate γ h , for h = 1,2. Infectious individuals in classes (I 1 i ) and (I i2 ) move to (T i1 ) and (T i2 ) at a constant rate τ p , for p = 1, . . . , 4, as consequence of treatment with ART. The infectious individuals on ART in class (T i1 ) progress to (T i2 ) at a rate ψ h as a result of decline in the immunological status. Similarly, infectious individuals on ART in class (T i2 ) progress to (T i1 ) at a rate ω h as a consequence of improvement in the immunological status. Infectious individuals in classes I i2 and T i2 may die as a consequence of infection, at a disease-induced death rate δ 1 and δ 2 , respectively.

Model assumptions
The following key assumptions have been made.
(i) 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. (ii) The standard HIV transmission incidence has been used to model the disease transmission. This is the form most commonly used for sexually transmitted diseases (Hethcote, 2000). (iii) There is homogeneous mixing and the transmission of HIV is assumed to be mainly through heterosexual means. (iv) 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. Furthermore, AIDS patients are usually too ill to remain sexually active and they are unable to transmit HIV through sexual activity. In 2016, Kenya issued full regulatory approval of PrEP, becoming the second country in Sub-Saharan Africa after South Africa, to make such approval (UNAIDS, 2016b). PrEP is used by people who do not have HIV but are at high risk of acquiring it to prevent HIV infection (AIDS, 2016; HIV/AIDS, 2016). The successes of PrEP use in Kenya 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, 2016b). The challenges posed by the continued occurrence of new infections call for a better understanding of the disease transmission and development of effective strategies for prevention and control of the spread of HIV. Therefore, we add control 0 ≤ φ ≤ 1 to measure the effectiveness of PrEP in the prevention measurements against acquiring HIV infection. Thus the infection terms given in (3) are, respectively, modified as follows ( 4 ) Figure 1 shows the schematic diagram for the compartmental model structure.

Model equations
The above assumptions lead to the following 10 non-linear system of ordinary differential equations.
The description of state variables and parameters of model (5) are given in Tables 1 and 2, respectively. The total population of both male and female sub-populations, respectively, evolve according to the following The solutions to each of the equations in system (7) are given by respectively. Here, N 0m and N 0f are the initial populations for the males and females.
The solution of all the equations from system (5) all remain non-negative for all t ≥ 0. The total population {N m ; N f } in each of the sub-populations is bounded by α /μ and ((1 − α) )/μ, respectively. The total population in both sub-populations is given by N = N m + N f such that the evolution of the total population over a specified period is given by It can easily be seen that lim sup t→∞ N ≤ μ .
Therefore, the phase space of the system (5) is given by The solutions in are all non-negative and bounded. Hence the domain of biological significance is positively invariant and attracting. Therefore, all solutions starting in remain in .

Equilibrium points
To better understand the dynamics of the proposed model, we begin by examining the model's behaviour about the steady states. This analysis is crucial for identifying the parameters of the model that help to achieve a HIV-free state. Since the rate of change in each compartment is constant at equilibrium, we set the right-hand side of equations in the system (5) to zero and solve for the state variables in terms of the infection terms λ mq and λ fq and obtain the following , , , , Substituting the expressions for the state variables (8) into infection terms in (4) then simplifying, the following polynomial is obtained. where The expression for R 0 is computed in (13). Note that the case λ * mq = 0 corresponds to the disease-free equilibrium which can be expressed as The disease-free equilibrium is refers to a scenario in which HIV infection is not present in a population. The solutions to the remaining part of (9), given by Equation (11) gives possible disease-persistent state (E 1 ).
The existence of disease-persitent state refers to a scenario in which HIV infection persists in the population. We note that A 0 > 0 if R 0 < 1 and expression (11)

The basic reproduction number
The basic reproduction number is defined as the average number of new infections generated by one infected individual, via heterosexual means, in a completely susceptible population (Van den Driessche, & Watmough, 2002). The basic reproduction number, R 0 , which is a measure of the infection on the susceptible populations, for the system (5) is obtained using the next generation method described in Van den Driessche, and Watmough (2002). This method has been explored in many papers, for instance, Feng, Ruan, Teng, and Wang (2015); Munz, Hudea, Imad, and Smith (2009);Nyabadza, Njagarah, and Smith (2013); Omondi, Nyabadza, and Smith (2018a); Omondi, Nyabadza, Bonyah, and Badu (2017). Following the explanation given in Van den Driessche, and Watmough (2002), we obtain F which is the matrix of new infections and V which is the matrix of transfers between compartments given by where If there is an existence of an infection in the population, then the corresponding threshold number for the persistence or eradication of HIV obtained from the spectral radius of FV −1 is given by where The components of the basic reproduction numbers, R 0i , for i = m,f, are explained as follows:-(i) R 0mj , for j = 1, . . . , 4, represents the basic reproductive number of the associated with each category of the male infected patients in compartments I mk and T mk , for k = 1,2, when introduced into a population, respectively. (ii) R 0fj , for j = 1, . . . , 4, represents the basic reproductive number of the associated with each category of the female infected patients in compartments I fk and T fk , for k = 1,2, when introduced into a population, respectively.
The term 1 = ψ 1 ω 1 /Q 3 Q 4 indicates the fraction of male individuals in compartments T m1 and T m2 who move from either the compartments due to either improvement or decline in their immunological status. The same explanation applies to the term 2 = ψ 2 ω 2 /Q 7 Q 8 in the female population. The values (1/Q 3 ), (1/Q 4 ) indicate the average time male individuals in compartments T m1 and T m2 stay in their respective compartments. A similar explanation is given for the values (1/Q 7 ), (1/Q 8 ) in the female population. Therefore, (1 − 1 ) and (1 − 2 ) refer to the fraction of individuals who do not cycle between compartments T i1 and T i2 in the male and female populations, respectively. Note that the square-root arises due to the fact that it takes two generations for infected individuals to produce new infections.
From Theorem 2 in Van den Driessche, and Watmough (2002), fundamental results about the equilibria analyses of the system (5) are given by the following propositions.
Proposition 3.1: The disease-free equilibrium point E 0 is locally asymptotically stable when R 0 < 1 and unstable otherwise.

Remark 3.1:
It is important to note that If R 0 < 1, then on average, an infected individual produces less than one new infected individual over the course of its infectious period and the infection cannot grow. Conversely, if R 0 > 1, then each infected individual produces, on average, more than one new infection. As a result, the infection will spread and become endemic in the population.

Proposition 3.2:
The disease-persistent equilibrium point E 1 given by the solution of expression (11) is locally asymptotically stable in when R 0 > 1 and unstable otherwise.

Proof:
The local stability of the disease-persistent equilibrium point E 1 is proved based on the centre manifold theory as described in Castillo-Chavez and Song (2004). We avoid restating the theorem and compute the components of a and b as explained in the theorem. Let ϑ = β m1 c f 1 be our bifurcation parameter so that for Furthermore, to linearize the system (5), we define S m = x 1 , I m1 = x 2 , I m2 = x 3 , T m1 = x 4 , T m2 = x 5 , S f = x 6 , I f 1 = x 7 , I f 2 = x 8 , T f 1 = x 9 , T f 2 = x 10 , and dx dt = f (x, ϑ), f : R 10 × R → R and f ∈ C 2 (R 10 × R).
Thus, we linearize the system (5) at disease-free equilibrium with the bifurcation parameter ϑ to obtain where It is important to note that with the bifurcation parameter expressed in (14), the linearised system (15) has a zero eigenvalue while the rest of the eigenvalues are negative. The left eigenvector of (15), V = (v 1 , v 2 , v 3 , v 4 , v 5 , v 6 , v 7 , v 8 , v 9 , v 10 ) and the right eigenvector W = (w 1 , w 2 , w 3 , w 4 , w 5 , w 6 , w 7 , w 8 , w 9 , w 10 ) T both associated to the eigenvalue zero are solutions of the linearised system (15) , v 10 = Q 5 Q 6 (ω 2 β f 3 c m3 + Q 7 β f 4 c m4 ) (Q 9 + Q 10 + Q 11 ) , while the associated right eigenvectors are given by , w 6 = −Q 1 Q 5 Q 6 (μω 2 + Q 7 (δ 2 + μ)) μ(1 − φ)(Q 9 + Q 10 + Q 11 ) , According Castillo-Chavez and Song (2004), the local dynamics of the system (5) around zero is governed by the the signs of a and b, where To compute a and b we, respectively, evaluate the non-zero second-order mixed derivatives with respect to variables and the non-zero partial derivatives with respect to bifurcation parameter. These are given by , Now, substituting the expressions into (16) and simplifying, we obtain From expressions in (18), it can be seen that a < 0 and b > 0 for all parameter values. Therefore, the disease-persistent equilibrium point E 1 is locally asymptotically stable close to R 0 = 1.

Parameter estimation
Some of the parameters ranges used in this paper have been estimated from previously published articles, while others are estimated intuitively. The baseline parameter values are obtained through curve fitting are presented under the caption of Figure 4. The full list of parameter ranges used in the simulation is given in Table 2.

HIV data
The data were obtained from Kenya Health Information System (KHIS) (KHIS, 2017). The data used represent new HIV infection for both males and females in Kenya. Data are collected routinely on a monthly basis and was retrieved for the period beginning January 2011 to December 2017. Only variables of interest were pulled out to excel spreadsheet. Data were stored in excel and thereafter analysed in R. The pictorial representation of the raw data is given in Figure 2.

Initial conditions
To fit and validate the model (5) to the data on reported cases of new HIV infections in Kenya, the initial conditions are set as follows. The total population of Kenya at the end of 2010 was estimated to be 41.35 million according to KNBS (2015). Of this 57.6% of the population was aged 15 years and above. Therefore, the total population (N) considered in this paper is 23.8176 million. The population of male (N m ) is estimated to be 49% and that of female (N f ) is estimated at 51%. The reported number of adults living with HIV in January 2011 was 1.2 million with approximated 661,515 on ART treatment (NACC, 2017). From here, the number of new HIV infections (I m1 and I fm ) have been estimated based on the data reported in January 2011 to be 500 and 530, respectively (KHIS, 2017). The susceptible population for each gender is estimated using the following relations The breakdown of the initial populations used in the curve fitting is given in Table 1. We therefore resort to curve fitting and numerical simulation to understand the effect of PrEP in the evolution of HIV.

Data input description
Input data were summarized using error bar and density plot to illustrate the extent of variability and presented in Figure 3. For process indicators, the mean number of infections are reported and the 95% confidence intervals are presented using the error bar in Figure 3(a). The results suggest that there appears a significant difference in HIV infections between males and females. This is supported by the non-overlapping standard deviation error bars. Furthermore, it can easily be seen that the mean number of new cases of infections in the female population is higher than the mean number of new cases of infection in the male population. From result in Figure 3(b), it can clearly be seen that the data for the two sets do not follow a normal distribution. Thus, there is need to perform further tests to establish any significance difference. To establish whether there is significant difference in the male and female infections, Mann-Whitney U-test is used. This is a non-parametric test that is used to compare means of two groups that do not follow a normal distribution as suggested by the results in Figure 3(b). The results in Table 3 show that there is significant difference in the mean infections between male and female given that the p-value = 3.686e − 10 < 0.05, performed at 95% confidence level.

Curve fitting
In this section, we fit system (5) to data to determine the trend of HIV in male and female populations. Curve fitting is a process that allows us to quantitatively estimate the trend of  the outcomes. The curve fitting process fits equations of approximating curves to the raw field data. However, for a given set of data, the fitting curves of a given type are generally not unique. Thus, a curve with a minimal deviation from all data points is desired. This bestfitting curve can be obtained by the method of least squares. In this method, the parameters not known are approximated through minimization of the sum of the squared deviations between the data and the model. It minimizes the sum of squared distances between the observed values and the model values. This can be mathematically expressed as where θ i = (y −ȳ) and n refers to the data points and RSS refers to the sum of square error estimate which is assumed to follow a normal distribution. The FME package (A flexible modelling environment for inverse modelling, sensitivity, identifiability and Monte Carlo Analysis) in R is used to fit the model to data. A R code is used in which, the parameters with unknown values are given lower and upper bounds from which the set of parameter values that produce the best fit are obtained. The following parameters were fixed at the following values: = 0.0014N, μ = 0.00139, δ 1 = 0.00915, δ 2 = 0.00725, c m1 = c f 1 = 3, c m2 = c m3 = c m4 = c f 2 = c f 3 = c f 4 = 2, α = 0.49. The parameter ranges/values in Table 2 Figure 5(b) show that infection will decline towards 2030 when PrEP uptake is maintained slightly above 40%. Note that the government of Kenya approved PrEP uptake in 2017 and much of its successes are yet to be reported. The effect of introduction of PrEP is seen in Figure 5(b) where there is a sudden fall in 2017 following the government approval. Thus, there is need to emphasise on preventive measures through educational campaigns and social programmes that ensure minimal or reduced infections. Although, the population shows a general fall in HIV infections, the female population will still continue to be disproportionately affected by the epidemic compared to the the male population as reflected in Figure 5(a). This finding agrees with the finding in the report by Kenya National AIDS Control Council (NACC) (Kenyan Ministry of Health [MOH], 2016;NACC, 2016) .
In order to establish the correlation between the parameters with respect to variables I m1 and I f 1 , we present pairs plot of the markov chain monte carlo (MCMC) samples for the model parameters. As seen Figure 6, the scatter plot matrix in the upper panel describes the pairwise relationship between parameters with corresponding correlation coefficients shown in the lower panel. The marginal distribution for each parameter is shown on the diagonal. The scatter in blue and green correspond to I m1 and I f 1 , respectively. It is shown in Figure 6 there is a negative correlation between β m1 and ω 1 as well as β f 1 and ω 2 . Implying that an increase in the drug efficacy (ART) results to a decrease in the infection terms. This is particular true with regards to treatment of HIV due to the fact that an increase in the drug efficacy results to an increase in viral load suppression which in turn lowers the chances of infection through heterosexual means.

Sensitivity analysis
We perform sensitivity analysis to examine the output's (basic reproduction number) response to the simultaneous variation of the parameter values within a range in the parameter space is described in Table 2. Following the work in Marino, Hogue, Ray, and Kirschner (2008), Wu, Dhingra, Gambhir, and Remais (2013), and Stein (1987), we use Latin Hypercube Sampling (LHS) to determine the Partial Rank Correlation Coefficients (PRCCs) with 5000 simulations per run. PRCC takes values between −1 and +1 in which the sign indicates how the model output is qualitatively related to each model parameter. The parameters are assumed to be random variables with uniform distribution with  their range values given in Table 2 and point values given under the caption of Figure 4. We observe from Figure 7 that the parameters with the greatest potential to increase the HIV infection are the effective person to person contact rates. Moreover, the uptake of PrEP, φ is the parameter with the greatest potential to make the epidemic better when increased. This is supported by the results in Figure 5(b).

Discussion
In this paper, we have analysed a sex-structured population model and studied the HIV infection trends in males and females. The model has assumed that the main mode of HIV transmission is heterosexual. The basic reproduction number and equilibrium points are computed. From our computation, it therefore suffices to deduce that the model exhibits two equilibria namely the disease-free equilibrium and the endemic equilibrium. The data representing new cases of HIV infection for both males and females has been extracted from Kenya Health Information System (KHIS) and analysed. The descriptive statistics of the data shows that there exists a higher number of cases of infection in females as opposed to males. This difference has been established to be significant through Mann-Whitney Utest. The model has been fitted to data using least squares method in R. The trend shows that the females are still disproportionately affected with HIV as compared to males. In order to establish the impact of the recent roll-out of PrEP, we investigated its role in limiting HIV infection. We fixed the rate of PrEP use to zero for the period before May 2017 when the PrEP use was launched in Kenya and after May 2017, the rate of PrEP use was fixed at 0.40 representing 40% coverage. From the trend projection to year 2030, it suffices to conclude that PrEP plays an important role in reducing the number of new cases of HIV. We notice that when the value of parameter (φ) is fixed at 0.40, the cases of infection declines towards 2030 to a near complete eradication of HIV. This implies that controlling and eventual eradication of HIV in Kenya requires aggressive campaigns by the Kenyan government in favour of PrEP use. Furthermore, sensitivity analysis has been carried out using Latin Hypercube Sampling (LHS) technique. It is seen that the model output (basic reproduction number) is highly sensitive to the effective contact rates suggesting that efforts made to reduce the contacts between uninfected individuals and the infected individuals will be most appropriate in limiting the occurrence of new infections.
In conclusion, we show that prevention of HIV infection still remains the most vital way of curbing further spread. HIV patients under ART treatment are possibly capable of aiding the eradication of HIV by convincing their sexual partners of the need to adhere to protection via use of PrEP or any other protection means and ART treatment. The model presented in this paper is a very simplified description of HIV infection in Kenya and therefore it has some cogent limitation. The model presented in this study does not take into account the full stages of HIV. Even though, the model recognizes the fact that there is need for immediate treatment once an individual is found to be positive in line with WHO regulations of 2015 based on viral load suppression, the model did not factor in the viral load levels. This limitation can be circumvented in various ways. First,there is a need to link the model to laboratory experiments for a clearer determination of parameter values based on the viral load of the patients. Second, the development of mathematical models elucidating all the HIV stages will greatly advance our understanding of HIV spread in Kenya. Despite the limitations highlighted, the model results have significant bearings on HIV dynamics and its treatment with ART. the manuscript. RW participated in the mathematical analysis and numerical analysis. LL conceived of the study, and participated in model formulation and helped to draft the manuscript. All authors read and approved the final manuscript.