# GLOBAL STABILITY OF VIRUS DYNAMICS MODEL WITH IMMUNE RESPONSE, CELLULAR INFECTION AND HOLLING TYPE-II

• ELAIW, A.M. (DEPARTMENT OF MATHEMATICS, KING ABDULAZIZ UNIVERSITY) ;
• GHALEB, SH.A. (DEPARTMENT OF MATHEMATICS, KING ABDULAZIZ UNIVERSITY)
• Accepted : 2019.03.11
• Published : 2019.03.25

#### Abstract

In this paper, we study the effect of Cytotoxic T Lymphocyte (CTL) and antibody immune responses on the virus dynamics with both virus-to-cell and cell-to-cell transmissions. The infection rate is given by Holling type-II. We first show that the model is biologically acceptable by showing that the solutions of the model are nonnegative and bounded. We find the equilibria of the model and investigate their global stability analysis. We derive five threshold parameters which fully determine the existence and stability of the five equilibria of the model. The global stability of all equilibria of the model is proven using Lyapunov method and applying LaSalle's invariance principle. To support our theoretical results we have performed some numerical simulations for the model. The results show the CTL and antibody immune response can control the disease progression.

# 1. INTRODUCTION

Mathematical modeling and analysis of within-host pathogen dynamics have attracted the interest of several mathematicians during the recent decades. The basic pathogen infection model has been proposed in [1] which describes the interaction between susceptible host cells, infected cells and pathogens and has been used to describe the dynamics of some types of viruses such as human immunodeficiency virus (HIV) and hepatitis B virus HBV and is given by

$\dot{s}$ = β − ˆδs − αsp,       (1.1)

$\dot{y}$ = αsp − εy,       (1.2)

$\dot{p}$ = my − γp,       (1.3)

where s, y and p denote the concentrations of susceptible cells, infected cells and pathogens, respectively. The susceptible (uninfected) cells are generated at rate β, die at rate $\hat{\delta}$s and become infected at rate αsp. The infected cells die at rate εy. Parameters m and γ represent, respectively, the generation and clearance rate constants of pathogens. Several modifications of model (1.1-1.3) have been done to take into account either Cytotoxic T Lymphocyte (CTL) immune response (see e.g. [2]-[7]) or humoral immune response (see e.g. [8]-[21]). Wodarz [22] has presented the following mathematical model to incorporate both humoral and CTL immunity into the pathogen dynamics:

$\dot{s}$ = β − ˆδs − αsp,       (1.4)

$\dot{y}$ = αsp − kwy − εy,       (1.5)

$\dot{p}$ = my − qzp − γp,       (1.6)

$\dot{z}$ = rzp − µz,       (1.7)

$\dot{w}$ = gwy − hw,       (1.8)

where, w and z denote the concentrations of CTL cells and B cells, respectively. Model (1.4)- (1.8) has been modified in several works (see e.g. [23]-[29]). In these works, only the pathogento-cell transmission has been considered. Mathematical models of pathogen dynamics with both pathogen-to-cell and cell-to-cell transmissions have been studied in several works (see e.g. [30]-[38]). However, both CTL and humoral immune responses have not been taken into account in these works.

The aim of this paper is to propose a pathogen infection model with both humoral and CTL immune responses. We have incorporated both pathogen-to-cell and cell-to-cell transmissions. The infection rate is given by Holling type-II. We first show that the solutions of the model are nonnegative and bounded. To investigate the global stability of the equilibria we construct Lyapunov functions using the method presented [39] and followed by [40]-[52]. Numerical simulations is performed to confirm our theoretical results.

# 2. MODEL WITH HOLLING-TYPE INCIDENCE

In this section, we propose the following pathogen dynamics model with Holling type-II incidence:

$\dot{s}=\beta-\hat{\delta} s-\frac{\alpha_{1} s p}{1+\eta s}-\frac{\alpha_{2} s y}{1+\eta s},$       (2.1)

$\dot{y}=\frac{\alpha_{1} s p}{1+\eta s}+\frac{\alpha_{2} s y}{1+\eta s}-\varepsilon y-k w y,$       (2.2)

$\dot{p}$ = my − γp − qzp,       (2.3)

$\dot{z}$ = rzp − µz,       (2.4)

$\dot{w}$ = gwy − hw,       (2.5)

where, η is positive constant, α1 and α2 are infection rates of pathogen-to-cell and cell-to-cell transmissions, respectively. The other variables are parameters have the same meaning as given in the previous section.

Proposition 2.1. There exist positive numbers Li , i = 1, 2, 3, 4 such that the compact set

Ω = { (s, y, p, z, w) ∈ $\mathbb{R}_{\geq 0}^{5}$ : 0 ≤ s, , y ≤ L1 , 0 ≤ p ≤ L2 , 0 ≤ z ≤ L3 , 0 ≤ w ≤ L4 }

is positively invariant.

Proof. From Eqs. (2.1)-(2.5) we have

$\dot{s}$ |s=0= β > 0,

$\dot{y}$  |y=0= $\frac{\alpha_{1} s p}{1+\eta s}$≥ 0,     for all s, p ≥ 0,

$\dot{p}$ |p=0= my ≥ 0,     for all y ≥ 0,

$\dot{z}$ |z=0= 0,

$\dot{w}$ |w=0= 0.

Hence, s(t) > 0, y(t) ≥ 0, p(t) ≥ 0, z(t) ≥ 0 and w(t) ≥ 0 for all t ≥ 0. Therefore model (2.1)-(2.5) is biologically acceptable in the sense that no population goes negative.

Next we show the boundedness of the solutions. Let Q(t) = s(t) + y(t) + $\frac{\varepsilon}{2 m}$p(t) +$\frac{\varepsilon q}{2 r m}$z(t) +$\frac{k}{g}$w(t), then

$\dot{Q}=\dot{s}+\dot{y}+\frac{\varepsilon}{2 m} \dot{p}+\frac{\varepsilon q}{2 r m} \dot{z}+\frac{k}{g} \dot{w}$

$=\beta-\hat{\delta} s-\frac{\varepsilon}{2} y-\frac{\varepsilon}{2 m} \gamma p-\frac{\varepsilon q}{2 r m} \mu z-\frac{k}{g} h w$

$\leq \beta-\sigma\left(s+y+\frac{\varepsilon}{2 m} p+\frac{\varepsilon q}{2 r m} z+\frac{k}{g} w\right)=\beta-\sigma Q,$

where σ = min{ $\hat{\delta}$,$\frac{\varepsilon}{2}$, γ, µ, h}. Then

$Q(t) \leq e^{-\sigma t}\left(Q(0)-\frac{\beta}{\sigma}\right)+\frac{\beta}{\sigma}$

Hence, 0 ≤ Q(t) ≤ L1, where L1 =$\frac{β}{σ}$. It follows that, 0 ≤ s(t), y(t) ≤ L1, 0 ≤ p(t) ≤ L2, 0≤ z(t) ≤ L3 and 0 ≤ w(t) ≤ L4 for all t ≥ 0 if s(0) + y(0) +$\frac{ε}{2m}$p(0) +$\frac{εq}{2rm}$z(0) +$\frac{k}{g}$w(0) ≤ L1, where L2 = $\frac{2mL_1}{ε}$, L3 =$\frac{2mrL_1}{εq}$ and L4 = $\frac{gL_1}{k}$. Therefore, s(t), y(t), p(t), z(t) and w(t) are bounded.

Lemma 2.2. For system (2.1)-(2.5) there exist five threshold parameters R0 > 0, $R_{1}^{z}$> 0, $R_{1}^{w}$> 0, $R_{2}^{w}$ > 0 and $R_{2}^{z}$ > 0, with $R_{1}^{w}$ < R0 such that

(i) if R0 ≤ 1, then there exists only one equilibrium Π0,

(ii) if $R_{1}^{z}$ ≤ 1, and $R_{1}^{w}$ ≤ 1 < R0, then there exist only two equilibria Π0 and Π1,

(iii) if $R_{1}^{z}$ > 1 and $R_{2}^{w}$ ≤ 1, then there exist only three equilibria Π0, Π1 and Π2,

(iv) if $R_{1}^{w}$ > 1 and $R_{2}^{z}$ ≤ 1, then there exist only three equilibria Π0, Π1 and Π3, and

(v) if $R_{2}^{w}$ > 1 and $R_{2}^{z}$ > 1, then there exist five equilibria Π0, Π1, Π2, Π3 and Π4.

Proof. The equilibria of system (2.6)-(2.10) satisfying

$\beta-\hat{\delta} s-\frac{\alpha_{1} s p}{1+\eta s}-\frac{\alpha_{2} s y}{1+\eta s}$ = 0,       (2.6)

$\frac{\alpha_{1} s p}{1+\eta s}+\frac{\alpha_{2} s y}{1+\eta s}$ − εy − kwy = 0,       (2.7)

my − γp − qzp = 0,       (2.8)

(rp − µ)z = 0,       (2.9)

(gy − h) w = 0.       (2.10)

We find that system (2.6)-(2.10) admits five equilibria.

(i) Infection-free equilibrium Π0 = (s0, 0, 0, 0, 0), where s0 = β/$\hat{\delta}$.

(ii) Chronic-infection equilibrium without immune response Π1 = (s1, y1, p1, 0, 0), where

s1 = $\frac{\beta}{(\beta \eta+\hat{\delta})\left(R_{0}-1\right)+\hat{\delta}},$

y1 = $\frac{\beta(\beta \eta+\hat{\delta})}{\varepsilon\left((\beta \eta+\hat{\delta})\left(R_{0}-1\right)+\hat{\delta}\right)}\left(R_{0}-1\right),$

p1 = $\frac{m \beta(\beta \eta+\hat{\delta})}{\varepsilon \gamma\left((\beta \eta+\hat{\delta})\left(R_{0}-1\right)+\hat{\delta}\right)}\left(R_{0}-1\right),$

and

R0 =$\frac{\beta\left(m \alpha_{1}+\gamma \alpha_{2}\right)}{\varepsilon \gamma(\beta \eta+\hat{\delta})}.$

Clearly Π1 exists if R0 > 1.

(iii) Chronic-infection equilibrium with only humoral immune response Π2 = (s2, y2, p2, z2, 0). Now we show that s2, y2, p2 and z2 which satisfy Eqs. (2.6)-(2.9) are positive. From Eq. (2.9) we have p2 = $µ\over r$. From Eqs. (2.6)-(2.7) we have

$\beta-\hat{\delta} s=\varepsilon y \Rightarrow y=\frac{\beta-\hat{\delta} s}{\varepsilon}.$       (2.11)

Substitute Eq. (2.11) in Eq. (2.6) and define a function G(s) as

$G(s)=\beta-\hat{\delta} s-\frac{\alpha_{1} s p_{2}}{1+\eta s}-\frac{\alpha_{2} s}{1+\eta s}(\beta-\hat{\delta} s)=0$       (2.12)

we have

G(0) = β > 0,

G(s0) = −$\frac{\alpha_{1} s_{0} p_{2}}{1+\eta s_{0}}$< 0.

Then, there exists s2 ∈ (0, s0) such as G(s2) = 0. Now

y = $\frac{\beta-\hat{\delta} s}{\varepsilon}=\frac{\hat{\delta}\left(\frac{\beta}{\delta}-s\right)}{\varepsilon}=\frac{\hat{\delta}}{\varepsilon}\left(s_{0}-s\right)$

⇒ y2 = $\frac{\hat{\delta}}{\varepsilon}$(s0 − s2) > 0.

And from Eq. (2.8) we have z2 = $\frac{\gamma}{q}\left(\frac{m y_{2}}{p_{2} \gamma}-1\right)$. Now we define

$R_{1}^{z}=\frac{m y_{2}}{p_{2} \gamma}=\frac{r m y_{2}}{\mu \gamma}.$       (2.13)

Then, z2 = $\frac{\gamma}{q}$($R_{1}^{z}$ -1) > 0 when $R_{1}^{z}$ > 1.

(iv) Chronic-infection equilibrium with only CTL immune response Π3 = (s3, y3, p3, 0, w3), where

$s_{3}=\frac{g \beta \gamma}{h\left(m \alpha_{1}+\gamma \alpha_{2}\right)+g \gamma \hat{\delta}}, \quad y_{3}=\frac{h}{g}$,

$p_{3}=\frac{h m}{g \gamma}, \quad w_{3}=\frac{\varepsilon}{k}\left(R_{1}^{w}-1\right)$,

and

$R_{1}^{w}=\frac{g \beta\left(m \alpha_{1}+\gamma \alpha_{2}\right)}{\varepsilon\left(h\left(m \alpha_{1}+\gamma \alpha_{2}\right)+g \gamma(\beta \eta+\hat{\delta})\right)}=\frac{R_{0}}{1+\frac{\varepsilon h}{g \beta} R_{0}}$

Clearly $R_{1}^{w}$< R0.

Hence, Π3 exists when $R_{1}^{w}$> 1.

(v) Chronic-infection equilibrium with both humoral and CTL immune responses Π4 = (s4, y4, p4, z4, w4), where

s4$\frac{\beta r g}{g \mu \alpha_{1}+h r \alpha_{2}+r g \hat{\delta}}, \quad y_{4}=\frac{h}{g}, \quad p_{4}=\frac{\mu}{r}$,

w4$\frac{\varepsilon}{k}\left(R_{2}^{w}-1\right), \quad z_{4}=\frac{\gamma}{q}\left(R_{2}^{z}-1\right)$,

and

$R_{2}^{w}=\frac{\beta g\left(\mu g \alpha_{1}+r h \alpha_{2}\right)}{\varepsilon h\left(g \mu \alpha_{1}+h r \alpha_{2}+r g(\beta \eta+\hat{\delta})\right)} \text { and } R_{2}^{z}=\frac{r}{\mu} p_{3}=\frac{m h r}{\gamma g \mu}.$

Hence, Π4 exists when $R_{2}^{w}$ > 1 and $R_{2}^{z}$ > 1.

# 3. GLOBAL STABILITY

We will use the arithmetic mean- geometric mean (AM-GM) inequality [38], through the paper. If θj ≥ 0, j = 1, 2, ..., n, it follows that

$\frac{1}{n} \sum_{j=1}^{n} \theta_{j} \geq \sqrt[n]{\prod_{j=1}^{n} \theta_{j}},$       (3.1)

with equality holding if and only if θ1 = θ2 = ... = θn. Define the function

H(ℓ) = ℓ − 1 − ln ℓ.       (3.2)

Clearly H(ℓ) ≥ 0 for any ℓ > 0 and H has the global minimum H(1) = 0.

Theorem 3.1. Let R0 ≤ 1, then Π0 is globally asymptotically stable (GAS).

Proof. Consider a Lyapunov functional

$L_{0}(s, y, p, z, w)=s-s_{0}-\int_{s_{0}}^{s} \frac{s_{0}(1+\eta \theta)}{\theta\left(1+\eta s_{0}\right)} d \theta+y+\frac{\alpha_{1} s_{0}}{\gamma\left(1+\eta s_{0}\right)} p+\frac{q \alpha_{1} s_{0}}{r \gamma\left(1+\eta s_{0}\right)} z+\frac{k}{g} w.$

We note L0(s, y, p, z, w) > 0 for all s, y, p, z, w > 0 and L0(s0, 0, 0, 0, 0) = 0. We calculate $\frac{d L_{0}}{d t}$ along the solutions of model (2.1)-(2.5) as:

$\frac{d L_{0}}{d t}=\left(1-\frac{s_{0}(1+\eta s)}{s\left(1+\eta s_{0}\right)}\right)\left(\beta-\hat{\delta} s-\frac{\alpha_{1} s p}{1+\eta s}-\frac{\alpha_{2} s y}{1+\eta s}\right)$

$+\frac{\alpha_{1} s p}{1+\eta s}+\frac{\alpha_{2} s y}{1+\eta s}-\varepsilon y-k w y+\frac{\alpha_{1} s_{0}}{\gamma\left(1+\eta s_{0}\right)}(m y-\gamma p-q z p)$

$+\frac{q \alpha_{1} s_{0}}{r \gamma\left(1+\eta s_{0}\right)}(r z p-\mu z)+\frac{k}{g}(g w y-h w).$       (3.3)

Collecting terms of Eq. (3.3) we get

\begin{aligned} \frac{d L_{0}}{d t} &=\left(1-\frac{s_{0}(1+\eta s)}{s\left(1+\eta s_{0}\right)}\right)(\beta-\hat{\delta} s)+\frac{\alpha_{2} s_{0} y}{1+\eta s_{0}}-\varepsilon y \\ &+\frac{\alpha_{1} s_{0}}{\gamma\left(1+\eta s_{0}\right)} m y-\frac{q \alpha_{1} s_{0}}{r \gamma\left(1+\eta s_{0}\right)} \mu z-\frac{k h}{g} w. \end{aligned}

Using the condition β = $\hat{\delta}$s0, we obtain

$\frac{d L_{0}}{d t}=\left(1-\frac{s_{0}(1+\eta s)}{s\left(1+\eta s_{0}\right)}\right)\left(\hat{\delta} s_{0}-\hat{\delta} s\right)+\frac{\alpha_{2} s_{0} y}{1+\eta s_{0}}-\varepsilon y$

$+\frac{\alpha_{1} s_{0}}{\gamma\left(1+\eta s_{0}\right)} m y-\frac{q \alpha_{1} s_{0}}{r \gamma\left(1+\eta s_{0}\right)} \mu z-\frac{k h}{g} w$

$=-\hat{\delta} \frac{\left(s-s_{0}\right)^{2}}{s\left(1+\eta s_{0}\right)}+\varepsilon\left(R_{0}-1\right) y-\frac{q \alpha_{1} s_{0}}{r \gamma\left(1+\eta s_{0}\right)} \mu z-\frac{k h}{g} w.$       (3.4)

Since R0 ≤ 1, then $\frac{d L_{0}}{d t}$ ≤ 0 for all s, y, z, w > 0. Moreover, $\frac{d L_{0}}{d t}$ = 0 when s = s0, y = 0, z = 0 and w = 0. Applying LaSalle’s invariance principle (LIP) we get that Π0 is GAS.

Theorem 3.2. Let $R_{1}^{z}$ ≤ 1 and $R_{1}^{w}$ ≤ 1 < R0, then Π1 is GAS.

Proof. Let us define a function L1(s, y, p, z, w) as:

L1 = $s-s_{1}-\int_{\mathrm{s} 1}^{s} \frac{s_{1}(1+\eta \theta)}{\theta\left(1+\eta s_{1}\right)} d \theta+y_{1} H\left(\frac{y}{y_{1}}\right)$       (3.5)

$+\frac{\alpha_{1} s_{1} p_{1}}{m y_{1}\left(1+\eta s_{1}\right)} p_{1} H\left(\frac{p}{p_{1}}\right)+\frac{q \alpha_{1} s_{1} p_{1}}{r m y_{1}\left(1+\eta s_{1}\right)} z+\frac{k}{g} w.$

Clearly, L1(s, y, p, z, w) > 0 for all s, y, p, z, w > 0 and L1(s1, y1, p1, 0, 0) = 0. Calculating $\frac{d L_{1}}{d t}$ along the trajectories of system (2.1)-(2.5), we obtain

\begin{aligned} \frac{d L_{1}}{d t} &=\left(1-\frac{s_{1}(1+\eta s)}{s\left(1+\eta s_{1}\right)}\right)\left(\beta-\hat{\delta} s-\frac{\alpha_{1} s p}{1+\eta s}-\frac{\alpha_{2} s y}{1+\eta s}\right) \\ &+\left(1-\frac{y_{1}}{y}\right)\left(\frac{\alpha_{1} s p}{1+\eta s}+\frac{\alpha_{2} s y}{1+\eta s}-\varepsilon y-k w y\right) \end{aligned}

$+\frac{\alpha_{1} s_{1} p_{1}}{m y_{1}\left(1+\eta s_{1}\right)}\left(1-\frac{p_{1}}{p}\right)(m y-\gamma p-q z p)$

$+\frac{q \alpha_{1} s_{1} p_{1}}{r m y_{1}\left(1+\eta s_{1}\right)}(r z p-\mu z)+\frac{k}{g}(g w y-h w)$       (3.6)

Collecting terms of Eq. (3.6) we get

\begin{aligned} \frac{d L_{1}}{d t} &=\left(1-\frac{s_{1}(1+\eta s)}{s\left(1+\eta s_{1}\right)}\right)(\beta-\hat{\delta} s)+\frac{\alpha_{1} s_{1} p}{1+\eta s_{1}}+\frac{\alpha_{2} s_{1} y}{1+\eta s_{1}}-\varepsilon y \\ &-\left(\frac{\alpha_{1} s p}{1+\eta s}+\frac{\alpha_{2} s y}{1+\eta s}\right) \frac{y_{1}}{y}+\varepsilon y_{1}+\frac{\alpha_{1} s_{1} p_{1}}{1+\eta s_{1}} \frac{y}{y_{1}}-\frac{\alpha_{1} s_{1} p_{1}}{m y_{1}\left(1+\eta s_{1}\right)} \gamma p \end{aligned}

$\begin{array}{l} -\frac{\alpha_{1} s_{1} p_{1}}{\left(1+\eta s_{1}\right)} \frac{y p_{1}}{y_{1} p}+\frac{\alpha_{1} s_{1} p_{1}}{m y_{1}\left(1+\eta s_{1}\right)} \gamma p_{1}+\frac{q \alpha_{1} s_{1} p_{1}}{m y_{1}\left(1+\eta s_{1}\right)} p_{1} z \\ -\frac{q \alpha_{1} s_{1} p_{1}}{m y_{1}\left(1+\eta s_{1}\right)} \frac{\mu}{r} z+k y_{1} w-\frac{k h}{g} w. \end{array}$

$=\left(1-\frac{s_{1}(1+\eta s)}{s\left(1+\eta s_{1}\right)}\right)(\beta-\hat{\delta} s)+\frac{\alpha_{1} s_{1} p_{1}}{1+\eta s_{1}} \frac{p}{p_{1}}+\frac{\alpha_{2} s_{1} y_{1}}{1+\eta s_{1}} \frac{y}{y_{1}}-\varepsilon y_{1} \frac{y}{y_{1}}$

$-\frac{\alpha_{1} s_{1} p_{1}}{1+\eta s_{1}} \frac{s p y_{1}\left(1+\eta s_{1}\right)}{s_{1} p_{1} y(1+\eta s)}-\frac{\alpha_{2} s_{1} y_{1}}{1+\eta s_{1}} \frac{s\left(1+\eta s_{1}\right)}{s_{1}(1+\eta s)}+\varepsilon y_{1}+\frac{\alpha_{1} s_{1} p_{1}}{1+\eta s_{1}} \frac{y}{y_{1}}$

$-\frac{\alpha_{1} s_{1} p_{1}}{m y_{1}\left(1+\eta s_{1}\right)} \gamma p_{1} \frac{p}{p_{1}}-\frac{\alpha_{1} s_{1} p_{1}}{\left(1+\eta s_{1}\right)} \frac{y p_{1}}{y_{1} p}+\frac{\alpha_{1} s_{1} p_{1}}{m y_{1}\left(1+\eta s_{1}\right)} \gamma p_{1}$

$+\frac{q \alpha_{1} s_{1} p_{1}}{m y_{1}\left(1+\eta s_{1}\right)}\left(p_{1}-\frac{\mu}{r}\right) z+k\left(y_{1}-\frac{h}{g}\right) w.$

Applying the equilibrium conditions for Π:

β = $\hat{\delta} s_{1}+\frac{\alpha_{1} s_{1} p_{1}}{1+\eta s_{1}}+\frac{\alpha_{2} s_{1} y_{1}}{1+\eta s_{1}},$

εy1$\frac{\alpha_{1} s_{1} p_{1}}{1+\eta s_{1}}+\frac{\alpha_{2} s_{1} y_{1}}{1+\eta s_{1}} \text { and } m y_{1}=\gamma p_{1},$

we obtain

$\left(\frac{\alpha_{1} s_{1} p_{1}}{1+\eta s_{1}}-\frac{\alpha_{1} s_{1} p_{1}}{m y_{1}\left(1+\eta s_{1}\right)} \gamma p_{1}\right) \frac{p}{p_{1}}$ = 0

and

$\left(\frac{\alpha_{2} s_{1} y_{1}}{1+\eta s_{1}}+\frac{\alpha_{1} s_{1} p_{1}}{1+\eta s_{1}}-\varepsilon y_{1}\right) \frac{y}{y_{1}}$ = 0.

Then

$\frac{d L_{1}}{d t}=-\hat{\delta} \frac{\left(s-s_{1}\right)^{2}}{s\left(1+\eta s_{1}\right)}+\left(\frac{\alpha_{1} s_{1} p_{1}}{1+\eta s_{1}}+\frac{\alpha_{2} s_{1} y_{1}}{1+\eta s_{1}}\right)\left(1-\frac{s_{1}(1+\eta s)}{s\left(1+\eta s_{1}\right)}\right)$

$-\frac{\alpha_{1} s_{1} p_{1}}{1+\eta s_{1}} \frac{s p y_{1}\left(1+\eta s_{1}\right)}{s_{1} p_{1} y(1+\eta s)}-\frac{\alpha_{2} s_{1} y_{1}}{1+\eta s_{1}} \frac{s\left(1+\eta s_{1}\right)}{s_{1}(1+\eta s)}+2 \frac{\alpha_{1} s_{1} p_{1}}{1+\eta s_{1}}+\frac{\alpha_{2} s_{1} y_{1}}{1+\eta s_{1}}$

$-\frac{\alpha_{1} s_{1} p_{1}}{1+\eta s_{1}} \frac{y p_{1}}{y_{1} p}+\frac{q \alpha_{1} s_{1} p_{1}}{m y_{1}\left(1+\eta s_{1}\right)}\left(p_{1}-\frac{\mu}{r}\right) z+k\left(y_{1}-\frac{h}{g}\right) w.$       (3.7)

Thus

$\frac{d L_{1}}{d t}=-\hat{\delta} \frac{\left(s-s_{1}\right)^{2}}{s\left(1+\eta s_{1}\right)}+\frac{\alpha_{1} s_{1} p_{1}}{1+\eta s_{1}}\left(3-\frac{s_{1}(1+\eta s)}{s\left(1+\eta s_{1}\right)}-\frac{s p y_{1}\left(1+\eta s_{1}\right)}{s_{1} p_{1} y(1+\eta s)}-\frac{y p_{1}}{y_{1} p}\right)$

$\begin{array}{l} +\frac{\alpha_{2} s_{1} y_{1}}{1+\eta s_{1}}\left(2-\frac{s_{1}(1+\eta s)}{s\left(1+\eta s_{1}\right)}-\frac{s\left(1+\eta s_{1}\right)}{s_{1}(1+\eta s)}\right)+\frac{q \alpha_{1} s_{1}}{\gamma\left(1+\eta s_{1}\right)}\left(p_{1}-p_{2}\right) z \\ \end{array}$

+ k (y1 − y3) w.       (3.8)

Thus if $R_{1}^{z}$ ≤ 1, then Π2 does not exist since z2$\frac{\gamma}{q}$ ($R_{1}^{z}$ - 1) ≤ 0. It follows that, $\frac{dz}{dt}$ = r ( p (t) − $\frac{\mu}{r}$) z (t) ≤ 0 for all z > 0. Then p (t) ≤ $\frac{\mu}{r}$ = p2 and p1 ≤ p2. If $R_{1}^{w}$ ≤ 1, then Π3 does not exist since w3$\frac{\varepsilon}{k}$ ($R_{1}^{w}$ − 1) ≤ 0. It follows that, $\frac{dw}{dt}$ = g ( y (t) −$\frac{h}{g}$) w (t) ≤ 0 for all w > 0. Hence, y (t) ≤ $\frac{h}{g}$ = y3 and y1 ≤ y3.

Using AM-GM inequality (3.1), with j = 2 and j = 3 we get $\frac{dL_1}{dt}$ ≤ 0 for all s, y, p, z, w > 0. Moreover, $\frac{dL_1}{dt}$ = 0 when s = s1, y = y1, p = p1, z = w = 0. The solutions of system (2.1)-(2.5) tend to $\grave{Γ}$the largest invariant subset of $\grave{Γ}$1  = { (s, y, p, z, w) : $\frac{dL_1}{dt}$ = 0}. Clearly $\grave{Γ}$= {Π1} . LIP implies that Π1 is GAS.

Theorem 3.3. Let $R_{1}^{z}$ > 1 and $R_{2}^{w}$ ≤ 1, then Π2 is GAS.

Proof. Define L2(s, y, p, z, w) as:

L2 = s - s-$\int_{s_{2}}^{s} \frac{s_{2}(1+\eta \theta)}{\theta\left(1+\eta s_{2}\right)} d \theta+y_{2} H\left(\frac{y}{y_{2}}\right)$

$+\frac{\alpha_{1} s_{2} p_{2}}{m y_{2}\left(1+\eta s_{2}\right)} p_{2} H\left(\frac{p}{p_{2}}\right)+\frac{q \alpha_{1} s_{2} p_{2}}{r m y_{2}\left(1+\eta s_{2}\right)} z_{2} H\left(\frac{z}{z_{2}}\right)+\frac{k}{g} w.$

We note that, L2(s, y, p, z, w) > 0 for all s, y, p, z, w > 0 and L2(s2, y2, p2, z2, 0) = 0. Calculating $\frac{dL_2}{dt}$ along the trajectories of system (2.1)-(2.5), we get

$\frac{d L_{2}}{d t}=\left(1-\frac{s_{2}(1+\eta s)}{s\left(1+\eta s_{2}\right)}\right)\left(\beta-\delta s-\frac{\alpha_{1} s p}{1+\eta s}-\frac{\alpha_{2} s y}{1+\eta s}\right)$

$+\left(1-\frac{y_{2}}{y}\right)\left(\frac{\alpha_{1} s p}{1+\eta s}+\frac{\alpha_{2} s y}{1+\eta s}-\varepsilon y-k w y\right)$

$+\frac{\alpha_{1} s_{2} p_{2}}{m y_{2}\left(1+\eta s_{2}\right)}\left(1-\frac{p_{2}}{p}\right)(m y-\gamma p-q z p)$

$+\frac{q \alpha_{1} s_{2} p_{2}}{r m y_{2}\left(1+\eta s_{2}\right)}\left(1-\frac{z_{2}}{z}\right)(r z p-\mu z)+\frac{k}{g}(g w y-h w).$       (3.9)

Collecting terms of Eq. (3.9) we get

$\frac{d L_{2}}{d t}=\left(1-\frac{s_{2}(1+\eta s)}{s\left(1+\eta s_{2}\right)}\right)(\beta-\hat{\delta} s)+\frac{\alpha_{1} s_{2} p}{1+\eta s_{2}}+\frac{\alpha_{2} s_{2} y}{1+\eta s_{2}}-\varepsilon y$

$-\left(\frac{\alpha_{1} s p}{1+\eta s}+\frac{\alpha_{2} s y}{1+\eta s}\right) \frac{y_{2}}{y}+\varepsilon y_{2}+\frac{\alpha_{1} s_{2} p_{2}}{\left(1+\eta s_{2}\right)} \frac{y}{y_{2}}-\frac{\alpha_{1} s_{2} p_{2}}{m y_{2}\left(1+\eta s_{2}\right)} \gamma p$

$-\frac{\alpha_{1} s_{2} p_{2}}{\left(1+\eta s_{2}\right)} \frac{y p_{2}}{y_{2} p}+\frac{\alpha_{1} s_{2} p_{2}}{m y_{2}\left(1+\eta s_{2}\right)} \gamma p_{2}+\frac{q \alpha_{1} s_{2} p_{2}}{m y_{2}\left(1+\eta s_{2}\right)} p_{2} z$

$-\frac{q \alpha_{1} s_{2} p_{2}}{m y_{2}\left(1+\eta s_{2}\right)} \frac{\mu}{r} z-\frac{q \alpha_{1} s_{2} p_{2}}{m y_{2}\left(1+\eta s_{2}\right)} p z_{2}+\frac{q \alpha_{1} s_{2} p_{2}}{m y_{2}\left(1+\eta s_{2}\right)} \frac{\mu}{r} z_{2}$

$+k y_{2} w-\frac{k h}{g} w$

$=\left(1-\frac{s_{2}(1+\eta s)}{s\left(1+\eta s_{2}\right)}\right)(\beta-\hat{\delta} s)+\frac{\alpha_{1} s_{2} p_{2}}{1+\eta s_{2}} \frac{p}{p_{2}}+\frac{\alpha_{2} s_{2} y_{2}}{1+\eta s_{2}} \frac{y}{y_{2}}-\varepsilon y_{2} \frac{y}{y_{2}}$

$-\frac{\alpha_{1} s_{2} p_{2}}{1+\eta s_{2}} \frac{s p y_{2}\left(1+\eta s_{2}\right)}{s_{2} p_{2} y(1+\eta s)}-\frac{\alpha_{2} s_{2} y_{2}}{1+\eta s} \frac{s\left(1+\eta s_{2}\right)}{s_{2}(1+\eta s)}+\varepsilon y_{2}+\frac{\alpha_{1} s_{2} p_{2}}{\left(1+\eta s_{2}\right)} \frac{y}{y_{2}}$

$-\frac{\alpha_{1} s_{2} p_{2}}{m y_{2}\left(1+\eta s_{2}\right)} \gamma p_{2} \frac{p}{p_{2}}-\frac{\alpha_{1} s_{2} p_{2}}{\left(1+\eta s_{2}\right)} \frac{y p_{2}}{y_{2} p}+\frac{\alpha_{1} s_{2} p_{2}}{m y_{2}\left(1+\eta s_{2}\right)} \gamma p_{2}$

$+\frac{q \alpha_{1} s_{2} p_{2}}{m y_{2}\left(1+\eta s_{2}\right)} p_{2} z-\frac{q \alpha_{1} s_{2} p_{2}}{m y_{2}\left(1+\eta s_{2}\right)} \frac{\mu}{r} z-\frac{q \alpha_{1} s_{2} p_{2}}{m y_{2}\left(1+\eta s_{2}\right)} p_{2} z_{2} \frac{p}{p_{2}}$

$+\frac{q \alpha_{1} s_{2} p_{2}}{m y_{2}\left(1+\eta s_{2}\right)} \frac{\mu}{r} z_{2}+\frac{k h}{g}\left(\frac{g}{h} y_{2}-1\right) w.$

Applying the equilibrium conditions for Π2:

\begin{aligned} \beta &=\hat{\delta} s_{2}+\frac{\alpha_{1} s_{2} p_{2}}{1+\eta s_{2}}+\frac{\alpha_{2} s_{2} y_{2}}{1+\eta s_{2}}, \quad \varepsilon y_{2}=\frac{\alpha_{1} s_{2} p_{2}}{1+\eta s_{2}}+\frac{\alpha_{2} s_{2} y_{2}}{1+\eta s_{2}} \\ m y_{2} &\left.=\gamma p_{2}+q p_{2} z_{2}\right), \quad p_{2}=\frac{\mu}{r} \end{aligned}

we get

$\left(\frac{\alpha_{1} s_{2} p_{2}}{1+\eta s_{2}}-\frac{\alpha_{1} s_{2} p_{2}}{m y_{2}\left(1+\eta s_{2}\right)} \gamma p_{2}-\frac{\alpha_{1} s_{2} p_{2}}{m y_{2}\left(1+\eta s_{2}\right)} q z_{2} p_{2}\right) \frac{p}{p_{2}}$

$=\frac{\alpha_{1} s_{2} p_{2}}{1+\eta s_{2}}\left(1-\frac{\gamma p_{2}+q z_{2} p_{2}}{m y_{2}}\right)=0$

and

$\left(\frac{\alpha_{2} s_{2} y_{2}}{1+\eta s_{2}}+\frac{\alpha_{1} s_{2} p_{2}}{1+\eta s_{2}}-\varepsilon y_{2}\right) \frac{y}{y_{2}}=0.$

Then

$\frac{d L_{2}}{d t}=-\hat{\delta} \frac{\left(s-s_{2}\right)^{2}}{s\left(1+\eta s_{2}\right)}+\left(\frac{\alpha_{1} s_{2} p_{2}}{1+\eta s_{2}}+\frac{\alpha_{2} s_{2} y_{2}}{1+\eta s_{2}}\right)\left(1-\frac{s_{2}(1+\eta s)}{s\left(1+\eta s_{2}\right)}\right)+2 \frac{\alpha_{1} s_{2} p_{2}}{1+\eta s_{2}}$

$+\frac{\alpha_{2} s_{2} y_{2}}{1+\eta s_{2}}-\frac{\alpha_{1} s_{2} p_{2}}{1+\eta s_{2}} \frac{s p y_{2}\left(1+\eta s_{2}\right)}{s_{2} p_{2} y(1+\eta s)}-\frac{\alpha_{2} s_{2} y_{2}}{1+\eta s} \frac{s\left(1+\eta s_{2}\right)}{s_{2}(1+\eta s)}$

$-\frac{\alpha_{1} s_{2} p_{2}}{\left(1+\eta s_{2}\right)} \frac{y p_{2}}{y_{2} p}+k\left(y_{2}-\frac{h}{g}\right) w$

Thus

$\frac{d L_{2}}{d t}=-\hat{\delta} \frac{\left(s-s_{2}\right)^{2}}{s\left(1+\eta s_{2}\right)}+\frac{\alpha_{1} s_{2} p_{2}}{1+\eta s_{2}}\left(3-\frac{s_{2}(1+\eta s)}{s\left(1+\eta s_{2}\right)}-\frac{s p y_{2}}{s_{2} p_{2} y} \frac{1+\eta s_{2}}{1+\eta s}-\frac{y}{y_{2}} \frac{p_{2}}{p}\right)$

$+\frac{\alpha_{2} s_{2} y_{2}}{1+\eta s_{2}}\left(2-\frac{s_{2}(1+\eta s)}{s\left(1+\eta s_{2}\right)}-\frac{s\left(1+\eta s_{2}\right)}{s_{2}(1+\eta s)}\right)+k\left(y_{2}-y_{3}\right) w.$

Thus, if $R_{2}^{w}$ ≤ 1, then Π4 does not exist since w4$\frac{\varepsilon}{k}$($R_{2}^{w}$ − 1) ≤ 0. It follow that $\frac{dw}{dt}$ = g ( y(t) − $\frac{h}{g}$)w(t) ≤ 0 for all w > 0. Then, y (t) ≤ $\frac{h}{g}$ and y2 ≤ $\frac{h}{g}$ = y3. Then by using (AM-GM) inequality (3.1), with j = 2 and j = 3 we have $\frac{dL_2}{dt}$ ≤ 0 for all s, y, p, w > 0. We have $\frac{dL_2}{dt}$ = 0 when s = s2, y = y2, p = p2 and w = 0. Let Γ2 = { (s, y, p, z, w) : $\frac{dL_2}{dt}$ = 0} and $\grave{Γ}$2 be the largest invariant subset of Γ2. The solutions of system (2.1)-(2.5) tend to  $\grave{Γ}$2. For each element of  $\grave{Γ}$we have p(t) = p2 and y(t) = y2. From Eq. (2.3) we get

$\dot{p}$(t) = 0 = my2 − γp2 − qz(t)p2.

z(t) = $\frac{\gamma}{g}\left(\frac{m y_{2}}{\gamma p_{2}}-1\right)$ = z2.

Then $\grave{Γ}$contains a single point that is {Π2} . LIP implies that Π2 is GAS.

Theorem 3.4. If $R_{1}^{w}$ > 1 and $R_{2}^{z}$ 2 ≤ 1, then Π3 is GAS in ∆.

Proof. Define a function L3(s, y, p, z, w) as:

L3 = s − s3 $-\int_{s_{3}}^{s} \frac{s_{3}(1+\eta \theta)}{\theta\left(1+\eta s_{3}\right)} d \theta+y_{3} H\left(\frac{y}{y_{3}}\right)+\frac{\alpha_{1} s_{3} p_{3}}{m y_{3}\left(1+\eta s_{3}\right)} p_{3} H\left(\frac{p}{p_{3}}\right)$

$+\frac{q \alpha_{1} s_{3} p_{3}}{r m y_{3}\left(1+\eta s_{3}\right)} z+\frac{k}{g} w_{3} H\left(\frac{w}{w_{3}}\right)$.

Clearly, L3(s, y, p, z, w) > 0 for all s, y, p, z, w > 0 and L3(s3, y3, p3, 0, w3) = 0. Calculating $\frac{dL_3}{dt}$ along the trajectories of system (2.1)-(2.5), we get

$\frac{d L_{3}}{d t}=\left(1-\frac{s_{3}(1+\eta s)}{s\left(1+\eta s_{3}\right)}\right)\left(\beta-\hat{\delta} s-\frac{\alpha_{1} s p}{1+\eta s}-\frac{\alpha_{2} s y}{1+\eta s}\right)$

$+\left(1-\frac{y_{3}}{y}\right)\left(\frac{\alpha_{1} s p}{1+\eta s}+\frac{\alpha_{2} s y}{1+\eta s}-\varepsilon y-k w y\right)$

$+\frac{\alpha_{1} s_{3} p_{3}}{m y_{3}\left(1+\eta s_{3}\right)}\left(1-\frac{p_{3}}{p}\right)(m y-\gamma p-q z p)$

$+\frac{q \alpha_{1} s_{3} p_{3}}{r m y_{3}\left(1+\eta s_{3}\right)}(r z p-\mu z)+\frac{k}{g}\left(1-\frac{w_{3}}{w}\right)(g w y-h w)$.       (3.10)

Collecting terms of Eq. (3.10) we get

$\frac{d L_{3}}{d t}=\left(1-\frac{s_{3}(1+\eta s)}{s\left(1+\eta s_{3}\right)}\right)(\beta-\hat{\delta} s)+\frac{\alpha_{1} s_{3} p}{1+\eta s_{3}}+\frac{\alpha_{2} s_{3} y}{1+\eta s_{3}}-\varepsilon y$

$-\left(\frac{\alpha_{1} s p}{1+\eta s}+\frac{\alpha_{2} s y}{1+\eta s}\right) \frac{y_{3}}{y}+\varepsilon y_{3}+\frac{\alpha_{2} s_{3} p_{3}}{\left(1+\eta s_{3}\right)} \frac{y}{y_{3}}-\frac{\alpha_{1} s_{3} p_{3}}{m y_{3}\left(1+\eta s_{3}\right)} \gamma p$

$-\frac{\alpha_{1} s_{3} p_{3}}{\left(1+\eta s_{3}\right)} \frac{y p_{3}}{y_{3} p}+\frac{\alpha_{1} s_{3} p_{3}}{m y_{3}\left(1+\eta s_{3}\right)} \gamma p_{3}+\frac{q \alpha_{1} s_{3} p_{3}}{m y_{3}\left(1+\eta s_{3}\right)} p_{3} z$

$-\frac{q \alpha_{1} s_{3} p_{3}}{m y_{3}\left(1+\eta s_{3}\right)} \frac{\mu}{r} z+k w y_{3}-\frac{k h}{g} w-k w_{3} y+\frac{k h}{g} w_{3}.$

$=\left(1-\frac{s_{3}(1+\eta s)}{s\left(1+\eta s_{3}\right)}\right)(\beta-\hat{\delta} s)+\frac{\alpha_{1} s_{3} p_{3}}{1+\eta s_{3}} \frac{p}{p_{3}}+\frac{\alpha_{2} s_{3} y_{3}}{1+\eta s_{3}} \frac{y}{y_{3}}-\varepsilon y_{3} \frac{y}{y_{3}}$

$-\frac{\alpha_{1} s_{3} p_{3}}{1+\eta s_{3}} \frac{s p y_{3}\left(1+\eta s_{3}\right)}{s_{3} p_{3} y(1+\eta s)}-\frac{\alpha_{2} s_{3} y_{3}}{1+\eta s_{3}} \frac{s\left(1+\eta s_{3}\right)}{s_{3}(1+\eta s)}+\varepsilon y_{3}+\frac{\alpha_{1} s_{3} p_{3}}{\left(1+\eta s_{3}\right)} \frac{y}{y_{3}}$

$-\frac{\alpha_{1} s_{3} p_{3}}{m y_{3}\left(1+\eta s_{3}\right)} \gamma p_{3} \frac{p}{p_{3}}-\frac{\alpha_{1} s_{3} p_{3}}{\left(1+\eta s_{3}\right)} \frac{y p_{3}}{y_{3} p}+\frac{\alpha_{1} s_{3} p_{3}}{m y_{3}\left(1+\eta s_{3}\right)} \gamma p_{3}$

$+\frac{q \alpha_{1} s_{3} p_{3}}{m y_{3}\left(1+\eta s_{3}\right)} \frac{\mu}{r}\left(\frac{r}{\mu} p_{3}-1\right) z+k w y_{3}-\frac{k h}{g} w-k w_{3} y_{3} \frac{y}{y_{3}}+\frac{k h}{g} w_{3}.$

Using the equilibrium conditions for Π3:

$\beta=\hat{\delta} s_{3}+\frac{\alpha_{1} s_{3} p_{3}}{1+\eta s_{3}}+\frac{\alpha_{2} s_{3} y_{3}}{1+\eta s_{3}}, \quad \gamma p_{3}=m y_{3},$

$y_{3}=\frac{h}{g}, \quad \frac{\alpha_{1} s_{3} p_{3}}{1+\eta s_{3}}+\frac{\alpha_{2} s_{3} y_{3}}{1+\eta s_{3}}=\varepsilon y_{3}+k y_{3} w_{3}.$

we obtain

$\left(\frac{\alpha_{1} s_{3} p_{3}}{1+\eta s_{3}}-\frac{\alpha_{1} s_{3} p_{3}}{m y_{3}\left(1+\eta s_{3}\right)} \gamma p_{3}\right) \frac{p}{p_{3}}=0$

and

$\left(\frac{\alpha_{2} s_{3} y_{3}}{1+\eta s_{3}}+\frac{\alpha_{1} s_{3} p_{3}}{1+\eta s_{3}}-\varepsilon y_{3}-k w_{3} y_{3}\right) \frac{y}{y_{3}}=0.$

Thus

$\frac{d L_{3}}{d t}=-\hat{\delta} \frac{\left(s-s_{3}\right)^{2}}{s\left(1+\eta s_{3}\right)}+\left(\frac{\alpha_{1} s_{3} p_{3}}{1+\eta s_{3}}+\frac{\alpha_{2} s_{3} y_{3}}{1+\eta s_{3}}\right)\left(1-\frac{s_{3}(1+\eta s)}{s\left(1+\eta s_{3}\right)}\right)+2 \frac{\alpha_{1} s_{3} p_{3}}{1+\eta s_{3}}$

$+\frac{\alpha_{2} s_{3} y_{3}}{1+\eta s_{3}}\left(2-\frac{s_{3}(1+\eta s)}{s\left(1+\eta s_{3}\right)}-\frac{s\left(1+\eta s_{3}\right)}{s_{3}(1+\eta s)}\right)+\frac{q \alpha_{1} s_{3} p_{3}}{m y_{3}\left(1+\eta s_{3}\right)} \frac{\mu}{r}\left(R_{2}^{z}-1\right) z.$

Therefore, if $R_{2}^{z}$ ≤ 1, then using (AM-GM) inequality (3.1), with j = 2 and j = 3 we have $\frac{dL_3}{dt}$ ≤ 0 for all s, y, p, z > 0 and $\frac{dL_3}{dt}$ = 0 when s = s3, y = y3, p = p3 and z = 0. Let Γ3 = { (s, y, p, z, w) : $\frac{dL_3}{dt}$ = 0 } and $\grave{Γ}$3 be the largest invariant subset of Γ3. The solutions of system (2.1)-(2.5) tend to  $\grave{Γ}$3. For each element of  $\grave{Γ}$we have y(t) = y3. From Eq. (2.2) we get

$\dot{y}(t)=0=\frac{\alpha_{2} s_{3} y_{3}}{1+\eta s_{3}}+\frac{\alpha_{1} s_{3} p_{3}}{1+\eta s_{3}}-\varepsilon y_{3}-k w(t) y_{3}.$

Hence

\begin{aligned} w(t) &=\frac{\varepsilon}{k}\left(\frac{\left(m \alpha_{1}+\gamma \alpha_{2}\right) s_{3}}{\varepsilon \gamma\left(1+\eta s_{3}\right)}-1\right) \\ &=\frac{\varepsilon}{k}\left(\frac{g \beta\left(m \alpha_{1}+\gamma \alpha_{2}\right)}{\varepsilon\left(h\left(m \alpha_{1}+\gamma \alpha_{2}\right)+g \gamma(\beta \eta+\hat{\delta})\right)}-1\right)=w_{3}. \end{aligned}

Then $\grave{Γ}$3 contains a single point that is {Π3} . Hence, global stability of Π3 follows from LIP.

Theorem 3.5. For system (2.1)-(2.5), suppose that $R_{2}^{w}$ > 1 $R_{2}^{z}$ > 1, then Π4 is GAS in $\stackrel{\circ}{\Delta}$.

Proof. We construct a function L4(s, y, p, z, w) as:

L4 = s − s4 $-\int_{s_{4}}^{s} \frac{s_{4}(1+\eta \theta)}{\theta\left(1+\eta s_{4}\right)} d \theta+y_{4} H\left(\frac{y}{y_{4}}\right)+\frac{\alpha_{1} s_{4} p_{4}}{m y_{4}\left(1+\eta s_{4}\right)} p_{4} H\left(\frac{p}{p_{4}}\right)$

$+\frac{q \alpha_{1} s_{4} p_{4}}{r m y_{4}\left(1+\eta s_{4}\right)} z_{4} H\left(\frac{z}{z_{4}}\right)+\frac{k}{g} w_{4} H\left(\frac{w}{w_{4}}\right).$

It is seen that, L4(s, y, p, z, w) > 0 for all s, y, p, z, w > 0 and L4(s4, y4, p4, z4, w4) = 0. We Calculating $\frac{dL_4}{dt}$ along the trajectories of system (2.1)-(2.5), we get

$\frac{d L_{4}}{d t}=\left(1-\frac{s_{4}(1+\eta s)}{s\left(1+\eta s_{4}\right)}\right)\left(\beta-\hat{\delta} s-\frac{\alpha_{1} s p}{1+\eta s}-\frac{\alpha_{2} s y}{1+\eta s}\right)$

$+\left(1-\frac{y_{4}}{y}\right)\left(\frac{\alpha_{1} s p}{1+\eta s}+\frac{\alpha_{2} s y}{1+\eta s}-\varepsilon y-k w y\right)$

$+\frac{\alpha_{1} s_{4} p_{4}}{m y_{4}\left(1+\eta s_{4}\right)}\left(1-\frac{p_{4}}{p}\right)(m y-\gamma p-q z p)$

$+\frac{\alpha_{1} s_{4} p_{4}}{m y_{4}\left(1+\eta s_{4}\right)}\left(1-\frac{p_{4}}{p}\right)(m y-\gamma p-q z p)$

$+\frac{q \alpha_{1} s_{4} p_{4}}{r m y_{4}\left(1+\eta s_{4}\right)}\left(1-\frac{z_{4}}{z}\right)(r z p-\mu z)$

$+\frac{k}{g}\left(1-\frac{w_{4}}{w}\right)(g w y-h w)$.       (3.12)

Collecting terms of Eq. (3.12) we get

$\frac{d L_{4}}{d t}=\left(1-\frac{s_{4}(1+\eta s)}{s\left(1+\eta s_{4}\right)}\right)(\beta-\hat{\delta} s)+\frac{\alpha_{1} s_{4} p_{4}}{1+\eta s_{4}} \frac{p}{p_{4}}+\frac{\alpha_{2} s_{4} y_{4}}{1+\eta s_{4}} \frac{y}{y_{4}}$

$-\varepsilon y_{4} \frac{y}{y_{4}} \frac{\alpha_{1} s_{4} p_{4}}{1+\eta s_{4}} \frac{s p y_{4}}{s_{4} p_{4} y} \frac{1+\eta s_{4}}{1+\eta s}-\frac{\alpha_{2} s_{4} y_{4}}{1+\eta s_{4}} \frac{s\left(1+\eta s_{4}\right)}{s_{4}(1+\eta s)}+\varepsilon y_{4}$

$+\frac{\alpha_{1} s_{4} p_{4}}{\left(1+\eta s_{4}\right)} \frac{y}{y_{4}}-\frac{\alpha_{1} s_{4} p_{4}}{m y_{4}\left(1+\eta s_{4}\right)} \gamma p_{4} \frac{p}{p_{4}}-\frac{\alpha_{1} s_{4} p_{4}}{\left(1+\eta s_{4}\right)} \frac{y p_{4}}{y_{4} p}$

$+\frac{\alpha_{1} s_{4} p_{4}}{m y_{4}\left(1+\eta s_{4}\right)} \gamma p_{4}+\frac{q \alpha_{1} s_{4} p_{4}}{m y_{4}\left(1+\eta s_{4}\right)} p_{4} z-\frac{q \alpha_{1} s_{4} p_{4}}{m y_{4}\left(1+\eta s_{4}\right)} \frac{\mu}{r} z.$

Applying the equilibrium conditions for Π4:

$\beta=\hat{\delta} s_{4}+\frac{\alpha_{1} s_{4} p_{4}}{1+\eta s_{4}}+\frac{\alpha_{2} s_{4} y_{4}}{1+\eta s_{4}}, \quad \varepsilon y_{4}=\frac{\alpha_{1} s_{4} p_{4}}{1+\eta s_{4}}+\frac{\alpha_{2} s_{4} y_{4}}{1+\eta s_{4}}-k y_{4} w_{4},$

$m y_{4}=\gamma p_{4}+q p_{4} z_{4}, \quad p_{4}=\frac{\mu}{r}, \quad y_{4}=\frac{h}{g}$

we obtain

$\frac{\alpha_{1} s_{4} p}{1+\eta s_{4}}-\frac{\alpha_{1} s_{4} p_{4}}{m y_{4}\left(1+\eta s_{4}\right)} \gamma p-\frac{\alpha_{1} s_{4} p_{4}}{m y_{4}\left(1+\eta s_{4}\right)} q z_{4} p$

$=\left(\frac{\alpha_{1} s_{4} p_{4}}{1+\eta s_{4}}-\frac{\alpha_{1} s_{4} p_{4}}{m y_{4}\left(1+\eta s_{4}\right)} \gamma p_{4}-\frac{\alpha_{1} s_{4} p_{4}}{m y_{4}\left(1+\eta s_{4}\right)} q z_{4} p_{4}\right) \frac{p}{p_{4}}=0$

and

$\frac{\alpha_{2} s_{4} y}{1+\eta s_{4}}-\varepsilon y+\frac{\alpha_{1} s_{4} p_{4}}{\left(1+\eta s_{4}\right) y_{4}} y-k w_{4} y$

$=\left(\frac{\alpha_{2} s_{4} y_{4}}{1+\eta s_{4}}-\varepsilon y_{4}+\frac{\alpha_{1} s_{4} p_{4}}{1+\eta s_{4}}-k w_{4} y_{4}\right) \frac{y}{y_{4}}=0$

Then

$\frac{d L_{4}}{d t}=-\hat{\delta} \frac{\left(s-s_{4}\right)^{2}}{s\left(1+\eta s_{4}\right)}+\left(\frac{\alpha_{1} s_{4} p_{4}}{1+\eta s_{4}}+\frac{\alpha_{2} s_{4} y_{4}}{1+\eta s_{4}}\right)\left(1-\frac{s_{4}(1+\eta s)}{s\left(1+\eta s_{4}\right)}\right)$

$+2 \frac{\alpha_{1} s_{4} p_{4}}{1+\eta s_{4}}+\frac{\alpha_{2} s_{4} y_{4}}{1+\eta s_{4}}-\frac{\alpha_{1} s_{4} p_{4}}{1+\eta s_{4}} \frac{s p y_{4}\left(1+\eta s_{4}\right)}{s_{4} p_{4} y(1+\eta s)}$

$-\frac{\alpha_{1} s_{4} y_{4}}{1+\eta s} \frac{s\left(1+\eta s_{4}\right)}{s_{4}(1+\eta s)}-\frac{\alpha_{1} s_{4} p_{4}}{\left(1+\eta s_{4}\right)} \frac{y p_{4}}{y_{4} p}.$       (3.13)

Eq. (3.13) can be simplified as:

$\frac{d L_{4}}{d t}=-\hat{\delta} \frac{\left(s-s_{4}\right)^{2}}{s\left(1+\eta s_{4}\right)}+\frac{\alpha_{1} s_{4} p_{4}}{1+\eta s_{4}}\left(3-\frac{s_{4}(1+\eta s)}{s\left(1+\eta s_{4}\right)}-\frac{s p y_{4}\left(1+\eta s_{4}\right)}{s_{4} p_{4} y(1+\eta s)}-\frac{y p_{4}}{y_{4} p}\right)$

$+\frac{\alpha_{2} s_{4} y_{4}}{1+\eta s_{4}}\left(2-\frac{s_{4}(1+\eta s)}{s\left(1+\eta s_{4}\right)}-\frac{s\left(1+\eta s_{4}\right)}{s_{4}(1+\eta s)}\right).$

Using (AM-GM) inequality (3.1), with j = 2 and j = 3 we have, $dL4 \over dt$ ≤ 0 for all s, y, p > 0 and $dL4 \over dt$ = 0 when s = s4, y = y4, p = p4. Let Γ4 = { (s, y, p, z, w) : $dL4 \over dt$= 0} and $\grave{Γ}$be the largest invariant subset of Γ4. The solutions of system (2.1)-(2.5) tend to$\grave{Γ}$4. For each element of $\grave{Γ}$4 we have y(t) = y4. From Eq. (2.2) we get

$\dot{y}(t)=0=\frac{\alpha_{2} s_{4} y_{4}}{1+\eta s_{4}}+\frac{\alpha_{1} s_{4} p_{4}}{1+\eta s_{4}}-\varepsilon y_{4}-k w(t) y_{4}$

Hence

$w(t)=\frac{\varepsilon}{k}\left(\frac{\left(\mu g \alpha_{1}+r h \alpha_{2}\right) s_{4}}{\varepsilon\left(1+\eta s_{3}\right)}-1\right)$

$=\frac{\varepsilon}{k}\left(\frac{\beta g\left(\mu g \alpha_{1}+r h \alpha_{2}\right)}{\varepsilon h\left(g \mu \alpha_{1}+h r \alpha_{2}+r g(\beta \eta+\hat{\delta})\right)}-1\right)=w_{4}$

From Eq. (2.3) we get

$\dot{p}$(t) = 0 = my4 − γp4 − qz(t)p4.

Hence

z(t) = $\frac{\gamma}{g}\left(\frac{m h r}{\gamma \mu g}-1\right)$ = z4.

Then $\grave{Γ}$contains a single point that is {Π4} . LIP implies that Π4 is GAS.

# 4. NUMERICAL SIMULATIONS

In this section we perform some numerical simulations for model (2.1)-(2.5), with parameters values given in Table 1. In the figures we show the evolution of the five states of the system s, y, p, z and w. We have used MATLAB for all computations. Now we investigate our

TABLE 1. Some parameters and their values of model (2.1)-(2.5).

theoretical results given in Theorems 3.1 - 3.5.

4.1. Stability of equilibria. In this subsection, we take η = 0 and chose three different initial conditions for model (2.1)-(2.5) as follows:

IC1 : (s(0), y(0), p(0), z(0), w(0)) = (600, 5, 0.5, 5, 0.5), (Solid lines in the figures),

IC2: (s(0), y(0), p(0), z(0), w(0)) = (400, 10, 3, 8, 1), (Dashed lines in the figures),

IC3: (s(0), y(0), p(0), z(0), w(0)) = (200, 15, 4, 14, 2). (Dotted lines in the figures),

Scenario 1: α1 = α2 = 0.0001, r = 0.01 and g = 0.01. For this set of parameters, we have R0 = 0.6667 < 1. Figure 1 shows that, the solutions of the system with IC1-IC3 converge to Π0 = (1000, 0, 0, 0, 0, 0). According to Theorem 3.1, Π0 is GAS.

FIGURE 1. The simulation of trajectories of system (2.1)-(2.5) for scenario 1.

Scenario 2: α1 = α2 = 0.0003, r = 0.001 and g = 0.001. With such choice we get, $R^z_1$= 0.33 < 1 and $R^w_1$ = 0.22 < 1 < R0 = 2 and Π1 exists with Π1 = (500, 12.5, 20.83, 0, 0). This result supports Lemma 1. Figure 2 support Theorem 3.2 that, Π1 is GAS.

FIGURE 2. The simulation of trajectories of system (2.1)-(2.5) for scenario 2.

Scenario 3: α1 = α2 = 0.0003, r = 0.006 and g = 0.001. Then, we calculate R0 = 2 > 1, $R^z_1$= 1.14 > 1 and $R^w_2$ = 0.19 < 1. Figure 3 shows that the solution of the system with different intial conditions reach the equilibrium Π2 = (542.57, 11.44, 16.67, 2.15, 0). This support Theorem 3.3.

FIGURE 3. The simulation of trajectories of system (2.1)-(2.5) for scenario 3.

Scenario 4: α1 = α2 = 0.0003, r = 0.001 and g = 0.01. Then, we calculate R0 = 2 > 1, $R^w_1$ = 1.11 > 1 and $R^z_2$ = 0.17 < 1. The results presented in Lemma 2.2 and Theorem 3.4 show that the equilibrium Π3 exists and it is GAS. Figure 4 supports the theoretical results of Theorem 3.4, where the states of the system reach the equilibrium Π3 = (555.56, 10, 16.67, 0, 0.44), for all initial conditions.

FIGURE 4. The simulation of trajectories of system (2.1)-(2.5) for scenario 4.

Scenario 5: α1 = α2 = 0.0004, r = 0.01 and g = 0.01. Then, we calculate R0 = 2 > 1 and $R^z_2$ = 1.67 > 1 and $R^w_2$ = 1.11 > 1. According to Lemma 2.2 and Theorem 3.5, Π4 exists and it is GAS. Figure 5, confirms the results of Theorem 3.5 where the states of the system starting with different initials converge to the equilibrium Π4 = (555.56, 10, 10, 10, 0.44).

FIGURE 5. The simulation of trajectories of system (2.1)-(2.5) for scenario 5.

4.2. Effect of the Holling type-II parameter η on the pathogen dynamics. Let us take the initial conditions (IC2). We choose the values α1 = α2 = 0.006, r = 0.01 and g = 0.01 and η is varied. Figure 6 shows the effect of the Holling type-II incidence η on the stability of the equilibria of the system. We observe that, as η is increased, both the virus-target and infectedtarget infection rates are decreased, and then the concentration of the uninfected (susceptible host) cells is increased, while the concentrations of the infected cells and free virus particles (pathogens), B cells and CTL cells are decreased. We note that R0 is a decreasing function of η. Now we compute ηcr such thatR0 = 1 = $\frac{\beta\left(m \alpha_{1}+\gamma \alpha_{2}\right)}{\varepsilon \gamma\left(\beta \eta_{c r}+\hat{\delta}\right)}$. Then ηcr = $\frac{m \alpha_{1}+\gamma \alpha_{2}}{\varepsilon \gamma} -\frac{1}{s_{0}}$. Therefore

FIGURE 6. The effect of holling rate constant η on the behaviour of all trajectories of system (2.1)-(2.5).

ηcr = max $\left\{0, \frac{m \alpha_{1}+\gamma \alpha_{2}}{\varepsilon \gamma}-\frac{1}{s_{0}}\right\}.$

It follows that, if η ≥ ηcr, then Π0 is GAS. For the choice of values of the parameters given above we found that ηcr = 0.039.

# 5. CONCLUSION

In this paper, we have studied a virus dynamics model with both virus-to-cell and cell-to-cell transmissions. The effect of both CTL and antibodies on the virus dynamics have been studied. The virus-uninfected and infected-uninfected incidence rates have been given by Holling type-II. We have shown that, the solutions of the model are nonnegative and bounded which ensure the well-posedness of the model. We have derived five threshold numbers which fully determines the existence and stability of the five equilibria of the model. We have investigated the global stability of the equilibria of the model by using Lyapunov method and LaSalle’s invariance principle. We have conducted numerical simulations and have shown that both the theoretical and numerical results are consistent. The results show that the CTL and antibodies can control the disease progression by reducing the concentration of the free virus particles and infected cells. Our proposed model can be extended by incorporating different types of time delay. Moreover, following the work of Gibelli et al. [54], viral infection models adaptive immune response and with a stochastic parameters dynamics can also be studied.

#### References

1. M. A. Nowak, R. Anderson, M. Boerlijst, S. Bonhoe er, R. May, and A. McMichael, HIV-1 evolution and disease progression, Science, 274 (1996), 1008-1010. https://doi.org/10.1126/science.274.5289.1008
2. M. Y. Li, and H. Shu, Global dynamics of a mathematical model for HTLV-I infection of CD4+ T cells with delayed CTL response, Nonlinear Anal. Real World Appl. 13 (2012), 1080-1092. https://doi.org/10.1016/j.nonrwa.2011.02.026
3. C. Lv, L. Huang, and Z. Yuan, Global stability for an HIV-1 infection model with Beddington-DeAngelis incidence rate and CTL immune response, Commun. Nonlinear Sci. Numer. Simul. 19 (2014) 121-127. https://doi.org/10.1016/j.cnsns.2013.06.025
4. A. M. Elaiw, A. A. Raezah and S. A. Azoz, Stability of delayed HIV dynamics models with two latent reservoirs and immune impairment, Adv. Differ. Equ. (2018).
5. X. Li and S. Fu, Global stability of a virus dynamics model with intracellular delay and CTL immune response, Math. Meth. Appl. Sci. 38 (2015), 420-430. https://doi.org/10.1002/mma.3078
6. D. Huang, X. Zhang, Y. Guo, and H. Wang, Analysis of an HIV infection model with treatments and delayed immune response, Appl. Math. Model. 40 (2016), 3081-3089. https://doi.org/10.1016/j.apm.2015.10.003
7. Y. Zhao and Z. Xu, Global dynamics for a delyed hepatitis C virus,infection model, Electron. J. Differ. Equ. 2014 (2014), 1-18. https://doi.org/10.1186/1687-1847-2014-1
8. A. Murase, T. Sasaki, and T. Kajiwara, Stability analysis of pathogen-immune interaction dynamics, J. Math. Biol. 51 (2005), 247-267. https://doi.org/10.1007/s00285-005-0321-y
9. A. M. Elaiw, E. Kh. Elnahary, E. Kh, A. M. Shehata and M. Abul-Ez, Stability of delay-distributed HIV infection models with multiple viral producer cells, J. Korean Soc. Ind. Appl. Math. 22 (2018), 29-62. https://doi.org/10.12941/JKSIAM.2018.22.029
10. A. M. Elaiw, E. Kh. Elnahary and A. A. Raezah, Effect of cellular reservoirs and delays on the global dynamics of HIV, Adv. Differ. Equ. (2018).
11. S.Wang, and D. Zou, Global stability of in host viral models with humoral immunity and intracellular delays, Appl. Math. Model. 36 (2012), 1313-1322. https://doi.org/10.1016/j.apm.2011.07.086
12. A. M. Elaiw, S. A. Ghaleb and A. Hobiny, Effect of time delay and antibodies on HCV dynamics with cure rate and two routes of infection. Journal of Applied Mathematics and Physics, 6 (2018), 1120. https://doi.org/10.4236/jamp.2018.65096
13. D.S. Callaway, and A.S. Perelson, HIV-1 infection and low steady state viral loads, Bull. Math. Biol. 64 (2002), 29-64. https://doi.org/10.1006/bulm.2001.0266
14. A. M. Elaiw and N. H. AlShamrani, Global stability of humoral immunity virus dynamics models with nonlinear infection rate and removal, Nonlinear Anal.-Real World Appl. 26, (2015), 161-190. https://doi.org/10.1016/j.nonrwa.2015.05.007
15. N. E. Tarfulea, A mathematical model for CTL effect on a latently infected cell inclusive HIV dynamics and treatment, AIP Conference Proceedings, 1895 (2017).
16. A. M. Elaiw and N. H. AlShamrani, Stability of a general delay-distributed virus dynamics model with multistaged infected progression and immune response, Math. Meth. Appl. Sci. 40 (2017), 699-719. https://doi.org/10.1002/mma.4002
17. Elaiw, A. M., A. A. Raezah, and A. S. Alofi, Effect of humoral immunity on HIV-1 dynamics with virus-totarget and infected-to-target infections, AIP Adv. 6 (2016).
18. T. Wang, Z. Hu, F. Liao and W. Ma, Global stability analysis for delayed virus infection model with general incidence rate and humoral immunity, Math. Comput. Simul. 89 (2013), 13-22. https://doi.org/10.1016/j.matcom.2013.03.004
19. L. Li and R. Xu, Global dynamics of an age-structured in-host viral infection model with humoral immunity, Adv. Differ. Equ. (2016).
20. J. Xu, Y. Zhou, Y. Li and Y. Yang, Global dynamics of a intracellular infection model with delays and humoral immunity, Math. Meth. Appl. Sci. 39 (2016), 5427-5435. https://doi.org/10.1002/mma.3927
21. H. Miao, Z. Teng, C. Kang, A. Muhammadhaji, Stability analysis of a virus infection model with humoral immunity response and two time delays, Math. Meth. Appl. Sci. 39 (2016), 3434-3449. https://doi.org/10.1002/mma.3790
22. D. Wodarz, Killer cell dynamics: mathematical and computational approaches to immunology. Springer Verlag, New York, (2007).
23. J. Wang, J. Pang, T. Kuniya and Y. Enatsu, Global threshold dynamics in a five-dimensional virus model with cell-mediated, humoral immune responses and distributed delays, Appl. Math. Comput. 241 (2014), 298-316. https://doi.org/10.1016/j.amc.2014.05.015
24. A. M. Elaiw and N. H. AlShamrani, Stability of an adaptive immunity pathogen dynamics model with latency and multiple delays, Math. Meth. Appl. Sci. 36 (2018), 125-142. https://doi.org/10.1002/mma.2576
25. A. M. Elaiw and N. H. AlShamrani, Stability of latent pathogen infection model with adaptive immunity and delays, J. Integr. Neurosci. 17 (2018), 547-576. https://doi.org/10.3233/JIN-180087
26. K. Hattaf and N. Yousfi, A class of delayed viral infection models with general incidence rate and adaptive immune response, International Journal of Dynamics and Control, 4 (2016), 254-265. https://doi.org/10.1007/s40435-015-0158-1
27. A. Rezounenko, Stability of a viral infection model with state-dependent delay, CTL and antibody immune responses, Discrete Contin. Dyn. Syst.-Ser. B. 22 (2016).
28. Y. Yan and W.Wang, Global stability of a five-dimensional model with immune responses and delay, Discrete Contin. Dyn. Syst.-Ser. B. 17 (2012), 401-416. https://doi.org/10.3934/dcdsb.2012.17.401
29. N. Yousfi, K. Hattaf, A. Tridane, Modeling the adaptive immune response in HBV infection, J. Math. Biol. 63 (2011), 933-957. https://doi.org/10.1007/s00285-010-0397-x
30. R.V. Culshaw, S. Ruan, G. Webb, A mathematical model of cell-to-cell spread of HIV-1 that includes a time delay, J. Math. Biol. 46 (2003), 425-444. https://doi.org/10.1007/s00285-002-0191-5
31. X. Lai and X. Zou, Modeling HIV-1 virus dynamics with both virus-to-cell infection and cell-to-cell transmission, SIAM Journal of Applied Mathematics, 74 (2014), 898-917. https://doi.org/10.1137/130930145
32. H. Pourbashash, S.S. Pilyugin, P. De Leenheer, C. McCluskey, Global analysis of within host virus models with cell-to-cell viral transmission, Discrete Contin. Dyn. Syst.-Ser. B. 10 (2014) 3341-3357.
33. J. Wang, J. Lang, X. Zou, Analysis of an age structured HIV infection model with virus-to-cell infection and cell-to-cell transmission, Nonlinear Anal.-Real World Appl. 34 (2017), 75-96. https://doi.org/10.1016/j.nonrwa.2016.08.001
34. S.-S. Chen, C.-Y. Cheng, Y. Takeuchi, Stability analysis in delayed within-host viral dynamics with both viral and cellular infections, J. Math. Anal. Appl. 442 (2016), 642-672. https://doi.org/10.1016/j.jmaa.2016.05.003
35. J. Lin, R. Xu, X. Tian, Threshold dynamics of an HIV-1 virus model with both virus-to-cell and cell-to-cell transmissions, intracellular delay, and humoral immunity, Appl. Math. Comput. 315 (2017), 516-530.
36. Y. Yang, L. Zou and S. Ruan, Global dynamics of a delayed within-host viral infection model with both virusto- cell and cell-to-cell transmissions, Math. Biosci. 270 (2015), 183-191. https://doi.org/10.1016/j.mbs.2015.05.001
37. A. M. Elaiw and A. A. Raezah, Stability of general virus dynamics models with both cellular and viral infections and delays, Math. Meth. Appl. Sci. 40 (2017), 5863-5880. https://doi.org/10.1002/mma.4436
38. S. Pan, S.P. Chakrabarty, Threshold dynamics of HCV model with cell-to-cell transmission and a non-cytolytic cure in the presence of humoral immunity, Commun. Nonlinear Sci. Numer. Simul. 61 (2018), 180-197. https://doi.org/10.1016/j.cnsns.2018.02.010
39. A. Korobeinikov, Global properties of basic virus dynamics models, Bull. Math. Biol. 66 (2004), 879-883. https://doi.org/10.1016/j.bulm.2004.02.001
40. A. M. Elaiw, A. A. Raezah and B. S. Alofi, Dynamics of delayed pathogen infection models with pathogenic and cellular infections and immune impairment, AIP Adv. 8 (2018).
41. A. M. Elaiw and S.A. Azoz, Global properties of a class of HIV infection models with Beddington-DeAngelis functional response, Math. Meth. Appl. Sci. 36 (2013), 383-394. https://doi.org/10.1002/mma.2596
42. A. M. Elaiw, Global properties of a class of HIV models, Nonlinear Anal.-RealWorld Appl. 11 (2010), 2253-2263. https://doi.org/10.1016/j.nonrwa.2009.07.001
43. H. Shu, L. Wang and J. Watmough, Global stability of a nonlinear viral infection model with infinitely distributed intracellular delays and CTL imune responses, SIAM Journal of Applied Mathematics, 73 (2013), 1280-1302. https://doi.org/10.1137/120896463
44. A. M. Elaiw and N. H. AlShamrani, Stability of a general delay-distributed virus dynamics model with multistaged infected progression and immune response, Math. Meth. Appl. Sci. 40 (2017), 699-719. https://doi.org/10.1002/mma.4002
45. A. M. Elaiw and N. A. Almuallem, Global dynamics of delay-distributed HIV infection models with differential drug efficacy in cocirculating target cells, Math. Meth. Appl. Sci. 39 (2016), 4-31. https://doi.org/10.1002/mma.3453
46. A. M. Elaiw, I. A. Hassanien and S. A. Azoz, Global stability of HIV infection models with intracellular delays, J. Korean. Math. Soc. 49 (2012), 779-794. https://doi.org/10.4134/JKMS.2012.49.4.779
47. A. M. Elaiw, Global dynamics of an HIV infection model with two classes of target cells and distributed delays, Discrete Dyn. Nat. Soc. 2012 (2012).
48. A. D. Hobiny, A. M. Elaiw, Almatrafi, Stability of delayed pathogen dynamics models with latency and two routes of infection, Adv. Differ. Equ. (2018).
49. A. M. Elaiw and A. A. Raezah, Stability of general virus dynamics models with both cellular and viral infections and delays, Math. Meth. Appl. Sci. 40 (2017), 5863-5880. https://doi.org/10.1002/mma.4436
50. A. M. Elaiw, T. O. Alade, S. M. Alsulami, Analysis of latent CHIKV dynamics models with general incidence rate and time delays. J. Biol. Dyn. 12 (2018), 700-730. https://doi.org/10.1080/17513758.2018.1503349
51. A. M. Elaiw, T. O. Alade, S. M. Alsulami, Analysis of within-host CHIKV dynamics models with general incidence rate. Int. J. Biomath. 11 (2018).
52. A.M. Elaiw, Global properties of a class of virus infection models with multitarget cells, Nonlinear Dyn. 69 (2012), 423-435. https://doi.org/10.1007/s11071-011-0275-0
53. J.K. Hale, and S. V. Lunel, Introduction to functional differential equations, Springer-Verlag, New York, (1993).
54. L. Gibelli, A. Elaiw, M.A. Alghamdi and A.M. Althiabi, Heterogeneous population dynamics of active particles: Progression, mutations, and selection dynamics, Math. Models Meth. Appl. Sci. 27 (2017), 617-640. https://doi.org/10.1142/S0218202517500117