# Stability of a Discrete-Time Pathogen Infection Model with Adaptive Immune Response

Authors: M. A. Alshaikh; A. M. Elaiw
Journal: Discrete Dynamics in Nature and Society
 (2020)

---

## Abstract

This paper studies the global stability of a discrete-time pathogen dynamic model with both cell-mediated and antibody immune responses. Both latently and actively infected cells are incorporated into the model. We discretize the continuous-time model by using the nonstandard finite difference (NSFD) method. We establish that NSFD preserves the nonnegativity and boundedness of the solutions of the model. We derive four threshold parameters which govern the existence and stability of the steady states. We establish by using the Lyapunov method, the global stability of the five steady states of the model. We illustrate our theoretical results by using numerical simulations.


---

## Body

## 1. Introduction

Studying and analysing the inward host dynamics of pathogens that infect the human body such as the human immunodeficiency virus (HIV), hepatitis C virus (HCV), hepatitis B virus (HBV), and human T-cell leukemia virus (HTLV) have received great attention from scientists and researchers (see, e.g., [1–12]). The dynamical behavior of pathogens can be understood by using mathematical modeling. The adaptive immune system plays an important role in controlling pathogenic infections. B cells and cytotoxic T lymphocytes (CTLs) are the main types of lymphocytes. B cells produce antibodies (antibody immune response) to neutralize the pathogens. CTLs (cell-mediated immune response) attack and kill cells infected by pathogen. In the literature, the impact of antibody immune response on pathogenic infections has been studied in several mathematical models (see, e.g., [13–20]). Further, the effect of cell-mediated immune response on pathogen dynamics has been modeled in several works (see, e.g., [21–27]). Wodarz [28] has combined both CTL and antibody immune responses in the following pathogen dynamic model:(1)dFdt=ω−ϑF−ϰFH,(2)dSdt=ϰFH−aS−δSZ,(3)dHdt=ηS−cH−dHY,(4)dYdt=κHY−λY,(5)dZdt=gSZ−μZ,where F=Ft, S=St, H=Ht, Y=Yt, and Z=Zt are the concentrations of healthy cells, infected cells, free pathogens, antibodies, and CTLs at time t, respectively; the healthy cells are generated, dead, and become infected at rates ω, ϑF, and ϰFH; the term aS represents the death rate of infected cells; the infected cells are killed by CTLs at a rate δSZ; the free pathogens are produced from infected cells at a rate ηS, die at a rate cH, and neutralized by the antibodies at a rate dHY; the antibodies and CTLs are proliferated at rates κHY and gSZ, respectively, while they die at rates λY and μZ, respectively.In models (1)–(5), it is assumed that cells generate mature pathogens directly after they are infected. However, it has been reported in [29] that a number of intracellular processes are needed to produce new mature pathogens. It has been shown in [30–32] that latent HIV reservoirs serve as a major barrier in HIV treatment. Pathogen dynamic models with latently infected cells and adaptive immunity have been studied in [33–35].In the previous mentioned works, the pathogen dynamics are modeled by systems of nonlinear differential equations which are difficult to be solved analytically. Therefore, we can use an approximate solution which is obtained by discretizing the continuous-time model. Nonstandard finite difference (NSFD) has become one of the famous discretization methods which has been utilized in solving different types of differential equations [36]. NSFD can preserve the basic and global properties of several types of nonlinear systems [37]. NSFD has been used to discretize different pathogen dynamic models (see, e.g., [37–40]). The effect of the cell-mediated immune response has been considered of discrete-time pathogen dynamic models in [41, 42]. In addition, Elaiw and Alshaikh [43–45] have studied discrete-time pathogen infection models with antibodies.The present paper aims to study the impact of both antibody and cell-mediated immune responses on the pathogen dynamic model. The model takes into account two types of infected cells, latently infected cells and actively infected cells. We use the NSFD method to discretize the continuous-time model. We show that the model is well posed. Further, we construct Lyapunov functions to prove the global stability of the steady states.


## 2. The Model

We consider the following continuous pathogen infection model with cell-mediated immune response, antibody immune response, and latency:(6)dFdt=ω−ϑF−ϰFH,(7)dKdt=1−εϰFH−α+mK,(8)dSdt=εϰFH+mK−aS−δSZ,(9)dHdt=ηS−cH−dHY,(10)dYdt=κHY−λY,(11)dZdt=gSZ−μZ,where K=Kt is the concentration of latently infected cells at time t; a fraction ε∈0,1 of healthy cells becomes actively infected cells, and the remaining fraction 1−ε will become latently infected cells; the terms mK and αK represent the activation and death rates of the latently infected cells, respectively.We apply the NSFD method [36] on systems (6)–(11), and we obtain(12)Fn+1−FnΠh=ω−ϑFn+1−ϰFn+1Hn,(13)Kn+1−KnΠh=1−εϰFn+1Hn−α+mKn+1,(14)Sn+1−SnΠh=εϰFn+1Hn+mKn+1−aSn+1−δSn+1Zn+1,(15)Hn+1−HnΠh=ηSn+1−cHn+1−dHn+1Yn+1,(16)Yn+1−YnΠh=κHn+1Yn+1−λYn+1,(17)Zn+1−ZnΠh=gSn+1Zn+1−μZn+1,where h>0 is the time step size and Fn,Kn,Sn,Hn,Yn,Zn are the approximations of the solution Ftn,Ktn,Stn,Htn,Ytn,Ztn of systems (6)–(11) at the discrete-time points tn=nh, n∈ℕ=0,1,2,…. The denominator function Πh is chosen such that Πh=h+Oh2 [46, 47]. Let us consider the following form of Π [47]:(18)Πh=1−e−ϑhϑ.The initial conditions of systems (12)–(17) are(19)F0,K0,S0,H0,Y0,Z0∈ℝ+6=F,K,S,H,Y,ZF>0,K>0,S>0,H>0,Y>0,Z>0.


### 2.1. Preliminaries

We define the set(20)Γ=F,K,S,H,Y,Z:0<F,K,S<δ1,0<H<δ3,0<Y<δ4,0<Z<δ2,where δ1=ω/σ, δ2=g/δδ1, δ3=2η/aδ1, δ4=κ/dδ3, and σ=minϑ,α,a/2,c,λ,μ.Lemma 1.
Any solutionFn,Kn,Sn,Hn,Yn,Zn of models (12)–(17) with initial condition (19) is nonnegative and ultimately bounded.Proof.
From equations (12)–(17), we obtain(21)Fn+1=Fn+Πhω1+Πhϑ+ϰHn,(22)Kn+1=Kn+Πh1−εϰFn+1Hn1+Πhα+m,(23)Sn+1=Sn+ΠhεϰFn+1Hn+ΠhmKn+11+Πha+δZn+1,(24)Hn+1=Hn+ΠhηSn+11+Πhc+dYn+1,(25)Yn+1=Yn1+Πhλ−κHn+1,(26)Zn+1=Zn1+Πhμ−gSn+1.
We show thatF1,K1,S1,H1,Y1,Z1 exists uniquely and is positive. If n=0, then from equations (21) and (22) and the initial condition (19), we get F1>0 and K1>0. From equations (23) and (26), we get(27)Z1=Z0+ΠhgS0+ΠhεϰF1H0+ΠhmK11+Πha+δZ1Z1−ΠhμZ1.
Then,(28)A1Z12+B1Z1+C1=0,where(29)A1=1+ΠhμΠhδ,B1=1+Πhμ1+Πha−ΠhδZ0−ΠhgS0+ΠhεϰF1H0+ΠhmK1,C1=−1+ΠhaZ0.
SinceA1>0 and C1<0, then B12−4A1C1>0, and hence equation (27) admits a unique positive root, Z1>0. It follows from equation (23) that(30)S1=S0+ΠhεϰF1H0+ΠhmK11+Πha+δZ1>0.
Next, we considerY1. It follows from equations (24) and (25) that(31)Y1=Y0+ΠhκH0+ΠhηS11+Πhc+dY1Y1−ΠhλY1.
Then, we get(32)A2Y12+B2Y1+C2=0,where(33)A2=1+ΠhλΠhd,B2=1+Πhλ1+Πhc−ΠhdY0−ΠhκH0+ΠhηS1,C2=−1+ΠhcY0.
SinceA2>0 and C2<0, then B22−4A2C2>0, and hence there exists a unique positive root of equation (32), Y1>0.
From equation (24), we get(34)H1=H0+ΠhηS11+Πhc+dY1>0.
Therefore, the solutionF1,K1,S1,H1,Y1,Z1 exists uniquely and is nonnegative.
Forn=1, repeat the above process; we can show that F2,K2,S2,H2,Y2,Z2 exists uniquely and is positive. Utilizing the mathematical induction, for all n∈ℕ, we know that Fn,Kn,Sn,Hn,Yn,Zn exists uniquely and is positive.
To investigate the boundedness of solutions, we define(35)Ωn=Fn+Kn+Sn+a2ηHn+ad2ηκYn+δgZn.
Using equations (12)–(17), we get(36)Ωn+1−Ωn=Fn+1−Fn+Kn+1−Kn+Sn+1−Sn+a2ηHn+1−Hn+ad2ηκYn+1−Yn+δgZn+1−Zn=Πhω−ϑFn+1−ϰFn+1Hn+Πh1−εϰFn+1Hn−α+mKn+1+ΠhεϰFn+1Hn+mKn+1−aSn+1−δSn+1Zn+1+Πha2ηηSn+1−cHn+1−dHn+1Yn+1+Πhad2ηκκHn+1Yn+1−λYn+1+ΠhδggSn+1Zn+1−μZn+1=Πhω−ϑFn+1−αKn+1−a2Sn+1−ac2ηHn+1−adλ2ηκYn+1−δμgZn+1.
Then,(37)Ωn+1−Ωn≤Πhω−ΠhσFn+1+Kn+1+Sn+1+a2ηHn+1+ad2ηκYn+1+δgZn+1=Πhω−ΠhσΩn+1.
Hence,(38)Ωn+1≤Ωn1+Πhσ+Πhω1+Πhσ.
Using lemma 2.2 in [48], we have(39)Ωn≤11+ΠhσnΩ0+ωσ1−11+Πhσn.
Then,limn⟶∞supΩn≤δ1, limn⟶∞supFn≤δ1, limn⟶∞supKn≤δ1, limn⟶∞supSn≤δ1, limn⟶∞supHn≤δ3, limn⟶∞supYn≤δ4, and limn⟶∞supZn≤δ2. Therefore, the solution Fn,Kn,Sn,Hn,Yn,Zn converges to Γ as n⟶∞.


### 2.2. Steady States

Systems (12)–(17) have five steady states:(i)
Pathogen-free steady stateΘ0F0,0,0,0,0,0, where F0=ω/ϑ.(ii)
Chronic pathogen steady state without immune responseΘ∗F∗,K∗,S∗,H∗,0,0, where(40)F∗=F0ℛ0,K∗=ϑac1−εϰηαε+mℛ0−1,S∗=ϑcϰηℛ0−1,H∗=ϑϰℛ0−1,in whichℛ0 is the basic reproduction number given as(41)ℛ0=ϰηωαε+mϑacα+m.Thus,Θ∗ exists if ℛ0>1.(iii)
Chronic pathogen steady state with only antibody immune responseΘ¯F¯,K¯,S¯,H¯,Y¯,0, where(42)F¯=κωκϑ+ϰλ,K¯=ωϰλ1−εα+mκϑ+ϰλ,S¯=ωϰλαε+maα+mκϑ+ϰλ,H¯=λκ,Y¯=cdℛ1Y−1,in which(43)ℛ1Y=ϰηωκαε+macα+mκϑ+ϰλ=κϑℛ0κϑ+ϰλ=ℛ01+ϰλ/κϑ.Clearly,ℛ1Y<ℛ0. If ℛ1Y>1, then Θ¯ exists.(iv)
Chronic pathogen steady state with only cell-mediated immune responseΘ^F^,K^,S^,H^,0,Z^, where(44)F^=cgωcgϑ+ηϰμ,K^=ωϰημ1−εα+mcgϑ+ηϰμ,S^=μg,H^=ημcg,Z^=aδℛ1Z−1,in which(45)ℛ1Z=gϰωηαε+maα+mcgϑ+ηϰμ.Therefore,Θ^ exists if ℛ1Z>1.(v)
Chronic pathogen steady state with both antibody and cell-mediated immune responsesΘ˜F˜,K˜,S˜,H˜,Y˜,Z˜, where(46)F˜=κωκϑ+ϰλ=F¯,K˜=ωϰλ1−εα+mκϑ+ϰλ=K¯,S˜=μg=S^,H˜=λκ=H¯,Y˜=cdℛ1Y/ℛ2Z−1,Z˜=aδℛ2Z−1,in which(47)ℛ2Z=λgϰωαε+maμα+mκϑ+ϰλ,ℛ1Yℛ2Z=ηκμcgλ.Thus,Θ˜ exists when ℛ2Z>1 and ℛ1Y/ℛ2Z>1.


### 2.3. Global Stability Analysis

In this section, we study the global stability of the steady states by using the Lyapunov method. We define a functionG as(48)Gθ=θ−1−lnθ,which yields(49)lnθ≤θ−1.Theorem 1.
Letℛ0≤1, then Θ0 of models (12)–(17) is globally asymptotically stable (GAS).Proof.
Define a functionℒnFn,Kn,Sn,Hn,Yn,Zn as(50)ℒn=1ΠhF0GFnF0+γ1Kn+γ2Sn+γ31+ΠhcHn+γ41+ΠhλYn+γ51+ΠhμZn,where(51)γ11−ε+γ2ε=1,γ1α+m=γ2m,γ2a=γ3η,γ3d=γ4κ,γ2δ=γ5g.
The solution of equation (51) is given by(52)γ1=mαε+m,γ2=α+mαε+m,γ3=aα+mηαε+m,γ4=daα+mκηαε+m,γ5=δα+mgαε+m.
Clearly,ℒn>0 for all Fn,Kn,Sn,Hn,Yn,Zn>0 and ℒn=0 at Θ0. We calculate Δℒn=ℒn+1−ℒn as(53)Δℒn=1ΠhF0GFn+1F0+γ1Kn+1+γ2Sn+1+γ31+ΠhcHn+1+γ41+ΠhλYn+1+γ51+ΠhμZn+1−1ΠhF0GFnF0+γ1Kn+γ2Sn+γ31+ΠhcHn+γ41+ΠhλYn+γ51+ΠhμZn=1ΠhF0Fn+1F0−FnF0+lnFnFn+1+γ1Kn+1−Kn+γ2Sn+1−Sn+γ31+ΠhcHn+1−Hn+γ41+ΠhλYn+1−Yn+γ51+ΠhμZn+1−Zn.
Using inequality (49), we get(54)Δℒn≤1ΠhFn+1−Fn+F0FnFn+1−1+γ1Kn+1−Kn+γ2Sn+1−Sn+γ31+ΠhcHn+1−Hn+γ41+ΠhλYn+1−Yn+γ51+ΠhμZn+1−Zn=1Πh1−F0Fn+1Fn+1−Fn+γ1Kn+1−Kn+γ2Sn+1−Sn+γ31+ΠhcHn+1−Hn+γ41+ΠhλYn+1−Yn+γ51+ΠhμZn+1−Zn.
From equations (12)–(17), we have(55)Δℒn≤1−F0Fn+1ω−ϑFn+1−ϰFn+1Hn+γ11−εϰFn+1Hn−α+mKn+1+γ2εϰFn+1Hn+mKn+1−aSn+1−δSn+1Zn+1+γ3ηSn+1−cHn+1−dHn+1Yn+1+cγ3Hn+1−Hn+γ4κHn+1Yn+1−λYn+1+λγ4Yn+1−Yn+γ5gSn+1Zn+1−μZn+1+μγ5Zn+1−Zn=1−F0Fn+1ω−ϑFn+1+ϰF0Hn−cγ3Hn−λγ4Yn−μγ5Zn=−ϑFn+1Fn+1−F02+cγ3ϰF0cγ3−1Hn−λγ4Yn−μγ5Zn=−ϑFn+1Fn+1−F02+cγ3ℛ0−1Hn−λγ4Yn−μγ5Zn.
Since,ℛ0≤1, then for all n≥0, we have Δℒn≤0. Hence, the sequence ℒn is monotonically decreasing. Since ℒn≥0, then limn⟶∞ℒn≥0, and hence limn⟶∞Δℒn=0. This yields limn⟶∞Fn=F0,limn⟶∞Yn=0, limn⟶∞Zn=0, and limn⟶∞ℛ0−1Hn=0. For the case ℛ0<1, we have limn⟶∞Fn=F0, limn⟶∞Hn=0, limn⟶∞Yn=0, and limn⟶∞Zn=0. From equations (14) and (15), we obtain limn⟶∞Sn=0 and limn⟶∞Kn=0. If ℛ0=1, we have limn⟶∞Fn=F0, limn⟶∞Yn=0, and limn⟶∞Zn=0. From equations (12), (14), and (15), we obtain limn⟶∞Hn=0, limn⟶∞Kn=0, and limn⟶∞Sn=0. Hence, the pathogen-free steady state Θ0 is GAS.Theorem 2.
Ifℛ1Y≤1<ℛ0 and ℛ1Z≤1<ℛ0, then Θ∗ is GAS.Proof.
DefineUnFn,Kn,Sn,Hn,Yn,Zn as(56)Un=1ΠhF∗GFnF∗+γ1K∗GKnK∗+γ2S∗GSnS∗+γ31+ΠhcH∗GHnH∗+γ4Yn+γ5Zn.
We haveUnFn,Kn,Sn,Hn,Yn,Zn>0 for all Fn,Kn,Sn,Hn,Yn,Zn>0, and moreover UnF∗,K∗,S∗,H∗,0,0=0. Computing ΔUn=Un+1−Un as(57)ΔUn=1ΠhF∗GFn+1F∗+γ1K∗GKn+1K∗+γ2S∗GSn+1S∗+γ31+ΠhcH∗GHn+1H∗+γ4Yn+1+γ5Zn+1−1ΠhF∗GFnF∗+γ1K∗GKnK∗+γ2S∗GSnS∗+γ31+ΠhcH∗GHnH∗+γ4Yn+γ5Zn=1ΠhF∗Fn+1F∗−FnF∗+lnFnFn+1+γ1K∗Kn+1K∗−KnK∗+lnKnKn+1+γ2S∗Sn+1S∗−SnS∗+lnSnSn+1+γ3H∗Hn+1H∗−HnH∗+lnHnHn+1+γ4Yn+1−Yn+γ5Zn+1−Zn+cγ3H∗GHn+1H∗−GHnH∗.
Using inequality (49), we get(58)ΔUn≤1ΠhFn+1−Fn+F∗FnFn+1−1+γ1Kn+1−Kn+K∗KnKn+1−1+γ2Sn+1−Sn+S∗SnSn+1−1+γ3Hn+1−Hn+H∗HnHn+1−1+γ4Yn+1−Yn+γ5Zn+1−Zn+cγ3H∗GHn+1H∗−GHnH∗=1Πh1−F∗Fn+1Fn+1−Fn+γ11−K∗Kn+1Kn+1−Kn+γ21−S∗Sn+1Sn+1−Sn+γ31−H∗Hn+1Hn+1−Hn+γ4Yn+1−Yn+γ5Zn+1−Zn+cγ3H∗Hn+1H∗−HnH∗+lnHnHn+1.
From equations (12)–(17), we have(59)ΔUn≤1−F∗Fn+1ω−ϑFn+1−ϰFn+1Hn+γ11−K∗Kn+11−εϰFn+1Hn−α+mKn+1+γ21−S∗Sn+1εϰFn+1Hn+mKn+1−aSn+1−δSn+1Zn+1+γ31−H∗Hn+1ηSn+1−cHn+1−dHn+1Yn+1+γ4κHn+1Yn+1−λYn+1+γ5gSn+1Zn+1−μZn+1+cγ3H∗Hn+1H∗−HnH∗+lnHnHn+1=1−F∗Fn+1ω−ϑFn+1+ϰF∗Hn−γ11−εϰK∗Fn+1HnKn+1+γ1α+mK∗−γ2εϰS∗Fn+1HnSn+1−γ2mS∗Kn+1Sn+1+γ2aS∗+γ2δS∗Zn+1−γ3ηH∗Sn+1Hn+1+γ3cH∗+γ3dH∗Yn+1−γ4λYn+1−γ5μZn+1−cγ3Hn+cγ3H∗lnHnHn+1.
From the steady state conditions ofΘ∗, we have(60)ω=ϑF∗+ϰF∗H∗,1−εϰF∗H∗=α+mK∗,εϰF∗H∗+mK∗=aS∗,ηS∗=cH∗,and we get(61)γ11−εϰF∗H∗=γ1α+mK∗=γ2mK∗,ϰF∗H∗=γ2aS∗=γ3ηS∗,ϰF∗=cγ3,ΔUn≤1−F∗Fn+1ϑF∗+ϰF∗H∗−ϑFn+1−γ11−εϰF∗H∗K∗Fn+1HnKn+1F∗H∗+γ11−εϰF∗H∗−γ2εϰF∗H∗S∗Fn+1HnSn+1F∗H∗−γ11−εϰF∗H∗S∗Kn+1Sn+1K∗+ϰF∗H∗+γ5gS∗−μgZn+1−ϰF∗H∗H∗Sn+1Hn+1S∗+ϰF∗H∗+γ4κH∗−λκYn+1+ϰF∗H∗lnHnHn+1=1−F∗Fn+1ϑF∗−ϑFn+1+γ11−εϰF∗H∗4−F∗Fn+1−K∗Fn+1HnKn+1F∗H∗−S∗Kn+1Sn+1K∗−H∗Sn+1Hn+1S∗+lnHnHn+1+γ2εϰF∗H∗3−F∗Fn+1−S∗Fn+1HnSn+1F∗H∗−H∗Sn+1Hn+1S∗+lnHnHn+1+γ5gS∗−S^Zn+1+γ4κH∗−H¯Yn+1=−ϑFn+1Fn+1−F∗2−γ11−εϰF∗H∗GF∗Fn+1+GK∗Fn+1HnKn+1F∗H∗+GS∗Kn+1Sn+1K∗+GH∗Sn+1Hn+1S∗−γ2εϰF∗H∗GF∗Fn+1+GS∗Fn+1HnSn+1F∗H∗+GH∗Sn+1Hn+1S∗+γ5cgϑ+ηϰμηϰℛ1Z−1Zn+1+γ4κϑ+ϰλϰℛ1Y−1Yn+1.
Sinceℛ1Z≤1 and ℛ1Y≤1, then Un is monotonically decreasing. We have Un≥0, and then limn⟶∞Un≥0 and hence limn⟶∞ΔUn=0, which implies limn⟶∞Fn=F∗, limn⟶∞Kn=K∗, limn⟶∞Sn=S∗,limn⟶∞Hn=H∗, limn⟶∞ℛ1Y−1Yn=0, and limn⟶∞ℛ1Z−1Zn=0. We have four cases:(i)
ℛ1Y=1 and ℛ1Z=1, and then from equation (14),(62)0=εϰF∗H∗+mK∗−aS∗−δS∗limn⟶∞Zn+1.Applying the steady state condition (60), we obtain limn⟶∞Zn=0. In addition, from equation (15), we get(63)0=ηS∗−cH∗−dH∗limn⟶∞Yn+1.From equation (60), we get limn⟶∞Yn=0.(ii)
ℛ1Y=1, ℛ1Z<1, and limn⟶∞Zn=0. From equation (63), we get limn⟶∞Yn=0.(iii)
ℛ1Y<1, ℛ1Z=1, and limn⟶∞Yn=0. From equation (62), we get limn⟶∞Zn=0.(iv)
ℛ1Y<1, ℛ1Z<1, limn⟶∞Yn=0, and limn⟶∞Zn=0. It follows that if ℛ1Z≤1 and ℛ1Y≤1, then limn⟶∞Fn=F∗, limn⟶∞Kn=K∗, limn⟶∞Sn=S∗,limn⟶∞Hn=H∗, limn⟶∞Yn=0, and limn⟶∞Zn=0. Then, Θ∗ is GAS.Theorem 3.
Letℛ1Y>1 and ℛ2Z≤1, then Θ¯ of models (12)–(17) is GAS.Proof.
DefineWnFn,Kn,Sn,Hn,Yn,Zn as(64)Wn=1ΠhF¯GFnF¯+γ1K¯GKnK¯+γ2S¯GSnS¯+γ31+Πhc+dY¯H¯GHnH¯+γ4Y¯GYnY¯+γ5Zn.
Clearly,WnFn,Kn,Sn,Hn,Yn,Zn>0 for all Fn,Kn,Sn,Hn,Yn,Zn>0 and WnF¯,K¯,S¯,H¯,Y¯,0=0. Computing ΔWn=Wn+1−Wn as(65)ΔWn=1ΠhF¯GFn+1F¯+γ1K¯GKn+1K¯+γ2S¯GSn+1S¯+γ31+Πhc+dY¯H¯GHn+1H¯+γ4Y¯GYn+1Y¯+γ5Zn+1−1ΠhF¯GFnF¯+γ1K¯GKnK¯+γ2S¯GSnS¯+γ31+Πhc+dY¯H¯GHnH¯+γ4Y¯GYnY¯+γ5Zn=1ΠhF¯Fn+1F¯−FnF¯+lnFnFn+1+γ1K¯Kn+1K¯−KnK¯+lnKnKn+1+γ2S¯Sn+1S¯−SnS¯+lnSnSn+1+γ3H¯Hn+1H¯−HnH¯+lnHnHn+1+γ4Y¯Yn+1Y¯−YnY¯+lnYnYn+1+γ5Zn+1−Zn+γ3c+dY¯H¯GHn+1H¯−GHnH¯.
Using inequality (49), we get(66)ΔWn≤1ΠhF¯Fn+1−FnF¯+FnFn+1−1+γ1K¯Kn+1−KnK¯+KnKn+1−1+γ2S¯Sn+1−SnS¯+SnSn+1−1+γ3H¯Hn+1−HnH¯+HnHn+1−1+γ4Y¯Yn+1−YnY¯+YnYn+1−1+γ5Zn+1−Zn+γ3c+dY¯H¯GHn+1H¯−GHnH¯=1Πh1−F¯Fn+1Fn+1−Fn+γ11−K¯Kn+1Kn+1−Kn+γ21−S¯Sn+1Sn+1−Sn+γ31−H¯Hn+1Hn+1−Hn+γ41−Y¯Yn+1Yn+1−Yn+γ5Zn+1−Zn+γ3c+dY¯H¯Hn+1H¯−HnH¯+lnHnHn+1.
From equations (12)–(17), we have(67)ΔWn≤1−F¯Fn+1ω−ϑFn+1−ϰFn+1Hn+γ11−K¯Kn+11−εϰFn+1Hn−α+mKn+1+γ21−S¯Sn+1εϰFn+1Hn+mKn+1−aSn+1−δSn+1Zn+1+γ31−H¯Hn+1ηSn+1−cHn+1−dHn+1Yn+1+γ41−Y¯Yn+1κHn+1Yn+1−λYn+1+γ5gSn+1Zn+1−μZn+1+γ3c+dY¯H¯Hn+1H¯−HnH¯+lnHnHn+1=1−F¯Fn+1ω−ϑFn+1+ϰF¯Hn−γ11−εϰK¯Fn+1HnKn+1+γ1α+mK¯−γ2εϰS¯Fn+1HnSn+1−γ2mS¯Kn+1Sn+1+γ2aS¯+γ2δS¯Zn+1−γ3ηH¯Sn+1Hn+1+γ3cH¯+γ3dH¯Yn+1−γ4λYn+1+γ4λY¯−γ5μZn+1−γ3c+dY¯Hn+γ3c+dY¯H¯lnHnHn+1.
We have(68)ω=ϑF¯+ϰF¯H¯,1−εϰF¯H¯=α+mK¯,εϰF¯H¯+mK¯=aS¯,ηS¯=cH¯+dH¯Y¯,λ=κH¯,and we get(69)γ11−εϰF¯H¯=γ1α+mK¯=γ2mK¯,ϰF¯H¯=γ2aS¯=γ3ηS¯,ϰF¯=γ3c+dY¯,ΔWn≤1−F¯Fn+1ϑF¯+ϰF¯H¯−ϑFn+1−γ11−εϰF¯H¯Fn+1HnK¯F¯H¯Kn+1+γ11−εϰF¯H¯−γ2εϰF¯H¯S¯Fn+1HnSn+1F¯H¯−γ11−εϰF¯H¯S¯Kn+1Sn+1K¯+ϰF¯H¯+γ5gS¯−μgZn+1−ϰF¯H¯H¯Sn+1Hn+1S¯+ϰF¯H¯+ϰF¯H¯lnHnHn+1=−ϑFn+1Fn+1−F¯2+γ11−εϰF¯H¯4−F¯Fn+1−Fn+1HnK¯F¯H¯Kn+1−S¯Kn+1Sn+1K¯−H¯Sn+1Hn+1S¯+lnHnHn+1+γ2εϰF¯H¯3−F¯Fn+1−S¯Fn+1HnSn+1F¯H¯−H¯Sn+1Hn+1S¯+lnHnHn+1+γ5gS¯−S^Zn+1=−ϑFn+1Fn+1−F¯2−γ11−εϰF¯H¯GF¯Fn+1+GFn+1HnK¯F¯H¯Kn+1+GS¯Kn+1Sn+1K¯+GH¯Sn+1Hn+1S¯−γ2εϰF¯H¯GF¯Fn+1+GS¯Fn+1HnSn+1F¯H¯+GH¯Sn+1Hn+1S¯+γ5μℛ2Z−1Zn+1.
Sinceℛ2Z≤1, then ΔWn≤0. Then, the sequence Wn is monotonically decreasing. Since Wn≥0, then limn⟶∞Wn≥0 and hence limn⟶∞ΔWn=0. Thus, limn⟶∞Fn=F¯, limn⟶∞Kn=K¯, limn⟶∞Sn=S¯, limn⟶∞Hn=H¯, and limn⟶∞ℛ2Z−1Zn=0. We have two cases:(i)
ℛ2Z=1, and then from equation (14),(70)0=εϰF¯H¯+mK¯−aS¯−δS¯limn⟶∞Zn+1,and we getlimn⟶∞Zn=0. From equation (15), we get(71)0=ηS¯−cH¯−dH¯limn⟶∞Yn+1.This giveslimn⟶∞Yn=Y¯.(ii)
ℛ2Z<1 and limn⟶∞Zn=0. From equation (71), we get limn⟶∞Yn=Y¯. Hence, Θ¯ is GAS.Theorem 4.
Letℛ1Z>1 and ℛ1Y/ℛ2Z≤1, then Θ^ of models (12)–(17) is GAS.Proof.
DefineℳnFn,Kn,Sn,Hn,Yn,Zn as(72)ℳn=1ΠhF^GFnF^+γ1K^GKnK^+γ2S^GSnS^+γ3a+δZ^a1+ΠhcH^GHnH^+γ4a+δZ^aYn+γ5Z^GZnZ^.
Clearly,ℳnFn,Kn,Sn,Hn,Yn,Zn>0 for all Fn,Kn,Sn,Hn,Yn,Zn>0 and ℳnF^,K^,S^,H^,0,Z^=0. We have Δℳn=ℳn+1−ℳn as(73)Δℳn=1ΠhF^GFn+1F^+γ1K^GKn+1K^+γ2S^GSn+1S^+γ3a+δZ^a1+ΠhcH^GHn+1H^+γ4a+δZ^aYn+1+γ5Z^GZn+1Z^−1ΠhF^GFnF^+γ1K^GKnK^+γ2S^GSnS^+γ3a+δZ^a1+ΠhcH^GHnH^+γ4a+δZ^aYn+γ5Z^GZnZ^=1ΠhF^Fn+1F^−FnF^+lnFnFn+1+γ1K^Kn+1K^−KnK^+lnKnKn+1+γ2S^Sn+1S^−SnS^+lnSnSn+1+γ3a+δZ^aH^Hn+1H^−HnH^+lnHnHn+1+γ4a+δZ^aYn+1−Yn+γ5Z^Zn+1Z^−ZnZ^+lnZnZn+1+γ3ca+δZ^aGHn+1H^−GHnH^.
Using inequality (49), we get(74)Δℳn≤1ΠhF^Fn+1−FnF^+FnFn+1−1+γ1K^Kn+1−KnK^+KnKn+1−1+γ2S^Sn+1−SnS^+SnSn+1−1+γ3a+δZ^aH^Hn+1−HnH^+HnHn+1−1+γ4a+δZ^aYn+1−Yn+γ5Z^Zn+1Z^−ZnZ^+ZnZn+1−1+γ3ca+δZ^aH^GHn+1H^−GHnH^=1Πh1−F^Fn+1Fn+1−Fn+γ11−K^Kn+1Kn+1−Kn+γ21−S^Sn+1Sn+1−Sn+γ3a+δZ^a1−H^Hn+1Hn+1−Hn+γ4a+δZ^aYn+1−Yn+γ51−Z^Zn+1Zn+1−Zn+γ3ca+δZ^aH^Hn+1H^−HnH^+lnHnHn+1.
From equations (12)–(17), we have(75)Δℳn≤1−F^Fn+1ω−ϑFn+1−ϰFn+1Hn+γ11−K^Kn+11−εϰFn+1Hn−α+mKn+1+γ21−S^Sn+1εϰFn+1Hn+mKn+1−aSn+1−δSn+1Zn+1+γ3a+δZ^a1−H^Hn+1ηSn+1−cHn+1−dHn+1Yn+1+γ4a+δZ^aκHn+1Yn+1−λYn+1+γ51−Z^Zn+1gSn+1Zn+1−μZn+1+γ3ca+δZ^aH^Hn+1H^−HnH^+lnHnHn+1=1−F^Fn+1ω−ϑFn+1+ϰF^Hn−γ11−εϰK^Fn+1HnKn+1+γ1α+mK^−γ2εϰS^Fn+1HnSn+1−γ2mS^Kn+1Sn+1−γ2aSn+1+γ2aS^+γ2δS^Zn+1+γ3a+δZ^aηSn+1−γ3a+δZ^ηaH^Sn+1Hn+1+γ3a+δZ^acH^+γ3a+δZ^adH^Yn+1−γ4a+δZ^aλYn+1−γ5gZ^Sn+1−γ5μZn+1+γ5μZ^−γ3ca+δZ^aHn+γ3ca+δZ^aH^lnHnHn+1.
We have(76)ω=ϑF^+ϰF^H^,1−εϰF^H^=α+mK^,εϰF^H^+mK^=a+δZ^S^,ηS^=cH^,S^=μg,and we get(77)γ11−εϰF^H^=γ1α+mK^=γ2mK^,ϰF^H^=γ2a+δZ^S^=γ3a+δZ^ηaS^,Δℳn≤1−F^Fn+1ϑF^+ϰF^H^−ϑFn+1−γ11−εϰF^H^Fn+1HnK^F^H^Kn+1+γ11−εϰF^H^−γ2εϰF^H^S^Fn+1HnSn+1F^H^−γ11−εϰF^H^S^Kn+1Sn+1K^+ϰF^H^−ϰF^H^H^Sn+1Hn+1S^+ϰF^H^+γ3da+δZ^aH^−λκYn+1+ϰF^H^lnHnHn+1=−ϑFn+1Fn+1−F^2+γ11−εϰF^H^4−F^Fn+1−Fn+1HnK^F^H^Kn+1−S^Kn+1Sn+1K^−H^Sn+1Hn+1S^+lnHnHn+1+γ2εϰF^H^3−F^Fn+1−S^Fn+1HnSn+1F^H^−H^Sn+1Hn+1S^+lnHnHn+1+γ3da+δZ^aH^−H˜Yn+1=−ϑFn+1Fn+1−F^2−γ11−εϰF^H^GF^Fn+1+GFn+1HnK^F^H^Kn+1+GS^Kn+1Sn+1K^+GH^Sn+1Hn+1S^−γ2εϰF^H^GF^Fn+1+GS^Fn+1HnSn+1F^H^+GH^Sn+1Hn+1S^+γ2a+δZ^dλκηℛ1Yℛ2Z−1Yn+1.
Sinceℛ1Z>1 and if ℛ1Y/ℛ2Z≤1, then Δℳn≤0. Hence, the sequence ℳn is monotonically decreasing. Since ℳn≥0, then limn⟶∞ℳn≥0 and hence limn⟶∞Δℳn=0. Thus, limn⟶∞Fn=F^, limn⟶∞Kn=K^, limn⟶∞Sn=S^, limn⟶∞Hn=H^, and limn⟶∞ℛ1Y/ℛ2Z−1Yn+1=0. We have two cases:(i)
ℛ1Y/ℛ2Z=1, and from equation (14),(78)0=εϰF^H^+mK^−aS^−δS^limn⟶∞Zn+1,and this giveslimn⟶∞Zn=Z^; moreover, from equation (15), we have(79)0=ηS^−cH^−dH^limn⟶∞Yn+1,then we getlimn⟶∞Yn=0.(ii)
ℛ1Y/ℛ2Z<1 and limn⟶∞Yn=0. From equation (78), we get limn⟶∞Zn=Z^ and hence Θ^ is GAS.Theorem 5.
Letℛ1Y>ℛ2Z>1, then Θ˜ of models (12)–(17) is GAS.Proof.
DefineVnFn,Kn,Sn,Hn,Yn,Zn as(80)Vn=1ΠhF˜GFnF˜+γ1K˜GKnK˜+γ2S˜GSnS˜+γ3a+δZ˜a1+Πhc+dY˜H˜GHnH˜+γ4a+δZ˜aY˜GYnY˜+γ5Z˜GZnZ˜.
Clearly,VnFn,Kn,Sn,Hn,Yn,Zn>0 for all Fn,Kn,Sn,Hn,Yn,Zn>0 and VnF˜,K˜,S˜,H˜,Y˜,Z˜=0. Computing ΔVn=Vn+1−Vn as(81)ΔVn=1ΠhF˜GFn+1F˜+γ1K˜GKn+1K˜+γ2S˜GSn+1S˜+γ3a+δZ˜a1+Πhc+dY˜H˜GHn+1H˜+γ4a+δZ˜aY˜GYn+1Y˜+γ5Z˜GZn+1Z˜−1ΠhF˜GFnF˜+γ1K˜GKnK˜+γ2S˜GSnS˜+γ3a+δZ˜a1+Πhc+dY˜H˜GHnH˜+γ4a+δZ˜aY˜GYnY˜+γ5Z˜GZnZ˜=1ΠhF˜Fn+1F˜−FnF˜+lnFnFn+1+γ1K˜Kn+1K˜−KnK˜+lnKnKn+1+γ2S˜Sn+1S˜−SnS˜+lnSnSn+1+γ3a+δZ˜aH˜Hn+1H˜−HnH˜+lnHnHn+1+γ4a+δZ˜aY˜Yn+1Y˜−YnY˜+lnYnYn+1+γ5Z˜Zn+1Z˜−ZnZ˜+lnZnZn+1+γ3c+dY˜a+δZ˜aH˜GHn+1H˜−GHnH˜.
Using inequality (49), we get(82)ΔVn≤1ΠhF˜Fn+1−FnF˜+FnFn+1−1+γ1K˜Kn+1−KnK˜+KnKn+1−1+γ2S˜Sn+1−SnS˜+SnSn+1−1+γ3a+δZ˜aH˜Hn+1−HnH˜+HnHn+1−1+γ4a+δZ˜aY˜Yn+1−YnY˜+YnYn+1−1+γ5Z˜Zn+1−ZnZ˜+ZnZn+1−1+γ3c+dY˜a+δZ˜aH˜GHn+1H˜−GHnH˜=1Πh1−F˜Fn+1Fn+1−Fn+γ11−K˜Kn+1Kn+1−Kn+γ21−S˜Sn+1Sn+1−Sn+γ3a+δZ˜a1−H˜Hn+1Hn+1−Hn+γ4a+δZ˜a1−Y˜Yn+1Yn+1−Yn+γ51−Z˜Zn+1Zn+1−Zn+γ3c+dY˜a+δZ˜aH˜Hn+1H˜−HnH˜+lnHnHn+1.
From equations (12)–(17), we have(83)ΔVn≤1−F˜Fn+1ω−ϑFn+1−ϰFn+1Hn+γ11−K˜Kn+11−εϰFn+1Hn−α+mKn+1+γ21−S˜Sn+1εϰFn+1Hn+mKn+1−aSn+1−δSn+1Zn+1+γ3a+δZ˜a1−H˜Hn+1ηSn+1−cHn+1−dHn+1Yn+1+γ4a+δZ˜a1−Y˜Yn+1κHn+1Yn+1−λYn+1+γ51−Z˜Zn+1gSn+1Zn+1−μZn+1+γ3c+dY˜a+δZ˜H˜aHn+1H˜−HnH˜+lnHnHn+1=1−F˜Fn+1ω−ϑFn+1+ϰF˜Hn−γ11−εϰK˜Fn+1HnKn+1+γ1α+mK˜−γ2εϰS˜Fn+1HnSn+1−mγ2S˜Kn+1Sn+1−aγ2Sn+1+aγ2S˜+δγ2S˜Zn+1+γ3a+δZ˜ηaSn+1−γ3a+δZ˜ηaH˜Sn+1Hn+1+γ3a+δZ˜acH˜+γ3da+δZ˜aH˜Yn+1−γ4a+δZ˜aλYn+1+γ4a+δZ˜aλY˜−γ5μZn+1−γ5gZ˜Sn+1+γ5μZ˜−γ3c+dY˜a+δZ˜aHn+γ3c+dY˜a+δZ˜aH˜lnHnHn+1.
We have(84)ω=ϑF˜+ϰF˜H˜,1−εϰF˜H˜=α+mK˜,εϰF˜H˜+mK˜=a+δZ˜S˜,ηS˜=c+dY˜H˜,S˜=μg,H˜=λκ,and we get(85)γ11−εϰF˜H˜=γ1α+mK˜=γ2mK˜,ϰF˜H˜=γ2a+δZ˜S˜=γ3a+δZ˜ηaS˜,ΔVn≤1−F˜Fn+1ϑF˜+ϰF˜H˜−ϑFn+1−γ11−εϰF˜H˜Fn+1HnK˜F˜H˜Kn+1+γ11−εϰF˜H˜−γ2εϰF˜H˜S˜Fn+1HnSn+1F˜H˜−γ11−εϰF˜H˜S˜Kn+1Sn+1K˜+ϰF˜H˜−ϰF˜H˜H˜Sn+1Hn+1S˜+ϰF˜H˜+ϰF˜H˜lnHnHn+1=−ϑFn+1Fn+1−F˜2+γ11−εϰF˜H˜4−F˜Fn+1−Fn+1HnK˜F˜H˜Kn+1−S˜Kn+1Sn+1K˜−H˜Sn+1Hn+1S˜+lnHnHn+1+γ2εϰF˜H˜3−F˜Fn+1−S˜Fn+1HnSn+1F˜H˜−H˜Sn+1Hn+1S˜+lnHnHn+1=−ϑFn+1Fn+1−F˜2−γ11−εϰF˜H˜GF˜Fn+1+GFn+1HnK˜F˜H˜Kn+1+GS˜Kn+1Sn+1K˜+GH˜Sn+1Hn+1S˜−γ2εϰF˜H˜GF˜Fn+1+GS˜Fn+1HnSn+1F˜H˜+GH˜Sn+1Hn+1S˜.
We note thatΔVn≤0. It follows that Vn is a monotone decreasing sequence. Since Vn≥0, then limn⟶∞Vn≥0 and hence limn⟶∞ΔVn=0. Thus, limn⟶∞Fn=F˜, limn⟶∞Kn=K˜, limn⟶∞Sn=S˜, and limn⟶∞Hn=H˜. From equations (14) and (15), we have(86)0=εϰF˜H˜+mK˜−aS˜−δS˜limn⟶∞Zn+1,0=ηS˜−cH˜−dH˜limn⟶∞Yn+1,then limn⟶∞Yn=Y˜ and limn⟶∞Zn=Z˜. Hence, Θ˜ is GAS.


## 2.1. Preliminaries

We define the set(20)Γ=F,K,S,H,Y,Z:0<F,K,S<δ1,0<H<δ3,0<Y<δ4,0<Z<δ2,where δ1=ω/σ, δ2=g/δδ1, δ3=2η/aδ1, δ4=κ/dδ3, and σ=minϑ,α,a/2,c,λ,μ.Lemma 1.
Any solutionFn,Kn,Sn,Hn,Yn,Zn of models (12)–(17) with initial condition (19) is nonnegative and ultimately bounded.Proof.
From equations (12)–(17), we obtain(21)Fn+1=Fn+Πhω1+Πhϑ+ϰHn,(22)Kn+1=Kn+Πh1−εϰFn+1Hn1+Πhα+m,(23)Sn+1=Sn+ΠhεϰFn+1Hn+ΠhmKn+11+Πha+δZn+1,(24)Hn+1=Hn+ΠhηSn+11+Πhc+dYn+1,(25)Yn+1=Yn1+Πhλ−κHn+1,(26)Zn+1=Zn1+Πhμ−gSn+1.
We show thatF1,K1,S1,H1,Y1,Z1 exists uniquely and is positive. If n=0, then from equations (21) and (22) and the initial condition (19), we get F1>0 and K1>0. From equations (23) and (26), we get(27)Z1=Z0+ΠhgS0+ΠhεϰF1H0+ΠhmK11+Πha+δZ1Z1−ΠhμZ1.
Then,(28)A1Z12+B1Z1+C1=0,where(29)A1=1+ΠhμΠhδ,B1=1+Πhμ1+Πha−ΠhδZ0−ΠhgS0+ΠhεϰF1H0+ΠhmK1,C1=−1+ΠhaZ0.
SinceA1>0 and C1<0, then B12−4A1C1>0, and hence equation (27) admits a unique positive root, Z1>0. It follows from equation (23) that(30)S1=S0+ΠhεϰF1H0+ΠhmK11+Πha+δZ1>0.
Next, we considerY1. It follows from equations (24) and (25) that(31)Y1=Y0+ΠhκH0+ΠhηS11+Πhc+dY1Y1−ΠhλY1.
Then, we get(32)A2Y12+B2Y1+C2=0,where(33)A2=1+ΠhλΠhd,B2=1+Πhλ1+Πhc−ΠhdY0−ΠhκH0+ΠhηS1,C2=−1+ΠhcY0.
SinceA2>0 and C2<0, then B22−4A2C2>0, and hence there exists a unique positive root of equation (32), Y1>0.
From equation (24), we get(34)H1=H0+ΠhηS11+Πhc+dY1>0.
Therefore, the solutionF1,K1,S1,H1,Y1,Z1 exists uniquely and is nonnegative.
Forn=1, repeat the above process; we can show that F2,K2,S2,H2,Y2,Z2 exists uniquely and is positive. Utilizing the mathematical induction, for all n∈ℕ, we know that Fn,Kn,Sn,Hn,Yn,Zn exists uniquely and is positive.
To investigate the boundedness of solutions, we define(35)Ωn=Fn+Kn+Sn+a2ηHn+ad2ηκYn+δgZn.
Using equations (12)–(17), we get(36)Ωn+1−Ωn=Fn+1−Fn+Kn+1−Kn+Sn+1−Sn+a2ηHn+1−Hn+ad2ηκYn+1−Yn+δgZn+1−Zn=Πhω−ϑFn+1−ϰFn+1Hn+Πh1−εϰFn+1Hn−α+mKn+1+ΠhεϰFn+1Hn+mKn+1−aSn+1−δSn+1Zn+1+Πha2ηηSn+1−cHn+1−dHn+1Yn+1+Πhad2ηκκHn+1Yn+1−λYn+1+ΠhδggSn+1Zn+1−μZn+1=Πhω−ϑFn+1−αKn+1−a2Sn+1−ac2ηHn+1−adλ2ηκYn+1−δμgZn+1.
Then,(37)Ωn+1−Ωn≤Πhω−ΠhσFn+1+Kn+1+Sn+1+a2ηHn+1+ad2ηκYn+1+δgZn+1=Πhω−ΠhσΩn+1.
Hence,(38)Ωn+1≤Ωn1+Πhσ+Πhω1+Πhσ.
Using lemma 2.2 in [48], we have(39)Ωn≤11+ΠhσnΩ0+ωσ1−11+Πhσn.
Then,limn⟶∞supΩn≤δ1, limn⟶∞supFn≤δ1, limn⟶∞supKn≤δ1, limn⟶∞supSn≤δ1, limn⟶∞supHn≤δ3, limn⟶∞supYn≤δ4, and limn⟶∞supZn≤δ2. Therefore, the solution Fn,Kn,Sn,Hn,Yn,Zn converges to Γ as n⟶∞.


## 2.2. Steady States

Systems (12)–(17) have five steady states:(i)
Pathogen-free steady stateΘ0F0,0,0,0,0,0, where F0=ω/ϑ.(ii)
Chronic pathogen steady state without immune responseΘ∗F∗,K∗,S∗,H∗,0,0, where(40)F∗=F0ℛ0,K∗=ϑac1−εϰηαε+mℛ0−1,S∗=ϑcϰηℛ0−1,H∗=ϑϰℛ0−1,in whichℛ0 is the basic reproduction number given as(41)ℛ0=ϰηωαε+mϑacα+m.Thus,Θ∗ exists if ℛ0>1.(iii)
Chronic pathogen steady state with only antibody immune responseΘ¯F¯,K¯,S¯,H¯,Y¯,0, where(42)F¯=κωκϑ+ϰλ,K¯=ωϰλ1−εα+mκϑ+ϰλ,S¯=ωϰλαε+maα+mκϑ+ϰλ,H¯=λκ,Y¯=cdℛ1Y−1,in which(43)ℛ1Y=ϰηωκαε+macα+mκϑ+ϰλ=κϑℛ0κϑ+ϰλ=ℛ01+ϰλ/κϑ.Clearly,ℛ1Y<ℛ0. If ℛ1Y>1, then Θ¯ exists.(iv)
Chronic pathogen steady state with only cell-mediated immune responseΘ^F^,K^,S^,H^,0,Z^, where(44)F^=cgωcgϑ+ηϰμ,K^=ωϰημ1−εα+mcgϑ+ηϰμ,S^=μg,H^=ημcg,Z^=aδℛ1Z−1,in which(45)ℛ1Z=gϰωηαε+maα+mcgϑ+ηϰμ.Therefore,Θ^ exists if ℛ1Z>1.(v)
Chronic pathogen steady state with both antibody and cell-mediated immune responsesΘ˜F˜,K˜,S˜,H˜,Y˜,Z˜, where(46)F˜=κωκϑ+ϰλ=F¯,K˜=ωϰλ1−εα+mκϑ+ϰλ=K¯,S˜=μg=S^,H˜=λκ=H¯,Y˜=cdℛ1Y/ℛ2Z−1,Z˜=aδℛ2Z−1,in which(47)ℛ2Z=λgϰωαε+maμα+mκϑ+ϰλ,ℛ1Yℛ2Z=ηκμcgλ.Thus,Θ˜ exists when ℛ2Z>1 and ℛ1Y/ℛ2Z>1.


## 2.3. Global Stability Analysis

In this section, we study the global stability of the steady states by using the Lyapunov method. We define a functionG as(48)Gθ=θ−1−lnθ,which yields(49)lnθ≤θ−1.Theorem 1.
Letℛ0≤1, then Θ0 of models (12)–(17) is globally asymptotically stable (GAS).Proof.
Define a functionℒnFn,Kn,Sn,Hn,Yn,Zn as(50)ℒn=1ΠhF0GFnF0+γ1Kn+γ2Sn+γ31+ΠhcHn+γ41+ΠhλYn+γ51+ΠhμZn,where(51)γ11−ε+γ2ε=1,γ1α+m=γ2m,γ2a=γ3η,γ3d=γ4κ,γ2δ=γ5g.
The solution of equation (51) is given by(52)γ1=mαε+m,γ2=α+mαε+m,γ3=aα+mηαε+m,γ4=daα+mκηαε+m,γ5=δα+mgαε+m.
Clearly,ℒn>0 for all Fn,Kn,Sn,Hn,Yn,Zn>0 and ℒn=0 at Θ0. We calculate Δℒn=ℒn+1−ℒn as(53)Δℒn=1ΠhF0GFn+1F0+γ1Kn+1+γ2Sn+1+γ31+ΠhcHn+1+γ41+ΠhλYn+1+γ51+ΠhμZn+1−1ΠhF0GFnF0+γ1Kn+γ2Sn+γ31+ΠhcHn+γ41+ΠhλYn+γ51+ΠhμZn=1ΠhF0Fn+1F0−FnF0+lnFnFn+1+γ1Kn+1−Kn+γ2Sn+1−Sn+γ31+ΠhcHn+1−Hn+γ41+ΠhλYn+1−Yn+γ51+ΠhμZn+1−Zn.
Using inequality (49), we get(54)Δℒn≤1ΠhFn+1−Fn+F0FnFn+1−1+γ1Kn+1−Kn+γ2Sn+1−Sn+γ31+ΠhcHn+1−Hn+γ41+ΠhλYn+1−Yn+γ51+ΠhμZn+1−Zn=1Πh1−F0Fn+1Fn+1−Fn+γ1Kn+1−Kn+γ2Sn+1−Sn+γ31+ΠhcHn+1−Hn+γ41+ΠhλYn+1−Yn+γ51+ΠhμZn+1−Zn.
From equations (12)–(17), we have(55)Δℒn≤1−F0Fn+1ω−ϑFn+1−ϰFn+1Hn+γ11−εϰFn+1Hn−α+mKn+1+γ2εϰFn+1Hn+mKn+1−aSn+1−δSn+1Zn+1+γ3ηSn+1−cHn+1−dHn+1Yn+1+cγ3Hn+1−Hn+γ4κHn+1Yn+1−λYn+1+λγ4Yn+1−Yn+γ5gSn+1Zn+1−μZn+1+μγ5Zn+1−Zn=1−F0Fn+1ω−ϑFn+1+ϰF0Hn−cγ3Hn−λγ4Yn−μγ5Zn=−ϑFn+1Fn+1−F02+cγ3ϰF0cγ3−1Hn−λγ4Yn−μγ5Zn=−ϑFn+1Fn+1−F02+cγ3ℛ0−1Hn−λγ4Yn−μγ5Zn.
Since,ℛ0≤1, then for all n≥0, we have Δℒn≤0. Hence, the sequence ℒn is monotonically decreasing. Since ℒn≥0, then limn⟶∞ℒn≥0, and hence limn⟶∞Δℒn=0. This yields limn⟶∞Fn=F0,limn⟶∞Yn=0, limn⟶∞Zn=0, and limn⟶∞ℛ0−1Hn=0. For the case ℛ0<1, we have limn⟶∞Fn=F0, limn⟶∞Hn=0, limn⟶∞Yn=0, and limn⟶∞Zn=0. From equations (14) and (15), we obtain limn⟶∞Sn=0 and limn⟶∞Kn=0. If ℛ0=1, we have limn⟶∞Fn=F0, limn⟶∞Yn=0, and limn⟶∞Zn=0. From equations (12), (14), and (15), we obtain limn⟶∞Hn=0, limn⟶∞Kn=0, and limn⟶∞Sn=0. Hence, the pathogen-free steady state Θ0 is GAS.Theorem 2.
Ifℛ1Y≤1<ℛ0 and ℛ1Z≤1<ℛ0, then Θ∗ is GAS.Proof.
DefineUnFn,Kn,Sn,Hn,Yn,Zn as(56)Un=1ΠhF∗GFnF∗+γ1K∗GKnK∗+γ2S∗GSnS∗+γ31+ΠhcH∗GHnH∗+γ4Yn+γ5Zn.
We haveUnFn,Kn,Sn,Hn,Yn,Zn>0 for all Fn,Kn,Sn,Hn,Yn,Zn>0, and moreover UnF∗,K∗,S∗,H∗,0,0=0. Computing ΔUn=Un+1−Un as(57)ΔUn=1ΠhF∗GFn+1F∗+γ1K∗GKn+1K∗+γ2S∗GSn+1S∗+γ31+ΠhcH∗GHn+1H∗+γ4Yn+1+γ5Zn+1−1ΠhF∗GFnF∗+γ1K∗GKnK∗+γ2S∗GSnS∗+γ31+ΠhcH∗GHnH∗+γ4Yn+γ5Zn=1ΠhF∗Fn+1F∗−FnF∗+lnFnFn+1+γ1K∗Kn+1K∗−KnK∗+lnKnKn+1+γ2S∗Sn+1S∗−SnS∗+lnSnSn+1+γ3H∗Hn+1H∗−HnH∗+lnHnHn+1+γ4Yn+1−Yn+γ5Zn+1−Zn+cγ3H∗GHn+1H∗−GHnH∗.
Using inequality (49), we get(58)ΔUn≤1ΠhFn+1−Fn+F∗FnFn+1−1+γ1Kn+1−Kn+K∗KnKn+1−1+γ2Sn+1−Sn+S∗SnSn+1−1+γ3Hn+1−Hn+H∗HnHn+1−1+γ4Yn+1−Yn+γ5Zn+1−Zn+cγ3H∗GHn+1H∗−GHnH∗=1Πh1−F∗Fn+1Fn+1−Fn+γ11−K∗Kn+1Kn+1−Kn+γ21−S∗Sn+1Sn+1−Sn+γ31−H∗Hn+1Hn+1−Hn+γ4Yn+1−Yn+γ5Zn+1−Zn+cγ3H∗Hn+1H∗−HnH∗+lnHnHn+1.
From equations (12)–(17), we have(59)ΔUn≤1−F∗Fn+1ω−ϑFn+1−ϰFn+1Hn+γ11−K∗Kn+11−εϰFn+1Hn−α+mKn+1+γ21−S∗Sn+1εϰFn+1Hn+mKn+1−aSn+1−δSn+1Zn+1+γ31−H∗Hn+1ηSn+1−cHn+1−dHn+1Yn+1+γ4κHn+1Yn+1−λYn+1+γ5gSn+1Zn+1−μZn+1+cγ3H∗Hn+1H∗−HnH∗+lnHnHn+1=1−F∗Fn+1ω−ϑFn+1+ϰF∗Hn−γ11−εϰK∗Fn+1HnKn+1+γ1α+mK∗−γ2εϰS∗Fn+1HnSn+1−γ2mS∗Kn+1Sn+1+γ2aS∗+γ2δS∗Zn+1−γ3ηH∗Sn+1Hn+1+γ3cH∗+γ3dH∗Yn+1−γ4λYn+1−γ5μZn+1−cγ3Hn+cγ3H∗lnHnHn+1.
From the steady state conditions ofΘ∗, we have(60)ω=ϑF∗+ϰF∗H∗,1−εϰF∗H∗=α+mK∗,εϰF∗H∗+mK∗=aS∗,ηS∗=cH∗,and we get(61)γ11−εϰF∗H∗=γ1α+mK∗=γ2mK∗,ϰF∗H∗=γ2aS∗=γ3ηS∗,ϰF∗=cγ3,ΔUn≤1−F∗Fn+1ϑF∗+ϰF∗H∗−ϑFn+1−γ11−εϰF∗H∗K∗Fn+1HnKn+1F∗H∗+γ11−εϰF∗H∗−γ2εϰF∗H∗S∗Fn+1HnSn+1F∗H∗−γ11−εϰF∗H∗S∗Kn+1Sn+1K∗+ϰF∗H∗+γ5gS∗−μgZn+1−ϰF∗H∗H∗Sn+1Hn+1S∗+ϰF∗H∗+γ4κH∗−λκYn+1+ϰF∗H∗lnHnHn+1=1−F∗Fn+1ϑF∗−ϑFn+1+γ11−εϰF∗H∗4−F∗Fn+1−K∗Fn+1HnKn+1F∗H∗−S∗Kn+1Sn+1K∗−H∗Sn+1Hn+1S∗+lnHnHn+1+γ2εϰF∗H∗3−F∗Fn+1−S∗Fn+1HnSn+1F∗H∗−H∗Sn+1Hn+1S∗+lnHnHn+1+γ5gS∗−S^Zn+1+γ4κH∗−H¯Yn+1=−ϑFn+1Fn+1−F∗2−γ11−εϰF∗H∗GF∗Fn+1+GK∗Fn+1HnKn+1F∗H∗+GS∗Kn+1Sn+1K∗+GH∗Sn+1Hn+1S∗−γ2εϰF∗H∗GF∗Fn+1+GS∗Fn+1HnSn+1F∗H∗+GH∗Sn+1Hn+1S∗+γ5cgϑ+ηϰμηϰℛ1Z−1Zn+1+γ4κϑ+ϰλϰℛ1Y−1Yn+1.
Sinceℛ1Z≤1 and ℛ1Y≤1, then Un is monotonically decreasing. We have Un≥0, and then limn⟶∞Un≥0 and hence limn⟶∞ΔUn=0, which implies limn⟶∞Fn=F∗, limn⟶∞Kn=K∗, limn⟶∞Sn=S∗,limn⟶∞Hn=H∗, limn⟶∞ℛ1Y−1Yn=0, and limn⟶∞ℛ1Z−1Zn=0. We have four cases:(i)
ℛ1Y=1 and ℛ1Z=1, and then from equation (14),(62)0=εϰF∗H∗+mK∗−aS∗−δS∗limn⟶∞Zn+1.Applying the steady state condition (60), we obtain limn⟶∞Zn=0. In addition, from equation (15), we get(63)0=ηS∗−cH∗−dH∗limn⟶∞Yn+1.From equation (60), we get limn⟶∞Yn=0.(ii)
ℛ1Y=1, ℛ1Z<1, and limn⟶∞Zn=0. From equation (63), we get limn⟶∞Yn=0.(iii)
ℛ1Y<1, ℛ1Z=1, and limn⟶∞Yn=0. From equation (62), we get limn⟶∞Zn=0.(iv)
ℛ1Y<1, ℛ1Z<1, limn⟶∞Yn=0, and limn⟶∞Zn=0. It follows that if ℛ1Z≤1 and ℛ1Y≤1, then limn⟶∞Fn=F∗, limn⟶∞Kn=K∗, limn⟶∞Sn=S∗,limn⟶∞Hn=H∗, limn⟶∞Yn=0, and limn⟶∞Zn=0. Then, Θ∗ is GAS.Theorem 3.
Letℛ1Y>1 and ℛ2Z≤1, then Θ¯ of models (12)–(17) is GAS.Proof.
DefineWnFn,Kn,Sn,Hn,Yn,Zn as(64)Wn=1ΠhF¯GFnF¯+γ1K¯GKnK¯+γ2S¯GSnS¯+γ31+Πhc+dY¯H¯GHnH¯+γ4Y¯GYnY¯+γ5Zn.
Clearly,WnFn,Kn,Sn,Hn,Yn,Zn>0 for all Fn,Kn,Sn,Hn,Yn,Zn>0 and WnF¯,K¯,S¯,H¯,Y¯,0=0. Computing ΔWn=Wn+1−Wn as(65)ΔWn=1ΠhF¯GFn+1F¯+γ1K¯GKn+1K¯+γ2S¯GSn+1S¯+γ31+Πhc+dY¯H¯GHn+1H¯+γ4Y¯GYn+1Y¯+γ5Zn+1−1ΠhF¯GFnF¯+γ1K¯GKnK¯+γ2S¯GSnS¯+γ31+Πhc+dY¯H¯GHnH¯+γ4Y¯GYnY¯+γ5Zn=1ΠhF¯Fn+1F¯−FnF¯+lnFnFn+1+γ1K¯Kn+1K¯−KnK¯+lnKnKn+1+γ2S¯Sn+1S¯−SnS¯+lnSnSn+1+γ3H¯Hn+1H¯−HnH¯+lnHnHn+1+γ4Y¯Yn+1Y¯−YnY¯+lnYnYn+1+γ5Zn+1−Zn+γ3c+dY¯H¯GHn+1H¯−GHnH¯.
Using inequality (49), we get(66)ΔWn≤1ΠhF¯Fn+1−FnF¯+FnFn+1−1+γ1K¯Kn+1−KnK¯+KnKn+1−1+γ2S¯Sn+1−SnS¯+SnSn+1−1+γ3H¯Hn+1−HnH¯+HnHn+1−1+γ4Y¯Yn+1−YnY¯+YnYn+1−1+γ5Zn+1−Zn+γ3c+dY¯H¯GHn+1H¯−GHnH¯=1Πh1−F¯Fn+1Fn+1−Fn+γ11−K¯Kn+1Kn+1−Kn+γ21−S¯Sn+1Sn+1−Sn+γ31−H¯Hn+1Hn+1−Hn+γ41−Y¯Yn+1Yn+1−Yn+γ5Zn+1−Zn+γ3c+dY¯H¯Hn+1H¯−HnH¯+lnHnHn+1.
From equations (12)–(17), we have(67)ΔWn≤1−F¯Fn+1ω−ϑFn+1−ϰFn+1Hn+γ11−K¯Kn+11−εϰFn+1Hn−α+mKn+1+γ21−S¯Sn+1εϰFn+1Hn+mKn+1−aSn+1−δSn+1Zn+1+γ31−H¯Hn+1ηSn+1−cHn+1−dHn+1Yn+1+γ41−Y¯Yn+1κHn+1Yn+1−λYn+1+γ5gSn+1Zn+1−μZn+1+γ3c+dY¯H¯Hn+1H¯−HnH¯+lnHnHn+1=1−F¯Fn+1ω−ϑFn+1+ϰF¯Hn−γ11−εϰK¯Fn+1HnKn+1+γ1α+mK¯−γ2εϰS¯Fn+1HnSn+1−γ2mS¯Kn+1Sn+1+γ2aS¯+γ2δS¯Zn+1−γ3ηH¯Sn+1Hn+1+γ3cH¯+γ3dH¯Yn+1−γ4λYn+1+γ4λY¯−γ5μZn+1−γ3c+dY¯Hn+γ3c+dY¯H¯lnHnHn+1.
We have(68)ω=ϑF¯+ϰF¯H¯,1−εϰF¯H¯=α+mK¯,εϰF¯H¯+mK¯=aS¯,ηS¯=cH¯+dH¯Y¯,λ=κH¯,and we get(69)γ11−εϰF¯H¯=γ1α+mK¯=γ2mK¯,ϰF¯H¯=γ2aS¯=γ3ηS¯,ϰF¯=γ3c+dY¯,ΔWn≤1−F¯Fn+1ϑF¯+ϰF¯H¯−ϑFn+1−γ11−εϰF¯H¯Fn+1HnK¯F¯H¯Kn+1+γ11−εϰF¯H¯−γ2εϰF¯H¯S¯Fn+1HnSn+1F¯H¯−γ11−εϰF¯H¯S¯Kn+1Sn+1K¯+ϰF¯H¯+γ5gS¯−μgZn+1−ϰF¯H¯H¯Sn+1Hn+1S¯+ϰF¯H¯+ϰF¯H¯lnHnHn+1=−ϑFn+1Fn+1−F¯2+γ11−εϰF¯H¯4−F¯Fn+1−Fn+1HnK¯F¯H¯Kn+1−S¯Kn+1Sn+1K¯−H¯Sn+1Hn+1S¯+lnHnHn+1+γ2εϰF¯H¯3−F¯Fn+1−S¯Fn+1HnSn+1F¯H¯−H¯Sn+1Hn+1S¯+lnHnHn+1+γ5gS¯−S^Zn+1=−ϑFn+1Fn+1−F¯2−γ11−εϰF¯H¯GF¯Fn+1+GFn+1HnK¯F¯H¯Kn+1+GS¯Kn+1Sn+1K¯+GH¯Sn+1Hn+1S¯−γ2εϰF¯H¯GF¯Fn+1+GS¯Fn+1HnSn+1F¯H¯+GH¯Sn+1Hn+1S¯+γ5μℛ2Z−1Zn+1.
Sinceℛ2Z≤1, then ΔWn≤0. Then, the sequence Wn is monotonically decreasing. Since Wn≥0, then limn⟶∞Wn≥0 and hence limn⟶∞ΔWn=0. Thus, limn⟶∞Fn=F¯, limn⟶∞Kn=K¯, limn⟶∞Sn=S¯, limn⟶∞Hn=H¯, and limn⟶∞ℛ2Z−1Zn=0. We have two cases:(i)
ℛ2Z=1, and then from equation (14),(70)0=εϰF¯H¯+mK¯−aS¯−δS¯limn⟶∞Zn+1,and we getlimn⟶∞Zn=0. From equation (15), we get(71)0=ηS¯−cH¯−dH¯limn⟶∞Yn+1.This giveslimn⟶∞Yn=Y¯.(ii)
ℛ2Z<1 and limn⟶∞Zn=0. From equation (71), we get limn⟶∞Yn=Y¯. Hence, Θ¯ is GAS.Theorem 4.
Letℛ1Z>1 and ℛ1Y/ℛ2Z≤1, then Θ^ of models (12)–(17) is GAS.Proof.
DefineℳnFn,Kn,Sn,Hn,Yn,Zn as(72)ℳn=1ΠhF^GFnF^+γ1K^GKnK^+γ2S^GSnS^+γ3a+δZ^a1+ΠhcH^GHnH^+γ4a+δZ^aYn+γ5Z^GZnZ^.
Clearly,ℳnFn,Kn,Sn,Hn,Yn,Zn>0 for all Fn,Kn,Sn,Hn,Yn,Zn>0 and ℳnF^,K^,S^,H^,0,Z^=0. We have Δℳn=ℳn+1−ℳn as(73)Δℳn=1ΠhF^GFn+1F^+γ1K^GKn+1K^+γ2S^GSn+1S^+γ3a+δZ^a1+ΠhcH^GHn+1H^+γ4a+δZ^aYn+1+γ5Z^GZn+1Z^−1ΠhF^GFnF^+γ1K^GKnK^+γ2S^GSnS^+γ3a+δZ^a1+ΠhcH^GHnH^+γ4a+δZ^aYn+γ5Z^GZnZ^=1ΠhF^Fn+1F^−FnF^+lnFnFn+1+γ1K^Kn+1K^−KnK^+lnKnKn+1+γ2S^Sn+1S^−SnS^+lnSnSn+1+γ3a+δZ^aH^Hn+1H^−HnH^+lnHnHn+1+γ4a+δZ^aYn+1−Yn+γ5Z^Zn+1Z^−ZnZ^+lnZnZn+1+γ3ca+δZ^aGHn+1H^−GHnH^.
Using inequality (49), we get(74)Δℳn≤1ΠhF^Fn+1−FnF^+FnFn+1−1+γ1K^Kn+1−KnK^+KnKn+1−1+γ2S^Sn+1−SnS^+SnSn+1−1+γ3a+δZ^aH^Hn+1−HnH^+HnHn+1−1+γ4a+δZ^aYn+1−Yn+γ5Z^Zn+1Z^−ZnZ^+ZnZn+1−1+γ3ca+δZ^aH^GHn+1H^−GHnH^=1Πh1−F^Fn+1Fn+1−Fn+γ11−K^Kn+1Kn+1−Kn+γ21−S^Sn+1Sn+1−Sn+γ3a+δZ^a1−H^Hn+1Hn+1−Hn+γ4a+δZ^aYn+1−Yn+γ51−Z^Zn+1Zn+1−Zn+γ3ca+δZ^aH^Hn+1H^−HnH^+lnHnHn+1.
From equations (12)–(17), we have(75)Δℳn≤1−F^Fn+1ω−ϑFn+1−ϰFn+1Hn+γ11−K^Kn+11−εϰFn+1Hn−α+mKn+1+γ21−S^Sn+1εϰFn+1Hn+mKn+1−aSn+1−δSn+1Zn+1+γ3a+δZ^a1−H^Hn+1ηSn+1−cHn+1−dHn+1Yn+1+γ4a+δZ^aκHn+1Yn+1−λYn+1+γ51−Z^Zn+1gSn+1Zn+1−μZn+1+γ3ca+δZ^aH^Hn+1H^−HnH^+lnHnHn+1=1−F^Fn+1ω−ϑFn+1+ϰF^Hn−γ11−εϰK^Fn+1HnKn+1+γ1α+mK^−γ2εϰS^Fn+1HnSn+1−γ2mS^Kn+1Sn+1−γ2aSn+1+γ2aS^+γ2δS^Zn+1+γ3a+δZ^aηSn+1−γ3a+δZ^ηaH^Sn+1Hn+1+γ3a+δZ^acH^+γ3a+δZ^adH^Yn+1−γ4a+δZ^aλYn+1−γ5gZ^Sn+1−γ5μZn+1+γ5μZ^−γ3ca+δZ^aHn+γ3ca+δZ^aH^lnHnHn+1.
We have(76)ω=ϑF^+ϰF^H^,1−εϰF^H^=α+mK^,εϰF^H^+mK^=a+δZ^S^,ηS^=cH^,S^=μg,and we get(77)γ11−εϰF^H^=γ1α+mK^=γ2mK^,ϰF^H^=γ2a+δZ^S^=γ3a+δZ^ηaS^,Δℳn≤1−F^Fn+1ϑF^+ϰF^H^−ϑFn+1−γ11−εϰF^H^Fn+1HnK^F^H^Kn+1+γ11−εϰF^H^−γ2εϰF^H^S^Fn+1HnSn+1F^H^−γ11−εϰF^H^S^Kn+1Sn+1K^+ϰF^H^−ϰF^H^H^Sn+1Hn+1S^+ϰF^H^+γ3da+δZ^aH^−λκYn+1+ϰF^H^lnHnHn+1=−ϑFn+1Fn+1−F^2+γ11−εϰF^H^4−F^Fn+1−Fn+1HnK^F^H^Kn+1−S^Kn+1Sn+1K^−H^Sn+1Hn+1S^+lnHnHn+1+γ2εϰF^H^3−F^Fn+1−S^Fn+1HnSn+1F^H^−H^Sn+1Hn+1S^+lnHnHn+1+γ3da+δZ^aH^−H˜Yn+1=−ϑFn+1Fn+1−F^2−γ11−εϰF^H^GF^Fn+1+GFn+1HnK^F^H^Kn+1+GS^Kn+1Sn+1K^+GH^Sn+1Hn+1S^−γ2εϰF^H^GF^Fn+1+GS^Fn+1HnSn+1F^H^+GH^Sn+1Hn+1S^+γ2a+δZ^dλκηℛ1Yℛ2Z−1Yn+1.
Sinceℛ1Z>1 and if ℛ1Y/ℛ2Z≤1, then Δℳn≤0. Hence, the sequence ℳn is monotonically decreasing. Since ℳn≥0, then limn⟶∞ℳn≥0 and hence limn⟶∞Δℳn=0. Thus, limn⟶∞Fn=F^, limn⟶∞Kn=K^, limn⟶∞Sn=S^, limn⟶∞Hn=H^, and limn⟶∞ℛ1Y/ℛ2Z−1Yn+1=0. We have two cases:(i)
ℛ1Y/ℛ2Z=1, and from equation (14),(78)0=εϰF^H^+mK^−aS^−δS^limn⟶∞Zn+1,and this giveslimn⟶∞Zn=Z^; moreover, from equation (15), we have(79)0=ηS^−cH^−dH^limn⟶∞Yn+1,then we getlimn⟶∞Yn=0.(ii)
ℛ1Y/ℛ2Z<1 and limn⟶∞Yn=0. From equation (78), we get limn⟶∞Zn=Z^ and hence Θ^ is GAS.Theorem 5.
Letℛ1Y>ℛ2Z>1, then Θ˜ of models (12)–(17) is GAS.Proof.
DefineVnFn,Kn,Sn,Hn,Yn,Zn as(80)Vn=1ΠhF˜GFnF˜+γ1K˜GKnK˜+γ2S˜GSnS˜+γ3a+δZ˜a1+Πhc+dY˜H˜GHnH˜+γ4a+δZ˜aY˜GYnY˜+γ5Z˜GZnZ˜.
Clearly,VnFn,Kn,Sn,Hn,Yn,Zn>0 for all Fn,Kn,Sn,Hn,Yn,Zn>0 and VnF˜,K˜,S˜,H˜,Y˜,Z˜=0. Computing ΔVn=Vn+1−Vn as(81)ΔVn=1ΠhF˜GFn+1F˜+γ1K˜GKn+1K˜+γ2S˜GSn+1S˜+γ3a+δZ˜a1+Πhc+dY˜H˜GHn+1H˜+γ4a+δZ˜aY˜GYn+1Y˜+γ5Z˜GZn+1Z˜−1ΠhF˜GFnF˜+γ1K˜GKnK˜+γ2S˜GSnS˜+γ3a+δZ˜a1+Πhc+dY˜H˜GHnH˜+γ4a+δZ˜aY˜GYnY˜+γ5Z˜GZnZ˜=1ΠhF˜Fn+1F˜−FnF˜+lnFnFn+1+γ1K˜Kn+1K˜−KnK˜+lnKnKn+1+γ2S˜Sn+1S˜−SnS˜+lnSnSn+1+γ3a+δZ˜aH˜Hn+1H˜−HnH˜+lnHnHn+1+γ4a+δZ˜aY˜Yn+1Y˜−YnY˜+lnYnYn+1+γ5Z˜Zn+1Z˜−ZnZ˜+lnZnZn+1+γ3c+dY˜a+δZ˜aH˜GHn+1H˜−GHnH˜.
Using inequality (49), we get(82)ΔVn≤1ΠhF˜Fn+1−FnF˜+FnFn+1−1+γ1K˜Kn+1−KnK˜+KnKn+1−1+γ2S˜Sn+1−SnS˜+SnSn+1−1+γ3a+δZ˜aH˜Hn+1−HnH˜+HnHn+1−1+γ4a+δZ˜aY˜Yn+1−YnY˜+YnYn+1−1+γ5Z˜Zn+1−ZnZ˜+ZnZn+1−1+γ3c+dY˜a+δZ˜aH˜GHn+1H˜−GHnH˜=1Πh1−F˜Fn+1Fn+1−Fn+γ11−K˜Kn+1Kn+1−Kn+γ21−S˜Sn+1Sn+1−Sn+γ3a+δZ˜a1−H˜Hn+1Hn+1−Hn+γ4a+δZ˜a1−Y˜Yn+1Yn+1−Yn+γ51−Z˜Zn+1Zn+1−Zn+γ3c+dY˜a+δZ˜aH˜Hn+1H˜−HnH˜+lnHnHn+1.
From equations (12)–(17), we have(83)ΔVn≤1−F˜Fn+1ω−ϑFn+1−ϰFn+1Hn+γ11−K˜Kn+11−εϰFn+1Hn−α+mKn+1+γ21−S˜Sn+1εϰFn+1Hn+mKn+1−aSn+1−δSn+1Zn+1+γ3a+δZ˜a1−H˜Hn+1ηSn+1−cHn+1−dHn+1Yn+1+γ4a+δZ˜a1−Y˜Yn+1κHn+1Yn+1−λYn+1+γ51−Z˜Zn+1gSn+1Zn+1−μZn+1+γ3c+dY˜a+δZ˜H˜aHn+1H˜−HnH˜+lnHnHn+1=1−F˜Fn+1ω−ϑFn+1+ϰF˜Hn−γ11−εϰK˜Fn+1HnKn+1+γ1α+mK˜−γ2εϰS˜Fn+1HnSn+1−mγ2S˜Kn+1Sn+1−aγ2Sn+1+aγ2S˜+δγ2S˜Zn+1+γ3a+δZ˜ηaSn+1−γ3a+δZ˜ηaH˜Sn+1Hn+1+γ3a+δZ˜acH˜+γ3da+δZ˜aH˜Yn+1−γ4a+δZ˜aλYn+1+γ4a+δZ˜aλY˜−γ5μZn+1−γ5gZ˜Sn+1+γ5μZ˜−γ3c+dY˜a+δZ˜aHn+γ3c+dY˜a+δZ˜aH˜lnHnHn+1.
We have(84)ω=ϑF˜+ϰF˜H˜,1−εϰF˜H˜=α+mK˜,εϰF˜H˜+mK˜=a+δZ˜S˜,ηS˜=c+dY˜H˜,S˜=μg,H˜=λκ,and we get(85)γ11−εϰF˜H˜=γ1α+mK˜=γ2mK˜,ϰF˜H˜=γ2a+δZ˜S˜=γ3a+δZ˜ηaS˜,ΔVn≤1−F˜Fn+1ϑF˜+ϰF˜H˜−ϑFn+1−γ11−εϰF˜H˜Fn+1HnK˜F˜H˜Kn+1+γ11−εϰF˜H˜−γ2εϰF˜H˜S˜Fn+1HnSn+1F˜H˜−γ11−εϰF˜H˜S˜Kn+1Sn+1K˜+ϰF˜H˜−ϰF˜H˜H˜Sn+1Hn+1S˜+ϰF˜H˜+ϰF˜H˜lnHnHn+1=−ϑFn+1Fn+1−F˜2+γ11−εϰF˜H˜4−F˜Fn+1−Fn+1HnK˜F˜H˜Kn+1−S˜Kn+1Sn+1K˜−H˜Sn+1Hn+1S˜+lnHnHn+1+γ2εϰF˜H˜3−F˜Fn+1−S˜Fn+1HnSn+1F˜H˜−H˜Sn+1Hn+1S˜+lnHnHn+1=−ϑFn+1Fn+1−F˜2−γ11−εϰF˜H˜GF˜Fn+1+GFn+1HnK˜F˜H˜Kn+1+GS˜Kn+1Sn+1K˜+GH˜Sn+1Hn+1S˜−γ2εϰF˜H˜GF˜Fn+1+GS˜Fn+1HnSn+1F˜H˜+GH˜Sn+1Hn+1S˜.
We note thatΔVn≤0. It follows that Vn is a monotone decreasing sequence. Since Vn≥0, then limn⟶∞Vn≥0 and hence limn⟶∞ΔVn=0. Thus, limn⟶∞Fn=F˜, limn⟶∞Kn=K˜, limn⟶∞Sn=S˜, and limn⟶∞Hn=H˜. From equations (14) and (15), we have(86)0=εϰF˜H˜+mK˜−aS˜−δS˜limn⟶∞Zn+1,0=ηS˜−cH˜−dH˜limn⟶∞Yn+1,then limn⟶∞Yn=Y˜ and limn⟶∞Zn=Z˜. Hence, Θ˜ is GAS.


## 3. Numerical Simulations

To conduct numerical simulations for models (12)–(17), we use the following data: ω=10, ϑ=0.01, ε=0.3, α=0.4, m=0.1, a=0.6, δ=0.1, η=1.2, c=4, λ=0.1, μ=0.03, d=0.8, and h=0.1. The remaining parameters will be chosen as below.


### 3.1. Stability of Steady States

We consider the following initial values:IV1:F0=800, K0=8, S0=1, H0=0.5, Y0=1, and Z0=0.5IV2:F0=600, K0=10, S0=3, H0=1, Y0=3, and Z0=1IV3:F0=400, K0=12, S0=5, H0=2, Y0=5, and Z0=2We choose variedϰ, κ, and g:Case (1)ϰ=0.003, κ=0.05, and g=0.002. This yields ℛ0=0.66<1. Figure 1 shows that the concentration of healthy cells increases and tends to the value F0=1000, while the concentrations of the other compartments decrease and tend to zero for all initial IV1–IV3. Therefore, Θ0 is GAS. This illustrates the result of Theorem 1.Case (2)ϰ=0.01, κ=0.05, and g=0.002. These values give ℛ1Y=0.7333<1<ℛ0=2.2 and ℛ1Z=0.4<1<ℛ0=2.2. From Figure 2, it can be seen that the solutions of the system tend to the steady state Θ∗=454.545,7.636,4,1.2,0,0, for all the three initial values IV1–IV3. As a result, Θ∗ exists, and based on Theorem 2, Θ∗ is GAS.Case (3)ϰ=0.01, κ=0.2, and g=0.002, and then ℛ1Y=1.46667>1 and ℛ1Z=0.4<1. Figure 3 shows that the trajectories of the system reach the steady state Θ¯=666.667,4.667,2.444,0.5,2.333,0 for all the initial values IV1–IV3. This shows the existence of Θ¯. According to Theorem 3, Θ¯ is GAS.Case (4)ϰ=0.01, κ=0.005, and g=0.01, and then ℛ1Z=1.15789>1 and ℛ1Y/ℛ2Z=0.045<1. Figure 4 shows that the solutions starting with initial values IV1–IV3 converge to the steady state Θ^=526.316,6.632,3,0.9,0,0.947. This shows that Θ^ exists, and it is GAS according to Theorem 4.Case (5)ϰ=0.01, κ=0.3, and g=0.02, and then ℛ0=2.2>1, ℛ1Y=1.65>1 and ℛ2Z=1.222>1. Figure 5 displays that the solutions with initial values IV1–IV3 approach the steady state Θ˜=750,3.5,1.5,0.333,1.75,1.333. Hence, Θ˜ exists; moreover, on the basis of Theorem 5, Θ˜ is GAS.Figure 1
The solutions of systems (12)–(17) for case ℛ0≤1.
(a)(b)(c)(d)(e)(f)Figure 2
The solutions of systems (12)–(17) for case ℛ1Y≤1<ℛ0 and ℛ1Z≤1<ℛ0.
(a)(b)(c)(d)(e)(f)Figure 3
The solutions of systems (12)–(17) for case ℛ1Y>1 and ℛ1Z≤1.
(a)(b)(c)(d)(e)(f)Figure 4
The solutions of systems (12)–(17) for case ℛ1Z>1 and ℛ1Y/ℛ2Z≤1.
(a)(b)(c)(d)(e)(f)Figure 5
The solutions of systems (12)–(17) for case ℛ1Y>ℛ2Z>1.
(a)(b)(c)(d)(e)(f)


### 3.2. Effect of Antipathogenic Drugs

In this section, we study the effect of antipathogenic drugs on the pathogen dynamics. We develop models (12)–(17) by incorporating two categories of antipathogenic drugs with efficacies ρ1∈0,1 and ρ2∈0,1 as(87)Fn+1−FnΠh=ω−ϑFn+1−1−ρ1ϰFn+1Hn,Kn+1−KnΠh=1−ρ11−εϰFn+1Hn−α+mKn+1,Sn+1−SnΠh=1−ρ1εϰFn+1Hn+mKn+1−aSn+1−δSn+1Zn+1,Hn+1−HnΠh=1−ρ2ηSn+1−cHn+1−dHn+1Yn+1,Yn+1−YnΠh=κHn+1Yn+1−λYn+1,Zn+1−ZnΠh=gSn+1Zn+1−μZn+1.The basic reproduction number for this system is given by(88)ℛ0ρ1,ρ2=1−ρ11−ρ2ϰηωαε+mϑacα+m.Without the loss of generality, we letρ=ρ1=ρ2, and then(89)ℛ0ρ=1−ρ2ϰηωαε+mϑacα+m.It follows that(90)ℛ0ρ=1−ρ2ℛ00.Since the target of antipathogenic drugs is to clear the pathogens from the body, we have to determine the minimum drug efficacyρmin such that ℛ0ρ≤1 for all ρmin≤ρ≤1. We can find the value of ρmin by solving the following algebraic equation:(91)1−ρ2ℛ00=1.The positive root of equation (91) is given by(92)ρ=1−1ℛ00.Therefore,(93)ρmin=max0,1−1ℛ00.We takeϰ=0.01, κ=0.3, and g=0.02. We get ρmin=0.3258. We consider the initial IV3 and choose different values of ρ. Figure 6 shows the effect of antipathogenic drugs ρ on the pathogen dynamics. It is observed that as ρ is increased, the concentration of the healthy cells is increased, while the concentrations of the other compartments are decreased. In addition, we have(i)
If0.3258≤ρ≤1, then ℛ0ρ≤1 and Θ0 is GAS. In this case, the antipathogenic drugs can remove the pathogens from the body.(ii)
Ifρ<0.3258, then ℛ0ρ>1 and Θ0 had lost its stability and the solution trajectories of the system converge to one of the other steady states of the system. This means that using the drug treatment with less efficacy will not able to clear the pathogens from the body.Figure 6
The solutions of systems (12)–(17) with different values of ρ.
(a)(b)(c)(d)(e)(f)


### 3.3. Comparison Results

In this section, we perform a comparison between the NSFD method and Runge–Kutta method with a variable step size presented in MATLAB (ODE45). We compare the solutions of NSFD and ODE45 for Case (5). The numerical results show that there is a very good agreement between the result of NSFD (h=0.1 and h=0.5) and that of ODE45 (see Figure 7). Of course, if one considers the Runge–Kutta method with a fixed numerical step size h, then the numerical error will increase by increasing h. Moreover, the global stability of the discrete-time model obtained by the Runge–Kutta method may not guaranteed. In the case of NSFD, we can see that the increasing of h does not affect the stability of the steady states. This is in agreement with the results of several papers [37–40] that NSFD preserves the essential qualitative features of these models independently of the chosen step size h.Figure 7
Comparison between the solutions of NSFD and ODE45.
(a)(b)(c)(d)(e)(f)


## 3.1. Stability of Steady States

We consider the following initial values:IV1:F0=800, K0=8, S0=1, H0=0.5, Y0=1, and Z0=0.5IV2:F0=600, K0=10, S0=3, H0=1, Y0=3, and Z0=1IV3:F0=400, K0=12, S0=5, H0=2, Y0=5, and Z0=2We choose variedϰ, κ, and g:Case (1)ϰ=0.003, κ=0.05, and g=0.002. This yields ℛ0=0.66<1. Figure 1 shows that the concentration of healthy cells increases and tends to the value F0=1000, while the concentrations of the other compartments decrease and tend to zero for all initial IV1–IV3. Therefore, Θ0 is GAS. This illustrates the result of Theorem 1.Case (2)ϰ=0.01, κ=0.05, and g=0.002. These values give ℛ1Y=0.7333<1<ℛ0=2.2 and ℛ1Z=0.4<1<ℛ0=2.2. From Figure 2, it can be seen that the solutions of the system tend to the steady state Θ∗=454.545,7.636,4,1.2,0,0, for all the three initial values IV1–IV3. As a result, Θ∗ exists, and based on Theorem 2, Θ∗ is GAS.Case (3)ϰ=0.01, κ=0.2, and g=0.002, and then ℛ1Y=1.46667>1 and ℛ1Z=0.4<1. Figure 3 shows that the trajectories of the system reach the steady state Θ¯=666.667,4.667,2.444,0.5,2.333,0 for all the initial values IV1–IV3. This shows the existence of Θ¯. According to Theorem 3, Θ¯ is GAS.Case (4)ϰ=0.01, κ=0.005, and g=0.01, and then ℛ1Z=1.15789>1 and ℛ1Y/ℛ2Z=0.045<1. Figure 4 shows that the solutions starting with initial values IV1–IV3 converge to the steady state Θ^=526.316,6.632,3,0.9,0,0.947. This shows that Θ^ exists, and it is GAS according to Theorem 4.Case (5)ϰ=0.01, κ=0.3, and g=0.02, and then ℛ0=2.2>1, ℛ1Y=1.65>1 and ℛ2Z=1.222>1. Figure 5 displays that the solutions with initial values IV1–IV3 approach the steady state Θ˜=750,3.5,1.5,0.333,1.75,1.333. Hence, Θ˜ exists; moreover, on the basis of Theorem 5, Θ˜ is GAS.Figure 1
The solutions of systems (12)–(17) for case ℛ0≤1.
(a)(b)(c)(d)(e)(f)Figure 2
The solutions of systems (12)–(17) for case ℛ1Y≤1<ℛ0 and ℛ1Z≤1<ℛ0.
(a)(b)(c)(d)(e)(f)Figure 3
The solutions of systems (12)–(17) for case ℛ1Y>1 and ℛ1Z≤1.
(a)(b)(c)(d)(e)(f)Figure 4
The solutions of systems (12)–(17) for case ℛ1Z>1 and ℛ1Y/ℛ2Z≤1.
(a)(b)(c)(d)(e)(f)Figure 5
The solutions of systems (12)–(17) for case ℛ1Y>ℛ2Z>1.
(a)(b)(c)(d)(e)(f)


## 3.2. Effect of Antipathogenic Drugs

In this section, we study the effect of antipathogenic drugs on the pathogen dynamics. We develop models (12)–(17) by incorporating two categories of antipathogenic drugs with efficacies ρ1∈0,1 and ρ2∈0,1 as(87)Fn+1−FnΠh=ω−ϑFn+1−1−ρ1ϰFn+1Hn,Kn+1−KnΠh=1−ρ11−εϰFn+1Hn−α+mKn+1,Sn+1−SnΠh=1−ρ1εϰFn+1Hn+mKn+1−aSn+1−δSn+1Zn+1,Hn+1−HnΠh=1−ρ2ηSn+1−cHn+1−dHn+1Yn+1,Yn+1−YnΠh=κHn+1Yn+1−λYn+1,Zn+1−ZnΠh=gSn+1Zn+1−μZn+1.The basic reproduction number for this system is given by(88)ℛ0ρ1,ρ2=1−ρ11−ρ2ϰηωαε+mϑacα+m.Without the loss of generality, we letρ=ρ1=ρ2, and then(89)ℛ0ρ=1−ρ2ϰηωαε+mϑacα+m.It follows that(90)ℛ0ρ=1−ρ2ℛ00.Since the target of antipathogenic drugs is to clear the pathogens from the body, we have to determine the minimum drug efficacyρmin such that ℛ0ρ≤1 for all ρmin≤ρ≤1. We can find the value of ρmin by solving the following algebraic equation:(91)1−ρ2ℛ00=1.The positive root of equation (91) is given by(92)ρ=1−1ℛ00.Therefore,(93)ρmin=max0,1−1ℛ00.We takeϰ=0.01, κ=0.3, and g=0.02. We get ρmin=0.3258. We consider the initial IV3 and choose different values of ρ. Figure 6 shows the effect of antipathogenic drugs ρ on the pathogen dynamics. It is observed that as ρ is increased, the concentration of the healthy cells is increased, while the concentrations of the other compartments are decreased. In addition, we have(i)
If0.3258≤ρ≤1, then ℛ0ρ≤1 and Θ0 is GAS. In this case, the antipathogenic drugs can remove the pathogens from the body.(ii)
Ifρ<0.3258, then ℛ0ρ>1 and Θ0 had lost its stability and the solution trajectories of the system converge to one of the other steady states of the system. This means that using the drug treatment with less efficacy will not able to clear the pathogens from the body.Figure 6
The solutions of systems (12)–(17) with different values of ρ.
(a)(b)(c)(d)(e)(f)


## 3.3. Comparison Results

In this section, we perform a comparison between the NSFD method and Runge–Kutta method with a variable step size presented in MATLAB (ODE45). We compare the solutions of NSFD and ODE45 for Case (5). The numerical results show that there is a very good agreement between the result of NSFD (h=0.1 and h=0.5) and that of ODE45 (see Figure 7). Of course, if one considers the Runge–Kutta method with a fixed numerical step size h, then the numerical error will increase by increasing h. Moreover, the global stability of the discrete-time model obtained by the Runge–Kutta method may not guaranteed. In the case of NSFD, we can see that the increasing of h does not affect the stability of the steady states. This is in agreement with the results of several papers [37–40] that NSFD preserves the essential qualitative features of these models independently of the chosen step size h.Figure 7
Comparison between the solutions of NSFD and ODE45.
(a)(b)(c)(d)(e)(f)


## 4. Conclusion

We studied a discrete-time pathogen dynamic model with adaptive immunity. The model incorporated two classes of infected cells, latently infected cells and actively infected cells. We proved that the solutions of the proposed model are positive and bounded. We showed that the model has five steady states, and their existence are governed by four threshold parameters. We constructed Lyapunov functions to prove the global stability of the steady states. In addition, we conducted numerical simulations to illustrate the global stability of the steady states. We studied the effect of antipathogenic treatment on the pathogen dynamics. We observed that asρ is increased, the concentration of the healthy cells is increased while the concentration of pathogens is decreased. Since the pathogen-free steady state Θ0 is the desired steady state to be stabilized, we determined the critical drug efficacy ρmin by solving the equation ℛ0ρmin=1 and showed that Θ0 is GAS when ρ≥ρmin. We presented a comparison between the NSFD method and the Runge–Kutta method. The numerical results showed that there is a very good agreement between the result of NSFD and that of the Runge–Kutta method.We mention that the proposed pathogen dynamic model is constructed at a supermacroscopic scale, namely, by means of population dynamics with deterministic interactions. This pathogen model can be extended to incorporate stochastic interactions (see [49]). Recently, many authors argued that the pathogens move freely in the body and follow the Fickian diffusion (see, e.g., [50–53]). Therefore, it is more reasonable to study reaction-diffusion versions of our discrete-time model. In certain circumstances, some pathogens such as HIV can suppress immune response or even destroy it especially when the concentration of the pathogens is too high. In the literature, several pathogen dynamic models were proposed by considering either CTL immune impairment (see, e.g., [54–56]) or antibody immune impairment (see, e.g., [57, 58]). However, modeling the impairment of both CTL and antibody immune responses has not been studied before. We leave these extensions as a possible future work.



---

*Source: 1264175-2020-08-25.xml*