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

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 L_i , 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 R₀ ≤ 1, then there exists only one equilibrium Π₀,

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

(iii) if \(R_{1}^{z}\) > 1 and \(R_{2}^{w}\) ≤ 1, then there exist only three equilibria Π₀, Π₁ and Π₂,

(iv) if \(R_{1}^{w}\) > 1 and \(R_{2}^{z}\) ≤ 1, then there exist only three equilibria Π₀, Π₁ and Π₃, and

(v) if \(R_{2}^{w}\) > 1 and \(R_{2}^{z}\) > 1, then there exist five equilibria Π₀, Π₁, Π₂, Π₃ and Π₄.

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 Π₁ = (s₁, y₁, p₁, 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 Π₁ exists if R₀ > 1.

(iii) Chronic-infection equilibrium with only humoral immune response Π₂ = (s₂, y₂, p₂, z₂, 0). Now we show that s₂, y₂, p₂ and z₂ 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 s₂ ∈ (0, s₀) such as G(s₂) = 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 Π₃ = (s₃, y₃, p₃, ₀, w₃), 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 Π₄ = (s₄, y₄, p₄, z₄, w₄), where

s₄ = \(\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}\),

w₄ = \(\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, Π₄ 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 θ₁ = θ₂ = ... = θ_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 R₀ ≤ 1, then Π₀ 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(s₀, 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}\)s₀, 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 = s₀, 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 < R₀, then Π₁ is GAS.

Proof. Let us define a function L₁(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, L₁(s, y, p, z, w) > 0 for all s, y, p, z, w > 0 and L1(s₁, y₁, p₁, 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}},\)

εy₁ = \(\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 z₂ = \(\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 Π₃ does not exist since w₃ = \(\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}\) = y₃ and y₁ ≤ y₃.

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 = s₁, y = y₁, p = p₁, z = w = 0. The solutions of system (2.1)-(2.5) tend to \(\grave{Γ}\)₁ the largest invariant subset of \(\grave{Γ}\)₁  = { (s, y, p, z, w) : \(\frac{dL_1}{dt}\) = 0}. Clearly \(\grave{Γ}\)₁ = {Π₁} . LIP implies that Π₁ is GAS.

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

Proof. Define L₂(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, L₂(s, y, p, z, w) > 0 for all s, y, p, z, w > 0 and L₂(s₂, y₂, p₂, z₂, 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 Π₂:

\(\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 w₄ = \(\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 y₂ ≤ \(\frac{h}{g}\) = y₃. 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{Γ}\)₂ be the largest invariant subset of Γ₂. The solutions of system (2.1)-(2.5) tend to  \(\grave{Γ}\)₂. 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)\) = z₂.

Then \(\grave{Γ}\)₂ contains a single point that is {Π₂} . LIP implies that Π₂ is GAS.

Theorem 3.4. If \(R_{1}^{w}\) > 1 and \(R_{2}^{z}\) 2 ≤ 1, then Π₃ is GAS in ∆.

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

L₃ = 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 L₃(s₃, y₃, p₃, 0, w₃) = 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 Π₃:

\(\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 = s₃, y = y₃, p = p₃ and z = 0. Let Γ₃ = { (s, y, p, z, w) : \(\frac{dL_3}{dt}\) = 0 } and \(\grave{Γ}\)₃ be the largest invariant subset of Γ₃. The solutions of system (2.1)-(2.5) tend to  \(\grave{Γ}\)₃. For each element of  \(\grave{Γ}\)₃ we have y(t) = y₃. 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{Γ}\)₃ contains a single point that is {Π₃} . Hence, global stability of Π₃ follows from LIP.

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

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

L₄ = s − s₄ \(-\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 L₄(s₄, y₄, p₄, z₄, w₄) = 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 Π₄:

\(\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{Γ}\)₄. For each element of \(\grave{Γ}\)₄ 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 {Π₄} . LIP implies that Π₄ 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: α₁ = α₂ = 0.0001, r = 0.01 and g = 0.01. For this set of parameters, we have R₀ = 0.6667 < 1. Figure 1 shows that, the solutions of the system with IC1-IC3 converge to Π₀ = (1000, 0, 0, 0, 0, 0). According to Theorem 3.1, Π₀ is GAS.

 

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

Scenario 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 < R₀ = 2 and Π₁ exists with Π₁ = (500, 12.5, 20.83, 0, 0). This result supports Lemma 1. Figure 2 support Theorem 3.2 that, Π₁ is GAS.

 

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

Scenario 3: α₁ = α₂ = 0.0003, r = 0.006 and g = 0.001. Then, we calculate R₀ = 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 Π₂ = (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: α₁ = α₂ = 0.0003, r = 0.001 and g = 0.01. Then, we calculate R₀ = 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 Π₃ exists and it is GAS. Figure 4 supports the theoretical results of Theorem 3.4, where the states of the system reach the equilibrium Π₃ = (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: α₁ = α₂ = 0.0004, r = 0.01 and g = 0.01. Then, we calculate R₀ = 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, Π₄ 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 Π₄ = (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 α₁ = α₂ = 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 thatR₀ = 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 Π₀ 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.