Skip to main content

ORIGINAL RESEARCH article

Front. Appl. Math. Stat., 14 August 2024
Sec. Mathematical Biology
This article is part of the Research Topic Dynamics and Control of Human and Animal Viral Diseases: Deterministic and Stochastic Models View all articles

Stability of generalized models for HIV-1 dynamics with impaired CTL immunity and three pathways of infection

  • 1Department of Mathematics and Statistics, College of Science, University of Jeddah, Jeddah, Saudi Arabia
  • 2Department of Mathematics, Faculty of Science, King Abdulaziz University, Jeddah, Saudi Arabia

Highly active antiretroviral therapy (HAART) stands out as the most effective treatment for human immunodeficiency virus type-1 (HIV-1). While total eradication is still difficult, HAART can dramatically lower the virus's plasma viral level below the detection threshold. The activation of latently infected cells is considered to be the main obstacle to the total eradication of HIV-1 infection. This study investigates the dynamic characteristics of two generalized HIV-1 infection models taking into account the impairment of cytotoxic T lymphocytes (CTLs). These models include CD4+T cells that are latently infected and equipped with the capability to engage in cell-to-cell infection and elude immune responses. We introduce models featuring three infection pathways: virus-to-cell (VTC), latent cell-to-cell (L-CTC), and active cell-to-cell (A-CTC). The three pathways' infection rates are characterized by general functions, which cover the many types of infection rates documented in the literature. The second model integrates three distinct types of distributed-time delays. We demonstrate the validity of the suggested models, through their well-posedness. We determine the basic reproduction ratio (ℜ0) of the systems. Lyapunov functions and LaSalle's invariance principle are employed to verify that the global stability of both the virus-free steady state (O0) and the virus-persistence steady state (O1). More precisely, O0 achieves global asymptotic stability when ℜ0 ≤ 1, whereas O1 attains global asymptotic stability when ℜ0 > 1. To demonstrate the impact of the parameter values on ℜ0, we examine the sensitivity analysis. It is illustrated that ℜ0 comprises three components, namely 0E, 0P, and 0K, corresponding to the transmissions of VTC, L-CTC, and A-CTC, respectively. Thus, if the L-CTC pathway is disregarded in the HIV-1 infection model, ℜ0 may be underestimated, which could lead to inadequate or erroneous medication therapy focused on eradicating HIV-1 within the body. To demonstrate the associated mathematical outcomes, we conduct numerical simulations through an illustrative example. Specifically, we delve into how the dynamics of HIV-1 are influenced by both immune impairment and time delay. Our findings suggest a significant role of reduced immunity in the progression of the infection. Furthermore, time delays possess the potential to markedly reduce ℜ0, thereby impeding the replication of HIV-1.

1 Introduction

In recent years, there has been a significant upsurge in interest regarding the mathematical modeling of human immunodeficiency virus type-1 (HIV-1) dynamics within the host [1]. Researchers have developed mathematical models that have the potential to substantially advance our knowledge of infection and the medication therapy approaches used to combat it. The dynamics of HIV-1 infection within-host have been explored through various mathematical models. The basic model, governing the time evolution (t) of the healthy CD4+T cell population (J), the infected CD4+T cell population (K), and HIV-1 particles (E), utilizes a set of ordinary differential equations [2]. This model has been extended and modified by using non-integer derivatives [35]. The concept described in reference [2] postulated that CD4+T cells get infected upon coming into contact with HIV-1 particles is known as the virus-to-cell (VTC) pathway. Numerous studies have demonstrated that HIV-1 may spread directly through the development of virological synapses from an infected cell to an uninfected cell, a process referred to as the cell-to-cell (CTC) pathway (see, e.g., [68]). The CTC mechanism can decrease the time required for HIV-1 particle production by 0.9 times and enhance HIV-1 fitness by 3.9 times [9].

A viral infection model with two infection pathways, VTC and CTC, can be formulated as outlined in Graw and Perelson [10]:

{J˙=θηJJ(ρ+ρKK)J,K˙=(ρ+ρKK)JηKK,˙=ϖKη,    (1)

where J=J(t), K=K(t) and E=E(t) represent the concentrations of healthy CD4+T cells, actively infected cells and HIV-1 particles at time t, respectively. The production rate of healthy CD4+T cells is θ, while their death rate is ηJJ. At rates of ρEJE and ρKJK, respectively, the healthy CD4+T cells get infected via VTC and CTC routes. HIV-1 particles are generated at a rate of ϖK from cells that are actively infected. The HIV-1 and actively infected cell death rates are represented by ηEE and ηKK, respectively. Because HIV-1 replication and transmission in lymphoid tissues are so complex, an increasing number of significant biological parameters (such as cell-to-cell transfer, latency, immune response, and time delay). Therefore, many authors have interested in extending model (1) by including different factors.

Cytotoxic T Lymphocyte (CTL) immunity plays an important role against viral infection. Cytotoxic T Lymphocyte (CTLs) kill the cells infected by viruses. Under the effect of CTL immune response, model (1) can be written as [1115]:

{J˙=θ-ηJJ-(ρEE+ρKK)J,K˙=(ρEE+ρKK)J-ηKK-ϕLK,E˙=ϖK-ηEE,L˙=ξLK-ηLL,    (2)

where, L=L(t) is the concentration of the CTLs at time t. CTLs are stimulated, die and kill infected cells at rate, ξLK, ηLL and ϕLK, respectively. It has been documented that HIV-1 can impair CTL immunity and inhibit CTL responsiveness [1622]. By considering the CTL impairment, the fourth equation of system (2) can be modified as [2325]:

L˙=ξK-βLK-ηLL,

where CTL immunological impairment rate is denoted by βLK and CTL immunity stimulation is represented by ξK.

Although effective combination therapy often lowers the HIV-1 plasma viral load to below the detection limit, complete virus elimination remains elusive [26]. Despite the advancements in antiviral medication, one of the major challenges in eradicating HIV-1 lies in the persistence of latently infected cells [26, 27]. Latently infected cells harbor HIV-1 virions but remain dormant until activated. In the HIV-1 infection model, which encompasses two infection pathways, CTL immune impairment and latently infected cells, the formulation can be expressed as: [24, 25]:

{J˙=θ-ηJJ-(ρEE+ρKK)J,P˙=υ(ρEE+ρKK)J-(ψ+ηP)P,K˙=(1-υ)(ρEE+ρKK)J+ψP-ηKK-ϕLK,E˙=ϖK-ηEE,L˙=ξK-βLK-ηLL,

where P=P(t) is the concentration of the latently infected cells at time t. Latently infected cells become active at a rate of ψP and undergo cell death at a rate of ηPP. Among all infections, a fraction υ becomes latent, while the remaining fraction 1 − υ becomes active.

According to the experimental data by Agosto et al. [28], it has been observed that latently infected cells can spread the infection through the CTC mechanism to uninfected cells. Consequently, in recent works [2931], viral infection models were developed by including three pathways of infection, (i) VTC, ρEEJ, (ii) latent cell-to-cell (L-CTC), ρPPJ and (iii) active cell-to-cell (A-CTC), ρKKJ. One crucial factor influencing the spread of viruses is the infection rate [32]. The following types of infection rates were taken into account in the viral dynamics models including CTL immune impairment and/or CTC transmission:

• Bilinear infection rate: ρEEJ [21, 33]; (ρEE+ρKK)J [25]; (ρEE+ρPP+ρKK)J, [34, 35]. This situation implies that the rate of infection is precisely proportional to the product of the concentrations, involving viruses or infected cells interacting with uninfected cells, a phenomenon known as the mass-action principle. However, practical applications reveal that this concept is not universally applicable. The law of mass-action, for instance, won't apply if the quantity of viruses or infected cells surpasses the quantity of uninfected cells. In such a scenario, an increase in virus or infected cell concentration won't result in a rise in infection.

• Saturated infection rate: ρEE1+ϵEEJ [18]; ρEE1+ϵEEJ+ρKK1+ϵKKJ [23, 36], where ϵE,ϵK0.

• Beddington-DeAngelis infection rate ρEE1+ϵEE+ϵJJJ [37], where ϵE,ϵJ0.

• General infection rate: ΥE(J,E) [19]; (ΥE(E)+ΥP(K))ΥK(J) [24]; ΥE(J,E)E+ΥK(J,K)K [38]; ΥE(J,E)E+ΥP(J,P)P+ΥK(J,K)K [29], [30]; ΥE(J,E)+ΥP(J,P)+ΥK(J,K) [31], where ΥE, ΥP and ΥK are general functions.

It should be noted that the model in Alofi and Azoz [24] disregarded the L-CTC route, but included consideration for CTL immune impairment. While the models in Wang et al. [29] and Hattaf and Dutta [30] incorporated three infection pathways (VTC, L-CTC, and A-CTC), they did not account for CTL immunity. In another model discussed in Elaiw and AlShamrani [31], the VTC, L-CTC, and A-CTC pathways were considered, yet the model omitted consideration for CTL immunological impairment.

The current study aims to investigate two HIV-1 infection models with CTL immunological impairment, taking into account three different infection pathways: VTC, L-CTC, and A-CTC which are modeled by general nonlinear functions. Additionally, three kinds of distributed-time delays are incorporated into the second model. We examine the fundamental characteristics of the model's solutions in addition to the steady states' global stability. We provide numerical simulations to validate the theoretical results. We wrap up by talking about the outcomes. Our proposed model can be seen as a generalization of several HIV-1 infection models with CTL impairment presented in the literature. Our system's basic reproduction ratio (ℜ0) is more accurate than models that do not account for the L-CTC route. As a result, the acquired (ℜ0) may result in appropriate drug therapy aimed at eliminating HIV-1 from the body.

2 Model with CTL immune impairment, general infection rate and three pathways of infection

2.1 System overview

We introduce a model for HIV-1 that incorporates CTL immune impairment while also accounting for the infection of healthy CD4+T cells through VTC, L-CTC, and A-CTC transmission modes. The model is outlined as follows:

{J˙=θ-ηJJ-ΥE(J,E)-ΥP(J,P)-ΥK(J,K),P˙=ΥE(J,E)+ΥP(J,P)+ΥK(J,K)-(ψ+ηP)P,K˙=ψP-ηKK-ϕLK,E˙=ϖK-ηEE,L˙=ξK-ηLL-βLK.    (3)

The general functions ΥE(J,E), ΥP(J,P) and ΥK(J,K) represent, respectively, the VTC, L-CTC and A-CTC routes of infection. The meanings assigned to all remaining parameters and variables remain in line with the explanations provided in Section 1. Model (3) is very general as it considers nonlinear incidences (ΥE(J,E), ΥP(J,P) and ΥK(J,K)). Additionally, it is crucial to highlight that model (3) encompasses a broad spectrum of pre-existing models, as seen in earlier research [24, 30, 31].

Define Z(J), Z{E,P,K} as:

E(J)=limE0+ΥE(J,E)E=ΥE(J,0)E,P(J)=limP0+ΥP(J,P)P=ΥP(J,0)P,K(J)=limK0+ΥK(J,K)K=ΥK(J,0)K.

The paper relies on the following conditions regarding the functions ΥZ(J,Z), Z{E,P,K} [39, 40]:

Condition C1. ΥZ(J,Z) is continuously differentiable, ΥZ(J,Z)>0 and ΥZ(0,Z)=ΥZ(J,0)=0 for all J>0 and Z{E,P,K}>0.

Condition C2. ΥZ(J,Z)J>0, ΥZ(J,Z)Z>0 for all J>0 and Z{E,P,K}>0.

Condition C3. Z(J)>0 and Z(J)>0 for all J>0, Z{E,P,K}.

Condition C4. ΥZ(J,Z)Z is non-increasing with respect to Z for all Z{E,P,K}>0.

2.2 Fundamental characteristics

2.2.1 Model's well-posedness

To guarantee the biological acceptability of our model, the following assertion establishes a constrained domain for the compartment concentrations. More specifically, the concentrations shouldn't go negative or unbounded.

Proposition 2.1. Assume that Condition C1 is satisfied, then the solutions of system (3) are nonnegative and bounded.

Proof. From system (3) we get the following:

J˙J=0=θ>0, P˙P=0=ΥE(J,E)+ΥK(J,K)0,for all J,E,K0,K˙K=0=ψP0,for all P0, E˙E=0=ϖK0,for all K0, L˙L=0=ξK0,for all K0.

Thus, according to Proposition B.7 of Smith and Waltman [41], we conclude that, for all t ≥ 0, we have (J(t),P(t),K(t),E(t),L(t))05, whenever (J(0),P(0),K(0),E(0),L(0))05. Let

Φ=J+P+K+ηK2ϖE+ηK4ξL.

Then, we have

Φ˙=J˙+P˙+K˙+ηK2ϖE˙+ηK4ξL˙   =θ-ηJJ-ηPP-ηK4K-(ϕ+ηKβ4ξ)LK   -ηKηE2ϖE-ηKηL4ξL   θ-ηJJ-ηPP-ηK4K-ηKηE2ϖE-ηKηL4ξL   θ-ε(J+P+K+ηK2ϖE+ηK4ξL)=θ-εΦ,

where ε=min{ηJ,ηP,ηK4,ηE,ηL}. Hence

Φ(t)e-εt(Φ(0)-θε)+θε.

This yields 0 ≤ Φ(t) ≤ Γ1 if Φ(0) ≤ Γ1, where Γ1=θε.

Considering that all state variables have non-negative values, 0J(t),P(t),K(t)Γ1, 0E(t)Γ2, and 0L(t)Γ3, for all t ≥ 0 if J(0)+P(0)+K(0)+ηK2ϖE(0)+ηK4ξL(0)Γ1, where Γ2=2ϖΓ1ηK and Γ3=4ξΓ1ηK. To sum up, the boundedness of J(t),P(t),K(t),E(t) and L(t) implies that the solutions of system (3) are bounded.

2.2.2 Derivation of reproduction ratios and steady states

Proposition 2.2. Assuming that Conditions C1C4 are met, there exists a positive basic reproduction ratio 0=ϖψE(J0)+ηEηKP(J0)+ψηEK(J0)ηEηK(ψ+ηP) for system (3) in a way that

 (i) ensures the system consistently maintains a virus-free steady state, denoted as O0, and

(ii) if ℜ0 > 1, the system additionally possesses a virus-persistence steady state, denoted as O1.

Proof. It is observed that a virus-free steady state labeled as O0=(J0,0,0,0,0) consistently exists in system (3), where J0=θηJ. In the upcoming analysis, we plan to calculate the basic reproduction ratio for system (3) using the next-generation matrix (NGM) approach, as introduced by Driessche and Watmough [42]. The model (3) features infected compartments, denoted as (P,K,E), and is defined by matrices that describe the nonlinear components related to the new infection Ω^1 and the outflow term Λ^1 within these compartments as given below.

Ω^1=(ΥE(J,E)+ΥP(J,P)+ΥK(J,K)00),Λ^1=((ψ+ηP)P-ψP+ηKK+ϕLK-ϖK+ηEE).

The derivatives of Ω^1 and Λ^1 at the steady state O0 are calculated, yielding the following matrices:

Ω1=(ΥP(J0,0)PΥK(J0,0)KΥE(J0,0)E000000),Λ1=(ψ+ηP00-ψηK00-ϖηE).

The form of the NGM is as follows:

Ω1Λ1-1=(ϖψE(J0)ηEηK(ψ+ηP)+P(J0)ψ+ηP+ψK(J0)ηK(ψ+ηP)ϖE(J0)ηEηK+K(J0)ηKE(J0)ηE000000).

The spectral radius of matrix Ω1Λ1-1, symbolized as ℜ0, defines the basic reproduction ratio. Its calculation is determined by the following expression:

0=ϖψE(J0)+ηEηKP(J0)+ψηEK(J0)ηEηK(ψ+ηP)   =0E+0P+0K,    (4)

where

0E=ϖψE(J0)ηEηK(ψ+ηP),    0P=P(J0)ψ+ηP,0K=ψK(J0)ηK(ψ+ηP).

The clinical relevance of the parameter ℜ0 holds significance as it plays a crucial role in determining whether the HIV-1 infection will progress into a chronic state. Within this framework, ℜ0Z symbolizes the average quantities of secondary infected cells generated from engagements with viruses and infected cells, where Z{E,P,K}. To extend our exploration beyond O0, we scrutinize (J,P,K,E,L) as a potential steady state governed by the following set of algebraic equations:

0=θ-ηJJ-ΥE(J,E)-ΥP(J,P)-ΥK(J,K),    (5)
0=ΥE(J,E)+ΥP(J,P)+ΥK(J,K)-(ψ+ηP)P,    (6)
0=ψP-ηKK-ϕLK,    (7)
0=ϖK-ηEE,    (8)
0=ξK-ηLL-βLK.    (9)

From Equations (8, 9), we get

E=Θ1(K)=ϖKηE, L=ξKηL+βK.    (10)

Clearly, Θ1(0) = 0. The following outcome is obtained after inserting Equation (10) into Equation (7):

P=Θ2(K)=ηLηKK+(ϕξ+ηKβ)K2ψ(ηL+βK).    (11)

It is clear that Θ2(0) = 0. From Equations (5, 6), we get

θ-ηJJ=(ψ+ηP)P.    (12)

The following outcome is obtained after inserting Equation (11) into Equation (12):

J=Θ3(K)   =1ηJ(θ-(ψ+ηP)(ηLηKK+(ϕξ+ηKβ)K2)ψ(ηL+βK)).        (13)

Note that Θ3(0)=J0. Upon substitution of Equations (10, 11, 13) into Equation (6), the result is as follows:

ΥE(Θ3(K),Θ1(K))+ΥP(Θ3(K),Θ2(K))+ΥK(Θ3(K),K)   -(ψ+ηP)(ηLηKK+(ϕξ+ηKβ)K2)ψ(ηL+βK)=0.    (14)

From Equation (14), we have

1. When K=0, the virus-free steady state O0 is derived from Equations (10, 11, 12, 13).

2. When K0, let us define a function Ψ(K) on [0, ∞) as:

Ψ(K)=ΥE(Θ3(K),Θ1(K))+ΥP(Θ3(K),Θ2(K))         +ΥK(Θ3(K),K)         -(ψ+ηP)(ηLηKK+(ϕξ+ηKβ)K2)ψ(ηL+βK).

We have Ψ(0) = 0. Let K^ be such that Θ3(K^)=0, i.e.,

J0-(ψ+ηP)(ηLηKK^+(ϕξ+ηKβ)K^2)ηJψ(ηL+βK^)=0,

which implies that

(ϕξ+ηKβ)(ψ+ηP)K^2+(ηLηK(ψ+ηP)-ηJψβJ0)K^-ηJψηLJ0=0.    (15)

Therefore, Equation (15) has a positive solution

K^=-B+B2-4AC2A,

where

A=(ϕξ+ηKβ)(ψ+ηP),B=ηLηK(ψ+ηP)-ηJψβJ0,C=-ηJψηLJ0.

We can see that

Ψ(K^)=ΥE(0,Θ1(K^))+ΥP(0,Θ2(K^))+ΥK(0,K^)      -(ψ+ηP)(ηLηKK^+(ϕξ+ηKβ)K^2)ψ(ηL+βK^)      =-(ψ+ηP)(ηLηKK^+(ϕξ+ηKβ)K^2)ψ(ηL+βK^)<0.

In addition

Ψ(K)=Θ3(K)ΥE(J,E)J+Θ1(K)ΥE(J,E)E         +Θ3(K)ΥP(J,P)J+Θ2(K)ΥP(J,P)P         +Θ3(K)ΥK(J,K)J+ΥK(J,K)K         -ψ+ηPψ(ηK+ϕξ(2ηL+βK)K(ηL+βK)2).

Condition C1 implies that ΥZ(J0,0)J=0, Z{E,P,K}. Then

Ψ(0)=Θ1(0)ΥE(J0,0)E+Θ2(0)ΥP(J0,0)P      +ΥK(J0,0)K-ηK(ψ+ηP)ψ      =ϖηEΥE(J0,0)E+ηKψΥP(J0,0)P+ΥK(J0,0)K      -ηK(ψ+ηP)ψ      =ϖE(J0)ηE+ηKP(J0)ψ+K(J0)-ηK(ψ+ηP)ψ      =ηK(ψ+ηP)ψ(0-1),

where ℜ0 is defined in Equation (4). Therefore, if ℜ0 > 1 then Ψ′(0) > 0 and there exists K1(0,K^) such that Ψ(K1)=0. Let K=K1 in Equation (5) and define

(J)=θ-ηJJ-ΥE(J,Θ1(K1))-ΥP(J,Θ2(K1))-ΥK(J,K1).

Subsequently, based on Condition C1, we obtain ℵ(0) = θ > 0 and

(J0)=-(ΥE(J0,Θ1(K1))+ΥP(J0,Θ2(K1))+ΥK(J0,K1))<0.

Under Condition C2, it follows that (J) decreases strictly as a function of J. Consequently, there is a unique value J1 within the interval (0,J0) for which (J1) equals zero. Moreover, considering Equations (10, 11), we find

P1=ηLηKK1+(ϕξ+βηK)K12ψ(ηL+βK1), E1=ϖK1ηE,L1=ξK1ηL+βK1.

The presence of the virus-persistence steady state O1=(J1,P1,K1,E1,L1) becomes evident when ℜ0 > 1.

2.2.3 Global stability of steady states

The forthcoming theorems explore the global asymptotic stability of both virus-free and virus-persistence steady states. The development of the Lyapunov function will adhere to the approach outlined in Korobeinikov [43].

In the remaining work, we focus on a function Ξi(J,P,K,E,L), and we characterize Mi  as the largest invariant subset of

Mi={(J,P,K,E,L):dΞidt=0}, i=1,...,4.

In order to examine whether the steady state O0, is globally asymptotically stable (G.A.S), we necessitate adherence to Condition C5 as outlined below [40]:

Condition C5. The supremum of Z(J)E(J) on (0,θ/ηJ] is attained at J=θηJ, where Z{P,K}.

Theorem 2.1. For system (3), let0 ≤ 1 and Conditions C1C5 are satisfied, then O0 is G.A.S.

Proof. Let's contemplate a potential Lyapunov function as:

Ξ1=J-J0-J0JE(J0)E(λ)dλ+P+ϖE(J0)+ηEK(J0)ηKηEK               +E(J0)ηEE+ϕ(ϖE(J0)+ηEK(J0))2ξηKηEL2.

Evidently, Ξ1(J,P,K,E,L)>0 for every J,P,K,E,L>0, as well as Ξ1(J0,0,0,0,0)=0. Computing dΞ1dt, we derive

dΞ1dt=(1(J0)(J))J˙+P˙+ϖ(J0)+ηK(J0)ηKηK˙         +(J0)η˙+ϕ(ϖ(J0)+ηK(J0))ξηKη˙         =(1(J0)(J))(θηJJΥ(J,)ΥP(J,P)         ΥK(J,K))+Υ(J,)+ΥP(J,P)+ΥK(J,K)         (ψ+ηP)P         +ϖ(J0)+ηK(J0)ηKη(ψPηKKϕK)         +(J0)η(ϖKη)         +ϕ(ϖ(J0)+ηK(J0))ξηKη(ξKηβK)         =(1(J0)(J))(θηJJ)         +((J0)(J)Υ(J,)(J0))         +((J0)(J)ΥP(J,P)P(ψ+ηP)         +ψ(ϖ(J0)+ηK(J0))ηKη)P         +((J0)(J)ΥK(J,K)KK(J0))K         ϕ(ϖ(J0)+ηK(J0))ξηKη(η+βK)2.    (16)

From Condition C4, we have

ΥE(J,E)ElimE0+ΥE(J,E)E=E(J),    for all J,E>0,ΥP(J,P)PlimP0+ΥP(J,P)P=P(J),    for all J,P>0,ΥK(J,K)KlimK0+ΥK(J,K)K=K(J),    for all J,K>0.    (17)

Then

E(J0)E(J)ΥE(J,E)E-E(J0)E(J0)E(J)E(J)-E(J0)=0.

Further, Conditions C4 and C5 in case of J0=θηJ imply that

E(J0)E(J)ΥP(J,P)P-(ψ+ηP)+ψ(ϖE(J0)+ηEK(J0))ηKηEE(J0)P(J)E(J)-(ψ+ηP)+ψ(ϖE(J0)+ηEK(J0))ηKηEE(J0)P(J0)E(J0)-(ψ+ηP)+ψ(ϖE(J0)+ηEK(J0))ηKηE=P(J0)-(ψ+ηP)+ψ(ϖE(J0)+ηEK(J0))ηKηE=ψ(ϖE(J0)+ηEK(J0))+ηKηEP(J0)ηKηE-(ψ+ηP)=(ψ+ηP)(0-1).

Furthermore

E(J0)E(J)ΥK(J,K)K-K(J0)E(J0)K(J)E(J)-K(J0)E(J0)K(J0)E(J0)-K(J0)=0.

Upon conducting a computation and substituting the value J0=θ/ηJ, Equation (16) transforms into the following expression:

dΞ1dtθ(1-E(J0)E(J))(1-JJ0)+(ψ+ηP)(0-1)P-ϕ(ϖE(J0)+ηEK(J0))ξηKηE(ηL+βK)L2.

From Condition C3, we have

(1-E(J0)E(J))(1-JJ0)0.

Clearly, dΞ1dt0 when ℜ0 ≤ 1 and dΞ1dt=0 when J=J0 and P=K=L=0. All solutions convey to M1 , wherein every element fulfills J(t)=J0 and P(t)=K(t)=L(t)=0 for all t [44]. Following that, the first equation of system (3) in conjunction with Condition C1 results in

0=J˙=θ-ηJJ0-ΥE(J0,E(t))ΥE(J0,E(t))=0E(t)=0,for all t.

Therefore, with M1 ={O0}, our conclusion asserts that O0 is G.A.S under the condition ℜ0 ≤ 1, in accordance with LaSalle's invariance principle (L.I.P) [45].

In investigating the global stability of O1, we define 𝐹(α) = α − 1 − ln(α). Besides, the following remarks and Condition will be essential.

Remark 2.1. Considering Conditions C2 and C4, it is established that, for every positive J,E,E*, we have

(ΥE(J,E)E-ΥE(J,E*)E*)(ΥE(J,E)-ΥE(J,E*))0,

which implies that

(ΥE(J,E)ΥE(J,E*)-EE*)(1-ΥE(J,E*)ΥE(J,E))0,    (18)

where E*{E1,E˜1}.

Condition C6. For any J within the interval (0,θ/ηJ), and positive values of P and K, the following condition is satisfied

1.     (ΥP(J,P)ΥE(J,E*)P-ΥP(J*,P*)ΥE(J*,E*)P*)           ×(ΥP(J,P)ΥE(J,E*)-ΥP(J*,P*)ΥE(J*,E*))0,
2.     (ΥK(J,K)ΥE(J,E*)K-ΥK(J*,K*)ΥE(J*,E*)K*)           ×(ΥK(J,K)ΥE(J,E*)-ΥK(J*,K*)ΥE(J*,E*))0,

where J*{J1,J~1}>0,P*{P1,P~1}>0,K*{K1,K~1}>0, and E*{E1,E˜1}>0.

Remark 2.2. From Condition C6, For any J within the interval (0,θ/ηJ), and positive values of P and K, we get

(ΥE(J*,E*)ΥP(J,P)ΥE(J,E*)ΥP(J*,P*)-PP*)×(1-ΥE(J,E*)ΥP(J*,P*)ΥE(J*,E*)ΥP(J,P))0,    (19)
(ΥE(J*,E*)ΥK(J,K)ΥE(J,E*)ΥK(J*,K*)-KK*)×(1-ΥE(J,E*)ΥK(J*,K*)ΥE(J*,E*)ΥK(J,K))0,    (20)

where J*{J1,J~1}>0,P*{P1,P~1}>0,K*{K1,K~1}>0, and E*{E1,E˜1}>0.

Theorem 2.2. Given ℜ0 > 1 and the fulfillment of Conditions C1-C4 and C6, it can be concluded that the steady state O1 for system (3) is guaranteed to be G.A.S.

Proof. We formulate Ξ2(J,P,K,E,L) given by Equation (21) as:

Ξ2=J-J1-J1JΥE(J1,E1)ΥE(λ,E1)dλ+P1𝐹(PP1)   +ϖK1ΥE(J1,E1)+ηEE1ΥK(J1,K1)ηEE1(ηK+ϕL1)𝐹(KK1)   +ΥE(J1,E1)ηEг(EE1)   +ϕ(ϖK1ΥE(J1,E1)+ηEE1ΥK(J1,K1))2ηEK1E1(ηK+ϕL1)(ξ-βL1)(L-L1)2.    (21)

Steady state condition Equation (9) guarantees that ξ-βL1=ηLL1K1>0. Clearly, Ξ2 is positive definite. Calculating dΞ2dt:

dΞ2dt=(1Υ(J1,1)Υ(J,1))J˙+(1P1P)P˙       +ϖK1Υ(J1,1)+η1ΥK(J1,K1)ηK11(ηK+ϕ1)(1K1K)K˙       +Υ(J1,1)η1(11)˙       +ϕ(ϖK1Υ(J1,1)+η1ΥK(J1,K1))ηK11(ηK+ϕ1)(ξβ1)(1)˙       =(1Υ(J1,1)Υ(J,1))(θηJJΥ(J,)ΥP(J,P)       ΥK(J,K))       +(1P1P)(Υ(J,)+ΥP(J,P)+ΥK(J,K)       (ψ+ηP)P)       +ϖK1Υ(J1,1)+η1ΥK(J1,K1)ηK11(ηK+ϕ1)(1K1K)       ×(ψPηKKϕK)       +Υ(J1,1)η1(11)(ϖKη)       +ϕ(ϖK1Υ(J1,1)+η1ΥK(J1,K1))ηK11(ηK+ϕ1)(ξβ1)(1)       ×(ξKηβK)       =(1Υ(J1,1)Υ(J,1))(θηJJ)+Υ(J1,1)Υ(J,1)       ×(Υ(J,)+ΥP(J,P)+ΥK(J,K))       (Υ(J,)+ΥP(J,P)+ΥK(J,K))P1P       (ψ+ηP)P+(ψ+ηP)P1       +ψ(ϖK1Υ(J1,1)+η1ΥK(J1,K1))ηK11(ηK+ϕ1)P       ψ(ϖK1Υ(J1,1)+η1ΥK(J1,K1))ηK11(ηK+ϕ1)PK1K       ηK(ϖK1Υ(J1,1)+η1ΥK(J1,K1))ηK11(ηK+ϕ1)(KK1)       ϕ(ϖK1Υ(J1,1)+η1ΥK(J1,K1))ηK11(ηK+ϕ1)(KK1)       +Υ(J1,1)ϖKη1Υ(J1,1)ϖKηΥ(J1,1)1       +Υ(J1,1)       +ϕ(ϖK1Υ(J1,1)+η1ΥK(J1,K1))ηK11(ηK+ϕ1)(ξβ1)(1)       ×(ξKηβK).    (22)

Using the following steady state conditions for O1

θ=ηJJ1+ΥE(J1,E1)+ΥP(J1,P1)+ΥK(J1,K1),ΥE(J1,E1)+ΥP(J1,P1)+ΥK(J1,K1)=(ψ+ηP)P1,K1=ψP1ηK+ϕL1,    E1=ϖK1ηE,  ξK1=ηLL1+βL1K1,

we get

ΥE(J1,E1)+ΥK(J1,K1)=ϖK1ΥE(J1,E1)+ηEE1ΥK(J1,K1)ηEE1=ψP1(ϖK1ΥE(J1,E1)+ηEE1ΥK(J1,K1))ηEK1E1(ηK+ϕL1).

Consequently, the representation of Equation (22) will be as follows:

dΞ2dt=(1Υ(J1,1)Υ(J,1))(ηJJ1ηJJ)       +(1Υ(J1,1)Υ(J,1))(Υ(J1,1)+ΥP(J1,P1)       +ΥK(J1,K1))+Υ(J1,1)Υ(J,)Υ(J,1)       +ΥP(J1,P1)Υ(J1,1)ΥP(J,P)Υ(J,1)ΥP(J1,P1)       +ΥK(J1,K1)Υ(J1,1)ΥK(J,K)Υ(J,1)ΥK(J1,K1)       Υ(J1,1)Υ(J,)P1Υ(J1,1)P       ΥP(J1,P1)ΥP(J,P)P1ΥP(J1,P1)P       ΥK(J1,K1)ΥK(J,K)P1ΥK(J1,K1)P(ψ+ηP)P1PP1       +ΥP(J1,P1)+ΥK(J1,K1)       +ψP1(ϖK1Υ(J1,1)+η1ΥK(J1,K1))ηK11(ηK+ϕ1)PP1       ψP1(ϖK1Υ(J1,1)+η1ΥK(J1,K1))ηK11(ηK+ϕ1)PK1P1K       ηK(ϖK1Υ(J1,1)+η1ΥK(J1,K1))ηK11(ηK+ϕ1)(KK1)       ϕ(ϖK1Υ(J1,1)+η1ΥK(J1,K1))ηK11(ηK+ϕ1)(KK1)       +ϕ(ϖK1Υ(J1,1)+η1ΥK(J1,K1))ηK11(ηK+ϕ1)1(KK1)       ϕ(ϖK1Υ(J1,1)+η1ΥK(J1,K1))ηK11(ηK+ϕ1)1(KK1)       +Υ(J1,1)KK1Υ(J1,1)K1K1Υ(J1,1)1       +Υ(J1,1)       +ϕ(ϖK1Υ(J1,1)+η1ΥK(J1,K1))ηK11(ηK+ϕ1)(ξβ1)(1)       ×(ξKηβK       ξK1+η1+β1K1+β1Kβ1K)       =ηJJ1(1Υ(J1,1)Υ(J,1))(1JJ1)       +(2Υ(J1,1)Υ(J,1))(Υ(J1,1)+ΥP(J1,P1)       +ΥK(J1,K1))+Υ(J1,1)Υ(J,)Υ(J,1)       +ΥP(J1,P1)Υ(J1,1)ΥP(J,P)Υ(J,1)ΥP(J1,P1)       +ΥK(J1,K1)Υ(J1,1)ΥK(J,K)Υ(J,1)ΥK(J1,K1)       Υ(J1,1)Υ(J,)P1Υ(J1,1)P       ΥP(J1,P1)ΥP(J,P)P1ΥP(J1,P1)PΥK(J1,K1)       ×ΥK(J,K)P1ΥK(J1,K1)P(Υ(J1,1)+ΥP(J1,P1)       +ΥK(J1,K1))PP1+(Υ(J1,1)+ΥK(J1,K1))PP1       (Υ(J1,1)+ΥK(J1,K1))PK1P1K       ϖK1Υ(J1,1)+η1ΥK(J1,K1)ηK11(ηK+ϕ1)(ηK+ϕ1)(KK1)       ϕ(ϖK1Υ(J1,1)+η1ΥK(J1,K1))ηK11(ηK+ϕ1)       ×(1)(KK1)+Υ(J1,1)KK1Υ(J1,1)K1K1       Υ(J1,1)1+Υ(J1,1)       +ϕ(ϖK1Υ(J1,1)+η1ΥK(J1,K1))ηK11(ηK+ϕ1)(ξβ1)       ×(ξβ1)(1)(KK1)       ϕ(ϖK1Υ(J1,1)+η1ΥK(J1,K1))ηK11(ηK+ϕ1)(ξβ1)       ×(η+βK)(1)2.

This implies that

dΞ2dt=ηJJ1(1-ΥE(J1,E1)ΥE(J,E1))(1-JJ1)   +(2-ΥE(J1,E1)ΥE(J,E1))(ΥE(J1,E1)+ΥP(J1,P1)   +ΥK(J1,K1))+ΥE(J1,E1)ΥE(J,E)ΥE(J,E1)   +ΥP(J1,P1)ΥE(J1,E1)ΥP(J,P)ΥE(J,E1)ΥP(J1,P1)   +ΥK(J1,K1)ΥE(J1,E1)ΥK(J,K)ΥE(J,E1)ΥK(J1,K1)   -ΥE(J1,E1)ΥE(J,E)P1ΥE(J1,E1)P   -ΥP(J1,P1)ΥP(J,P)P1ΥP(J1,P1)P   -ΥK(J1,K1)ΥK(J,K)P1ΥK(J1,K1)P-ΥP(J1,P1)PP1   -(ΥE(J1,E1)+ΥK(J1,K1))PK1P1K   -ϖK1ΥE(J1,E1)+ηEE1ΥK(J1,K1)ηEE1(KK1-1)   +ΥE(J1,E1)KK1-ΥE(J1,E1)KE1K1E   -ΥE(J1,E1)EE1+ΥE(J1,E1)   -ϕ(ΥE(J1,E1)+ΥK(J1,K1))K1(ηK+ϕL1)(ξ-βL1)(ηL+βK)(L-L1)2   =ηJJ1(1-ΥE(J1,E1)ΥE(J,E1))(1-JJ1)   +(2-ΥE(J1,E1)ΥE(J,E1))(ΥE(J1,E1)+ΥP(J1,P1)   +ΥK(J1,K1))+ΥE(J1,E1)ΥE(J,E)ΥE(J,E1)   +ΥP(J1,P1)ΥE(J1,E1)ΥP(J,P)ΥE(J,E1)ΥP(J1,P1)   +ΥK(J1,K1)ΥE(J1,E1)ΥK(J,K)ΥE(J,E1)ΥK(J1,K1)   -ΥE(J1,E1)ΥE(J,E)P1ΥE(J1,E1)P-ΥP(J1,P1)ΥP(J,P)P1ΥP(J1,P1)P   -ΥK(J1,K1)ΥK(J,K)P1ΥK(J1,K1)P-ΥP(J1,P1)PP1   -(ΥE(J1,E1)+ΥK(J1,K1))PK1P1K-(ΥE(J1,E1)   +ΥK(J1,K1))(KK1-1)   +ΥE(J1,E1)KK1-ΥE(J1,E1)KE1K1E-ΥE(J1,E1)EE1   +ΥE(J1,E1)   -ϕ(ΥE(J1,E1)+ΥK(J1,K1))K1(ηK+ϕL1)(ξ-βL1)(ηL+βK)(L-L1)2   =ηJJ1(1-ΥE(J1,E1)ΥE(J,E1))(1-JJ1)   +(2-ΥE(J1,E1)ΥE(J,E1))(ΥE(J1,E1)+ΥP(J1,P1)   +ΥK(J1,K1))+ΥE(J1,E1)ΥE(J,E)ΥE(J,E1)   +ΥP(J1,P1)ΥE(J1,E1)ΥP(J,P)ΥE(J,E1)ΥP(J1,P1)   +ΥK(J1,K1)ΥE(J1,E1)ΥK(J,K)ΥE(J,E1)ΥK(J1,K1)   -ΥE(J1,E1)ΥE(J,E)P1ΥE(J1,E1)P-ΥP(J1,P1)ΥP(J,P)P1ΥP(J1,P1)P   -ΥK(J1,K1)ΥK(J,K)P1ΥK(J1,K1)P-ΥP(J1,P1)PP1   -(ΥE(J1,E1)+ΥK(J1,K1))PK1P1K   -ΥK(J1,K1)KK1+2ΥE(J1,E1)+ΥK(J1,K1)   -ΥE(J1,E1)KE1K1E-ΥE(J1,E1)EE1   -ϕ(ΥE(J1,E1)+ΥK(J1,K1))K1(ηK+ϕL1)(ξ-βL1)(ηL+βK)(L-L1)2   =ηJJ1(1-ΥE(J1,E1)ΥE(J,E1))(1-JJ1)   +ΥE(J1,E1)(ΥE(J,E)ΥE(J,E1)-EE1)   +ΥE(J1,E1)(4-ΥE(J1,E1)ΥE(J,E1)-ΥE(J,E)P1ΥE(J1,E1)P   -PK1P1K-KE1K1E)   +ΥP(J1,P1)(ΥE(J1,E1)ΥP(J,P)ΥE(J,E1)ΥP(J1,P1)-PP1)   +ΥP(J1,P1)(2-ΥE(J1,E1)ΥE(J,E1)-ΥP(J,P)P1ΥP(J1,P1)P)   +ΥK(J1,K1)(ΥE(J1,E1)ΥK(J,K)ΥE(J,E1)ΥK(J1,K1)-KK1)   +ΥK(J1,K1)(3-ΥE(J1,E1)ΥE(J,E1)-ΥK(J,K)P1ΥK(J1,K1)P   -PK1P1K)-ϕ(ΥE(J1,E1)+ΥK(J1,K1))K1(ηK+ϕL1)(ξ-βL1)   ×(ηL+βK)(L-L1)2.

Then

dΞ2dt=ηJJ1(1-ΥE(J1,E1)ΥE(J,E1))(1-JJ1)   +ΥE(J1,E1)(ΥE(J,E)ΥE(J,E1)-EE1-1+ΥE(J,E1)EΥE(J,E)E1)   +ΥE(J1,E1)(5-ΥE(J1,E1)ΥE(J,E1)-ΥE(J,E)P1ΥE(J1,E1)P   -PK1P1K-KE1K1E-ΥE(J,E1)EΥE(J,E)E1)   +ΥP(J1,P1)(ΥE(J1,E1)ΥP(J,P)ΥE(J,E1)ΥP(J1,P1)-PP1-1   +ΥE(J,E1)ΥP(J1,P1)PΥE(J1,E1)ΥP(J,P)P1)   +ΥP(J1,P1)(3-ΥE(J1,E1)ΥE(J,E1)-ΥP(J,P)P1ΥP(J1,P1)P   -ΥE(J,E1)ΥP(J1,P1)PΥE(J1,E1)ΥP(J,P)P1)   +ΥK(J1,K1)(ΥE(J1,E1)ΥK(J,K)ΥE(J,E1)ΥK(J1,K1)-KK1-1   +ΥE(J,E1)ΥK(J1,K1)KΥE(J1,E1)ΥK(J,K)K1)   +ΥK(J1,K1)(4-ΥE(J1,E1)ΥE(J,E1)-ΥK(J,K)P1ΥK(J1,K1)P   -PK1P1K-ΥE(J,E1)ΥK(J1,K1)KΥE(J1,E1)ΥK(J,K)K1)   -ϕ(ΥE(J1,E1)+ΥK(J1,K1))K1(ηK+ϕL1)(ξ-βL1)(ηL+βK)(L-L1)2   =ηJJ1(1-ΥE(J1,E1)ΥE(J,E1))(1-JJ1)   +ΥE(J1,E1)(ΥE(J,E)ΥE(J,E1)-EE1)(1-ΥE(J,E1)ΥE(J,E))   +ΥE(J1,E1)(5-ΥE(J1,E1)ΥE(J,E1)-ΥE(J,E)P1ΥE(J1,E1)P   -PK1P1K-KE1K1E-ΥE(J,E1)EΥE(J,E)E1)   +ΥP(J1,P1)(ΥE(J1,E1)ΥP(J,P)ΥE(J,E1)ΥP(J1,P1)-PP1)   ×(1-ΥE(J,E1)ΥP(J1,P1)ΥE(J1,E1)ΥP(J,P))   +ΥP(J1,P1)(3-ΥE(J1,E1)ΥE(J,E1)-ΥP(J,P)P1ΥP(J1,P1)P   -ΥE(J,E1)ΥP(J1,P1)PΥE(J1,E1)ΥP(J,P)P1)   +ΥK(J1,K1)(ΥE(J1,E1)ΥK(J,K)ΥE(J,E1)ΥK(J1,K1)-KK1)   ×(1-ΥE(J,E1)ΥK(J1,K1)ΥE(J1,E1)ΥK(J,K))   +ΥK(J1,K1)(4-ΥE(J1,E1)ΥE(J,E1)-ΥK(J,K)P1ΥK(J1,K1)P   -PK1P1K-ΥE(J,E1)ΥK(J1,K1)KΥE(J1,E1)ΥK(J,K)K1)   -ϕ(ΥE(J1,E1)+ΥK(J1,K1))K1(ηK+ϕL1)(ξ-βL1)(ηL+βK)(L-L1)2.

From Condition C2, we have

(1-ΥE(J1,E1)ΥE(J,E1))(1-JJ1)0.

Under the AM-GM inequality, we discern that

5ΥE(J1,E1)ΥE(J,E1)+ΥE(J,E)P1ΥE(J1,E1)P+PK1P1K+KE1K1E   +ΥE(J,E1)EΥE(J,E)E1,4ΥE(J1,E1)ΥE(J,E1)+ΥK(J,K)P1ΥK(J1,K1)P+PK1P1K   +ΥE(J,E1)ΥK(J1,K1)KΥE(J1,E1)ΥK(J,K)K1,3ΥE(J1,E1)ΥE(J,E1)+ΥP(J,P)P1ΥP(J1,P1)P+ΥE(J,E1)ΥP(J1,P1)PΥE(J1,E1)ΥP(J,P)P1.

Moreover, from Equations (1820), we conclude

(ΥE(J,E)ΥE(J,E1)-EE1)(1-ΥE(J,E1)ΥE(J,E))0,(ΥE(J1,E1)ΥP(J,P)ΥE(J,E1)ΥP(J1,P1)-PP1)×(1-ΥE(J,E1)ΥP(J1,P1)ΥE(J1,E1)ΥP(J,P))0,(ΥE(J1,E1)ΥK(J,K)ΥE(J,E1)ΥK(J1,K1)-KK1)×(1-ΥE(J,E1)ΥK(J1,K1)ΥE(J1,E1)ΥK(J,K))0.

Therefore, when ℜ0 > 1, it implies that dΞ2dt0 for all J,P,K,E,L>0. Moreover, dΞ2dt=0 when J=J1,P=P1,K=K1,E=E1, and L=L1. Consequently, M2 ={O1}. Applying the L.I.P allows us to infer that if ℜ0 > 1, the steady state O1 is being G.A.S [44, 45].

3 HIV-1 model with distributed delays

According to estimates, it takes HIV-1 0.9 days to penetrate a target cell and start creating new HIV-1 particles. Perelson et al. [46]. Numerous studies have delved into models of HIV-1 infection, considering time delays and dual infection routes (see e.g., [4751]). This section enhances the previously introduced model by integrating three varieties of distributed time delays.

3.1 System overview

The system of delay differential equations (DDEs) presented below will be studied.

{J˙=θ-ηJJ-ΥE(J,E)-ΥP(J,P)-ΥK(J,K),P˙=0f1π1(r)e-b1r(ΥE(J(t-r),E(t-r))       +ΥP(J(t-r),P(t-r))+ΥK(J(t-r),      K(t-r)))dr-(ψ+ηP)P,K˙=ψ0f2π2(r)e-b2rP(t-r)dr-ηKK-ϕLK,E˙=ϖ0f3π3(r)e-b3rK(t-r)dr-ηEE,L˙=ξK-ηLL-βLK.    (23)

Here, system (23) includes the following assumptions:

• At the moment t, healthy cells are encountering either HIV-1 particles or infected cells transition to a latent infection state r units of time later. The emergence of latently infected cells at time t depends on the count of cells recently contacted at time tr, which persists until time t.

• After a latency period of r units of time after infection, cells that were initially latent and became actively infected. The emergence of actively infected cells at the time t depends on the count of cells recently transitioned into latent infection at time tr, which persist until time t.

• Newly formed mature HIV-1 particles emerge from actively infected cells r units of time after infection. The occurrence of HIV-1 particle production at time t relies on the count of cells that recently transitioned into actively infected status at time tr that remain viable until time t.

As such, the likelihood of surviving the time frame from tr to t is represented by πi(r)e-bir, where bi are positive constants and i = 1, 2, 3. Besides, we choose the delay parameter, r, at random. It is drawn from a probability distribution function πi(r), which spans the interval [0, fi], in which fi refers to the maximum delay duration. The functions πi(r), for i = 1, 2, 3 adhere to and fulfill the following specified conditions:

πi(r)>0,   0fiπi(r)dr=1,   and   0fiπi(r)e-μrdr<,where   μ>0.

Assuming that

Π̄i(r)=πi(r)e-bir, Πi=0fiΠ̄i(r)dr,i=1,2,3,

implies the fact that 0 < Πi ≤ 1, i = 1, 2, 3.

System (23) has the outlined initial conditions below:

{(J(v),P(v),K(v),E(v),L(v))=(k1(v),k2(v),k3(v),k4(v),k5(v)),ki(v)0,i=1,2,...,5,v[-f,0],f=max{f1,f2,f3},    (24)

where ki(v) ∈ C([−f, 0], ℝ≥0), i = 1, 2, ..., 5 and C= C([−f, 0], ℝ≥0) is the Banach space comprising continuous functions, and it is equipped with the norm ki=sup-fζ0|ki(ζ)| for all kiC. Thus, a unique solution for system (23) under the specified initial conditions (24) is guaranteed, as affirmed by references [44] and [52]. The meanings assigned to all remaining parameters and variables remain in line with the explanations provided earlier. We presume that functions ΥZ(J,Z), Z{E,P,K}, satisfy Conditions C1C6 as presented previously.

3.2 Fundamental characteristics

3.2.1 Model's well-posedness

Proposition 3.1. Assuming that Condition C1 is fulfilled, then the solutions of system (23) under initial conditions (24) are nonnegative and ultimately bounded.

Proof. Since J˙J=0=θ>0, then it can be inferred that J(t) is always positive for every t ≥ 0. Furthermore, other equations within system (23) can be expressed as:

P˙+(ψ+ηP)P=0f1Π̄1(r)(ΥE(J(t-r),E(t-r))                               +ΥP(J(t-r),P(t-r))                               +ΥK(J(t-r),K(t-r)))drP(t)=k2(0)e-(ψ+ηP)t+0te-(ψ+ηP)(t-λ)×0f1Π̄1(r)(ΥE(J(λ-r),E(λ-r))+ΥP(J(λ-r),P(λ-r))+ΥK(J(λ-r),K(λ-r))))drdλ0.K˙+(ηK+ϕL(t))K=ψ0f2Π̄2(r)P(t-r)drK(t)=k3(0)e-0t(ηK+ϕL(u))du+ψ0te-λt(ηK+ϕL(u))du×0f2Π̄2(r)P(λ-r)drdλ0.E˙+ηEE=ϖ0f3Π̄3(r)K(t-r)drE(t)=k4(0)e-ηEt+ϖ0te-ηE(t-λ)×0f3Π̄3(r)K(λ-r)drdλ0.L˙+(ηL+βK)L=ξKL(t)=k5(0)e-0t(ηL+βK(u))du+ξ0te-λt(ηL+βK(u))duK(r)dr0,

for every t ∈ [0, f]. By employing a recursive argumentation approach, it can be shown that J(t),P(t),K(t),E(t), and L(t) remain nonnegative for every t ≥ 0. As a result, system (23) only admits solutions where (J(t),P(t),K(t),E(t),L(t))05, for every t ≥ 0.

Obviously, limtsupJ(t)θηJ, as deduced from the first equation in system (23). Following that, we proceed with defining Φ1 as follows:

Φ1=0f1Π̄1(r)J(t-r)dr+P.

Therefore

Φ˙1=0f1Π̄1(r)dJ(t-r)dtdr+P˙   =0f1Π̄1(r)(θ-ηJJ(t-r))dr-(ψ+ηP)P   =θΠ1-ηJ0f1Π̄1(r)J(t-r)dr-(ψ+ηP)P   θ-ε1(0f1Π̄1(r)J(t-r)dr+P)=θ-ε1Φ1,

where ε1=min{ηJ,ψ+ηP}. This indicates that limtsupΦ1(t)θε1=Γ^1. Based on the nonnegativity of 0f1Π̄1(r)J(t-r)dr and P, then limtsupP(t)Γ^1 is confirmed. Furthermore, we let

Φ2=K+ηK2ξL.

This produces

Φ˙2=K˙+ηK2ξL˙   =ψ0f2Π̄2(r)P(t-r)dr-ηKK-ϕLK   +ηK2ξ(ξK-ηLL-βLK)   =ψ0f2Π̄2(r)P(t-r)dr-ηK2K   -ηKηL2ξL-(ϕ+ηKβ2ξ)LK   ψ0f2Π̄2(r)P(t-r)dr-ηK2K-ηKηL2ξL   ψΓ^1-ε2(K+ηK2ξL)=ψΓ^1-ε2Φ2,

where ε2=min{ηK2,ηL}. Thus, limtsupΦ2(t)ψΓ^1ε2=Γ^2. Based on the nonnegativity of K and L, one can guarantee that limtsupK(t)Γ^2, and limtsupL(t)2ξΓ^2ηK=Γ^3. As a result of the fourth equation in system (23), one can acquire

E˙=ϖ0f3Π̄3(r)K(t-r)dr-ηEEϖΠ3Γ^2-ηEEϖΓ^2-ηEE.

Then, limtsupE(t)ϖΓ^2ηE=Γ^4. Consequently, it can be inferred that J(t), P(t), K(t), E(t), and L(t) are being ultimately bounded.

3.2.2 Derivation of reproduction ratios and steady states

Proposition 3.2. Assuming that Conditions C1C4 are met, there exists a positive basic reproduction ratio ~0=Π1(ϖψΠ2Π3E(J~0)+ηEηKP(J~0)+ψηEΠ2K(J~0))ηEηK(ψ+ηP) for system (23) in a way that

1. ensures the system consistently maintains a virus-free steady state, denoted as Õ0, and

2. if ~0>1, the system additionally possesses a virus-persistence steady state, denoted as Õ1.

Proof. It is observed that a virus-free steady state labeled as Õ0=(J~0,0,0,0,0) consistently exists in system (23), where J~0=θηJ. Using NGM, the matrices below characterize the nonlinear terms, denoted as Ω^2 and responsible for new infections, as well as Λ^2, representing the outflow.

Ω^2=(Π1(ΥE(J,E)+ΥP(J,P)+ΥK(J,K))00),Λ^2=((ψ+ηP)P-ψΠ2P+ηKK+ϕLK-ϖΠ3K+ηEE).

We compute the derivatives of Ω^2 and Λ^2 at the steady state Õ0, leading to the generation of the subsequent matrices:

Ω2=(Π1ΥP(J~0,0)PΠ1ΥK(J~0,0)KΠ1ΥE(J~0,0)E000000),Λ2=(ψ+ηP00-ψΠ2ηK00-ϖΠ3ηE).

The form of the NGM is as follows:

Ω2Λ2-1=(ϖψΠ1Π2Π3E(J~0)ηEηK(ψ+ηP)+Π1P(J~0)ψ+ηP+ψΠ1Π2K(J~0)ηK(ψ+ηP)ϖΠ1Π3E(J~0)ηEηK+Π1K(J~0)ηKΠ1E(J~0)ηE000000).

The expression representing the basic reproduction ratio ~0 is outlined as follows:

~0=Π1(ϖψΠ2Π3E(J~0)+ηEηKP(J~0)+ψηEΠ2K(J~0))ηEηK(ψ+ηP)=~0E+~0P+~0K,    (25)

where

~0E=ϖψΠ1Π2Π3E(J~0)ηEηK(ψ+ηP),   ~0P=Π1P(J~0)ψ+ηP, ~0K=ψΠ1Π2K(J~0)ηK(ψ+ηP).

The parameters ~0Z, Z{E,P,K} share the identical biological meaning as those elucidated in Section 2. To extend our exploration beyond Õ0, we scrutinize (J,P,K,E,L) as a potential steady state governed by the following set of algebraic equations:

0=θ-ηJJ-ΥE(J,E)-ΥP(J,P)-ΥK(J,K),    (26)
0=Π1(ΥE(J,E)+ΥP(J,P)+ΥK(J,K))-(ψ+ηP)P,    (27)
0=ψΠ2P-ηKK-ϕLK,    (28)
0=ϖΠ3K-ηEE,    (29)
0=ξK-ηLL-βLK.    (30)

From Equations (29, 30), we get

E=Θ~1(K)=ϖΠ3KηE, L=ξKηL+βK.    (31)

Clearly, Θ~1(0)=0. The following outcome is obtained after inserting Equation (31) into Equation (28):

P=Θ~2(K)=ηLηKK+(ϕξ+ηKβ)K2ψΠ2(ηL+βK).    (32)

It is clear that Θ~2(0)=0. From Equations (26, 27), we get

θ-ηJJ=ψ+ηPΠ1P.    (33)

The following outcome is obtained after inserting Equation (32) into Equation (33):

J=Θ~3(K) =1ηJ(θ-(ψ+ηP)(ηLηKK+(ϕξ+ηKβ)K2)ψΠ1Π2(ηL+βK)).       (34)

Note that Θ~3(0)=J~0. Upon substitution of Equations (31, 32, 34) into Equation (27), the result is as follows:

ΥE(Θ~3(K),Θ~1(K))+ΥP(Θ~3(K),Θ~2(K))+ΥK(Θ~3(K),K)-(ψ+ηP)(ηLηKK+(ϕξ+ηKβ)K2)ψΠ1Π2(ηL+βK)=0.    (35)

From Equation (35), we have

1. When K=0, the virus-free steady state Õ0 is derived from Equations (31, 32, 33, 34).

2. When K0, let us define a function Ψ~(K) on [0, ∞) as:

Ψ~(K)=ΥE(Θ~3(K),Θ~1(K))+ΥP(Θ~3(K),Θ~2(K))   +ΥK(Θ~3(K),K)   -(ψ+ηP)(ηLηKK+(ϕξ+ηKβ)K2)ψΠ1Π2(ηL+βK).

We have Ψ~(0)=0. Let Ǩ be such that Θ~3(Ǩ)=0, i.e.,

J~0-(ψ+ηP)(ηLηKǨ+(ϕξ+ηKβ)Ǩ2)ηJψΠ1Π2(ηL+βǨ)=0,

which implies that

(ϕξ+ηKβ)(ψ+ηP)Ǩ2   +(ηLηK(ψ+ηP)-ηJψβΠ1Π2J~0)Ǩ-ηJψηLΠ1Π2J~0=0.    (36)

Therefore, Equation (36) has a positive solution

Ǩ=-B~+B~2-4ÃC~2Ã,

where

Ã=(ϕξ+ηKβ)(ψ+ηP),   B~=(ηLηK(ψ+ηP)-ηJψβΠ1Π2J~0),   C~=-ηJψηLΠ1Π2J~0.

We can see that

Ψ~(Ǩ)=ΥE(0,Θ~1(Ǩ))+ΥP(0,Θ~2(Ǩ))+ΥK(0,Ǩ)      -(ψ+ηP)(ηLηKǨ+(ϕξ+ηKβ)Ǩ2)ψΠ1Π2(ηL+βǨ)      =-(ψ+ηP)(ηLηKǨ+(ϕξ+ηKβ)Ǩ2)ψΠ1Π2(ηL+βǨ)<0.

In addition

Ψ~(K)=Θ~3(K)ΥE(J,E)J+Θ~1(K)ΥE(J,E)E      +Θ~3(K)ΥP(J,P)J+Θ~2(K)ΥP(J,P)P      +Θ~3(K)ΥK(J,K)J+ΥK(J,K)K      -ψ+ηPψΠ1Π2(ηK+ϕξ(2ηL+βK)K(ηL+βK)2).

Condition C1 implies that ΥZ(J~0,0)J=0, Z{E,P,K}. Then

Ψ~(0)=Θ~1(0)ΥE(J~0,0)E+Θ~2(0)ΥP(J~0,0)P  +ΥK(J~0,0)K-ηK(ψ+ηP)ψΠ1Π2  =ϖΠ3ηEΥE(J~0,0)E+ηKψΠ2ΥP(J~0,0)P  +ΥK(J~0,0)K-ηK(ψ+ηP)ψΠ1Π2  =ϖΠ3E(J~0)ηE+ηKP(J~0)ψΠ2+K(J~0)  -ηK(ψ+ηP)ψΠ1Π2  =ηK(ψ+ηP)ψΠ1Π2(~0-1),

where ~0 is defined in Equation (25). Therefore, if ~0>1 then Ψ~(0)>0 and there exists K~1(0,Ǩ) such that Ψ~(K~1)=0. Let K=K~1 in Equation (26) and define

~(J)=θ-ηJJ-ΥE(J,Θ~1(K~1))-ΥP(J,Θ~2(K~1))      -ΥK(J,K~1).

Subsequently, based on Condition C1, we obtain ~(0)=θ>0 and

~(J~0)=-(ΥE(J~0,Θ~1(K~1))+ΥP(J~0,Θ~2(K~1))      +ΥK(J~0,K~1))<0.

Under Condition C2, it follows that ~(J) decreases strictly as a function of J. Consequently, there is a unique value J~1 within the interval (0,J~0) for which ~(J~1) equals zero. Moreover, considering Equations (31, 32), we find

P~1=ηLηKK~1+(ϕξ+ηKβ)K~12ψΠ2(ηL+βK~1), E˜1=ϖΠ3K~1ηE,      L~1=ξK~1ηL+βK~1.

The presence of the virus-persistence steady state Õ1=(J~1,P~1,K~1,E˜1,L~1) becomes evident when ~0 being > 1.

3.2.3 Global stability of steady states

The forthcoming theorems explore the global asymptotic stability of both virus-free and virus-persistence steady states. For clarity's sake, let's demonstrate (J(t-r),P(t-r),P(t-r),E(t-r)) by (Jr,Pr,Pr,Er).

Theorem 3.1. For system (23), let ~01 and Conditions C1-C5 satisfied, then Õ0 is G.A.S.

Proof. Let's contemplate a potential Lyapunov function, Ξ3(J,P,K,E,L), in the following manner:

Ξ3=J-J~0-J~0JE(J~0)E(λ)dλ+1Π1P      +ϖΠ3E(J~0)+ηEK(J~0)ηKηEK+E(J~0)ηEE      +ϕ(ϖΠ3E(J~0)+ηEK(J~0))2ξηKηEL2      +1Π10f1Π̄1(r)t-rt(ΥE(J(λ),E(λ))      +ΥP(J(λ),P(λ))+ΥK(J(λ),K(λ)))dλdr      +ψ(ϖΠ3E(J~0)+ηEK(J~0))ηKηE0f2Π̄2(r)      ×t-rtP(λ)dλdr+ϖE(J~0)ηE0f3Π̄3(r)t-rtK(λ)dλdr.

Evidently, Ξ3(J,P,K,E,L)>0 for every J,P,K,E,L>0, as well as Ξ3(J~0,0,0,0,0)=0. Further, dΞ3dt is given by:

dΞ3dt=(1-E(J~0)E(J))J˙+1Π1P˙      +ϖΠ3E(J~0)+ηEK(J~0)ηKηEK˙+E(J~0)ηEE˙      +ϕ(ϖΠ3E(J~0)+ηEK(J~0))LξηKηEL˙      +ΥE(J,E)+ΥP(J,P)+ΥK(J,K)      -1Π10f1Π̄1(r)(ΥE(Jr,Er)+ΥP(Jr,Pr)      +ΥK(Jr,Kr))dr      +ψΠ2(ϖΠ3E(J~0)+ηEK(J~0))ηKηEP      -ψ(ϖΠ3E(J~0)+ηEK(J~0))ηKηE0f2Π̄2(r)Prdr      +ϖΠ3E(J~0)ηEK      -ϖE(J~0)ηE0f3Π̄3(r)Krdr      =(1-E(J~0)E(J))(θ-ηJJ-ΥE(J,E)      -ΥP(J,P)-ΥK(J,K))      +1Π10f1Π̄1(r)(ΥE(Jr,Er)+ΥP(Jr,Pr)      +ΥK(Jr,Kr))dr      -ψ+ηPΠ1P+ϖΠ3E(J~0)+ηEK(J~0)ηKηE      ×(ψ0f2Π̄2(r)Prdr-ηKK-ϕLK)      +E(J~0)ηE(ϖ0f3Π̄3(r)Krdr-ηEE)      +ϕ(ϖΠ3E(J~0)+ηEK(J~0))LξηKηE(ξK      -ηLL-βLK)+ΥE(J,E)+ΥP(J,P)+ΥK(J,K)      -1Π10f1Π̄1(r)(ΥE(Jr,Er)+ΥP(Jr,Pr)      +ΥK(Jr,Kr))dr      +ψΠ2(ϖΠ3E(J~0)+ηEK(J~0))ηKηEP      -ψ(ϖΠ3E(J~0)+ηEK(J~0))ηKηE0f2Π̄2(r)Prdr      +ϖΠ3E(J~0)ηEK-ϖE(J~0)ηE0f3Π̄3(r)Krdr      =(1-E(J~0)E(J))(θ-ηJJ)      +(E(J~0)E(J)ΥE(J,E)E-E(J~0))E      +(E(J~0)E(J)ΥP(J,P)P-ψ+ηPΠ1      +ψΠ2(ϖΠ3E(J~0)+ηEK(J~0))ηKηE)P      +(E(J~0)E(J)ΥK(J,K)K-K(J~0))K      -ϕ(ϖΠ3E(J~0)+ηEK(J~0))ξηKηE(ηL+βK)L2.    (37)

Using inequality (17), we obtain

E(J~0)E(J)ΥE(J,E)E-E(J~0)E(J~0)E(J)E(J)-E(J~0)=0.

Further, Condition C4 and C5 in case of J~0=θηJ imply that

E(J~0)E(J)ΥP(J,P)P-ψ+ηPΠ1+ψΠ2(ϖΠ3E(J~0)+ηEK(J~0))ηKηEE(J~0)P(J)E(J)-ψ+ηPΠ1+ψΠ2(ϖΠ3E(J~0)+ηEK(J~0))ηKηEE(J~0)P(J~0)E(J~0)-ψ+ηPΠ1+ψΠ2(ϖΠ3E(J~0)+ηEK(J~0))ηKηE=P(J~0)-ψ+ηPΠ1+ψΠ2(ϖΠ3E(J~0)+ηEK(J~0))ηKηE=ψΠ2(ϖΠ3E(J~0)+ηEK(J~0))+ηKηEP(J~0)ηKηE-ψ+ηPΠ1=ψ+ηPΠ1(~0-1).

Furthermore

E(J~0)E(J)ΥK(J,K)K-K(J~0)E(J~0)K(J)E(J)-K(J~0)      E(J~0)K(J~0)E(J~0)-K(J~0)=0.

Upon conducting a straightforward computation and incorporating the value J~0=θ/ηJ, Equation (37) will manifest in the subsequent manner:

dΞ3dtθ(1-E(J~0)E(J))(1-JJ~0)+ψ+ηPΠ1(~0-1)P      -ϕ(ϖΠ3E(J~0)+ηEK(J~0))ξηKηE(ηL+βK)L2.

From Condition C3, we have

(1-E(J~0)E(J))(1-JJ~0)0.

Evidently, for ~01, the rate dΞ3dt remains non-positive. Furthermore, dΞ3dt equals zero at J=J~0, with P=K=L=0. It follows that all solutions convey to M3 . Within M3 , every element satisfy J(t)=J~0 and P(t)=K(t)=L(t)=0 for all t. Following this, the first equation of model (23) along with Condition C1 results in:

0=J˙=θ-ηJJ~0-ΥE(J~0,E(t))ΥE(J~0,E(t))=0      E(t)=0,for all t.

Then, M3 ={Õ0}. Consequently, according to the L.I.P, the inference can be drawn that Õ0 exhibits global asymptotic stability whenever ~01 [44, 45].

Theorem 3.2. Given ~0>1 and the fulfillment of Conditions C1-C4 and C6, it can be concluded that the steady state Õ1 for system (23) is guaranteed to be G.A.S.

Proof. Building up Ξ4(J,P,K,E,L) according to Equation (38) as follows:

Ξ4=J-J~1-J~1JΥE(J~1,E˜1)ΥE(λ,E˜1)dλ+P~1Π1𝐹(PP~1)   +ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1)ηEE˜1(ηK+ϕL~1)𝐹(KK~1)   +ΥE(J~1,E˜1)ηE𝐹(EE˜1)   +ϕ(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))2ηEK~1E˜1(ηK+ϕL~1)(ξ-βL~1)(L-L~1)2   +ΥE(J~1,E˜1)Π10f1Π̄1(r)t-rt𝐹(ΥE(J(λ),E(λ))ΥE(J~1,E˜1))dλdr   +ΥP(J~1,P~1)Π10f1Π̄1(r)t-rt𝐹(ΥP(J(λ),P(λ))ΥP(J~1,P~1))dλdr   +ΥK(J~1,K~1)Π10f1Π̄1(r)t-rt𝐹(ΥK(J(λ),K(λ))ΥK(J~1,K~1))dλdr   +ψP~1(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))ηEK~1E˜1(ηK+ϕL~1)   0f2Π̄2(r)t-rt𝐹(P(λ)P~1)dλdr   +ϖK~1ΥE(J~1,E˜1)ηEE˜10f3Π̄3(r)t-rt𝐹(K(λ)K~1)dλdr.    (38)

We can deduce from the steady state condition Equation (30) that ξ-βL~1=ηLL~1K~1>0. The positivity definiteness of Ξ4 becomes evident. When calculating dΞ4dt along the trajectories of the model described in Equation (23), we obtain

dΞ4dt=(1-ΥE(J~1,E˜1)ΥE(J,E˜1))J˙+1Π1(1-P~1P)P˙+ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1)ηEK~1E˜1(ηK+ϕL~1)(1-K~1K)K˙+ΥE(J~1,E˜1)ηEE˜1(1-E˜1E)E˙+ϕ(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))ηEK~1E˜1(ηK+ϕL~1)(ξ-βL~1)×(L-L~1)L˙+ΥE(J~1,E˜1)Π10f1Π̄1(r)(ΥE(J,E)ΥE(J~1,E˜1)-ΥE(Jr,Er)ΥE(J~1,E˜1)+ln(ΥE(Jr,Er)ΥE(J,E)))dr+ΥP(J~1,P~1)Π10f1Π̄1(r)(ΥP(J,P)ΥP(J~1,P~1)-ΥP(Jr,Pr)ΥP(J~1,P~1)+ln(ΥP(Jr,Pr)ΥP(J,P)))dr+ΥK(J~1,K~1)Π10f1Π̄1(r)(ΥK(J,K)ΥK(J~1,K~1)-ΥK(Jr,Kr)ΥK(J~1,K~1)+ln(ΥK(Jr,Kr)ΥK(J,K)))dr+ψP~1(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))ηEK~1E˜1(ηK+ϕL~1)×0f2Π̄2(r)(PP~1-PrP~1+ln(PrP))dr+ϖK~1ΥE(J~1,E˜1)ηEE˜10f3Π̄3(r)(KK~1-KrK~1+ln(KrK))dr=(1-ΥE(J~1,E˜1)ΥE(J,E˜1))(θ-ηJJ-ΥE(J,E)-ΥP(J,P)-ΥK(J,K))+1Π1(1-P~1P)(0f1Π̄1(r)(ΥE(Jr,Er)+ΥP(Jr,Pr)+ΥK(Jr,Kr))dr-(ψ+ηP)P)+ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1)ηEK~1E˜1(ηK+ϕL~1)(1-K~1K)×(ψ0f2Π̄2(r)Prdr-ηKK-ϕLK)+ΥE(J~1,E˜1)ηEE˜1(1-E˜1E)(ϖ0f3Π̄3(r)Krdr-ηEE)+ϕ(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))ηEK~1E˜1(ηK+ϕL~1)(ξ-βL~1)×(L-L~1)(ξK-ηLL-βLK)+ΥE(J,E)-1Π10f1Π̄1(r)ΥE(Jr,Er)dr+ΥE(J~1,E˜1)Π10f1Π̄1(r)ln(ΥE(Jr,Er)ΥE(J,E))dr+ΥP(J,P)-1Π10f1Π̄1(r)ΥP(Jr,Pr)dr+ΥP(J~1,P~1)Π10f1Π̄1(r)ln(ΥP(Jr,Pr)ΥP(J,P))dr+ΥK(J,K)-1Π10f1Π̄1(r)ΥK(Jr,Kr)dr+ΥK(J~1,K~1)Π10f1Π̄1(r)ln(ΥK(Jr,Kr)ΥK(J,K))dr+ψΠ2P~1(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))ηEK~1E˜1(ηK+ϕL~1)PP~1-ψ(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))ηEK~1E˜1(ηK+ϕL~1)×0f2Π̄2(r)Prdr+ψP~1(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))ηEK~1E˜1(ηK+ϕL~1)×0f2Π̄2(r)ln(PrP)dr+ϖΠ3K~1ΥE(J~1,E˜1)ηEE˜1KK~1-ϖΥE(J~1,E˜1)ηEE˜1×0f3Π̄3(r)Krdr+ϖK~1ΥE(J~1,E˜1)ηEE˜10f3Π̄3(r)ln(KrK)dr.

Thus

dΞ4dt=(1-ΥE(J~1,E˜1)ΥE(J,E˜1))(θ-ηJJ)   +ΥE(J~1,E˜1)ΥE(J,E˜1)(ΥE(J,E)+ΥP(J,P)+ΥK(J,K))   -1Π10f1Π̄1(r)ΥE(Jr,Er)P~1Pdr   -1Π10f1Π̄1(r)ΥP(Jr,Pr)P~1Pdr   -1Π10f1Π̄1(r)ΥK(Jr,Kr)P~1Pdr   -ψ+ηPΠ1P+ψ+ηPΠ1P~1   -ψ(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))ηEK~1E˜1(ηK+ϕL~1)   ×0f2Π̄2(r)PrK~1Kdr   -ηK(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))ηEK~1E˜1(ηK+ϕL~1)(K-K~1)   -ϕ(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))ηEK~1E˜1(ηK+ϕL~1)L(K-K~1)   -ϖΥE(J~1,E˜1)ηE0f3Π̄3(r)KrEdr-ΥE(J~1,E˜1)EE˜1   +ΥE(J~1,E˜1)   +ϕ(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))ηEK~1E˜1(ηK+ϕL~1)(ξ-βL~1)   ×(L-L~1)(ξK-ηLL-βLK)   +ΥE(J~1,E˜1)Π10f1Π̄1(r)ln(ΥE(Jr,Er)ΥE(J,E))dr   +ΥP(J~1,P~1)Π10f1Π̄1(r)ln(ΥP(Jr,Pr)ΥP(J,P))dr   +ΥK(J~1,K~1)Π10f1Π̄1(r)ln(ΥK(Jr,Kr)ΥK(J,K))dr   +ψΠ2P~1(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))ηEK~1E˜1(ηK+ϕL~1)PP~1   +ψP~1(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))ηEK~1E˜1(ηK+ϕL~1)   ×0f2Π̄2(r)ln(PrP)dr+ϖΠ3K~1ΥE(J~1,E˜1)ηEE˜1KK~1    +ϖK~1ΥE(J~1,E˜1)ηEE˜10f3Π̄3(r)ln(KrK)dr.    (39)

Implementing the specified conditions for steady state Õ1

θ=ηJJ~1+ΥE(J~1,E˜1)+ΥP(J~1,P~1)+ΥK(J~1,K~1),ΥE(J~1,E˜1)+ΥP(J~1,P~1)+ΥK(J~1,K~1)=(ψ+ηP)P~1Π1,K~1=ψΠ2P~1ηK+ϕL~1,   E˜1=ϖΠ3K~1ηE,   ξK~1=ηLL~1+βL~1K~1,

we derive

ΥE(J~1,E˜1)+ΥK(J~1,K~1)=ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1)ηEE˜1, =ψΠ2P~1(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))ηEK~1E˜1(ηK+ϕL~1).

Therefore, Equation (39) will be expressed as follows:

dΞ4dt=(1-ΥE(J~1,E˜1)ΥE(J,E˜1))(ηJJ~1-ηJJ)   +(2-ΥE(J~1,E˜1)ΥE(J,E˜1))(ΥE(J~1,E˜1)   +ΥP(J~1,P~1)+ΥK(J~1,K~1))+ΥE(J~1,E˜1)   ×ΥE(J,E)ΥE(J,E˜1)   +ΥP(J~1,P~1)ΥE(J~1,E˜1)ΥP(J,P)ΥE(J,E˜1)ΥP(J~1,P~1)   +ΥK(J~1,K~1)ΥE(J~1,E˜1)ΥK(J,K)ΥE(J,E˜1)ΥK(J~1,K~1)   -ΥE(J~1,E˜1)Π10f1Π̄1(r)ΥE(Jr,Er)P~1ΥE(J~1,E˜1)Pdr   -ΥP(J~1,P~1)Π10f1Π̄1(r)ΥP(Jr,Pr)P~1ΥP(J~1,P~1)Pdr   -ΥK(J~1,K~1)Π10f1Π̄1(r)ΥK(Jr,Kr)P~1ΥK(J~1,K~1)Pdr   -(ψ+ηP)P~1Π1PP~1   -ψP~1(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))ηEK~1E˜1(ηK+ϕL~1)   ×0f2Π̄2(r)PrK~1P~1Kdr   -ηK(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))ηEK~1E˜1(ηK+ϕL~1)(K-K~1)   -ϕ(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))ηEK~1E˜1(ηK+ϕL~1)L(K-K~1)   +ϕ(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))ηEK~1E˜1(ηK+ϕL~1)L~1(K-K~1)   -ϕ(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))ηEK~1E˜1(ηK+ϕL~1)L~1(K-K~1)   -ϖK~1ΥE(J~1,E˜1)ηE0f3Π̄3(r)KrK~1Edr-ΥE(J~1,E˜1)EE˜1   +ΥE(J~1,E˜1)   +ϕ(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))ηEK~1E˜1(ηK+ϕL~1)(ξ-βL~1)(L-L~1)   ×(ξK-ηLL-βLK-ξK~1+ηLL~1+βL~1K~1   -βL~1K+βL~1K)   +ΥE(J~1,E˜1)Π10f1Π̄1(r)ln(ΥE(Jr,Er)ΥE(J,E))dr   +ΥP(J~1,P~1)Π10f1Π̄1(r)ln(ΥP(Jr,Pr)ΥP(J,P))dr   +ΥK(J~1,K~1)Π10f1Π̄1(r)ln(ΥK(Jr,Kr)ΥK(J,K))dr   +(ΥE(J~1,E˜1)+ΥK(J~1,K~1))PP~1   +ΥE(J~1,E˜1)+ΥK(J~1,K~1)Π20f2Π̄2(r)ln(PrP)dr   +ΥE(J~1,E˜1)KK~1+ΥE(J~1,E˜1)Π30f3Π̄3(r)ln(KrK)dr.

Consequently, the form of the expression becomes as follows:

dΞ4dt=ηJJ~1(1-ΥE(J~1,E˜1)ΥE(J,E˜1))(1-JJ~1)   +(2-ΥE(J~1,E˜1)ΥE(J,E˜1))(ΥE(J~1,E˜1)+ΥP(J~1,P~1)   +ΥK(J~1,K~1))+ΥE(J~1,E˜1)ΥE(J,E)ΥE(J,E˜1)   +ΥP(J~1,P~1)ΥE(J~1,E˜1)ΥP(J,P)ΥE(J,E˜1)ΥP(J~1,P~1)   +ΥK(J~1,K~1)ΥE(J~1,E˜1)ΥK(J,K)ΥE(J,E˜1)ΥK(J~1,K~1)   -ΥE(J~1,E˜1)Π10f1Π̄1(r)ΥE(Jr,Er)P~1ΥE(J~1,E˜1)Pdr   -ΥP(J~1,P~1)Π10f1Π̄1(r)ΥP(Jr,Pr)P~1ΥP(J~1,P~1)Pdr   -ΥK(J~1,K~1)Π10f1Π̄1(r)ΥK(Jr,Kr)P~1ΥK(J~1,K~1)Pdr   -(ΥE(J~1,E˜1)+ΥP(J~1,P~1)   +ΥK(J~1,K~1))PP~1-ΥE(J~1,E˜1)+ΥK(J~1,K~1)Π2   ×0f2Π̄2(r)PrK~1P~1Kdr   -ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1)ηEK~1E˜1(ηK+ϕL~1)   ×(ηK+ϕL~1)(K-K~1)   -ϕ(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))ηEK~1E˜1(ηK+ϕL~1)   ×(L-L~1)(K-K~1)   -ΥE(J~1,E˜1)Π30f3Π̄3(r)KrE˜1K~1Edr-ΥE(J~1,E˜1)EE˜1   +ΥE(J~1,E˜1)   +ϕ(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))ηEK~1E˜1(ηK+ϕL~1)(ξ-βL~1)   ×(ξ-βL~1)(L-L~1)(K-K~1)   -ϕ(ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1))(ηL+βK)ηEK~1E˜1(ηK+ϕL~1)(ξ-βL~1)   ×(L-L~1)2+ΥE(J~1,E˜1)Π10f1Π̄1(r)ln(ΥE(Jr,Er)ΥE(J,E))dr   +ΥP(J~1,P~1)Π10f1Π̄1(r)ln(ΥP(Jr,Pr)ΥP(J,P))dr   +ΥK(J~1,K~1)Π10f1Π̄1(r)ln(ΥK(Jr,Kr)ΥK(J,K))dr   +(ΥE(J~1,E˜1)   +ΥK(J~1,K~1))PP~1+ΥE(J~1,E˜1)+ΥK(J~1,K~1)Π2   ×0f2Π̄2(r)ln(PrP)dr   +ΥE(J~1,E˜1)KK~1+ΥE(J~1,E˜1)Π30f3Π̄3(r)ln(KrK)dr   =ηJJ~1(1-ΥE(J~1,E˜1)ΥE(J,E˜1))(1-JJ~1)   +(2-ΥE(J~1,E˜1)ΥE(J,E˜1))(ΥE(J~1,E˜1)+ΥP(J~1,P~1)   +ΥK(J~1,K~1))+ΥE(J~1,E˜1)ΥE(J,E)ΥE(J,E˜1)   +ΥP(J~1,P~1)ΥE(J~1,E˜1)ΥP(J,P)ΥE(J,E˜1)ΥP(J~1,P~1)   +ΥK(J~1,K~1)ΥE(J~1,E˜1)ΥK(J,K)ΥE(J,E˜1)ΥK(J~1,K~1)   -ΥE(J~1,E˜1)Π10f1Π̄1(r)ΥE(Jr,Er)P~1ΥE(J~1,E˜1)Pdr   -ΥP(J~1,P~1)Π10f1Π̄1(r)ΥP(Jr,Pr)P~1ΥP(J~1,P~1)Pdr   -ΥK(J~1,K~1)Π10f1Π̄1(r)ΥK(Jr,Kr)P~1ΥK(J~1,K~1)Pdr   -ΥP(J~1,P~1)PP~1   -ΥE(J~1,E˜1)+ΥK(J~1,K~1)Π20f2Π̄2(r)PrK~1P~1Kdr   -ϖΠ3K~1ΥE(J~1,E˜1)+ηEE˜1ΥK(J~1,K~1)ηEE˜1(KK~1-1)   -ΥE(J~1,E˜1)Π30f3Π̄3(r)KrE˜1K~1Edr-ΥE(J~1,E˜1)EE˜1   +ΥE(J~1,E˜1)   -ϕ(ΥE(J~1,E˜1)+ΥK(J~1,K~1))K~1(ηK+ϕL~1)(ξ-βL~1)(ηL+βK)(L-L~1)2   +ΥE(J~1,E˜1)Π10f1Π̄1(r)ln(ΥE(Jr,Er)ΥE(J,E))dr   +ΥP(J~1,P~1)Π10f1Π̄1(r)ln(ΥP(Jr,Pr)ΥP(J,P))dr   +ΥK(J~1,K~1)Π10f1Π̄1(r)ln(ΥK(Jr,Kr)ΥK(J,K))dr   +ΥE(J~1,E˜1)+ΥK(J~1,K~1)Π20f2Π̄2(r)ln(PrP)dr   +ΥE(J~1,E˜1)KK~1+ΥE(J~1,E˜1)Π30f3Π̄3(r)ln(KrK)dr.

Thus

dΞ4dt=ηJJ˜1(1Υ(J˜1,˜1)Υ(J,˜1))(1JJ˜1)     +(2Υ(J˜1,˜1)Υ(J,˜1))(Υ(J˜1,˜1)+ΥP(J˜1,P˜1)     +ΥK(J˜1,K˜1))+Υ(J˜1,˜1)Υ(J,)Υ(J,˜1)     +ΥP(J˜1,P˜1)Υ(J˜1,˜1)ΥP(J,P)Υ(J,˜1)ΥP(J˜1,P˜1)     +ΥK(J˜1,K˜1)Υ(J˜1,˜1)ΥK(J,K)Υ(J,˜1)ΥK(J˜1,K˜1)     Υ(J˜1,˜1)π10f1π¯1(r)Υ(Jr,r)P˜1Υ(J˜1,˜1)Pdr     ΥP(J˜1,P˜1)π10f1π¯1(r)ΥP(Jr,Pr)P˜1ΥP(J˜1,P˜1)Pdr     ΥK(J˜1,K˜1)π10f1π¯1(r)ΥK(Jr,Kr)P˜1ΥK(J˜1,K˜1)Pdr     ΥP(J˜1,P˜1)PP˜1     Υ(J˜1,˜1)+ΥK(J˜1,K˜1)π20f2π¯2(r)PrK˜1P˜1Kdr     (Υ(J˜1,˜1)+ΥK(J˜1,K˜1))(KK˜11)     Υ(J˜1,˜1)π30f3π¯3(r)Kr˜1K˜1dr     Υ(J˜1,˜1)˜1     +Υ(J˜1,˜1)ϕ(Υ(J˜1,˜1)+ΥK(J˜1,K˜1))K˜1(ηK+ϕ˜1)(ξβ˜1)(η+βK)(˜1)2     +Υ(J˜1,˜1)π10f1π¯1(r)ln(Υ(Jr,r)Υ(J,))dr     +ΥP(J˜1,P˜1)π10f1π¯1(r)ln(ΥP(Jr,Pr)ΥP(J,P))dr     +ΥK(J˜1,K˜1)π10f1π¯1(r)ln(ΥK(Jr,Kr)ΥK(J,K))dr     +Υ(J˜1,˜1)+ΥK(J˜1,K˜1)π20f2π¯2(r)ln(PrP)dr     +Υ(J˜1,˜1)KK˜1+Υ(J˜1,˜1)π30f3π¯3(r)ln(KrK)dr     =ηJJ˜1(1Υ(J˜1,˜1)Υ(J,˜1))(1JJ˜1)     +(2Υ(J˜1,˜1)Υ(J,˜1))(Υ(J˜1,˜1)+ΥP(J˜1,P˜1)     +ΥK(J˜1,K˜1))+Υ(J˜1,˜1)Υ(J,)Υ(J,˜1)     +ΥP(J˜1,P˜1)Υ(J˜1,˜1)ΥP(J,P)Υ(J,˜1)ΥP(J˜1,P˜1)     +ΥK(J˜1,K˜1)Υ(J˜1,˜1)ΥK(J,K)Υ(J,˜1)ΥK(J˜1,K˜1)     Υ(J˜1,˜1)π10f1π¯1(r)Υ(Jr,r)P˜1Υ(J˜1,˜1)Pdr     ΥP(J˜1,P˜1)π10f1π¯1(r)ΥP(Jr,Pr)P˜1ΥP(J˜1,P˜1)Pdr     ΥK(J˜1,K˜1)π10f1π¯1(r)ΥK(Jr,Kr)P˜1ΥK(J˜1,K˜1)Pdr     ΥP(J˜1,P˜1)PP˜1     Υ(J˜1,˜1)+ΥK(J˜1,K˜1)π20f2π¯2(r)PrK˜1P˜1Kdr     ΥK(J˜1,K˜1)KK˜1+ΥK(J˜1,K˜1)     Υ(J˜1,˜1)π30f3π¯3(r)Kr˜1K˜1drΥ(J˜1,˜1)˜1     +2Υ(J˜1,˜1)     ϕ(Υ(J˜1,˜1)+ΥK(J˜1,K˜1))K˜1(ηK+ϕ˜1)(ξβ˜1)(η+βK)(˜1)2     +Υ(J˜1,˜1)π10f1π¯1(r)ln(Υ(Jr,r)Υ(J,))dr     +ΥP(J˜1,P˜1)π10f1π¯1(r)ln(ΥP(Jr,Pr)ΥP(J,P))dr     +ΥK(J˜1,K˜1)π10f1π¯1(r)ln(ΥK(Jr,Kr)ΥK(J,K))dr     +Υ(J˜1,˜1)+ΥK(J˜1,K˜1)π20f2π¯2(r)ln(PrP)dr     +Υ(J˜1,˜1)π30f3π¯3(r)ln(KrK)dr.

Additionally, we have the following equalities:

ln(ΥE(Jr,Er)ΥE(J,E))=ln(ΥE(Jr,Er)P~1ΥE(J~1,E˜1)P)+ln(ΥE(J~1,E˜1)ΥE(J,E˜1))+ln(ΥE(J,E˜1)EΥE(J,E)E˜1)+ln(PK~1P~1K)+ln(KE˜1K~1E),ln(ΥP(Jr,Pr)ΥP(J,P))=ln(ΥP(Jr,Pr)P~1ΥP(J~1,P~1)P)+ln(ΥE(J~1,E˜1)ΥE(J,E˜1))+ln(ΥE(J,E˜1)ΥP(J~1,P~1)PΥE(J~1,E˜1)ΥP(J,P)P~1),ln(ΥK(Jr,Kr)ΥK(J,K))=ln(ΥK(Jr,Kr)P~1ΥK(J~1,K~1)P)+ln(ΥE(J~1,E˜1)ΥE(J,E˜1))+ln(ΥE(J,E˜1)ΥK(J~1,K~1)KΥE(J~1,E˜1)ΥK(J,K)K~1)+ln(PK~1P~1K),            ln(PrP)=ln(PrK~1P~1K)+ln(P~1KPK~1),ln(KrK)=ln(KrE˜1K~1E)+ln(K~1EKE˜1).

Therefore, dΞ4dt will be

dΞ4dt=ηJJ~1(1-ΥE(J~1,E˜1)ΥE(J,E˜1))(1-JJ~1)   +ΥE(J~1,E˜1)(ΥE(J,E)ΥE(J,E˜1)-EE˜1)   +ΥP(J~1,P~1)(ΥE(J~1,E˜1)ΥP(J,P)ΥE(J,E˜1)ΥP(J~1,P~1)-PP~1)   +ΥK(J~1,K~1)(ΥE(J~1,E˜1)ΥK(J,K)ΥE(J,E˜1)ΥK(J~1,K~1)-KK~1)   +(ΥE(J~1,E˜1)+ΥP(J~1,P~1)+ΥK(J~1,K~1))   (2-ΥE(J~1,E˜1)ΥE(J,E˜1))   -ΥE(J~1,E˜1)Π10f1Π̄1(r)ΥE(Jr,Er)P~1ΥE(J~1,E˜1)Pdr   -ΥP(J~1,P~1)Π10f1Π̄1(r)ΥP(Jr,Pr)P~1ΥP(J~1,P~1)Pdr   -ΥK(J~1,K~1)Π10f1Π̄1(r)ΥK(Jr,Kr)P~1ΥK(J~1,K~1)Pdr   -ΥE(J~1,E˜1)+ΥK(J~1,K~1)Π20f2Π̄2(r)PrK~1P~1Kdr   +2ΥE(J~1,E˜1)   +ΥK(J~1,K~1)-ΥE(J~1,E˜1)Π30f3Π̄3(r)KrE˜1K~1Edr   +ΥE(J~1,E˜1)Π10f1Π̄1(r)[ln(ΥE(Jr,Er)P~1ΥE(J~1,E˜1)P)   +ln(ΥE(J~1,E˜1)ΥE(J,E˜1))   +ln(ΥE(J,E˜1)EΥE(J,E)E˜1)+ln(PK~1P~1K)+ln(KE˜1K~1E)]dr   +ΥP(J~1,P~1)Π10f1Π̄1(r)[ln(ΥP(Jr,Pr)P~1ΥP(J~1,P~1)P)   +ln(ΥE(J~1,E˜1)ΥE(J,E˜1))+ln(ΥE(J,E˜1)ΥP(J~1,P~1)PΥE(J~1,E˜1)ΥP(J,P)P~1)]dr   +ΥK(J~1,K~1)Π10f1Π̄1(r)[ln(ΥK(Jr,Kr)P~1ΥK(J~1,K~1)P)   +ln(ΥE(J~1,E˜1)ΥE(J,E˜1))   +ln(ΥE(J,E˜1)ΥK(J~1,K~1)KΥE(J~1,E˜1)ΥK(J,K)K~1)+ln(PK~1P~1K)]dr   +ΥE(J~1,E˜1)+ΥK(J~1,K~1)Π20f2Π̄2(r)   ×[ln(PrK~1P~1K)+ln(P~1KPK~1)]dr   +ΥE(J~1,E˜1)Π30f3Π̄3(r)[ln(KrE˜1K~1E)+ln(K~1EKE˜1)]dr   -ϕ(ΥE(J~1,E˜1)+ΥK(J~1,K~1))K~1(ηK+ϕL~1)(ξ-βL~1)(ηL+βK)(L-L~1)2.

Then

dΞ4dt=ηJJ~1(1-ΥE(J~1,E˜1)ΥE(J,E˜1))(1-JJ~1)   +ΥE(J~1,E˜1)(ΥE(J,E)ΥE(J,E˜1)-EE˜1)   +ΥP(J~1,P~1)(ΥE(J~1,E˜1)ΥP(J,P)ΥE(J,E˜1)ΥP(J~1,P~1)-PP~1)   +ΥK(J~1,K~1)(ΥE(J~1,E˜1)ΥK(J,K)ΥE(J,E˜1)ΥK(J~1,K~1)-KK~1)   -(ΥE(J~1,E˜1)+ΥP(J~1,P~1)+ΥK(J~1,K~1))   ×(ΥE(J~1,E˜1)ΥE(J,E˜1)-1-ln(ΥE(J~1,E˜1)ΥE(J,E˜1)))   -ΥE(J~1,E˜1)Π10f1Π̄1(r)[ΥE(Jr,Er)P~1ΥE(J~1,E˜1)P-1   -ln(ΥE(Jr,Er)P~1ΥE(J~1,E˜1)P)]dr   -ΥP(J~1,P~1)Π10f1Π̄1(r)[ΥP(Jr,Pr)P~1ΥP(J~1,P~1)P-1   -ln(ΥP(Jr,Pr)P~1ΥP(J~1,P~1)P)]dr   -ΥK(J~1,K~1)Π10f1Π̄1(r)[ΥK(Jr,Kr)P~1ΥK(J~1,K~1)P-1      -ln(ΥK(Jr,Kr)P~1ΥK(J~1,K~1)P)]dr      -ΥE(J~1,E˜1)+ΥK(J~1,K~1)Π20f2Π̄2(r)      ×[PrK~1P~1K-1-ln(PrK~1P~1K)]dr      -ΥE(J~1,E˜1)Π30f3Π̄3(r)[KrE˜1K~1E-1-ln(KrE˜1K~1E)]dr      +ΥE(J~1,E˜1)ln(ΥE(J,E˜1)EΥE(J,E)E˜1)      +ΥP(J~1,P~1)ln(ΥE(J,E˜1)ΥP(J~1,P~1)PΥE(J~1,E˜1)ΥP(J,P)P~1)      +ΥK(J~1,K~1)ln(ΥE(J,E˜1)ΥK(J~1,K~1)KΥE(J~1,E˜1)ΥK(J,K)K~1)      -ϕ(ΥE(J~1,E˜1)+ΥK(J~1,K~1))K~1(ηK+ϕL~1)(ξ-βL~1)(ηL+βK)(L-L~1)2.

By simplifying the previous result, we arrive at

dΞ4dt=ηJJ~1(1-ΥE(J~1,E˜1)ΥE(J,E˜1))(1-JJ~1)   +ΥE(J~1,E˜1)(ΥE(J,E)ΥE(J,E˜1)-EE˜1-1+ΥE(J,E˜1)EΥE(J,E)E˜1)   +ΥP(J~1,P~1)(ΥE(J~1,E˜1)ΥP(J,P)ΥE(J,E˜1)ΥP(J~1,P~1)-PP~1-1   +ΥE(J,E˜1)ΥP(J~1,P~1)PΥE(J~1,E˜1)ΥP(J,P)P~1)      +ΥK(J~1,K~1)(ΥE(J~1,E˜1)ΥK(J,K)ΥE(J,E˜1)ΥK(J~1,K~1)-KK~1-1   +ΥE(J,E˜1)ΥK(J~1,K~1)KΥE(J~1,E˜1)ΥK(J,K)K~1)   -ΥE(J~1,E˜1)Π10f1Π̄1(r)[𝐹(ΥE(Jr,E(t-r))P~1ΥE(J~1,E˜1)P)   +𝐹(ΥE(J~1,E˜1)ΥE(J,E˜1))]dr   -ΥP(J~1,P~1)Π10f1Π̄1(r)[𝐹(ΥP(Jr,Pr)P~1ΥP(J~1,P~1)P)   +𝐹(ΥE(J~1,E˜1)ΥE(J,E˜1))]dr   -ΥK(J~1,K~1)Π10f1Π̄1(r)[𝐹(ΥK(Jr,Kr)P~1ΥK(J~1,K~1)P)   +𝐹(ΥE(J~1,E˜1)ΥE(J,E˜1))]dr   -ΥE(J~1,E˜1)+ΥK(J~1,K~1)Π20f2Π̄2(r)𝐹(PrK~1P~1K)dr   -ΥE(J~1,E˜1)Π30f3Π̄3(r)𝐹(KrE˜1K~1E)dr   -ΥE(J~1,E˜1)[ΥE(J,E˜1)EΥE(J,E)E˜1-1-ln(ΥE(J,E˜1)EΥE(J,E)E˜1)]   -ΥP(J~1,P~1)[ΥE(J,E˜1)ΥP(J~1,P~1)PΥE(J~1,E˜1)ΥP(J,P)P~1-1   -ln(ΥE(J,E˜1)ΥP(J~1,P~1)PΥE(J~1,E˜1)ΥP(J,P)P~1)]   -ΥK(J~1,K~1)[ΥE(J,E˜1)ΥK(J~1,K~1)KΥE(J~1,E˜1)ΥK(J,K)K~1-1   -ln(ΥE(J,E˜1)ΥK(J~1,K~1)KΥE(J~1,E˜1)ΥK(J,K)K~1)]   -ϕ(ΥE(J~1,E˜1)+ΥK(J~1,K~1))K~1(ηK+ϕL~1)(ξ-βL~1)(ηL+βK)(L-L~1)2   =ηJJ~1(1-ΥE(J~1,E˜1)ΥE(J,E˜1))(1-JJ~1)+ΥE(J~1,E˜1)   ×(ΥE(J,E)ΥE(J,E˜1)-EE˜1)(1-ΥE(J,E˜1)ΥE(J,E))   +ΥP(J~1,P~1)(ΥE(J~1,E˜1)ΥP(J,P)ΥE(J,E˜1)ΥP(J~1,P~1)-PP~1)   ×(1-ΥE(J,E˜1)ΥP(J~1,P~1)ΥE(J~1,E˜1)ΥP(J,P))   +ΥK(J~1,K~1)(ΥE(J~1,E˜1)ΥK(J,K)ΥE(J,E˜1)ΥK(J~1,K~1)-KK~1)   ×(1-ΥE(J,E˜1)ΥK(J~1,K~1)ΥE(J~1,E˜1)ΥK(J,K))   -ΥE(J~1,E˜1)Π10f1Π̄1(r)[𝐹(ΥE(Jr,E(t-r))P~1ΥE(J~1,E˜1)P)   +𝐹(ΥE(J~1,E˜1)ΥE(J,E˜1))   +𝐹(ΥE(J,E˜1)EΥE(J,E)E˜1)]dr-ΥP(J~1,P~1)Π10f1Π̄1(r)   ×[𝐹(ΥP(Jr,Pr)P~1ΥP(J~1,P~1)P)+𝐹(ΥE(J~1,E˜1)ΥE(J,E˜1))   +𝐹(ΥE(J,E˜1)ΥP(J~1,P~1)PΥE(J~1,E˜1)ΥP(J,P)P~1)]dr   -ΥK(J~1,K~1)Π10f1Π̄1(r)[𝐹(ΥK(Jr,Kr)P~1ΥK(J~1,K~1)P)   +𝐹(ΥE(J~1,E˜1)ΥE(J,E˜1))   +𝐹(ΥE(J,E˜1)ΥK(J~1,K~1)KΥE(J~1,E˜1)ΥK(J,K)K~1)]dr   -ΥE(J~1,E˜1)+ΥK(J~1,K~1)Π20f2Π̄2(r)   ×𝐹(PrK~1P~1K)dr-ΥE(J~1,E˜1)Π30f3Π̄3(r)𝐹(KrE˜1K~1E)dr   -ϕ(ΥE(J~1,E˜1)+ΥK(J~1,K~1))K~1(ηK+ϕL~1)(ξ-βL~1)   ×(ηL+βK)(L-L~1)2.

From Condition C2, we have

(1-ΥE(J~1,E˜1)ΥE(J,E˜1))(1-JJ~1)0.

In addition, from Equations (1820), we conclude

(ΥE(J,E)ΥE(J,E˜1)-EE˜1)(1-ΥE(J,E˜1)ΥE(J,E))0,(ΥE(J~1,E˜1)ΥP(J,P)ΥE(J,E˜1)ΥP(J~1,P~1)-PP~1)×(1-ΥE(J,E˜1)ΥP(J~1,P~1)ΥE(J~1,E˜1)ΥP(J,P))0,(ΥE(J~1,E˜1)ΥK(J,K)ΥE(J,E˜1)ΥK(J~1,K~1)-KK~1)×(1-ΥE(J,E˜1)ΥK(J~1,K~1)ΥE(J~1,E˜1)ΥK(J,K))0.

Consequently, when ~0>1, it is established that dΞ4dt0 holds true for all positive values of J,P,K,E and L. Additionally, dΞ4dt=0 is realized when J=J~1,P=P~1,K=K~1,E=E˜1, and L=L~1. As a result, with M4 ={Õ1}, we can affirm that the steady state Õ1 is G.A.S under the condition ~0>1 in accordance of L.I.P [44, 45].

4 Illustrative examples and numerical simulations

Dive into the real world! This section brings our theoretical predictions to life through numerical simulations. We use the MATLAB solvers' ode45 and dde23 to explore the dynamics of systems (3) and (23), respectively. We utilize saturated incidence functions for the incidence rates in the system (3), while choosing Holling type-II incidence for system (23). Additionally, within the framework of model (3), we investigate the dynamic impact of CTL immune impairment and saturation effects. Furthermore, our analysis extends to exploring the influences of time delays and the Holling effect on the dynamics of the system (23). We utilize sensitivity analysis to reveal how basic reproduction ratios respond to variations in parameter values.

4.1 Numerical simulation for model (3)

The subsequent saturation incidence rates will be selected for system (3):

ΥE(J,E)=ρEJE1+ϵEE,   ΥP(J,P)=ρPJP1+ϵPP,ΥK(J,K)=ρKJK1+ϵKK,

where ρZ > 0, Z{E,P,K}, account for the infection rate constants. Parameters ϵZ>0,Z{E,P,K}, represent the saturation infection rate constants. Hence, the structure of system (3) will be as follows:

{J˙=θ-ηJJ-ρEJE1+ϵEE-ρPJP1+ϵPP-ρKJK1+ϵKK,P˙=ρEJE1+ϵEE+ρPJP1+ϵPP+ρKJK1+ϵKK-(ψ+ηP)P,K˙=ψP-ηKK-ϕLK,E˙=ϖK-ηEE,L˙=ξK-ηLL-βLK.    (40)

Clearly, ΥZ(J,Z), Z{E,P,K} are continuously differentiable functions. In addition, ΥZ(J,Z), Z{E,P,K} satisfying the following: ΥZ(J,Z)>0 and ΥZ(0,Z)=ΥZ(J,0)=0 for all J>0 and Z > 0. Hence, Condition C1 is fulfilled. Moreover, we have

ΥZ(J,Z)J=ρZZ1+ϵZZ>0,   for all Z>0,ΥZ(J,Z)Z=ρZJ(1+ϵZZ)2>0,   for all J,Z>0,

where Z{E,P,K}. This means that Condition C2 is valid. It is obvious that

Z(J)=ΥZ(J,0)Z=ρZJ>0,   for all J>0, Z{E,P,K}.

Furthermore

Z(J)=ddJ(ΥZ(J,0)Z)=ρZ>0,    Z{E,P,K},

which confirms that Condition C3 is also met. Moreover, we have

Z(ΥZ(J,Z)Z)=-ρZϵZJ(1+ϵZZ)2<0,   for all J,Z>0,Z{E,P,K}.

Then Condition C4 is verified. It is observable that Z(J)E(J)=ρZρE, Z{P,K}. Hence, Condition C5 is satisfied. Furthermore, we have

ΥP(J,P)ΥE(J,E1)=ρPP(1+ϵEE1)ρEE1(1+ϵPP),ΥP(J1,P1)ΥE(J1,E1)=ρPP1(1+ϵEE1)ρEE1(1+ϵPP1),ΥK(J,K)ΥE(J,E1)=ρKK(1+ϵEE1)ρEE1(1+ϵKK),ΥK(J1,K1)ΥE(J1,E1)=ρKK1(1+ϵEE1)ρEE1(1+ϵKK1),

and

(ΥP(J,P)ΥE(J,E1)P-ΥP(J1,P1)ΥE(J1,E1)P1)×(ΥP(J,P)ΥE(J,E1)-ΥP(J1,P1)ΥE(J1,E1))=-ϵPρP2(1+ϵEE1)2(P-P1)2ρE2E12(1+ϵPP)2(1+ϵPP1)20,(ΥK(J,K)ΥE(J,E1)K-ΥK(J1,K1)ΥE(J1,E1)K1)×(ΥK(J,K)ΥE(J,E1)-ΥK(J1,K1)ΥE(J1,E1))=-ϵKρK2(1+ϵEE1)2(K-K1)2ρE2E12(1+ϵKK)2(1+ϵKK1)20,

for all P,K>0,J(0,J0). Consequently, Condition C6 is also satisfied. Given that Conditions C1-C6 are fulfilled, the conclusions regarding global stability, as outlined in Theorems 2.1 and 2.2, persist. System (40) is associated with the basic reproduction ratio, provided as

0(40)=J0(ϖψρE+ηEηKρP+ψηEρK)ηEηK(ψ+ηP)      =0E(40)+0P(40)+0K(40),    (41)

where

0E(40)=ϖψJ0ρEηEηK(ψ+ηP), 0P(40)=J0ρPψ+ηP,0K(40)=ψJ0ρKηK(ψ+ηP).

4.1.1 Stability of steady states

We conduct numerical simulations for the system described in Equation (40), employing the parameter with values outlined in Table 1. Several parameters have been drawn from previous research, complemented by specific choices (ρZ, β, and ϵ = ϵZ for Z{E,P,K}) made for this study. For the assessment of steady state stability in system the (40), we commence simulations using three varied sets of initial states:

Table 1
www.frontiersin.org

Table 1. Summary of model parameters and references.

I.S.1: (J(0),P(0),K(0),E(0),L(0))=(405,4,2,1,0.2),

I.S.2: (J(0),P(0),K(0),E(0),L(0))=(350,10,2.85,3.5,0.5),

I.S.3: (J(0),P(0),K(0),E(0),L(0))=(500,20,4,4,0.6).

Utilizing the infection rates ρZ for Z{E,P,K}, the stability of steady states is governed by ℜ0(40), as defined in Equation (41). To scrutinize the impact of variations in these parameters, we explore two distinct scenarios:

Stability of O0(40). For the parameters ρE=0.0003, ρP=0.0001, ρK=0.00002, ϵ = 0.02, and β = 0.001, the computed value of ℜ0(40) = 0.503 is below 1. As illustrated in Figure 1, the trajectories originating I.S.1-I.S.3 approach the steady state O0(40) = (1000, 0, 0, 0, 0). This observation reinforces the assertion that O0(40) is being G.A.S, aligning with the outcome established in Theorem 2.1. Biologically speaking, this situation implies the complete eradication of the infection, signifying the successful elimination of the pathogen by the human body.

Figure 1
www.frontiersin.org

Figure 1. Solution patterns of the dynamical system (40) in the state of ℜ0(40) ≤ 1. (A) Healthy CD4+T cells. (B) Latently infected cells. (C) Actively infected cells. (D) HIV-1 particles. (E) CTLs.

Stability of O1(40). Choosing ρE=0.003, ρP=0.0002, ρK=0.0001, ϵ = 0.02, and β = 0.001, yields ℜ0(40) = 2.804>1. Consequently, the steady state O1(40) exists when ℜ0(40) > 1, with specific values given by O1(40) = (408.467, 15.987, 3.903, 4.229, 0.479). In Figure 2, the numerical results align with the findings of Theorem 2.2, indicating the convergence of solutions for system (40) to O1(40) when ℜ0(40) > 1 across all I.S.1-I.S.3. From a biological standpoint, this situation implies the coexistence of both HIV-1 particles and CTLs within the host organism.

Figure 2
www.frontiersin.org

Figure 2. Solution patterns of the dynamical system (40) in the state of ℜ0(40) > 1. (A) Healthy CD4+T cells. (B) Latently infected cells. (C) Actively infected cells. (D) HIV-1 particles. (E) CTLs.

4.1.2 Impact of impaired CTL immunity

In this given context, we introduce fluctuations in the β value while maintaining specific values for ρE=0.003, ρP=0.0002, ρK=0.0001, and ϵ = ϵZ = 0.02, where Z{E,P,K}. We aim to explore how CTL immune impairment impacts the dynamics of the system (40), by obtaining numerical solutions with different β values as provided in Table 2. Under these circumstances, we employ the subsequent initial state:

Table 2
www.frontiersin.org

Table 2. Effect of the immune impairment parameter of CTLs on system (40).

I.S.4: (J(0),P(0),K(0),E(0),L(0))=(405,16.06,4,4.3,0.3).

As β increases, we notice a decline in the quantity of CTLs, as shown in Table 2. This decrease is accompanied by a greater quantity of infected cells, whether latent or active, along with a higher quantity of HIV-1 particles. Due to this, there is a reduction in the count of healthy cells. Notably, Figure 3 reveals that the stability criteria of the steady states remain unaffected by CTL immune impairment. This is a consequence of the fact that the parameter ℜ0(40) stays constant regardless of varying β values.

Figure 3
www.frontiersin.org

Figure 3. Solution patterns of the dynamical system (40) along different β values. (A) Healthy CD4+T cells. (B) Latently infected cells. (C) Actively infected cells. (D) HIV-1 particles. (E) CTLs.

4.1.3 The saturation infection rates effect on the system dynamics

In this instance, we fix specific values for ρE=0.003, ρP=0.0002, ρK=0.0001, and β = 0.001. We numerically solve system (40) with various values of ϵ = ϵZ for Z{E,P,K}, accompanied by the following initial state:

I.S.5: (J(0),P(0),K(0),E(0),L(0))=(500,12,4,4.3,0.3).

Upon examining the results in Table 3, a positive correlation is apparent between ϵ and healthy cells, meaning that when ϵ becomes higher, the concentration of healthy cells becomes higher too. In contrast, the counts of other compartments become lower. Notably, the parameter ℜ0(40) remains independent of ϵ. Consequently, Figure 4 illustrates that altering ϵ does not have an impact on the steady states' stability.

Table 3
www.frontiersin.org

Table 3. Exploring system (40) dynamics under the influence of saturation infection rates.

Figure 4
www.frontiersin.org

Figure 4. The influence of saturation infection rates on the solution patterns in system (40). (A) Healthy CD4+T cells. (B) Latently infected cells. (C) Actively infected cells. (D) HIV-1 particles. (E) CTLs.

4.2 Numerical simulation for model (23)

In the context of numerical analysis, we adopt the following specific form regarding the probability distribution functions πi(r), for i = 1, 2, 3:

πi(r)=β*(r-ri), ri[0,fi], i=1,2,3,

where β*(.) is the Dirac delta function. When fi → ∞, we observe

0πi(r)dr=1, i=1,2,3.

Furthermore, we obtain

Πi=0β*(r-ri)e-birdr=e-biri,      i=1,2,3.

Additionally, we will select the Holling type-II incidence rates of infection as follows:

ΥE(J,E)=ρEJE1+ϑJ,   ΥP(J,P)=ρPJP1+ϑJ,ΥK(J,K)=ρKJK1+ϑJ,

Here, the positive constants ρZ for Z{E,P,K} denote the infection rate constants, while ϑ > 0 represents the Holling infection rate constant. Consequently, Based on the aforementioned reasoning, system (23) was modified into a system with a discrete time delay, as illustrated below:

{J˙=θ-ηJJ-ρEJE1+ϑJ-ρPJP1+ϑJ-ρKJK1+ϑJ,P˙=e-b1r1(ρEJ(t-r1)E(t-r1)1+ϑJ(t-r1)   +ρPJ(t-r1)P(t-r1)1+ϑJ(t-r1)    +ρKJ(t-r1)K(t-r1)1+ϑJ(t-r1))-(ψ+ηP)P,K˙=ψe-b2r2P(t-r2)-ηKK-ϕLK,E˙=ϖe-b3r3K(t-r3)-ηEE,L˙=ξK-ηLL-βLK.    (42)

The functions ΥZ(J,Z) for Z{E,P,K} are evidently continuously differentiable. Additionally, for each Z{E,P,K}, these functions satisfy the conditions: ΥZ(J,Z)>0 and ΥZ(0,Z)=ΥZ(J,0)=0 for all J>0 and Z > 0. Thus, the fulfillment of Condition C1 is evident. Furthermore, we observe

ΥZ(J,Z)J=ρZZ(1+ϑJ)2>0,   for all J,Z>0,ΥZ(J,Z)Z=ρZJ1+ϑJ>0,   for all J>0,

where Z{E,P,K}. This can be expressed by stating that Condition C2 is valid. It is obvious that

Z(J)=ΥZ(J,0)Z=ρZJ1+ϑJ>0,for all J>0, Z{E,P,K}.

Furthermore

Z(J)=ddJ(ΥZ(J,0)Z)=ρZ(1+ϑJ)2>0,   Z{E,P,K},

which confirms that Condition C3 is also met. Moreover, we have

Z(ΥZ(J,Z)Z)=0,   Z{E,P,K}.

Hence, the verification of Condition C4 is straightforward. Additionally, it is evident that Z(J)E(J)=ρZρE for Z{P,K}. Consequently, Condition C5 is satisfied. Moreover, we have:

ΥP(J,P)ΥE(J,E˜1)=ρPPρEE˜1,    ΥP(J~1,P~1)ΥE(J~1,E˜1)=ρPP~1ρEE˜1,ΥK(J,K)ΥE(J,E˜1)=ρKKρEE˜1,    ΥK(J~1,K~1)ΥE(J~1,E˜1)=ρKK~1ρEE˜1,

and

(ΥP(J,P)ΥE(J,E˜1)P-ΥP(J~1,P~1)ΥE(J~1,E˜1)P~1)×(ΥP(J,P)ΥE(J,E˜1)-ΥP(J~1,P~1)ΥE(J~1,E˜1))=0,(ΥK(J,K)ΥE(J,E˜1)K-ΥK(J~1,K~1)ΥE(J~1,E˜1)K~1)×(ΥK(J,K)ΥE(J,E˜1)-ΥK(J~1,K~1)ΥE(J~1,E˜1))=0,

for all P,K>0,J(0,J~0). Hence, Condition C6 is fulfilled as well. With the fulfillment of Conditions C1C6, the global stability results stated in Theorems 3.1 and 3.2 persist. Accordingly, the determination of the basic reproduction ratio for system (42) is outlined below:

~0(42)=J~0e-b1r1(ψe-b2r2(ρEϖe-b3r3+ρKηE)+ρPηEηK)ηEηK(ψ+ηP)(1+ϑJ~0)      =~0E(42)+~0P(42)+~0K(42),    (43)

where

~0E(42)=ϖψJ~0ρEe-b1r1-b2r2-b3r3ηEηK(ψ+ηP)(1+ϑJ~0),~0P(42)=J~0ρPe-b1r1(ψ+ηP)(1+ϑJ~0),~0K(42)=ψJ~0ρKe-b1r1-b2r2ηK(ψ+ηP)(1+ϑJ~0).

4.2.1 The stability implications of time delays on steady states

To assess how time delay parameters affect the solutions of system (42), we maintain constant values for parameters ρE=0.003, ρP=0.0001, ρK=0.0004, β = 0.001, b1 = 0.1, b2 = 0.2, b3 = 0.3, and set ϑ = 0.0001. Furthermore, you can refer to Table 1 for the values of the other parameters. Next, we examine how the dynamics are affected by experimenting with different variations of the delay parameters ri, where i takes values from 1 to 3. The stability of steady states is highly sensitive to changes in the parameters ri due to the dependence of ~0(42) on them [as indicated in Equation (43)]. Let's delve into the following choices for the delay parameters:

D.P.1: r1 = 0.003, r2 = 0.001, r3 = 0.002.

D.P.2: r1 = 0.01, r2 = 0.02, r3 = 0.03.

D.P.3: r1 = 0.8, r2 = 0.6, r3 = 0.7.

D.P.4: r1 = 1.189, r2 = 1.5, r3 = 2.5.

D.P.5: r1 = 3, r2 = 4, r3 = 5.

D.P.6: r1 = 5, r2 = 6, r3 = 7. We address the initial state for solving system (42) as follows:

I.S.6: (J(v),P(v),K(v),E(v),L(v))=(600,5,2,1,0.2), v∈[−r, 0], r = max{r1, r2, r3}. Table 4 displays the values of ~0(42) corresponding to different ri values, where i = 1, 2, 3. Significantly reducing ~0(42) is noteworthy when the ri parameters increase. Figure 5 visually depicts the numerical solutions derived from the system, emphasizing the significant impact of the included time-delays. Specifically, a positive correlation exists between ri and the count of healthy cells, indicating that their values increase simultaneously, while a decrease is observed in the counts of other compartments. This suggests the potential creation of a new category of medication aimed at extending delay times and ultimately mitigating HIV-1 multiplication.

Table 4
www.frontiersin.org

Table 4. Analyzing the influence of delay parameters, ri, on System (42).

Figure 5
www.frontiersin.org

Figure 5. The influence of ri on the solution patterns in system (42). (A) Healthy CD4+T cells. (B) Latently infected cells. (C) Actively infected cells. (D) HIV-1 particles. (E) CTLs.

4.2.2 The Holling type-II infection rates effect on the system dynamics

In this scenario, we set ρE=0.003, ρP=0.0001, ρK=0.0004, b1 = 0.1, b2 = 0.2, b3 = 0.3, β = 0.001, and ri = 0.002 for i = 1, 2, 3. We consider different values of ϑ and numerically solve system (42) with the initial state I.S.6 to analyze the impact of Holling infection rates on the dynamics of the system (42). Analyzing Table 5 reveals a correlation wherein an augmentation in ϑ leads to an increase in the count of healthy cells, while there is a reduction in the counts of other compartments. Clearly, the parameter ~0(42) depends on ϑ. Therefore, Figure 6 confirms that the Holling infection alters the stability properties of steady states. Furthermore, we derive the subsequent observations:

1. Õ1(42) is being G.A.S under the condition 0 < ϑ < 0.0017336.

2. Õ0(42) is being G.A.S under the condition ϑ≥0.0017336.

Table 5
www.frontiersin.org

Table 5. The dynamics of system (42) affected by Holling infection rates.

Figure 6
www.frontiersin.org

Figure 6. The influence of Holling infection rates on the solution patterns in system (42). (A) Healthy CD4+T cells. (B) Latently infected cells. (C) Actively infected cells. (D) HIV-1 particles. (E) CTLs.

4.3 Analysis of parameter sensitivity

4.3.1 Investigating the sensitivity of parameters for model (40)

The principal goal of analyzing parameter sensitivity is to determine the variable that has the most substantial impact on ℜ0, which should be the focus of control strategies. To achieve this, we will employ direct differentiation, a method that calculates sensitivity indices by considering changes in variables based on parameters. We can define the normalized forward sensitivity index of ℜ0, which depends on the differentiability with respect to a parameter ϱ, as follows:

Sϱ=ϱ00ϱ.    (44)

For each parameter belonging to ℜ0(40), we utilized Equation (44) to compute the sensitivity index values. For this purpose, we considered the parameters specified in Table 1, while also incorporating supplementary values for ρE=0.003, ρP=0.0002, and ρK=0.0001. Sensitivity index values for ℜ0(40) are showcased in Table 6 and Figure 7, revealing positive index values for θ, ρE, ρP, ρK, ϖ, and ψ. In this context, higher values of these parameters correspond to a higher prevalence of HIV-1 infection within this population. In contrast, the remaining indices show negativity, implying that an increase in ηJ, ηP, ηK, and ηE corresponds to a reduction in the ℜ0(40) value. Notably, the parameters θ, ρE, and ϖ emerge as the pivotal factors in terms of sensitivity, whereas ρP, ρK, and ψ are considered less critical. Interestingly, the parameter ξ, representing CTLs responsiveness, demonstrates no impact on ℜ0(40).

Table 6
www.frontiersin.org

Table 6. Examining sensitivity in parameters regarding ℜ0(40).

Figure 7
www.frontiersin.org

Figure 7. Analyzing sensitivity of parameters for ℜ0(40) in model (40).

4.3.2 Investigating the sensitivity of parameters for model (42)

Here, Equation (44) will be applied to calculate sensitivity indices concerning ~0(42), considering each parameter contributing to the basic reproduction ratio. The calculations were performed with the parameters specified in Table 1, supplemented by the inclusion of the additional parameters: ρE=0.003, ρP=0.0001, ρK=0.0004, b1 = 0.1, b2 = 0.2, b3 = 0.3, r1 = 0.8, r2 = 0.6, r3 = 0.7, and ϑ = 0.0001. The sensitivity index values for ~0(42) are presented in Table 7, and Figure 8 provides a graphical representation. Positive indices for θ, ρE, ρP, ρK, ϖ, and ψ suggest that an increase in these values will elevate the prevalence of HIV-1 disease. Conversely, negative indices, including ηJ, ηP, ηK, ηE, b1, b2, b3, r1, r2, r3, and ϑ, indicate that an increase in these values will lower the parameter ~0(42). Notably, θ, ρE, and ϖ emerge as the most influential parameters, while ρP, ρK, and ψ are considered less critical. Furthermore, the parameter representing CTLs responsiveness, ξ, demonstrates no discernible effect on ~0(42).

Table 7
www.frontiersin.org

Table 7. Examining sensitivity in parameters regarding ~0(42).

Figure 8
www.frontiersin.org

Figure 8. Analyzing sensitivity of parameters for ~0(42) in model (42).

5 Discussion

To emphasize the importance of integrating L-CTC spread into our investigated models, we analyze model (3) in the context of various antiviral interventions, outlined as follows:

{J˙=θ-ηJJ-(1-ς1)ΥE(J,E)-(1-ς2)ΥP(J,P)   -(1-ς3)ΥK(J,K),P˙=(1-ς1)ΥE(J,E)+(1-ς2)ΥP(J,P)   +(1-ς3)ΥK(J,K)-(ψ+ηP)P,K˙=ψP-ηKK-ϕLK,E˙=ϖK-ηEE,L˙=ξK-ηLL-βLK.    (45)

The parameter 0 ≤ ς1 ≤ 1 represents the effectiveness of antiviral therapy in inhibiting VTC transmission, while 0 ≤ ς2, ς3 ≤ 1 denote the therapeutic efficacies in blocking L-CTC and A-CTC transmissions, respectively [62].

For system (45), the following establishes the basic reproduction ratio:

0=(1-ς1)ϖψηEηK(ψ+ηP)ΥE(J0,0)E+(1-ς2)ψ+ηPΥP(J0,0)P+(1-ς3)ψηK(ψ+ηP)ΥK(J0,0)K.

Assuming ς to be equal to ς1, ς2, and ς3, we obtain:

0ς=(1-ς)[ϖψηEηK(ψ+ηP)ΥE(J0,0)E   +1ψ+ηPΥP(J0,0)P+ψηK(ψ+ηP)ΥK(J0,0)K]   =(1-ς)0.

Currently, we compute the drug efficacy, ς, which results in 0ς<1 and stabilizes O0 in system (45) as follows:

1ς>ς~min=max{0,1-10}.    (46)

Neglecting the spread of L-CTC in model (45) yields

{J˙=θ-ηJJ-(1-ς)ΥE(J,E)-(1-ς)ΥK(J,K),P˙=(1-ς)ΥE(J,E)+(1-ς)ΥK(J,K)-(ψ+ηP)P,K˙=ψP-ηKK-ϕLK,E˙=ϖK-ηEE,L˙=ξK-ηLL-βLK,    (47)

and the definition of the basic reproduction ratio of model (47) will be

^0ς=(1-ς)[ϖψηEηK(ψ+ηP)ΥE(J0,0)E   +ψηK(ψ+ηP)ΥK(J0,0)K]=(1-ς)^0.

We identify the drug efficacy, denoted as ς, that ensures ^0ς<1, thereby stabilizing the steady state O0 of system (47). The condition is expressed as:

1ς>ς^min=max{0,1-1^0}.    (48)

It is evident that ^0<0. Interestingly, comparing Equations (46, 48) reveal that ς^min is always smaller than ς~min. Therefore, applying drugs with efficacy ς in case of ς^minς<ς~min will guarantee that ^0ς<1, ensuring the global asymptotic stability of O0 in the system (47). However, 0ς>1, indicating the instability of O0 in system (45). As a result, drug therapies are designed without considering the spread of L-CTC may be inaccurate or insufficient for achieving viral eradication from the body. Therefore, compared to the model given in Alofi and Azoz [24] and Zhang and Xu [25], our suggested models are more pertinent in explaining HIV-1 infections.

6 Conclusion

This work formulates two HIV-1 dynamics models with CTL immune impairment. These models comprise five compartments: healthy CD4+T cells, latently infected cells, actively infected cells, free HIV-1 and CTLs. To enhance realism, we introduced a scenario in which healthy cells acquire susceptibility to infection upon exposure to infectious compartments (HIV-1 particles, latently infected cells and actively infected cells). The second model takes the first model a step further by incorporating multiple distributed-time delays. Both models incorporate general functions to represent the incidence rates. Notably, the solutions produced by these models display properties of nonnegativity and boundedness. Within this framework, we identified two crucial steady states: the virus-free steady state, O00), and the virus-persistence steady state, O11). We calculated the basic reproduction ratios, referred to as ℜ0 (~0), which are essential for determining whether the steady states mentioned earlier exist and remain stable over time. It is noteworthy that ℜ0 (~0) consists of three separate components, representing the contributions from VTC, L-CTC, and A-CTC. The Lyapunov function technique and L.I.P. are employed to comprehend the system's behavior over time by studying the global asymptotic stability of the steady states. Two key findings emerged from our investigation: firstly, the virus-free steady state O00) is being G.A.S whenever ℜ0 < 1 (~0<1), ultimately leading to the disappearance of the infection. In contrast, the steady state O00) loses its stability and the virus-persistence steady state O11) is being G.A.S whenever ℜ0 > 1 (~0>1), indicating a long-term infection. For experimental verification of our theoretical derivations, we employed numerical simulations, which yielded results consistent with our analytical solutions. Our exploration encompassed the influence of CTL immune impairment, time delays, and L-CTC on HIV-1 dynamics, revealing reduced immune function as a key factor contributing significantly to disease progression. Moreover, time delays significantly reduced the basic reproduction ratio (~0), leading to suppressed HIV-1 reproduction. Therefore, eliminating HIV-1 requires prioritizing strategies that keep ~0 below 1. Our study also highlights the critical role of L-CTC spread in HIV-1 dynamics. Additionally, a sensitivity analysis revealed key factors influencing the system's behavior, deepening our understanding of HIV-1 dynamics.

We are optimistic that the strategy outlined in our article may be improved upon. The modeling of the in-host dynamics with an emphasis on the immune system vs invading particles would be an intriguing viewpoint. Recently, active-particle approaches have been utilized to describe the immunological rivalry in detail at the cellular level, allowing for the modeling of epidemics (see [63]).

It seems like an interesting avenue to investigate the memory impact on the dynamics of our model with fractional differential equations (FDEs) [6466]. Non-local and memory-dependent effects, which are frequently essential in virological systems, can be captured by FDEs. In our proposed modeled we assumed that the production rate of healthy CD4+T cells is constant. Nonetheless, HIV-1 infection possesses the capability to infect and disrupt the thymus's supply of these cells. As a result, rather than being constant, the thymus's supply of new cells has been seen to be changeable [55, 67]. These topics are left for later research.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

AE: Conceptualization, Formal analysis, Methodology, Writing – original draft, Writing – review & editing. NA: Conceptualization, Formal analysis, Investigation, Methodology, Writing – original draft, Writing – review & editing. RH: Formal analysis, Investigation, Methodology, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare that no financial support was received for the research, authorship, and/or publication of this article.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Nowak MA, May RM. Virus Dynamics. New York, NY: Oxford University Press (2000). doi: 10.1093/oso/9780198504184.001.0001

Crossref Full Text | Google Scholar

2. Nowak MA, Bangham CRM. Population dynamics of immune responses to persistent viruses. Science. (1996) 272:74–9. doi: 10.1126/science.272.5258.74

PubMed Abstract | Crossref Full Text | Google Scholar

3. Jan R, Khan A, Boulaaras S, Zubair SA. Dynamical behaviour and chaotic phenomena of HIV infection through fractional calculus. Discrete Dyn Nat Soc. (2022) 2022:5937420. doi: 10.1155/2022/5937420

Crossref Full Text | Google Scholar

4. Jan A, Srivastava HM, Khan A, Mohammed PO, Jan R, Hamed YS, et al. In vivo HIV dynamics, modeling the interaction of HIV and immune system via non-integer derivatives. Fractal Fract. (2023) 7:361. doi: 10.3390/fractalfract7050361

Crossref Full Text | Google Scholar

5. Jan R, Boulaaras S, Shah SAA. Fractional-calculus analysis of human immunodeficiency virus and CD4+ T-cells with control interventions. Commun Theor Phys. (2022) 74:105001. doi: 10.1088/1572-9494/ac7e2b

Crossref Full Text | Google Scholar

6. Jolly C, Sattentau Q. Retroviral spread by induction of virological synapses. Traffic. (2004) 5:643–50. doi: 10.1111/j.1600-0854.2004.00209.x

PubMed Abstract | Crossref Full Text | Google Scholar

7. Komarova NL, Wodarz D. Virus dynamics in the presence of synaptic transmission. Math Biosci. (2013) 242:161–71. doi: 10.1016/j.mbs.2013.01.003

PubMed Abstract | Crossref Full Text | Google Scholar

8. Sigal A, Kim JT, Balazs AB, Dekel E, Mayo A, Milo R, et al. Cell-to-cell spread of HIV permits ongoing replication despite antiretroviral therapy. Nature. (2011) 477:95–8. doi: 10.1038/nature10347

PubMed Abstract | Crossref Full Text | Google Scholar

9. Gao Y, Wang J. Threshold dynamics of a delayed nonlocal reaction-diffusion HIV infection model with both cell-free and cell-to-cell transmissions. J Math Anal Appl. (2020) 488:124047. doi: 10.1016/j.jmaa.2020.124047

Crossref Full Text | Google Scholar

10. Graw F, Perelson AS. Modeling viral spread. Ann Rev Viro. (2016) 3:555–72. doi: 10.1146/annurev-virology-110615-042249

PubMed Abstract | Crossref Full Text | Google Scholar

11. Guo T, Qiu Z. The effects of CTL immune response on HIV infection model with potent therapy, latently infected cells and cell-to-cell viral transmission. Math Biosci Eng. (2019) 16:6822–41. doi: 10.3934/mbe.2019341

PubMed Abstract | Crossref Full Text | Google Scholar

12. Wang J, Guo M, Liu X, Zhao Z. Threshold dynamics of HIV-1 virus model with cell-to-cell transmission, cell-mediated immune responses and distributed delay. Appl Math Comput. (2016) 291:149–61. doi: 10.1016/j.amc.2016.06.032

Crossref Full Text | Google Scholar

13. Yang Y, Zhang T, Xu Y, Zhou J. A delayed virus infection model with cell-to-cell transmission and CTL immune response. Int J Bifurcat Chaos. (2017) 27:1750150. doi: 10.1142/S0218127417501504

Crossref Full Text | Google Scholar

14. Cervantes-Perez AG, Avila-Vales E. Dynamical analysis of multipathways and multidelays of general virus dynamics model. Int J Bifurcat Chaos. (2019) 29:1950031–30. doi: 10.1142/S0218127419500317

Crossref Full Text | Google Scholar

15. Elaiw AM, AlShamrani NH. Global stability of a delayed adaptive immunity viral infection with two routes of infection and multi-stages of infected cells. Commun Nonlinear Sci Numer Simul. (2020) 86:105259. doi: 10.1016/j.cnsns.2020.105259

PubMed Abstract | Crossref Full Text | Google Scholar

16. Regoes R, Wodarz D, Nowak MA. Virus dynamics: the effect to target cell limitation and immune responses on virus evolution. J Theor Biol. (1998) 191:451–62. doi: 10.1006/jtbi.1997.0617

PubMed Abstract | Crossref Full Text | Google Scholar

17. Hu Z, Zhang J, Wang H, Ma W, Liao F. Dynamics analysis of a delayed viral infection model with logistic growth and immune impairment. Appl Math Model. (2014) 38:524–34. doi: 10.1016/j.apm.2013.06.041

PubMed Abstract | Crossref Full Text | Google Scholar

18. Krishnapriya P, Pitchaimani M. Modeling and bifurcation analysis of a viral infection with time delay and immune impairment. Jpn J Ind Appl Math. (2017) 34:99–139. doi: 10.1007/s13160-017-0240-5

PubMed Abstract | Crossref Full Text | Google Scholar

19. Eric AV, Noe CC, Gerardo GA. Analysis of a viral infection model with immune impairment, intracellular delay and general non-linear incidence rate. Chaos Solitons Fractals. (2014) 69:1–9. doi: 10.1016/j.chaos.2014.08.009

Crossref Full Text | Google Scholar

20. Wang S, Song X, Ge Z. Dynamics analysis of a delayed viral infection model with immune impairment. Appl Math Model. (2011) 35:4877–85. doi: 10.1016/j.apm.2011.03.043

Crossref Full Text | Google Scholar

21. Krishnapriya P, Pitchaimani M. Analysis of time delay in viral infection model with immune impairment. J Appl Math Comput. (2017) 55:421–53. doi: 10.1007/s12190-016-1044-5

Crossref Full Text | Google Scholar

22. Wang ZP, Liu XN. A chronic viral infection model with immune impairment. J Theor Biol. (2007) 249:532–42. doi: 10.1016/j.jtbi.2007.08.017

PubMed Abstract | Crossref Full Text | Google Scholar

23. Elaiw AM, Raezah A, Alofi BS. Dynamics of delayed pathogen infection models with pathogenic and cellular infections and immune impairment. AIP Adv. (2018) 8:025323. doi: 10.1063/1.5023752

Crossref Full Text | Google Scholar

24. Alofi BS, Azoz SA. Stability of general pathogen dynamic models with two types of infectious transmission with immune impairment. AIMS Mathem. (2020) 6:114–40. doi: 10.3934/math.2021009

Crossref Full Text | Google Scholar

25. Zhang L, Xu R. Dynamics analysis of an HIV infection modelwith latent reservoir, delayed CTL immune response and immune impairment. Nonlinear Anal Model Control. (2023) 28:1–19. doi: 10.15388/namc.2023.28.29615

PubMed Abstract | Crossref Full Text | Google Scholar

26. Xu R, Song C. Dynamics of an HIV infection model with virus diffusion and latently infected cell activation. Nonlinear Anal Real World Appl. (2022) 67:103618. doi: 10.1016/j.nonrwa.2022.103618

Crossref Full Text | Google Scholar

27. Nath BJ, Sadri K, Sarmah HK, Hosseini K. An optimal combination of antiretroviral treatment and immunotherapy for controlling HIV infection. Math Comput Simul. (2024) 217:226–43. doi: 10.1016/j.matcom.2023.10.012

Crossref Full Text | Google Scholar

28. Agosto L, Herring M, Mothes W, Henderson A. HIV-1-infected CD4+ T cells facilitate latent infection of resting CD4+ T cells through cell-cell contact. Cell. (2018) 24:2088–100. doi: 10.1016/j.celrep.2018.07.079

PubMed Abstract | Crossref Full Text | Google Scholar

29. Wang W, Wang X, Guo K, Ma W. Global analysis of a diffusive viral model with cell-to-cell infection and incubation period. Math Methods Appl Sci. (2020) 43:5963–78. doi: 10.1002/mma.6339

Crossref Full Text | Google Scholar

30. Hattaf K, Dutta H. Modeling the dynamics of viral infections in presence of latently infected cells. Chaos Solitons Fractals. (2020) 136:109916. doi: 10.1016/j.chaos.2020.109916

PubMed Abstract | Crossref Full Text | Google Scholar

31. Elaiw AM, AlShamrani NH. Stability of a general CTL-mediated immunity HIV infection model with silent infected cell-to-cell spread. Adv Diff Equ. (2020) 2020:355. doi: 10.1186/s13662-020-02818-3

Crossref Full Text | Google Scholar

32. Li J, Wang X, Chen Y. Analysis of an age-structured HIV infection model with cell-to-cell transmission. Eur Phys J Plus. (2024) 139:78. doi: 10.1140/epjp/s13360-024-04873-1

PubMed Abstract | Crossref Full Text | Google Scholar

33. Bai N, Xu R. Mathematical analysis of an HIV model with latent reservoir, delayed CTL immune response and immune impairment. Math Biosc Eng. (2021) 18:1689–707. doi: 10.3934/mbe.2021087

PubMed Abstract | Crossref Full Text | Google Scholar

34. AlShamrani NH, Halawani RH, Shammakh W, Elaiw AM. Global properties of HIV-1 dynamics models with CTL immune impairment and latent cell-to-cell spread. Mathematics. (2023) 11:3743. doi: 10.3390/math11173743

Crossref Full Text | Google Scholar

35. AlShamrani NH, Halawani RH, Elaiw AM. Effect of impaired B-cell and CTL functions on HIV-1 dynamics. Mathematics. (2023) 11:4385. doi: 10.3390/math11204385

Crossref Full Text | Google Scholar

36. Ma Y, Yu X. The effect of environmental noise on threshold dynamics for a stochastic viral infection model with two modes of transmission and immune impairment. Chaos Solitons Fractals. (2020) 134:109699. doi: 10.1016/j.chaos.2020.109699

Crossref Full Text | Google Scholar

37. Yang Y, Xu R. Mathematical analysis of a delayed HIV infection model with saturated CTL immune response and immune impairment. J Appl Math Comput. (2022) 68:2365–80. doi: 10.1007/s12190-021-01621-x

Crossref Full Text | Google Scholar

38. Yaagoub Z, Sadki M, Allali K. A generalized fractional hepatitis B virus infection modelwith both cell-to-cell and virus-to-cell transmissions. Nonlinear Dyn. (2024). doi: 10.21203/rs.3.rs-3958680/v1

PubMed Abstract | Crossref Full Text | Google Scholar

39. Elaiw AM, AlShamrani NH. Global stability of humoral immunity virus dynamics models with nonlinear infection rate and removal. Nonlinear Anal: Real World Appl. (2015) 26:161–90. doi: 10.1016/j.nonrwa.2015.05.007

Crossref Full Text | Google Scholar

40. Shu H, Chen Y, Wang L. Impacts of the cell-free and cell-to-cell infection modes on viral dynamics. J Dynam Differential Equations. (2018) 30:1817–36. doi: 10.1007/s10884-017-9622-2

Crossref Full Text | Google Scholar

41. Smith HL, Waltman P. The Theory of the Chemostat: Dynamics of Microbial Competition. Cambridge: Cambridge University Press (1995). doi: 10.1017/CBO9780511530043

Crossref Full Text | Google Scholar

42. van den Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math Biosci. (2002) 180:29–48. doi: 10.1016/S0025-5564(02)00108-6

PubMed Abstract | Crossref Full Text | Google Scholar

43. Korobeinikov A. Global properties of basic virus dynamics models. Bull Math Biol. (2004) 66:879–83. doi: 10.1016/j.bulm.2004.02.001

PubMed Abstract | Crossref Full Text | Google Scholar

44. Hale JK, Lunel SMV. Introduction to Functional Differential Equations. New York, NY: Springer-Verlag (1993). doi: 10.1007/978-1-4612-4342-7

Crossref Full Text | Google Scholar

45. Khalil HK. Nonlinear Systems, 3rd Edn. Upper Saddle River, NJ: Prentice Hall. (2002).

Google Scholar

46. Perelson A, Neumann A, Markowitz M, Leonard J, Ho D. HIV-1 dynamics in vivo: virion clearance rate, infected cell life-span, and viral generation time. Science. (1996) 271:1582–6. doi: 10.1126/science.271.5255.1582

PubMed Abstract | Crossref Full Text | Google Scholar

47. Culshaw RV, Ruan S, Webb G. A mathematical model of cell-to-cell spread of HIV-1 that includes a time delay. J Math Biol. (2003) 46:425–44. doi: 10.1007/s00285-002-0191-5

PubMed Abstract | Crossref Full Text | Google Scholar

48. Lai X, Zou X. Modelling HIV-1 virus dynamics with both virus-to-cell infection and cell-to-cell transmission. SIAM J Appl Math. (2014) 74:898–917. doi: 10.1137/130930145

Crossref Full Text | Google Scholar

49. Elaiw AM, Raezah AA. Stability of general virus dynamics models with both cellular and viral infections and delays. Math Methods Appl Sci. (2017) 40:5863–80. doi: 10.1002/mma.4436

Crossref Full Text | Google Scholar

50. Adak D, Bairagi N. Analysis and computation of multi-pathways and multi-delays HIV-1 infection model. Appl Math Model. (2018) 54:517–36. doi: 10.1016/j.apm.2017.09.051

Crossref Full Text | Google Scholar

51. Yang Y, Zou L, Ruan S. Global dynamics of a delayed within-host viral infection model with both virus-to-cell and cell-to-cell transmissions. Math Biosci. (2015) 270:183–91. doi: 10.1016/j.mbs.2015.05.001

PubMed Abstract | Crossref Full Text | Google Scholar

52. Kuang Y. Delay Differential Equations with Applications in Population Dynamics. San Diego: Academic Press (1993).

Google Scholar

53. Sahani SK, Yashi. Effects of eclipse phase and delay on the dynamics of HIV infection. J Biol Syst. (2018) 26:421–54. doi: 10.1142/S0218339018500195

Crossref Full Text | Google Scholar

54. Li F, Ma W. Dynamics analysis of an HTLV-1 infection model with mitotic division of actively infected cells and delayed CTL immune response. Math Methods Appl Sci. (2018) 41:3000–17. doi: 10.1002/mma.4797

Crossref Full Text | Google Scholar

55. Perelson AS, Kirschner DE, de Boer R. Dynamics of HIV infection of CD4+ T cells. Math Biosci. (1993) 114:81–125. doi: 10.1016/0025-5564(93)90043-A

PubMed Abstract | Crossref Full Text | Google Scholar

56. Callaway DS, Perelson AS. HIV-1 infection and low steady state viral loads. Bull Math Biol. (2002) 64:29–64. doi: 10.1006/bulm.2001.0266

PubMed Abstract | Crossref Full Text | Google Scholar

57. Mohri H, Bonhoeffer S, Monard S, Perelson AS, Ho DD. Rapid turnover of T lymphocytes in SIV-infected rhesus macaques. Science. (1998) 279:1223–7. doi: 10.1126/science.279.5354.1223

PubMed Abstract | Crossref Full Text | Google Scholar

58. Wang Y, Liu J, Liu L. Viral dynamics of an HIV model with latent infection incorporating antiretroviral therapy. Adv Differ Equ. (2016) 2016:225. doi: 10.1186/s13662-016-0952-x

Crossref Full Text | Google Scholar

59. Yan C, Wang W. Modeling HIV dynamics under combination therapy with inducers and antibodies. Bull Math Biol. (2019) 81:2625–48. doi: 10.1007/s11538-019-00621-0

PubMed Abstract | Crossref Full Text | Google Scholar

60. Allali K, Danane J, Kuang Y. Global analysis for an HIV infection model with CTL immune response and infected cells in eclipse phase. Appl Sci. (2017) 7:861. doi: 10.3390/app7080861

Crossref Full Text | Google Scholar

61. Sun Q, Min L, Kuang Y. Global stability of infection-free state and endemic infection state of a modified human immunodeficiency virus infection model. IET Syst Biol. (2015) 9:95–103. doi: 10.1049/iet-syb.2014.0046

PubMed Abstract | Crossref Full Text | Google Scholar

62. Wang X, Rong L. HIV low viral load persistence under treatment: insights from a model of cell-to-cell viral transmission. Appl Math Lett. (2019) 94:44–51. doi: 10.1016/j.aml.2019.02.019

Crossref Full Text | Google Scholar

63. Burini D, Knopoff D. Epidemics and society-a multiscale vision from the small world to the globally interconnected world. Math Models Methods Appl Sci. (2024) 34:1564–96. doi: 10.1142/S0218202524500295

Crossref Full Text | Google Scholar

64. Boulaaras S, Rehman ZU, Abdullah FA, Jan R, Abdullah M, Jan A, et al. Coronavirus dynamics, infections and preventive interventions using fractional-calculus analysis. AIMS Math. (2023) 8:8680–701. doi: 10.3934/math.2023436

Crossref Full Text | Google Scholar

65. Tang TQ, Shah Z, Jan R, Deebani W, Shutaywi M. A robust study to conceptualize the interactions of CD4+ T-cells and human immunodeficiency virus via fractional-calculus. Phys Scr. (2021) 96:125231. doi: 10.1088/1402-4896/ac2d7b

Crossref Full Text | Google Scholar

66. Hattaf K. A new mixed fractional derivative with applications in computational biology. Computation. (2024) 12:7. doi: 10.3390/computation12010007

Crossref Full Text | Google Scholar

67. Attaullah R, Jan S, YüzbaşI Ş. Dynamical behaviour of HIV Infection with the influence of variable source term through Galerkin method. Chaos Solitons Fractals. (2021) 152:11429. doi: 10.1016/j.chaos.2021.111429

Crossref Full Text | Google Scholar

Keywords: HIV-1, immune impairment, global stability, delays, Lyapunov stability, sensitivity analysis

Citation: AlShamrani NH, Halawani RH and Elaiw AM (2024) Stability of generalized models for HIV-1 dynamics with impaired CTL immunity and three pathways of infection. Front. Appl. Math. Stat. 10:1412357. doi: 10.3389/fams.2024.1412357

Received: 04 April 2024; Accepted: 15 July 2024;
Published: 14 August 2024.

Edited by:

Fahad Al Basir, Asansol Girls' College, India

Reviewed by:

Rashid Jan, University of Swabi, Pakistan
Zakaria Yaagoub, University of Hassan II Casablanca, Morocco
Nicola Bellomo, Polytechnic University of Turin, Italy

Copyright © 2024 AlShamrani, Halawani and Elaiw. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Noura H. AlShamrani, bmhhbHNoYW1yYW5pJiN4MDAwNDA7dWouZWR1LnNh

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.