# Global Dynamics of a Within-Host COVID-19/AIDS Coinfection Model with Distributed Delays

Authors: S. A. Azoz; A. M. Elaiw; E. Ramadan; A. D. Al Agha; Aeshah A. Raezah
Journal: Journal of Mathematics
 (2022)

---

## Abstract

Acquired immunodeficiency syndrome (AIDS) is a spectrum of conditions caused by infection with the human immunodeficiency virus (HIV). Among people with AIDS, cases of COVID-19 have been reported in many countries. COVID-19 (coronavirus disease 2019) is caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). In this manuscript, we are going to present a within-host COVID-19/AIDS coinfection model to study the dynamics and influence of the coinfection between COVID-19 and AIDS. The model is a six-dimensional delay differential equation that describes the interaction between uninfected epithelial cells, infected epithelial cells, free SARS-CoV-2 particles, uninfected CD4+ T cells, infected CD4+ T cells, and free HIV-1 particles. We demonstrated that the proposed model is biologically acceptable by proving the positivity and boundedness of the model solutions. The global stability analysis of the model is carried out in terms of the basic reproduction number. Numerical simulations are carried out to investigate that if COVID-19/AIDS coinfected individuals have a poor immune response or a low number of CD4+ T cells, then the viral load of SARS-CoV-2 and the number of infected epithelial cells will rise. On the contrary, the existence of time delays can rise the number of uninfected CD4+ T cells and uninfected epithelial cells, thus reducing the viral load within the host.


---

## Body

## 1. Introduction

In December 2019, the first case of the emergence of the severe acute respiratory syndrome coronavirus 2 (COVID-19) occurred in Wuhan, China. In March 2020, the World Health Organization (WHO) declared COVID-19 a worldwide epidemic. Globally, as of 27 April 2022, over 500 million people were infected with COVID-19, including 6 million deaths [1]. Old age and its accompanying symptoms such as diabetes, heart disease, and high blood pressure are considered risk factors for developing severe COVID-19 infection and are associated with a high death rate [2, 3]. Some other risk factors are associated when infection with COVID-19 occurs in people with chronic diseases such as acquired immunodeficiency syndrome (AIDS) [4]. In 2020, there were 37.7 million persons living with HIV-1 (PLWH) worldwide; HIV-1 causes acquired immunodeficiency syndrome (AIDS) with 680,000 of them dying from HIV-1-related diseases, and only 73% of them were on antiretroviral medication (ART) [5]. Because their immune systems are impaired, PLWH who do not receive ART or whose condition is poorly managed could be more susceptible to developing COVID-19. If infected with COVID-19, such people are at a greater risk of developing acute symptoms and dying. The coinfection cases are challenging due to the scarcity of data on the outcomes and consequences of SARS-CoV-2 infection in HIV-1 positive individuals [3, 6, 7].HIV-1 and SARS-CoV-2 are both RNA viruses. SARS-CoV-2 attacks upper respiratory epithelial cells, and the virus generated by infected cells goes down to the lower airway, infecting bronchial and alveolar epithelial cells [4, 8]. On the other hand, HIV-1 targets CD4+ lymphocytes, which are the immune system’s most plentiful white blood cells (referred to as CD4+ T cells). A great effort is being made in many areas of the world to create measures to battle these viruses and study their biological and immunological features and clinical outcomes. Some of these studies indicate that COVID-19 pandemic has caused disruptions in HIV-1 care facilities in many countries [9, 10]. However, it is unclear whether people infected with HIV-1 having an increased incidence of COVID-19 and significant clinical signs, despite a controversial suggestion that antiretroviral therapy or HIV-1-related immunosuppression could protect HIV-1 infected people from severe COVID-19. A number of HIV-1 and SARS-CoV-2 coinfection cases have been documented throughout the world [11, 12]. Most studies of COVID-19/AIDS coinfection reported that there is a lack of clarity on what constitutes the primary illness and what constitutes comorbidity in the context of coinfection. Few studies inferred that SARS-CoV-2 infection does not increase the course of HIV-1 infection in PLWH nor does HIV-1 infection have an impact on COVID-19 infection course in PLWH [13–15]. However, Wang et al. [16] published a case report of an HIV-1/COVID-19 patient with such a lower CD4+ T cell number, and as a result, the patient had a prolonged COVID-19 course and decreased antibody levels. Moreover, COVID-19/AIDS coinfection has been observed to cause pneumonia problems more frequently than COVID-19 alone [17]. This study aims to give a comprehensive picture of SARS-CoV-2 infection in persons having HIV-1/AIDS.Mathematical models that consist of a system of differential equations have proven their effectiveness in studying the interactions between viruses and their hosts and the common interactions between diseases (see e.g., [12, 18–24]). HIV within-host models have been widely investigated and great results have been reached [18, 19, 25–28]. On the other side, SARS-CoV-2 within-host modeling has received less attention ([24, 29–32]). Some coinfection models between SARS-CoV-2 and other viruses have been developed. For example, Pinky and Dobrovolny [33] used a within-host model to investigate SARS-CoV-2 coinfections with several viruses types such as influenza A virus (IAV), parainfluenza virus (PIV), and human rhinovirus (HRV). In fact, the models of coinfection are essential to grasp the coinfection dynamics between SARS-CoV-2 and HIV, to assist the experimental studies and save time, and to develop effective treatments for coinfected people. Ahmed et al. [34] created a fractional epidemiological model to analyze the pandemic scenario in numerous HIV and COVID-19 affected countries, including South Africa and Brazil. Then, to the best of our knowledge, the first ordinary differential within-host SARS-CoV-2/HIV coinfection system is presented by Al Agha et al. [20]. The formulation of their model is based on Nowak and Bangham’s model that was used widely to model HIV monoinfection and SARS-CoV-2 monoinfection. Al Agha et al. used the same principals to model SARS-CoV-2/HIV coinfection and connect the two infections together. The model is formulated as follows:(1)X˙t=ρ−αXt−ηXtVt,Y˙t=ηXtVt−kYt−μYtSt,V˙t=aYt−δVt,S˙t=ξ+uYtSt−γSt−θStHt,W˙t=θStHt−βWt,H˙t=λWt−ωHt,where Xt, Yt, and Vt represent the healthy epithelial cells, infected epithelial cells, and free SARS-CoV-2 particles, respectively, whilst St, Wt, and Ht depict healthy CD4+ T cells, infected CD4+ T cells, and free HIV particles concentrations at time t, respectively. Epithelial cells are recruited at rate ρ and turned into infected cells at pace ηXtVt. Infected produce SARS-CoV-2 particles at rate aYt. CD4+ T cells are recruited at rate ξ, eliminate infected epithelial cells at a proportion μYtSt, and proliferate at rate uYtSt. HIV particles infect CD4+ T cells at rate θStHt. The infected cells produces HIV at rate λWt. The components Xt, Yt, Yt, Yt, Wt, and Ht die at rates αXt, kYt, δVt, γSt, βWt, and ωHt, respectively. Then, Elaiw et al. [21] adopted the same previous model with the addition of the effect of latent cells, and then, they made a comprehensive study of the proposed model. Ringa et al. [22] presented a new mathematical model for COVID-19 and HIV/AIDS. The dynamics of the full model is driven by that of its submodels. Also, they studied the impact of intervention measures by incorporating it into the model using time-dependent controls.Most of the previous publications are the assumptions that cells produce viruses immediately after they are infected. It is commonly observed that in many biological processes, a time delay is inevitable. For HIV-1 infection, it roughly takes about one day for a newly infected cell to become productive and then to be able to produce new virus particles. Therefore, mathematicians have frequently used different types of delays to make biological models more realistic. In [26, 28, 35], HIV models with time delay were introduced, whilst modeling and analysis of COVID-19 based on a time delay dynamic model are presented in [12, 23]. Although there are some publications that combine the coinfection between viruses in the presence of time delay, there are still no models of coinfection between SARS-CoV-2 and HIV with time delay. Due to the decisive role of time delays in dynamic systems, the objective of this work is to expand model (1) to accommodate distributed delays. This can help comprehend the coinfection dynamics between SARS-CoV-2 and HIV-1 from a different perspective. A continuous distribution function is used to represent the delay in case of distributed time delay. This makes distributed delays more realistic than discrete time delays which presume that each individual in the population has the same delay period. Thus, we have investigated a model with six delay differential equations, and we have established the solutions nonnegativity and boundedness, listed the prospective equilibrium points and the conditions of existence, discussed the global stability of the equilibria, and examined time delay impact on the model’s dynamics. The document includes the following sections: the model is presented in Section 2. Section 3 confirms the basic properties of the model. Section 4 exhibits the global properties of the model. Section 5 lists the numerical simulations. Finally, Section 6 debates the results and some potential next directions.


## 2. COVID-19/AIDS Coinfection Model with Distributed Delay

In this section, we extend model (1) by considering a variety of distributed time delays as follows:(2)X˙t=ρ−αXt−ηXtVt,Y˙t=η∫0∞g1ϵe−m1ϵXt−ϵVt−ϵdϵ−kYt−υYtSt,V˙t=a∫0∞g2ϵe−m2ϵYt−ϵdϵ−℘Vt,S˙t=ξ+uYtSt−γSt−ℑStHt,W˙t=ℑ∫0∞g3ϵe−m3ϵSt−ϵHt−ϵdϵ−βWt,H˙t=λ∫0∞g4ϵe−m4ϵWt−ϵdϵ−ωHt.Thus, we have a system of six delay differential equations whereXt, Yt, Vt, St, Wt, and Ht stand for the concentrations of uninfected epithelial cells, infected epithelial cells, free SARS-CoV-2 particles, uninfected CD4+ T cells, infected CD4+ T cells, and HIV-1 particles at time t, respectively. A scheme describing the coinfection between SARS-CoV-2 and HIV in host without time delay is shown in Figure 1. The factor g1εe−m1ε designates the likelihood that uninfected epithelial cells were in touch with SARS-CoV-2 particles at time t−ϵ survived ϵ time units, and infection occurs at time t. The term g2ϵe−m2ϵ simulates the probability of new immature SARS-CoV-2 particles at time t−ϵ survived ϵ time units and mature at time t. Moreover, the factor g3ϵe−m3ϵ symbolizes the probability that uninfected CD4+ T cells contacted by HIV-1 particles at time t−ϵ survived ϵ time units and become infected at time t. The term g4ϵe−m4ϵ represents the probability that new immature HIV-1 particles at time t−ϵ persisted ϵ time units and mature at time t, where miandi=1,2,3,4, and are the positive constants. The delay parameter ϵ is a random variable picked from probability distribution functions giϵ during time interval 0,∞. The functions giϵ (i=1,2,3,4) satisfy giϵ>0 and(3)∫0∞giϵdϵ=1,∫0∞giϵe−nϵdϵ<∞,where n>0. Let us denote the following model:(4)ℒ¯iϵ=giϵe−miϵ,Li=∫0∞ℒ¯iϵdϵ,where i=1,2,3,4. This implies that 0<Li≤1. The initial conditions of model (2) are specified as follows:(5)Xϖ=φ1ϖ,Yϖ=φ2ϖ,Vϖ=φ3ϖ,Sϖ=φ4ϖ,Wϖ=φ5ϖ,Hϖ=φ6ϖ,φjϖ≥0,ϖ∈−∞,0,j=1,2,…6,where φjϖ∈ℂ=ζ∈C−∞,0,ℝ:ζθeαθis uniformly continuous forθ∈−∞,0,ζ<∞, and ζ=supθ≤0ζθeαθ such that α is a positive constant. Here, ℂ is the Banach space of fading memory type [36]. Therefore, using the standard theory of differential equations with infinitely distributed delays [37, 38], model (2) with initial constraints (3) has a single solution.Figure 1
Scheme describing the coinfection between SARS-CoV-2 and HIV.


## 3. Basic Characteristics

This section proves that model (2) solutions are non-negative and ultimately bounded. Additionally, it computes whole potential equilibria and the threshold numbers.


### 3.1. Non-Negativity and Boundedness

Proposition 1.
All of model (2) solutions with beginning conditions (3) are non-negative and eventually bounded.Proof.
Starting with model (2) first equation, we obtain X˙t|X=0=ρ>0, which yields that Xt>0 for all t≥0. From fourth equation of the model, we get S˙t|S=0=ξ>0; then, St>0 for all t≥0. Furthermore, the rest of the model equations give us the following model:(6)Yt=φ20e−∫0tk+υsqdq+η∫0te−∫ℓtk+υsqdq∫0∞ℒ¯1ϵXℓ−ϵVℓ−ϵdϵdℓ≥0,Vt=φ30e−℘t+a∫0te−℘t−ℓ∫0∞ℒ¯2ϵYℓ−ϵdϵdℓ≥0,Wt=φ50e−βt+ℑ∫0te−βt−ℓ∫0∞ℒ¯3ϵSℓ−ϵHℓ−ϵdϵdℓ≥0,Ht=φ60e−ωt+λ∫0te−ωt−ℓ∫0∞ℒ¯4ϵWℓ−ϵdϵdℓ≥0.
For allt∈0,∞, as a result of the recursive argument, we obtain Xt,Yt,Vt,St,Wt,Ht≥0 for all t≥0. Hence, system (2) solutions with initial conditions (3) realize Xt,Yt,Vt,St,Wt,Ht∈ℝ≥06 for all non-negative values of t.
Now, we establish the boundedness of the model’s solutions. Based on model (2) first equation, we gain limt⟶∞supXt≤Ω1, where Ω1=ρ/α. We define the following model:(7)Ψ1t=∫0∞ℒ¯1εXt−εdε+Yt+υuSt.
Then, we get the following model:(8)Ψ˙1t=∫0∞ℒ¯1ερ−αXt−ε−ηXt−εVt−εdε+η∫0∞ℒ¯1εXt−εVt−εdε−kYt−υYtSt+υuξ+uYtSt−γSt−ℑStHt≤ρ∫0∞ℒ¯1εdε+υξu−α∫0∞ℒ¯1εXt−εdε−kYt−υγuSt≤ρ+υξu−ϕ1∫0∞ℒ¯1εXt−εdε+Yt+υuSt=ρ+υξu−ϕ1Ψ1t,where ϕ1=minα,k,γ. This implies that limt⟶∞supΨ1t≤Ω2, where Ω2=ρ/ϕ1+υξ/uϕ1. Since Yt and St are non-negative, then limt⟶∞supYt≤Ω2 and limt⟶∞supSt≤Ω3, where Ω3=uΩ2/υ. Using model (2) third equation, we get the following model:(9)V˙t=a∫0∞ℒ¯2ϵYt−ϵdϵ−℘Vt≤aL2Ω2−℘Vt≤aΩ2−℘Vt.
This implies thatlimt⟶∞supVt≤Ω4, where Ω4=aΩ2/℘.
We define the following model:(10)Ψ2t=∫0∞ℒ¯3εSt−εdε+Wt.
Following that, we get the following model:(11)Ψ˙2t=∫0∞ℒ¯3εξ+uYt−εSt−ε−γSt−ε−ℑSt−εHt−εdε+ℑ∫0∞ℒ¯3εSt−εHt−εdε−βWt≤ξ+uΩ2Ω3−ϕ2Ψ2t,where ϕ2=minγ,β. Thus, we have limt⟶∞supWt≤Ω5, whereΩ5=ξ+uΩ2Ω3/ϕ2. Finally, the last equation of model (2) gives the following model:(12)H˙t=λ∫0∞ℒ¯4ϵWt−ϵdϵ−ωHt≤λΩ5−ωHt.
Thus,limt⟶∞supHt≤Ω6, where Ω6=λΩ5/ω. Accordingly, the following region is positively invariant with regard to model (2).(13)Ξ=X,Y,V,S,W,H∈C≥06:X≤Ω1,Y≤Ω2,V≤Ω4,S≤Ω3,W≤Ω5,H≤Ω6.


### 3.2. Equilibrium Points

This subsection displays all of model (2) possible equilibrium states and deduces four threshold parameters that determine the equilibria existence.We solve the following set of algebraic equations to calculate the following model equilibrium points:(14)0=ρ−αX−ηXV,0=ηL1XV−kY−υYS,0=aL2Y−℘V,0=ξ+uYS−γS−ℑSH,0=ℑL3SH−βW,0=λL4W−ωH.From last equation of model (4), we have W=ωH/λL4. Then, we substitute in the fifth equation, and we get the following model:(15)0=ℑL3S−βωλL4H.So, we have two possibilities:(16)H=0orS=βωλℑL3L4.Then, doing the same for second and third equations, we obtain another two possibilities as follows:(17)V=0orηL1X−kρaL2−υρSaL2=0.Equations (16) and (17) provide us with four possibilities. Accordingly, model (4) has four equilibrium points:(i)
Uninfected equilibriumEP0=X0,0,0,S0,0,0, where X0=ρ/α and S0=ξ/γ(ii)
The HIV-1 monoinfection equilibriumEPH=X1,0,0,S1,W1,H1, where(18)X1=ρα,S1=βωℑλL3L4,W1=γωℑλL4ξℑλL3L4βγω−1,H1=γℑξℑλL3L4βγω−1.It follows thatW1>0 and H1>0 only when ξℑλL3L4/βγω>1. Thus, we have the following model:(19)X1=X0,S1=S0R1,W1=γωℑλL4R1−1,H1=γℑR1−1.Here,R1=ξℑλL3L4/βγω. Here, R1 is the basic reproduction number for HIV-1 infection. It sets start of HIV-1 infection in host body. We note that W1>0 and H1>0 if R1>1. Therefore, EPH exists when R1>1.(iii)
SARS-CoV-2 monoinfection equilibriumEPV=X2,Y2,V2,S2,0,0, where(20)X2=Y2k+S2Y2υηV2L1,Y2=℘V2aL2,S2=ξγ−uY2.V2 satisfies the following equation:(21)u℘2ηkV22+uα℘2k−aγ℘ηkL2−a℘ηυξL2−au℘ηρL1L2V2−aαγ℘kL2−aα℘υξL2+a2γηρL1L22aηL1L2aγL2−u℘V2=0.To prove that equation (21) has a positive root, we introduce a function BV as follows:(22)BV=u℘2ηkV2+uα℘2k−aγ℘ηkL2−a℘ηυξL2−au℘ηρL1L2V−aαγ℘kL2−aα℘υξL2+a2γηρL1L22aηL1L2aγL2−u℘V.Then, we have the following equation:(23)B0=−aαγ℘kL2−aα℘υξL2+a2γηρL1L22a2γηL1L22=αγ℘k+α℘υξaγηL1L2R2−1.Here,R2=aγηρL1L2/α℘γk+υξ. This implies that B0>0 when R2>1. In addition, we find that(24)limV⟶aγL2u℘−BV=−∞.It follows that there exists0<V2<aγL2/u℘ such that BV2=0. From equation (20), we get Y2>0, S2>0, and X2>0. As a result, we deduce that EPV exists when R2>1. Here, R2 is the basic reproduction number for SARS-CoV-2 infection. It defines start of SARS-CoV-2 infection in host body.(iv)
COVID-19/AIDS coinfection equilibriumEPVH=X3,Y3,V3,S3,W3,H3, where(25)X3=℘kℑλL3L4+βυωaηℑλL1L2L3L4,Y3=α℘aηL2aηℑλρL1L2L3L4α℘ℑkλL3L4+βυω−1,V3=αηaηℑλρL1L2L3L4α℘ℑkλL3L4+βυω−1,S3=βωℑλL3L4,W3=ωuα℘+aγηL2aηℑλL2L4λξβω+uλρL1ℑkλL3L4+βυωaηℑL2L3L4uα℘+aγηL2−1,H3=uα℘+aγηL2aηℑL2λξβω+uλρL1ℑkλL3L4+βυωaηℑL2L3L4uα℘+aγηL2−1.It follows thatW3>0 and H3>0 only when λξ/βω+uλρL1/ℑkλL3L4+βυωaηℑL2L3L4/uα℘+aγηL2>1. On the other hand, Y3>0 and V3>0 only when aηℑλρL1L2L3L4/α℘ℑkλL3L4+βυω>1.Thus, we can rewrite the components ofEPVH as follows:(26)X3=X0R4,Y3=α℘aηL2R4−1,V3=αηR4−1,S3=βωℑλL3L4,W3=ωuα℘+aγηL2aηℑλL2L4R3−1,H3=uα℘+aγηL2aηℑL2R3−1,where(27)R3=λξβω+uλρL1ℑkλL3L4+βυωaηℑL2L3L4uα℘+aγηL2,R4=aηℑλρL1L2L3L4α℘ℑkλL3L4+βυω.Therefore,EPVH exists when R3>1 and R4>1. Here, R3 and R4 are the threshold parameters that mark the COVID-19/AIDS coinfection incidence.The threshold parameters are defined as follows:(28)R1=ξℑλL3L4βγω,R2=aγηρL1L2α℘γk+υξ,R3=λξβω+uλρL1ℑkλL3L4+βυωaηℑL2L3L4uα℘+aγηL2,R4=aηℑλρL1L2L3L4α℘ℑkλL3L4+βυω.For simplicity, the contractions listed will be used in the parts that follow(29)Xt≡X,Yt≡Y,Vt≡V,St≡S,Wt≡W,Ht≡H,and(30)Xt−ϵ≡Xϵ,Yt−ϵ≡Yϵ,Vt−ϵ≡Vϵ,St−ϵ≡Sϵ,Wt−ϵ≡Wϵ,Ht−ϵ≡Hϵ.


## 3.1. Non-Negativity and Boundedness

Proposition 1.
All of model (2) solutions with beginning conditions (3) are non-negative and eventually bounded.Proof.
Starting with model (2) first equation, we obtain X˙t|X=0=ρ>0, which yields that Xt>0 for all t≥0. From fourth equation of the model, we get S˙t|S=0=ξ>0; then, St>0 for all t≥0. Furthermore, the rest of the model equations give us the following model:(6)Yt=φ20e−∫0tk+υsqdq+η∫0te−∫ℓtk+υsqdq∫0∞ℒ¯1ϵXℓ−ϵVℓ−ϵdϵdℓ≥0,Vt=φ30e−℘t+a∫0te−℘t−ℓ∫0∞ℒ¯2ϵYℓ−ϵdϵdℓ≥0,Wt=φ50e−βt+ℑ∫0te−βt−ℓ∫0∞ℒ¯3ϵSℓ−ϵHℓ−ϵdϵdℓ≥0,Ht=φ60e−ωt+λ∫0te−ωt−ℓ∫0∞ℒ¯4ϵWℓ−ϵdϵdℓ≥0.
For allt∈0,∞, as a result of the recursive argument, we obtain Xt,Yt,Vt,St,Wt,Ht≥0 for all t≥0. Hence, system (2) solutions with initial conditions (3) realize Xt,Yt,Vt,St,Wt,Ht∈ℝ≥06 for all non-negative values of t.
Now, we establish the boundedness of the model’s solutions. Based on model (2) first equation, we gain limt⟶∞supXt≤Ω1, where Ω1=ρ/α. We define the following model:(7)Ψ1t=∫0∞ℒ¯1εXt−εdε+Yt+υuSt.
Then, we get the following model:(8)Ψ˙1t=∫0∞ℒ¯1ερ−αXt−ε−ηXt−εVt−εdε+η∫0∞ℒ¯1εXt−εVt−εdε−kYt−υYtSt+υuξ+uYtSt−γSt−ℑStHt≤ρ∫0∞ℒ¯1εdε+υξu−α∫0∞ℒ¯1εXt−εdε−kYt−υγuSt≤ρ+υξu−ϕ1∫0∞ℒ¯1εXt−εdε+Yt+υuSt=ρ+υξu−ϕ1Ψ1t,where ϕ1=minα,k,γ. This implies that limt⟶∞supΨ1t≤Ω2, where Ω2=ρ/ϕ1+υξ/uϕ1. Since Yt and St are non-negative, then limt⟶∞supYt≤Ω2 and limt⟶∞supSt≤Ω3, where Ω3=uΩ2/υ. Using model (2) third equation, we get the following model:(9)V˙t=a∫0∞ℒ¯2ϵYt−ϵdϵ−℘Vt≤aL2Ω2−℘Vt≤aΩ2−℘Vt.
This implies thatlimt⟶∞supVt≤Ω4, where Ω4=aΩ2/℘.
We define the following model:(10)Ψ2t=∫0∞ℒ¯3εSt−εdε+Wt.
Following that, we get the following model:(11)Ψ˙2t=∫0∞ℒ¯3εξ+uYt−εSt−ε−γSt−ε−ℑSt−εHt−εdε+ℑ∫0∞ℒ¯3εSt−εHt−εdε−βWt≤ξ+uΩ2Ω3−ϕ2Ψ2t,where ϕ2=minγ,β. Thus, we have limt⟶∞supWt≤Ω5, whereΩ5=ξ+uΩ2Ω3/ϕ2. Finally, the last equation of model (2) gives the following model:(12)H˙t=λ∫0∞ℒ¯4ϵWt−ϵdϵ−ωHt≤λΩ5−ωHt.
Thus,limt⟶∞supHt≤Ω6, where Ω6=λΩ5/ω. Accordingly, the following region is positively invariant with regard to model (2).(13)Ξ=X,Y,V,S,W,H∈C≥06:X≤Ω1,Y≤Ω2,V≤Ω4,S≤Ω3,W≤Ω5,H≤Ω6.


## 3.2. Equilibrium Points

This subsection displays all of model (2) possible equilibrium states and deduces four threshold parameters that determine the equilibria existence.We solve the following set of algebraic equations to calculate the following model equilibrium points:(14)0=ρ−αX−ηXV,0=ηL1XV−kY−υYS,0=aL2Y−℘V,0=ξ+uYS−γS−ℑSH,0=ℑL3SH−βW,0=λL4W−ωH.From last equation of model (4), we have W=ωH/λL4. Then, we substitute in the fifth equation, and we get the following model:(15)0=ℑL3S−βωλL4H.So, we have two possibilities:(16)H=0orS=βωλℑL3L4.Then, doing the same for second and third equations, we obtain another two possibilities as follows:(17)V=0orηL1X−kρaL2−υρSaL2=0.Equations (16) and (17) provide us with four possibilities. Accordingly, model (4) has four equilibrium points:(i)
Uninfected equilibriumEP0=X0,0,0,S0,0,0, where X0=ρ/α and S0=ξ/γ(ii)
The HIV-1 monoinfection equilibriumEPH=X1,0,0,S1,W1,H1, where(18)X1=ρα,S1=βωℑλL3L4,W1=γωℑλL4ξℑλL3L4βγω−1,H1=γℑξℑλL3L4βγω−1.It follows thatW1>0 and H1>0 only when ξℑλL3L4/βγω>1. Thus, we have the following model:(19)X1=X0,S1=S0R1,W1=γωℑλL4R1−1,H1=γℑR1−1.Here,R1=ξℑλL3L4/βγω. Here, R1 is the basic reproduction number for HIV-1 infection. It sets start of HIV-1 infection in host body. We note that W1>0 and H1>0 if R1>1. Therefore, EPH exists when R1>1.(iii)
SARS-CoV-2 monoinfection equilibriumEPV=X2,Y2,V2,S2,0,0, where(20)X2=Y2k+S2Y2υηV2L1,Y2=℘V2aL2,S2=ξγ−uY2.V2 satisfies the following equation:(21)u℘2ηkV22+uα℘2k−aγ℘ηkL2−a℘ηυξL2−au℘ηρL1L2V2−aαγ℘kL2−aα℘υξL2+a2γηρL1L22aηL1L2aγL2−u℘V2=0.To prove that equation (21) has a positive root, we introduce a function BV as follows:(22)BV=u℘2ηkV2+uα℘2k−aγ℘ηkL2−a℘ηυξL2−au℘ηρL1L2V−aαγ℘kL2−aα℘υξL2+a2γηρL1L22aηL1L2aγL2−u℘V.Then, we have the following equation:(23)B0=−aαγ℘kL2−aα℘υξL2+a2γηρL1L22a2γηL1L22=αγ℘k+α℘υξaγηL1L2R2−1.Here,R2=aγηρL1L2/α℘γk+υξ. This implies that B0>0 when R2>1. In addition, we find that(24)limV⟶aγL2u℘−BV=−∞.It follows that there exists0<V2<aγL2/u℘ such that BV2=0. From equation (20), we get Y2>0, S2>0, and X2>0. As a result, we deduce that EPV exists when R2>1. Here, R2 is the basic reproduction number for SARS-CoV-2 infection. It defines start of SARS-CoV-2 infection in host body.(iv)
COVID-19/AIDS coinfection equilibriumEPVH=X3,Y3,V3,S3,W3,H3, where(25)X3=℘kℑλL3L4+βυωaηℑλL1L2L3L4,Y3=α℘aηL2aηℑλρL1L2L3L4α℘ℑkλL3L4+βυω−1,V3=αηaηℑλρL1L2L3L4α℘ℑkλL3L4+βυω−1,S3=βωℑλL3L4,W3=ωuα℘+aγηL2aηℑλL2L4λξβω+uλρL1ℑkλL3L4+βυωaηℑL2L3L4uα℘+aγηL2−1,H3=uα℘+aγηL2aηℑL2λξβω+uλρL1ℑkλL3L4+βυωaηℑL2L3L4uα℘+aγηL2−1.It follows thatW3>0 and H3>0 only when λξ/βω+uλρL1/ℑkλL3L4+βυωaηℑL2L3L4/uα℘+aγηL2>1. On the other hand, Y3>0 and V3>0 only when aηℑλρL1L2L3L4/α℘ℑkλL3L4+βυω>1.Thus, we can rewrite the components ofEPVH as follows:(26)X3=X0R4,Y3=α℘aηL2R4−1,V3=αηR4−1,S3=βωℑλL3L4,W3=ωuα℘+aγηL2aηℑλL2L4R3−1,H3=uα℘+aγηL2aηℑL2R3−1,where(27)R3=λξβω+uλρL1ℑkλL3L4+βυωaηℑL2L3L4uα℘+aγηL2,R4=aηℑλρL1L2L3L4α℘ℑkλL3L4+βυω.Therefore,EPVH exists when R3>1 and R4>1. Here, R3 and R4 are the threshold parameters that mark the COVID-19/AIDS coinfection incidence.The threshold parameters are defined as follows:(28)R1=ξℑλL3L4βγω,R2=aγηρL1L2α℘γk+υξ,R3=λξβω+uλρL1ℑkλL3L4+βυωaηℑL2L3L4uα℘+aγηL2,R4=aηℑλρL1L2L3L4α℘ℑkλL3L4+βυω.For simplicity, the contractions listed will be used in the parts that follow(29)Xt≡X,Yt≡Y,Vt≡V,St≡S,Wt≡W,Ht≡H,and(30)Xt−ϵ≡Xϵ,Yt−ϵ≡Yϵ,Vt−ϵ≡Vϵ,St−ϵ≡Sϵ,Wt−ϵ≡Wϵ,Ht−ϵ≡Hϵ.


## 4. Global Properties

We demonstrate the global asymptotic stability of all equilibria in this section by building Lyapunov functions using the approach described in [39]. We define FΔ=Δ−1−lnΔ, where Δ can be any variable for the model.Theorem 1.
Globally asymptotically stable (G.A.S) of equilibriumEP0 is satisfied when R1≤1 and R2≤1.Proof.
Take a Lyapunov functionϑ0X,Y,V,S,W,H as follows:(31)ϑ0=X0FXX0+1L1Y+ηX0℘V+υuL1S0FSS0+υuL1L3W+υβuλL1L3L4H+ηL1∫0∞ℒ¯1ε∫t−εtXℓVℓdℓdε+aηX0℘∫0∞ℒ¯2ε∫t−εtYℓdℓdε+υℑuL1L3∫0∞ℒ¯3ε∫t−εtSℓHℓdℓdε+υβuL1L3L4∫0∞ℒ¯4ℓ∫t−εtWℓdℓdε.
Clearly,ϑ0X,Y,V,S,W,H>0 for all X,Y,V,S,W,H>0 and ϑ0X0,0,0,S0,0,0=0. Calculating dϑ0/dt along the solutions of system (2) gives the following equation:(32)dϑ0dt=1−X0Xρ−αX−ηXV+1L1η∫0∞ℒ¯1εXεVεdε−kY−υYS+ηX0℘a∫0∞ℒ¯2εYεdε−℘V+υuL11−S0Sξ+uYS−γS−ℑSH+υuL1L3ℑ∫0∞ℒ¯3εSεHεdε−βW+υβuλL1L3L4λ∫0∞ℒ¯4εWεdε−ωH+ηL1∫0∞ℒ¯1εXV−XεVεdε+aηX0℘∫0∞ℒ¯2εY−Yεdε+υℑuL1L3∫0∞ℒ¯3εSH−SεHεdε+υβuL1L3L4∫0∞ℒ¯4εW−Wεdε.
Adding up the terms in equation (32), we obtain the following equation:(33)dϑ0dt=1−X0Xρ−αX−ηXV+ηX0V−1L1kY−1L1υYS−ηX0V+υuL11−S0Sξ−γS+1L1υYS−1L1υYS0−υuL1ℑSH+υuL1ℑS0H−υβuL1L3W−υβωuλL1L3L4H+ηXV+aηX0℘L2Y+υuL1ℑSH+υβuL1L3W.
Usingρ=αX0 and ξ=γS0, we obtain the following equation:(34)dϑ0dt=−αXX−X02−υγuL1SS−S02+ηX0℘aL2−1L1k−1L1υS0Y+υuL1ℑS0−βωλL3L4H=−αXX−X02−υγuL1SS−S02+γk+υξγL1aγηρL1L2α℘γk+υξ−1Y+υβωuλL1L3L4ξℑλL3L4βγω−1H=−αXX−X02−υγuL1SS−S02+γk+υξγL1R2−1Y+υβωuλL1L3L4R1−1H.
SinceR1≤1 and R2≤1, we get dϑ0/dt≤0 for all X,Y,V,S,W,H>0. Also, dϑ0/dt=0 when X=X0, S=S0, and Y=H=0. Set T0=X,Y,V,S,W,H:dϑ0/dt=0 and the largest invariant subset (L.I.S) of T0 by T0′. Then, the model solutions converge to T0′. The set T0′ contains elements with Xt=X0, St=S0, and Yt=Ht=0, and hence, Y˙t=H˙t=0. The second and last equations of the model (2) give the following equation:(35)0=Y˙t=η∫0∞ℒ¯1ϵX0Vϵdϵ,0=H˙t=λ∫0∞ℒ¯4ϵWϵdϵ.
Thus, we getVt=Wt=0 for all t. Then, T0′=EP0, and using Lyapunov–LaSalle asymptotic stability theorem [40–42], EP0 is G.A.S.In the following theorems, we need to use the equalities:(36)lnSεHεSH=lnSεHεWjSjHjW+lnSjS+lnWHjWjH,lnWεW=lnWεHjWjH+lnWjHWHj,j=1,3.Furthermore,(37)lnXϵVϵXV=lnXϵVϵYiXiViY+lnXiX+lnYViYiV,lnYϵY=lnYϵViYiV+lnYiVYVi,i=2,3.Theorem 2.
IfR1>1 and R4≤1, the equilibrium EPH is G.A.S.Proof.
Consider a Lyapunov functionϑ1X,Y,V,S,W,H as follows:(38)ϑ1=X1FXX1+1L1Y+ηX1℘V+υuL1S1FSS1+υuL1L3W1FWW1+υβuλL1L3L4H1FHH1+ηL1∫0∞ℒ¯1ε∫t−εtXℓVℓdℓdε+aηX1℘∫0∞ℒ¯2ε∫t−εtYldℓdε+υℑS1H1uL1L3∫0∞ℒ¯3ε∫t−εtFSℓHℓS1H1dℓdε+υβW1uL1L3L4∫0∞ℒ¯4ε∫t−εtFWℓW1dℓdε.
Differentiatingϑ1, we obtain the following equation:(39)dϑ1dt=1−X1Xρ−αX−ηXV+1L1η∫0∞ℒ¯1εXεVεdε−kY−υYS+ηX1℘a∫0∞ℒ¯2εYεdε−℘V+υuL11−S1Sξ+uYS−γS−ℑSH+υuL1L31−W1Wℑ∫0∞ℒ¯3εSεHεdε−βW+υβuλL1L3L41−H1Hλ∫0∞ℒ¯4εWεdε−ωH+ηL1∫0∞ℒ¯1εXV−XεVεdε+aηX1℘∫0∞ℒ¯2εY−Yεdε+υℑS1H1uL1L3∫0∞ℒ¯3εSHS1H1−SεHεS1H1+lnSεHεSHdε+υβW1uL1L3L4∫0∞ℒ¯4εWW1−WεW1+lnWεWdε.
Summing the terms of equation (39), we obtain the following equation:(40)dϑ1dt=1−X1Xρ−αX−1L1kY+υuL11−S1Sξ−γS−1L1υYS1+υuL1ℑS1H−υuL1L3ℑ∫0∞ℒ¯3εSεHεW1Wdε+υuL1L3βW1−υβuL1L3L4∫0∞ℒ¯4εWεH1Hdε−υβωuλL1L3L4H+υβωuλL1L3L4H1+aηX1℘L2Y+υℑS1H1uL1L3∫0∞ℒ¯3εlnSεHεSHdε+υβW1uL1L3L4∫0∞ℒ¯4εlnWεWdε.
By utilizing the equilibrium conditions forEPH, we get the following equation:(41)ρ=αX1,ξ=γS1+ℑH1S1,ℑL3H1S1=βW1,λL4W1=ωH1.
Then, we obtain the following equation:(42)dϑ1dt=−αXX−X12+ηX1℘aL2−1L1k−1L1υS1Y−υγuL1SS−S12+υuL1ℑS1H11−S1S+υuL1ℑS1H−υℑS1H1uL1L3∫0∞ℒ¯3εSεHεW1S1H1Wdε+υuL1ℑS1H1−υβW1uL1L3L4∫0∞ℒ¯4εWεH1W1Hdε−υuL1ℑS1H+υuL1ℑS1H1+υℑS1H1uL1L3∫0∞ℒ¯3εlnSεHεSHdε+υℑS1H1uL1L4∫0∞ℒ¯4εlnWεWdε.
Using the equalities given by equation (36) in case of j=1, we get the following equation:(43)dϑ1dt=−αXX−X12−υγuL1SS−S12+ℑkλL3L4+βυωℑλL1L3L4aηℑλρL1L2L3L4α℘ℑkλL3L4+βυω−1Y−υuL1ℑS1H1S1S−1−lnS1S−υℑS1H1uL1L3∫0∞ℒ¯3εSεHεW1S1H1W−1−lnSεHεW1S1H1Wdε−υℑS1H1uL1L4∫0∞ℒ¯4εWεH1W1H−1−lnWεH1W1Hdε.
Therefore, equation (43) becomes(44)dϑ1dt=−αXX−X12−υγuL1SS−S12+ℑkλL3L4+βυωℑλL1L3L4R4−1Y−υuL1ℑS1H1FS1S−υℑS1H1uL1L3∫0∞ℒ¯3εFSεHεW1S1H1Wdε−υℑS1H1uL1L4∫0∞ℒ¯4εFWεH1W1Hdε.
SinceR4≤1, we find that dϑ1/dt≤0 for all X,Y,V,S,W,H>0. Also, dϑ1/dt=0 when X=X1, S=S1, Y=0, W=W1, and H=H1. Model (2) solutions converge to T1′ the L.I.S of T1=X,Y,V,S,W,H:dϑ1/dt=0. The set T1′ contains elements with Yt=0; then, Yt=0. Second equation of system (2) implies(45)0=Y˙t=η∫0∞ℒ¯1ϵX1Vϵdϵ,which gives Vt=0 for all t. Therefore, T1′=EPH and EPH is G.A.S according to Lyapunov–LaSalle asymptotic stability theorem [40–42].Theorem 3.
IfR2>1 and R3≤1, then the equilibrium EPV is G.A.S.Proof.
We introduce a Lyapunov functionϑ2X,Y,V,S,W,H as follows:(46)ϑ2=X2FXX2+1L1Y2FYY2+ηX2℘V2FVV2+υuL1S2FSS2+υuL1L3W+υβuλL1L3L4H+ηX2V2L1∫0∞ℒ¯1ε∫t−εtFXℓVℓX2V2dℓdε+aηX2Y2℘∫0∞ℒ¯2ε∫t−εtFYℓY2dℓdε+υuL1L3ℑ∫0∞ℒ¯3ε∫t−εtSℓHℓdℓdε+υuL1L3L4β∫0∞ℒ¯4ε∫t−εtWℓdℓdε.
Differentiatingϑ2, we obtain the following equation:(47)dϑ2dt=1−X2Xρ−αX−ηXV+1L11−Y2Yη∫0∞ℒ¯1εXεVεdε−kY−υYS+ηX2℘1−V2V×a∫0∞ℒ¯2εYεdε−℘V+υuL11−S2Sξ+uYS−γS−ℑSH+υuL1L3ℑ∫0∞ℒ¯3εSεHεdε−βW+υβuλL1L3L4λ∫0∞ℒ¯4εWεdε−ωH+ηX2V2L1∫0∞ℒ¯1εXVX2V2−XεVεX2V2+lnXεVεXVdε+aηX2Y2℘∫0∞ℒ¯2εYY2−YεY2+lnYεYdε+υuL1L3ℑ∫0∞ℒ¯3εSH−SεHεdε+υβuL1L3L4∫0∞ℒ¯4εW−Wεdε.
By collecting the terms of equation (47), we have the following equation:(48)dϑ2dt=1−X2Xρ−αX−ηL1∫0∞ℒ¯1εXεVεY2Ydε−1L1kY+1L1kY2+1L1υY2S−ηX2℘a∫0∞ℒ¯2εYεV2Vdε+ηX2V2+υuL11−S2Sξ−γS−1L1υYS2+υuL1ℑS2H−υβωuλL1L3L4H+ηX2V2L1∫0∞ℒ¯1εlnXεVεXVdε+aηX2℘L2Y+aηX2Y2℘∫0∞ℒ¯2εlnYεYdε.
By using the equilibrium conditions forEPV,(49)ρ=αX2+ηX2V2,ηL1X2V2=kY2+υY2S2,aL2Y2=℘V2,ξ=γS2−uY2S2,we obtain the following equation:(50)dϑ2dt=−αXX−X22+ηV2X21−X2X+ηX2℘aL2−1L1k−1L1υS2Y−ηX2V2L1∫0∞ℒ¯1εXεVεY2X2V2Ydε+ηX2V2−1L1υY2S2+1L1υY2S−aηX2Y2℘∫0∞ℒ¯2εYεV2Y2Vdε+ηX2V2−υγuL1SS−S22−1L1υY2S21−S2S+υuL1ℑS2−βωλL3L4H+ηX2V2L1∫0∞ℒ¯1εlnXεVεXVdε+ηX2V2L2∫0∞ℒ¯2εlnYεYdε.
Using the equalities given by equation (37) in case of i=2, we get the following equation:(51)dϑ2dt=−αXX−X22−υγuL1SS−S22−1L1υY2S22−S2S−SS2−ηX2V2X2X−1−lnX2X+υuL1ℑS2−βωλL3L4H−ηX2V2L1∫0∞ℒ¯1εXεVεY2X2V2Y−1−lnXεVεY2X2V2Ydε−ηX2V2L2∫0∞ℒ¯2εYεV2Y2V−1−lnYεV2Y2Vdε.
Therefore, equation (51) becomes(52)dϑ2dt=−αXX−X22−υξuL1SS2S−S22−ηX2V2FX2X+υuL1ℑS2−βωλL3L4H−ηX2V2L1∫0∞ℒ¯1εFXεVεY2X2V2Ydε−ηX2V2L2∫0∞ℒ¯2εFYεV2Y2Vdε.
IfR3≤1, then EpVH does not exist since H3≤0 and W3≤0. This implies that(53)H˙t=λL4W−ωH≤0,W˙t=ℑL3SH−βW≤0.
Therefore, we getℑS−βω/λL3L4H≤0 for all H,S>0. Hence, we have ℑS2−βω/λL3L4≤0, and therefore, dϑ2/dt≤0 for all X,Y,V,S,W,H>0. In addition, dϑ2/dt=0 when X=X2, S=S2, H=0, Y=Y2, and V=V2. Solutions of the model (2) that converge to T2′ is the L.I.S of T2=X,Y,V,S,W,H:dϑ2/dt=0. The set T2′ has elements with Ht=0, and thus, H˙t=0. Using system (2) last equation, we get the following equation:(54)0=H˙t=λ∫0∞ℒ¯4ϵWϵdϵ.
YieldWt=0 for all values of t. Therefore, T2′=EPV and EPV is G.A.S according to Lyapunov–LaSalle asymptotic stability theorem [40–42].Theorem 4.
IfR4>1 and 1<R3≤1+aηℑλξL2L3L4/βωuα℘+aγηL2, then the equilibrium EPVH is G.A.S.Proof.
We consider a Lyapunov functionϑ3X,Y,V,S,W,H as follows:(55)ϑ3=X3FXX3+1L1Y3FYY3+ηX3℘V3FVV3+υuL1S3FSS3+υuL1L3W3FWW3+υβuλL1L3L4H3FHH3+ηX3V3L1∫0∞ℒ¯1ε∫t−εtFXℓVℓX3V3dℓdε+aηX3Y3℘∫0∞ℒ¯2ε∫t−εtFYℓY3dℓdε+υℑS3H3uL1L3∫0∞ℒ¯3ε∫t−εtFSℓHℓS3H3dℓdε+υβW3uL1L3L4∫0∞ℒ¯4ε∫t−εtFWℓW3dℓdε.
By differentiatingϑ3, we obtain the following equation:(56)dϑ3dt=1−X3Xρ−αX−ηXV+1L11−Y3Yη∫0∞ℒ¯1εXεVεdε−kY−υYS+ηX3℘1−V3V×a∫0∞ℒ¯2εYεdε−℘V+υuL11−S3Sξ+uYS−γS−ℑSH+υuL1L31−W3W×ℑ∫0∞ℒ¯3εSεHεdε−βW+υβuλL1L3L41−H3Hλ∫0∞ℒ¯4εWεdε−ωH+ηX3V3L1∫0∞ℒ¯1εXVX3V3−XεVεX3V3+lnXεVεXVdε+aηX3Y3℘∫0∞ℒ¯2εYY3−YεY3+lnYεYdε+υℑS3H3uL1L3∫0∞ℒ¯3εSHS3H3−SεHεS3H3+lnSεHεSHdε+υβW3uL1L3L4∫0∞ℒ¯4εWW3−WεW3+lnWεWdε.
Collecting terms of equation (56) gives the following equation:(57)dϑ3dt=1−X3Xρ−αX−ηL1∫0∞ℒ¯1εXεVεY3Ydε−1L1kY+1L1kY3+1L1υY3S−ηX3℘a∫0∞ℒ¯2εYεV3Vdε+ηX3V3+υuL11−S3Sξ−γS−1L1υYS3+υuL1ℑS3H−υuL1L3ℑ∫0∞ℒ¯3εSεHεW3Wdε+υuL1L3βW3−υβuL1L3L4∫0∞ℒ¯4εWεH3Hdε−υβuλL1L3L4ωH+υβuλL1L3L4ωH3+ηX3V3L1∫0∞ℒ¯1εlnXεVεXVdε+aηX3℘L2Y+aηX3Y3℘∫0∞ℒ¯2εlnYεYdε+υℑS3H3uL1L3∫0∞ℒ¯3εlnSεHεSHdε+υβW3uL1L3L4∫0∞ℒ¯4εlnWεWdε.
By using the equilibrium conditions forEPVH,(58)ρ=αX3+ηV3X3,ηL1V3X3=kY3+υY3S3,aL2Y3=℘V3,ξ=γS3+ℑH3S3−uY3S3,ℑL3H3S3=βW3,λL4W3=ωH3,we get the following equation:(59)dϑ3dt=−αXX−X32+ηV3X31−X3X+ηX3℘aL2−1L1k−1L1υS3Y−ηX3V3L1∫0∞ℒ¯1εXεVεY3X3V3Ydε+ηX3V3−1L1υY3S3+1L1υY3S−ηX3Y3℘a∫0∞ℒ¯2εYεV3Y3Vdε+ηX3V3−υγuL1SS−S32−1L1υY3S31−S3S+υuL1ℑS3H31−S3S+υuL1ℑS3H−υℑS3H3uL1L3∫0∞ℒ¯3εSεHεW3S3H3Wdε+υuL1ℑS3H3−υβW3uL1L3L4∫0∞ℒ¯4εWεH3W3Hdε−υuL1ℑS3H+υuL1ℑS3H3+ηX3V3L1∫0∞ℒ¯1εlnXεVεXVdε+ηX3V3L2∫0∞ℒ¯2εlnYεYdε+υℑS3H3uL1L3∫0∞ℒ¯3εlnSεHεSHdε+υℑS3H3uL1L4∫0∞ℒ¯4εlnWεWdε.
Using the equalities given by equations (36) and (37) in case of i,j=3, we get the following equation:(60)dϑ3dt=−αXX−X32−υγuL1SS−S32−1L1υY3S32−S3S−SS3−ηX3V3X3X−1−lnX3X−υℑS3H3uL1S3S−1−lnS3S−ηX3V3L1∫0∞ℒ¯1εXεVεY3X3V3Y−1−lnXεVεY3X3V3Ydε−ηX3V3L2a∫0∞ℒ¯2εYεV3Y3V−1−lnYεV3Y3Vdε−υℑS3H3uL1L3∫0∞ℒ¯3εSεHεW3S3H3W−1−lnSεHεW3S3H3Wdε−υℑS3H3uL1L4∫0∞ℒ¯4εWεH3W3H−1−lnWεH3W3Hdε.
Therefore, equation (60) becomes the following equation:(61)dϑ3dt=−αXX−X32+υuα℘+aγηL2aηuSL1L2S−S32aηℑuλρL1L2L3L4ℑkλL3L4+βυωuα℘+aγηL2−1−ηX3V3FX3X−υℑS3H3uL1FS3S−ηX3V3L1∫0∞ℒ¯1εFXεVεY3X3V3Ydε−ηX3V3L2∫0∞ℒ¯2εFYεV3Y3Vdε−υℑS3H3uL1L3∫0∞ℒ¯3εFSεHεW3S3H3Wdε−υℑS3H3uL1L4∫0∞ℒ¯4εFWεH3W3Hdε.
Since1<R3≤1+aηℑλξL2L3L4/βωuα℘+aγηL2, then dϑ3/dt≤0 for the positive values of X,Y,V,S,W,andH. Moreover, dϑ3/dt=0 when X=X3, S=S3, Yt=Y3, Vt=V3, Wt=W3, and Ht=H3. The model trajectories that converge to T3′ be the L.I.S of T3=X,Y,V,S,W,H:dϑ3/dt=0. Hence, T3′=EPVH and EPVH is G.A.S according to Lyapunov–LaSalle stability theorem.All equilibria of model (2) with the existence conditions and global stability constraints are summarized in Table 1.Table 1
Equilibrium points of model (2), existence conditions, and global stability conditions.
EquilibriumExistence conditionsGlobal stability conditionsEP0=X0,0,0,S0,0,0NoneR1≤1 and R2≤1EPH=X1,0,0,S1,W1,H1R1>1R1>1 and R4≤1EPV=X2,Y2,V2,S2,0,0R2>1R2>1 and R3≤1EPVH=X3,Y3,V3,S3,W3,H3R3>1 and R4>1R4>1 and 1<R3≤1+aηℑλξL2L3L4/βωuα℘+aγηL2


## 5. Numerical Simulations

We execute numerical simulations in this part to enhance the outcomes of Theorems1–4. Moreover, the impact of time delays on system dynamical behavior will be tested. To transform a model with distributed time delay (2) to a discrete one, we choose a Dirac delta function D. as a specific formula of kernel gi. as follows:(62)giν=Dν−εi,εi∈0,∞,i=1,2,3,4.Then, we get the following equation:(63)Lj=∫0∞Dς−ϵje−mjςdς=e−mjϵj,j=1,2,3,4.Thus, model (2) is reduced as follows:(64)X˙t=ρ−αXt−ηXtVt,Y˙t=ηe−m1ϵ1Xt−ϵ1Vt−ϵ1−kYt−υYtSt,V˙t=ae−m2ϵ2Yt−ϵ2−℘Vt,S˙t=ξ+uYtSt−γSt−ℑStHt,W˙t=ℑe−m3ϵ3St−ϵ3Ht−ϵ3−βWt,H˙t=λe−m4ϵ4Wt−ϵ4−ωHt.For model (64), the threshold parameters are given by the following equation:(65)R1=ξℑλe−m3ϵ3+m4ϵ4βγω,R2=aγηρe−m1ϵ1+m2ϵ2α℘γk+υξ,R3=λξβω+uλρe−m1ϵ1ℑkλe−m3ϵ3+m4ϵ4+βυωaηℑe−m2ϵ2+m3ϵ3+m4ϵ4uα℘+aγηe−m2ϵ2,R4=aηℑλρe−m1ϵ1+m2ϵ2+m3ϵ3+m4ϵ4α℘ℑkλe−m3ϵ3+m4ϵ4+βυω.To solve system (18) numerically, we change some parameters values whilst assigning fixed estimate to the rest parameters (Table 2). We modify the parameters η,υ,℘, and ℑ to test the conclusions of Theorems 1–4. Furthermore, to test the impact of the time delays upon COVID-19/AIDS dynamics, delays parameters ϵ1,ϵ2,ϵ3, and ϵ4 have been changed.Table 2
Model (64) parameter values.
ParametersValueReferenceρ0.02241[24]α10−3[24]ηVaried—k0.11[24]υVaried—a0.24[24]℘Varied—ξ10[19]u0.1[44]γ0.01[45]ℑVaried—β0.5[43]λ5[35]ω2[35]m10.1—m20.2—m30.2—m40.3—ϵi,i=1,…,4Varied—


### 5.1. Stability of Equilibrium Points

During this part, we choose delay parameters as follows:ϵ1=1, ϵ2=0.8, ϵ3=1, and ϵ4=0.8. Additionally, we select three distinct starting conditions of the model (18):Initial-1:Xϵ=5, Yϵ=0.0001, Vϵ=0.0002, Sϵ=100, Wϵ=5, and Hϵ=10,Initial-2:Xϵ=10, Yϵ=0.001, Vϵ=0.002, Sϵ=200, Wϵ=10, and Hϵ=15,Initial-3:Xϵ=15, Yϵ=0.002, Vϵ=0.003, Sϵ=300, Wϵ=15, and Hϵ=20.Here,ϵ∈−maxϵ1,ϵ2,ϵ3,ϵ4,0 and it is optional to pick these values. Moreover, the initial conditions are split into three groups to provide global stability for any starting conditions. To dissolve system (18), we utilize MATLAB solver dde23. Based on equilibrium points EP0, EPH, EPV, and EPVH global stability explained in Theorems 1–4, the simulations are divided into four cases. In these instances, we change values of η, υ, ℘, and ℑ of system (18). Other parameters values are set and recorded in Table 2. The four scenarios are detailed as follows:(i)
Case 1 (stability ofEP0): we take η=0.006, υ=0.01, ℘=0.3, and ℑ=0.0001. The thresholds in this case are given by R1=0.322<1 and R2=0.008<1. In harmony with Theorem 1, the equilibrium EP0=22.41,0,0,1000,0,0 is G.A.S (Figure 2). This is the best case scenario when the person is free of SARS-CoV-2 and HIV-1 infection.(ii)
Case 2 (stability ofEPH): we get η=0.0006, υ=0.01, ℘=0.3, and ℑ=0.0016. This provides us with R1=5.15>1 and R4=0.004<1. According to Theorem 2, the equilibrium EPH=22.41,0,0,186.478,13.3211,27.266 is G.A.S (Figure 3). This simulates the situation in which a person has HIV-1 infection with depressed CD4+ T cell levels, but SARS-CoV-2 infection is not present.(iii)
Case 3 (stability ofEPV): we select η=2.9, υ=0.002, ℘=0.1, and ℑ=0.0001. This gives R2=56.997>1 and R3=0.3535<1. In this situation, the system solutions converge globally to equilibrium EPV=0.4285,0.0087,0.018,1094.69,0,0. This result accords with Theorem 3 (Figure 4). This scenario simulates the case of a person infected with SARS-CoV-2 but not HIV-1 infection.(iv)
Case 4 (stability ofEPVH): we consider η=2.9, υ=0.02, ℘=0.1, and ℑ=0.0016. This implies that R3=5.19433>1, R3<1+aηℑλξL2L3L4/βωuα℘+aγηL2=6.1436, and R4=30.1279>1. In agreement with Theorem 4, the equilibrium EPVH=0.7155,0.005,0.0105,186.478,13.48,27.59 is G.A.S (Figure 5). In this case, COVID-19/AIDS coinfection occurs, where an HIV-1 patient gets infected with SARS-CoV-2. CD4+ T cells, which are the main target of HIV-1, are recruited to eliminate SARS-CoV-2 infection from the body. However, if the patient has low CD4+ T cell counts, the clearance of SARS-CoV-2 may not be achieved. This can cause severe infection and death.Figure 2
The numerical simulations of model (64) for η=0.006, υ=0.01, ℘=0.3, and ℑ=0.0001 using three different initial conditions sets. Uninfected equilibrium EP0=22.41,0,0,1000,0,0 is G.A.S. (a) Uninfected epithelial cells. (b) Infected epithelial cells. (c) SARS-CoV-2. (d) Uninfected CD4+ T cells. (e) Infected CD4+ T cells. (f) HIV-1.
(a)(b)(c)(d)(e)(f)Figure 3
The numerical simulations of model (64) for η=0.0006, υ=0.01, ℘=0.3, and ℑ=0.0016 using three different initial conditions sets. HIV-1 monoinfection equilibrium EPH=22.41,0,0,186.478,13.3211,27.266 is G.A.S. (a) Uninfected epithelial cells. (b) Infected epithelial cells. (c) SARS-CoV-2. (d) Uninfected CD4+ T cells. (e) Infected CD4+ T cells. (f) HIV-1.
(a)(b)(c)(d)(e)(f)Figure 4
The numerical simulations of model (64) for η=2.9, υ=0.002, ℘=0.1, and ℑ=0.0001 using three different initial conditions sets. SARS-CoV-2 monoinfection equilibrium EPV=0.4285,0.0087,0.018,1094.69,0,0 is G.A.S. (a) Uninfected epithelial cells. (b) Infected epithelial cell. (c) SARS-CoV-2. (d) Uninfected CD4+ T cells. (e) Infected CD4+ T cells. (f) HIV-1.
(a)(b)(c)(d)(e)(f)Figure 5
Model (64) numerical simulations for η=2.9, υ=0.02, ℘=0.1, and ℑ=0.0016 with three different initial conditions sets. COVID-19/AIDS coinfection equilibrium EPVH=0.7155,0.005,0.0105,186.478,13.48,27.59 is G.A.S. (a) Uninfected epithelial cells. (b) Infected epithelial cells. (c) SARS-CoV-2. (d) Uninfected CD4+ T cells. (e) Infected CD4+ T cells. (f) HIV-1.
(a)(b)(c)(d)(e)(f)


### 5.2. Impact of Time Delays on COVID-19/AIDS Dynamics

Here, we adjust parameters of delayϵi,i=1,2,…,4 and set the parameters values η=2.9,υ=0.02,℘=0.1, and ℑ=0.0016. Since R1,R2,R3, and R4 offered by equation (65) rely on ϵi,i=1,2,…,4, varying parameters ϵi will convert stability of the equilibria. We consider the following cases:(D.P.S1)ϵ1=ϵ2=ϵ3=ϵ4=0(D.P.S2)ϵ1=0.3,ϵ2=0.4,ϵ3=0.5, and ϵ4=0.6(D.P.S3)ϵ1=10,ϵ2=11,ϵ3=12, and ϵ4=13With the above values, we solve model (64) with given initial conditions:Initial-3: (Xϵ,Yϵ,Vϵ,Sϵ,Wϵ,Hϵ) = (15, 0.002, 0.003, 300, 15, 20).The inclusion of time delays can increase the number of uninfected epithelial and CD4+ T cells while diminish the number of other compartments, as shown in Figure6. Table 3 shows the values R1 and R2 for selected values of ϵi, i=1,2,…,4. Clearly, R1 and R2 decrease when ϵi are increased, and accordingly, the stability of EP0 can be changed. Let us compute the critical value of the time delay that changes the stability of EP0. Without loss of generality, we let the parameters ϵ3=ϵ4=ϵ34 and ϵ1=ϵ2=ϵ12, and write R1 and R2 as functions of ϵ34 and ϵ12, respectively, as follows:(66)R1ϵ34=ξℑλe−m3+m4ϵ34βγω,R2ϵ12=aγηρe−m1+m2ϵ12α℘γk+υξ.Figure 6
Model (64) numerical simulations for η=2.9, υ=0.02, ℘=0.1, and ℑ=0.0016 with three different sets of delay parameters. (a) Uninfected epithelial cells. (b) Infected epithelial cells. (c) SARS-CoV-2. (d) Uninfected CD4+ T cells. (e) Infected CD4+ T cells. (f) HIV-1.
(a)(b)(c)(d)(e)(f)Table 3
The variation ofR1 and R2 with respect to the delay parameters.
Delay parametersR1R2ϵ12=ϵ34=087.75602ϵ12=0.4 and ϵ34=0.27.246.88ϵ12=0.6 and ϵ34=0.46.556.48ϵ12=3 and ϵ34=14.853.15ϵ12=5 and ϵ34=31.791.731ϵ12=6.82823 and ϵ34=4.1588811ϵ12=7 and ϵ34=50.660.95ϵ12=8 and ϵ34=60.39830.704ϵ12=14 and ϵ34=120.01980.1163ϵ12=18 and ϵ34=160.00270.035ϵ12=24 and ϵ34=200.00040.0058To compel basic reproduction numbersR1 and R2 to verify R1ϵ34≤1 and R2ϵ12≤1, respectively, we choose the following equations:(67)ϵ34≥ϵ34min,whereϵ34min=max0,1m3+m4lnξℑλβγω.And(68)ϵ12≥ϵ12min,whereϵ12min=max0,1m1+m2lnaγηρα℘γk+υξ.Therefore, ifϵ34≥ϵ34min and ϵ12≥ϵ12min, then EP0 is G.A.S. Computing ϵ34 and ϵ12 gives ϵ34=4.15888 and ϵ12=6.82823, respectively. It follows(i)
Ifϵ34≥4.15888 and ϵ12≥6.82823, then R1ϵ34≤1, R2ϵ12≤1, and EP0 is G.A.S.(ii)
Ifϵ34<4.15888 or ϵ12<6.82823, then R1ϵ34>1, R2ϵ12>1, and EP0 will lose its stability.


## 5.1. Stability of Equilibrium Points

During this part, we choose delay parameters as follows:ϵ1=1, ϵ2=0.8, ϵ3=1, and ϵ4=0.8. Additionally, we select three distinct starting conditions of the model (18):Initial-1:Xϵ=5, Yϵ=0.0001, Vϵ=0.0002, Sϵ=100, Wϵ=5, and Hϵ=10,Initial-2:Xϵ=10, Yϵ=0.001, Vϵ=0.002, Sϵ=200, Wϵ=10, and Hϵ=15,Initial-3:Xϵ=15, Yϵ=0.002, Vϵ=0.003, Sϵ=300, Wϵ=15, and Hϵ=20.Here,ϵ∈−maxϵ1,ϵ2,ϵ3,ϵ4,0 and it is optional to pick these values. Moreover, the initial conditions are split into three groups to provide global stability for any starting conditions. To dissolve system (18), we utilize MATLAB solver dde23. Based on equilibrium points EP0, EPH, EPV, and EPVH global stability explained in Theorems 1–4, the simulations are divided into four cases. In these instances, we change values of η, υ, ℘, and ℑ of system (18). Other parameters values are set and recorded in Table 2. The four scenarios are detailed as follows:(i)
Case 1 (stability ofEP0): we take η=0.006, υ=0.01, ℘=0.3, and ℑ=0.0001. The thresholds in this case are given by R1=0.322<1 and R2=0.008<1. In harmony with Theorem 1, the equilibrium EP0=22.41,0,0,1000,0,0 is G.A.S (Figure 2). This is the best case scenario when the person is free of SARS-CoV-2 and HIV-1 infection.(ii)
Case 2 (stability ofEPH): we get η=0.0006, υ=0.01, ℘=0.3, and ℑ=0.0016. This provides us with R1=5.15>1 and R4=0.004<1. According to Theorem 2, the equilibrium EPH=22.41,0,0,186.478,13.3211,27.266 is G.A.S (Figure 3). This simulates the situation in which a person has HIV-1 infection with depressed CD4+ T cell levels, but SARS-CoV-2 infection is not present.(iii)
Case 3 (stability ofEPV): we select η=2.9, υ=0.002, ℘=0.1, and ℑ=0.0001. This gives R2=56.997>1 and R3=0.3535<1. In this situation, the system solutions converge globally to equilibrium EPV=0.4285,0.0087,0.018,1094.69,0,0. This result accords with Theorem 3 (Figure 4). This scenario simulates the case of a person infected with SARS-CoV-2 but not HIV-1 infection.(iv)
Case 4 (stability ofEPVH): we consider η=2.9, υ=0.02, ℘=0.1, and ℑ=0.0016. This implies that R3=5.19433>1, R3<1+aηℑλξL2L3L4/βωuα℘+aγηL2=6.1436, and R4=30.1279>1. In agreement with Theorem 4, the equilibrium EPVH=0.7155,0.005,0.0105,186.478,13.48,27.59 is G.A.S (Figure 5). In this case, COVID-19/AIDS coinfection occurs, where an HIV-1 patient gets infected with SARS-CoV-2. CD4+ T cells, which are the main target of HIV-1, are recruited to eliminate SARS-CoV-2 infection from the body. However, if the patient has low CD4+ T cell counts, the clearance of SARS-CoV-2 may not be achieved. This can cause severe infection and death.Figure 2
The numerical simulations of model (64) for η=0.006, υ=0.01, ℘=0.3, and ℑ=0.0001 using three different initial conditions sets. Uninfected equilibrium EP0=22.41,0,0,1000,0,0 is G.A.S. (a) Uninfected epithelial cells. (b) Infected epithelial cells. (c) SARS-CoV-2. (d) Uninfected CD4+ T cells. (e) Infected CD4+ T cells. (f) HIV-1.
(a)(b)(c)(d)(e)(f)Figure 3
The numerical simulations of model (64) for η=0.0006, υ=0.01, ℘=0.3, and ℑ=0.0016 using three different initial conditions sets. HIV-1 monoinfection equilibrium EPH=22.41,0,0,186.478,13.3211,27.266 is G.A.S. (a) Uninfected epithelial cells. (b) Infected epithelial cells. (c) SARS-CoV-2. (d) Uninfected CD4+ T cells. (e) Infected CD4+ T cells. (f) HIV-1.
(a)(b)(c)(d)(e)(f)Figure 4
The numerical simulations of model (64) for η=2.9, υ=0.002, ℘=0.1, and ℑ=0.0001 using three different initial conditions sets. SARS-CoV-2 monoinfection equilibrium EPV=0.4285,0.0087,0.018,1094.69,0,0 is G.A.S. (a) Uninfected epithelial cells. (b) Infected epithelial cell. (c) SARS-CoV-2. (d) Uninfected CD4+ T cells. (e) Infected CD4+ T cells. (f) HIV-1.
(a)(b)(c)(d)(e)(f)Figure 5
Model (64) numerical simulations for η=2.9, υ=0.02, ℘=0.1, and ℑ=0.0016 with three different initial conditions sets. COVID-19/AIDS coinfection equilibrium EPVH=0.7155,0.005,0.0105,186.478,13.48,27.59 is G.A.S. (a) Uninfected epithelial cells. (b) Infected epithelial cells. (c) SARS-CoV-2. (d) Uninfected CD4+ T cells. (e) Infected CD4+ T cells. (f) HIV-1.
(a)(b)(c)(d)(e)(f)


## 5.2. Impact of Time Delays on COVID-19/AIDS Dynamics

Here, we adjust parameters of delayϵi,i=1,2,…,4 and set the parameters values η=2.9,υ=0.02,℘=0.1, and ℑ=0.0016. Since R1,R2,R3, and R4 offered by equation (65) rely on ϵi,i=1,2,…,4, varying parameters ϵi will convert stability of the equilibria. We consider the following cases:(D.P.S1)ϵ1=ϵ2=ϵ3=ϵ4=0(D.P.S2)ϵ1=0.3,ϵ2=0.4,ϵ3=0.5, and ϵ4=0.6(D.P.S3)ϵ1=10,ϵ2=11,ϵ3=12, and ϵ4=13With the above values, we solve model (64) with given initial conditions:Initial-3: (Xϵ,Yϵ,Vϵ,Sϵ,Wϵ,Hϵ) = (15, 0.002, 0.003, 300, 15, 20).The inclusion of time delays can increase the number of uninfected epithelial and CD4+ T cells while diminish the number of other compartments, as shown in Figure6. Table 3 shows the values R1 and R2 for selected values of ϵi, i=1,2,…,4. Clearly, R1 and R2 decrease when ϵi are increased, and accordingly, the stability of EP0 can be changed. Let us compute the critical value of the time delay that changes the stability of EP0. Without loss of generality, we let the parameters ϵ3=ϵ4=ϵ34 and ϵ1=ϵ2=ϵ12, and write R1 and R2 as functions of ϵ34 and ϵ12, respectively, as follows:(66)R1ϵ34=ξℑλe−m3+m4ϵ34βγω,R2ϵ12=aγηρe−m1+m2ϵ12α℘γk+υξ.Figure 6
Model (64) numerical simulations for η=2.9, υ=0.02, ℘=0.1, and ℑ=0.0016 with three different sets of delay parameters. (a) Uninfected epithelial cells. (b) Infected epithelial cells. (c) SARS-CoV-2. (d) Uninfected CD4+ T cells. (e) Infected CD4+ T cells. (f) HIV-1.
(a)(b)(c)(d)(e)(f)Table 3
The variation ofR1 and R2 with respect to the delay parameters.
Delay parametersR1R2ϵ12=ϵ34=087.75602ϵ12=0.4 and ϵ34=0.27.246.88ϵ12=0.6 and ϵ34=0.46.556.48ϵ12=3 and ϵ34=14.853.15ϵ12=5 and ϵ34=31.791.731ϵ12=6.82823 and ϵ34=4.1588811ϵ12=7 and ϵ34=50.660.95ϵ12=8 and ϵ34=60.39830.704ϵ12=14 and ϵ34=120.01980.1163ϵ12=18 and ϵ34=160.00270.035ϵ12=24 and ϵ34=200.00040.0058To compel basic reproduction numbersR1 and R2 to verify R1ϵ34≤1 and R2ϵ12≤1, respectively, we choose the following equations:(67)ϵ34≥ϵ34min,whereϵ34min=max0,1m3+m4lnξℑλβγω.And(68)ϵ12≥ϵ12min,whereϵ12min=max0,1m1+m2lnaγηρα℘γk+υξ.Therefore, ifϵ34≥ϵ34min and ϵ12≥ϵ12min, then EP0 is G.A.S. Computing ϵ34 and ϵ12 gives ϵ34=4.15888 and ϵ12=6.82823, respectively. It follows(i)
Ifϵ34≥4.15888 and ϵ12≥6.82823, then R1ϵ34≤1, R2ϵ12≤1, and EP0 is G.A.S.(ii)
Ifϵ34<4.15888 or ϵ12<6.82823, then R1ϵ34>1, R2ϵ12>1, and EP0 will lose its stability.


## 6. Discussion

Coinfection between COVID-19/AIDS has become a serious problem during COVID-19 pandemic. Mathematical modeling represents a main tool in helping experimental studies understand new diseases. We studied a within-host COVID-19/AIDS coinfection model with distributed delays in this paper. The model explores the contacts between healthy epithelial cells, infected epithelial cells, free SARS-CoV-2 particles, uninfected CD4+ T cells, infected CD4+ T cells, and free HIV-1 particles. There are four equilibrium points for the model with the following listed properties:(a)
Uninfected equilibriumEP0: its existence is permanent and it is G.A.S if R1≤1 and R2≤1. This represents the situation of a person without SARS-CoV-2 or HIV-1 infections.(b)
The HIV-1 monoinfection equilibriumEPH exists if R1>1, and it is G.A.S if R4≤1. At this point, the person has only HIV-1 infection, but he is not infected by SARS-CoV-2.(c)
SARS-CoV-2 monoinfection equilibriumEPV is appeared when R2>1, and if R3≤1, then it is G.A.S. It is the instance of a person who is suffering from SARS-CoV-2 infection only.(d)
COVID-19/AIDS coinfection equilibriumEPVH exists and G.A.S if R4>1 and 1<R3≤1+aηℑλξL2L3L4/βωuα℘+aγηL2. In this case, the patient suffers from COVID-19/AIDS coinfection.The numerical and theoretical results were found to be in agreement. The time delays increase the concentrations of uninfected epithelial and CD4+ T cells, while they decrease concentrations of free SARS-CoV-2 and HIV-1 particles. Thus, parameters of delay can be examined and used in developing effective treatments for COVID-19/AIDS coinfected patients. Moreover, the model with distributed delays confirmed the effect observed in [20] that low numbers of CD4+ T cells can increase the risk of severe SARS-CoV-2 infection in coinfected patient. Thus, our model can be used to estimate the parameters required to get rid of SARS-CoV-2 in HIV-1 patients. Also, a bifurcation analysis can be executed in order to get a deeper understanding of the stability changes. Furthermore, the work can be developed by finding a better approximation of all parameters in model (2) through fitting with real data. We will keep these points in mind for future projects. [46].



---

*Source: 9129187-2022-12-22.xml*