Two strains and drug adherence: An HIV model in the paradigm of community transmission.

A two-strain model, comprising of drug-sensitive and drug-resistant strains, is proposed for the dynamics of Human Immunodeficiency Virus (HIV) spread in a community. A treatment model is introduced by taking drug adherence into account. The treatment-free model is analyzed for the effect of treatment availability and drug adherence on disease dynamics. The analysis revealed that for the treatment-free model, at least one strain faces competitive exclusion, and co-existence of both strains is not possible. On the contrary, both strains may co-exist in presence of treatment. The analysis carried out was both local, as well as global. A comprehensive bifurcation analysis showed periodic behaviour and all solutions approached a stable limit cycle for a wide range of parametric values. Overall, we concluded that the treatment availability and drug adherence play a significant role in determining the dynamics of HIV spread. Numerical simulations are performed to validate the analytical results using MATLAB.

Introduction

From the epidemiological perspective, Human Immunodeficiency Virus (HIV) is a RNA retrovirus that progressively weakens the hosts’ immune system by targeting the helper CD4 T cells [1, 2]. The entry of HIV into CD4 T cells happens through the translation of RNA into DNA, using a viral enzyme called reverse transcriptase. The level of severity of HIV infection can be measured by the viral load or the CD4 T cell count. In particular, an HIV-infected patient is classified as having Acquired Immunodeficiency Syndrome (AIDS), if the CD4 T cell count decreases from the normal level (of around 1000 ) to 200 or below [2]. In the absence of any treatment, HIV infection is almost certainly fatal, leading to the death of an infected person within 5–10 years [3]. Despite the continuous progress in the field of treatment and prevention, there is no cure or vaccine for HIV infection yet. However, there are therapeutic interventions like antiretroviral therapy (ART), which reduces the replication of HIV significantly, thereby improving the prognosis for the infected person [4]. Clinical trials and observational studies have shown a substantial reduction in mortality and morbidity, leading to increased life expectancy in HIV-infected patients treated with ART [4-7]. With the enhanced accessibility to medical facilities, around of individuals living with HIV had started ART by the end of 2020 [8].

Drug adherence usually refers to the degree to which a patient adheres to the prescribed medications, and plays a vital role in the control of chronic diseases, in general. The viral suppression of HIV replication is significantly associated with the adherence to ART regimens [9]. Optimal adherence to ART helps in preventing onward HIV transmission to others [10, 11], minimizing the emergence of drug resistance [12, 13] and brings down the HIV related mortality [14]. There are several factors, namely, personal attributes, institutional management, treatment related factors and psychological factors that lead to ART non-adherence, in developing as well as developed countries [15]. The fear of disclosure of HIV diagnosis and subsequent stigma in the society [16-18], low self-efficacy and reading ability [19], as well as mental illness [17, 20] are some major concerning factors related to ART adherence. The situation gets exacerbated when two major diseases collide. In this context, the COVID-19 pandemic has severely affected the ART adherence in HIV patients. On one hand, COVID-19 has interrupted the ART drug supply due to transportation difficulties (as most cities have implemented some form of lockdown), while on the other hand, overstretched healthcare system under COVID-19 resulted in lower quality of clinical care for HIV patients, in addition to suspension of HIV testing [21-23]. It has also increased the cases of mental health issues such as depression and anxiety [24]. In July 2020, WHO announced that 73 countries have warned of being at risk of stock of ART medicines running out because of the COVID-19 pandemic, and 24 countries reported having either an extremely low stock of ARTs or disruptions in the supply of these vital medicines [25]. Therefore, it becomes very crucial to study the impact of drug adherence on the dynamics of HIV infection in a population.

The phenomenon of HIV drug resistance (HIVDR) continuously remains a major clinical and public health issue, because resistant HIV not only increases the vulnerability of the infected patients, but also the likelihood of transmission to others [26, 27]. The drug-resistant mutant generally develops via inadequate drug concentration during the treatment. There are lot of reasons of inadequate drug level such as poor adherence, treatment interruptions, host genetics, or use of sub-optimal drug combinations [13]. HIVDR can be classified as primary (transmitted) and secondary (acquired). Primary HIVDR exists even before the initiation of ART. This type of resistance is usually transmitted from other drug-resistant infected person or acquired during previous treatment to eliminate vertical transmission of HIV [28]. This leads to the failure of first-line regimens, especially when not identified at the time of ART initiation [29]. On the other hand, secondary HIVDR evolves through selection process between virus mutants under sub-optimal drug adherence [13]. The development of drug resistance can be easily represented by the following schematic diagram.

Among people failing ART, the levels of resistance reached up to [30]. ART is highly exposed to poor or sub-optimal drug adherence due to the above discussed factors. As a result, the risk of emergence of drug resistance is increasing, concurrently, with the improvement in the coverage of ART. To control the development of HIVDR, it is extremely important to study the impact of ART and its adherence to the dynamics of drug resistance in HIV infected population.

Mathematical modelling has been useful in epidemiology over the years to determine the major contributing factors to the disease dynamics. Several researchers have analyzed the within-host [2, 31–33] and between-host [3, 34–40] dynamics of HIV spread. The initial efforts towards modelling the HIV dynamics in a community were carried out in [34, 35] and subsequently various refinements have been added to the modelling frameworks. In [3], Cai et al. considered an HIV/AIDS treatment model by dividing the period of infection into the asymptomatic and the symptomatic phases and showed the dependency of disease status on different treatment levels. Naresh et al. [37] proposed a mathematical model for the spread of HIV/AIDS in a population of varying size with immigration of infectives and deduced that the restriction on a direct inflow of infective population can slow down the spread of infection. A different HIV/AIDS epidemic treatment model with nonlinear incidence rate showed that appreciable change in the susceptible individual’s sexual habits reduces both incidence and prevalence of the disease [39]. Recent studies [36, 41–43] have shown the dynamics of multiple pathogen strains of various infectious diseases including HIV. Sharomi et al. [36] developed a two-strain HIV model of six compartments with the inclusion of drug-sensitive and drug-resistant HIV-infected populations. In this study, the authors analyzed the impact of ART on HIV dynamics and concluded that the widespread use of ART, despite the risk of the development and transmission of drug-resistant strain, can lead to a significant reduction in disease burden or even eradicate the HIV infection from a community, under certain conditions. Further, Kuddus et al. [43] investigated a two-strain model for general infectious diseases that have a protracted infectious period with treatment. They suggested that the emergence of drug resistance could be reduced if the treatment rate is sufficient to eliminate the drug-sensitive strain from the population. We note that none of the related prior studies have explored the role of drug adherence in the presence of treatment. Also, the effect of drug adherence on the drug-resistant infected population has not been studied yet. In this study, we develop a novel mathematical model that includes two separate compartments for the infected population, namely, the drug-sensitive and the drug-resistant, and analyze the impact of treatment, and its adherence on the community transmission of HIV. In particular, we focus on the role of treatment availability and its adherence, in the control of transmission and emergence of drug-resistant strain in a community. The key features of relevant models in literature along with the novelty of the model proposed here are summarized in Table 1.

Table 1

Comparison between present work and previously existed studies based on various factors

HIV/AIDS dynamics	With treatment class	Multiple strain study	Incorporation of drug adherence	With drug resistant infected class	References
Yes	No	No	No	No	[37]
Yes	Yes	No	No	No	[3, 38, 39]
Yes	No	Yes	No	Yes	[40]
No	No	Yes	No	No	[41]
No	Yes	Yes	No	Yes	[42]
No	No	Yes	Yes	Yes	[43]
Yes	Yes	Yes	Yes	Yes	Present work

To the best of the knowledge of the authors, a rigorous study of the dynamics of HIV spread under treatment in presence of drug-sensitive and drug-resistant infected population has not be conducted using mathematical modelling. Keeping all these facts in mind, we presented a mathematical model in Sect. 2 and rest of this article is organized as follows. In Sect. 3, we analyzed the treatment-free model, followed by the treatment model. In this section, we discussed non-negativity and boundedness of the solutions, existence of equilibrium points and their local and global stability followed by bifurcation analysis. Section 4 is devoted to the estimation of all the parameters and interpretation of all analytic results with the help of numerical simulations. Finally, in Sect. 5, we concluded this article by discussing various possible biological interpretations of obtained results.

Formulation of mathematical model

The model divides the total sexually active population into four mutually exclusive groups, namely, susceptible (S), infected individuals with the drug-sensitive strain (), infected individuals under treatment (T) and infected individuals with the drug-resistant strain (). We assume that the individuals are homogeneously distributed in the given population. The susceptible population increases through the process of recruitment of individuals, who enter into the sexually active class, at a constant rate . The susceptibles become infected after an effective contact either with the drug-sensitive infected or the drug-resistant infected individuals with the incidence rates of and , respectively. Here and represent the effective contact rates of susceptible with the drug-sensitive infected and the drug-resistant infected individuals, respectively. Generally, is larger than , since most mutant strains are less fit than the wild-type strains. We assume that an effective contact transmits the same type of strain in the newly infected person. The primary reason for the drug resistance in the infected population is its acquisition from the failure of treatment, and not the transmission [44]. Therefore, we do not consider the transmission of resistance from drug-sensitive infected population. Further, individuals in each of the groups die naturally at the rate . The poor economic condition of infected patients, in addition to limited medical infrastructure, restricts the reach of treatment to only a fraction () of the whole drug-sensitive infected population. Moreover, we assume that only a fraction () of total population under treatment adheres to it, as directed. The remaining fraction of the treatment class is non-adherent to the prescribed drugs, which causes this population to develop a drug-resistant strain of the virus. As a result, individuals of the treatment class will transfer to the drug-resistant infected class at the rate of . Also, these remaining fraction of the population in the drug-sensitive infected and treatment compartments die with disease induced death rate of . This population also develops drug resistance because of the sub-optimal adherence to the treatment. Here we are considering only the sub-optimal adherence to the treatment as the reason for the development of drug resistance. The drug-resistant infected population dies with the disease induced death rate of . We also assume that a drug-resistant infected individual is at drug resistance level. In other words, the treatment does not work at all on drug-resistant infected population. A schematic representation for this type of HIV spread is presented in Fig. 1.

Fig. 1

Schematic representation of the model (2.1)

Taking into account the aforementioned assumptions, the population in each group is determined by the following deterministic system of coupled nonlinear ordinary differential equations:where all the parameters are as defined in the preceding description. For simplicity, we define,

Model analysis

Since the state variables of the model represents the human population, the solution of system (2.1) needs to be non-negative (for all non-negative initial conditions) and bounded. We begin with the non-negativity part. Let denote the time at which the susceptible population S becomes extinct. Then and . Therefore any such that and . Hence and . Now, from the second equation of (2.1), we have, for . In an analogous manner, from the third and fourth equation of (2.1), we have, and , respectively.

For the boundedness part, we define . Then,which implies that . This shows that W is bounded by . Therefore S, , T and are also bounded by . From these results, we can now state the following Theorem (Fig. 2).

Fig. 2

Schematic representation of the model (3.1)

Theorem 1

The biologically feasible region defined by:where denotes the non-negative cone and its lower dimensional faces, is positively invariant for the system (2.1) with non-negative initial conditions.

Treatment-free model

We choose and in model (2.1) to analyze the dynamics of HIV spread in the absence of the treatment. Accordingly, the resultant model is given by:where all the parameters have same meaning as in system (2.1). From Theorem 1, we can conclude that for the treatment-free model (3.1), the feasible region defined by:is positively invariant with non-negative initial conditions. Therefore, it is sufficient to consider the dynamics of model system (3.1) in the region .

Equilibrium points and local stability analysis

For the system given by equation (3.1), there exists three equilibrium points, as enumerated below. Note that the interior equilibrium point does not exists, unless , in which case there are infinitely many solutions for different initial conditions.

Disease-free equilibrium point: .

Drug-sensitive equilibrium point: , which exists provided .

Drug-resistant equilibrium point: , which exists provided .

Now, in order to determine the “Basic Reproduction Number”, for system (3.1), we consider, the next generation matrix as,Hence the “Basic Reproduction Number” for the drug-sensitive strain and the drug-resistant strain are given by and , respectively. We denote the ratio of and , as .

Finally, we present the stability analysis for system (3.1), for which we consider the following Jacobian,

The eigenvalues of the Jacobin matrix, evaluated at are , and . So the disease-free equilibrium is locally stable if, and , i.e, and , respectively. This is equivalent to the condition . Note that under these conditions and do not exist.

One of the eigenvalues of the Jacobin matrix, evaluated at is given by . The remaining two eigenvalues are the solutions of the characteristic equation , which have negative real parts provided (this condition is already satisfied from the existential condition of ). Note that the first eigenvalue is negative, provided,

. Therefore the equilibrium point is locally stable if and only if .

One of the eigenvalues of the Jacobin matrix, evaluated at is given by . The remaining two eigenvalues are the solution of the characteristic equation , which have negative real parts provided (this condition is already satisfied from the existential condition of ). Note that the first eigenvalue is negative, provided,

. Therefore the equilibrium point is locally stable if and only if .

The summary of existence and stability conditions for different equilibrium points are presented in Table 2 and Fig. 3.

Table 2

Existence and stability conditions for different equilibrium points

Equilibrium Point	Existence Condition	Stability Condition
λμ,0,0	Always	maxR0(S),R0(R)<1
μ+μ_Sα,λα-μ(μ+μ_S)α(μ+μ_S),0	R0(S)>1	R0(SR)>1
μ+μ_Rβ,0,λβ-μ(μ+μ_R)β(μ+μ_R)	R0(R)>1	R0(SR)<1

Fig. 3

Two-parameter bifurcation diagram showing existence and stability regions for different equilibrium points of model (3.1). The dashed boxes contain total existing equilibrium point while solid boxes contain only stable equilibrium points. The dashed blue line represents

Existent equilibrium points and their stability in different cases are presented in Table 3, for which we define , and

Table 3

Existent equilibrium points and their stability in different cases

Case	In terms of basic reproduction number	Existence of equilibria	Stable (S)/unstable (US)
P_1<P_2<P_3	R0(S)>1, R0(R)>1 and R0(SR)>1	E0(1),E1(1),E2(1)	S={E1(1)} and US={E0(1),E0(2)}
P_1<P_3<P_2	R0(S)>1, R0(R)<1 and R0(SR)>1	E0(1),E1(1)	S={E1(1)} and US={E0(1)}
P_2<P_1<P_3	R0(S)>1, R0(R)>1 and R0(SR)<1	E0(1),E1(1),E2(1)	S={E2(1)} and US={E0(1),E1(2)}
P_2<P_3<P_1	R0(S)<1, R0(R)>1 and R0(SR)<1	E0(1),E2(1)	S={E2(1)} and US={E0(1)}
P_3<P_1<P_2	R0(S)<1, R0(R)<1 and R0(SR)>1	E0(1)	S={E0(1)} and US=Φ
P_3<P_2<P_1	R0(S)<1, R0(R)<1 and R0(SR)<1	E0(1)	S={E0(1)} and US=Φ

Global stability analysis

Theorem 2

The equilibrium point is globally asymptotically stable with

Proof

We define the Lyapunov function as:Then and The derivative of along the solution of system (3.1) givesTherefore, by LaSalle invariance principle, is globally asymptotically stable if

Theorem 3

The equilibrium point is globally asymptotically stable whenever it exists if .

We consider the Lyapunov function as:The derivative of along the solution of system (3.1) givesTherefore, by LaSalle invariance principle, is globally asymptotically stable if

Theorem 4

The equilibrium point is globally asymptotically stable whenever it exists if .

We define the Lyapunov function asThe remaining proof is similar to the proof in Theorem 3.

Bifurcation analysis

Theorem 5

The system (3.1) undergoes through transcritical bifurcation at if . Here, equilibrium points and exchange their stability at .

Let , thenandNotice that one eigenvalue of matrix become zero at and remaining two eigenvalues are negative only if . Now, we choose and as the eigenvectors corresponding to the zero eigenvalues of the matrices and , respectively. We rewrite the system (3.1) as and . Then,andAlso, So,Now, from equations (3.2), (3.3) and (3.4), we can conclude that, Therefore, F satisfies all the transversality conditions (from Sotomayor’s theorem) for transcritical bifurcation at .

,

,

Theorem 6

The system (3.1) undergoes through transcritical bifurcation at if . Here, equilibrium points and exchange their stability at .

The proof of this theorem is similar to the proof in Theorem 5.

In the absence of treatment, there is no development of drug resistance. The only way of getting drug resistance is transmission from already drug-resistant infected population. There is a strong competition between both strains as they share their common resource, namely, the susceptible population. The strongest and the fittest one will win this battle and survive for long time. The basic reproduction number is a good indicator of the overall fitness of a strain. Co-existence of both the mutants is not possible, due to the competitive selection process. This model lacks a range of biological scenarios due to its apparent simplicity. By excluding the possibility of development of drug resistance, only three scenarios of competitive exclusion are possible: (a) both the types of infected population die out, (b) only the drug-sensitive infected population survives and, (c) only the drug-resistant infected population survives, which makes co-existence impossible. But, there are more possible outcomes between the two strains, beyond the competitive exclusion. Thus, it becomes extremely important to study the model system (2.1) that considers the fact that the drug-resistant infected population increases as a result of development of drug resistance, due to sub-optimal drug adherence, along with the transmission from the already existent drug-resistant infected population.

Treatment model

Recall that the model system (2.1) shows the dynamics of HIV spread with treatment. We analyze this mathematical model by determining the equilibrium states, and their stability, followed by the bifurcation analysis.

For the system given by equation (2.1), there exist three equilibrium points as enumerated below: where the parameters a, b and c are as defined in Sect. 2.

Disease-free equilibrium point: .

Planar equilibrium point:

Interior equilibrium point: , where, , , and ,

Now, we determine the “Basic Reproduction Number”, for system (2.1), by considering the next generation matrix as,Hence the “Basic Reproduction Number” for the sensitive strain and the drug-resistant strain are given by and , respectively. In particular, the “Basic Reproduction Number” for the system (2.1) is . We denote the ratio of and , as .

For the presentation of the stability analysis, we consider the following Jacobian,The stability analysis for each of the equilibria is enumerated below: The summary of existence and stability conditions for different equilibrium points of model (2.1) are presented in Table 4 and Fig. 4.

Table 4

Existence and stability conditions for different equilibrium points of model (2.1)

Equilibrium point	Existence condition	Stability condition
E0(2)λμ,0,0,0	Always	max{R¯0(S),R¯0(R)}<1
E1(2)(S¯_1,0,0,I¯_R_1)	R¯0(R)>1	R¯0(SR)<1
E∗(2)(S_∗,I_S_∗,T_∗,I_R_∗)	R¯0(S)>1 and R¯0(SR)>1	A_∗B_∗C_∗-A∗2D_∗-C∗2>0

Fig. 4

Two-parameter bifurcation diagram showing existence and stability regions for different equilibrium points of model (2.1). The dashed boxes contain total existing equilibrium point while solid boxes contain only stable equilibrium points. The dashed blue line represents . The interior equilibrium point is not stable in the whole region showed in above figure, instead it is stable only in a sub-region of above stability region where condition (C) holds

The eigenvalues of the Jacobin matrix, evaluated at are , , and . Therefore, the disease-free equilibrium is locally stable, provided, and , i.e, and , respectively.

The eigenvalues of the Jacobin matrix, evaluated at are given by , and , where and . All the eigenvalues are either negative or have negative real part provided and , which are equivalent to , and the existential condition for , namely . Thus, the equilibrium point is stable, provided .

Finally, the eigenvalues of the interior equilibrium point are given by the roots of the following characteristic equation of the Jacobian matrix, evaluated at : where By the Routh-Hurwitz criterion, all the eigenvalues will have negative real part if, The coefficients of the characteristic equation, in terms of the Basic Reproduction Numbers are given by, The condition (A) has already been satisfied since the interior equilibrium point exists if and only if and . For condition (B), we have which shows that condition (B) also holds whenever exists. So, we are left only with condition (C) to obtain a locally stable interior equilibrium point.

,

,

.

Theorem 7

The equilibrium point is globally asymptotically stable with , where .

From the system (2.1), we havewhich can be integrated to giveAfter using the upper bound of S(t), which is , we obtainIt follows then that as , if Therefore, the hyperplane attracts all solutions of the system (2.1) originating in domain .

Since as for , from the system (2.1), we haveIt follows then that as , if and the hyperplane attracts all solution of the system (2.1) originating in domain . Similarly, if , as , since as and the hyperplane attracts all solution of the system (2.1) originating in domain . Now, it is simple to show that if , then . Therefore, the equilibrium point is globally asymptotically stable when .

Theorem 8

The equilibrium point is globally asymptotically stable whenever it exists if

From the second and last equations of the system (2.1), we have:Now, we divide equations (3.6) and (3.7) by and , respectively, and obtainHere, equations (3.8) and (3.9) leads to the following relation:and immediately obtain the following inequality:Now, integrating both sides of the above inequality yieldswhich we rearrange in the following formFinally, as we take the limit as , since both and are bounded, we get:Therefore, all solutions of the system (2.1) converge to the hyperplane when . Consequently, all solutions of this system converge to the hyperplane also, as proved in Theorem 7.

Next, we construct the Lyapunov function to show the globally asymptotically stable behaviour of the equilibrium point on the hyperplanes and . We consider the Lyapunov function asThen and The derivative of along the solution of system (2.1) givesNow, at the equilibrium point , we haveReplacing this in equation (3.10), we obtainSince the arithmetic mean is always greater than or equal to the geometric mean, we obtain . Therefore, by LaSalle invariance principle, it follows that the equilibrium point is globally asymptotically stable whenever it exists if

Remark

The global behaviour of the interior equilibrium point is uncertain, since we are unable to construct any suitable Lyapunov function for it. However, we have numerically shown that for a given set of parametric values, this equilibrium point is an attractor for solutions starting from a wide range of initial conditions (see Fig. 8).

Fig. 8

Phase portrait for model (2.1) with various initial conditions assuring the global stability of existing equilibrium points (except ). Other parametric values are

Theorem 9

The system (2.1) undergoes through transcritical bifurcation at if . Here, equilibrium points and exchange their stability at .

Since the equilibrium point becomes non-hyperbolic at , therefore linearization is inconclusive. We use the center manifold theory to study the local behaviour of this non-hyperbolic equilibrium point at . It also describes the existence of another equilibrium point , which is bifurcated from the non-hyperbolic equilibrium. The Jacobian matrix of model system (2.1) around the equilibrium evaluated at , represented by , gives one simple zero eigenvalue and all the other eigenvalues have negative real parts if . Therefore, we can apply center manifold theory. We calculate the right eigenvector u and left eigenvector v corresponding to the zero eigenvalue of the matrix , which are given by and . Now, we rewrite the system (2.1) as and , where and . LetFurther, by calculating the partial derivatives of the functions associated with model system (2.1) at equilibrium point and , we get and . Therefore, from Theorem 4.1(iv) of [45], we conclude that the model system (2.1) undergoes transcritical bifurcation at .

Theorem 10

The equilibrium point changes its stability through Hopf-bifurcation when the real part of eigenvalues of intersects the imaginary axis at , where is a root of the equation .

The characteristic equation of the matrix has purely imaginary roots if: and the system becomes structurally unstable in the neighborhood of the equilibrium point , at where is a root of equation . Further, this equilibrium point will change its stability at through Hopf-bifurcation if conditions (A), (B) and (C) are satisfied along with the transversality condition [46], where x is the root of characteristic equation (3.5). Now differentiating equation (3.5) with respect to , we have:Therefore,Since x is purely imaginary at (say ), we get:where and are as defined in previous section. In order to satisfy the transversality condition, the numerator (say ) must be non-zero at . It is extremely difficult to simplify this expression furthermore, therefore we choose a set of parametric values and then investigate the existence of the Hopf-bifurcation numerically. This analysis has been done in the next section.

,

,

Parameter estimation and numerical simulation

In this Section, we shall estimate the parametric values for the models (3.1) and (2.1) based on the real data from the Indian population. Further, we will also verify the analytic results on existence and stability of the equilibrium points and different bifurcations, through several numerical illustrations, in case of both the models.

The recruitment rate () can be estimated as the sum of new individuals, who enter in the sexually active class and number of net migrated individuals, during the whole year, at the initial time. In India, we assume that most of the adults become sexually active during the age interval of 18–30 years. We set the base year to be 2019, for the purpose of our simulation. The number of individuals recruited in the susceptible class during base year is equal to the average net births during the period 1989–2001, in addition to the net migration in the base year. The net births are total number of births adjusted with the infant deaths in a year. According to United Nations- World Population Prospectus [47], the average birth rate and average infant mortality rate in India during time period 1989–2001 was 29.38 per thousand population per year and 78.87 per thousand live births per year, respectively. Further, the average total population of India during this time period was 955.29 million. Also, the net migration into India in 2019 was million.

Therefore, the constant recruitment rate can be estimated as follows:According to United Nations- World Population Prospectus [47], in India, around 7.2 persons died for every 1000 individuals during 2019. Therefore, the natural death rate is 0.007 per year in 2019. The disease induced death rates can be estimated by using the life expectancy of an HIV-infected person. In [48], the authors have concluded that the life expectancy of an HIV-infected patient reduces about 38 years from the general population at the exact age of 20. In India, the life expectancy from birth was 69.50 years in 2019 [47]. From this we can conclude that the life expectancy of an HIV infected person without any treatment would be around 11 years. Therefore, 1 in every 11 HIV-infected person will die every year. Thus, we estimate that, per year.

The transmission rate represents the rate of change of drug-sensitive infected population in presence of unit susceptible and unit drug-sensitive infected population. According to India HIV Estimate 2019 Report, NACO [49], there were 69200 new HIV cases reported with total 2.35 million cases. The population in age group 20–50 can be considered as total susceptible population since it includes most of the sexually active people. According to [47], in 2019, India had 625.19 million people in the age group of 20–50. Therefore, the transmission rate is given by,We have seen that due to the enhanced accessibility to medical facilities, around of individuals living with HIV had started ART by the end of 2020 [8]. Therefore, we assume for most of the numerical simulations carried out. Further, we assume that of the population in the treatment class follows the prescribed treatment properly. The remaining parameters related to the drug-resistant infected class would be varied within a suitable range, considering the corresponding parametric values associated with the drug-sensitive infected class. For each of the parameters, the unit, estimated value and range of variation has been included in Table 5.

Table 5

Units and values (estimate and range) of various parameters

Parameter	Unit	Estimated value	Range of variation
λ	million year^-1	25.31	Fixed
α	million^-1 year^-1	47.134×10^-6	0.000025–0.0003
β	million^-1 year^-1	Assumed	0.00001–0.00007
μ	year^-1	0.007	0.004–0.007
μ_S	year^-1	0.0909	Fixed
μ_R	year^-1	Assumed	0.06–0.2
η	Unitless	0.7	0–1
ϵ	Unitless	0.75	0–1

For numerical representation of all previous discussed analytic results, we fix for model (3.1) and for model (2.1) and vary other parameters to illustrate all the possible existence and stability scenarios for both models. We choose the initial population as 625 million for susceptible, 2.2 million for drug-sensitive infected, 1.5 million for under treatment and 1 million (assumption) for drug-resistant infected individuals, , which is based on the number of individuals in different classes in the base year. We also choose two other initial conditions to validate the global stability results: and for the model (3.1) and and for the model (2.1).

For model (3.1), first we consider and for which and . The solutions are represented in Fig. 5a which points out the local asymptotic stable behaviour of the disease-free equilibrium point , as expected, since and . The population level of susceptible individuals for different initial values converges to while the drug-sensitive and drug-resistant population decrease continuously and eventually get extinct (see Fig. 6a). This ensures the global asymptotic stability of equilibrium point and supports the theoretical result in Theorem 2. In order to illustrate the case and , we choose and which results in and . The dynamics for these parametric values is represented in Fig. 5b which indicates the persistence of susceptible and drug-sensitive infected population and extinction of drug-resistant infected population after a certain time. It is seen from Fig. 6b, that the population levels of susceptible and drug-sensitive infected individuals for different initial conditions exhibit oscillatory behaviour before stabilizing at the level 2082.98 and 109.59, respectively, while the drug-resistant infected population becomes extinct in a very short period of time. This also suggests the global asymptotic stability of drug-sensitive mutant dominant equilibrium point and is in agreement with the theoretical result of Theorem 3. Note that the basic reproduction number of drug-resistant infected population is greater than 1 but still this population gets wiped out from the system. This is because of a higher basic reproduction number of drug-sensitive mutant which makes it dominant in comparison to the drug-resistant mutant. Further, in order to illustrate the case and , we choose and which leads to and . The solution trajectories are drawn in Fig. 5c which shows that the susceptible and drug-resistant infected population saturate to a positive level and the drug-sensitive infected population dies out after some time. The global stability of drug-resistant mutant dominant equilibrium point , as already obtained theoretically in Theorem 4, can be confirmed from Fig. 6c which also depicts an oscillatory behaviour of the susceptible and drug-resistant infected population before progressing to a steady state. The drug-resistant infected population dominates the drug-sensitive infected population due to its higher basic reproduction number. However, the basic reproduction number for both infected populations is greater than one. This analysis indicates that both the infected populations can be controlled by their transmission rate and disease induced death rate. A higher transmission rate of an infected population is in favor of its persistence at a higher level, while a higher disease induced death rate can cause its extinction from the system. We also observed that the coexistence of both infected populations is not possible because they share the only resource available to them, that is, the susceptible population.

Fig. 5

Time series plot for each population of the model system (3.1) with other parametric values as and initial condition [625, 2.2, 1]

Fig. 6

Phase portrait for model (3.1) with various initial conditions assuring the global stability of each existing equilibrium point. Other parametric values as

On the other hand, we vary the parameters to investigate the dynamics of model (2.1). We consider and and get the corresponding basic reproduction numbers as and . The dynamics for these parametric values is represented in Fig. 7a. Clearly, the compartment of infected and under treatment individuals become empty at a very early stage, while the susceptible population increase continuously and eventually reach to its highest level . The phase portraits in Fig. 8a and b confirms the global stability of disease-free equilibrium point . Next, we choose and to illustrate the case when and . In this case, the drug-resistant mutant dominant equilibrium point behaves like a sink and only the susceptible and drug-resistant infected populations are able to survive for long run. The dynamics for this case is illustrated in Fig. 7b and the global stability of is shown in Fig. 8c, d. Note that the drug-sensitive infected population is wiped out from the system after some time, although it has basic reproduction number greater than one. This happens since the drug-resistant infected population has a higher basic reproduction number than drug-sensitive infected population. In other words we can say that the transmission rate of drug-sensitive infected population is not enough for its survival. Finally, we choose and . This set of parametric values results in , and . Figure 7c represents the long time dynamics of the model (2.1) and suggests that each population converges to the interior equilibrium point after some initial oscillations. Further, the phase portraits in Fig. 8e, f show that the phase point converges to the interior equilibrium point for various initial conditions. It also indicates that the equilibrium point is an attractor for a large region around it.

Fig. 7

Time series plot for each population of the model system (2.1) with other parametric values as and initial condition [625, 2.2, 1.5, 1].

Now we choose the parameters as and to illustrate the transcritical bifurcation when equilibrium point exchanges its stability with equilibrium point . Also, we vary the parameter to show the exchange of stability between equilibrium point and , we vary the parameter . In Fig. 9a, we observe that the equilibrium point bifurcates and losses its stability at and a new stable equilibrium point comes into the picture. We also demonstrated this bifurcation behaviour between these equilibrium points with respect to basic reproduction number . For this, we fix and vary the basic reproduction number and by changing common parameters and . Figure 9b depicts that if , only one equilibrium point exists, that is , which is locally asymptotically stable. On , a new equilibrium point emerges and exchanges its stability with this new equilibrium point provided that . We find similar figures of the transcritical bifurcation for model system (2.1) with same parametric values as discussed above, when disease-free equilibrium point changes its stability with axial equilibrium point . These numerical results support Theorems 5, 6 and 9.

Fig. 9

Transcritical bifurcation diagrams with respect to parameters a and b

For the set of parametric values: and , we observe that condition is not satisfied. Therefore, the interior equilibrium point of model system (2.1) is not locally asymptotically stable for this set of parametric values. Our numerical simulation shows that the solution curves of this system admits periodic behaviour and converge to a stable limit cycle. The time series plot of each population for the above set of parametric values is shown in Fig. 10a, which assures the existence of the periodic nature of the solutions and predicts that the interior equilibrium point is unstable. In Fig. 10b, solutions starting from different initial values are converging to a stable limit cycle.

Fig. 10

a Time series plot for each population of the model system (2.1) showing periodic nature of solutions. b Phase portrait for the model system (2.1) with two different initial conditions assuring the existence of a stable limit cycle for the parametric values: and

Considering as a bifurcation parameter, we note that the solutions of equation for are and , which lies in the interval [0, 1]. The set of eigenvalues for these parametric values corresponding to Jacobian matrices and are and , respectively. These sets consist only either negative real or almost pure imaginary eigenvalues. Further, we have and , which is in agreement with the transversality condition (3.11). Here is as discussed in Theorem 10. So, the equilibrium point changes its stability through Hopf-bifurcation at bifurcation points and . Geometrically, Figs. 11 and 12 confirm the occurrence of Hopf-bifurcation at points and . In Fig. 11, the curves before the bifurcation point and after the bifurcation point show the population levels for different values of the parameter . The limit cycles occur between these two bifurcation points and the red coloured dots represent the lowest population level, while the blue coloured dots represent the highest population level, at a specific value of the parameter . In addition, Fig. 12 shows the periodic nature of solutions of system (2.1) for different values of parameter . We can conclude from Fig. 10b, 12 and stability of coexistence equilibrium point that solutions of model system (2.1) achieve steady state at for and while converge to a stable limit cycle for . This confirms that Hopf-bifurcation is supercritical in nature.

Fig. 11

The bifurcation diagram of each population of the model system (2.1) with respect to the Hopf-bifurcation parameter . The other parameters are

Fig. 12

The existence of limit cycle for the model system (2.1) corresponding to variation in Hopf-bifurcation parameter . The other parameters are

Further, we have plotted density plots for each population on the plane to know the best possible scenario to minimize the infected populations and maximize the under treatment population. We used parametric values: for these density plots. Figure 13a suggests that lower treatment rate will minimize the susceptible population. This happens because lower treatment levels lead to a major increase in total infected population which we can observe in Fig. 13b. It is observed from Fig. 13c, d that perfect adherence with medium treatment rate is an optimal scenario for maximizing the population under treatment. On the other hand, Fig. 13e depicts that a medium treatment rate with imperfect adherence lead to a greater risk of emergence of drug resistance, while extremely low or high treatment rate minimizes this risk. In addition, a perfect drug adherence can also hold the drug-resistant infected population at lower levels.

Fig. 13

Density Plot for a S, b , c T, d T and e on plane. The other parameters are

Conclusion

In this study, we formulated a mathematical model to understand the dynamics of HIV spread under treatment in presence of drug-sensitive and drug-resistant infected population. We considered that the drug resistance transmits from drug-resistant infected population to susceptible via direct contact. Also, drug resistance can develop through improper drug adherence. Therefore, transmission is the only source for the newly drug-resistant infected population, in the absence of treatment. For the treatment-free model, we obtained the basic reproduction number for drug-sensitive and drug-resistant strains with the help of the next generation matrix. The existence of equilibrium points and their stability is totally controlled by these basic reproduction numbers. The stability analysis of equilibrium points provides us with three possible scenarios: In all of the above cases, we observe that the co-existence of all populations is impossible. This happens because both infected populations with drug-sensitive and drug-resistant strains share their common source (susceptible population) of growth. The bifurcation analysis shows that the transmission rate and disease induced death rate of both strains are the determining factors for the elimination of the disease from the system. The lower transmission rate or higher disease induced death rate from certain threshold values makes the system disease free after some time. On the other hand, a greater transmission rate or a lower disease induced death rate of a strain from these threshold values is directly related to the increased fitness of that strain.

Both type of infected population dies out when and .

Only drug-sensitive infected population survives when and .

Only drug-resistant infected population survives when and .

In presence of treatment to the infected population, however, we expect some more potential outcomes between two strains, beyond the competitive exclusion. We calculated the basic reproduction number for drug-sensitive and drug-resistant strains for treatment model to analyze it. We analyzed all existing equilibrium points and their stability. The treatment model also shows three different possible scenarios for infected population: From this, we can conclude that the system becomes disease free at lower reproductive capacities of both the strains. At the same time, a higher basic reproduction number is not enough for the survival of the drug-sensitive infected population, since the drug-resistant infected population out-competes the former, if . This could be a result of high transmission rate or low disease induced death rate of drug-resistant infected class and a higher number of infected individuals in treatment class. These results are consistent with observations in a prior two-strain disease model [43]. We also found that the population in the drug-resistant infected class faces survival risk only at their lower reproductive capacity, in the absence of drug-sensitive infected population. We observed that the drug resistance infected population remains in the system, even at its lower reproductive capacity, if the number of individuals in drug-sensitive infected and treatment class are at a positive level. This shows that a better treatment availability is responsible for the emergence of drug resistance. Therefore, one could question the sustainability of the treatment of the drug-sensitive infected population at the expense of an increment in the drug-resistant individuals in the community [50]. However, during the initial phase, such resistant strains have low fitness levels and are less transmissible than their drug-susceptible counterparts [51]. So an early awareness campaign in non-adherent patients may reduce the risk of emergence of the drug-resistant strain. From the existence condition of the interior equilibrium point, we conclude that the co-existence of both the infected populations, under better treatment availability, demands a high transmission rate in drug-sensitive infected population compare to the other one.

Both type of infected population dies out when and .

Only drug-resistant infected population survives when and .

drug-sensitive and drug-resistant infected population co-exist in the system when and .

Finally, we have performed an extensive analysis of possible bifurcations for treatment model. We found that the system undergoes a transcritical bifurcation at a critical transmission rate of drug-resistant infected population. The system remains at disease-free state if the transmission rate remains below this threshold value, without which the system goes to the only drug-resistant strain endemic equilibrium point. A lower disease induced death rate of drug-resistant infected population also contributes in this transformation from disease-free to endemic state of the system. We also found that the co-existence endemic equilibrium point changes its stability through Hopf-bifurcation at a critical value of parameter and the system shows periodic behaviour. Numerical simulation suggests that all solutions converge to a stable limit cycle for a range of parametric values of . Extremely low or high treatment availability makes the co-existence endemic equilibrium point stable. On the other hand, a large range of medium treatment availability forces the system to show a periodic behaviour, eventually approaching to a limit cycle. In this case, it is extremely difficult for policy makers to reach any concluding strategy to control the spread of the disease if the ratio of number of individuals in each class changes continuously. Therefore, a medium range of treatment availability is not an efficient way to control the disease. We also observed that both basic reproduction numbers are independent of the proportion of population which properly adheres to the treatment, but the stability condition of co-existence endemic equilibrium point depends on it. Sharomi et al. [36] concluded that the widespread use of ART could significantly reduce the disease burden, despite the risk of the development and transmission of drug-resistant strain. But, our study suggests that the drug-adherence to the treatment also plays a vital role in determining the dynamics and is responsible for periodic behavior. The density plot for treatment model shows that a medium level of treatment availability with sub-optimal drug adherence is suitable template for emergence of the drug resistance. This result can be a stepping stone for public health debates [50] regarding the potential benefits and dangers of providing ART to the infected population at risk for the emergence of drug-resistance in patients.

In summary, the present study shows that without treatment, there is no emergence of drug resistance which results in competitive exclusion of one of the infected population. Co-existence of both type of infected population is impossible. On the contrary, when treatment is provided to the infected population, the drug-sensitive infected population either face competitive exclusion or co-exist with drug-resistant infected population. The number of infected patients under treatment and patients who follow the treatment properly are the main determining factors of the dynamics of disease spread. Also, a higher transmission rate of drug-sensitive strain with enough treatment availability and its improper adherence ultimately supports the growth of the drug-resistant infected population. Our analysis shows that policy makers should focus on providing better conditions to the patients to adhere their treatment properly. Otherwise, a higher treatment availability with sub-optimal drug adherence could lead us to a danger of increased number of drug-resistant infected population, or an endemic of drug-resistant infection in the worst of conditions.