A numerical method for a delayed viral infection model with general incidence rate

Journal of King Saud University – Science (2015) xxx, xxx–xxx

King Saud University

Journal of King Saud University – Science www.ksu.edu.sa www.sciencedirect.com

ORIGINAL ARTICLE

A numerical method for a delayed viral infection model with general incidence rate Khalid Hattaf a,b,*, Noura Yousfi b a

Centre Re´gional des Me´tiers de l’Education et de la Formation (CRMEF), 20340 Derb Ghalef, Casablanca, Morocco Department of Mathematics and Computer Science, Faculty of Sciences Ben M’sik, Hassan II University, P.O Box 7955 Sidi Othman, Casablanca, Morocco b

Received 26 March 2015; accepted 16 October 2015

KEYWORDS Delay difference equations; General incidence rate; Mixed Euler method; Global stability

Abstract In this paper, we construct a numerical method for a delayed viral infection model with general incidence rate. We prove that the obtained discrete model has the same dynamics as the corresponding continuous model, such as positivity, boundedness and global behaviors of solutions with no restriction on the time step size. Furthermore, numerical simulations are given to illustrate and confirm our main analytical results. Ó 2015 The Authors. Production and hosting by Elsevier B.V. on behalf of King Saud University. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

1. Introduction In recent years, several authors are interested in the study of the dynamics of viral infections by proposing the continuous mathematical models with delays and different forms of incidence rate, such as mass action process (Zhu and Zou, 2008; Li and Shu, 2010; Hattaf and Yousfi, 2011; Vargas-De-Leo´n, 2012), standard incidence function (Gourley et al., 2008; Eikenberry et al., 2009; Tian and Xu, 2010), saturated mass action (Li and Ma, 2007; Xu, 2011), Beddington–DeAngelis functional response (Huang et al., 2011; Xiang et al., 2013) and Crowley–Martin functional response (Zhou and Cui,

* Corresponding author. E-mail address: [email protected] (K. Hattaf). Peer review under responsibility of King Saud University.

Production and hosting by Elsevier

2011). In 2013, the authors (Hattaf et al. (2013)) have generalized all previous forms by proposing the following model: 8 _ ¼ k  dxðtÞ  fðxðtÞ; yðtÞ; vðtÞÞvðtÞ; xðtÞ > > > < yðtÞ _ ¼ fðxðt  s1 Þ; yðt  s1 Þ; vðt  s1 ÞÞvðt  s1 Þea1 s1  ayðtÞ; > _ ¼ kyðt  s2 Þea2 s2  uvðtÞ; vðtÞ > > : ð1Þ where xðtÞ; yðtÞ and vðtÞ denote the concentration of uninfected cells, infected cells and free virus particles at time t, respectively. The parameter k is the recruited rate of uninfected cells, k is the production rate of free virus by infected cells, d, a and u are, respectively, the death rates of uninfected cells, infected cells and free virus. The first delay s1 represents the time needed for infected cells to produce virions after viral entry and the factor ea1 s1 accounts for the probability of surviving from time t  s1 to time t, where a1 is the death rate for infected but not yet virus-producing cells. The second delay s2 denotes the time necessary for the newly produced virions to become mature and then infectious particles. The probability

http://dx.doi.org/10.1016/j.jksus.2015.10.003 1018-3647 Ó 2015 The Authors. Production and hosting by Elsevier B.V. on behalf of King Saud University. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Please cite this article in press as: Hattaf, K., Yousfi, N. A numerical method for a delayed viral infection model with general incidence rate. Journal of King Saud University – Science (2015), http://dx.doi.org/10.1016/j.jksus.2015.10.003

2

K. Hattaf, N. Yousfi

of survival of immature virions is given by ea2 s2 and the average life time of an immature virus is given by a12 . The incidence function fðx; y; vÞ is assumed to be continuously differentiable in the interior of R3þ and satisfies the three fundamental hypotheses given in Hattaf et al. (2012) and used in Hattaf et al. (2014), Hattaf and Yousfi (2014) and Wang et al. (2013), that are: fð0; y; vÞ ¼ 0; for all y P 0 and v P 0; @f ðx; y; vÞ > 0; for all x > 0; y P 0 and v P 0; @x @f @f ðx; y; vÞ 6 0 and ðx; y; vÞ 6 0; for all x P 0; @y @v y P 0 and v P 0:

ðH1Þ ðH2Þ

2. Delayed discrete model and preliminaries Let h be a time step size. Assume that there exist integers ðm1 ; m2 Þ 2 N with s1 ¼ m1 h and s2 ¼ m2 h. The grids points are tn ¼ nh for n 2 N. By applying both forward and backward Euler methods and using the approximations xðtn Þ  xn ; yðtn Þ  yn and vðtn Þ  vn , we obtain the following delayed discrete model 8 xnþ1 ¼ xn þ hðk  dxnþ1  fðxnþ1 ; yn ; vn Þvn Þ; > > <   ynþ1 ¼ yn þ h fðxnm1 þ1 ; ynm1 ; vnm1 Þvnm1 ea1 s1  aynþ1 ; > >   : vnþ1 ¼ vn þ h kynm2 þ1 ea2 s2  uvnþ1 : ð2Þ

ðH3Þ

From the biological point of view, the three hypotheses are reasonable. Indeed, the first means that the incidence rate is equal to zero if there are no susceptible cells. The second one signifies that the incidence rate is increasing when the numbers of infected cells and virus are constant and the number of susceptible cells increases. Hence, the second hypothesis means the more the amount of susceptible cells, the more the average number of cells which are infected by each virus in the unit time will occur. Similarly, the third assumption means the more the amount of infected cells or virus, the less the average number of cells which are infected by each virus in the unit time will be. On the other hand, the infectious process is not instantaneous. For this reason, we choose to use delay differential equations in order to take into account the time needed for infected cells to produce new virions after viral entry and the time necessary for the newly produced virions to become mature and infectious. In Hattaf et al. (2013), Hattaf et al. proved the positivity and boundedness of solutions. Also, they identified the basic reproduction number of model (1) as follows   k k R0 ¼ f ; 0; 0 ea1 s1 a2 s2 : au d Moreover, they established the global stability of equilibria. In reality, scientists often collect the data and analyze the results at discrete times. In addition, the numerical simulations of continuous models are obtained by discretizing these models. For these reasons, we will discretize the model (1) by using ‘mixed’ Euler method which is a mixture of both forward and backward Euler methods. Furthermore, we will show that the discrete model obtained by the mixed Euler method maintains essential dynamical properties, such as positivity, boundedness and global behaviors of solutions with no restriction on the time step size. The remainder of this paper is organized as follows. In the next section, we introduce our discrete virus dynamics model with general and two delays, and establish some preliminary results. The stability of the disease-free equilibrium and the chronic infection equilibrium of the new delayed discrete model is investigated in Sections 3 and 4. Numerical simulations are given to verify the main theoretical results in Section 5. The paper ends with a conclusion in Section 6.

The sequences fxn g; fyn g and fvn g represent the concentrations of cells and free virus at time n. Biologically, these concentrations are positives and bounded. For this, we assume that the initial values of model (2) satisfy: xðsÞ P 0;

yðsÞ P 0;

vðsÞ P 0;

for all

s ¼ m; m þ 1; . . . ; 0;

ð3Þ

where m ¼ maxðm1 ; m2 Þ. The following result establishes the positivity and boundedness of solutions of the discrete model (2). Proposition 2.1. All solutions of system (2) subject to condition (3) remain nonnegative and bounded for all n 2 N. Proof. First, we prove the positivity of solutions by using mathematical induction. When n ¼ 0, we have ð1 þ hdÞx1 þ hfðx1 ; y0 ; v0 Þv0 ¼ x0 þ hk: Let gðxÞ ¼ ð1 þ hdÞx þ hfðx; y0 ; v0 Þv0  x0  hk. Clearly, g(x) is a continuous function for x; gðRÞ ¼ R and g0 ðxÞ ¼ 1 þ hd þ hv0

@f ðx; y0 ; v0 Þ > 0: @x

Then gðxÞ is monotonically increasing on R. Hence, the equation gðxÞ ¼ 0 has a unique solution on R. Since     x0 þ hk x0 þ hk gð0Þ ¼ x0  hk < 0 and g ¼ hf ;y ;v0 v0 > 0; 1 þ hd 1 þ hd 0   0 þhk we have x 2 0; x1þhd . Hence, x ¼ x1 > 0. From (2), we obtain y0 þ hfðxm1 þ1 ; ym1 ; vm1 Þvm1 ea1 s1 ; 1 þ ah v0 þ hky1 v1 ¼ : 1 þ uh

y1 ¼

Then y1 P 0 and v1 P 0. Therefore, by using the induction, we get xn P 0; yn P 0 and vn P 0 for all n P 0. This proves the positivity of solutions. Next, we prove the boundedness of solutions. Let Tn ¼ xn þ yn þ h

n1 X

fðxjþ1 ; yj ; vj Þvj eha1 ðnjÞ :

j¼nm1

Please cite this article in press as: Hattaf, K., Yousfi, N. A numerical method for a delayed viral infection model with general incidence rate. Journal of King Saud University – Science (2015), http://dx.doi.org/10.1016/j.jksus.2015.10.003

A numerical method for a delayed viral infection model with general incidence rate Then, Tnþ1  Tn ¼ hðk  dxnþ1  fðxnþ1 ; yn ; vn Þvn Þ þ hðfðxnm1 þ1 ; ynm1 ; vnm1 Þvnm1 ea1 s1  aynþ1 Þ n X þh fðxjþ1 ; yj ; vj Þvj eha1 ðnþ1jÞ j¼nm1 þ1

h

n1 X

fðxjþ1 ; yj ; vj Þvj eha1 ðnjÞ

¼ hðk  dxnþ1  aynþ1 Þ n X fðxjþ1 ; yj ; vj Þvj eha1 ðnþ1jÞ þ hð1  eha1 Þ j¼nm1 þ1

Tnþ1 6

3. Stability of the disease-free equilibrium

Theorem 3.1. The disease-free equilibrium Ef of system (2) is globally asymptotically stable whenever R0 6 1, and unstable otherwise.

6 hðk  dTnþ1 Þ

Proof. We construct a discrete Lyapunov functional as follows

ha1 1

h

2. If R0 > 1, then the system (2) has a unique chronic infection equilibrium of the form E ðx ; y  ; v Þ besides Ef , where ky  kdx  x 2 ðdk ; 0Þ; y  ¼ ae a1 s1 and v ¼ uea2 s2 .

In this section, we investigate the stability of the disease-free equilibrium.

j¼nm1

where d ¼ minfd; a; e

3

g. Hence,

Z

fðx0 ; 0; 0Þ a ds þ ea1 s1 yn þ ea1 s1 þa2 s2 ð1 þ uhÞvn fðs; 0; 0Þ k x0 n1 n1 X X þh fðxjþ1 ; yj ; vj Þvj þ aea1 s1 h yjþ1 ;

1 hk Tn þ : 1 þ hd 1 þ hd

Ln ¼ xn  x0 

By using the induction, we get the following inequality   n  n 1 k 1 Tn 6 : T0 þ 1  1 þ hd d 1 þ hd

xn

j¼nm1

j¼nm2

Then,

Calculating the first difference of Ln along the positive solution of system (2), we have

k lim sup Tn 6 : d n!þ1

DLn ¼ Lnþ1  Ln

This implies that fTn g is bounded. Therefore, fxn g and fyn g are also bounded. By the third equation of (2), we obtain vnþ1 ¼

1 hkea2 s2 vn þ y : 1 þ hu 1 þ hu nm2 þ1

As fyn g is bounded, then there exists a M such that yn 6 M for all n 2 fm2 ; m2 þ 1; . . . ; 0; 1; . . .g. Thus, vnþ1 6

1 hkea2 s2 vn þ M: 1 þ hu 1 þ hu

By induction, we get  n  n  1 kM a2 s2 1 kM a2 s2 e e v0 þ 1 ; 6 v0 þ vn 6 1 þ hu u 1 þ hu u Therefore, fvn g is bounded. This completes the proof.

h

If in addition we assume that x0 > 0; y0 > 0 and v0 > 0, it is easy to get the following result. Remark 2.2. If x0 > 0; y0 > 0 and v0 > 0, then all solutions of system (2) subject to condition (3) are positive for any n > 0. In this case, the mixed Euler numerical scheme for the system (1) is called unconditionally positive. Next, we find the steady states of system (2). By a simple computation, it is easy to see that the system (2) has the same steady states as those of the corresponding continuous system (1). k

 k Theorem 2.3. Let us define R0 ¼ au f d ; 0; 0 ea1 s1 a2 s2 . 1. If R0 6 1, then the system  (2) has a unique disease-free equilibrium of the form Ef dk ; 0; 0 .

¼ xnþ1  xn 

Z

xnþ1

xn

fðx0 ; 0; 0Þ ds þ ea1 s1 ðynþ1  yn Þ fðs; 0; 0Þ

a þ ð1 þ uhÞea1 s1 þa2 s2 ðvnþ1  vn Þ k  þ h fðxnþ1 ; yn ; vn Þvn  fðxnm1 þ1 ; ynm1 ; vnm1 Þvnm1   þ aea1 s1 h ynþ1  ynm2 þ1   fðx0 ; 0; 0Þ ðxnþ1  xn Þ þ ea1 s1 ðynþ1  yn Þ 6 1 fðxnþ1 ; 0; 0Þ a þ ð1 þ uhÞea1 s1 þa2 s2 ðvnþ1  vn Þ k  þ h fðxnþ1 ; yn ; vn Þvn  fðxnm1 þ1 ; ynm1 ; vnm1 Þvnm1   þ aea1 s1 h ynþ1  ynm2 þ1 : Using the equality k ¼ dx0 , we get    xnþ1  fðx0 ; 0; 0Þ 0 1 DLn 6 hdx 1  0 fðxnþ1 ; 0; 0Þ x   auh a1 s1 þa2 s2 fðxnþ1 ; yn ; vn Þ e R0  1 vn þ k fðxnþ1 ; 0; 0Þ     xnþ1 fðx0 ; 0; 0Þ 1 6 hdx0 1  0 fðxnþ1 ; 0; 0Þ x auh a1 s1 þa2 s2 e þ ðR0  1Þvn : k Since fðx; y; vÞ is strictly monotonically increasing with respect to x, we have    xnþ1  fðx0 ; 0; 0Þ 1 0 1 6 0: x fðxnþ1 ; 0; 0Þ Since R0  1, we have DLn ¼ Lnþ1  Ln 6 0 for all n P 0. Then, fLn g is a monotonically decreasing sequence. Since Ln P 0, we have limn!þ1 Ln P 0. Hence limn!þ1 ðLnþ1  Ln Þ ¼ 0, from which we obtain limn!þ1 xnþ1 ¼ x0 and limn!þ1 ððR0  1Þvn Þ ¼ 0. We discuss two cases:

Please cite this article in press as: Hattaf, K., Yousfi, N. A numerical method for a delayed viral infection model with general incidence rate. Journal of King Saud University – Science (2015), http://dx.doi.org/10.1016/j.jksus.2015.10.003

4

K. Hattaf, N. Yousfi

 If R0 < 1, then limn!þ1 vn ¼ 0. By the third equation of (2), we get limn!þ1 y n ¼ 0.  If R0 ¼ 1. Using limn!þ1 xn ¼ x0 and the first equation of (2), it is not hard to have limn!þ1 vn ¼ 0. By the above discussion, we deduce that Ef is globally asymptotically stable if R0  1. Now, we prove that the disease-free equilibrium Ef is unstable when R0 > 1. Calculating the linearization system of model (2) at equilibrium Ef , we get a new system of the form 8   k > < Xnþ1 ¼ Xn þ hdXnþ1  fðd ; 0; 0ÞVn ;  ð4Þ Ynþ1 ¼ Yn þ h ea1 s1 fðkd ; 0; 0ÞVnm1  aYnþ1 ; > : a2 s2 Ynm2 þ1  uVnþ1 Þ; Vnþ1 ¼ Vn þ hðke where

Xn ¼ xn  kd ; Yn ¼ yn

Vn ¼ vn .

and

Let

Zn ¼ ðXn ; Yn ; Vn ÞT . Then, system (4) is equivalent to Znþ1 ¼ AZn þ BZnm1 þ CZnm2 þ DZnm1 m2 ; where 0

1

0

hfð Þ 1 0 0 B 1þhd 0  1þhd C B B C B 1 A¼B 0 C; B ¼ @ 0 0 0 1þha @ A 1 0 0 0 0 1þhu 0 0 1 0 0 0 0 B B C B0 0 0C C¼B @0 A; D ¼ @ hk a2 s2 0 0 ð1þhaÞð1þhuÞ e 0 k;0;0 d

0

ð5Þ

hfðkd;0;0Þ a s C e 1 1C A; 1þha

0 0

0

0

0

h2 kfðkd;0;0Þ ð1þhaÞð1þhuÞ

Remark 4.1. The assumption (H4) is satisfied by various types of the incidence rate including the mass action when bx fðx; y; vÞ ¼ bx, the saturation incidence when fðx; y; vÞ ¼ 1þav , the incidence function was used in Zhuo (2012) and Sun and bx Min (2014) when fðx; y; vÞ ¼ xþv , Beddington–DeAngelis response when fðx; y; vÞ ¼ 1þa1bx xþa2 v, Crowley-Martin response when fðx; y; vÞ ¼ 1þa1 xþabx and the more generalized 2 vþa1 a2 xv response introduced by Hattaf et al. (see Section 5 in Hattaf bx et al. (2013)) when fðx; y; vÞ ¼ 1þa1 xþa , where b is a positive 2 vþa3 xv constant rate describing the infection process, a; a1 ; a2 and a3 are nonnegative constants. Further, the fourth hypothesis given in Wang et al. (2013) on the incidence function depending only of x and v is a particular case of the assumption (H4). The following theorem establish the global stability of E . Theorem 4.2. Assume R0 > 1 and (H4) hold. Then the chronic infection equilibrium E of system (2) is globally asymptotically stable.

1

0

First, we give the following important remark.

1 C C: A

ea1 s1 a2 s2

Hence, the characteristic equation corresponding to linearized system (5) is given by detðA  nI þ nm1 B þ nm2 C þ nm1 m2 DÞ ¼ 0: Therefore, this characteristic equation can be rewritten as   1  n PðnÞ ¼ 0; 1 þ hd where PðnÞ ¼ ð1 þ haÞð1 þ huÞnm1 þm2 þ1  ð2 þ ha þ huÞnm1 þm2

Proof. We define a discrete Lyapunov functional as follows   Z xn fðx ; y ; v Þ yn a1 s1  ds þ e Wn ¼ xn  x  y /   fðs; y ; v Þ y x v  a n þ ð1 þ uhÞea1 s1 þa2 s2 v /  k v   n1 X fðx jþ1 ; yj ; vj Þvj / þ fðx ; y ; v Þv h fðx ; y ; v Þv j¼nm1   n1 X yjþ1 ; / þ aea1 s1 y h y j¼nm2 where /ðxÞ ¼ x  1  ln x; x 2 Rþ . Clearly, / : Rþ ! Rþ attains its strict global minimum at x ¼ 1 and /ð1Þ ¼ 0. R x  ;y ;v Þ The function w : x # x  x  x fðx fðs;y ;v Þ ds has the global   minimum at x ¼ x and wðx Þ ¼ 0. So, wðxÞ P 0 for all x > 0. Thus, Wn P 0 with equality holding if and only if yn xn vn x ¼ y ¼ v ¼ 1 for all n P 0. The first difference of Wn satisfies DWn ¼ Wnþ1  Wn

k þ nm1 þm2 1  h2 kfð ; 0; 0Þea1 s1 a2 s2 : d 1 is the root of this equation. The remaining Clearly, n ¼ 1þhd roots are given by the solutions of the equation PðnÞ ¼ 0. We have Pð1Þ ¼ h2 auð1  R0 Þ < 0; limn!þ1 PðnÞ ¼ þ1 and P is a continuous function on interval 2 ½1; þ1Þ. Then there exists a n 2 ð1; þ1Þ such that PðnÞ ¼ 0. Therefore, Ef is unstable when R0 > 1. This completes the proof. h

4. Stability of the chronic infection equilibrium In this section, we establish the global stability of the chronic infection equilibrium E , by assuming that R0 > 1 and the function f satisfies the following hypothesis    fðx; y; vÞ fðx; y ; v Þ v 6 0; for all x; y; v > 0:  1 fðx; y ; v Þ fðx; y; vÞ v ðH4Þ

Z xnþ1 fðx ; y ; v Þ ds ¼ xnþ1  xn  fðs; y ;v Þ xn    yn þ ea1 s1 ynþ1  yn þ y ln ynþ1    a vn a1 s1 þa2 s2 þ ð1 þ uhÞe þ fðx ; y ; v Þv h vnþ1  vn þ v ln vnþ1 k      fðxnm1 þ1 ; ynm1 ; vnm1 Þvnm1 fðxnþ1 ; yn ; vn Þvn /  /         fðx ; y ; v Þv fðx ; y ; v Þv      ynm2 þ1 ynþ1 a1 s1  / þ ae y h / y y         fðx ; y ;v Þ yn a 1 s1  x ð  x Þ þ e y  y þ y ln 6 1 nþ1 n nþ1 n ynþ1 fðxnþ1 ; y ; v Þ    a v n þ ð1 þ uhÞea1 s1 þa2 s2 vnþ1  vn þ v ln vnþ1 k þ hfðxnþ1 ; yn ; vn Þvn  hfðxnm1 þ1 ; ynm1 ; vnm1 Þvnm1   fðxnm1 þ1 ; ynm1 ; vnm1 Þvnm1 þ fðx ; y ; v Þv h ln fðxnþ1 ; yn ; vn Þvn    ynm2 þ1 a1 s1 : þ ae h ynþ1  ynm2 þ1 þ y ln ynþ1

Please cite this article in press as: Hattaf, K., Yousfi, N. A numerical method for a delayed viral infection model with general incidence rate. Journal of King Saud University – Science (2015), http://dx.doi.org/10.1016/j.jksus.2015.10.003

A numerical method for a delayed viral infection model with general incidence rate

5

obtain limn!þ1 vn ¼ v and limn!þ1 yn ¼ y . Therefore, we conclude that E is globally asymptotically stable. h

By the inequality ln x 6 x  1, we get fðx ; y ; v Þ Þðxnþ1  xn Þ fðxnþ1 ; y ; v Þ y Þðy  yn Þ þ ea1 s1 ð1  ynþ1 nþ1 a v þ ea1 s1 þa2 s2 ð1  Þðvnþ1  vn Þ k vnþ1    auh a1 s1 þa2 s2 vn e vnþ1  vn þ v ln þ k vnþ1

DWn 6 ð1 

5. Numerical simulations In this section, we will confirm and illustrate our previous theoretical results by numerical simulations. For this, we consider the following delayed discrete model for HIV infection: 8 xnþ1 ¼ xn þ hðk  dxnþ1  bxnþ1 vn Þ; > > <   ynþ1 ¼ yn þ h ea1 s1 bxnm1 þ1 vnm1  aynþ1 ; ð6Þ > >  a2 s2  : vnþ1 ¼ vn þ h ke ynm2 þ1  uvnþ1 :

þ hfðxnþ1 ; yn ; vn Þvn  hfðxnm1 þ1 ; ynm1 ; vnm1 Þvnm1   fðxnm1 þ1 ; ynm1 ; vnm1 Þvnm1 þ fðx ; y ; v Þv h ln fðxnþ1 ; yn ; vn Þvn    ynm2 þ1 a1 s1 : þ ae h ynþ1  ynm2 þ1 þ y ln ynþ1 Using k ¼ dx þ ay ea1 s1 ; fðx ; y ; v Þv ¼ ay ea1 s1 y u a 2 s2 e ¼ v , we obtain k

The model (6) is a special case of (2) with fðx; y; vÞ ¼ bx. The basic reproduction number of (2) is given by and

R0 ¼

kbk a1 s1 a2 s2 e : dau

ð7Þ

    xnþ1  fðx ; y ; v Þ fðx ; y ; v Þ vn fðxnþ1 ; yn ; vn Þ  a1 s1 þ þ hay e DWn 6 hdx 1   1 1 fðxnþ1 ; y ; v Þ fðxnþ1 ; y ; v Þ v fðxnþ1 ; y ; v Þ x     y vnm1 fðxnm1 þ1 ; ynm1 ; vnm1 Þ vn ynm2 þ1 v  a1 s1 þ hay e 1   þ hay ea1 s1 1  ynþ1 v v y vnþ1 fðx ; y ; v Þ       fðxnm1 þ1 ; ynm1 ; vnm1 Þvnm1 ynm2 þ1 xnþ1 fðx ; y ; v Þ 1 þ hay ea1 s1 ln ¼ hdx ea1 s1 1   fðxnþ1 ; yn ; vn Þvnþ1 ynþ1 x fðxnþ1 ; y ; v Þ " # fðx ; y ; v Þ y vn fðxnþ1 ; yn ; vn Þ ynm2 þ1 v fðxnþ1 ; y ; v Þ  a 1 s1    þ hay e 4  fðxnþ1 ; y ; v Þ ynm2 þ1 v fðx ; y ; v Þ fðxnþ1 ; yn ; vn Þ y vnþ1     vn fðxnþ1 ; y ; v Þ vn fðxnþ1 ; yn ; vn Þ þ þ hay ea1 s1 1   þ fðxnþ1 ; yn ; vn Þ v fðxnþ1 ; y ; v Þ v ! y vn fðxnþ1 ; yn ; vn Þ y vnm1 fðxnm1 þ1 ; ynm1 ; vnm1 Þ  a 1 s1  þ hay e ynm2 þ1 v fðx ; y ; v Þ ynþ1 v fðx ; y ; v Þ   fðxnm1 þ1 ; ynm1 ; vnm1 Þvnm1 ynm2 þ1 þ hay ea1 s1 ln fðxnþ1 ; yn ; vn Þvnþ1 ynþ1      xnþ1  fðx ; y ; v Þ vn fðxnþ1 ; y ; v Þ vn fðxnþ1 ; yn ; vn Þ  a1 s1  a1 s1 þ hay þ 1 1  e 1  þ ¼ hdx e fðxnþ1 ; y ; v Þ x v fðxnþ1 ; yn ; vn Þ v fðxnþ1 ; y ; v Þ              ynm2 þ1 v fðx ; y ; v Þ fðxnþ1 ; y ; v Þ y vnm1 fðxnm1 þ1 ; ynm1 ; vnm1 Þ þ / þ/ :  hay ea1 s1 / þ / y vnþ1 ynþ1 v fðx ; y ; v Þ fðxnþ1 ; y ; v Þ fðxnþ1 ; yn ; vn Þ 



Since fðx; y; vÞ is strictly monotonically increasing with respect to x, we obtain that    xnþ1  fðx ; y ; v Þ 1 6 0: 1  fðxnþ1 ; y ; v Þ x According to (H4), we have vn fðxnþ1 ; y ; v Þ vn fðxnþ1 ; yn ; vn Þ 1   þ þ fðxnþ1 ; yn ; vn Þ v fðxnþ1 ; y ; v Þ v    fðxnþ1 ; yn ; vn Þ fðxnþ1 ; y ; v Þ vn ¼ 1  6 0: fðxnþ1 ; y ; v Þ fðxnþ1 ; yn ; vn Þ v Since /ðxÞ P 0 for x > 0, we deduce that DWn 6 0 for all n P 0. Then, fWn g is a monotonically decreasing sequence. Since WðnÞ P 0, it follows that limn!þ1 WðnÞ P 0. Hence, we obtain that limn!þ1 ðWnþ1  Wn Þ ¼ 0, which implies that  y 2 þ1 ¼ yv . From system (2), we limn!þ1 xn ¼ x and limn!þ1 nm vnþ1

In addition, the hypotheses (H1), (H2), (H3) and (H4) are satisfied. Firstly, we simulate the model (6) by using the following parameter values: k = 10 cells mm3 day1 (Perelson et al., 1993), d = 0.02 day1 (Perelson et al., 1993), a = 0.5 day1 (Perelson et al., 1996), u = 3 day1 (Perelson et al., 1996), b = 0.000024 mm3 virion1 day1 (Perelson et al., 1996; Stafford et al., 2000), a1 ¼ d (Perelson et al., 1993), k = 600 virions cell1 day1 (Hattaf and Yousfi, 2012) a2 = 0.65 day1, s1 = 3.5 days, s2 = 2.5 days and h = 0.1 days. By calculating, we have R0 ¼ 0:8813 < 1. By Theorem 3.1, we deduce that the disease-free equilibrium Ef ð500; 0; 0Þ of (6) is globally asymptotically stable, which means that the virus is cleared and the infection dies out. Fig. 1 validates the above analysis.

Please cite this article in press as: Hattaf, K., Yousfi, N. A numerical method for a delayed viral infection model with general incidence rate. Journal of King Saud University – Science (2015), http://dx.doi.org/10.1016/j.jksus.2015.10.003

6

K. Hattaf, N. Yousfi 140 120

400

Infected cells

Uninfected cells

500

300 200

100 80 60 40

100 0

20 0

500

1000

1500

0

2000

0

500

Time

1000

1500

2000

Time Phase diagram

5000 6000 4000

n

3000

v

Virus

4000

2000

2000

Ef

0 200

1000

600 100

0

0

500

1000

1500

2000

0

xn

Plot demonstrates the global stability of Ef .

200

35 30

Infected cells

Uninfected cells

0

yn

Time Figure 1

400 200

150

100

50

25 20 15 10 5

0

0

500

1000

1500

0

2000

0

500

Time

1000

1500

2000

Time Phase diagram

1400 1200

1500

E*

1000

n

800

v

Virus

1000

600

500

400

0 40

200

150 20

0

0

500

1000

1500

2000

Time Figure 2

yn

100 0

50 0

x

n

Plot demonstrates the global stability of E .

Secondly, we choose b = 0.00024 mm3 virion1 day1 (Perelson et al., 1996; Stafford et al., 2000) and the other parameter values are the same as above. The reason to just modify the parameter b is based on the fact that R0 is an increasing function with respect to b (see the explicit formula (7) for R0 ). By calculating, we have R0 ¼ 8:8128 > 1. Then,

system (6) has a unique chronic infection equilibrium E ð56:7359; 16:5319; 651:0636Þ. By applying Theorems 3.1 and 4.2, we see that Ef becomes unstable and E is globally asymptotically stable. In this case, the virus persists in the host and the infection becomes chronic. Fig. 2 confirms this observation.

Please cite this article in press as: Hattaf, K., Yousfi, N. A numerical method for a delayed viral infection model with general incidence rate. Journal of King Saud University – Science (2015), http://dx.doi.org/10.1016/j.jksus.2015.10.003

A numerical method for a delayed viral infection model with general incidence rate

7

References

Figure 3 Plot of the basic reproduction number R0 as a function of the time delays s1 and s2 .

According to the above, we deduce a strategy to control the viral infection. This strategy is based on reducing the value of R0 and making it less than or equal to one. From the explicit expression of R0 in (7), it is clear that with the increase in time delays s1 and s2 , the value of R0 decreases which is demonstrated in Fig. 3. 6. Conclusion In this work, we have proposed a discrete mathematical model with two delays to describe the dynamics of viral infection, such as human immunodeficiency virus (HIV), the hepatitis B virus (HBV) and the hepatitis C virus (HCV). The discrete model is derived from the continuous system (1) by using a mixed Euler method. Also, the infection transmission process is modeled by a general incidence function that includes various types of incidence rate existing in the literature. We have proved that the proposed mixed Euler method is unconditionally positive. Furthermore, the dynamical behaviors of the the delayed discrete model are investigated by linearization method and by constructing suitable discrete Lyapunov functionals. More precisely, we have proved that the disease-free equilibrium Ef is globally asymptotically stable if the basic reproduction number satisfies R0 6 1, which means that the virus is cleared and the infection dies out. When R0 > 1; Ef becomes unstable and the chronic infection equilibrium E is globally asymptotically stable. In this case, the virus persists in the host and the infection becomes chronic. Therefore, we conclude that the discrete model has the same qualitative properties as the corresponding continuous viral infection model (1) with no restriction on the time step size. Acknowledgement We would like to thank the editor and anonymous referees for their very helpful comments and suggestions that greatly improve the presentation of this work.

Eikenberry, S., Hews, S., Nagy, J.D., Kuang, Y., 2009. The dynamics of a delay model of HBV infection with logistic hepatocyte growth. Math. Biosci. Eng. 6, 283–299. Gourley, S.A., Kuang, Y., Nagy, J.D., 2008. Dynamics of a delay differential model of hepatitis B virus. J. Biol. Dyn. 2, 140–153. Hattaf, K., Yousfi, N., 2011. A delay differential equation model of HIV with therapy and cure rate. Int. J. Nonlinear Sci. 12, 503–512. Hattaf, K., Yousfi, N., 2012. Two optimal treatments of HIV infection model. World J. Model. Simul. 8 (1), 27–35. Hattaf, K., Yousfi, N., 2014. Global stability of a virus dynamics model with cure rate and absorption. J. Egypt. Math. Soc. 22 (2014), 386–389. Hattaf, K., Yousfi, N., Tridane, A., 2012. Mathematical analysis of a virus dynamics model with general incidence rate and cure rate. Nonlinear Anal. RWA 13, 1866–1872. Hattaf, K., Yousfi, N., Tridane, A., 2013. Stability analysis of a virus dynamics model with general incidence rate and two delays. Appl. Math. Comput. 221, 514–521. Hattaf, K., Yousfi, N., Tridane, A., 2014. A Delay virus dynamics model with general incidence rate. Differ. Equ. Dyn. Syst. 22 (2), 181–190. Huang, G., Ma, W., Takeuchi, Y., 2011. Global analysis for delay virus dynamics model with Beddington–DeAngelis functional response. Appl. Math. Lett. 24 (7), 1199–1203. Li, D., Ma, W., 2007. Asymptotic properties of an HIV-1 infection model with time delay. J. Math. Anal. Appl. 335, 683–691. Li, M.Y., Shu, H., 2010. Global dynamics of an in-host viral model with intracellular delay. Bull. Math. Biol. 72, 1492–1505. Perelson, A.S., Kirschner, D.E., De Boer, R., 1993. Dynamics of HIV infection of CD4+ T cells. Math. Biosci. 114 (1), 81–125. Perelson, A.S., Neumann, A., Markowitz, M., Leonard, J., Ho, D.D., 1996. HIV-1 dynamics in vivo: virion clearance rate, infected cell life-span, and viral generation time. Science 271, 1582–1586. Stafford, M.A., Corey, L., Cao, Y., Daar, E.S., Ho, D.D., Perelson, A. S., 2000. Modeling plasma virus concentration during primary HIV infection. J. Theor. Biol. 203 (3), 285–301. Sun, Q., Min, L., 2014. Dynamics analysis and simulation of a modified HIV infection model with a saturated infection rate. Computat. Math. Methods Med. 2014. http://dx.doi.org/10.1155/ 2014/145162. Tian, X., Xu, R., 2010. Asymptotic properties of a hepatitis B virus infection model with time delay. Disc. Dyn. Nat. Soc. http://dx.doi. org/10.1155/2010/182340, 21 pages. Vargas-De-Leo´n, C., 2012. Stability analysis of a model for HBV infection with cure of infected cells and intracellular delay. Appl. Math. Computat. 219, 389–398. Wang, T., Hu, Z., Liao, F., Ma, W., 2013. Global stability analysis for delayed virus infection model with general incidence rate and humoral immunity. Math. Comput. Simul. 89, 13–22. Xiang, H., Feng, L.-X., Huo, H.-F., 2013. Stability of the virus dynamics model with Beddington–DeAngelis functional response and delays. Appl. Math. Model. 37, 5414–5423. Xu, R., 2011. Global stability of an HIV-1 infection model with saturation infection and intracellular delay. J. Math. Anal. Appl. 375, 75–81. Zhou, X., Cui, J., 2011. Global stability of the viral dynamics with Crowley–Martin functional response. Bull. Korean Math. Soc. 48, 555–574. Zhu, H., Zou, X., 2008. Impact of delays in cell infection and virus production on HIV-1 dynamics. Math. Med. Biol. 25, 99–112. Zhuo, X., 2012. Analysis of a HBV infection model with non-cytolytic cure process. IEEE 6th Int. Conf. Syst. Biol., 148–151

Please cite this article in press as: Hattaf, K., Yousfi, N. A numerical method for a delayed viral infection model with general incidence rate. Journal of King Saud University – Science (2015), http://dx.doi.org/10.1016/j.jksus.2015.10.003