We study the global stability of within-host Chikungunya virus (CHIKV) infection models with antibodies. We incorporate two modes of infections, attaching a CHIKV to a host monocyte, and contacting an infected monocyte with an uninfected monocyte. The CHIKV-monocyte and infected-monocyte incidence rates are given by saturation. In the second model we consider two classes of infected monocytes, latently infected monocytes and actively infected monocytes. The global stability analysis of the equilibria are established using Lyapunov method. We support our theoretical results by numerical simulations.

Topics
Stability theory, Mathematical modeling, Insects, Diseases and conditions, Antibody, Antigen, Blood cells, Viruses

Chikungunya is one of the mosquito-borne diseases caused by Chikungunya virus (CHIKV). This kind of viruses is transmitted to humans by infected Aedes albopictus and Aedes agypti mosquito. CHIKV causes severe joint and muscle pain, fever, rash, headache, nausea and fatigue. In the literature, many researchers have constructed and analyzed mathematical models for the transmission of the CHIKV to human population (see e.g. Refs. 1–8). In other words, within-host CHIKV dynamics model has been presented in a recent paper9 as:

(1)
(2)
(3)
(4)

where, s, y, p and x are the concentrations of uninfected monocytes, infected monocytes, CHIKV particles and antibodies, respectively. New uninfected monocytes are generated with rate β and die with rate δs. The CHIKV-monocyte incidence rate is given by ηsy, where η is constant. Constants ϵ, c and m represent, respectively, the death rate of the infected monocytes, CHIKV and antibodies. Constant π is the generation rate of the CHIKV from actively infected monocytes. Antibodies attack the CHIKV at rate rxp. Once antigen is encountered, the antibodies expand at a constant rate λ and proliferate at rate ρxp.

Model (1)–(4) has been extended in Refs. 10 and 11 by considering general CHIKV-monocyte incidence rate. In Refs. 9–11 it has been assumed that the uninfected monocyte becomes infected by contacting with CHIKV(CHIKV-to-monocyte transmission). Long and Heise12 reported that the CHIKV can also spread by infected-to-monocyte transmission. Mathematical models of different viruses with both cellular and viral infections have been studied in several works.13–24 However, the dynamics of CHIKV with two routes of infection did not studied before.

The aim of the present paper is to propose and analyze two CHIKV dynamics models with two routes of infection. We consider saturated incidence rate which modifies the bilinear incidence presented in model (1)–(4). In the second model we incorporate two classes of infected monocytes, actively infected monocytes and latently infected monocytes. We calculate the basic reproduction number R0 which determines the existence and stability of the equilibria. To investigate the global stability of the equilibria we construct Lyapunov functions using the method presented25 and followed by Refs. 26–41.

Our proposed models can be extended by incorporating different types of time delay. Moreover, following the work of Gibelli et al.,42 CHIKV models with a stochastic parameters dynamics can also be studied.

We investigate the CHIKV infection model with saturated CHIKV-monocyte and infected-monocyte incidence

(5)
(6)
(7)
(8)

where the terms η1sp1+γ1p and η2sy1+γ2y represent the CHIKV-monocyte and infected-monocyte incidence rates, respectively, and γ1 and γ2 are the saturation constants.

Let M1, M2, M3 > 0 and define

Γ1={(s,y,p,x)∈R≥04:0≤s,y≤M1,0≤p≤M2,0≤x≤M3}.

Lemma 1.

For system (5)–(8), the compact set Γ1 is positively invariant.

Proof.
We have
ṡ∣s=0=β>0,ẏ∣y=0=η1sp1+γ1p≥0,for alls,p≥0,ṗ∣p=0=πy≥0,for ally≥0,ẋ∣x=0=λ>0.
Thus R≥04 is positively invariant with respect to system (5)–(8). Let us define
F1(t)=s(t)+y(t),F2(t)=p(t)+rρx(t).
Then from Eqs. (5)–(8) we get
Ḟ1(t)=β−δs(t)−ϵy(t)≤β−σ1(s(t)+y(t))=β−σ1F1(t),
where, σ1 = min{δ, ϵ}. Hence F1(t) ≤ M1, if F1(0) ≤ M1, where M1=βσ1. It follows that 0 ≤ s(t), y(t) ≤ M1 when 0 ≤ s(0) + y(0) ≤ M1. Moreover, we have
Ḟ2(t)=πy(t)−cp(t)+rρλ−mrρx(t)≤πM1+rρλ−σ2p(t)+rρx(t)=πM1+rρλ−σ2F2(t),
where, σ2 = min{c, m}. Hence F2(t) ≤ M2, if F2(0) ≤ M2, where M2=πM1+rρλσ2. Then, 0 ≤ p(t) ≤ M2 and 0 ≤ x(t) ≤ M3 if 0≤p(0)+rρx(0)≤M2, where M3=ρM2r.□

We define a threshold parameter

R0=(η1πm+η2cm+η2rλ)βϵδ(cm+rλ).

Lemma 2.

For system (5)–(8) if R0≤1, then there exists only one equilibrium E0 ∈ Γ1, and if R0>1, then there exist two equilibria E0 ∈ Γ1 and E1∈Γ1°, where Γ1° is the interior of Γ1.

Proof.
Let E(s, y, p, x) be any equilibrium satisfying
(9)
(10)
(11)
(12)
By solving Eqs. (9)–(12) we get two equilibria a CHIKV-free equilibrium E0 = (s0, 0, 0, x0), where s0=βδ and x0=λm and an infected equilibrium satisfying
C1p4+C2p3+C3p2+C4p+C5=0,
where
C1=− c2ϵ(γ1γ2δ+γ1 η1+γ2η2)ρ2,C2=− c2m2γ1γ2δϵ−c2m2γ2ϵη1−c2m2γ1ϵη2−2cmrγ1γ2δϵλ−2cmrγ2ϵη1λ−2cmrγ1ϵη2λ−r2γ1γ2δϵλ2 − r2γ2ϵη1λ2−r2γ1ϵη2λ2+2cmπγ1δϵρ+2c2mγ2δϵρ−2cmπβγ2η1ρ+2cmπϵη1ρ−2cmπβγ1η2ρ + 2c2mϵη2ρ+πrγ1δϵλρ+2crγ2δϵλρ−πrβγ2η1λρ+πrϵη1λρ−πrβγ1η2λρ+2crϵη2λρ−cπδϵρ2 + π2βη1ρ2+cπβη2ρ2,
C3=− c2m2γ1γ2δϵ−c2m2γ2ϵη1−c2m2γ1ϵη2−2cmrγ1γ2δϵλ−2cmrγ2ϵη1λ−2cmrγ1ϵη2λ−r2γ1γ2δϵλ2 − r2γ2ϵη1λ2−r2γ1ϵη2λ2+2cmπγ1δϵρ+2c2mγ2δϵρ−2cmπβγ2η1ρ+2cmπϵη1ρ−2cmπβγ1η2ρ + 2c2mϵη2ρ+πrγ1δϵλρ+2crγ2δϵλρ−πrβγ2η1λρ+πrϵη1λρ−πrβγ1η2λρ+2crϵη2λρ − cπδϵρ2+π2βη1ρ2+cπβη2ρ2,C4=− cm2πγ1δϵ−c2m2γ2δϵ+cm2πβγ2η1−cm2πϵη1+cm2πβγ1η2−c2m2ϵη2−mπrγ1δϵλ − 2cmrγ2δϵλ+mπrβγ2η1λ−mπrϵη1λ+mπrβγ1η2λ−2cmrϵη2λ−r2γ2δϵλ2−r2ϵη2λ2 + 2cmπδϵρ−2mπ2βη1ρ−2cmπβη2ρ+πrδϵλρ−πrβη2λρ,C5=− cm2πδϵ+m2π2βη1+cm2πβη2−mπrδϵλ+mπrβη2λ.
Let define a function X(p) as:
X(p)=C1p4+C2p3+C3p2+C4p+C5=0.
Then we obtain
X(0)=C5,Xmρ=− mr2ϵλ2mγ2η1+mγ1γ2δ+η2+γ2δ+η2ρρ2<0.
The constant C5 can be written as:
C5=mπδϵ(cm+rλ)R0−1
Hence if R0>1 then C5 > 0 and there exists p1∈(0,mρ) such that X(p1) = 0. Therefore, if R0>1⁠, then
x1=λm−ρp1>0,y1=(c+rx1)p1π>0,s1=y1(1+p1γ1)(1+y1γ2)ϵ(1+γ2y1)η1p1+η2y1(1+γ1p1)>0.
Then an infected equilibrium E1 = (s1, y1, p1, x1) exists when R0>1⁠.
Now we show that E0 ∈ Γ1 and E1∈Γ1°. Clearly, E0 ∈ Γ1. From the equilibrium conditions of E1 we have
β=δs1+η1s1p11+γ1p1+η2s1y11+γ2y1⇒δs1+ϵy1=β⇒0<s1<βδ≤M1,0<y1<βϵ≤M1.
Moreover, from Eqs. (11) and (12) we have
πy1−cp1−rx1p1+rρλ+ρx1p1−mx1=0
⇒cp1+mrρx1=πy1+rρλ<πM1+rρλ⇒0<p1<πM1+rρλc≤M2,0<x1<ρrπM1+rρλm≤ρM2r=M3.
It follows that, E1∈Γ°.□

Define a function G(z) = z − 1 − ln z.

Theorem 1.

The equilibrium E0 of system (5)–(8) is globally asymptotically stable in Γ1 when R0≤1.

Proof.
Letting R0≤1 and constructing a Lyapunov function U0(s, y, p, x) as:
U0(s,y,p,x)=s0Gss0+y+η1s0c+rx0p+rη1s0ρ(c+rx0)x0Gxx0.
Calculating dU0dt along system (5)–(8) we obtain
dU0dt=1−s0sβ−δs−η1sp1+γ1p−η2sy1+γ2y+η1sp1+γ1p+η2sy1+γ2y−ϵy+η1s0c+rx0πy−cp−rxp+rη1s0ρ(c+rx0)1−x0xλ+ρxp−mx=1−s0sβ−δs+η1s0p1+γ1p+η2s0y1+γ2y−ϵy+η1s0πc+rx0y−η1s0c+rx0cp−rη1s0x0c+rx0p+rη1s0ρ(c+rx0)1−x0xλ−mx=1−s0sβ−δs−η1s0γ1p21+γ1p−η2s0γ2y21+γ2y+η2s0y−ϵy+η1s0πc+rx0y+rη1s0ρ(c+rx0)1−x0xλ−mx.
Then
(13)
If R0≤1, then dU0dt≤0 for all s, y, p, x > 0, and dU0dt=0 if s = s0, x = x0, y = p = 0. Applying LaSalle’s Invariance Principle (LIP), we get E0 is globally asymptotically stable. □

Theorem 2.

The equilibrium E1 of system (5)–(8) is globally asymptotically stable in Γ1° when R0>1⁠.

Proof.
Define U1(s, y, p, x) as:
U1=s1Gss1+y1Gyy1+η1s1p1(1+γ1p1)πy1p1Gpp1+ rη1s1p1ρ(1+γ1p1)πy1x1Gxx1.
Calculating dU1dt along the trajectories of (5)–(8) we obtain
dU1dt=1−s1sβ−δs−η1sp1+γ1p−η2sy1+γ2y+1−y1yη1sp1+γ1p+η2sy1+γ2y−ϵy + η1s1p1(1+γ1p1)πy11−p1pπy−cp−rxp+rη1s1p1ρ(1+γ1p1)πy11−x1xλ+ρxp−mx=1−s1sβ−δs+η1s1p1+γ1p+η2s1y1+γ2y−η1sp1+γ1py1y−η2sy1+γ2yy1y−ϵy+ϵy1+η1s1p11+γ1p1yy1−η1s1p11+γ1p1p1ypy1 − η1s1p1(1+γ1p1)πy1cp+η1s1p1(1+γ1p1)πy1cp1+η1s1p1(1+γ1p1)πy1rxp1−rη1s1p1(1+γ1p1)πy1x1p+rη1s1p1ρ(1+γ1p1)πy11−x1xλ−mx.
Applying the equilibrium conditions for E1
β=δs1+η1s1p11+γ1p1+η2s1y11+γ2y1, ϵy1=η1s1p11+γ1p1+η2s1y11+γ2y1,cp1=πy1−rx1p1, λ=mx1−ρx1p1.
we get
dU1dt=− δ(s−s1)2s−η1s1p11+γ1p1γ1(p−p1)2(1+γ1p)(1+γ1p1)p1 − η2s1y11+γ2y1γ2(y−y1)2(1+γ2y)(1+γ2y1)y1 + η2s1y11+γ2y13−s1s−s(1+γ2y1)s1(1+γ2y)−1+γ2y1+γ2y1 + η1s1p11+γ1p14−s1s−spy1(1+γ1p1)s1p1y(1+γ1p)−p1ypy1−1+γ1p1+γ1p1 − η1s1p1(1+γ1p1)πy1rλρx1(x−x1)2x.
Using the rule
(14)
we get
14s1s+spy1(1+γ1p1)s1p1y(1+γ1p)+p1ypy1+1+γ1p1+γ1p1≥1,
13s1s+s(1+γ2y1)s1(1+γ2y)+1+γ2y1+γ2y1≥1.
Therefore, dU1dt≤0 for all s, y, p, x > 0 and dU1dt=0 when s = s1, y = y1, p = p1 and x = x1. The global stability of E1 is induced from LIP. □

In this section, we consider two classes of infected monocytes, actively infected monocytes (y) and latently infected monocytes (w). We propose the following CHIKV model:

(15)
(16)
(17)
(18)
(19)

where 0<n<1⁠, bw and dw are the activation and dearth rates of the latently infected monocytes, respectively.

Let M1L,M2L,M3L>0 and define the set

Γ2=(s,w,y,p,x)∈R≥05:0≤s,w,y≤M1L,   0≤p≤M2L,0≤x≤M3L .

Lemma 3.

The compact set Γ2 is positively invariant for system (15)–(19)

Proof.
We have
ṡ∣s=0=β>0,ẇ∣w=0=(1−n)η1sp1+γ1p+η2sy1+γ2y≥0,for alls,p≥0,ẏ∣y=0=nη1sp1+γ1p+dw≥0,for alls,p,w≥0,ṗ∣p=0=πy≥0,for ally≥0,ẋ∣x=0=λ>0.
Then, R≥05 is positively invariant for system (15)–(19). We let
H1(t)=s(t)+w(t)+y(t),H2(t)=p(t)+rρx(t),
then
Ḣ1(t)=β−δs(t)−bw(t)−ϵy(t)≤β−σ1L(s(t)+w(t)+y(t))=β−σ1LH1(t),
where, σ1L=min{δ,b,ϵ}. Hence H1(t)≤M1L, if H1(0)≤M1L, where M1L=βσ1L. Hence, 0≤s(t),w(t),y(t)≤M1L if 0≤s(0)+w(0)+y(0)≤M1L. Moreover, we have
Ḣ2(t)=πy(t)−cp(t)+rρλ−mrρx(t)≤πM1L+rρλ−σ2p(t)+rρx(t)=πM1L+rρλ−σ2LH2(t),
where, σ2L=σ2⁠. Hence H2(t)≤M2L, if H2(0)≤M2L, where M2L=πM1L+rρλσ2L. Thus, 0≤p(t)≤M2L and x(t)≤M3L if 0≤p(0)+rρx(0)≤M2L, where M3L=ρM2Lr.□

We define a threshold parameter for system (15)–(19) as:

R0L=β(d+bn)(η1πm+η2cm+η2rλ)ϵδ(cm+rλ)(b+d).

Lemma 4.

For the system (15)–(19) if R0L≤1, then there exists only one equilibrium E0 ∈ Γ2 and if R0L>1, then there exist two equilibria E0 ∈ Γ2 and E1∈Γ2°.

Proof.
The equilibria of system (15)–(19) satisfying
(20)
(21)
(22)
(23)
(24)
Solving Eqs. (20)–(24) we get a CHIKV-free equilibrium E0 = (s0, 0, 0, 0, x0), where s0=βδ and x0=λm. From Eqs. (20)–(24) we have
(25)
(26)
(27)
(28)
Substituting from Eqs. (25)–(28) into Eq. (22) we get
D1p4+D2p3+D3p2+D4p+D5=0,
where
D1=− c2(b+d)ϵ(γ1γ2δ+γ1η1+γ2η1)ρ2,D2=2bc2mγ1γ2δϵρ+2c2dmγ1γ2δϵρ+2bc2mγ1ϵη1ρ+2c2dmγ1ϵη1ρ+2bc2mγ2ϵη1ρ+2c2dmγ2ϵη1ρ + 2bcrγ1γ2δϵλρ+2cdrγ1γ2δϵλρ+2bcrγ1ϵη1λρ+2cdrγ1ϵη1λρ+2bcrγ2ϵη1λρ+2cdrγ2ϵη1λρ − bcπγ1δϵρ2−cdπγ1δϵρ2−bc2γ2δϵρ2−c2dγ2δϵρ2+cdπβγ1η1ρ2+bcnπβγ1η1ρ2 + cdπβγ2η1ρ2+bcnπβγ2η1ρ2−bc2ϵη1ρ2−c2dϵη1ρ2−bcπϵη1ρ2−cdπϵη1ρ2,
D3=− bc2m2γ1γ2δϵ−c2dm2γ1γ2δϵ−bc2m2γ1ϵη1−c2dm2γ1ϵη1−bc2m2γ2ϵη1−c2dm2γ2ϵη1 − 2bcmrγ1γ2δϵλ−2cdmrγ1γ2δϵλ−2bcmrγ1ϵη1λ−2cdmrγ1ϵη1λ−2bcmrγ2ϵη1λ−2cdmrγ2ϵη1λ − br2γ1γ2δϵλ2−dr2γ1γ2δϵλ2−br2γ1ϵη1λ2−dr2γ1ϵη1λ2−br2γ2ϵη1λ2−dr2γ2ϵη1λ2 + 2bcmπγ1δϵρ+2cdmπγ1δϵρ+2bc2mγ2δϵρ+2c2dmγ2δϵρ−2cdmπβγ1η1ρ−2bcmnπβγ1η1ρ − 2cdmπβγ2η1ρ−2bcmnπβγ2η1ρ+2bc2mϵη1ρ+2c2dmϵη1ρ+2bcmπϵη1ρ+2cdmπϵη1ρ + bπrγ1δϵλρ+dπrγ1δϵλρ+2bcrγ2δϵλρ+2cdrγ2δϵλρ−dπrβγ1η1λρ−bnπrβγ1η1λρ − dπrβγ2η1λρ−bnπrβγ2η1λρ+2bcrϵη1λρ+2cdrϵη1λρ+bπrϵη1λρ+dπrϵη1λρ − bcπδϵρ2−cdπδϵρ2+cdπβη1ρ2+bcnπβη1ρ2+dπ2βη1ρ2+bnπ2βη1ρ2,D4=− bcm2πγ1δϵ−cdm2πγ1δϵ−bc2m2γ2δϵ−c2dm2γ2δϵ+cdm2πβγ1η1+bcm2nπβγ1η1 + cdm2πβγ2η1+bcm2nπβγ2η1−bc2m2ϵη1−c2dm2ϵη1−bcm2πϵη1−cdm2πϵη1 − bmπrγ1δϵλ−dmπrγ1δϵλ−2bcmrγ2δϵλ−2cdmrγ2δϵλ+dmπrβγ1η1λ+bmnπrβγ1η1λ + dmπrβγ2η1λ+bmnπrβγ2η1λ−2bcmrϵη1λ−2cdmrϵη1λ−bmπrϵη1λ−dmπrϵη1λ − br2γ2δϵλ2−dr2γ2δϵλ2−br2ϵη1λ2−dr2ϵη1λ2+2bcmπδϵρ+2cdmπδϵρ−2cdmπβη1ρ − 2bcmnπβη1ρ−2dmπ2βη1ρ−2bmnπ2βη1ρ+bπrδϵλρ+dπrδϵλρ−dπrβη1λρ−bnπrβη1λρ.D5=− bcm2πδϵ−cdm2πδϵ+cdm2πβη1+bcm2nπβη1+dm2π2βη1+bm2nπ2βη1−bmπrδϵλ − dmπrδϵλ+dmπrβη1λ+bmnπrβη1λ.
Let
XL(p)=D1p4+D2p3+D3p2+D4p+D5=0.
Then
XL(0)=D5,XLmρ=−(b+d)mr2ϵλ2(mγ2η1+mγ1(γ2δ+η2)+(γ2δ+η2)ρ)ρ2<0.
The term D5 can be written as:
D5=mπ(bδϵ+dδϵ)(cm+rλ)R0L−1.
Therefore if R0L>1, then D5 > 0 and there exists p1∈(0,mρ) such that XL(p1) = 0. If R0L>1, then system (15)–(19) has an infected equilibrium E1 = (s1, y1, p1, x1), where
x1=λm−ρp1>0,y1=(c+rx1)p1π>0,s1=β(1+γ1p1)(1+γ2y1)δ(1+γ1p1)(1+γ2y1)+η1p1(1+γ2y1)+η2y1(1+γ1p1)>0,w1=(1−n)s1(p1η1+p1y1γ2η1+y1η2+p1y1γ1η2)(b+d)(1+p1γ1)(1+y1γ2)>0.
□

Theorem 3.

The equilibrium E0 of system (15)–(19), is globally asymptotically stable in Γ2 when if R0L≤1⁠.

Proof.
Define V0(s, w, y, p, x) as:
V0=s0Gss0+dbn+dw+b+dbn+dy+η1s0c+rx0p+rη1s0ρ(c+rx0)x0Gxx0.
Calculating dV0dt along system (15)–(19) we obtain
dV0dt=1−s0sβ−δs−η1sp1+γ1p−η2sy1+γ2y+dbn+d(1−n)η1sp1+γ1p+η2sy1+γ2y−(b+d)w + b+dbn+dnη1sp1+γ1p+η2sy1+γ2y+dw−ϵy+η1s0c+rx0πy−cp−rxp+rη1s0ρ(c+rx0)1−x0xλ+ρxp−mx=1−s0sβ−δs+η1sp1+γ1p+η2sy1+γ2y−b+dbn+dϵy+η1s0πc+rx0y−η1s0c+rx0cp−rη1s0ρ(c+rx0)p+rη1s0ρ(c+rx0)1−x0xλ−mx=1−s0sβ−δs−η1s0γ1p21+γ1p−η2s0γ2y21+γ2y+η2s0y−b+dbn+dϵy+η1s0πc+rx0y+rη1s0ρ(c+rx0)1−x0xλ−mx.
Then
(29)
Therefore if R0L≤1, then dV0dt≤0 for all s, w, y, p, x > 0 and dV0dt=0 at E0. LIP implies that E0 is globally asymptotically stable when R0L≤1.□

Theorem 4.

The equilibrium E1 of system (15)–(19), is globally asymptotically stable in Γ2° when if R0L>1⁠.

Proof.
Let
V1(s,w,y,p,x)=s1Gss1+dbn+dw1Gww1+b+dbn+dy1Gyy1+η1s1p1(1+γ1p1)πy1p1Gpp1+rη1s1p1ρ(1+γ1p1)πy1x1Gxx1.
Then
dV1dt=1−s1sβ−δs−η1sp1+γ1p−η2sy1+γ2y+dbn+d1−w1w(1−n)η1sp1+γ1p+η2sy1+γ2y−(b+d)w + b+dbn+d1−y1ynη1sp1+γ1p+η2sy1+γ2y+dw−ϵy+η1s1p1(1+γ1p1)πy11−p1pπy−cp−rxp+rη1s1p1ρ(1+γ1p1)πy11−x1xλ+ρxp−mx=1−s1sβ−δs+η1s1p1+γ1p+η2s1y1+γ2y−η1d(1−n)bn+dspw1(1+γ1p)w−η2d(1−n)bn+dsyw1(1+γ2y)w + d(b+d)bn+dw1−nη1(b+d)bn+dspy1(1+γ1p)y−nη2(b+d)bn+dsyy1(1+γ2y)y−d(b+d)bn+dwy1y−b+dbn+dϵy + b+dbn+dϵy1+η1s1p11+γ1p1yy1−η1s1p11+γ1p1yp1y1p−η1s1p1(1+γ1p1)πy1cp+η1s1p1(1+γ1p1)πy1cp1 + η1s1p1(1+γ1p1)πy1rxp1−rη1s1p1(1+γ1p1)πy1x1p+rη1s1p1ρ(1+γ1p1)πy11−x1xλ−mx.
Applying the equilibria conditions for E1
β=δs1+η1s1p11+γ1p1+η2s1y11+γ2y1,(b+d)w1=(1−n)η1s1p11+γ1p1+η2s1y11+γ2y1,b+dbn+dϵy1=η1s1p11+γ1p1+η2s1y11+γ2y1,cp1=πy1−rx1p1,λ=mx1−ρx1p1,
we get
dV1dt=− δ(s−s1)2s−η1s1p11+γ1p1γ1(p−p1)2(1+γ1p)(1+γ1p1)p1−η2s1y11+γ2y1γ2(y−y1)2(1+γ2y)(1+γ2y1)y1−η1s1p1rλ(1+γ1p1)πy1ρx1(x−x1)2x + d(1−n)bn+dη1s1p11+γ1p15−s1s−spw1(1+γ1p1)s1p1w(1+γ1p)−wy1w1y−yp1y1p−1+γ1p1+γ1p1+d(1−n)bn+dη2s1y11+γ2y14−s1s−syw1(1+γ2y1)s1y1w(1+γ2y)−y1wyw1−1+γ2y1+γ2y1 + n(b+d)bn+dη1s1p11+γ1p14−s1s−spy1(1+γ1p1)s1p1y(1+γ1p)−yp1y1p−1+γ1p1+γ1p1+n(b+d)bn+dη2s1y11+γ2y13−s1s−s(1+γ2y1)s1(1+γ2y)−1+γ2y1+γ2y1.
Using rule (14) we find that dV1dt≤0 for all s, w, y, p, x > 0 and dU1dt=0 at E1. The global stability of E1 is induced from LIP. □

In this section, we will perform numerical simulations for system (5)–(8) and (15)–(19) using MATLAB.

We will use the values of the parameters given in Table I. Moreover, we simulate system (5)–(8) with three different initial values as:

  • IV1:

    s(0) = 14.0, y(0) = 1.0, p(0) = 1.5, and x(0) = 1.5,

  • IV2:

    s(0) = 8.0, y(0) = 2.0, p(0) = 3.0, and x(0) = 4.0,

  • IV3:

    s(0) = 4.0, y(0) = 3.5, p(0) = 6.0, and x(0) = 7.0.

TABLE I.

The value of the parameters of model (5)–(8).

ParameterValueParameterValue
β 2 δ 0.1 
η1, η2 varied ϵ 0.5 
π 4 c 0.1 
r 0.5 λ 1.4 
m 1 ρ 0.2 
γ1 varied γ2 varied 
ParameterValueParameterValue
β 2 δ 0.1 
η1, η2 varied ϵ 0.5 
π 4 c 0.1 
r 0.5 λ 1.4 
m 1 ρ 0.2 
γ1 varied γ2 varied 
View Large

Then we consider two cases:

  • Set (I):

    We take η1 = η2 = 0.001 and γ1 = γ2 = 0.09. The value of R0 is computed as R0=0.2400<1⁠. Figure 1 shows that, the solutions s(t) and x(t) with the initial values IV1-IV3 reach the values s0=βδ=20 and x0=λm=1.4⁠, respectively, while y(t) and p(t) approach zero. This shows that, E0 is globally asymptotically stable which agrees with Theorem 1.

  • Set (II):

    We take η1 = η2 = 0.05 and γ1 = γ2 = 0.09. Then, we calculate R0=12.0>1. We compute the equilibrium E1 = (5.61, 2.87, 3.80, 5.84). Figure 1 shows that when R0>1⁠, the states of the system tend to E1 for all the three initial values IV1-IV3. This confirms that the validity of Theorem 2.

FIG. 1.
FIG. 1. The simulation of trajectories of system (5)–(8). (a) Uninfected monocytes. (b) Infected monocytes. (c) Free CHIKV particles. (d) Antibodies.
View largeDownload slide

The simulation of trajectories of system (5)–(8). (a) Uninfected monocytes. (b) Infected monocytes. (c) Free CHIKV particles. (d) Antibodies.

FIG. 1.
FIG. 1. The simulation of trajectories of system (5)–(8). (a) Uninfected monocytes. (b) Infected monocytes. (c) Free CHIKV particles. (d) Antibodies.
View largeDownload slide

The simulation of trajectories of system (5)–(8). (a) Uninfected monocytes. (b) Infected monocytes. (c) Free CHIKV particles. (d) Antibodies.

Close modal

1. Effect of the saturation parameters for system (5)–(8)

We fix the value η1 = η2 = 0.05. Let us consider γ = γ1 = γ2 and the initial s(0) = 12, y(0) = 1.0, p(0) = 3.0, and x(0) = 4.0. The evolution of the system’s states with different values of γ is shown in Figure 2. It is shown that, as γ is increased, the concentration of the uninfected monocytes are increased while the concentrations of the other compartments are decreased.

FIG. 2.
FIG. 2. The simulation of trajectories of system (5)–(8) with different values of γ. (a) Uninfected monocytes. (b) Infected monocytes. (c) Free CHIKV particles. (d) Antibodies.
View largeDownload slide

The simulation of trajectories of system (5)–(8) with different values of γ. (a) Uninfected monocytes. (b) Infected monocytes. (c) Free CHIKV particles. (d) Antibodies.

FIG. 2.
FIG. 2. The simulation of trajectories of system (5)–(8) with different values of γ. (a) Uninfected monocytes. (b) Infected monocytes. (c) Free CHIKV particles. (d) Antibodies.
View largeDownload slide

The simulation of trajectories of system (5)–(8) with different values of γ. (a) Uninfected monocytes. (b) Infected monocytes. (c) Free CHIKV particles. (d) Antibodies.

Close modal

In this case we select n = 0.5, d = 0.5 and b = 0.2. The other parameters are the same as given in Table I. Let us choose three initial values as:

  • IV4:

    s(0) = 14.0, w(0) = 1.0, y(0) = 1.0, p(0) = 1.5, and x(0) = 1.5,

  • IV5:

    s(0) = 8.0, w(0) = 2.0, y(0) = 2.0, p(0) = 3.0, and x(0) = 4.0,

  • IV6:

    s(0) = 4.0, w(0) = 4.0, y(0) = 3.5, p(0) = 6.0, and x(0) = 7.0.

Consider the two cases:

  • Set (I):

    Take η1 = η2 = 0.001 and γ1 = γ2 = 0.09, then R0=0.2057<1. From Figure 3 we can see that, the concentrations of the uninfected monocytes and antibodies return to their values s0=βδ=20 and x0=λm=1.4⁠, respectively. Figure 3 establishes that, the solutions s(t) and x(t) with the initial values IV4-IV6 reach the values s0=βδ=20 and x0=λm=1.4⁠, respectively, while w(t), y(t) and p(t) approach zero. This confirms the global asymptotic stability result of Theorem 3.

  • Set (II):

    We take η1 = η2 = 0.05 and γ1 = γ2 = 0.09. Then, we calculate R0=10.2857>1, and E1 = (5.95, 1.00, 2.40, 3.62, 5.10). Figure 3 shows that when R0>1⁠, the solutions of the system starting from the initials IV4-IV6 will tend to E1. This displays the agreement between the numerical and theoretical results of Theorem 4.

FIG. 3.
FIG. 3. The simulation of trajectories of system (15)–(19). (a) Uninfected monocytes. (b) Latently infected monocytes. (c) Actively infected monocytes. (d) Free CHIKV particles. (e) Antibodies.
View largeDownload slide

The simulation of trajectories of system (15)–(19). (a) Uninfected monocytes. (b) Latently infected monocytes. (c) Actively infected monocytes. (d) Free CHIKV particles. (e) Antibodies.

FIG. 3.
FIG. 3. The simulation of trajectories of system (15)–(19). (a) Uninfected monocytes. (b) Latently infected monocytes. (c) Actively infected monocytes. (d) Free CHIKV particles. (e) Antibodies.
View largeDownload slide

The simulation of trajectories of system (15)–(19). (a) Uninfected monocytes. (b) Latently infected monocytes. (c) Actively infected monocytes. (d) Free CHIKV particles. (e) Antibodies.

Close modal

This project was funded by the Deanship of Scientific Research (DSR), King Abdulaziz University, Jeddah, Saudi Arabia under grant no. (KEP-MSc-13-130-38). The authors, therefore, acknowledge with thanks DSR technical and financial support.

1.
Y.
 
Dumont
 and
F.
 
Chiroleu
, “
Vector control for the chikungunya disease
,”
Mathematical Biosciences and Engineering
 
7
,
313
–
345
(
2010
).
Google Scholar
PubMed
 
2.
Y.
 
Dumont
 and
J. M.
 
Tchuenche
, “
Mathematical studies on the sterile insect technique for the chikungunya disease and Aedes albopictus
,”
Journal of Mathematical Biology
 
65
(
5
),
809
–
854
(
2012
).
https://doi.org/10.1007/s00285-011-0477-6
Google Scholar
Crossref
Search ADS
PubMed
 
3.
Y.
 
Dumont
,
F.
 
Chiroleu
, and
C.
 
Domerg
, “
On a temporal model for the chikungunya disease: Modeling, theory and numerics
,”
Mathematical Biosciences
 
213
,
80
–
91
(
2008
).
https://doi.org/10.1016/j.mbs.2008.02.008
Google Scholar
Crossref
Search ADS
PubMed
 
4.
D.
 
Moulay
,
M.
 
Aziz-Alaoui
, and
M.
 
Cadivel
, “
The chikungunya disease: Modeling, vector and transmission global dynamics
,”
Mathematical Biosciences
 
229
,
50
–
63
(
2011
).
https://doi.org/10.1016/j.mbs.2010.10.008
Google Scholar
Crossref
Search ADS
PubMed
 
5.
D.
 
Moulay
,
M.
 
Aziz-Alaoui
, and
H. D.
 
Kwon
, “
Optimal control of chikungunya disease: Larvae reduction, treatment and prevention
,”
Mathematical Biosciences and Engineering
 
9
,
369
–
392
(
2012
).
https://doi.org/10.3934/mbe.2012.9.369
Google Scholar
Crossref
Search ADS
PubMed
 
6.
C. A.
 
Manore
,
K. S.
 
Hickmann
,
S.
 
Xu
,
H. J.
 
Wearing
, and
J. M.
 
Hyman
, “
Comparing dengue and chikungunya emergence and endemic transmission in A. aegypti and A. albopictus
,”
Journal of Theoretical Biology
 
356
,
174
–
191
(
2014
).
https://doi.org/10.1016/j.jtbi.2014.04.033
Google Scholar
Crossref
Search ADS
PubMed
 
7.
L.
 
Yakob
 and
A. C.
 
Clements
, “
A mathematical model of chikungunya dynamics and control: The major epidemic on Reunion Island
,”
PLoS One
 
8
,
e57448
(
2013
).
https://doi.org/10.1371/journal.pone.0057448
Google Scholar
Crossref
Search ADS
PubMed
 
8.
X.
 
Liu
 and
P.
 
Stechlinski
, “
Application of control strategies to a seasonal model of chikungunya disease
,”
Applied Mathematical Modelling
 
39
,
3194
–
3220
(
2015
).
https://doi.org/10.1016/j.apm.2014.10.035
Google Scholar
Crossref
Search ADS
 
9.
Y.
 
Wang
 and
X.
 
Liu
, “
Stability and Hopf bifurcation of a within-host chikungunya virus infection model with two delays
,”
Mathematics and Computers in Simulation
 
138
,
31
–
48
(
2017
).
https://doi.org/10.1016/j.matcom.2016.12.011
Google Scholar
Crossref
Search ADS
 
10.
A. M.
 
Elaiw
,
T. O.
 
Alade
, and
S. M.
 
Alsulami
, “
Analysis of latent CHIKV dynamics models with general incidence rate and time delays
,”
Journal of Biological Dynamics
 
12
(
1
),
700
–
730
(
2018
).
https://doi.org/10.1080/17513758.2018.1503349
Google Scholar
Crossref
Search ADS
PubMed
 
11.
A. M.
 
Elaiw
,
T. O.
 
Alade
, and
S. M.
 
Alsulami
, “
Analysis of within-host CHIKV dynamics models with general incidence rate
,”
International Journal of Biomathematics
 
11
(
5
),
1850062
(
2018
).
https://doi.org/10.1142/s1793524518500626
Google Scholar
Crossref
Search ADS
 
12.
K. M.
 
Long
 and
M. T.
 
Heise
, “
Protective and pathogenic responses to chikungunya virus infection
,”
Curr Trop Med Rep
 
2
(
1
),
13
–
21
(
2015
).
https://doi.org/10.1007/s40475-015-0037-z
Google Scholar
Crossref
Search ADS
PubMed
 
13.
J.
 
Wang
,
J.
 
Lang
, and
X.
 
Zou
, “
Analysis of an age structured HIV infection model with virus-to-cell infection and cell-to-cell transmission
,”
Nonlinear Analysis: Real World Applications
 
34
,
75
–
96
(
2017
).
https://doi.org/10.1016/j.nonrwa.2016.08.001
Google Scholar
Crossref
Search ADS
 
14.
F.
 
Li
 and
J.
 
Wang
, “
Analysis of an HIV infection model with logistic target cell growth and cell-to-cell transmission
,”
Chaos, Solitons and Fractals
 
81
,
136
–
145
(
2015
).
https://doi.org/10.1016/j.chaos.2015.09.003
Google Scholar
Crossref
Search ADS
 
15.
A. M.
 
Elaiw
,
A.
 
Raezah
, and
A.
 
Alofi
, “
Stability of a general delayed virus dynamics model with humoral immunity and cellular infection
,”
AIP Advances
 
7
(
6
),
065210
(
2017
).
https://doi.org/10.1063/1.4989569
Google Scholar
Crossref
Search ADS
 
16.
A. M.
 
Elaiw
,
A.
 
Raezah
, and
A. S.
 
Alofi
, “
Effect of humoral immunity on HIV-1 dynamics with virus-to-target and infected-to-target infections
,”
AIP Advances
 
6
(
8
),
085204
(
2016
).
https://doi.org/10.1063/1.4960987
Google Scholar
Crossref
Search ADS
 
17.
A. M.
 
Elaiw
,
A. A.
 
Almatrafi
, and
A. D.
 
Hobiny
, “
Effect of antibodies on pathogen dynamics with delays and two routes of infection
,”
AIP Advances
 
8
,
065104
(
2018
).
https://doi.org/10.1063/1.5029483
Google Scholar
Crossref
Search ADS
 
18.
H.
 
Shu
,
Y.
 
Chen
, and
L.
 
Wang
, “
Impacts of the cell-free and cell-to-cell infection modes on viral dynamics
,”
Journal of Dynamics and Differential Equations
 
30
(
4
),
1817
–
1836
(
2018
).
https://doi.org/10.1007/s10884-017-9622-2
Google Scholar
Crossref
Search ADS
 
19.
X.
 
Lai
 and
X.
 
Zou
, “
Modeling cell-to-cell spread of HIV-1 with logistic target cell growth
,”
Journal of Mathematical Analysis and Applications
 
426
,
563
–
584
(
2015
).
https://doi.org/10.1016/j.jmaa.2014.10.086
Google Scholar
Crossref
Search ADS
 
20.
A. M.
 
Elaiw
,
A. A.
 
Almatrafi
,
A. D.
 
Hobiny
, and
I. A.
 
Abbas
, “
Stability of latent pathogen infection model with CTL immune response and saturated cellular infection
,”
AIP Advances
 
8
,
125021
(
2018
).
https://doi.org/10.1063/1.5079402
Google Scholar
Crossref
Search ADS
 
21.
A. M.
 
Elaiw
,
A. A.
 
Raezah
, and
K.
 
Hattaf
, “
Stability of HIV-1 infection with saturated virus-target and infected-target incidences and CTL immune response
,”
International Journal of Biomathematics
 
10
(
5
),
1750070
(
2017
).
Google Scholar
Crossref
Search ADS
 
22.
A. M.
 
Elaiw
 and
A. A.
 
Raezah
, “
Stability of general virus dynamics models with both cellular and viral infections and delays
,”
Mathematical Methods in the Applied Sciences
 
40
(
16
),
5863
–
5880
(
2017
).
https://doi.org/10.1002/mma.4436
Google Scholar
Crossref
Search ADS
 
23.
X.
 
Lai
 and
X.
 
Zou
, “
Modelling HIV-1 virus dynamics with both virus-to-cell infection and cell-to-cell transmission
,”
SIAM Journal of Applied Mathematics
 
74
,
898
–
917
(
2014
).
https://doi.org/10.1137/130930145
Google Scholar
Crossref
Search ADS
 
24.
Y.
 
Yang
,
L.
 
Zou
, and
S.
 
Ruanc
, “
Global dynamics of a delayed within-host viral infection model with both virus-to-cell and cell-to-cell transmissions
,”
Mathematical Biosciences
 
270
,
183
–
191
(
2015
).
https://doi.org/10.1016/j.mbs.2015.05.001
Google Scholar
Crossref
Search ADS
PubMed
 
25.
A.
 
Korobeinikov
, “
Global properties of basic virus dynamics models
,”
Bulletin of Mathematical Biology
 
66
(
4
),
879
–
883
, (
2004
).
https://doi.org/10.1016/j.bulm.2004.02.001
Google Scholar
Crossref
Search ADS
PubMed
 
26.
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 Advances
 
8
,
025323
(
2018
).
https://doi.org/10.1063/1.5023752
Google Scholar
Crossref
Search ADS
 
27.
A. M.
 
Elaiw
 and
S. A.
 
Azoz
, “
Global properties of a class of HIV infection models with Beddington-DeAngelis functional response
,”
Mathematical Methods in the Applied Sciences
 
36
,
383
–
394
(
2013
).
https://doi.org/10.1002/mma.2596
Google Scholar
Crossref
Search ADS
 
28.
A. M.
 
Elaiw
, “
Global properties of a class of HIV models
,”
Nonlinear Analysis: Real World Applications
 
11
,
2253
–
2263
(
2010
).
https://doi.org/10.1016/j.nonrwa.2009.07.001
Google Scholar
Crossref
Search ADS
 
29.
H.
 
Shu
,
L.
 
Wang
, and
J.
 
Watmough
, “
Global stability of a nonlinear viral infection model with infinitely distributed intracellular delays and CTL immune responses
,”
SIAM Journal of Applied Mathematics
 
73
(
3
),
1280
–
1302
(
2013
).
https://doi.org/10.1137/120896463
Google Scholar
Crossref
Search ADS
 
30.
A. M.
 
Elaiw
 and
N. H.
 
AlShamrani
, “
Global stability of humoral immunity virus dynamics models with nonlinear infection rate and removal
,”
Nonlinear Analysis: Real World Applications
 
26
,
161
–
190
(
2015
).
https://doi.org/10.1016/j.nonrwa.2015.05.007
Google Scholar
Crossref
Search ADS
 
31.
A. M.
 
Elaiw
 and
N. H.
 
AlShamrani
, “
Stability of a general delay-distributed virus dynamics model with multi-staged infected progression and immune response
,”
Mathematical Methods in the Applied Sciences
 
40
(
3
),
699
–
719
(
2017
).
https://doi.org/10.1002/mma.4002
Google Scholar
Crossref
Search ADS
 
32.
A. M.
 
Elaiw
 and
N. A.
 
Almuallem
, “
Global dynamics of delay-distributed HIV infection models with differential drug efficacy in cocirculating target cells
,”
Mathematical Methods in the Applied Sciences
 
39
,
4
–
31
(
2016
).
https://doi.org/10.1002/mma.3453
Google Scholar
Crossref
Search ADS
 
33.
G.
 
Huang
,
Y.
 
Takeuchi
, and
W.
 
Ma
, “
Lyapunov functionals for delay differential equations model of viral infections
,”
SIAM Journal of applied Mathematics
 
70
(
7
),
2693
–
2708
(
2010
).
https://doi.org/10.1137/090780821
Google Scholar
Crossref
Search ADS
 
34.
A. M.
 
Elaiw
,
E. Kh.
 
Elnahary
, and
A. A.
 
Raezah
, “
Effect of cellular reservoirs and delays on the global dynamics of HIV
,”
Advances in Difference Equations
 
2018
,
85
.
https://doi.org/10.1186/s13662-018-1523-0
Crossref
Search ADS
 
35.
A. M.
 
Elaiw
 and
N. H.
 
AlShamrani
, “
Stability of an adaptive immunity pathogen dynamics model with latency and multiple delays
,”
Mathematical Methods in the Applied Sciences
 
36
,
125
–
142
(
2018
).
Google Scholar
 
36.
A. M.
 
Elaiw
,
I. A.
 
Hassanien
, and
S. A.
 
Azoz
, “
Global stability of HIV infection models with intracellular delays
,”
Journal of the Korean Mathematical Society
 
49
(
4
),
779
–
794
(
2012
).
https://doi.org/10.4134/jkms.2012.49.4.779
Google Scholar
Crossref
Search ADS
 
37.
C.
 
Lv
,
L.
 
Huang
, and
Z.
 
Yuan
, “
Global stability for an HIV-1 infection model with Beddington–DeAngelis incidence rate and CTL immune response
,”
Communications in Nonlinear Science and Numerical Simulation
 
19
,
121
–
127
(
2014
).
https://doi.org/10.1016/j.cnsns.2013.06.025
Google Scholar
Crossref
Search ADS
 
38.
A. M.
 
Elaiw
,
A. A.
 
Raezah
, and
S. A.
 
Azoz
, “
Stability of delayed HIV dynamics models with two latent reservoirs and immune impairment
,”
Advances in Difference Equations
 
2018
,
414
.
https://doi.org/10.1186/s13662-018-1869-3
Crossref
Search ADS
 
39.
J.
 
Xu
,
Y.
 
Zhou
,
Y.
 
Li
, and
Y.
 
Yang
, “
Global dynamics of a intracellular infection model with delays and humoral immunity
,”
Mathematical Methods in the Applied Sciences
 
39
(
18
),
5427
–
5435
(
2016
).
https://doi.org/10.1002/mma.3927
Google Scholar
Crossref
Search ADS
 
40.
A. M.
 
Elaiw
, “
Global dynamics of an HIV infection model with two classes of target cells and distributed delays
,”
Discrete Dynamics in Nature and Society
 
2012
,
253703
.
41.
A. M.
 
Elaiw
, “
Global properties of a class of virus infection models with multitarget cells
,”
Nonlinear Dynamics
 
69
(
1-2
),
423
–
435
(
2012
).
https://doi.org/10.1007/s11071-011-0275-0
Google Scholar
Crossref
Search ADS
 
42.
L.
 
Gibelli
,
A.
 
Elaiw
,
M. A.
 
Alghamdi
, and
A. M.
 
Althiabi
, “
Heterogeneous population dynamics of active particles: Progression, mutations, and selection dynamics
,”
Mathematical Models and Methods in Applied Sciences
 
27
,
617
–
640
(
2017
).
https://doi.org/10.1142/s0218202517500117
Google Scholar
Crossref
Search ADS
 
© 2019 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
2019
Author(s)