Stability of HIV-1 Dynamics Models with Viral and Cellular Infections in the Presence of Macrophages

Abstract

In this research work, we suggest two mathematical models that take into account (i) two categories of target cells, CD$4^{+}$T cells and macrophages, and (ii) two modes of infection transmissions, the direct virus-to-cell (VTC) method and cell-to-cell (CTC) infection transmission, where CTC is an effective method of spreading human immunodeficiency virus type-1 (HIV-1), as with the VTC method. The second model incorporates four time delays. In both models, the presence of a bounded and positive solution of the biological model is investigated. The existence conditions of all equilibria are established. The basic reproduction number ${\mathcal{R}}_0$ that identifies a disease index is obtained. Lyapunov functions are utilized to verify the global stability of all equilibria. The theoretical findings are verified through numerical simulations. According to the outcomes, the trajectories of the solutions approach the infection-free equilibrium and infection-present equilibrium when ${\mathcal{R}}_0\le 1$ and ${\mathcal{R}}_0>1$, respectively. Further, we study the sensitivity analysis to investigate how the values of all the parameters of the suggested model affect ${\mathcal{R}}_0$ for given data. We discuss the impact of the time delay on HIV-1 progression. We find that a longer time delay results in suppression of the HIV-1 infection and vice versa.

1. Introduction

Human immunodeficiency virus type-1 (HIV-1) has been an area of attention in many studies for a long time; HIV-1 causes Acquired Immunodeficiency Syndrome (AIDS). World HIV-1 prevalence reached 38.4 million people by the end of 2021, with approximately 680,000 of people dying from HIV-1-related diseases and an estimated 1.5 million people acquiring HIV-1 [1]. HIV-1 is an RNA virus that targets the CD$4^{+}$T cells of the immune system. Despite efforts by scientists to create measures to battle this virus as well as research on its immunological and biological features and clinical outcomes, the study of HIV-1 dynamics and the identification of factors that control the infection process are still of interest.

Various mathematical models of numerous viral infections have confirmed their effectiveness in describing virus dynamics within hosts (see, e.g., [2,3,4,5,6]). Wide research on HIV infection has been carried out to understand virus dynamics between hosts [7,8,9,10]. Furthermore, there are other studies that have examined the dynamics of HIV-1, the progression of infection, and the dynamics of virus and immune system interaction within hosts [11,12,13,14]. Modeling transmissible diseases mathematically has become an important tool in analyzing and controlling infectious diseases. Models and simulations can help construct and test hypotheses, evaluate quantitative hypotheses, provide specific answers, establish the impact of parameter changes, and provide parameter estimations [15,16].

The traditional method of HIV-1 transmission is through virus-to-cell (VTC) contact, but recent studies indicate that HIV-1 can also be efficiently spread via cell-to-cell contact. It is reported that the CTC transmission mechanism can result in roughly $60\%$ of the total virus infection and that this infection mode shortens the period of virus creation by $0.9$ times [17]. Therefore, mathematical modeling should take into account the CTC viral infection method, which may be crucial to the virus’s ability to spread. Culshaw et al. [18] built a two-dimensional HIV-1 model with the CTC mode. Then, researchers started looking into HIV-1 mathematical models that account for both VTC and CTC transmission modes [19,20,21]. Although the mathematical analysis of models with CTC transmission is more difficult, the theoretical conclusions are more plentiful. These models can more accurately depict the biological process of HIV-1 infection.

The standard HIV-1 dynamics model assumes only that the virus interacts with one type of cell [22]; this model is based on ordinary differential equations and has been extended to incorporate the existing time delay that occurs between when a virus enters a cell and when new particles are created [12,23,24,25,26]. It was reported in [11] that the HIV-1 virus is capable of infecting both macrophages and CD4${}^{+}$T cells. An HIV-1 infection model with two target cells was proposed in [27]; Ref. [28] contains five variables describing the densities of uninfected CD4${}^{+}$T cells ($x_1$), infected CD4${}^{+}$T cells ($y_1$), uninfected macrophages ($x_2$), infected macrophages ($y_2$), and free HIV-1 particles (v), where $(x_1,y_1,x_2,y_2,v)=$$(x_1\left(t\right),y_1\left(t\right),x_2\left(t\right),y_2\left(t\right),v\left(t\right))$ and t is the time. The model is formulated as:

$$\begin{array}{cc}{\dot{x}}_1\hfill & =\underset{\mathrm{production}\mathrm{of}{\mathrm{CD}4}^{+}\mathrm{T}\mathrm{cells}}{\underbrace{{\alpha }_1}}-\underset{\mathrm{death}}{\underbrace{{\gamma }_1{x}_1}}-\underset{\mathrm{HIV}-1\mathrm{infectious}\mathrm{transmission}}{\underbrace{{\beta }_1{x}_{1}v}},\hfill \end{array}$$

(1)

$$\begin{array}{cc}{\dot{y}}_1\hfill & =\underset{\mathrm{HIV}-1\mathrm{infectious}\mathrm{transmission}}{\underbrace{{\beta }_1{x}_{1}v}}-\underset{\mathrm{death}}{\underbrace{{\vartheta }_1{y}_1}},\hfill \end{array}$$

(2)

$$\begin{array}{cc}{\dot{x}}_2\hfill & =\underset{\mathrm{production}\mathrm{of}\mathrm{macrophages}}{\underbrace{{\alpha }_2}}-\underset{\mathrm{death}}{\underbrace{{\gamma }_2{x}_2}}-\underset{\mathrm{HIV}-1\mathrm{infectious}\mathrm{transmission}}{\underbrace{{\beta }_2{x}_{2}v}},\hfill \end{array}$$

(3)

$$\begin{array}{cc}{\dot{y}}_2\hfill & =\underset{\mathrm{HIV}-1\mathrm{infectious}\mathrm{transmission}}{\underbrace{{\beta }_2{x}_{2}v}}-\underset{\mathrm{death}}{\underbrace{{\vartheta }_2{y}_2}},\hfill \end{array}$$

(4)

$$\begin{array}{cc}\dot{v}\hfill & =\underset{\mathrm{production}\mathrm{of}\mathrm{HIV}-1\mathrm{from}\mathrm{infected}{\mathrm{CD}4}^{+}\mathrm{T}\mathrm{cells}\mathrm{and}\mathrm{macrophages}}{\underbrace{{\lambda }_1{\vartheta }_1{y}_1+{\lambda }_2{\vartheta }_2{y}_2}}-\underset{\mathrm{death}}{\underbrace{\phi v}}.\hfill \end{array}$$

(5)

Elaiw [29] provided an analysis study of the mathematical model (1)–(5). The model has been extended in different directions by including an eclipse phase [27,29], immune response [30,31], time delay [30,32], and optimal drug schedules [33,34,35].

In the literature, all HIV-1 studies with two target cells have considered only VTC infection as the mode of infection. Therefore, the purpose of this research work is to develop and analyze two dynamical HIV-1 models with two categories of target cells and two modes of infection: VTC and CTC. The second model is an extension of the first one by considering the effects of four types of time delays. We aim to determine the conditions in which the equilibrium points reach the state of stability, the effect of including the delay time on the behavior of solutions of the introduced model, and also the extent to which the parameters of the model influence the basic reproduction number ${\mathcal{R}}_0$ for given data.

The paper is structured as follows: The introduction is the first section. In Section 2, we introduce the proposed HIV-1 model with two types of target cells and two modes of infection. In Section 3, we incorporate into the proposed first model the discrete time delays. In both Section 2 and Section 3 and for the proposed models; the basic properties are established—such as the solutions’ nonnegativity and boundedness—and calculations are performed to determine the parameter ${\mathcal{R}}_0$ along with the probable equilibrium points and the conditions under which they exist, and we examine the global stability of the equilibria by using appropriate Lyapunov functions and applying the LaSalle invariance principle. In Section 4, some numerical simulations are carried out to ensure the theoretical results. In addition to analyzing sensitivity, we examine the stability of equilibria and test key parameters and time delays to determine how they affect the dynamics of the model. Finally, conclusions are provided in Section 5.

2. The Model

In this section, the authors reformulate and analyze the HIV-1 dynamics model to depict the dispersion of the HIV-1 virus in the host through two ways of transmission of the infection (VTC and CTC) and in the presence of two types of target cells (CD$4^{+}$T cells and macrophages).

2.1. Model Formulation

The proposed HIV-1 infection model incorporates both VTC and CTC transmission and two categories of target cells: CD$4^{+}$T cells and macrophages. The schematic diagram for the interaction is presented in Figure 1, and the model is as follows:

$$\begin{array}{cc}\hfill {\dot{x}}_{\ell }& ={\alpha }_{\ell }-{\gamma }_{\ell }{x}_{\ell }-{\beta }_{\ell }{x}_{\ell }v-{\overline{\beta }}_{\ell }{x}_{\ell }{y}_{\ell },\ell =1,2,\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill {\dot{y}}_{\ell }& ={\beta }_{\ell }{x}_{\ell }v+{\overline{\beta }}_{\ell }{x}_{\ell }{y}_{\ell }-{\vartheta }_{\ell }{y}_{\ell },\ell =1,2,\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \dot{v}& =\underset{\ell =1}{\sum ^2}{\lambda }_{\ell }{\vartheta }_{\ell }{y}_{\ell }-\phi v,\hfill \end{array}$$

(8)

where ${\overline{\beta }}_1{x}_1{y}_1$ and ${\overline{\beta }}_2{x}_2{y}_2$ are the CTC infection rates of the CD$4^{+}$T cells and macrophages, respectively. When the CTC infection mode is neglected, then Model (6)–(8) leads to Model (1)–(5).

2.2. Invariant Region

Through the following lemma, we prove that Model (6)–(8) is biologically well-defined.

 Lemma 1.

Let ${\zeta }_{\ell }>0,\ell =1,2,3$; then, there exists a positively invariant compact set for Model (6)–(8):

$$\Omega =\left\{\left(x_1,y_1,x_2,y_2,v\right)\in {\mathbb{R}}_{\ge 0}^5:0\le x_1,y_1\le {\zeta }_1,0\le x_2,y_2\le {\zeta }_2,0\le v\le {\zeta }_3\right\}.$$

 Proof.

We have

$$\begin{array}{cc}\hfill {\left.{\dot{x}}_{\ell }\right|}_{\left(x_{\ell }=0\right)}& ={\alpha }_{\ell }>0,\ell =1,2,\hfill \\ \hfill {\left.{\dot{y}}_{\ell }\right|}_{(y_{\ell }=0)}& ={\beta }_{\ell }{x}_{\ell }v\ge 0,\mathrm{for}\mathrm{all}{x}_{\ell },v\ge 0,\ell =1,2,\hfill \\ \hfill {\left.\dot{v}\right|}_{(v=0)}& =\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{y}_{\ell }\ge 0,\mathrm{for}\mathrm{all}{y}_{\ell }\ge 0,\ell =1,2;\hfill \end{array}$$

thus, the positively invariant property of ${\mathbb{R}}_{\ge 0}^5$ with respect to system (6)–(8) has been proofed.

Next, let $F_{\ell }=x_{\ell }+y_{\ell }.$ Then

$${\dot{F}}_{\ell }={\dot{x}}_{\ell }+{\dot{y}}_{\ell }={\alpha }_{\ell }-{\gamma }_{\ell }{x}_{\ell }-{\vartheta }_{\ell }{y}_{\ell }\le {\alpha }_{\ell }-{\delta }_{\ell }{F}_{\ell },$$

where, ${\delta }_{\ell }=min\{{\gamma }_{\ell },{\vartheta }_{\ell }\}$, $\ell =1,2$. Then

$$F_{\ell }\left(t\right)\le e^{-{\delta }_{\ell }t}\left(F_{\ell }\left(0\right)-\frac{{\alpha }_{\ell }}{{\delta }_{\ell }}\right)+\frac{{\alpha }_{\ell }}{{\delta }_{\ell }}.$$

This shows that $0\le F_{\ell }\left(t\right)\le {\zeta }_{\ell }$ for all $t\ge 0$ if $F_{\ell }\left(0\right)\le {\zeta }_{\ell },$ where ${\zeta }_{\ell }=\frac{{\alpha }_{\ell }}{{\delta }_{\ell }}$. Consequently, $0\le x_{\ell }\left(t\right),y_{\ell }\left(t\right)\le {\zeta }_{\ell }$ for all $t\ge 0$ if $x_{\ell }\left(0\right)+y_{\ell }\left(0\right)\le {\zeta }_{\ell },\ell =1,2$. Further,

$$\dot{v}\left(t\right)=\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{y}_{\ell }\left(t\right)-\phi v\left(t\right)\le \stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{\zeta }_{\ell }-\phi v\left(t\right);$$

thus, $0\le v\left(t\right)\le$${\zeta }_3$ for all $t\ge 0$ if $v\left(0\right)\le$ ${\zeta }_3$, where ${\zeta }_3={\displaystyle \sum _{\ell =1}^2}{\displaystyle \frac{{\lambda }_{\ell }{\vartheta }_{\ell }{\zeta }_{\ell }}{\phi }}.$ This guarantees the boundedness of $x_{\ell }\left(t\right),y_{\ell }\left(t\right)$ and $v\left(t\right),\ell =1,2$. □

2.3. Equilibria

To calculate the equilibrium points of System (6)–(8), we solve the following system:

$$\begin{array}{cc}0\hfill & ={\alpha }_{\ell }-{\gamma }_{\ell }{x}_{\ell }-{\beta }_{\ell }{x}_{\ell }v-{\overline{\beta }}_{\ell }{x}_{\ell }{y}_{\ell },\ell =1,2,\hfill \end{array}$$

(9)

$$\begin{array}{cc}0\hfill & ={\beta }_{\ell }{x}_{\ell }v+{\overline{\beta }}_{\ell }{x}_{\ell }{y}_{\ell }-{\vartheta }_{\ell }{y}_{\ell },\ell =1,2,\hfill \end{array}$$

(10)

$$\begin{array}{cc}0\hfill & =\underset{\ell =1}{\sum ^2}{\lambda }_{\ell }{\vartheta }_{\ell }{y}_{\ell }-\phi v;\hfill \end{array}$$

(11)

then, we obtain that System (6)–(8) has two equilibrium points:

(i)

   The infection-free equilibrium (IFE) ${\Pi }_0=(x_1^0,0,x_2^0,0,0),$ where $x_{\ell }^0=\frac{{\alpha }_{\ell }}{{\gamma }_{\ell }},\ell =1,2,$

(ii)

  The infection-present equilibrium (IPE) ${\Pi }^{*}=(x_1^{*},y_1^{*},x_2^{*},y_2^{*},v^{*}),$ where

$$x_{\ell }^{*}=\frac{{\alpha }_{\ell }}{{\gamma }_{\ell }+{\beta }_{\ell }{v}^{*}+{\overline{\beta }}_{\ell }{y}_{\ell }^{*}},y_{\ell }^{*}=\frac{-B_{\ell }+\sqrt{B_{\ell }^2+4A_{\ell }{C}_{\ell }}}{2A_{\ell }},$$

(12)

and

$$A_{\ell }={\vartheta }_{\ell }{\overline{\beta }}_{\ell },B_{\ell }=v^{*}{\vartheta }_{\ell }{\beta }_{\ell }+{\vartheta }_{\ell }{\gamma }_{\ell }-{\overline{\beta }}_{\ell }{\alpha }_{\ell },C_{\ell }=v^{*}{\beta }_{\ell }{\alpha }_{\ell }$$

(13)

for all $\ell =1,2$, in which $v^{*}$ satisfies this equation:

$$\phi v^{*}=\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{y}_{\ell }^{*}.$$

(14)

The basic reproduction number ${\mathcal{R}}_0$: In infectious disease studies, the basic reproduction number, abbreviated as ${\mathcal{R}}_0$, is a measure used to describe the contagiousness of transmissibility of an infectious disease. ${\mathcal{R}}_0$ signifies the average number of cells that are going to be infected by a cell that is already infected and is introduced into a susceptible cell. Utilizing the next-generation-matrix method approach described in [36], ${\mathcal{R}}_0$ of System (6)–(8) can be determined as follows:

The Jacobian of the matrix of new infection terms at IFE is given by:

$$\mathcal{F}=\left[\begin{array}{ccc}{\overline{\beta }}_1{x}_1^0& 0& {\beta }_1{x}_1^0\\ 0& {\overline{\beta }}_2{x}_2^0& {\beta }_2{x}_2^0\\ 0& 0& 0\end{array}\right],$$

and the Jacobian of the matrix of the other terms at IFE is as follows:

$$\nabla =\left[\begin{array}{ccc}{\vartheta }_1& 0& 0\\ 0& {\vartheta }_2& 0\\ -{\lambda }_1{\vartheta }_1& -{\lambda }_2{\vartheta }_2& \phi \end{array}\right].$$

The basic reproduction number is calculated as:

$${\mathcal{R}}_0=\rho \left(\mathcal{F}{\nabla }^{-1}\right)$$

where the matrix $\mathcal{F}{\nabla }^{-1}$ is obtained as:

$$\mathcal{F}{\nabla }^{-1}=\left[\begin{array}{ccc}\frac{{\overline{\beta }}_1{x}_1^0}{{\vartheta }_1}+\frac{{\beta }_1{x}_1^0{\lambda }_1}{\phi }& \frac{{\beta }_1{x}_1^0{\lambda }_2}{\phi }& \frac{{\beta }_1{x}_1^0}{\phi }\\ \frac{{\beta }_2{x}_2^0{\lambda }_1}{\phi }& \frac{{\overline{\beta }}_2{x}_2^0}{{\vartheta }_2}+\frac{{\beta }_2{x}_2^0{\lambda }_2}{\phi }& \frac{{\beta }_2{x}_2^0}{\phi }\\ 0& 0& 0\end{array}\right],$$

and $\rho \left(\mathcal{F}{\nabla }^{-1}\right)$ refers to the spectral radius of $\mathcal{F}{\nabla }^{-1}$. Hence,

$${\mathcal{R}}_0=\frac{1}{2}\left(\overline{\Phi }+\overline{\Psi }+\sqrt{{\left(\overline{\Phi }-\overline{\Psi }\right)}^2+4{\overline{\varphi }}_b{\overline{\psi }}_b}\right),$$

(15)

where

$$\begin{array}{cc}\hfill \overline{\Phi }& ={\overline{\varphi }}_a+{\overline{\varphi }}_b,\overline{\Psi }={\overline{\psi }}_a+{\overline{\psi }}_b,\hfill \\ \hfill {\overline{\varphi }}_a& =\frac{{\overline{\beta }}_1{x}_1^0}{{\vartheta }_1},{\overline{\varphi }}_b=\frac{{\beta }_1{x}_1^0{\lambda }_1}{\phi },{\overline{\psi }}_a=\frac{{\overline{\beta }}_2{x}_2^0}{{\vartheta }_2},{\overline{\psi }}_b=\frac{{\beta }_2{x}_2^0{\lambda }_2}{\phi }.\hfill \end{array}$$

(16)

In the following Lemmas, we show the condition when $v^{*}$ is positive.

 Lemma 2.

Suppose that ${\mathcal{R}}_0>1$. If $\overline{\Phi }<1,\overline{\Psi }<1,$ then $\overline{M}=\frac{{\overline{\varphi }}_b}{1-{\overline{\varphi }}_a}+\frac{{\overline{\psi }}_b}{1-{\overline{\psi }}_a}>1$.

 Proof.

Let ${\mathcal{R}}_0>1,$ then

$$\frac{1}{2}\left(\overline{\Phi }+\overline{\Psi }+\sqrt{{\left(\overline{\Phi }-\overline{\Psi }\right)}^2+4{\overline{\varphi }}_b{\overline{\psi }}_b}\right)>1$$

$$\begin{array}{cc}\hfill & ⟹\sqrt{{\left(\overline{\Phi }-\overline{\Psi }\right)}^2+4{\overline{\varphi }}_b{\overline{\psi }}_b}>2-\left(\overline{\Phi }+\overline{\Psi }\right)>0\hfill \\ \hfill & ⟹{\left(\overline{\Phi }-\overline{\Psi }\right)}^2+4{\overline{\varphi }}_b{\overline{\psi }}_b>4-4\left(\overline{\Phi }+\overline{\Psi }\right)+{\left(\overline{\Phi }+\overline{\Psi }\right)}^2>0\hfill \\ \hfill & ⟹{\overline{\varphi }}_b+{\overline{\psi }}_b-{\overline{\varphi }}_b{\overline{\psi }}_a-{\overline{\varphi }}_a{\overline{\psi }}_b>\left(1-{\overline{\varphi }}_a\right)\left(1-{\overline{\psi }}_a\right).\hfill \end{array}$$

(17)

Since $\overline{\Phi }<1,\overline{\Psi }<1,$ then $1-{\overline{\varphi }}_a>0$ and $1-{\overline{\psi }}_a>0.$ Therefore, (17) can be written as follows:

$$\frac{{\overline{\varphi }}_b+{\overline{\psi }}_b-{\overline{\varphi }}_b{\overline{\psi }}_a-{\overline{\varphi }}_a{\overline{\psi }}_b}{\left(1-{\overline{\varphi }}_a\right)\left(1-{\overline{\psi }}_a\right)}=\frac{{\overline{\varphi }}_b}{\left(1-{\overline{\varphi }}_a\right)}+\frac{{\overline{\psi }}_b}{\left(1-{\overline{\psi }}_a\right)}=\overline{M}>1.$$

□

 Lemma 3.

Suppose that ${\mathcal{R}}_0>1$, then the IFE ${\Pi }^{*}$ exists.

 Proof.

First, we have that any equilibrium exists that satisfies Equations (9)–(11). In the case of the equilibrium ${\Pi }^{*}$, we have $v\ne 0;$ then, from Equation (11), we find

$$\stackrel{2}{\sum _{\ell =1}}\frac{{\lambda }_{\ell }{\vartheta }_{\ell }{y}_{\ell }}{\phi }-v=0.$$

(18)

Adding Equations (9) and (10), we get

$${\vartheta }_{\ell }{y}_{\ell }={\alpha }_{\ell }-{\gamma }_{\ell }{x}_{\ell },\ell =1,2.$$

(19)

Substituting from Equation (19) into Equation (18), we obtain

$$\stackrel{2}{\sum _{\ell =1}}\frac{{\lambda }_{\ell }{\alpha }_{\ell }}{\phi }-\stackrel{2}{\sum _{\ell =1}}\frac{{\lambda }_{\ell }{\gamma }_{\ell }{x}_{\ell }}{\phi }-v=0.$$

Since $x_{\ell }=x_{\ell }\left(v\right)$, $y_{\ell }=y_{\ell }\left(v\right),\ell =1,2,$ then we can define a function $\mathcal{G}\left(v\right)$ as:

$$\mathcal{G}\left(v\right)=\stackrel{2}{\sum _{\ell =1}}\frac{{\lambda }_{\ell }{\alpha }_{\ell }}{\phi }-\stackrel{2}{\sum _{\ell =1}}\frac{{\lambda }_{\ell }{\gamma }_{\ell }{x}_{\ell }\left(v\right)}{\phi }-v,$$

(20)

in which $x_{\ell },y_{\ell }$ satisfy Equations (9)–(11) for $\ell =1,2$.

Now, we need to show that ∃$v^{*}>0$ such that $\mathcal{G}\left(v^{*}\right)=0$ as follows:

If $v=\check{v}={\displaystyle \sum _{\ell =1}^2}\frac{{\lambda }_{\ell }{\alpha }_{\ell }}{\phi }>0,$ then $x_{\ell }\left(\check{v}\right)>0,y_{\ell }\left(\check{v}\right)>0,\ell =1,2$ and $\mathcal{G}\left(\check{v}\right)=-{\displaystyle \sum _{\ell =1}^2}\frac{{\lambda }_{\ell }{\gamma }_{\ell }{x}_{\ell }\left(\check{v}\right)}{\phi }<0.$ Next, by calculating $\mathcal{G}\left(0\right)$ and ${\mathcal{G}}^{\prime }\left(0\right)$, we get:

$$\begin{array}{cc}\hfill \mathcal{G}\left(0\right)& =\stackrel{2}{\sum _{\ell =1}}\frac{{\lambda }_{\ell }{\alpha }_{\ell }}{\phi }-\stackrel{2}{\sum _{\ell =1}}\frac{{\lambda }_{\ell }{\gamma }_{\ell }{x}_{\ell }\left(0\right)}{\phi }\hfill \\ \hfill & =\frac{{\lambda }_1{\gamma }_1{x}_1^0}{\phi }\left(1-\frac{2}{2-\left(1-{\overline{\varphi }}_a\right)+\sqrt{{\left(1-{\overline{\varphi }}_a\right)}^2}}\right)+\frac{{\lambda }_2{\gamma }_2{x}_2^0}{\phi }\left(1-\frac{2}{2-\left(1-{\overline{\psi }}_a\right)+\sqrt{{\left(1-{\overline{\psi }}_a\right)}^2}}\right),\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\mathcal{G}}^{\prime }\left(0\right)& =-\stackrel{2}{\sum _{\ell =1}}\frac{{\lambda }_{\ell }{\gamma }_{\ell }}{\phi }{x}_{\ell }^{\prime }\left(0\right)-1\hfill \\ \hfill & =\frac{{\lambda }_1{x}_1^0{\beta }_1}{2\phi }\frac{\sqrt{{\left(1-{\overline{\varphi }}_a\right)}^2}+\left(1+{\overline{\varphi }}_a\right)}{{\left(1+\frac{\sqrt{{\left(1-{\overline{\varphi }}_a\right)}^2}-\left(1-{\overline{\varphi }}_a\right)}{2}\right)}^2\sqrt{{\left(1-{\overline{\varphi }}_a\right)}^2}}+\frac{{\lambda }_2{x}_2^0{\beta }_2}{2\phi }\frac{\sqrt{{\left(1-{\overline{\psi }}_a\right)}^2}+\left(1+{\overline{\psi }}_a\right)}{{\left(1+\frac{\sqrt{{\left(1-{\overline{\psi }}_a\right)}^2}-\left(1-{\overline{\psi }}_a\right)}{2}\right)}^2\sqrt{{\left(1-{\overline{\psi }}_a\right)}^2}}-1.\hfill \end{array}$$

In the next step, we calculate all possible cases of the two functions $\mathcal{G}\left(0\right),$${\mathcal{G}}^{\prime }\left(0\right),$ and the results are provided by

Table 1. Functions $\mathcal{G}\left(0\right)$ and ${\mathcal{G}}^{\prime }\left(0\right)$ and their corresponding values for different conditions.

Case	Conditions	$\mathcal{G}\left(0\right)$	${\mathcal{G}}^{\prime }\left(0\right)$
1	${\overline{\varphi }}_a=1,{\overline{\psi }}_a\le 1$	$0$	$+\infty$
2	${\overline{\varphi }}_a\le 1,{\overline{\psi }}_a=1$	$0$	$+\infty$
3	${\overline{\varphi }}_a<1,{\overline{\psi }}_a<1$	$0$	$\overline{M}-1>0$(from Lemma 2)
4	${\overline{\varphi }}_a\le 1,{\overline{\psi }}_a>1$	${\displaystyle \frac{{\lambda }_2{\gamma }_2{x}_2^0}{\phi }}\left({\displaystyle \frac{{\overline{\psi }}_a-1}{{\overline{\psi }}_a}}\right)>0$	−
5	${\overline{\varphi }}_a>1,{\overline{\psi }}_a\le 1$	${\displaystyle \frac{{\lambda }_1{\gamma }_1{x}_1^0}{\phi }}\left({\displaystyle \frac{{\overline{\varphi }}_a-1}{{\overline{\varphi }}_a}}\right)>0$	−
6	${\overline{\varphi }}_a>1,{\overline{\psi }}_a>1$	${\displaystyle \frac{{\lambda }_1{\gamma }_1{x}_1^0}{\phi }}\left({\displaystyle \frac{{\overline{\varphi }}_a-1}{{\overline{\varphi }}_a}}\right)+{\displaystyle \frac{{\lambda }_2{\gamma }_2{x}_2^0}{\phi }}\left({\displaystyle \frac{{\overline{\psi }}_{ac}-1}{{\overline{\psi }}_a}}\right)>0$	−

. As shown in Table 1, the function $\mathcal{G}\left(v\right)$ in Cases (1–3) is strictly increasing at $v=0,$ and $\mathcal{G}\left(0\right)>0$ in Cases (4–6), while $\mathcal{G}\left(\check{v}\right)$ has a negative value. This means that in all possible cases ∃$v^{*}\in (0,\check{v})$ satisfying $\mathcal{G}\left(v^{*}\right)=0$ if the condition ${\mathcal{R}}_0>1.$ Therefore, from Equations (12)–(13), we have $x_{\ell }^{*}>0,y_{\ell }^{*}>0,v^{*}>0,\ell =1,2$. Thus, the endemic equilibrium ${\Pi }^{*}$ exists when ${\mathcal{R}}_0>1.$ □

We conclude that for System (6)–(8):

(i)

   If ${\mathcal{R}}_0\le 1,$ then there will be only one equilibrium: ${\Pi }_0$;

(ii)

  If ${\mathcal{R}}_0>1,$ then there will be two equilibria: ${\Pi }_0$ and ${\Pi }^{*}$.

2.4. Global Properties

The stability of a system is an important characteristic of its qualitative behavior. Several stability theories are used in the study of dynamic systems. Among the most important stability concepts is Lyapunov stability. Using this concept, we can maintain the future states of a system arbitrarily close to the equilibrium by simply taking the initial condition sufficiently close to the equilibrium.

A major role is played by Lyapunov functions when studying the stability of dynamical systems. Many researchers have applied Lyapunov’s method to analyze the stability of nonlinear systems over the last century [37], and it has recently been widely used in many studies since this method has many advantages [38,39], including (i) Lyapunov theory is a useful tool when studying uncertain (especially nonlinear) systems with time-varying parameters, (ii) the theory provides procedures that are both efficient and insightful, (iii) in the case of important classes of problems and specific classes of functions, the theory is underpinned by suitable numerical methods such as those based on linear matrix inequalities (LMIs), and (iv) this method does not need to find the actual solution.

In the following, and based on Lyapunov’s method and LaSalle’s invariance principle (LIP) [40,41,42], we prove the global stability of System (6)–(8) at the infection-free equilibrium (IFE) ${\Pi }_0$ and at the disease infection point (IPE) ${\Pi }^{*}$. Let ${\Gamma }_{\ell }^{{}^{\prime }}$ be the largest invariant subset of ${\Gamma }_{\ell }=\left\{(x_1,x_{2,}{y}_1,y_2,v):\frac{d{\Xi }_{\ell }}{dt}=0\right\}$, where $\ell =1,2.$ Further, we define the function $\Theta \left(\xi \right)=\xi -1-ln\xi$. For the purpose of proving global stability, we utilize the following Lemma:

 Lemma 4.

Suppose that ${\mathcal{R}}_0\le 1;$ then,

 (i)

   ${\overline{\varphi }}_a\le 1,{\overline{\varphi }}_b\le 1,$${\overline{\psi }}_a\le 1,$ and ${\overline{\psi }}_b\le 1$

 (ii)

  If $M=\overline{\Phi }+\overline{\Psi }-\overline{\Phi }\overline{\Psi }+{\overline{\varphi }}_b{\overline{\psi }}_b;$ thus, $0<M\le 1$.

 Proof.

 (i)  Let ${\mathcal{R}}_0\le 1,$ then

$$\begin{array}{cc}\hfill 1& \ge \frac{1}{2}\left(\overline{\Phi }+\overline{\Psi }+\sqrt{{\left(\overline{\Phi }-\overline{\Psi }\right)}^2+4{\overline{\varphi }}_b{\overline{\psi }}_b}\right)\hfill \\ \hfill & \ge \frac{1}{2}\left(\overline{\Phi }+\overline{\Psi }+\sqrt{{\left(\overline{\Phi }-\overline{\Psi }\right)}^2}\right).\hfill \end{array}$$

• If $\overline{\Phi }\ge \overline{\Psi }$, then

  $$\begin{array}{cc}\hfill 1& \ge \frac{1}{2}\left(\overline{\Phi }+\overline{\Psi }+\left(\overline{\Phi }-\overline{\Psi }\right)\right)\hfill \\ \hfill & =\overline{\Phi };\hfill \end{array}$$

  that is ${\overline{\varphi }}_a\le 1,$ ${\overline{\varphi }}_b\le 1,$ and then $\overline{\Psi }\le 1.$
•  Similarly, if $\overline{\Psi }\ge \overline{\Phi }$, then

  $$\begin{array}{cc}\hfill 1& \ge \frac{1}{2}\left(\overline{\Phi }+\overline{\Psi }+\left(\overline{\Psi }-\overline{\Phi }\right)\right)\hfill \\ \hfill & =\overline{\Psi };\hfill \end{array}$$

  thus, we have ${\overline{\psi }}_a\le 1$,${\overline{\psi }}_b\le 1,$ and then, $\overline{\Phi }\le 1.$

(ii)

  Since $\overline{\Phi }\le 1$ and $\overline{\Psi }\le 1,$ then $\overline{\Phi }+\overline{\Psi }\le 2.$ Additionally, we have ${\mathcal{R}}_0\le 1$ that is

$$\begin{array}{cc}\hfill 0& <\frac{1}{2}\left(\overline{\Phi }+\overline{\Psi }+\sqrt{{\left(\overline{\Phi }-\overline{\Psi }\right)}^2+4{\overline{\varphi }}_b{\overline{\psi }}_b}\right)\le 1\hfill \\ \hfill & ⟹4\left(\overline{\Phi }+\overline{\Psi }\right)-{\left(\overline{\Phi }+\overline{\Psi }\right)}^2<4\left(\overline{\Phi }+\overline{\Psi }\right)-{\left(\overline{\Phi }+\overline{\Psi }\right)}^2+{\left(\overline{\Phi }-\overline{\Psi }\right)}^2+4{\overline{\varphi }}_b{\overline{\psi }}_b\le 4\hfill \\ \hfill & ⟹\left(\overline{\Phi }+\overline{\Psi }\right)\left[1-\frac{1}{4}\left(\overline{\Phi }+\overline{\Psi }\right)\right]<\left(\overline{\Phi }+\overline{\Psi }\right)-\overline{\Phi }\overline{\Psi }+{\overline{\varphi }}_b{\overline{\psi }}_b\le 1;\hfill \end{array}$$

here, $\overline{\Phi }+\overline{\Psi }>\overline{\Phi }\overline{\Psi }$ since $\overline{\Phi }\le 1,\overline{\Psi }\le 1.$ Consequently,

$$\begin{array}{cc}\hfill & ⟹0<\overline{\Phi }+\overline{\Psi }-\overline{\Phi }\overline{\Psi }+{\overline{\varphi }}_b{\overline{\psi }}_b\le 1\hfill \\ \hfill & ⟹0<M\le 1.\hfill \end{array}$$

This completes the proof. □

 Theorem 1.

For System (6)–(8), if the value of ${\mathcal{R}}_0$ is less than or equal to one $({\mathcal{R}}_0\le 1)$, then ${\Pi }_0$ is globally asymptotically stable (GAS).

 Proof.

Let ${\mathcal{R}}_0\le 1$, and construct a function ${\Xi }_0(x_1,y_1,x_2,y_2,v)$ as:

$${\Xi }_0=\stackrel{2}{\sum _{\ell =1}}{\eta }_{\ell }\left[x_{\ell }^0\Theta \left(\frac{x_{\ell }}{x_{\ell }^0}\right)+y_{\ell }\right]+{\eta }_{3}v,$$

where

$$\begin{array}{cc}\hfill {\eta }_1& ={\vartheta }_1{\vartheta }_2{\lambda }_1(1-{\overline{\psi }}_a),\hfill \\ \hfill {\eta }_2& ={\vartheta }_1{\vartheta }_2{\lambda }_2(1-{\overline{\varphi }}_a),\hfill \\ \hfill {\eta }_3& ={\vartheta }_1{\vartheta }_2(1-{\overline{\varphi }}_a)(1-{\overline{\psi }}_a).\hfill \end{array}$$

Clearly, ${\Xi }_0(x_1,y_1,x_2,y_2,v)>0$ for all $x_1,y_1,x_2,y_2,v>0$, and ${\Xi }_0(x_1^0,0,x_2^0,0,0)=0$. Calculating ${\displaystyle \frac{d{\Xi }_0}{dt}}$ along System (6)–(8), we obtain

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\Xi }_0}{dt}}& =\stackrel{2}{\sum _{\ell =1}}{\eta }_{\ell }\left(1-\frac{x_{\ell }^0}{x_{\ell }}\right){\dot{x}}_{\ell }+\stackrel{2}{\sum _{\ell =1}}{\eta }_{\ell }{\dot{y}}_{\ell }+{\eta }_3\dot{v}\hfill \\ \hfill & =\stackrel{2}{\sum _{\ell =1}}{\eta }_{\ell }\left(1-\frac{x_{\ell }^0}{x_{\ell }}\right)\left({\alpha }_{\ell }-{\gamma }_{\ell }{x}_{\ell }-{\beta }_{\ell }{x}_{\ell }v-{\overline{\beta }}_{\ell }{x}_{\ell }{y}_{\ell }\right)+\stackrel{2}{\sum _{\ell =1}}{\eta }_{\ell }\left({\beta }_{\ell }{x}_{\ell }v+{{\overline{\beta }}_{\ell }x}_{\ell }{y}_{\ell }-{\vartheta }_{\ell }{y}_{\ell }\right)\hfill \\ \hfill & +{\eta }_3\left(\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{y}_{\ell }-\phi v\right)\hfill \\ \hfill & =-\stackrel{2}{\sum _{\ell =1}}{\eta }_{\ell }{\gamma }_{\ell }\frac{{\left(x_{\ell }-x_{\ell }^0\right)}^2}{x_{\ell }}+\stackrel{2}{\sum _{\ell =1}}\left({\eta }_{\ell }{\overline{\beta }}_{\ell }{x}_{\ell }^0+{\vartheta }_{\ell }{\eta }_3{\lambda }_{\ell }-{\eta }_{\ell }{\vartheta }_{\ell }\right)y_{\ell }+\left(\stackrel{2}{\sum _{\ell =1}}{\eta }_{\ell }{\beta }_{\ell }{x}_{\ell }^0-{\eta }_3\phi \right)v\hfill \\ \hfill & =-\stackrel{2}{\sum _{\ell =1}}{\eta }_{\ell }{\gamma }_{\ell }\frac{{\left(x_{\ell }-x_{\ell }^0\right)}^2}{x_{\ell }}+\left({\vartheta }_1{\vartheta }_2{\lambda }_1(1-{\overline{\psi }}_a){\beta }_1{x}_1^0+{\vartheta }_1{\vartheta }_2{\lambda }_2(1-{\overline{\varphi }}_a){\beta }_2{x}_2^0\right.\hfill \\ \hfill & \left.-{\vartheta }_1{\vartheta }_2(1-{\overline{\varphi }}_a)(1-{\overline{\psi }}_a)\phi \right)v.\hfill \end{array}$$

(21)

Simplifying Equation (21), we get

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\Xi }_0}{dt}}& =-\stackrel{2}{\sum _{\ell =1}}{\eta }_{\ell }{\gamma }_{\ell }\frac{{\left(x_{\ell }-x_{\ell }^0\right)}^2}{x_{\ell }}+{\vartheta }_1{\vartheta }_2\phi \left((1-{\overline{\psi }}_a){\overline{\varphi }}_b+(1-{\overline{\varphi }}_a){\overline{\psi }}_b\right.\hfill \\ \hfill & \left.-(1-{\overline{\varphi }}_a)(1-{\overline{\psi }}_a)\right)v\hfill \\ \hfill & =-\stackrel{2}{\sum _{\ell =1}}{\eta }_{\ell }{\gamma }_{\ell }\frac{{\left(x_{\ell }-x_{\ell }^0\right)}^2}{x_{\ell }}-{\vartheta }_1{\vartheta }_2\phi \left(1-M\right)v.\hfill \end{array}$$

As a result, $\frac{d{\Xi }_0}{dt}\le 0$ if ${\mathcal{R}}_0\le 1$ for $x_{\ell },y_{\ell },v\in (0,\infty ),\ell =1,2$. Moreover, $\frac{d{\Xi }_0}{dt}=0$ when $x_{\ell }\left(t\right)=x_{\ell }^0$ and $v\left(t\right)=0,\ell =1,2,$ for all t. The solutions of System (6)–(8) tend to ${\Gamma }_0^{{}^{\prime }}$, which has elements with $v\left(t\right)=0,$ so $\dot{v}\left(t\right)=0.$ It follows from Equation (8) that

$$0=\dot{v}\left(t\right)=\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{y}_{\ell }\left(t\right)⟹y_{\ell }\left(t\right)=0,\ell =1,2.$$

Hence, ${\Gamma }_0^{{}^{\prime }}=$$\left\{{\Pi }_0\right\}$, and by applying LIP, we get that ${\Pi }_0$ is globally asymptotically stable (GAS). □

 Theorem 2.

For System (6)–(8), if ${\mathcal{R}}_0$ exceeds one (${\mathcal{R}}_0>1)$, then ${\Pi }^{*}$ is GAS.

 Proof.

Formulate a function ${\Xi }_1(x_1,y_1,x_2,y_2,v)$ as:

$${\Xi }_1=\stackrel{2}{\sum _{\ell =1}}{\check{\eta }}_{\ell }\left[x_{\ell }^{*}\Theta \left(\frac{x_{\ell }}{x_{\ell }^{*}}\right)+y_{\ell }^{*}\Theta \left(\frac{y_{\ell }}{y_{\ell }^{*}}\right)\right]+v_2^{*}\Theta \left(\frac{v_2}{v_2^{*}}\right),$$

where

$${\check{\eta }}_{\ell }=\frac{{\lambda }_{\ell }{\vartheta }_{\ell }{y}_{\ell }^{*}}{{\beta }_{\ell }{x}_{\ell }^{*}{v}^{*}},\ell =1,2.$$

(22)

Clearly, ${\Xi }_1(x_1,y_1,x_2,y_2,v)>0$ for all $x_1,y_1,x_2,y_2,v>0$, and ${\Xi }_1(x_1^{*},y_1^{*},x_2^{*},y_2^{*},v^{*})=0$. Calculating $\frac{d{\Xi }_1}{dt}$ along the trajectories of (6)–(8), we get

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\Xi }_1}{dt}}& =\stackrel{2}{\sum _{\ell =1}}{\check{\eta }}_{\ell }\left[\left(1-\frac{x_{\ell }^{*}}{x_{\ell }}\right){\dot{x}}_{\ell }+\left(1-\frac{y_{\ell }^{*}}{y_{\ell }}\right){\dot{y}}_{\ell }\right]+\left(1-\frac{v^{*}}{v}\right)\dot{v}\hfill \\ \hfill & =\stackrel{2}{\sum _{\ell =1}}{\check{\eta }}_{\ell }\left[\left(1-\frac{x_{\ell }^{*}}{x_{\ell }}\right)\left({\alpha }_{\ell }-{\gamma }_{\ell }{x}_{\ell }-{\beta }_{\ell }{x}_{\ell }v-{\overline{\beta }}_{\ell }{x}_{\ell }{y}_{\ell }\right)+\left(1-\frac{y_{\ell }^{*}}{y_{\ell }}\right)\left({\beta }_{\ell }{x}_{\ell }v+{{\overline{\beta }}_{\ell }x}_{\ell }{y}_{\ell }-{\vartheta }_{\ell }{y}_{\ell }\right)\right]\hfill \\ \hfill & +\left(1-\frac{v^{*}}{v}\right)\left(\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{y}_{\ell }-\phi v\right).\hfill \end{array}$$

(23)

Simplify Equation (23), and in the following steps we apply the following conditions for ${\Pi }^{*}$:

$$\begin{array}{cc}{\alpha }_{\ell }\hfill & ={\gamma }_{\ell }{x}_{\ell }^{*}+{\beta }_{\ell }{x}_{\ell }^{*}{v}^{*}+{\overline{\beta }}_{\ell }{x}_{\ell }^{*}{y}_{\ell }^{*},\ell =1,2,\hfill \end{array}$$

(24)

$$\begin{array}{cc}{\vartheta }_{\ell }{y}_{\ell }^{*}\hfill & ={\beta }_{\ell }{x}_{\ell }^{*}{v}^{*}+{\overline{\beta }}_{\ell }{x}_{\ell }^{*}{y}_{\ell }^{*},\hfill \end{array}$$

(25)

$$\begin{array}{cc}\phi v^{*}\hfill & =\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{y}_{\ell }^{*},\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill & =\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }\left({\beta }_{\ell }{x}_{\ell }^{*}{v}^{*}+{\overline{\beta }}_{\ell }{x}_{\ell }^{*}{y}_{\ell }^{*}\right).\hfill \end{array}$$

(27)

We get

$$\begin{array}{cc}{\displaystyle \frac{d{\Xi }_1}{dt}}\hfill & =-\stackrel{2}{\sum _{\ell =1}}{\gamma }_{\ell }{\check{\eta }}_{\ell }\frac{{\left(x_{\ell }-x_{\ell }^{*}\right)}^2}{x_{\ell }}+\stackrel{2}{\sum _{\ell =1}}{\check{\eta }}_{\ell }\left[\left(1-\frac{x_{\ell }^{*}}{x_{\ell }}\right)\left({\beta }_{\ell }{x}_{\ell }^{*}{v}^{*}+{\overline{\beta }}_{\ell }{x}_{\ell }^{*}{y}_{\ell }^{*}-{\beta }_{\ell }{x}_{\ell }v-{\overline{\beta }}_{\ell }{x}_{\ell }{y}_{\ell }\right)\right.\hfill \\ \hfill & \left.+\left(1-\frac{y_{\ell }^{*}}{y_{\ell }}\right)\left({\beta }_{\ell }{x}_{\ell }v+{{\overline{\beta }}_{\ell }x}_{\ell }{y}_{\ell }-{\vartheta }_{\ell }{y}_{\ell }\right)\right]+\left(1-\frac{v^{*}}{v}\right)\left(\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{y}_{\ell }-\phi v\right).\hfill \end{array}$$

Collecting the terms of the last equation, we have

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\Xi }_1}{dt}}& =\stackrel{2}{\underset{\ell =1}{-\sum }}{\gamma }_{\ell }{\check{\eta }}_{\ell }\frac{{\left(x_{\ell }-x_{\ell }^{*}\right)}^2}{x_{\ell }}+\stackrel{2}{\sum _{\ell =1}}{\check{\eta }}_{\ell }\left[{\beta }_{\ell }{x}_{\ell }^{*}{v}^{*}+{\overline{\beta }}_{\ell }{x}_{\ell }^{*}{y}_{\ell }^{*}-{\beta }_{\ell }{x}_{\ell }^{*}{v}^{*}\frac{x_{\ell }^{*}}{x_{\ell }}-{\overline{\beta }}_{\ell }{x}_{\ell }^{*}{y}_{\ell }^{*}\frac{x_{\ell }^{*}}{x_{\ell }}\right.\hfill \\ \hfill & \left.+{\beta }_{\ell }{x}_{\ell }^{*}v+{\overline{\beta }}_{\ell }{x}_{\ell }^{*}{y}_{\ell }-{\vartheta }_{\ell }{y}_{\ell }-{\beta }_{\ell }{x}_{\ell }v\frac{y_{\ell }^{*}}{y_{\ell }}-{{\overline{\beta }}_{\ell }x}_{\ell }{y}_{\ell }^{*}+{\vartheta }_{\ell }{y}_{\ell }^{*}\right]+\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{y}_{\ell }\hfill \\ \hfill & -\phi v-\frac{v^{*}}{v}\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{y}_{\ell }+\phi v^{*}.\hfill \end{array}$$

Then,

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\Xi }_1}{dt}}& =\stackrel{2}{\underset{\ell =1}{-\sum }}{\gamma }_{\ell }{\check{\eta }}_{\ell }\frac{{\left(x_{\ell }-x_{\ell }^{*}\right)}^2}{x_{\ell }}+\stackrel{2}{\sum _{\ell =1}}{\check{\eta }}_{\ell }\left[2{\beta }_{\ell }{x}_{\ell }^{*}{v}^{*}-{\beta }_{\ell }{x}_{\ell }^{*}{v}^{*}\frac{x_{\ell }^{*}}{x_{\ell }}+{\beta }_{\ell }{x}_{\ell }^{*}v+{\overline{\beta }}_{\ell }{x}_{\ell }^{*}{y}_{\ell }-{\beta }_{\ell }{x}_{\ell }v\frac{y_{\ell }^{*}}{y_{\ell }}\right]\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\check{\eta }}_{\ell }\left[2{\overline{\beta }}_{\ell }{x}_{\ell }^{*}{y}_{\ell }^{*}-{\overline{\beta }}_{\ell }{x}_{\ell }^{*}{y}_{\ell }^{*}\frac{x_{\ell }^{*}}{x_{\ell }}-{{\overline{\beta }}_{\ell }x}_{\ell }{y}_{\ell }^{*}\right]+\stackrel{2}{\sum _{\ell =1}}\left[{\lambda }_{\ell }-{\check{\eta }}_{\ell }\right]{\vartheta }_{\ell }{y}_{\ell }\hfill \\ \hfill & -\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{y}_{\ell }^{*}\frac{v}{v^{*}}-\frac{v^{*}}{v}\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{y}_{\ell }+\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }\left({\beta }_{\ell }{x}_{\ell }^{*}{v}^{*}+{\overline{\beta }}_{\ell }{x}_{\ell }^{*}{y}_{\ell }^{*}\right).\hfill \end{array}$$

(28)

Simplifying Equation (28), we obtain

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\Xi }_1}{dt}}& =-\stackrel{2}{\sum _{\ell =1}}{\gamma }_{\ell }{\check{\eta }}_{\ell }\frac{{\left(x_{\ell }-x_{\ell }^{*}\right)}^2}{x_{\ell }}+\stackrel{2}{\sum _{\ell =1}}{\check{\eta }}_{\ell }{\beta }_{\ell }{x}_{\ell }^{*}{v}^{*}\left[2-\frac{x_{\ell }^{*}}{x_{\ell }}+\frac{v}{v^{*}}-\frac{y_{\ell }}{y_{\ell }^{*}}-\frac{x_{\ell }v{y}_{\ell }^{*}}{x_{\ell }^{*}{v}^{*}{y}_{\ell }}\right]\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\check{\eta }}_{\ell }{\overline{\beta }}_{\ell }{x}_{\ell }^{*}{y}_{\ell }^{*}\left[2-\frac{x_{\ell }^{*}}{x_{\ell }}-\frac{x_{\ell }}{x_{\ell }^{*}}\right]\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\beta }_{\ell }{x}_{\ell }^{*}{v}^{*}\left[1+\frac{y_{\ell }}{y_{\ell }^{*}}-\frac{v}{v^{*}}-\frac{y_{\ell }}{y_{\ell }^{*}}\frac{v^{*}}{v}\right]+\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\overline{\beta }}_{\ell }{x}_{\ell }^{*}{y}_{\ell }^{*}\left[1+\frac{y_{\ell }}{y_{\ell }^{*}}-\frac{v}{v^{*}}-\frac{y_{\ell }{v}^{*}}{y_{\ell }^{*}v}\right].\hfill \end{array}$$

The last equation can be reduced to

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\Xi }_1}{dt}}& =-\stackrel{2}{\sum _{\ell =1}}{\gamma }_{\ell }{\check{\eta }}_{\ell }\frac{{\left(x_{\ell }-x_{\ell }^{*}\right)}^2}{x_{\ell }}+\stackrel{2}{\sum _{\ell =1}}{\check{\eta }}_{\ell }{\beta }_{\ell }{x}_{\ell }^{*}{v}^{*}\left[3-\frac{x_{\ell }^{*}}{x_{\ell }}-\frac{x_{\ell }v{y}_{\ell }^{*}}{x_{\ell }^{*}{v}^{*}{y}_{\ell }}-\frac{y_{\ell }}{y_{\ell }^{*}}\frac{v^{*}}{v}\right]\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\check{\eta }}_{\ell }{\overline{\beta }}_{\ell }{x}_{\ell }^{*}{y}_{\ell }^{*}\left[2-\frac{x_{\ell }^{*}}{x_{\ell }}-\frac{x_{\ell }}{x_{\ell }^{*}}\right]\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}\left({\lambda }_{\ell }{\beta }_{\ell }{x}_{\ell }^{*}{v}^{*}+{\lambda }_{\ell }{\overline{\beta }}_{\ell }{x}_{\ell }^{*}{y}_{\ell }^{*}\right)\left[1+\frac{y_{\ell }}{y_{\ell }^{*}}-\frac{v}{v^{*}}-\frac{y_{\ell }{v}^{*}}{y_{\ell }^{*}v}\right]\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\check{\eta }}_{\ell }{\beta }_{\ell }{x}_{\ell }^{*}{v}^{*}\left[\frac{v}{v^{*}}-\frac{y_{\ell }}{y_{\ell }^{*}}+\frac{y_{\ell }{v}^{*}}{y_{\ell }^{*}v}-1\right].\hfill \end{array}$$

Now, from Equations (22) and (25), we have ${\check{\eta }}_{\ell }{\beta }_{\ell }$$x_{\ell }^{*}{v}^{*}={\lambda }_{\ell }\left[{\beta }_{\ell }{x}_{\ell }^{*}{v}^{*}+{\overline{\beta }}_{\ell }{x}_{\ell }^{*}{y}_{\ell }^{*}\right],\ell =1,2$; then, we get

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\Xi }_1}{dt}}& =-\stackrel{2}{\sum _{\ell =1}}{\gamma }_{\ell }{\check{\eta }}_{\ell }\frac{{\left(x_{\ell }-x_{\ell }^{*}\right)}^2}{x_{\ell }}+\stackrel{2}{\sum _{\ell =1}}{\check{\eta }}_{\ell }{\beta }_{\ell }{x}_{\ell }^{*}{v}^{*}\left[3-\frac{x_{\ell }^{*}}{x_{\ell }}-\frac{x_{\ell }v{y}_{\ell }^{*}}{x_{\ell }^{*}{v}^{*}{y}_{\ell }}-\frac{y_{\ell }}{y_{\ell }^{*}}\frac{v^{*}}{v}\right]\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\check{\eta }}_{\ell }{\overline{\beta }}_{\ell }{x}_{\ell }^{*}{y}_{\ell }^{*}\left[2-\frac{x_{\ell }^{*}}{x_{\ell }}-\frac{x_{\ell }}{x_{\ell }^{*}}\right]+\stackrel{2}{\sum _{\ell =1}}{\check{\eta }}_{\ell }{\beta }_{\ell }{x}_{\ell }^{*}{v}^{*}\left[1+\frac{y_{\ell }}{y_{\ell }^{*}}-\frac{v}{v^{*}}-\frac{y_{\ell }{v}^{*}}{y_{\ell }^{*}v}\right]\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\check{\eta }}_{\ell }{\beta }_{\ell }{x}_{\ell }^{*}{v}^{*}\left[\frac{v}{v^{*}}-\frac{y_{\ell }}{y_{\ell }^{*}}+\frac{y_{\ell }{v}^{*}}{y_{\ell }^{*}v}-1\right];\hfill \end{array}$$

hence,

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\Xi }_1}{dt}}& =-\stackrel{2}{\sum _{\ell =1}}{\gamma }_{\ell }{\check{\eta }}_{\ell }\frac{{\left(x_{\ell }-x_{\ell }^{*}\right)}^2}{x_{\ell }}+\stackrel{2}{\sum _{\ell =1}}{\check{\eta }}_{\ell }{\beta }_{\ell }{x}_{\ell }^{*}{v}^{*}\left[3-\frac{x_{\ell }^{*}}{x_{\ell }}-\frac{x_{\ell }v{y}_{\ell }^{*}}{x_{\ell }^{*}{v}^{*}{y}_{\ell }}-\frac{y_{\ell }}{y_{\ell }^{*}}\frac{v^{*}}{v}\right]\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\check{\eta }}_{\ell }{\overline{\beta }}_{\ell }{x}_{\ell }^{*}{y}_{\ell }^{*}\left[2-\frac{x_{\ell }^{*}}{x_{\ell }}-\frac{x_{\ell }}{x_{\ell }^{*}}\right],\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill 3& \le \frac{x_{\ell }^{*}}{x_{\ell }}+\frac{x_{\ell }v{y}_{\ell }^{*}}{x_{\ell }^{*}{v}^{*}{y}_{\ell }}+\frac{y_{\ell }{v}^{*}}{y_{\ell }^{*}v},\ell =1,2,\hfill \\ \hfill 2& \le \frac{x_{\ell }^{*}}{x_{\ell }}+\frac{x_{\ell }}{x_{\ell }^{*}},\ell =1,2.\hfill \end{array}$$

In this case, $x_{\ell }^{*},y_{\ell }^{*},v^{*}>0$ if ${\mathcal{R}}_0>1;$ then ${\displaystyle \frac{d{\Xi }_1}{dt}}\le 0$ for all $x_{\ell },y_{\ell }$,$v>0,\ell =1,2$. It can be seen that ${\displaystyle \frac{d{\Xi }_1}{dt}}=0$ if and only if $(x_1\left(t\right),y_1\left(t\right),x_2\left(t\right),y_2\left(t\right),v\left(t\right))=(x_1^{*},y_1^{*},x_2^{*},y_2^{*},v^{*})$. Therefore, the solutions of System (6)–(8) tend to ${\Gamma }_1^{{}^{\prime }}=\left\{{\Pi }^{*}\right\}$ and ${\Pi }^{*}$ is GAS when ${\mathcal{R}}_0>1$ according to LIP. □

3. The Model with Delay

We incorporate into System (6)–(8) the discrete time delays. The model takes the form:

$$\begin{array}{cc}{\dot{x}}_{\ell }\hfill & ={\alpha }_{\ell }-{\gamma }_{\ell }{x}_{\ell }-{\beta }_{\ell }{x}_{\ell }v-{\overline{\beta }}_{\ell }{x}_{\ell }{y}_{\ell },\ell =1,2,\hfill \end{array}$$

(29)

$$\begin{array}{cc}{\dot{y}}_{\ell }\hfill & ={\beta }_{\ell }{e}^{-{ϵ}_{\ell }{\sigma }_{\ell }}{x}_{\ell }(t-{\sigma }_{\ell })v(t-{\sigma }_{\ell })+{\overline{\beta }}_{\ell }{e}^{-{ϵ}_{\ell }{\sigma }_{\ell }}{x}_{\ell }(t-{\sigma }_{\ell })y_{\ell }(t-{\sigma }_{\ell })-{\vartheta }_{\ell }{y}_{\ell },\ell =1,2,\hfill \end{array}$$

(30)

$$\begin{array}{cc}\dot{v}\hfill & =\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}{y}_{\ell }(t-{\varpi }_{\ell })-\phi v.\hfill \end{array}$$

(31)

In this model, the loss of the cells during the delay period $[t-{\sigma }_{\ell },t]$ is given by $e^{-{ϵ}_{\ell }{\sigma }_{\ell }}$, where ${ϵ}_{\ell }>0,\ell =1,2$. The factor $e^{-{\epsilon }_{\ell }{\varpi }_{\ell }}$ represents the loss of the infected cells during the delay period $[t-{\varpi }_{\ell },t]$, where ${\epsilon }_{\ell }>0,\ell =1,2$.

3.1. Properties of Solutions

Let $\mathrm{{\rm Y}}=max\{{\sigma }_1,{\varpi }_1,{\sigma }_2,{\varpi }_2\},$ and denote by $\mathbb{C}$ the Banach space of continuous functions mapping the interval $[-\mathrm{{\rm Y}},0]$ into ${\mathbb{R}}_{\ge 0}^5$. For System (29)–(31), we consider the initial conditions

$$\begin{array}{cc}\hfill x_{\ell }\left(\xi \right)& ={\chi }_{\ell }\left(\xi \right),y_{\ell }\left(\xi \right)={\chi }_{\ell +2}\left(\xi \right),v\left(\xi \right)={\chi }_5\left(\xi \right),\ell =1,2,\hfill \\ \hfill {\chi }_j\left(\xi \right)& \ge 0,j=1,2,\dots ,5,\xi \in [-\mathrm{{\rm Y}},0],\hfill \end{array}$$

(32)

where $({\chi }_1\left(\xi \right),\dots ,{\chi }_5\left(\xi \right))\in \mathbb{C}([-\mathrm{{\rm Y}},0],{\mathbb{R}}_{\ge 0}^5)$. Thus, there exists a unique solution for System (29)–(31) with initial conditions (32) (see [43]).

 Lemma 5.

Solutions of System (29)–(31) that satisfy the initial conditions (32) are non-negative and are ultimately bounded by $t\ge 0$.

 Proof.

We have ${\dot{x}}_{\ell }{\mid }_{x_{\ell }=0}={\alpha }_{\ell }>0$; hence, $x_{\ell }\left(t\right)>0$. From Equations (29)–(31), we have:

$$\begin{array}{cc}\hfill y_{\ell }\left(t\right)& ={\chi }_{\ell +2}\left(0\right)e^{-{\vartheta }_{\ell }t}+\underset{0}{\overset{t}{\int }}\left(\beta {}_{\ell }{e}^{-{\vartheta }_{\ell }(t-\xi )}{x}_{\ell }\left(\xi -{\sigma }_{\ell }\right)v\left(\xi -{\sigma }_{\ell }\right)\overline{\beta }{}_{\ell }{e}^{-{\vartheta }_{\ell }(t-\xi )}{x}_{\ell }\left(\xi -{\sigma }_{\ell }\right)y_{\ell }\left(\xi -{\sigma }_{\ell }\right)\right)d\xi ,\ell =1,2,\hfill \\ \hfill v\left(t\right)& ={\chi }_5\left(0\right)e^{-\phi t}+\sum _{\ell =1}^2\underset{0}{\overset{t}{\int }}{e}^{-\phi \left(t-\xi \right)}{\lambda }_{\ell }{\vartheta }_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}{y}_{\ell }\left(\xi -{\varpi }_{\ell }\right)d\xi .\hfill \end{array}$$

These show that $y_{\ell }\left(t\right)\ge 0$, $\ell =1,2$, and $v\left(t\right)\ge 0$ for all $t\ge 0$.

Define ${\tilde{F}}_{\ell }\left(t\right)=e^{-{ϵ}_{\ell }{\sigma }_{\ell }}{x}_{\ell }(t-{\sigma }_{\ell })+y_{\ell }\left(t\right)$, $\ell =1,2,$ and from Equation (29), we have ${\displaystyle \underset{t\to \infty }{lim}}sup{x}_{\ell }\left(t\right)\le \frac{{\alpha }_{\ell }}{{\gamma }_{\ell }},$$\ell =1,2$. As a result,

$$\begin{array}{cc}\hfill \frac{d{\tilde{F}}_{\ell }\left(t\right)}{dt}& ={\alpha }_{\ell }{e}^{-{ϵ}_{\ell }{\sigma }_{\ell }}-{\gamma }_{\ell }{e}^{-{ϵ}_{\ell }{\sigma }_{\ell }}{x}_{\ell }(t-{\sigma }_{\ell })-{\vartheta }_{\ell }{y}_{\ell }\left(t\right)\hfill \\ \hfill & \le {\alpha }_{\ell }{e}^{-{ϵ}_{\ell }{\sigma }_{\ell }}-{\tilde{\delta }}_{\ell }\left(e^{-{ϵ}_{\ell }{\sigma }_{\ell }}{x}_{\ell }(t-{\sigma }_{\ell })+y_{\ell }\left(t\right)\right)\le {\alpha }_{\ell }-{\tilde{\delta }}_{\ell }{\tilde{F}}_{\ell }\left(t\right),\hfill \end{array}$$

where ${\tilde{\delta }}_{\ell }=min\{{\gamma }_{\ell },{\vartheta }_{\ell }\},\ell =1,2$. Hence, ${\displaystyle \underset{t\to \infty }{lim}}sup{\tilde{F}}_{\ell }\left(t\right)\le {\tilde{\zeta }}_{\ell }$, where ${\tilde{\zeta }}_{\ell }={\displaystyle \frac{{\alpha }_{\ell }}{{\tilde{\delta }}_{\ell }}},\ell =1,2.$ Thus, we get ${\displaystyle \underset{t\to \infty }{lim}}sup{y}_{\ell }\left(t\right)\le$${\tilde{\zeta }}_{\ell }$ for all $t\ge 0$. On the other hand,

$$\dot{v}\left(t\right)=\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}{y}_{\ell }(t-{\varpi }_{\ell })-\phi v\left(t\right)\le \stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}{\tilde{\zeta }}_{\ell }-\phi v\left(t\right)\le \stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{\tilde{\zeta }}_{\ell }-\phi v\left(t\right);$$

then $0\le v\left(t\right)\le$${\tilde{\zeta }}_3$ for all $t\ge 0$ if $v\left(0\right)\le$ ${\tilde{\zeta }}_3$, where ${\tilde{\zeta }}_3={\displaystyle \sum _{\ell =1}^2}{\displaystyle \frac{{\lambda }_{\ell }{\vartheta }_{\ell }{\tilde{\zeta }}_{\ell }}{\phi }}$. Consequently, $x_{\ell }\left(t\right),y_{\ell }\left(t\right),$$\ell =1,2$, and $v\left(t\right)$ are ultimately bounded. □

3.2. Equilibria

To calculate the equilibrium points of the System (29)–(31), we solve the following system:

$$\begin{array}{cc}0\hfill & ={\alpha }_{\ell }-{\gamma }_{\ell }{x}_{\ell }-{\beta }_{\ell }{x}_{\ell }v-{\overline{\beta }}_{\ell }{x}_{\ell }{y}_{\ell },\ell =1,2,\hfill \end{array}$$

(33)

$$\begin{array}{cc}0\hfill & ={\beta }_{\ell }{e}^{-{ϵ}_{\ell }{\sigma }_{\ell }}{x}_{\ell }v+{\overline{\beta }}_{\ell }{e}^{-{ϵ}_{\ell }{\sigma }_{\ell }}{x}_{\ell }{y}_{\ell }-{\vartheta }_{\ell }{y}_{\ell },\ell =1,2,\hfill \end{array}$$

(34)

$$\begin{array}{cc}0\hfill & =\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}{y}_{\ell }-\phi v.\hfill \end{array}$$

(35)

As a result of the calculations, two equilibrium points can be found:

(i)

   Infection-free equilibrium (IFE) ${\Pi }_0=(x_1^0,0,x_2^0,0,0),$ where $x_{\ell }^0=\frac{{\alpha }_{\ell }}{{\gamma }_{\ell }},\ell =1,2.$

(ii)

  Infection-present equilibrium (IPE) $\tilde{\Pi }=({\tilde{x}}_1,{\tilde{y}}_1,{\tilde{x}}_2,{\tilde{y}}_2,\tilde{v})$ with the following definitions of each component:

$${\tilde{x}}_{\ell }=\frac{{\alpha }_{\ell }}{{\gamma }_{\ell }+{\beta }_{\ell }\tilde{v}+{\overline{\beta }}_{\ell }{\tilde{y}}_{\ell }},{\tilde{y}}_{\ell }=\frac{-{\tilde{B}}_{\ell }+\sqrt{{\tilde{B}}_{\ell }^2+4{\tilde{A}}_{\ell }{\tilde{C}}_{\ell }}}{2{\tilde{A}}_{\ell }},$$

(36)

and

$${\tilde{A}}_{\ell }={\vartheta }_{\ell }{\overline{\beta }}_{\ell },{\tilde{B}}_{\ell }={\vartheta }_{\ell }{\beta }_{\ell }\tilde{v}+{\vartheta }_{\ell }{\gamma }_{\ell }-{\overline{\beta }}_{\ell }{\alpha }_{\ell }{e}^{-{ϵ}_{\ell }{\sigma }_{\ell }},{\tilde{C}}_{\ell }={\beta }_{\ell }{\alpha }_{\ell }{e}^{-{ϵ}_{\ell }{\sigma }_{\ell }}\tilde{v},$$

(37)

where $\ell =1,2$ and $\tilde{v}$ satisfies the following equation

$$\phi \tilde{v}=\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}{\tilde{y}}_{\ell }.$$

(38)

The basic reproduction number ${\mathcal{R}}_0$: As in the previous method in Section 2.3, we calculate the basic reproduction number ${\mathcal{R}}_0$ of System (29)–(31) by implementing the next-generation-matrix method [36] as follows:

$${\mathcal{R}}_0=\rho \left(\tilde{\mathcal{F}}{\tilde{\nabla }}^{-1}\right),$$

where the Jacobian of the matrix of new infection terms and the Jacobian of the matrix of the other terms at IFE are given, respectively, by

$$\begin{array}{cc}\hfill \tilde{\mathcal{F}}& =\left[\begin{array}{ccc}{\overline{\beta }}_1{e}^{-{ϵ}_1{\sigma }_1}{x}_1^0& 0& {\beta }_1{e}^{-{ϵ}_1{\sigma }_1}{x}_1^0\\ 0& {\overline{\beta }}_2{e}^{-{ϵ}_2{\sigma }_2}{x}_2^0& {\beta }_2{e}^{-{ϵ}_2{\sigma }_2}{x}_2^0\\ 0& 0& 0\end{array}\right],\hfill \\ \hfill \tilde{\nabla }& =\left[\begin{array}{ccc}{\vartheta }_1& 0& 0\\ 0& {\vartheta }_2& 0\\ -{\lambda }_1{\vartheta }_1{e}^{-{\epsilon }_1{\varpi }_1}& -{\lambda }_2{\vartheta }_2{e}^{-{\epsilon }_2{\varpi }_2}& \phi \end{array}\right].\hfill \end{array}$$

Thus,

$${\mathcal{R}}_0=\frac{1}{2}\left(\widehat{\Phi }+\widehat{\Psi }+\sqrt{{\left(\widehat{\Phi }-\widehat{\Psi }\right)}^2+4{\widehat{\varphi }}_b{\widehat{\psi }}_b}\right),$$

(39)

where

$$\begin{array}{cc}\hfill \widehat{\Phi }& ={\widehat{\varphi }}_a+{\widehat{\varphi }}_b,\widehat{\Psi }={\widehat{\psi }}_a+{\widehat{\psi }}_b,\hfill \\ \hfill {\widehat{\varphi }}_a& =\frac{{\overline{\beta }}_1{e}^{-{ϵ}_1{\sigma }_1}{x}_1^0}{{\vartheta }_1},{\widehat{\varphi }}_b=\frac{{\beta }_1{e}^{-({ϵ}_1{\sigma }_1+{\epsilon }_1{\varpi }_1)}{x}_1^0{\lambda }_1}{\phi },{\widehat{\psi }}_a=\frac{{\overline{\beta }}_2{e}^{-{ϵ}_2{\sigma }_2}{x}_2^0}{{\vartheta }_2},{\widehat{\psi }}_b=\frac{{\beta }_2{e}^{-({ϵ}_2{\sigma }_2+{\epsilon }_2{\varpi }_2)}{x}_2^0{\lambda }_2}{\phi }.\hfill \end{array}$$

(40)

In the following Lemmas, we show the condition when there is a positive value for $\tilde{v}$.

 Lemma 6.

Suppose that ${\mathcal{R}}_0>1.$$If$ $\widehat{\Phi }<1,\widehat{\Psi }<1,$ then $\tilde{M}=\frac{{\widehat{\varphi }}_b}{1-{\widehat{\varphi }}_a}+\frac{{\widehat{\psi }}_b}{1-{\widehat{\psi }}_a}>1$.

 Proof.

Similar to the proof of Lemma 2. □

 Lemma 7.

Suppose that ${\mathcal{R}}_0>1$, then the IPE $\tilde{\Pi }$ exists.

 Proof.

First, we have that there exists any equilibrium satisfying Equations (33)–(34). In case of the equilibrium $\tilde{\Pi }$, we have $v\ne 0;$ then, from Equation (35), we find

$$\stackrel{2}{\sum _{\ell =1}}\frac{{\lambda }_{\ell }{\vartheta }_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}{y}_{\ell }}{\phi }-v=0.$$

Substituting from Equations (33)–(35) into the last equation, we get

$$\stackrel{2}{\sum _{\ell =1}}{e}^{-\left({ϵ}_{\ell }{\sigma }_{\ell }+{\epsilon }_{\ell }{\varpi }_{\ell }\right)}\left(\frac{{\lambda }_{\ell }{\alpha }_{\ell }}{\phi }-\frac{{\lambda }_{\ell }{\gamma }_{\ell }{x}_{\ell }}{\phi }\right)-v=0.$$

Since $x_{\ell }=x_{\ell }\left(v\right)$, $y_{\ell }=y_{\ell }\left(v\right),\ell =1,2,$ then we can define a function $\tilde{\mathcal{G}}\left(v\right)$ as:

$$\tilde{\mathcal{G}}\left(v\right)=\stackrel{2}{\sum _{\ell =1}}{e}^{-\left({ϵ}_{\ell }{\sigma }_{\ell }+{\epsilon }_{\ell }{\varpi }_{\ell }\right)}\left(\frac{{\lambda }_{\ell }{\alpha }_{\ell }}{\phi }-\frac{{\lambda }_{\ell }{\gamma }_{\ell }{x}_{\ell }}{\phi }\right)-v,$$

in which $x_{\ell }$, $y_{\ell }$ satisfy Equations (33)–(35) for $\ell =1,2$.

Now we need to show that ∃$\tilde{v}>0$ such that $\tilde{\mathcal{G}}\left(\tilde{v}\right)=0$ as follow:

If $v={\widehat{v}}^{*}={\displaystyle \sum _{\ell =1}^2}2\frac{{\lambda }_{\ell }{\alpha }_{\ell }{e}^{-\left({ϵ}_{\ell }{\sigma }_{\ell }+{\epsilon }_{\ell }{\varpi }_{\ell }\right)}}{\phi }>0,$ then $x_{\ell }\left({\widehat{v}}^{*}\right)>0,y_{\ell }\left({\widehat{v}}^{*}\right)>0,\ell =1,2,$ and $\tilde{\mathcal{G}}\left({\widehat{v}}^{*}\right)=-{\displaystyle \sum _{\ell =1}^2}\frac{{\lambda }_{\ell }{\gamma }_{\ell }{x}_{\ell }\left({\widehat{v}}^{*}\right)e^{-\left({ϵ}_{\ell }{\sigma }_{\ell }+{\epsilon }_{\ell }{\varpi }_{\ell }\right)}}{\phi }<0.$

Next, by calculating $\tilde{\mathcal{G}}\left(0\right)$ and ${\tilde{\mathcal{G}}}^{\prime }\left(0\right)$, we get:

$$\begin{array}{cc}\hfill \tilde{\mathcal{G}}\left(0\right)& =\stackrel{2}{\sum _{\ell =1}}{e}^{-\left({ϵ}_{\ell }{\sigma }_{\ell }+{\epsilon }_{\ell }{\varpi }_{\ell }\right)}\left(\frac{{\lambda }_{\ell }{\gamma }_{\ell }{x}_{\ell }^0}{\phi }-\frac{{\lambda }_{\ell }{\gamma }_{\ell }}{\phi }{x}_{\ell }\left(0\right)\right)\hfill \\ \hfill & ={\displaystyle \frac{{\widehat{\varphi }}_b{\gamma }_1}{{\beta }_1}}\left(1-\frac{2}{2-\left(1-{\widehat{\varphi }}_a\right)+\sqrt{{\left(1-{\widehat{\varphi }}_a\right)}^2}}\right)+{\displaystyle \frac{{\widehat{\psi }}_b{\gamma }_2}{{\beta }_2}}\left(1-\frac{2}{2-\left(1-{\widehat{\psi }}_a\right)+\sqrt{{\left(1-{\widehat{\psi }}_a\right)}^2}}\right),\hfill \\ \hfill {\tilde{\mathcal{G}}}^{\prime }\left(0\right)& =\frac{{\widehat{\varphi }}_b}{2{\widehat{\varphi }}_a}\left(\frac{{\widehat{\varphi }}_a+1}{\sqrt{{\left({\widehat{\varphi }}_a-1\right)}^2}}-1\right)+\frac{{\widehat{\psi }}_b}{2{\widehat{\psi }}_a}\left(\frac{{\widehat{\psi }}_a+1}{\sqrt{{\left({\widehat{\psi }}_a-1\right)}^2}}-1\right)-1.\hfill \end{array}$$

We calculate all possible cases of functions $\tilde{\mathcal{G}}\left(0\right),$${\tilde{\mathcal{G}}}^{\prime }\left(0\right)$, and the results are provided by

Table 2. Functions $\tilde{\mathcal{G}}\left(0\right)$ and ${\tilde{\mathcal{G}}}^{\prime }\left(0\right)$ and their corresponding values for different conditions.

Case	Conditions	$\tilde{\mathcal{G}}\left(0\right)$	${\tilde{\mathcal{G}}}^{\prime }\left(0\right)$
1	${\widehat{\varphi }}_a=1,{\widehat{\psi }}_a\le 1$	$0$	$+\infty$
2	${\widehat{\varphi }}_a\le 1,{\widehat{\psi }}_a=1$	$0$	$+\infty$
3	${\widehat{\varphi }}_a<1,{\widehat{\psi }}_a<1$	$0$	$\tilde{M}-1>0$(from Lemma 6)
4	${\widehat{\varphi }}_a\le 1,{\widehat{\psi }}_a>1$	${\displaystyle \frac{{\widehat{\psi }}_b{\gamma }_2}{{\beta }_2}}\left({\displaystyle \frac{{\widehat{\psi }}_a-1}{{\widehat{\psi }}_a}}\right)>0$	−
5	${\widehat{\varphi }}_a>1,{\widehat{\psi }}_a\le 1$	${\displaystyle \frac{{\widehat{\varphi }}_b{\gamma }_1}{{\beta }_1}}\left({\displaystyle \frac{{\widehat{\varphi }}_a-1}{{\widehat{\varphi }}_a}}\right)>0$	−
6	${\widehat{\varphi }}_a>1,{\widehat{\psi }}_a>1$	${\displaystyle \frac{{\widehat{\varphi }}_b{\gamma }_1}{{\beta }_1}}\left({\displaystyle \frac{{\widehat{\varphi }}_a-1}{{\widehat{\varphi }}_a}}\right)+{\displaystyle \frac{{\widehat{\psi }}_b{\gamma }_2}{{\beta }_2}}\left({\displaystyle \frac{{\widehat{\psi }}_a-1}{{\widehat{\psi }}_a}}\right)>0$	−

. As shown in Table 2, the function $\tilde{\mathcal{G}}\left(v\right)$ in Cases (1–3) is strictly increasing at the point $v=0,$ and $\tilde{\mathcal{G}}\left(0\right)>0$ in Cases (4–6), while $\tilde{\mathcal{G}}\left(\check{v}\right)$ has a negative value. This means that in all possible cases ∃$\tilde{v}\in (0,{\widehat{v}}^{*})$ satisfying $\tilde{\mathcal{G}}\left(\tilde{v}\right)=0$ if the condition ${\mathcal{R}}_0>1.$ Therefore, from Equations (36)–(37), we have ${\tilde{x}}_{\ell }>0$, ${\tilde{y}}_{\ell }>0$, $\tilde{v}>0,\ell =1,2$. Thus, the disease equilibrium $\tilde{\Pi }$ exists when ${\mathcal{R}}_0>1.$ □

From the above, we obtain:

(i)

   If ${\mathcal{R}}_0\le 1,$ then there will be only one equilibrium ${\Pi }_0$;

(ii)

  If ${\mathcal{R}}_0>1,$ then there will be two equilibria ${\Pi }_0$ and $\tilde{\Pi }$.

3.3. Global Properties

For System (29)–(31), as in Section 2.4, we verify the global asymptotic stability of both ${\Pi }_0$ and $\tilde{\Pi }$.

Let ${\tilde{\Gamma }}_{\ell }^{{}^{\prime }}$ be the largest invariant subset of ${\tilde{\Gamma }}_{\ell }=\left\{(x_1,x_{2,}{y}_1,y_2,v):\frac{d{\tilde{\Xi }}_{\ell }}{dt}=0\right\}$, where $\ell =1,2$. For the purpose of proving global stability, we utilize the following Lemma:

 Lemma 8.

Suppose that${\mathcal{R}}_0\le 1$; then,

 (i)

   ${\widehat{\varphi }}_a\le 1,{\widehat{\varphi }}_b\le 1,$${\widehat{\psi }}_a\le 1,$ and ${\widehat{\psi }}_b\le 1$;

 (ii)

  If $\widehat{M}=\widehat{\Phi }+\widehat{\Psi }-\widehat{\Phi }\widehat{\Psi }+{\widehat{\varphi }}_b{\widehat{\psi }}_b;$ thus, $0<\widehat{M}\le 1.$

 Proof.

Similar to the proof of Lemma 4. □

 Theorem 3.

For System (29)–(31), if ${\mathcal{R}}_0\le 1$, then ${\Pi }_0$ is GAS.

Proof. 

Let ${\mathcal{R}}_0\le 1$ and construct a function ${\tilde{\Xi }}_0(x_1,y_1,x_2,y_2,v)$ as:

$$\begin{array}{cc}\hfill {\tilde{\Xi }}_0& =\stackrel{2}{\sum _{\ell =1}}{\tilde{\eta }}_{\ell }\left[x_{\ell }^0\Theta \left(\frac{x_{\ell }}{x_{\ell }^0}\right)+e^{{ϵ}_{\ell }{\sigma }_{\ell }}{y}_{\ell }\right]\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\tilde{\eta }}_{\ell }\underset{0}{\overset{{\sigma }_{\ell }}{\int }}\left({\beta }_{\ell }{x}_{\ell }\left(t-\xi \right)v\left(t-\xi \right)+{\overline{\beta }}_{\ell }{x}_{\ell }\left(t-\xi \right)y_{\ell }\left(t-\xi \right)\right)d\xi \hfill \\ \hfill & +{\tilde{\eta }}_3\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}\underset{0}{\overset{{\varpi }_{\ell }}{\int }}{y}_{\ell }\left(t-\xi \right)d\xi +{\tilde{\eta }}_{3}v,\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill {\tilde{\eta }}_1& ={\vartheta }_1{\vartheta }_2{\lambda }_1{e}^{-({ϵ}_1{\sigma }_1+{\epsilon }_1{\varpi }_1)}(1-{\widehat{\psi }}_a),\hfill \\ \hfill {\tilde{\eta }}_2& ={\vartheta }_1{\vartheta }_2{\lambda }_2{e}^{-({ϵ}_2{\sigma }_2+{\epsilon }_2{\varpi }_2)}(1-{\widehat{\varphi }}_a),\hfill \\ \hfill {\tilde{\eta }}_3& ={\vartheta }_1{\vartheta }_2(1-{\widehat{\varphi }}_a)(1-{\widehat{\psi }}_a).\hfill \end{array}$$

(41)

Clearly, ${\tilde{\Xi }}_0(x_1,y_1,x_2,y_2,v)>0$ for all $x_1,y_1,x_2,y_2,v>0$, and ${\tilde{\Xi }}_0(x_1^0,0,x_2^0,0,0)=0$. We calculate ${\displaystyle \frac{d{\tilde{\Xi }}_0}{dt}}$ as:

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\tilde{\Xi }}_0}{dt}}& =\stackrel{2}{\sum _{\ell =1}}{\tilde{\eta }}_{\ell }\left(1-\frac{x_{\ell }^0}{x_{\ell }}\right){\dot{x}}_{\ell }+\stackrel{2}{\sum _{\ell =1}}{\tilde{\eta }}_{\ell }{e}^{{ϵ}_{\ell }{\sigma }_{\ell }}{\dot{y}}_{\ell }\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\tilde{\eta }}_{\ell }\left({\beta }_{\ell }{x}_{\ell }v-{\beta }_{\ell }{x}_2\left(t-{\sigma }_{\ell }\right)v\left(t-{\sigma }_{\ell }\right)+{\overline{\beta }}_{\ell }{x}_{\ell }{y}_{\ell }-{\overline{\beta }}_{\ell }{x}_{\ell }\left(t-{\sigma }_{\ell }\right)y_{\ell }\left(t-{\sigma }_{\ell }\right)\right)\hfill \\ \hfill & +{\tilde{\eta }}_3\stackrel{2}{\sum _{\ell =1}}{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}{\lambda }_{\ell }{\vartheta }_{\ell }\left(y_{\ell }-y_{\ell }\left(t-{\varpi }_{\ell }\right)\right)+{\tilde{\eta }}_3\dot{v}.\hfill \end{array}$$

Using Equations (29)–(31), we find

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\tilde{\Xi }}_0}{dt}}& =\stackrel{2}{\sum _{\ell =1}}{\tilde{\eta }}_{\ell }\left(1-\frac{x_{\ell }^0}{x_{\ell }}\right)\left({\alpha }_{\ell }-{\gamma }_{\ell }{x}_{\ell }-{\beta }_{\ell }{x}_{\ell }v-{\overline{\beta }}_{\ell }{x}_{\ell }{y}_{\ell }\right)\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\tilde{\eta }}_{\ell }{e}^{{ϵ}_{\ell }{\sigma }_{\ell }}\left({\beta }_{\ell }{e}^{-{ϵ}_{\ell }{\sigma }_{\ell }}{x}_{\ell }(t-{\sigma }_{\ell })v(t-{\sigma }_{\ell })+{\overline{\beta }}_{\ell }{e}^{-{ϵ}_{\ell }{\sigma }_{\ell }}{x}_{\ell }(t-{\sigma }_{\ell })y_{\ell }(t-{\sigma }_{\ell })-{\vartheta }_{\ell }{y}_{\ell }\right)\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\tilde{\eta }}_{\ell }\left({\beta }_{\ell }{x}_{\ell }v-{\beta }_{\ell }{x}_2\left(t-{\sigma }_{\ell }\right)v\left(t-{\sigma }_{\ell }\right)+{\overline{\beta }}_{\ell }{x}_{\ell }{y}_{\ell }-{\overline{\beta }}_{\ell }{x}_{\ell }\left(t-{\sigma }_{\ell }\right)y_{\ell }\left(t-{\sigma }_{\ell }\right)\right)\hfill \\ \hfill & +{\tilde{\eta }}_3\stackrel{2}{\sum _{\ell =1}}{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}{\lambda }_{\ell }{\vartheta }_{\ell }\left(y_{\ell }-y_{\ell }\left(t-{\varpi }_{\ell }\right)\right)+{\tilde{\eta }}_3\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}{y}_{\ell }\left(t-{\varpi }_{\ell }\right)-\phi {\tilde{\eta }}_{3}v.\hfill \end{array}$$

(42)

Collecting the terms of Equation (42), we obtain

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\tilde{\Xi }}_0}{dt}}& =-\stackrel{2}{\sum _{\ell =1}}{\tilde{\eta }}_{\ell }{\gamma }_{\ell }\frac{{\left(x_{\ell }-x_{\ell }^0\right)}^2}{x_{\ell }}+{\vartheta }_1{\vartheta }_2\phi \left({\lambda }_1{e}^{-({ϵ}_1{\sigma }_1+{\epsilon }_1{\varpi }_1)}(1-{\widehat{\psi }}_a)\frac{{\beta }_1{x}_1^0}{\phi }\right.\hfill \\ \hfill & \left.+{\lambda }_2{e}^{-({ϵ}_2{\sigma }_2+{\epsilon }_2{\varpi }_2)}(1-{\widehat{\varphi }}_a)\frac{{\beta }_2{x}_2^0}{\phi }-(1-{\widehat{\varphi }}_a)(1-{\widehat{\psi }}_a)\right)v.\hfill \end{array}$$

(43)

Simplifying (43), we get

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\tilde{\Xi }}_0}{dt}}& =-\stackrel{2}{\sum _{\ell =1}}{\tilde{\eta }}_{\ell }{\gamma }_{\ell }\frac{{\left(x_{\ell }-x_{\ell }^0\right)}^2}{x_{\ell }}+{\vartheta }_1{\vartheta }_2\phi \left(\widehat{\Phi }+\widehat{\Psi }-{\widehat{\varphi }}_b{\widehat{\psi }}_a-{\widehat{\varphi }}_a{\widehat{\psi }}_b-{\widehat{\varphi }}_a{\widehat{\psi }}_a-1\right)v\hfill \\ \hfill & =-\stackrel{2}{\sum _{\ell =1}}{\tilde{\eta }}_{\ell }{\gamma }_{\ell }\frac{{\left(x_{\ell }-x_{\ell }^0\right)}^2}{x_{\ell }}-{\vartheta }_1{\vartheta }_2\phi \left(1-\widehat{M}\right)v.\hfill \end{array}$$

As a result, $\frac{d{\tilde{\Xi }}_0}{dt}\le 0$ if ${\mathcal{R}}_0\le 1$ for $x_{\ell },y_{\ell },v\in (0,\infty ),\ell =1,2$. Moreover, $\frac{d{\tilde{\Xi }}_0}{dt}\le 0=0$ when $x_{\ell }\left(t\right)=x_{\ell }^0$ and $v\left(t\right)=0,\ell =1,2,$ for all t. The solutions of System (29)–(31) tend to ${\tilde{\Gamma }}_0^{{}^{\prime }}$, which has elements with $v\left(t\right)=0,$ so $\dot{v}\left(t\right)=0.$ Hence, from Equation (35), we get

$$0=\dot{v}\left(t\right)=\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}{y}_{\ell }(t-{\varpi }_{\ell })⟹y_{\ell }\left(t\right)=0,\ell =1,2.$$

Hence, ${\tilde{\Gamma }}_0^{{}^{\prime }}=$$\left\{{\Pi }_0\right\}$, and by applying LIP, we get that ${\Pi }_0$ is GAS. □

 Theorem 4.

For System (29)–(31), if ${\mathcal{R}}_0>1$, then $\tilde{\Pi }$is GAS.

 Proof.

Define ${\tilde{\Xi }}_1(x_1,y_1,x_2,y_2,v)$ as:

$$\begin{array}{cc}\hfill {\tilde{\Xi }}_1& =\stackrel{2}{\sum _{\ell =1}}{\overline{\eta }}_{\ell }\left[{\tilde{x}}_{\ell }\Theta \left(\frac{x_{\ell }}{{\tilde{x}}_{\ell }}\right)+e^{{ϵ}_{\ell }{\sigma }_{\ell }}{\tilde{y}}_{\ell }\Theta \left(\frac{y_{\ell }}{{\tilde{y}}_{\ell }}\right)\right]\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\overline{\eta }}_{\ell }{\beta }_{\ell }{\tilde{x}}_{\ell }\tilde{v}\underset{0}{\overset{{\sigma }_{\ell }}{\int }}\Theta \left(\frac{x_{\ell }\left(t-\xi \right)v\left(t-\xi \right)}{{\tilde{x}}_{\ell }\tilde{v}}\right)d\xi \hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\overline{\eta }}_{\ell }{\overline{\beta }}_{\ell }{\tilde{x}}_{\ell }{\tilde{y}}_{\ell }\underset{0}{\overset{{\sigma }_{\ell }}{\int }}\Theta \left(\frac{x_{\ell }\left(t-\xi \right)y_{\ell }\left(t-\xi \right)}{{\tilde{x}}_{\ell }{\tilde{y}}_{\ell }}\right)d\xi \hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{\tilde{y}}_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}\underset{0}{\overset{{\varpi }_{\ell }}{\int }}\Theta \left(\frac{y_{\ell }\left(t-\xi \right)}{{\tilde{y}}_{\ell }}\right)d\xi +\tilde{v}\Theta \left(\frac{v}{\tilde{v}}\right),\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill {\overline{\eta }}_1& =\frac{{\lambda }_1{e}^{-{\epsilon }_1{\varpi }_1}{\vartheta }_1{\tilde{y}}_1}{{\beta }_1{\tilde{x}}_1\tilde{v}}=\frac{{\lambda }_1{e}^{-\left({ϵ}_1{\sigma }_1+{\epsilon }_1{\varpi }_1\right)}\left({\beta }_1{\tilde{x}}_1\tilde{v}+{\overline{\beta }}_1{\tilde{x}}_1{\tilde{y}}_1\right)}{{\beta }_1{\tilde{x}}_1\tilde{v}},\hfill \\ \hfill {\overline{\eta }}_2& =\frac{{\lambda }_2{e}^{-{\epsilon }_2{\varpi }_2}{\vartheta }_2{\tilde{y}}_2}{{\beta }_2{\tilde{x}}_2\tilde{v}}=\frac{{\lambda }_2{e}^{-\left({ϵ}_2{\sigma }_2+{\epsilon }_2{\varpi }_2\right)}\left({\beta }_2{\tilde{x}}_2\tilde{v}+{\overline{\beta }}_2{\tilde{x}}_2{\tilde{y}}_2\right)}{{\beta }_2{\tilde{x}}_2\tilde{v}}.\hfill \end{array}$$

(44)

Clearly, ${\tilde{\Xi }}_1(x_1,y_1,x_2,y_2,v)>0$ for all $x_1,y_1,x_2,y_2,v>0$, and ${\tilde{\Xi }}_1({\tilde{x}}_1,{\tilde{y}}_1,{\tilde{x}}_2,{\tilde{y}}_2,\tilde{v})=0$. Calculating $\frac{d{\tilde{\Xi }}_1}{dt}$ along the trajectories of (29)–(31), we get

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\tilde{\Xi }}_1}{dt}}& =\stackrel{2}{\sum _{\ell =1}}{\overline{\eta }}_{\ell }\left(1-\frac{{\tilde{x}}_{\ell }}{x_{\ell }}\right)\left({\alpha }_{\ell }-{\gamma }_{\ell }{x}_{\ell }-{\beta }_{\ell }{x}_{\ell }v-{\overline{\beta }}_{\ell }{x}_{\ell }{y}_{\ell }\right)\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\overline{\eta }}_{\ell }{e}^{{ϵ}_{\ell }{\sigma }_{\ell }}\left(1-\frac{{\tilde{y}}_{\ell }}{y_{\ell }}\right)\left({\beta }_{\ell }{e}^{-{ϵ}_{\ell }{\sigma }_{\ell }}{x}_{\ell }(t-{\sigma }_{\ell })v(t-{\sigma }_{\ell })+{\overline{\beta }}_{\ell }{e}^{-{ϵ}_{\ell }{\sigma }_{\ell }}{x}_{\ell }(t-{\sigma }_{\ell })y_{\ell }(t-{\sigma }_{\ell })-{\vartheta }_{\ell }{y}_{\ell }\right)\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\overline{\eta }}_{\ell }{\beta }_{\ell }{\tilde{x}}_{\ell }\tilde{v}\left(\frac{x_{\ell }v}{{\tilde{x}}_{\ell }\tilde{v}}-\frac{x_{\ell }\left(t-{\sigma }_{\ell }\right)v\left(t-{\sigma }_{\ell }\right)}{{\tilde{x}}_{\ell }\tilde{v}}+ln\left(\frac{x_{\ell }\left(t-{\sigma }_{\ell }\right)v\left(t-{\sigma }_{\ell }\right)}{x_{\ell }v}\right)\right)\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\overline{\eta }}_{\ell }{\overline{\beta }}_{\ell }{\tilde{x}}_{\ell }{\tilde{y}}_{\ell }\left(\frac{x_{\ell }{y}_{\ell }}{{\tilde{x}}_{\ell }{\tilde{y}}_{\ell }}-\frac{x_{\ell }\left(t-{\sigma }_{\ell }\right)y_{\ell }\left(t-{\sigma }_{\ell }\right)}{{\tilde{x}}_{\ell }{\tilde{y}}_{\ell }}+ln\left(\frac{x_{\ell }\left(t-{\sigma }_{\ell }\right)y_{\ell }\left(t-{\sigma }_{\ell }\right)}{x_{\ell }{y}_{\ell }}\right)\right)\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{\tilde{y}}_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}\left(\frac{y_{\ell }}{{\tilde{y}}_{\ell }}-\frac{y_{\ell }\left(t-{\varpi }_{\ell }\right)}{{\tilde{y}}_{\ell }}+ln\left(\frac{y_{\ell }\left(t-{\varpi }_{\ell }\right)}{y_{\ell }}\right)\right)\hfill \\ \hfill & +\left(1-\frac{\tilde{v}}{v}\right)\left(\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}{y}_{\ell }\left(t-{\varpi }_{\ell }\right)-\phi v\right).\hfill \end{array}$$

Thus,

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\tilde{\Xi }}_1}{dt}}& =\stackrel{2}{\sum _{\ell =1}}{\overline{\eta }}_{\ell }\left(1-\frac{{\tilde{x}}_{\ell }}{x_{\ell }}\right)\left({\alpha }_{\ell }-{\gamma }_{\ell }{x}_{\ell }\right)\stackrel{2}{\underset{\ell =1}{-\sum }}{\overline{\eta }}_{\ell }\left({\beta }_{\ell }{x}_{\ell }v+{\overline{\beta }}_{\ell }{x}_{\ell }{y}_{\ell }-{\beta }_{\ell }{x}_{\ell }v\frac{{\tilde{x}}_{\ell }}{x_{\ell }}-{\overline{\beta }}_{\ell }{x}_{\ell }{y}_{\ell }\frac{{\tilde{x}}_{\ell }}{x_{\ell }}\right)\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\overline{\eta }}_{\ell }\left({\beta }_{\ell }{x}_{\ell }(t-{\sigma }_{\ell })v(t-{\sigma }_{\ell })+{\overline{\beta }}_{\ell }{x}_{\ell }(t-{\sigma }_{\ell })y_{\ell }(t-{\sigma }_{\ell })-e^{{ϵ}_{\ell }{\sigma }_{\ell }}{\vartheta }_{\ell }{y}_{\ell }\right)\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\overline{\eta }}_{\ell }\left(-{\beta }_{\ell }{x}_{\ell }(t-{\sigma }_{\ell })v(t-{\sigma }_{\ell })\frac{{\tilde{y}}_{\ell }}{y_{\ell }}-{\overline{\beta }}_{\ell }{x}_{\ell }(t-{\sigma }_{\ell })y_{\ell }(t-{\sigma }_{\ell })\frac{{\tilde{y}}_{\ell }}{y_{\ell }}+e^{{ϵ}_{\ell }{\sigma }_{\ell }}{\vartheta }_{\ell }{y}_{\ell }\frac{{\tilde{y}}_{\ell }}{y_{\ell }}\right)\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\overline{\eta }}_{\ell }\left({\beta }_{\ell }{x}_{\ell }v-{\beta }_{\ell }{x}_{\ell }\left(t-{\sigma }_{\ell }\right)v\left(t-{\sigma }_{\ell }\right)+{\beta }_{\ell }{\tilde{x}}_{\ell }\tilde{v}ln\left(\frac{x_{\ell }\left(t-{\sigma }_{\ell }\right)v\left(t-{\sigma }_{\ell }\right)}{x_{\ell }v}\right)\right)\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}\overline{\eta }\left({\overline{\beta }}_{\ell }{x}_{\ell }{y}_{\ell }-{\overline{\beta }}_{\ell }{x}_{\ell }\left(t-{\sigma }_{\ell }\right)y_{\ell }\left(t-{\sigma }_{\ell }\right)+{\overline{\beta }}_{\ell }{\tilde{x}}_{\ell }{\tilde{y}}_{\ell }ln\left(\frac{x_{\ell }\left(t-{\sigma }_{\ell }\right)y_{\ell }\left(t-{\sigma }_{\ell }\right)}{x_{\ell }{y}_{\ell }}\right)\right)\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}\left({\lambda }_{\ell }{\vartheta }_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}{y}_{\ell }-{\lambda }_{\ell }{\vartheta }_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}{y}_{\ell }\left(t-{\varpi }_{\ell }\right)+{\lambda }_{\ell }{\vartheta }_{\ell }{\tilde{y}}_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}ln\left(\frac{y_{\ell }\left(t-{\varpi }_{\ell }\right)}{y_{\ell }}\right)\right)\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}{y}_{\ell }\left(t-{\varpi }_{\ell }\right)-\phi v-\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}{y}_{\ell }\left(t-{\varpi }_{\ell }\right)\frac{\tilde{v}}{v}+\phi \tilde{v}.\hfill \end{array}$$

(45)

Simplifying Equation (45), we then have

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\tilde{\Xi }}_1}{dt}}& =\stackrel{2}{\sum _{\ell =1}}{\overline{\eta }}_{\ell }\left(1-\frac{{\tilde{x}}_{\ell }}{x_{\ell }}\right)\left({\alpha }_{\ell }-{\gamma }_{\ell }{x}_{\ell }\right)\stackrel{2}{\underset{\ell =1}{+\sum }}{\overline{\eta }}_{\ell }\left({\beta }_{\ell }{x}_{\ell }v\frac{{\tilde{x}}_{\ell }}{x_{\ell }}+{\overline{\beta }}_{\ell }{x}_{\ell }{y}_{\ell }\frac{{\tilde{x}}_{\ell }}{x_{\ell }}\right)\hfill \\ \hfill & -\stackrel{2}{\sum _{\ell =1}}{\overline{\eta }}_{\ell }{e}^{{ϵ}_{\ell }{\sigma }_{\ell }}{\vartheta }_{\ell }{y}_{\ell }+\stackrel{2}{\sum _{\ell =1}}{\overline{\eta }}_{\ell }\left(-{\beta }_{\ell }{x}_{\ell }(t-{\sigma }_{\ell })v(t-{\sigma }_{\ell })\frac{{\tilde{y}}_{\ell }}{y_{\ell }}-{\overline{\beta }}_{\ell }{x}_{\ell }(t-{\sigma }_{\ell })y_{\ell }(t-{\sigma }_{\ell })\frac{{\tilde{y}}_{\ell }}{y_{\ell }}+e^{{ϵ}_{\ell }{\sigma }_{\ell }}{\vartheta }_{\ell }{\tilde{y}}_{\ell }\right)\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\overline{\eta }}_{\ell }{\beta }_{\ell }{\tilde{x}}_{\ell }\tilde{v}ln\left(\frac{x_{\ell }\left(t-{\sigma }_{\ell }\right)v\left(t-{\sigma }_{\ell }\right)}{x_{\ell }v}\right)+\stackrel{2}{\sum _{\ell =1}}\overline{\eta }{\overline{\beta }}_{\ell }{\tilde{x}}_{\ell }{\tilde{y}}_{\ell }ln\left(\frac{x_{\ell }\left(t-{\sigma }_{\ell }\right)y_{\ell }\left(t-{\sigma }_{\ell }\right)}{x_{\ell }{y}_{\ell }}\right)\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}\left({\lambda }_{\ell }{\vartheta }_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}{y}_{\ell }-{\lambda }_{\ell }{\vartheta }_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}{y}_{\ell }\left(t-{\varpi }_{\ell }\right)+{\lambda }_{\ell }{\vartheta }_{\ell }{\tilde{y}}_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}ln\left(\frac{y_{\ell }\left(t-{\varpi }_{\ell }\right)}{y_{\ell }}\right)\right)\hfill \\ \hfill & +\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}{y}_{\ell }\left(t-{\varpi }_{\ell }\right)-\phi v-\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}{y}_{\ell }\left(t-{\varpi }_{\ell }\right)\frac{\tilde{v}}{v}+\phi \tilde{v}.\hfill \end{array}$$

(46)

Collecting the terms of Equation (46) and applying the equilibrium conditions for $\tilde{\Pi }$:

$$\begin{array}{cc}\hfill {\alpha }_{\ell }& ={\gamma }_{\ell }{\tilde{x}}_{\ell }+{\beta }_{\ell }{\tilde{x}}_{\ell }\tilde{v}+{\overline{\beta }}_{\ell }{\tilde{x}}_{\ell }{\tilde{y}}_{\ell },\ell =1,2,\hfill \\ \hfill {\vartheta }_{\ell }{\tilde{y}}_{\ell }& ={\beta }_{\ell }{e}^{-{ϵ}_{\ell }{\sigma }_{\ell }}{\tilde{x}}_{\ell }\tilde{v}+{\overline{\beta }}_{\ell }{e}^{-{ϵ}_{\ell }{\sigma }_{\ell }}{\tilde{x}}_{\ell }{\tilde{y}}_{\ell },\ell =1,2,\hfill \\ \hfill \phi \tilde{v}& =\stackrel{2}{\sum _{\ell =1}}{\lambda }_{\ell }{\vartheta }_{\ell }{e}^{-{\epsilon }_{\ell }{\varpi }_{\ell }}{\tilde{y}}_{\ell }=\stackrel{2}{\sum _{\ell =1}}{\overline{\eta }}_{\ell }{\beta }_{\ell }{\tilde{x}}_{\ell }\tilde{v},\hfill \end{array}$$

we get

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\tilde{\Xi }}_1}{dt}}& =\stackrel{2}{\sum _{\ell =1}}\left[-{\overline{\eta }}_{\ell }{\gamma }_{\ell }\frac{{\left(x_{\ell }-{\tilde{x}}_{\ell }\right)}^2}{x_{\ell }}+{\overline{\eta }}_{\ell }\left({\beta }_{\ell }{\tilde{x}}_{\ell }\tilde{v}+{\overline{\beta }}_{\ell }{\tilde{x}}_{\ell }{\tilde{y}}_{\ell }\right)\left(1-\frac{{\tilde{x}}_{\ell }}{x_{\ell }}\right)\right.\hfill \\ \hfill & -{\overline{\eta }}_{\ell }{\beta }_{\ell }{\tilde{x}}_{\ell }\tilde{v}\frac{x_{\ell }(t-{\sigma }_{\ell })v(t-{\sigma }_{\ell }){\tilde{y}}_{\ell }}{{\tilde{x}}_{\ell }\tilde{v}{y}_{\ell }}-{\overline{\eta }}_{\ell }{\overline{\beta }}_{\ell }{\tilde{x}}_{\ell }{\tilde{y}}_{\ell }\frac{x_{\ell }(t-{\sigma }_{\ell })y_{\ell }(t-{\sigma }_{\ell })}{{\tilde{x}}_{\ell }{y}_{\ell }}\hfill \\ \hfill & +{\overline{\eta }}_{\ell }{\beta }_{\ell }{\tilde{x}}_{\ell }\tilde{v}+{\overline{\eta }}_{\ell }{\overline{\beta }}_{\ell }{\tilde{x}}_{\ell }{\tilde{y}}_{\ell }+{\overline{\eta }}_{\ell }{\beta }_{\ell }{\tilde{x}}_{\ell }\tilde{v}ln\left(\frac{x_{\ell }\left(t-{\sigma }_{\ell }\right)v\left(t-{\sigma }_{\ell }\right)}{x_{\ell }v}\right)\hfill \\ \hfill & +{\overline{\eta }}_{\ell }{\overline{\beta }}_{\ell }{\tilde{x}}_{\ell }{\tilde{y}}_{\ell }ln\left(\frac{x_{\ell }\left(t-{\sigma }_{\ell }\right)y_{\ell }\left(t-{\sigma }_{\ell }\right)}{x_{\ell }{y}_{\ell }}\right)+{\overline{\eta }}_{\ell }{\beta }_{\ell }{\tilde{x}}_{\ell }\tilde{v}ln\left(\frac{y_{\ell }\left(t-{\varpi }_{\ell }\right)}{y_{\ell }}\right)\hfill \\ \hfill & \left.-{\overline{\eta }}_{\ell }{\beta }_{\ell }{\tilde{x}}_{\ell }\tilde{v}\frac{y_{\ell }\left(t-{\varpi }_{\ell }\right)\tilde{v}}{{\tilde{y}}_{\ell }v}\right]+\stackrel{2}{\sum _{\ell =1}}{\overline{\eta }}_{\ell }{\beta }_{\ell }{\tilde{x}}_{\ell }\tilde{v}.\hfill \end{array}$$

(47)

Using the following equalities

$$\begin{array}{cc}\hfill ln\left(\frac{x_{\ell }(t-{\sigma }_{\ell })v(t-{\sigma }_{\ell })}{x_{\ell }v}\right)& =ln\left(\frac{{\tilde{x}}_{\ell }}{x_{\ell }}\right)+ln\left(\frac{x_{\ell }(t-{\sigma }_{\ell })v(t-{\sigma }_{\ell }){\tilde{y}}_{\ell }}{y_{\ell }{\tilde{x}}_{\ell }\tilde{v}}\right)+ln\left(\frac{y_{\ell }\tilde{v}}{{\tilde{y}}_{\ell }v}\right),\hfill \\ \hfill ln\left(\frac{x_{\ell }(t-{\sigma }_{\ell })y_{\ell }(t-{\sigma }_{\ell })}{x_{\ell }{y}_{\ell }}\right)& =ln\left(\frac{{\tilde{x}}_{\ell }}{x_{\ell }}\right)+ln\left(\frac{x_{\ell }(t-{\sigma }_{\ell })y_{\ell }(t-{\sigma }_{\ell })}{{\tilde{x}}_{\ell }{y}_{\ell }}\right),\hfill \\ \hfill ln\left(\frac{y_{\ell }(t-{\varpi }_{\ell })}{y_{\ell }}\right)& =ln\left(\frac{y_{\ell }(t-{\varpi }_{\ell })\tilde{v}}{{\tilde{y}}_{\ell }v}\right)+ln\left(\frac{{\tilde{y}}_{\ell }v}{y_{\ell }\tilde{v}}\right),\ell =1,2,\hfill \end{array}$$

(48)

Equation (47) becomes

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\tilde{\Xi }}_1}{dt}}& =-\stackrel{2}{\sum _{\ell =1}}\left[{\overline{\eta }}_{\ell }{\gamma }_{\ell }\frac{{\left(x_{\ell }-{\tilde{x}}_{\ell }\right)}^2}{x_{\ell }}+{\overline{\eta }}_{\ell }\left({\beta }_{\ell }{\tilde{x}}_{\ell }\tilde{v}+{\overline{\beta }}_{\ell }{\tilde{x}}_{\ell }{\tilde{y}}_{\ell }\right)\Theta \left(\frac{{\tilde{x}}_{\ell }}{x_{\ell }}\right)\right.\hfill \\ \hfill & +{\overline{\eta }}_{\ell }{\beta }_{\ell }{\tilde{x}}_{\ell }\tilde{v}\left(\Theta \left(\frac{x_{\ell }(t-{\sigma }_{\ell })v(t-{\sigma }_{\ell }){\tilde{y}}_{\ell }}{{\tilde{x}}_{\ell }\tilde{v}{y}_{\ell }}\right)+\Theta \left(\frac{y_{\ell }\left(t-{\varpi }_{\ell }\right)\tilde{v}}{{\tilde{y}}_{\ell }v}\right)\right)\hfill \\ \hfill & \left.+{\overline{\eta }}_{\ell }{\overline{\beta }}_{\ell }{\tilde{x}}_{\ell }{\tilde{y}}_{\ell }\Theta \left(\frac{x_{\ell }(t-{\sigma }_{\ell })y_{\ell }(t-{\sigma }_{\ell })}{{\tilde{x}}_{\ell }{y}_{\ell }}\right)\right].\hfill \end{array}$$

Therefore, ${\displaystyle \frac{d{\Xi }_1}{dt}}\le 0$ for all $x_{\ell },y_{\ell }$,$v>0,\ell =1,2$. It can be seen that ${\displaystyle \frac{d{\tilde{\Xi }}_1}{dt}}=0$ if and only if $(x_1\left(t\right),y_1\left(t\right),x_2\left(t\right),y_2\left(t\right),v\left(t\right))=({\tilde{x}}_1,{\tilde{y}}_1,{\tilde{x}}_2,{\tilde{y}}_2,\tilde{v}).$ Therefore, the solutions of System (29)–(31) tend to ${\tilde{\Gamma }}_1^{{}^{\prime }}=\left\{\tilde{\Pi }\right\}$ and $\tilde{\Pi }$ is GAS when ${\mathcal{R}}_0>1$ according to LIP. □

4. Numerical Simulations

System (29)–(31) is analyzed numerically in this section, with a thorough discussion of numerical sensitivity analysis and illustrating the results of Theorems 3 and 4 numerically. Furthermore, we investigate the impact of temporal delays on the system’s dynamic behavior. To solve Systems (29)–(31) numerically, we rely on our calculations based on fixed parameters obtained from the modeling literature; see

Table 3. Parameters and their corresponding values (29)–(31).

Parameter	Value	References	Parameter	Value	References	Parameter	Value	Source
${\alpha }_1$	10	[44,45,46]	${\alpha }_2$	$0.03198$	[31]	${ϵ}_1$	$0.2$	[47]
${\gamma }_1$	$0.01$	[28,46,48]	${\gamma }_2$	$0.01$	[33,34]	${ϵ}_2$	1	Assumed
${\vartheta }_1$	$0.5$	[24,27,49]	${\vartheta }_2$	$0.1$	[27]	${\epsilon }_1$	1	[32]
${\lambda }_1$	6	[32]	${\lambda }_2$	6	[32]	${\epsilon }_2$	1	[32]
$\phi$	2	[46]

.

4.1. Sensitivity Analysis

Especially in pathology and epidemiology, sensitivity analysis plays a crucial role in modeling complicated interactions. By analyzing sensitivities, we can determine what parameters work best to curb the spread of the disease or the crime. Calculating the sensitivity indices can be done in three different ways: directly through direct differentiation, by a Latin hypercube sampling method, or by linearizing the model and solving the obtained linear algebraic equations.

In this study, direct differentiation is used since the indices can be expressed analytically. In cases where variables vary with respect to parameters, then the sensitivity index can be determined by partial derivatives [50]. As a function of a parameter, the normalized forward sensitivity index of ${\mathcal{R}}_0$ is expressed as follows:

$$S_{\varkappa }^{{\mathcal{R}}_0}=\frac{\varkappa }{{\mathcal{R}}_0}\frac{\partial {\mathcal{R}}_0}{\partial \varkappa },$$

(49)

where $\varkappa$ is a given parameter. The sensitivity indices for each parameter included in ${\mathcal{R}}_0$ are calculated using Equation (49). For instance, the sensitivity index of a parameter value with respect to ${\beta }_1$ is computed as

$$S_{{\beta }_1}^{{\mathcal{R}}_0}=\frac{{\beta }_1}{{\mathcal{R}}_0}\frac{\partial {\mathcal{R}}_0}{\partial {\beta }_1}=\frac{1}{2}\frac{1}{{\mathcal{R}}_0}\left({\widehat{\varphi }}_b+\frac{\left(\widehat{\Phi }-\widehat{\Psi }\right){\widehat{\varphi }}_b+2{\widehat{\varphi }}_b{\widehat{\psi }}_b}{\sqrt{{\left(\widehat{\Phi }-\widehat{\Psi }\right)}^2+4{\widehat{\varphi }}_b{\widehat{\psi }}_b}}\right).$$

The sensitivity index of ${\mathcal{R}}_0$ is shown in Figure 2 and

Table 4. Sensitivity index of ${\mathcal{R}}_0$.

Parameter	Sensitivity Index	Parameter	Sensitivity Index	Parameter	Sensitivity Index
${\alpha }_1$	$0.999$	${\alpha }_2$	$9.141\times {10}^{-6}$	${ϵ}_2$	$-6.399\times {10}^{-6}$
${\gamma }_1$	$-0.999$	${\gamma }_2$	$-9.141\times {10}^{-6}$	${\epsilon }_1$	$-1.373\times {10}^{-2}$
${\beta }_1$	$6.864\times {10}^{-2}$	${\beta }_2$	$9.140\times {10}^{-6}$	${\epsilon }_2$	$-1.828\times {10}^{-6}$
${\overline{\beta }}_1$	$9.313\times {10}^{-1}$	${\overline{\beta }}_2$	$1.651\times {10}^{-9}$	${\sigma }_1$	$-2.000\times {10}^{-1}$
${\vartheta }_1$	$-9.313\times {10}^{-1}$	${\vartheta }_2$	$-1.651\times {10}^{-9}$	${\sigma }_2$	$-6.399\times {10}^{-6}$
${\lambda }_1$	$6.864\times {10}^{-2}$	${\lambda }_2$	$9.140\times {10}^{-6}$	${\varpi }_1$	$-1.373\times {10}^{-2}$
$\phi$	$-6.865\times {10}^{-2}$	${ϵ}_1$	$-2.000\times {10}^{-1}$	${\varpi }_2$	$-1.828\times {10}^{-6}$

based on the parameter values in Table 3 and the values ${\sigma }_1=1,{\sigma }_2=0.7,{\varpi }_1={\varpi }_2=0.2,{\beta }_1=0.00003$, ${\beta }_2=0.00003$, ${\overline{\beta }}_1=0.0005$, and ${\overline{\beta }}_2=0.00001.$

Clearly, ${\alpha }_1$, ${\alpha }_2$, ${\beta }_1,{\beta }_2,{\overline{\beta }}_1,{\overline{\beta }}_2,{\lambda }_1,$ and ${\lambda }_2$ have positive indices. In terms of sensitivity, ${\alpha }_1$ is the most important parameter and ${\overline{\beta }}_2$ is the least important. In this case, there is a positive relationship between the persistence of the infection of the disease and the increase in the values of the parameters ${\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,{\overline{\beta }}_1,{\overline{\beta }}_2,{\lambda }_1,$ and ${\lambda }_2$ while keeping other parameters constant. The remaining indices are negative, i.e., the value of ${\mathcal{R}}_0$ decreases as ${\gamma }_1,{\gamma }_2,{\vartheta }_1,{\vartheta }_2$,$\phi ,{ϵ}_1,{ϵ}_2,{\epsilon }_1,{\epsilon }_2,{\sigma }_1,{\sigma }_2,{\varpi }_1,$ and ${\varpi }_2$ values increase.

4.2. Stability of the Equilibria

We analyze in this part the dynamic behavior of Model (29)–(31) numerically using MATLAB with the dde23 solver. In order to carry out the numerical simulation, we take the following considerations:

-

  The chosen delay parameters are ${\sigma }_1=1,{\sigma }_2=0.7,{\varpi }_1={\varpi }_2=0.2,$

-

  We pick three different initial conditions for System (29)–(31):

•     I.1: $(x_1\left(\xi \right),y_1\left(\xi \right),x_2\left(\xi \right),y_2\left(\xi \right),v\left(\xi \right))=(700,5,1.5,0.005,15)$;
•     I.2: $(x_1\left(\xi \right),y_1\left(\xi \right),x_2\left(\xi \right),y_2\left(\xi \right),v\left(\xi \right))=(500,3,1,0.02,12)$;
•     I.3: $(x_1\left(\xi \right),y_1\left(\xi \right),x_2\left(\xi \right),y_2\left(\xi \right),v\left(\xi \right))=(300,1.5,0.5,0.03,8)$, where $\xi \in [-max$$\{{\sigma }_1,{\sigma }_2,{\varpi }_1,{\varpi }_2\},0].$

By choosing different parameter values of the infection rates ${\beta }_1$, ${\beta }_2$, ${\overline{\beta }}_1$, and ${\overline{\beta }}_2$, we have the following outcomes:

Scenario 1 (Stability of $\Pi$${}_0$): We select ${\beta }_1=0.00003$, ${\beta }_2=0.0005$, ${\overline{\beta }}_1=0.00003$, ${\overline{\beta }}_2=0.00001$. This gives ${\mathcal{R}}_0=0.1105<1,$ and the solution of system (29)–(31) converges asymptotically to the IFE ${\Pi }_0=\left(1000,0,3.198,0,0\right).$ Figure 3 shows that the concentrations of both uninfected CD$4^{+}$T cells and uninfected macrophages increase and reach healthy values $x_1=1000,x_2=3.198$, while the concentrations of other compartments rapidly decline and reach zero. This affirms the global stability of ${\Pi }_0$, and the numerical results in this scenario coincide with the result of Theorem 3. In this case, HIV-1 is cleared, and the case simulates the state of a healthy person without HIV-1.

Scenario 2 (Stability of $\tilde{\Pi }$): We consider ${\beta }_1=0.0005$,${\beta }_2=0.005$, ${\overline{\beta }}_1=0.0002$, ${\overline{\beta }}_2=0.0001$ to obtain ${\mathcal{R}}_0=1.3478>1$. The solution of System (29)–(31) converges asymptotically to the IPE $\tilde{\Pi }=\left(747.095,4.14,0.899,0.114,5.110\right)$. In Figure 4, we can observe the existence and the stability of the equilibrium $\tilde{\Pi }$ that is proved in Lemma 7 and Theorem 4. This point represents the situation of an HIV-1 patient; that is, HIV-1 infection will persist. Accordingly, the infection rate is one of the main factors in disease control during HIV infection.

Effect of the time delay on the stability of the equilibria: The parameter ${\mathcal{R}}_0,$ provided by Equations (39) and (40), depends on the time delay parameters ${\sigma }_{\ell }$ and ${\varpi }_{\ell }$, $\ell =1,2$, which can change the stability of equilibria. We select ${\beta }_1=0.0005$, ${\beta }_2=0.005$, ${\overline{\beta }}_1=0.0002$, ${\overline{\beta }}_2=0.0001$ and ${\sigma }_1,{\sigma }_2,{\varpi }_1,$ and ${\varpi }_2$ will be varied to show the impact of time delay parameters. We examine the following cases:

• Case $\mathbf{I}$: ${\sigma }_1={\sigma }_2={\varpi }_1=$ ${\varpi }_2=0$;
• Case $\mathbf{II}$: ${\sigma }_1=0.3,{\sigma }_2=0.5,{\varpi }_1=0.4,$ ${\varpi }_2=0.2$;
• Case $\mathbf{III}$: ${\sigma }_1=0.8,{\sigma }_2=0.1,{\varpi }_1=0.3,$ ${\varpi }_2=1.3$;
• Case $\mathbf{IV}$: ${\sigma }_1=1.6,{\sigma }_2=1.2,{\varpi }_1=1.4,$ ${\varpi }_2=1.1.$

In Case $\mathbf{I}$, ${\mathcal{R}}_0$ becomes the form given by Equations (15),(16), and System (29)–(31) is reduced to System (6)–(31). With the above values, we solve Model (29)–(31) with the initial conditions I.1. Furthermore, the basic number ${\mathcal{R}}_0$ is calculated, and the values are equal to $\{1.9381,1.3408,1.2965,0.5614\}$ for each case, respectively. The impact of the time delay on the solution of our system is shown in Figure 5; it is apparent that when time delays were included, the number of uninfected cells of both CD$4^{+}$T cells and macrophages increased, while the number of other categories decreased.

Without loss of generality, suppose $\tau ={\sigma }_1={\varpi }_1$$={\sigma }_2={\varpi }_2$; then ${\mathcal{R}}_0$ can be expressed as follows:

$${\mathcal{R}}_0\left(\tau \right)=\frac{1}{2}\left(\tilde{\Phi }+\tilde{\Psi }+\sqrt{{\left(\tilde{\Phi }-\tilde{\Psi }\right)}^2+4{\tilde{\varphi }}_b{\tilde{\psi }}_b}\right),$$

(50)

where

$$\begin{array}{cc}\hfill \tilde{\Phi }& ={\tilde{\varphi }}_a+{\tilde{\varphi }}_b,\tilde{\Psi }={\tilde{\psi }}_a+{\tilde{\psi }}_b,\hfill \\ \hfill {\tilde{\varphi }}_a& =\frac{{\overline{\beta }}_1{e}^{-{ϵ}_1\tau }{x}_1^0}{{\vartheta }_1},{\tilde{\varphi }}_b=\frac{{\beta }_1{e}^{-({ϵ}_1+{\epsilon }_1)\tau }{x}_1^0{\lambda }_1}{\phi },{\tilde{\psi }}_a=\frac{{\overline{\beta }}_2{e}^{-{ϵ}_2\tau }{x}_2^0}{{\vartheta }_2},{\tilde{\psi }}_b=\frac{{\beta }_2{e}^{-({ϵ}_2+{\epsilon }_2)\tau }{x}_2^0{\lambda }_2}{\phi }.\hfill \end{array}$$

(51)

Using the values of the parameters given in Table 3, we obtain the following:

(i)

   If $\tau \ge 0.7035$, then ${\mathcal{R}}_0\left(\tau \right)\le 1$, and the IFE ${\Pi }_0$ is GAS;

(ii)

  If $\tau <0.7035$, then ${\mathcal{R}}_0\left(\tau \right)>1$, and ${\Pi }_0$ will become unstable.

The values of ${\mathcal{R}}_0$ and the stability case of ${\Pi }_0$ are shown in

Table 5. Values of ${\mathcal{R}}_0\left(\tau \right)$ for System (29)–(31) with different values of the time delay $\tau$.

$\mathit{\tau }$	${\mathcal{R}}_0\left(\mathit{\tau }\right)$	Stability of ${\mathbf{\Pi }}_0$
$0.0$	$1.9381$	unstable
$0.15$	$1.6684$	unstable
$0.3$	$1.4427$	unstable
$0.45$	$1.2535$	unstable
$0.6$	$1.0947$	unstable
$0.7035$	1	stable
$0.85$	$0.8838$	stable
$1.0$	$0.7831$	stable
$1.15$	$0.6978$	stable
$2.0$	$0.4045$	stable
$5.0$	$0.1509$	stable

for selected values of time delay. Clearly, with increasing time delays, ${\mathcal{R}}_0$ decreases, which accords with the findings of Section 4.1. As a result, the number of infected cells increases with decreasing time delay and vice versa. For more investigation, with a small time delay, the contact between virus or infected cells with healthy cells is faster and produces more infected cells.

5. Conclusions

In this paper, two HIV-1 infection models with two types of target cells CD$4^{+}$T cells and macrophages cells and two modes of transmission VTC and CTC were considered. The second model is a modification of the first one but incorporates four time delays. We proved that the two proposed models’ solutions are non-negative and bounded. We proved that each model has two possible equilibrium points: infection-free equilibrium (IFE), which is always present, and infection-present equilibrium (IPE), which is present if the basic reproduction number ${\mathcal{R}}_0>1$. The basic number ${\mathcal{R}}_0$ governs the dynamic behavior of the model: if ${\mathcal{R}}_0\le 1$, then the IFE point is GAS, and if ${\mathcal{R}}_0>1$, then the IPE point is GAS. In order to verify the theoretical results and investigate how the time delay affects the model solutions and the system’s dynamic behavior, we conducted some numerical computations. From those computations, it was noticed that:

(i)

  The trajectory diagrams tend towards the IFE when the reproduction number ${\mathcal{R}}_0\le 1$, as shown in Figure 3. One significant finding from these figures is that for different initial conditions assumed for the model categories, their trajectories still point towards the IFE over the passage of time. These findings also confirm the global asymptotic stability analysis results, which were presented in Section 3.3.

(ii)

  The trajectory diagrams tend towards the IPE for different initial conditions when the reproduction number ${\mathcal{R}}_0>1$, as shown in Figure 4, which confirms that the point IPE is GAS when ${\mathcal{R}}_0>1$. Consequently, the model leads to an outcome in which the person is infected with HIV-1.

(iii)

  From Figure 5 and Table 5, increasing the time delay causes a decrease in the reproduction number, resulting in an increase in uninfected CD$4^{+}$ T cells, resulting in a decrease in viral load. That is, time delay contributes a very significant effect in governing the dynamic behavior of the system and should not be neglected in HIV-1 modeling.

We also studied the sensitivity analysis to show how the values of all the parameters of the suggested model affect ${\mathcal{R}}_0$ for given data, and we saw that the persistence of the disease and the increase in the values of the parameters ${\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,{\overline{\beta }}_1,{\overline{\beta }}_2,{\lambda }_1,$ and ${\lambda }_2$ are positively correlated.

We assumed in this paper that: (i) the growth rate of uninfected cells is given in the form $\mathcal{U}\left(x\right)=\alpha$ without proliferation of uninfected cells. In many recent studies, it has been suggested that the growth rate of uninfected cells can take several forms such as: (a) growth rate with simple proliferation $\mathcal{U}\left(x\right)=\alpha +sx\left(1-\frac{x}{x_{max}}\right),$ where s is the rate of growth and $x_{max}$ is the maximum capacity of uninfected cells in the human body [27,51,52], and (b) growth rate with a full proliferation $\mathcal{U}(x,y)=\alpha +sx\left(1-\frac{x+y}{x_{max}}\right)$ [53]. On the other hand, it has been reported in [44] that HIV may be able to infect cells in the thymus and bone marrow and thus lead to reduced production of new uninfected cells. In this case, the production of uninfected cells $\alpha$ is a decreasing function of the viral load as $\alpha \left(v\right)=\frac{k\alpha }{k+v}$, where k is a constant [44]. These forms add some difficulties to the analysis of the model; therefore, we leave them to future work. We mention that our model can be extended in different directions by (i) considering the mutations of HIV, (ii) including the stochastic interaction, and (iii) considering the diffusion of cells and viruses.