Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 2 (2024) https://internationalpubls.com 460 Understanding Drug Resistance of HIV Infection with Multidrug Treatment and Immune Response Dr. Deepmala Kamboj1*, Dr. Kapil Kumar2, Dr. Vijay Chawla3 1*Associate Professor Mathematics, MLN College, Yamuna Nagar, Haryana, India e-mail: dmala.math@mlncollegeynr.ac.in 2MD Community Medicine, ASMO, Mukand Lal Distt. Civil Hospital, Yamuna Nagar, Haryana, India, e-mail:drkambojkapil84@gmail.com 3Assistant Professor Mathematics, Maharaja Agrasen Mahavidyalaya, Jagadhri, Haryana, India e-mail: macmathsdepartment@gmail.com Article History: Received: 05-06-2024 Revised: 03-10-2024 Accepted: 12-11-2024 Abstract: In the present study, a twin-strain mathematical model comprising drug sensitive (wild type) and drug resistant (mutant) strains is proposed for the dynamics of HIV (Human Immunodeficiency Virus). The purpose is to investigate strategies in the multidrug treatment of HIV infection in the presence of drug resistant strains. The HIV infection dynamics is described by a system of nonlinear differential equations, which governs the interaction of uninfected CD4+ T-cells with free virus. The division of infected cells into pre and post-RT classes has been incorporated into the system to explain the biological steps between the viral infection of CD4+ T-cells and production of HIV virions. Further, a combined drug therapy consisting of Fusion Inhibitor (FI), Reverse Tanscriptase Inhibitor (RTI), and Protease Inhibitor (PI) is introduced into the system so as to reduce viral load and thus increase the T-cell population. Continuous viral replication in the presence of drug therapy results in the emergence of variants of drug resistant virus. Thus, there would not be a complete eradication of virus which enhances the risk of the progression of the disease towards AIDS. The system takes into account this fact by introducing the two types of viral strains- drug sensitive and drug resistant strain. The impact of immune response is also considered on this twin-strain model with multidrug treatment. The stability of the steady states emerging in the system is analysed. Conditions are obtained for stability and existence of uninfected and infected steady states. Results from numerical simulations are exhibited to illustrate the dynamic relationship between multidrug therapy administration, the prevalence of drug resistance, the total level of viral production, and the strength of immune response. Keywords: CD4+ T-cells; twin-strain; HIV infection; drug sensitive virus; drug resistant virus; mutation; efficacy; immune response 1 Introduction Human Immunodeficiency Virus (HIV) is a ribonucleic acid (RNA) virus whose replication cycle begins with the binding of the gp 120 protein of HIV on the CD4 molecule (i.e., on the host cell surface). CD4 molecule is found on the surface of dendritic cells, monocytes/macrophages on a subset of T-lymphocytes (also known as CD4+ T-cells) which are responsible for defense function in the immune system. The CD4+ T-cell membrane fuses with the HIV envelope and allows the HIV virus to enter into CD4+ T-cell. After this HIV transfuses the viral RNA (genetic material) into the host cell. Then, on entering the host cell, the viral RNA with the help of an enzyme reverse transcriptase forms deoxyribonucleic acid (DNA). This single-stranded DNA in the host cell converts into double- stranded DNA (viral DNA). HIV inserts this viral DNA into DNA of CD4+ T-cells with the help of an enzyme called integrase. After integrating DNA into CD4+ T-cells, HIV starts producing long chains of HIV proteins. Protease enzymes cut protein chains into smaller pieces to form the structure of a new virus. The newly formed copies of virus exist out of the host cell through budding process and proceed further to infect the new cells and this process continues. Consequently, after detecting the invasion of virus the human body stimulates CD4+ T-cells, which further stimulate CTLs. These CTLs by proliferation and surrounding kill the infected CD4+ T-cells. The count of CD4+ T-cells of an infected individual when reaches below 200 π‘šπ‘šβˆ’3 cells, the stage is then characterized as the mailto:dmala.math@mlncollegeynr.ac.in Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 2 (2024) https://internationalpubls.com 461 onset of AIDS (Acquired Immuno Deficiency Syndrome). Mathematical models played a significant role in developing a better understanding of the internal dynamics of HIV/AIDS, drug therapies and immune responses [1-7]. In 1989, Perelson [1] developed a simple model to explain the interaction between the human immune system and HIV. To identify the behaviour of viral dynamics this model [2] has been extended mathematically. This extended model successfully explained many of the symptoms of HIV/AIDS. Various studies [8-15] used these mathematical models to understand HIV dynamics and to devise drug treatment strategies to counter HIV infection. Since the replication rate of HIV is extremely high therefore its treatment demands simultaneous administration of two or more antiretroviral drugs [16-26]. Antiretroviral drugs interrupt the activities of those enzymes which are essential to complete the different stages of HIV replication cycle. For example, Fusion Inhibitors (FI) prevent the fusion of HIV envelope with the host cell membrane. Integrase Inhibitors block the activity of enzyme integrase that inserts the HIV DNA into DNA of the host cell. Reverse Transcriptase Inhibitors (RTI) directly block the action of reverse transcriptase enzyme and prevent HIV virus replication. Protease Inhibitors (PI) prevent immature HIV from becoming a mature virus by blocking the activity of protease enzymes.Thus, new copies of HIV will not be able to infect new cells. To control the extremely high replication rate of HIV, it is always preferable to devise drug therapies using simultaneous administration of two or more antiretroviral drugs. Most host and viral factors such as nonadherence to the treatment protocol, deleterious side effects, poor drug absorption, etc., are some of the main reasons for drug therapy regimen failure. Out of these, the major factor is found to be the presence or emergence of drug resistant viral strains. Due to HIV infection, infected cells can generate billions of viral particles everyday [7]. The chances of occurrence of mutations are quite high as the process by which RNA genome is reverse transcribed into proviral DNA is highly error-prone [27- 28]. Due to the occurrence of single or combined mutations, there is always a reasonable chance of the generation of drug resistant virus even before the initiation of drug treatment for HIV [16]. Several mathematical models have been designed to analyze the evolution of mutant strains and dynamics of HIV with antiretroviral therapies [16-26]. These models analyzed that the treatment with antiretroviral therapies failed due to preexistence of drug resistant virus. Bonhoeffer et al. [17] suggested that when there is inherited drug resistant virus, then a very effective drug therapy would be able to reduce the HIV viral load at the initial stage. Krischner and Webb [29] obtained an increase in the level of drug resistant viral load during monotherapy treatment of HIV infection. A comparison treatment outcome with drug therapy initiated at different T-cell levels has been made by them. Mclean and Nowak [18] showed that during the course of multidrug treatment, the resistant virus would dominate the wild type virus. Riberio et al. [30] suggested the preexistence of drug resistant virus by calculating the drug resistant viral load before the initiation of drug therapy. Nowak et al. [13] compared the clinical data available on drug resistant virus development in patients with the results of the twin-strain mathematical model. Bonhoeffer and Nowak [31] showed the significance of the existence of drug resistant virus and discussed whether it exists before the onset of therapy or produced by replication of virus during treatment. Keeping in view the fact that the process of reverse transcription takes place in the early stage of infection before an infected T-cell produces virus particles. The classification of infected cells is done in two subclasses: pre-RT and post-RT classes. Srivastava et al. [32] proposed a twin-strain model using the above classification of infected cells to study the effect of RTI and PI drug on the emergence of drug resistant HIV virus. Further, a healthy immune response plays an important role in delaying the progression of HIV infection [33]. Kamboj and Sharma [34] discussed the importance of coupling between multidrug therapy and the immune response of a host person in HIV infection dynamics. Thus, in the present study, along with pre-RT and post-RT classification, full logistic term representing the proliferation of T-cells, the immune response of the body, and the administration of three drugs FI, RTI and PI are incorporated in the mathematical model. The above model is Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 2 (2024) https://internationalpubls.com 462 biologically more realistic, depicts a clearer picture of HIV infection, and has not been studied earlier in any available study by incorporating all the above factors using mathematical model. The motive of this study is to analyze the effectiveness of drugs FI, RTI and PI over the two strains of virus i.e., drug resistant and drug sensitive virus. Further, the impact of immune response on drug resistant virus strain is explored in the presence of drug therapies. The ultimate expectation is to find out any possibility of complete eradication of the virus in the presence of drug therapy and immune response. In the present study, a variable that represents cell population is considered as a continuous, differentiable variable, and the exact value of population is approximated through the nearest integer value of the corresponding variable. Considering all the above mentioned facts, a mathematical model is presented in section 2. In section 3, the model is analysed for non negative and bounded solution. In section 4 and 5, existence of steady states and their stabilities have been discussed. Section 6 has been devoted to the interpretation of all results with the help of numerical simulations. Finally, in section 7, the study is concluded by discussing various biological interpretations of obtained results. 2 Model formulation In the present model, a patient is considered under multidrug treatment with healthy cells 𝑇(𝑑) and infected cells 𝐼(𝑑), which are infected with free virus 𝑉(𝑑).The Fusion Inhibitor (FI) of efficacy 𝑓 ∈ [0,1) when applied prevent the entry of free virus into healthy cells. Kamboj and Sharma [11], modified the model of Srivastava et al. [9] by dividing the population of infected cells 𝐼(𝑑) into two categories: 𝑇1(𝑑) for pre-RT class and 𝑇2(𝑑) for post-RT class. The cells in pre-RT class (i.e., 𝑇1(𝑑)) will proceed to post-RT class to complete the HIV replication life cycle at a rate 𝛼 . But, on the application of the RTI drug therapy with efficacy πœ‚ ∈ (0,1), all the cells in pre-RT class will not be able to complete reverse transcription process. Therefore, a fraction (πœ‚π›Όπ‘‡1) of them will revert back to uninfected class and the remaining will proceed to post-RT class and turn into productively infected virus. The Protease Inhibitor (PI) drug with efficacy 𝛾 ∈ [0,1), prevents the post-RT cells to produce non-infectious virions with rate 𝛾𝑁 . Then, (1 βˆ’ 𝛾)𝑁 measures the concentration of infectious virions, where 𝑁 denotes the average number of viral particles produced by an infected cell. That means, the effect of PI restricts only to the infectious (𝑉) virions, which constitutes a part of the newly produced virions. Since the replication rate of HIV virus is exponentially high therefore, the process of reverse transcription of viral RNA to proviral DNA is highly error prone. Consequently, the probability of occurrence of mutations is very high. For example, the average number of changes per genome is 0.3 per replication cycle, i.e., after reverse transcription about 22 percent of infected cells should carry proviral genomes with one mutation [35]. Hence, in presence of multidrug therapy, two strains of the virus, i.e., drug sensitive strain and drug resistant strain, are to be incorporated in the model. Now, the infected cells in pre-RT class are to be divided into two categories, i.e., 𝑇1 = 𝑇1 𝑠 + 𝑇1 π‘Ÿ; infected either by drug sensitive virus (𝑇1 𝑠) or drug resistant virus (𝑇1 π‘Ÿ). Similarly, 𝑇2 𝑠 and 𝑇2 π‘Ÿ, the parts of 𝑇2 cells, are infected by drug sensitive and drug resistant virus respectively. 𝑉𝑠 and π‘‰π‘Ÿ be the population of infectious virus, which are drug sensitive and drug resistant respectively. To incorporate the response of the immune system, the CTLs/ immune cell population (𝐸) present in the body is to be included in the model. Since, after reverse transcription, CTLs attack only productively infected (post-RT) cells. It means, the other infected cells, which either revert back to uninfected class or in which reverse transcription has not been completed (i.e., pre-RT cells) do not have the ability to express HIV and cannot invite CTLs for immunity support. Therefore, the intensity of the immune response should depend on the concentration of post-RT cells (𝑇2 𝑠 and 𝑇2 π‘Ÿ ). The mathematical model representing the above dynamics is written as follows: 𝑑𝑇 𝑑𝑑 = 𝑠 βˆ’ πœ‡π‘‡ + π‘Ÿπ‘‡ (1 βˆ’ 𝑇 π‘‡π‘šπ‘Žπ‘₯ ) βˆ’ (1 βˆ’ 𝑓𝑠)π‘˜π‘‰π‘ π‘‡ βˆ’ (1 βˆ’ π‘“π‘Ÿ)π‘˜π‘‰π‘Ÿπ‘‡ + Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 2 (2024) https://internationalpubls.com 463 (𝑏𝑠 + πœ‚ 𝑠𝛼𝑠)𝑇1 𝑠 + (π‘π‘Ÿ + πœ‚ π‘Ÿπ›Όπ‘Ÿ)𝑇1 π‘Ÿ, (2.1) 𝑑𝑇1 𝑠 𝑑𝑑 = (1 βˆ’ 𝑓𝑠)π‘˜π‘‰π‘ π‘‡ βˆ’ (πœ‡1 + 𝛼𝑠 + 𝑏𝑠)𝑇1 𝑠, (2.2) 𝑑𝑇2 𝑠 𝑑𝑑 = (1 βˆ’ πœ‡π‘š)(1 βˆ’ πœ‚ 𝑠)𝛼𝑠𝑇1 𝑠 βˆ’ 𝛿𝑠𝑇2 𝑠 βˆ’ 𝑑π‘₯𝐸𝑇2 𝑠, (2.3) 𝑑𝑉𝑠 𝑑𝑑 = 𝑁𝛿𝑠(1 βˆ’ 𝛾 𝑠)𝑇2 𝑠 βˆ’ πœ‡π‘£π‘‰π‘ , (2.4) 𝑑𝑇1 π‘Ÿ 𝑑𝑑 = (1 βˆ’ π‘“π‘Ÿ)π‘˜π‘‡π‘‰π‘Ÿ βˆ’ (πœ‡1 + π›Όπ‘Ÿ + π‘π‘Ÿ)𝑇1 π‘Ÿ , (2.5) 𝑑𝑇2 π‘Ÿ 𝑑𝑑 = πœ‡π‘š(1 βˆ’ πœ‚ 𝑠)𝛼𝑠𝑇1 𝑠 + π›Όπ‘Ÿ(1 βˆ’ πœ‚ π‘Ÿ)𝑇1 π‘Ÿ βˆ’ π›Ώπ‘Ÿπ‘‡2 π‘Ÿ βˆ’ 𝑑π‘₯𝐸𝑇2 π‘Ÿ, (2.6) π‘‘π‘‰π‘Ÿ 𝑑𝑑 = π‘π›Ώπ‘Ÿ(1 βˆ’ 𝛾 π‘Ÿ)𝑇2 π‘Ÿ βˆ’ πœ‡π‘£π‘‰π‘Ÿ, (2.7) 𝑑𝐸 𝑑𝑑 = 𝑝(𝑇2 𝑠 + 𝑇2 π‘Ÿ) βˆ’ 𝑑𝐸𝐸, (2.8) with 𝑇(0) = 𝑇00, 𝐸(0) = 𝐸0, 𝑇1 𝑠(0) = 𝑇10, 𝑇2 𝑠(0) = 𝑇20, 𝑉𝑠(0) = 𝑉10, 𝑇1 π‘Ÿ(0) = 𝑇11, 𝑇2 π‘Ÿ(0) = 𝑇21, π‘‰π‘Ÿ(0) = 𝑉20. In equation (2.1), 𝑠 represents the rate at which new T-cells are created from sources within the body, such as thymus. The natural decay of these cells with time is given by πœ‡π‘‡. The logistic expression π‘Ÿπ‘‡(1 βˆ’ 𝑇 π‘‡π‘šπ‘Žπ‘₯ ) represents the T-cells created by the proliferation of existing T-cells, in the presence of infection. Detailed discussion on the role of logistic term is found in [4, 36-37]. The parameter π‘˜ represents the interaction-infection rate of T-cells with the virus, assumed to be same for both strains and πœ‡1 is the death rate of infected cells in pre-RT class. 𝑏𝑠 and π‘π‘Ÿ are the reverting rates of infected T-cells to uninfected class due to the non-completion of reverse transcription for respective strains whereas 𝛼𝑠 and π›Όπ‘Ÿ denotes the rate of transition of T cells from pre-RT class to post-RT class. In equation (2.3) and (2.6) 𝛿𝑠 and π›Ώπ‘Ÿ denote the death rate of actively infected cells in post-RT class for respective strains. In equations (2.4) and (2.7) πœ‡π‘£ denotes the clearance rate of virus which is assumed to be the same for both strains. 𝑁 in equations (2.4) and (2.7) represents the average number of viral particles produced by an infected cell, assumed to be the same for both strains. The parameters 𝑓𝑠, πœ‚π‘  and 𝛾𝑠 and π‘“π‘Ÿ, πœ‚π‘Ÿ and π›Ύπ‘Ÿ in [0,1) represent the efficacy of drug FI, RTI and PI corresponding to drug sensitive and drug resistant virus strains, respectively. The parameter πœ‡π‘š in equations (2.3) and (2.6) represent the rate at which cells infected by drug sensitive virus mutate and become drug resistant virus during the process of reverse transcription. The backward mutation from drug resistant to drug sensitive strain has not been considered in the present study. The parameter 𝑑π‘₯ denotes the rate of clearance of infected cells (𝑇2 𝑠 and 𝑇2 π‘Ÿ) by CTLs. Therefore, term 𝑑π‘₯𝐸𝑇2 π‘Ÿ and 𝑑π‘₯𝐸𝑇2 𝑠 in equations (2.3) and (2.6) represent the loss of infected cells (𝑇2 𝑠 and 𝑇2 π‘Ÿ) by CTLs. For simplicity, it is assumed that the immune cell population (𝐸) are produced at the same constant proliferation rate 𝑝 and the death rate 𝑑𝐸 whether they are produced as a result of the presence of either kind of productively infected cells (𝑇2 𝑠 and 𝑇2 π‘Ÿ). 3 Analysis The variables of the mathematical model (2.1-2.8) considered in the previous section represent the populations and for the model to be biologically realistic, it does not allow the cell populations to grow unbounded or get a negative value for all time. For the positivity of the solutions of the model, a non-negative orthant, 𝑅+ 8 = {π‘₯ ∈ 𝑅8|π‘₯ β‰₯ 0}, is defined to contain, forever, any trajectory that starts in it. For this model, we have 𝑑𝑇 𝑑𝑑 |𝑇=0 = 𝑠 + (𝑏𝑠 + πœ‚ 𝑠𝛼𝑠)𝑇1 𝑠 + (π‘π‘Ÿ + πœ‚ π‘Ÿπ›Όπ‘Ÿ)𝑇1 π‘Ÿ β‰₯ 0, 𝑑𝑇1 𝑠 𝑑𝑑 |𝑇1𝑠=0 = (1 βˆ’ 𝑓𝑠)π‘˜π‘‰π‘ π‘‡ β‰₯ 0, Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 2 (2024) https://internationalpubls.com 464 𝑑𝑇2 𝑠 𝑑𝑑 |𝑇2𝑠=0 = (1 βˆ’ πœ‡π‘š)(1 βˆ’ πœ‚ 𝑠)𝛼𝑠𝑇1 𝑠 β‰₯ 0, 𝑑𝑉𝑠 𝑑𝑑 |𝑉𝑠=0 = 𝑁(1 βˆ’ 𝛾 𝑠)𝛿𝑠𝑇2 𝑠 β‰₯ 0, 𝑑𝑇1 π‘Ÿ 𝑑𝑑 |𝑇1π‘Ÿ=0 = (1 βˆ’ π‘“π‘Ÿ)π‘˜π‘‰π‘Ÿπ‘‡ β‰₯ 0, 𝑑𝑇2 π‘Ÿ 𝑑𝑑 |𝑇2π‘Ÿ=0 = πœ‡π‘š(1 βˆ’ πœ‚ 𝑠)𝛼𝑠𝑇1 𝑠 + π›Όπ‘Ÿ(1 βˆ’ πœ‚ π‘Ÿ)𝑇1 π‘Ÿ β‰₯ 0, π‘‘π‘‰π‘Ÿ 𝑑𝑑 |π‘‰π‘Ÿ=0 = 𝑁(1 βˆ’ 𝛾 π‘Ÿ)π›Ώπ‘Ÿπ‘‡2 π‘Ÿ β‰₯ 0, 𝑑𝐸 𝑑𝑑 |𝐸=0 = 𝑝(𝑇2 𝑠 + 𝑇2 π‘Ÿ) β‰₯ 0. This shows that the vector field (𝑇, 𝑇1 𝑠, 𝑇2 𝑠, 𝑉𝑠, 𝑇1 π‘Ÿ, 𝑇2 π‘Ÿ , π‘‰π‘Ÿ, 𝐸), on each bounding hyperplane of 𝑅+ 8 , is pointing to the inward direction of 𝑅+ 8 . That means, all the solution trajectories initiating in 𝑅+ 8 , will remain inside 𝑅+ 8 for all 𝑑. Hence, the positivity of the solutions initiating in the interior of 𝑅+ 8 is guaranteed. Further on adding the equations (2.1), (2.2), (2.3), (2.5) and (2.6), we have, 𝑑 𝑑𝑑 (𝑇 + 𝑇1 𝑠 + 𝑇2 𝑠 + 𝑇1 π‘Ÿ + 𝑇2 π‘Ÿ) = 𝑠 βˆ’ πœ‡π‘‡ + π‘Ÿπ‘‡ (1 βˆ’ 𝑇 π‘‡π‘šπ‘Žπ‘₯ ) βˆ’ πœ‡1𝑇1 𝑠 βˆ’ πœ‡1𝑇1 π‘Ÿ βˆ’ 𝛿𝑠𝑇2 𝑠 βˆ’ π›Ώπ‘Ÿπ‘‡2 π‘Ÿ βˆ’ 𝑑π‘₯𝐸(𝑇2 𝑠 + 𝑇2 π‘Ÿ), ≀ 𝑠 + π‘Ÿπ‘‡(1 βˆ’ 𝑇 π‘‡π‘šπ‘Žπ‘₯ ) βˆ’ πœ‡(𝑇 + 𝑇1 𝑠 + 𝑇2 𝑠 + 𝑇1 π‘Ÿ + 𝑇2 π‘Ÿ), (since 𝛿𝑠 > π›Ώπ‘Ÿ > πœ‡1 > πœ‡). Let us denote 𝐢 = max {𝑠 + π‘Ÿπ‘‡(1 βˆ’ 𝑇 π‘‡π‘šπ‘Žπ‘₯ )} , for 𝑇 ∈ (0, 𝑇0] , where 𝑇0 = (π‘Ÿβˆ’πœ‡)+√(π‘Ÿβˆ’πœ‡)2+4π‘Ÿ1𝑠 2π‘Ÿ1 with π‘Ÿ1 = π‘Ÿ π‘‡π‘šπ‘Žπ‘₯ , obtained in next section from equation 2.1 so that π‘‘β†’βˆž lim sup(𝑇 + 𝑇1 𝑠 + 𝑇2 𝑠 + 𝑇1 π‘Ÿ + 𝑇2 π‘Ÿ) ≀ 𝐢 πœ‡ . Therefore, without any loss of generality, it can be assumed that π‘‘β†’βˆž lim sup𝑇(𝑑) ≀ 𝐢 πœ‡ , π‘‘β†’βˆž lim sup𝑇1 𝑠(𝑑) ≀ 𝐢 πœ‡ π‘‘β†’βˆž lim sup𝑇2 𝑠(𝑑) ≀ 𝐢 πœ‡ , π‘‘β†’βˆž lim sup𝑇1 π‘Ÿ(𝑑) ≀ 𝐢 πœ‡ , π‘‘β†’βˆž lim sup𝑇2 π‘Ÿ(𝑑) ≀ 𝐢 πœ‡ . Now, this bound for 𝑇2 𝑠 and 𝑇2 π‘Ÿ enables to find the bounds for 𝑉𝑠(𝑑) and π‘‰π‘Ÿ(𝑑) and 𝐸(𝑑) from the equations (2.8) and (2.9), respectively. So, finally, we have a bounded set 𝑆 = {(𝑇, 𝑇1 𝑠, 𝑇2 𝑠, 𝑉𝑠, 𝑇1 π‘Ÿ , 𝑇2 π‘Ÿ, π‘‰π‘Ÿ, 𝐸) ∈ 𝑅+ 8 ; 0 ≀ 𝑇, 𝑇1 𝑠, 𝑇2 𝑠, 𝑇1 π‘Ÿ, 𝑇2 π‘Ÿ ≀ 𝐢 πœ‡ , 0 ≀ 𝑉𝑠 ≀ 𝑁(1βˆ’π›Ύπ‘ )𝛿𝑠𝐢 πœ‡π‘£πœ‡ , 0 ≀ π‘‰π‘Ÿ ≀ 𝑁(1βˆ’π›Ύπ‘Ÿ)π›Ώπ‘ŸπΆ πœ‡π‘£πœ‡ , 0 ≀ 𝐸 ≀ 2𝐢𝑝 π‘‘πΈπœ‡ }. Then, any solution trajectory, which initiates from an interior point of 𝑅+ 8 , enters 𝑆 and remains there forever. 4 Steady states The model system (2.1)-(2.8) has three steady states: (a) The infection free steady state 𝐸0 = (𝑇0, 0,0,0,0,0,0,0), where, for a new parameter π‘Ÿ1 = π‘Ÿ π‘‡π‘šπ‘Žπ‘₯ , the equation (2.1) is solved to get 𝑇0 = (π‘Ÿβˆ’πœ‡)+√(π‘Ÿβˆ’πœ‡)2+4π‘Ÿ1𝑠 2π‘Ÿ1 . (b) The infected steady state πΈπ‘Ÿ = (𝑇, 0,0,0, 𝑇1 π‘Ÿ , 𝑇2 π‘Ÿ, π‘‰π‘Ÿ, 𝐸), with only drug resistant viral strain, where 𝑇 = (π‘Ÿβˆ’πœ‡+𝛽2)+√(π‘Ÿβˆ’πœ‡+𝛽2) 2+4(𝛽1+π‘Ÿ1)𝑠 2(𝛽1+π‘Ÿ1) , 𝑇1 π‘Ÿ = (1βˆ’π‘“π‘Ÿ)π‘˜π‘‡π‘‰π‘Ÿ πœ‡1+π›Όπ‘Ÿ+π›½π‘Ÿ ,𝑇2 π‘Ÿ = πœ‡π‘£π‘‰π‘Ÿ π‘π›Ώπ‘Ÿ(1βˆ’π›Ύ π‘Ÿ) , π‘‰π‘Ÿ = π‘˜1𝑇 βˆ’ π‘˜2 , 𝐸 = π‘πœ‡π‘£π‘‰π‘Ÿ π‘‘πΈπ‘π›Ώπ‘Ÿ(1βˆ’π›Ύ π‘Ÿ) , π‘˜1 = 𝑁2π›Όπ‘Ÿπ›Ώπ‘Ÿ 2π‘‘πΈπ‘˜(1βˆ’πœ‚ π‘Ÿ)(1βˆ’π‘“π‘Ÿ)(1βˆ’π›Ύπ‘Ÿ)2 𝑑π‘₯π‘πœ‡π‘£ 2(πœ‡1+π›Όπ‘Ÿ+π›½π‘Ÿ) , π‘˜2 = π‘π‘‘πΈπ›Ώπ‘Ÿ 2(1βˆ’π›Ύπ‘Ÿ) 𝑑π‘₯π‘πœ‡π‘£ , 𝛽1 = (1βˆ’π‘“π‘Ÿ)((1βˆ’πœ‚π‘Ÿ)π›Όπ‘Ÿ+πœ‡1)π‘˜π‘˜1 πœ‡1+π›Όπ‘Ÿ+π›½π‘Ÿ , 𝛽2 = π‘˜2𝛽1 π‘˜1 . (c) The infected steady state with both drug sensitive and drug resistant strains present, is given by πΈπ‘š = (οΏ½ΜƒοΏ½, οΏ½ΜƒοΏ½1 𝑠, οΏ½ΜƒοΏ½2 𝑠, �̃�𝑠, οΏ½ΜƒοΏ½1 π‘Ÿ , οΏ½ΜƒοΏ½2 π‘Ÿ, οΏ½ΜƒοΏ½π‘Ÿ, οΏ½ΜƒοΏ½), where �̃�𝑠 = π‘˜7οΏ½ΜƒοΏ½ βˆ’ π‘˜8 βˆ’ π‘˜9οΏ½ΜƒοΏ½π‘Ÿ, οΏ½ΜƒοΏ½π‘Ÿ = 𝛼1οΏ½ΜƒοΏ½ 2βˆ’π›½3𝑇 𝛾1𝑇+𝛿1 , οΏ½ΜƒοΏ½1 𝑠 = (1βˆ’π‘“π‘ )π‘˜οΏ½ΜƒοΏ½π‘‰π‘  πœ‡1+𝛼𝑠+𝑏𝑠 , οΏ½ΜƒοΏ½2 𝑠 = πœ‡π‘£π‘‰π‘  𝑁𝛿𝑠(1βˆ’π›Ύ 𝑠) , οΏ½ΜƒοΏ½1 π‘Ÿ = (1βˆ’π‘“π‘Ÿ)π‘˜οΏ½ΜƒοΏ½π‘‰π‘Ÿ πœ‡1+π›Όπ‘Ÿ+π‘π‘Ÿ , οΏ½ΜƒοΏ½2 π‘Ÿ = πœ‡π‘£π‘‰π‘Ÿ π‘π›Ώπ‘Ÿ(1βˆ’π›Ύ π‘Ÿ) , οΏ½ΜƒοΏ½ = 𝑝(οΏ½ΜƒοΏ½2 𝑠+οΏ½ΜƒοΏ½2 π‘Ÿ) 𝑑𝐸 , Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 2 (2024) https://internationalpubls.com 465 π‘˜7 = π‘˜π‘‘πΈπ‘ 2𝛿𝑠 2(1βˆ’π›Ύπ‘ )2(1βˆ’πœ‚π‘ )(1βˆ’πœ‡π‘š)(1βˆ’π‘“ 𝑠)𝛼𝑠 𝑑π‘₯π‘πœ‡π‘£ 2(πœ‡1+𝛼𝑠+𝑏𝑠) , π‘˜8 = 𝑑𝐸𝑁𝛿𝑠 2(1βˆ’π›Ύπ‘ ) 𝑑π‘₯π‘πœ‡π‘£ , π‘˜9 = 𝛿𝑠(1βˆ’π›Ύ 𝑠) π›Ώπ‘Ÿ(1βˆ’π›Ύ π‘Ÿ) , 𝛼1 = π‘˜7π‘˜10, 𝛽3 = π‘˜8π‘˜10, 𝛾1 = π‘˜9π‘˜10 + π‘˜7π‘˜14 βˆ’ π‘˜11, 𝛿1 = π‘˜12 βˆ’ π‘˜14π‘˜8. οΏ½ΜƒοΏ½ is obtained as the root of the cubic equation, given by 𝛼2οΏ½ΜƒοΏ½ 3 + 𝛼3οΏ½ΜƒοΏ½ 2 + 𝛼4οΏ½ΜƒοΏ½ + 𝛼5 = 0, where 𝛼2 = ((1 βˆ’ 𝑓 𝑠)π‘˜π‘˜9 βˆ’ (1 βˆ’ 𝑓 π‘Ÿ)π‘˜ βˆ’ (1βˆ’π‘“π‘ )π‘˜π‘˜9(𝑏𝑠+πœ‚ 𝑠𝛼𝑠) πœ‡1+𝛼𝑠+𝑏𝑠 + (1βˆ’π‘“π‘Ÿ)(π‘π‘Ÿ+πœ‚ π‘Ÿπ›Όπ‘Ÿ)π‘˜ πœ‡1+π›Όπ‘Ÿ+π‘π‘Ÿ )𝛼1 + ( (1βˆ’π‘“π‘ )(𝑏𝑠+πœ‚ 𝑠𝛼𝑠)π‘˜π‘˜7 πœ‡1+𝛼𝑠+𝑏𝑠 βˆ’ (1 βˆ’ 𝑓𝑠)π‘˜π‘˜7 βˆ’ π‘Ÿ1)𝛾1, 𝛼3 = 𝛾1(βˆ’πœ‡ + (1 βˆ’ 𝑓 𝑠)π‘˜π‘˜8 βˆ’ (1βˆ’π‘“π‘ )(𝑏𝑠+πœ‚ 𝑠𝛼𝑠)π‘˜π‘˜8 πœ‡1+𝛼𝑠+𝑏𝑠 + π‘Ÿ) +𝛿1(βˆ’(1 βˆ’ 𝑓 𝑠)π‘˜π‘˜7 + (1βˆ’π‘“π‘ )(𝑏𝑠+πœ‚ 𝑠𝛼𝑠)π‘˜π‘˜7 πœ‡1+𝛼𝑠+𝑏𝑠 βˆ’ π‘Ÿ1) + 𝛽1((1 βˆ’ π‘“π‘Ÿ)π‘˜ βˆ’ (1 βˆ’ 𝑓𝑠)π‘˜π‘˜9 + (1βˆ’π‘“π‘ )(𝑏𝑠+πœ‚ 𝑠𝛼𝑠)π‘˜π‘˜9 πœ‡1+𝛼𝑠+𝑏𝑠 βˆ’ (1βˆ’π‘“π‘Ÿ)(π‘π‘Ÿ+πœ‚ π‘Ÿπ›Όπ‘Ÿ)π‘˜ πœ‡1+π›Όπ‘Ÿ+π‘π‘Ÿ ), 𝛼4 = 𝑠𝛾1 + ((1 βˆ’ 𝑓 𝑠)π‘˜π‘˜8 βˆ’ πœ‡ βˆ’ (1βˆ’π‘“π‘ )(𝑏𝑠+πœ‚ 𝑠𝛼𝑠) πœ‡1+𝛼𝑠+𝑏𝑠 π‘˜π‘˜8 + π‘Ÿ)𝛿1, 𝛼5 = 𝑠𝛿1. Further, it is noted that π‘‰π‘Ÿ = 𝑇1 π‘Ÿ = 𝑇2 π‘Ÿ = 0 if πœ‡π‘š = 0. Thus, the steady state πΈπ‘š reduces to steady state with sensitive virus only, (say, 𝐸𝑠). 5 Stability of steady states The asymptotic stability of a steady state is decided by the eigenvalues of the Jacobian matrix. In the present problem, the system (2.1)-(2.8) is linearized around a steady state, and the corresponding Jacobian matrix 𝑱 is obtained as follows: 𝑱 = ( βˆ’π‘€ 𝑏𝑠 + πœ‚ 𝑠𝛼𝑠 0 βˆ’(1 βˆ’ 𝑓𝑠)π‘˜π‘‡ π‘π‘Ÿ + πœ‚ π‘Ÿπ›Όπ‘Ÿ 0 βˆ’(1 βˆ’ π‘“π‘Ÿ)π‘˜π‘‡ 0 (1 βˆ’ 𝑓𝑠)π‘˜π‘‰π‘  βˆ’(πœ‡1+ 𝛼𝑠 + 𝑏𝑠) 0 (1 βˆ’ 𝑓𝑠)π‘˜π‘‡ 0 0 0 0 0 (1 βˆ’ πœ‚π‘ )(1 βˆ’ πœ‡π‘š)𝛼𝑠 βˆ’(𝛿𝑠 + 𝑑π‘₯𝐸) 0 0 0 0 βˆ’π‘‘π‘₯𝑇2 𝑠 0 0 𝑁(1 βˆ’ 𝛾𝑠)𝛿𝑠 βˆ’πœ‡π‘£ 0 0 0 0 (1 βˆ’ π‘“π‘Ÿ)π‘˜π‘‰π‘Ÿ 0 0 0 βˆ’(πœ‡1+ π›Όπ‘Ÿ + π‘π‘Ÿ) 0 (1 βˆ’ π‘“π‘Ÿ)π‘˜π‘‡ 0 0 πœ‡π‘š(1 βˆ’ πœ‚ 𝑠)𝛼𝑠 0 0 π›Όπ‘Ÿ(1 βˆ’ πœ‚ π‘Ÿ) βˆ’(π›Ώπ‘Ÿ + 𝑑π‘₯𝐸) 0 βˆ’π‘‘π‘₯𝑇2 π‘Ÿ 0 0 0 0 0 𝑁(1 βˆ’ π›Ύπ‘Ÿ)π›Ώπ‘Ÿ βˆ’πœ‡π‘£ 0 0 0 𝑝 0 0 𝑝 0 βˆ’π‘‘πΈ ) where 𝑀 = πœ‡ βˆ’ π‘Ÿ + 2π‘Ÿ1𝑇 + (1 βˆ’ 𝑓 π‘Ÿ)π‘˜π‘‰π‘Ÿ + (1 βˆ’ 𝑓 𝑠)π‘˜π‘‰π‘  is a positive value. 5.1 Stability of uninfected steady state π‘¬πŸŽ At steady state 𝐸0, the corresponding Jacobian matrix (say, π‘±πŸŽ) is obtained by substituting 𝑇 = 𝑇0 and 𝑇1 𝑠 = 𝑇2 𝑠 = 𝑉𝑠 = 𝑇1 π‘Ÿ = 𝑇2 π‘Ÿ = π‘‰π‘Ÿ = 𝐸 = 0 in the Jacobian matrix 𝑱. π‘±πŸŽ = Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 2 (2024) https://internationalpubls.com 466 ( βˆ’π‘€0 𝑏𝑠 + πœ‚ 𝑠𝛼𝑠 0 βˆ’(1 βˆ’ 𝑓𝑠)π‘˜π‘‡0 π‘π‘Ÿ + πœ‚ π‘Ÿπ›Όπ‘Ÿ 0 βˆ’(1 βˆ’ π‘“π‘Ÿ)π‘˜π‘‡0 0 0 βˆ’(πœ‡1 + 𝛼𝑠 + 𝑏𝑠) 0 (1 βˆ’ 𝑓𝑠)π‘˜π‘‡0 0 0 0 0 0 (1 βˆ’ πœ‚π‘ )(1 βˆ’ πœ‡π‘š)𝛼𝑠 βˆ’π›Ώπ‘  0 0 0 0 βˆ’π‘‘π‘₯𝑇2 𝑠 0 0 𝑁(1 βˆ’ 𝛾𝑠)𝛿𝑠 βˆ’πœ‡π‘£ 0 0 0 0 0 0 0 0 βˆ’(πœ‡1 + π›Όπ‘Ÿ + π‘π‘Ÿ) 0 (1 βˆ’ π‘“π‘Ÿ)π‘˜π‘‡0 0 0 πœ‡π‘š(1 βˆ’ πœ‚ 𝑠)𝛼𝑠 0 0 π›Όπ‘Ÿ(1 βˆ’ πœ‚ π‘Ÿ) βˆ’(π›Ώπ‘Ÿ + 𝑑π‘₯𝐸) 0 βˆ’π‘‘π‘₯𝑇2 π‘Ÿ 0 0 0 0 0 𝑁(1 βˆ’ π›Ύπ‘Ÿ)π›Ώπ‘Ÿ βˆ’πœ‡π‘£ 0 0 0 𝑝 0 0 𝑝 0 βˆ’π‘‘πΈ ) where 𝑀0 = πœ‡ βˆ’ π‘Ÿ + 2π‘Ÿ1𝑇0 is a positive value. The characteristic equation |π‘±πŸŽ βˆ’ πœ†π‘°| = 0 corresponding to the Jacobian matrix π‘±πŸŽ is given by (πœ† + 𝑀0)(πœ† + 𝑑𝐸)𝑃0(πœ†)𝑄0(πœ†) = 0; 𝑀0 = πœ‡ βˆ’ π‘Ÿ + 2π‘Ÿ1𝑇0 > 0, (5.1) where 𝑃0(πœ†) = πœ† 3 + 𝐴0πœ† 2 + 𝐡0πœ† + 𝐢0 and 𝑄0(πœ†) = πœ† 3 + 𝐴1πœ† 2 + 𝐡1πœ† + 𝐢1. The coefficients of these polynomials are expressed as follows: 𝐴0 = πœ‡1 + 𝛼𝑠 + 𝑏𝑠 + 𝛿𝑠 + πœ‡π‘£ , 𝐡0 = (πœ‡1 + 𝛼𝑠 + 𝑏𝑠)(𝛿𝑠 + πœ‡π‘£) + π›Ώπ‘ πœ‡π‘£ , 𝐢0 = π‘˜π›Όπ‘ π›Ώπ‘ (1 βˆ’ πœ‚ 𝑠)(1 βˆ’ 𝛾𝑠)(1 βˆ’ 𝑓𝑠)(1 βˆ’ πœ‡π‘š)(𝑁01 βˆ’π‘)𝑇0; 𝑁01 = (πœ‡1+𝛼𝑠+𝑏𝑠)πœ‡π‘£ π‘˜π‘‡0𝛼𝑠(1βˆ’πœ‚ 𝑠)(1βˆ’π›Ύπ‘ )(1βˆ’π‘“π‘ )(1βˆ’πœ‡π‘š) , 𝐴1 = πœ‡1 + π›Όπ‘Ÿ + π‘π‘Ÿ + π›Ώπ‘Ÿ + πœ‡π‘£ , 𝐡1 = (πœ‡1 + π›Όπ‘Ÿ + π‘π‘Ÿ)(π›Ώπ‘Ÿ + πœ‡π‘£) + π›Ώπ‘Ÿπœ‡π‘£ , 𝐢1 = π‘˜π›Όπ‘Ÿπ›Ώπ‘Ÿ(1 βˆ’ πœ‚ π‘Ÿ)(1 βˆ’ π›Ύπ‘Ÿ)(1 βˆ’ π‘“π‘Ÿ)(𝑁02 βˆ’π‘); 𝑁02 = (πœ‡1+π›Όπ‘Ÿ+π‘π‘Ÿ)πœ‡π‘£ π‘˜π‘‡0π›Όπ‘Ÿ(1βˆ’πœ‚ π‘Ÿ)(1βˆ’π‘“π‘Ÿ)(1βˆ’π›Ύπ‘Ÿ) . The characteristic equation (5.1) provides πœ† = βˆ’π‘€0 and βˆ’π‘‘πΈ as two eigenvalues of the Jacobian matrix 𝐽0. The remaining six eigenvalues are obtained from the roots of 𝑃0(πœ†) = 0 and 𝑄0(πœ†) = 0. The stability of uninfected steady state 𝐸0 is ensured through the negative real parts of all the eight eigenvalues of 𝐽0. Obviously, the eigenvalues βˆ’π‘€0 and βˆ’π‘‘πΈ meet this requirement. But, for other six eigenvalues, the roots of 𝑃0(πœ†) = 0 and 𝑄0(πœ†) = 0 are to be checked. According to Routh-Hurwitz criterion [38], all the roots of 𝑃0(πœ†) = 0 and 𝑄0(πœ†) = 0 will have negative real parts if and only if 𝐴0, 𝐡0, 𝐢0, 𝐴1, 𝐡1, 𝐢1 are all positive and 𝐴0𝐡0 > 𝐢0, 𝐴1𝐡1 > 𝐢1. It is noted that 𝐴0, 𝐡0, 𝐴0𝐡0 βˆ’ 𝐢0, 𝐴1, 𝐡1, 𝐴1𝐡1 βˆ’ 𝐢1 are all positive, and therefore, the onus of deciding the asymptotic stability of 𝐸0 stays with the value of 𝐢0 and 𝐢1 only. The coefficients 𝐢0 and 𝐢1 are positive if 𝑁01 > 𝑁 and 𝑁02 > 𝑁 respectively. That means, the asymptotically stability of the uninfected state 𝐸0 is ensured with 𝑁 < 𝑁01 or 𝑁 < 𝑁02. On the other hand, for 𝑁 > 𝑁01 or 𝑁 > 𝑁02, 𝐢0 and 𝐢1 become negative, which implies a sign change in the coefficients of the cubic equations 𝑃0(πœ†) = 0 and 𝑄0(πœ†) = 0. Then, according to the Descartes’ rule of signs, one positive root of the equation implies a positive eigenvalue for 𝐽0. That means, the uninfected state 𝐸0 cannot be stable for 𝑁 > 𝑁01 and 𝑁 > 𝑁02. Also for 𝑁 = 𝑁01 or 𝑁 = 𝑁02, the cubic equation 𝑃0(πœ†) = 0 or 𝑄0(πœ†) = 0 yields a zero eigenvalue and the reduced quadratic equation will have roots with negative real parts. Thus, according to Routh-Hurwitz conditions, the state 𝐸0 becomes neutrally stable, when 𝑁 = 𝑁01 or 𝑁 = 𝑁02, . Proposition 1. The uninfected steady state 𝐸0 is locally asymptotically stable if 𝑁01 and 𝑁02 are greater than 𝑁. 5.2 Stability of steady state (𝑬𝒓) infected with only drug resistant viral strain The Jacobian matrix 𝑱𝒓 evaluated at steady state πΈπ‘Ÿ is obtained from the Jacobian matrix 𝑱 , by substituting 𝑇 = 𝑇, 𝑇1 𝑠 = 𝑇2 𝑠 = 𝑉𝑠 = 0, 𝑇1 π‘Ÿ = 𝑇1 π‘Ÿ , 𝑇2 π‘Ÿ = 𝑇2 π‘Ÿ , π‘‰π‘Ÿ = π‘‰π‘Ÿ and 𝐸 = 𝐸. Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 2 (2024) https://internationalpubls.com 467 𝑱𝒓 = ( βˆ’π‘€1 𝑏𝑠 + πœ‚ 𝑠𝛼𝑠 0 βˆ’(1 βˆ’ 𝑓𝑠)π‘˜π‘‡ π‘π‘Ÿ + πœ‚ π‘Ÿπ›Όπ‘Ÿ 0 βˆ’(1 βˆ’ π‘“π‘Ÿ)π‘˜π‘‡ 0 π‘œ βˆ’(πœ‡1 + 𝛼𝑠 + 𝑏𝑠) 0 (1 βˆ’ 𝑓𝑠)π‘˜π‘‡ 0 0 0 0 0 (1 βˆ’ πœ‚π‘ )(1 βˆ’ πœ‡π‘š)𝛼𝑠 βˆ’(𝛿𝑠 + 𝑑π‘₯𝐸) 0 0 0 0 βˆ’π‘‘π‘₯𝑇2 𝑠 0 0 𝑁(1 βˆ’ 𝛾𝑠)𝛿𝑠 βˆ’πœ‡π‘£ 0 0 0 0 (1 βˆ’ π‘“π‘Ÿ)π‘˜π‘‰π‘Ÿ 0 0 0 βˆ’(πœ‡1 + π›Όπ‘Ÿ + π‘π‘Ÿ) 0 (1 βˆ’ π‘“π‘Ÿ)π‘˜π‘‡ 0 0 πœ‡π‘š(1 βˆ’ πœ‚ 𝑠)𝛼𝑠 0 0 π›Όπ‘Ÿ(1 βˆ’ πœ‚ π‘Ÿ) βˆ’(π›Ώπ‘Ÿ + 𝑑π‘₯𝐸) 0 βˆ’π‘‘π‘₯𝑇2 π‘Ÿ 0 0 0 0 0 𝑁(1 βˆ’ π›Ύπ‘Ÿ)π›Ώπ‘Ÿ βˆ’πœ‡π‘£ 0 0 0 𝑝 0 0 𝑝 0 βˆ’π‘‘πΈ ) where 𝑀1 = πœ‡ βˆ’ π‘Ÿ + 2π‘Ÿ1𝑇 + (1 βˆ’ 𝑓 π‘Ÿ)π‘˜π‘‰π‘Ÿ is a positive value. The corresponding characteristic equation |𝑱𝒓 βˆ’ πœ†π‘°| = 0 is expressed as 𝑃1(πœ†)𝑄1(πœ†) = 0, (5.2) where 𝑃1(πœ†) = πœ† 5 + 𝐴3πœ† 4 + 𝐡3πœ† 3 + 𝐢3πœ† 2 + 𝐷3πœ† + 𝐸3 and 𝑄1(πœ†) = πœ† 3 + 𝐴4πœ† 2 + 𝐡4πœ† + 𝐢4. The coefficients of these polynomials are expressed as follows: 𝐴3 = π‘˜4 + π‘˜5 + π‘˜6 +𝑀1; 𝐡3 = πœ‡π‘£π‘‘πΈ + (π‘˜4 + π‘˜5)π‘˜6 + π‘˜4π‘˜5 + 𝑝𝑑π‘₯𝑇2 π‘Ÿ + (π‘˜4 + π‘˜5 + π‘˜6)𝑀1 βˆ’ (π‘π‘Ÿ + πœ‚ π‘Ÿπ›Όπ‘Ÿ)(1 βˆ’ π‘“π‘Ÿ)π‘˜π‘‰π‘Ÿ, 𝐢3 = (π‘˜4 + π‘˜5)πœ‡π‘£π‘‘πΈ + π‘˜4π‘˜5π‘˜6 + 𝑑π‘₯𝑝(πœ‡π‘£ + π‘˜4)𝑇2 π‘Ÿ βˆ’ π‘˜π›Όπ‘Ÿπ›Ώπ‘Ÿ(1 βˆ’ πœ‚ π‘Ÿ)(1 βˆ’ 𝑓𝑠)(1 βˆ’ π›Ύπ‘Ÿ)𝑁𝑇 + (πœ‡π‘£π‘‘πΈ + (π‘˜4 + π‘˜5)π‘˜6 + π‘˜4π‘˜5 + 𝑝𝑑π‘₯𝑇2 π‘Ÿ )𝑀1 βˆ’ π‘˜(1 βˆ’ 𝑓 π‘Ÿ)(π‘π‘Ÿ + πœ‚π‘Ÿπ›Όπ‘Ÿ)(π‘˜5π‘˜6)π‘‰π‘Ÿ, 𝐷3 = πœ‡π‘£π‘‘πΈπ‘˜4π‘˜5 + 𝑑π‘₯π‘π‘˜4πœ‡π‘£π‘‡2 π‘Ÿ βˆ’ π‘˜π›Όπ‘Ÿπ›Ώπ‘Ÿ(1 βˆ’ πœ‚ π‘Ÿ)(1 βˆ’ π›Ύπ‘Ÿ)(𝑑𝐸 +𝑀1)𝑁𝑇 + (π‘˜4 + π‘˜5)πœ‡π‘£π‘‘πΈπ‘€1 + π‘˜4π‘˜5π‘˜6𝑀1 + 𝑝𝑑π‘₯(πœ‡π‘£ + π‘˜4)𝑀1𝑇2 π‘Ÿ βˆ’ π‘˜(1 βˆ’ π‘“π‘Ÿ)(π‘π‘Ÿ + πœ‚ π‘Ÿπ›Όπ‘Ÿ)(πœ‡π‘£π‘‘πΈ + π‘˜5π‘˜6 + 𝑑π‘₯𝑝𝑇2 π‘Ÿ )π‘‰π‘Ÿ + π›Όπ‘Ÿπ›Ώπ‘Ÿπ‘˜ 2(1 βˆ’ πœ‚π‘Ÿ)(1 βˆ’ π›Ύπ‘Ÿ)(1 βˆ’ π‘“π‘Ÿ)(1 βˆ’ 𝑓𝑠)π‘π‘‡π‘‰π‘Ÿ βˆ’ π‘˜π‘(π‘π‘Ÿ + πœ‚ π‘Ÿπ›Όπ‘Ÿ)𝑑π‘₯π‘‰π‘Ÿπ‘‡2 π‘Ÿ , 𝐸3 = π‘˜4π‘˜5πœ‡π‘£π‘‘πΈπ‘€1 + π‘˜4πœ‡π‘£π‘‘π‘₯𝑝𝑀1𝑇2 π‘Ÿ βˆ’ π‘˜(1 βˆ’ 𝑓𝑠)(1 βˆ’ πœ‚π‘Ÿ)(1 βˆ’ π›Ύπ‘Ÿ)π›Όπ‘Ÿπ›Ώπ‘Ÿπ‘‘πΈπ‘π‘€1𝑇 βˆ’ π‘˜(1 βˆ’ π‘“π‘Ÿ)(π‘π‘Ÿ + πœ‚ π‘Ÿπ›Όπ‘Ÿ) (π‘˜5πœ‡π‘£π‘‘πΈ + 𝑑π‘₯π‘πœ‡π‘£π‘‡2 π‘Ÿ ) π‘‰π‘Ÿ + π‘˜2π›Όπ‘Ÿπ›Ώπ‘Ÿ(1 βˆ’ πœ‚ π‘Ÿ)(1 βˆ’ π›Ύπ‘Ÿ)(1 βˆ’ 𝑓𝑠)(1 βˆ’ π‘“π‘Ÿ)π‘‘πΈπ‘π‘‡π‘‰π‘Ÿ, 𝐴4 = 𝑒1 + 𝛼𝑠 + 𝑏𝑠 + πœ‡π‘£ + 𝛿𝑠 + 𝑑π‘₯𝐸, 𝐡4 = πœ‡π‘£(πœ‡1 + 𝛼𝑠 + 𝑏𝑠) + (𝛿𝑠 + 𝑑π‘₯𝐸)(πœ‡1 + 𝛼𝑠 + 𝑏𝑠 + πœ‡π‘£), 𝐢4 = (𝛿𝑠 + 𝑑π‘₯𝐸)(πœ‡1 + 𝛼𝑠 + 𝑏𝑠)πœ‡π‘£ βˆ’ π‘˜π‘(1 βˆ’ πœ‡π‘š)𝛼𝑠𝛿𝑠(1 βˆ’ πœ‚ 𝑠)(1 βˆ’ 𝛾𝑠), where π‘˜4 = πœ‡1 + π›Όπ‘Ÿ + π‘π‘Ÿ, π‘˜5 = π›Ώπ‘Ÿ + 𝑑π‘₯𝐸, π‘˜6 = πœ‡π‘£ + 𝑑𝐸. The eigenvalues of π½π‘Ÿ will have negative real parts if the roots of 𝑃1(πœ†) = 0 and 𝑄1(πœ†) = 0 have Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 2 (2024) https://internationalpubls.com 468 negative real parts. In this case, the infected steady state πΈπ‘Ÿ, if exists, becomes asymptotically stable. According to Routh-Hurwitz criterion, the equation 𝑃1(πœ†) = 0 will have roots with negative real parts if 𝐴3 > 0, 𝐡3 > 0 , 𝐢3 > 0, 𝐷3 > 0 , 𝐸3 > 0 , 𝐴3𝐡3𝐢3 > 𝐢3 2 + 𝐴3 2𝐷3 and (𝐴3𝐷3 βˆ’ 𝐸3)(𝐴3𝐡3𝐢3 βˆ’ 𝐢3 2 βˆ’ 𝐴3 2𝐷3) > 𝐸3(𝐴3𝐡3 βˆ’ 𝐢3) 2 + 𝐴3𝐸3 2 . In an analogous manner, 𝑄1(πœ†) = 0 will have roots with negative real parts if 𝐴4 > 0, 𝐡4 > 0, 𝐢4 > 0 and 𝐴4𝐡4 βˆ’ 𝐢4 > 0. Proposition 2. The steady state πΈπ‘Ÿ infected with only drug resistant viral strain, if exists, will be asymptotically stable if the following conditions are satisfied i) 𝐴3 > 0, 𝐡3 > 0, 𝐢3 > 0,𝐷3 > 0, 𝐸3 > 0, 𝐴3𝐡3𝐢3 > 𝐢3 2 + 𝐴3 2𝐷3 and (𝐴3𝐷3 βˆ’ 𝐸3)(𝐴3𝐡3𝐢3 βˆ’ 𝐢3 2 βˆ’ 𝐴3 2𝐷3) > 𝐸3(𝐴3𝐡3 βˆ’ 𝐢3) 2 + 𝐴3𝐸3 2 ii) 𝐴4 > 0,𝐡4 > 0, 𝐢4 > 0 and 𝐴4𝐡4 βˆ’ 𝐢4 > 0. 5.3 Stability of steady state (π‘¬π’Ž) infected with both viral strains The Jacobian matrix π‘±π’Ž for the infected steady state πΈπ‘š is obtained by substituting 𝑇 = οΏ½ΜƒοΏ½, 𝑇1 𝑠 = οΏ½ΜƒοΏ½1 𝑠, 𝑇2 𝑠 = οΏ½ΜƒοΏ½2 𝑠, 𝑉𝑠 = �̃�𝑠, 𝑇1 π‘Ÿ = οΏ½ΜƒοΏ½1 π‘Ÿ, 𝑇2 π‘Ÿ = οΏ½ΜƒοΏ½2 π‘Ÿ, π‘‰π‘Ÿ = οΏ½ΜƒοΏ½π‘Ÿ, 𝐸 = οΏ½ΜƒοΏ½ in the Jacobian matrix 𝑱. It is noted that the corresponding characteristic equation |π‘±π’Ž βˆ’ πœ†π‘°| = 0, (5.3) is an eighth degree equation. This matrix π‘±π’Ž could not be divided into blocks so as to get a smaller degree characteristic equations, as in the previous cases. Thus, it is difficult to find the nature of roots for this eighth degree equation analytically. Hence, the nature of the roots of equation (5.3) is checked numerically, whenever required. Proposition 3. The infected steady state πΈπ‘š, if exists, will be asymptotically stable if determinant (5.3) will have all the roots with negative real parts. 6 Numerical example The system (2.1-2.8) of nonlinear ordinary differential equations is solved numerically using MATLAB for the following values of various parameters [1,32,39]. 𝑁 = 1000, π‘‡π‘šπ‘Žπ‘₯ = 1500mm βˆ’3, (𝑠, π‘˜) = (10, 0.000024)mm βˆ’3day βˆ’1; (𝑏𝑠, π‘π‘Ÿ, 𝛼𝑠, π›Όπ‘Ÿ, 𝛿𝑠 , π›Ώπ‘Ÿ) = (0.1, ,0.06, 7, 2, 0.26, 0.16)day βˆ’1 and (π‘Ÿ, πœ‡1, πœ‡π‘£, 𝑝, 𝑑π‘₯, 𝑑𝐸, πœ‡π‘š) = (0.3, 0.015, 2.4, 1.02, 0.01, 0.1, 0.3)day βˆ’1. Initial conditions are chosen as 𝑇(0) = 300 mm βˆ’3 , 𝑇1 𝑠(0) = 𝑇2 𝑠(0) = 𝑉𝑠(0) = 𝑇1 π‘Ÿ(0) = 𝑇2 π‘Ÿ(0) = π‘‰π‘Ÿ(0) = 10 mm βˆ’3 , and 𝐸(0) = 1mm βˆ’3. Numerical example is solved for different combinations of efficacies (𝑓𝑠, πœ‚π‘ , 𝛾𝑠, π‘“π‘Ÿ πœ‚π‘Ÿ , π›Ύπ‘Ÿ) in drug therapy. Case 1. Without drug therapy (i.e., 𝑓𝑠 = πœ‚π‘  = 𝛾𝑠 = π‘“π‘Ÿ = πœ‚π‘Ÿ = π›Ύπ‘Ÿ = 0) For the values of parameters mentioned above, all the roots of the equation (5.3) have negative real parts. Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 2 (2024) https://internationalpubls.com 469 Figure 1: Variations of T-Cell Population (𝑻), Sensitive Viral load (𝑽𝒔), Resistant Viral load (𝑽𝒓) and Total Viral load (𝑽) with time without drug therapy Thus, Proposition 3 implies the existence of an infected steady state πΈπ‘š with both drug sensitive and drug resistant strains, given by (𝑇, 𝑇1 𝑠, 𝑇2 𝑠, 𝑉𝑠, 𝑇1 π‘Ÿ, 𝑇2 π‘Ÿ , π‘‰π‘Ÿ, 𝐸) = (553.75,2.01,9.918,1074.4,0.42,1.01,66.89,111.39) mm βˆ’3 . Both the sensitive (𝑉𝑠 ) and resistant virus strains (π‘‰π‘Ÿ) coexist but as the process of reverse transcription is highly error-prone and as the number of changes per genome is 0.3 per replication cycle, therefore the chance of mutation is quite high. Thus, the drug resistant viral load exists before the initiation of therapy [16] as observed in figure 1. Case 2. Combination of FI and RTI drug therapy The conditions given in Proposition 2 are satisfied for the given parameters and set of efficacies (a, b, c, d) as a.𝑓𝑠 = 0.65, πœ‚π‘  = 0.7, 𝛾𝑠 = 0, π‘“π‘Ÿ = 0.4, πœ‚π‘Ÿ = 0.5, π›Ύπ‘Ÿ = 0 b.𝑓𝑠 = 0.5, πœ‚π‘  = 0.6, 𝛾𝑠 = 0, π‘“π‘Ÿ = 0.25, πœ‚π‘Ÿ = 0.4, π›Ύπ‘Ÿ = 0 c.𝑓𝑠 = 0.4, πœ‚π‘  = 0.5, 𝛾𝑠 = 0, π‘“π‘Ÿ = 0.2, πœ‚π‘Ÿ = 0.3, π›Ύπ‘Ÿ = 0 d.𝑓𝑠 = 0.2, πœ‚π‘  = 0.3, 𝛾𝑠 = 0, π‘“π‘Ÿ = 0.1, πœ‚π‘Ÿ = 0.2, π›Ύπ‘Ÿ = 0 Therefore, for the above values of efficacies (a, b, c, d), the infected steady state with only drug resistant viral strain exists. In each case, the resistant virus dominates the sensitive virus. The sensitive viral load decreases and it vanishes in about 10-20 days as shown in figure 2(2). It is observed in figure 2(3), that the resistant viral load increases with a decrease in efficacy. Thus, the total viral load increases with a decrease in efficacy shown in figure 2(4). Consequently, the T-cell population Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 2 (2024) https://internationalpubls.com 470 decreases with decrease in efficacy as shown in figure 2(1). It is also noted that the system could never reach or attain uninfected steady state for any values of efficacies between 0 to 1. Therefore, because of the presence of resistant viral load, there will never be complete eradication of viral load in spite of the vanishing of sensitive viral load. Figure 2: Variations of T-Cell Population (𝑻), Sensitive Viral load (𝑽𝒔), Resistant Viral load (𝑽𝒓) and Total Viral load (𝑽) with time with combination of FI and RTI drug therapy Case 3. Combination of RTI and PI drug therapy Analogous to previous case, the conditions of Proposition 2 are satisfied for the set of efficacies (a, b, c, d) as a.𝑓𝑠 = 0, πœ‚π‘  = 0.65, 𝛾𝑠 = 0.7, π‘“π‘Ÿ = 0, πœ‚π‘Ÿ = 0.4, π›Ύπ‘Ÿ = 0.5 b.𝑓𝑠 = 0, πœ‚π‘  = 0.5, 𝛾𝑠 = 0.6, π‘“π‘Ÿ = 0, πœ‚π‘Ÿ = 0.25, π›Ύπ‘Ÿ = 0.4 c.𝑓𝑠 = 0, πœ‚π‘  = 0.4, 𝛾𝑠 = 0.5, π‘“π‘Ÿ = 0, πœ‚π‘Ÿ = 0.2, π›Ύπ‘Ÿ = 0.3 d.𝑓𝑠 = 0, πœ‚π‘  = 0.2, 𝛾𝑠 = 0.3, π‘“π‘Ÿ = 0, πœ‚π‘Ÿ = 0.1, π›Ύπ‘Ÿ = 0.2 Thus, the infected steady state πΈπ‘Ÿ with only drug resistant viral strain exists. For each set of efficacy, again as obtained in the previous case the resistant virus dominates the sensitive virus. The sensitive viral load decreases and vanishes in about 10 days, as shown in figure 3(2). The resistant viral load increases with the decrease in efficacies. Consequently, the T-cell population decreases with a decrease in efficacies, and here the total viral load increases as shown in figure 3(1) and 3(4) respectively. It is observed that the total viral load obtained in this case is lower than as obtained in previous case i.e., with FI and RTI drug therapy. Thus, the combination of RTI and PI is more effective than FI and RTI drug therapy. Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 2 (2024) https://internationalpubls.com 471 Figure 3: Variations of T-Cell Population (𝑻), Sensitive Viral load (𝑽𝒔), Resistant Viral load (𝑽𝒓) and Total Viral load (𝑽) with time with combination of RTI and PI drug therapy Case 4. Combination of FI and PI drug therapy In the case of FI and PI, again conditions in Proposition 2 are satisfied for the given set of efficacies (a,b,c,d) a.𝑓𝑠 = 0.65, πœ‚π‘  = 0, 𝛾𝑠 = 0.7, π‘“π‘Ÿ = 0.4, πœ‚π‘Ÿ = 0, π›Ύπ‘Ÿ = 0.5 b.𝑓𝑠 = 0.5, πœ‚π‘  = 0, 𝛾𝑠 = 0.6, π‘“π‘Ÿ = 0.25, πœ‚π‘Ÿ = 0, π›Ύπ‘Ÿ = 0.4 c.𝑓𝑠 = 0.4, πœ‚π‘  = 0, 𝛾𝑠 = 0.5, π‘“π‘Ÿ = 0.2, πœ‚π‘Ÿ = 0, π›Ύπ‘Ÿ = 0.3 d.𝑓𝑠 = 0.2, πœ‚π‘  = 0, 𝛾𝑠 = 0.3, π‘“π‘Ÿ = 0.1, πœ‚π‘Ÿ = 0, π›Ύπ‘Ÿ = 0.2 Thus, the infected steady state πΈπ‘Ÿ with only drug resistant virus strain exists. The results obtained in this case are analogous to that obtained in the case of RTI and PI combination as shown in figure 4. The only difference observed in this case is that the viral load obtained remains much higher in the initial days (i.e., in first 50 days) of infection than that obtained in the previous case as shown in figure 4(4). Again as discussed in case 2 and 3 in this case also there will never be the complete eradication of viral load because of the presence of drug resitant viral load. Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 2 (2024) https://internationalpubls.com 472 Figure 4: Variations of T-Cell Population (𝑻), Sensitive Viral load (𝑽𝒔), Resistant Viral load (𝑽𝒓) and Total Viral load (𝑽) with time with combination of FI and PI drug therapy Case 5. Combination of FI, RTI and PI drug therapy For the set of values of efficacies (a, b, c, d), a.𝑓𝑠 = 0.6, πœ‚π‘  = 0.7, 𝛾𝑠 = 0.65, π‘“π‘Ÿ = 0.5, πœ‚π‘Ÿ = 0.55, π›Ύπ‘Ÿ = 0.45 b.𝑓𝑠 = 0.5, πœ‚π‘  = 0.65, 𝛾𝑠 = 0.6, π‘“π‘Ÿ = 0.4, πœ‚π‘Ÿ = 0.50, π›Ύπ‘Ÿ = 0.35 c.𝑓𝑠 = 0.45, πœ‚π‘  = 0.6, 𝛾𝑠 = 0.5, π‘“π‘Ÿ = 0.2, πœ‚π‘Ÿ = 0.35, π›Ύπ‘Ÿ = 0.3 d.𝑓𝑠 = 0.35, πœ‚π‘  = 0.5, 𝛾𝑠 = 0.4, π‘“π‘Ÿ = 0.15, πœ‚π‘Ÿ = 0.2, π›Ύπ‘Ÿ = 0.25 the conditions given in Proposition 2 are satisfied. Thus, the infected steady state πΈπ‘Ÿ exists. For each of these sets of efficacies, the sensitive viral load vanishes in about 10 days as shown in figure 5(2). The resistant viral load increases with a decrease in efficacies as shown in figure 5(3). Consequently, the T-cell population decreases with a decrease in efficacy, and total viral load increases as observed in figure 5(1) and 5(4) respectively. It is observed that the combined drug therapy may eradicate the sensitive virus but may not be able to eradicate the resistant viral load. Therefore, due to the presence of resistant viral load, combined drug therapy of a very high efficacy fails to eradicate the virus completely. As observed from figure 5, the T-cell population obtained is higher and total viral load obtained is very low in comparison to the case of combined FI and RTI, RTI and PI, and PI and FI drug therapies. Thus, combination of three drugs is very effective to reduce the viral load as compared to the combination of drugs in pairs. Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 2 (2024) https://internationalpubls.com 473 Figure 5: Variations of T-Cell Population (𝑻), Sensitive Viral load (𝑽𝒔), Resistant Viral load (𝑽𝒓) and Total Viral load (𝑽) with time with combination of FI, RTI and PI drug therapy Combined drug therapy and Immune response Figure 6 shows the variation of T-cell population (𝑇), Sensitive Viral load (𝑉𝑠), Resistant Viral load (π‘‰π‘Ÿ) and Total Viral load (𝑉) in presence of combined drug therapy without immune response. It is observed that the T-cell population with the same set of efficacies is very less as compared to the above case i.e., in the presence of combined drug therapy with immune response. Consequently, the resistant viral load strains obtained in this case are also very large as compared to all the cases discussed above with immune response. If the values of the parameters related to the immune system are doubled then it is observed from figure 7 that the T-cell population obtained is more as compared to the above cases with the same set of efficacies as in case 5. Consequently, the resistant viral load as well as the total viral load is very low. Again on triplicating the values of the parameters related to the immune system, it is observed in figure 8 that the values of resistant viral load are very low as compared to the above cases. This shows the importance of a strong immune response with combined drug therapy in reducing the resistant viral load. Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 2 (2024) https://internationalpubls.com 474 Figure 6: Variations of T-Cell Population (𝑻), Sensitive Viral load (𝑽𝒔), Resistant Viral load (𝑽𝒓) and Total Viral load (𝑽) with time with combination of FI, RTI and PI drug therapy and without immune response Figure 7: Variations of T-Cell Population (𝑻), Sensitive Viral load (𝑽𝒔), Resistant Viral load (𝑽𝒓) and Total Viral load (𝑽) with time with combination of FI, RTI and PI drug therapy and doubling the parameters of immune system Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 2 (2024) https://internationalpubls.com 475 Figure 8: Variations of T-Cell Population (𝑻), Sensitive Viral load (𝑽𝒔), Resistant Viral load (𝑽𝒓) and Total Viral load (𝑽) with time with combination of FI, RTI and PI drug therapy and triplicating the parameters of immune system Mutations and Resistant Viral load The presence of mutant virus is observed before the initiation of antiretroviral therapy as discussed in case 1. The variations of T-cell population, Sensitive Viral load 𝑉𝑠, Resistant Viral load π‘‰π‘Ÿ and the Total Viral load (𝑉) as shown in figure 9 for different values of πœ‡π‘š. In this figure, with the increase in the value of πœ‡π‘š , the sensitive viral load decreases, and the resistant viral load increases. Consequently, the total viral load decreases which results in an increase of T-cell population. Figure 9: Variations of T-Cell Population (𝑻), Sensitive Viral load (𝑽𝒔), Resistant Viral load (𝑽𝒓) and Total Viral load (𝑽) with time without drug therapy with different values of ππ’Ž Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 2 (2024) https://internationalpubls.com 476 7 Conclusions The mechanism of the emergence of HIV resistant virus strain under multidrug treatment consisting of Fusion Inhibitor, Reverse Transcriptase Inhibitor, and Protease Inhibitor with an active immune response system is studied by considering a mathematical model of nonlinear differential equations. The present study analyses how a combination of drug therapies (FI and RTI), (RTI and PI), (FI and PI), and (FI, RTI, and PI) in presence of immune response could become more effective drug therapy for both strains of virus. The study also analysed that the combination (RTI and PI) and (FI and PI) are equally effective and also more effective than the (FI and RTI) drug therapy combination. It is further observed that the combined drug therapy (FI, RTI and PI) works very effectively than the combination of therapies in pairs. It also increases the T-cell population to a desired level, which is essential to reduce the risk of disease progression. However, it fails to eradicate the resistant viral load completely. Thus, the drug regimen fails to eradicate the virus completely. It is analyzed that the drug resistant viral load can be reduced by strengthening the immune response system. This interprets that the drug of higher efficacies alone may not be able to eradicate the virus completely because of presence of drug resistant viral load. Whereas with the support of active immune response, the resistant viral load can be reduced and the progression of disease towards AIDS may be prevented even with the moderate efficacies drug combination of FI, RTI and PI. Since the action of immune response of the body is not instant therefore the considered model may be modified further with the introduction of time delay associated with immune response of the body. Also, the drug efficacies are assumed to be constant whereas the concentration of any drug in blood varies continuously due to various factors. Therefore, to get more realistic picture of HIV infection dynamics with viral strains, the model in the present study can be modified further by incorporating the time-dependent drug efficacies. 7.1 Conflict of interest On behalf of all authors, there is no conflict of interest. Compliance with Ethical Standards does not apply to the manuscript. No funding is taken for the manuscript. Data availability statement not applicable. References [1] A.S. Perelson, Modelling the interaction of the immune system with HIV, in: C. Castillo-Chavez (Ed.), Mathematical and Statistical Approaches to AIDS Epidemiology, Springer, Berlin, 1989, p.350. [2] A.S. Perelson, D.E. Kirschner, R. De Boer, Dynamics of HIV Infection of CD4+ T-cells, Math. Biosci. 114 (1993) 81-125. [3] A.S. Perelson, P.W. Nelson, Mathematical Analysis of HIV-1 Dynamics in vivo, SIAM Rev. 41 (1999) 3-44. [4] P. De Leenheer, H.L. Smith, Virus dynamics: a global analysis, SIAM J. Appl. Math. 63 (2003) 1313-1327. [5] L. Wang, M. Li, Mathematical analysis of the global dynamics of a model for HIV infection of CD4+ T-cells, Math. Biosci. 200 (2006) 44-57. [6] D. Wodarz, D.H. Hamer, Infection dynamic in HIV-specific CD4+ T-cells, Math. Biosci. 209 (2007) 14-29. [7] A.S. Perelson, A.U. Neumann, M. Markowitz, J.M. Leonard, D.D. Ho, HIV-1 dynamics in vivo: virion clearance rate, infected cell life span, and viral generation time, Sci. 271 (1996) 1582- 1586. [8] L. Rong, Z. Feng, A.S. Perelson, Mathematical analysis of Age-Structured HIV-1 dynamics with combination antiretroviral therapy, SIAM J. Appl. Math. 67 (2007) 731-756. Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 2 (2024) https://internationalpubls.com 477 [9] P.K. Srivastava, M. Banerjee, P. Chandra, Modeling the drug therapy for HIV infection, J. Biol. Sys. 17 (2009) 213-223. [10] S. Bonhoeffer, R.M. May, G.M. Shaw, M.A. Nowak, Virus dynamics and drug therapy, Proc. Nat. Acad. Sci. USA 94 (1997) 6971-6976. [11] D. Kamboj, M.D. Sharma, Effects of combined drug therapy on HIV-1 infection dynamics, Int. J. Biomath. DOI: 10.1142/S1793524516500650. [12] N.M. Dixit, A.S. Perelson, Complex patterns of viral load decay under antiretroviral therapy: influence of pharmacokinetics and intracellular delay, J. Theor. Biol. 226 (2004) 95-109. [13] M.A. Nowak, S. Bonhoeffer, G.M. Shaw, R.M. May, Anti-viral drug treatment: dynamics of resistance in free virus and infected cell populations, J. Theor. Biol. 184 (1997) 203-217. [14] P.W. Nelson, J. Mittler, A.S. Perelson, Effect of drug efficacy and the eclipse phase of the viral life cycle on estimates of HIV-1 viral dynamic parameters, J. Acquir. Immune. Defic. Syndr. 26 (2001) 405-412. [15] D.E. Kirschner, G.F. Webb, A Mathematical Model of Combined Drug Therapy of HIV Infection, J. Theor. Med. 1 (1997) 25-34. [16] R. Riberio, S. Bonhoeffer, Production of resistant HIV mutant during antiretroviral therapy, Proc. Natl. Acad. Sci. USA 97 (2000) 7861-7866. [17] S. Bonhoeffer, M.A. Nowak, Pre-existence and emergence of drug resistance in HIV-1 infection, Proc. Roy. Soc. Lon. B 264 (1997) 631-637. [18] A.R. McLean, M.A. Nowak, Competition between zidovudine sensitive and zidovudine resistant strains of HIV, AIDS 6 (1992) 71-79. [19] R. Riberio, S. Bonhoeffer, M. Nowak, The frequency of resistant mutant bfore antiviral therapy, AIDS 12 (1998) 461-465. [20] L. Rong, Z. Feng, A.S. Perelson, Emergence of HIV-1 drug resistance during antiretroviral treatment, Bull. Math. Bio. 69 (2007) 2027-2060. [21] C. P. Bhunu, W. Garira, and G. Magombedze, Mathematical Analysis of a Two Strain HIV/AIDS Model with Antiretroviral Treatment. Acta Biotheor 57 (2009) 361–381. [22] L. Rong, M. A. Gilchrist, Z. Feng, A.S. Perelson, Modeling within host HIV-1 dynamics and the evolution of drug resistance: Trade-offs between viral enzyme function and drug susceptibility, J. Theor. Bio. 247 (2007) 804-818. [23] M. Maziane, E. M. Lotfi, K. Hattaf, K. et al. Dynamics of a Class of HIV Infection Models with Cure of Infected Cells in Eclipse Stage. Acta Biotheor 63 (2015) 363–380. [24] M. Rabiu, R. Willie, N. Parumasur, Optimal Control Strategies and Sensitivity Analysis of an HIV/AIDS-Resistant Model with Behavior Change. Acta Biotheor 69 (2021) 543–589. [25] B. Seidu, O. D. Makinde, C.S. Bornaa, Mathematical Analysis of an Industrial HIV/AIDS Model that Incorporates Carefree Attitude Towards Sex. Acta Biotheor 69 (2021) 257–276. [26] A. Poonia, S.P. Chakrabarty, Two strains and drug adherence: An HIV model in the paradigm of community transmission, Nonlinear Dyn 108 (2022), 2767–2792. [27] A.S. Perelson, P. Essunger, D.D. Ho, Dynamics of HIV-1 and CD4+ lymphocytes in vivo, AIDS (Sup A) 11(1997)S17-S24. [28] A.S. Perelson, P.W. Nelson, Modeling viral infections, Proc. Symp. Appl. Math. 59 (2002) 139- 172. [29] D.E. Krischner, G.F. Webb, Understanding drug resistance for monotherapy treatment of HIV infection, Bull. Math. Bio. 59 (1997) 763-785. [30] R.M. Riberio, S. Bonhoeffer, M.A. Nowak, The frequency of resistant mutant virus before antiviral therapy, AIDS. 12 (1998) 461-465. [31] S. Bonhoeffer, M.A. Nowak, R.M. May, G.M. Shaw, Virus dynamics and drug dynamics, Proc. Natl. Acad. Sci. USA. 94 (1997) 6971-6976. [32] P.K. Srivastava, M. Banerjee, P. Chandra, Dynamical model of inhost HIV infection: with drug Communications on Applied Nonlinear Analysis ISSN: 1074-133X Vol 31 No. 2 (2024) https://internationalpubls.com 478 therapy and multi viral strains, J. Biol. Sys. 20 (2012) 303-325. [33] J. W. Mellors, C. R. Rinaldo, P. Gupta, R. M. White, J. A. Todd, L. A. Kingsley, Prognosis in HIV-1 infection predicted by the quantity of virus in plasma, Sci. 272 (1996) 1167-1170. [34] D. Kamboj, M.D. Sharma, Multidrug therapy for HIV infection: dynamics of immune system, Acta biotheoretica 67(2019)129-147. [35] A.S. Perelson, P.W. Nelson, Modeling viral infections, Proc. Sym. Appl. Math. 59 (2002) 139- 172. [36] Y. Wang, Y. Zhou, J. Wu, J. Heffernan, Oscillatory viral dynamics in a delayed HIV pathogenesis model, Math. Biosci. 219 (2009) 104-112. [37] L. Wang, S. Ellermeyer, HIV infection and CD4+ T-cell dynamics, Discrete, Conti. Dyn. Syst. B6 (2006) 1417-1430. [38] J. L. Willems, Stability Theory of Dynamical Systems, Wiley, New York, 1970. [39] P. K. Srivastava, P. Chandra, Hopf bifurcation and periodic solutions in model for the dynamics of HIV and immune response, Diff. Equa. Dyn. Sys. 16 (2008) 77-100.