Analysis of a stochastic distributed delay epidemic model with relapse and Gamma distribution kernel

Chaos, Solitons and Fractals 133 (2020) 109643

Contents lists available at ScienceDirect

Chaos, Solitons and Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos

Analysis of a stochastic distributed delay epidemic model with relapse and Gamma distribution kernel Tomás Caraballo a,∗, Mohamed El Fatini b, Mohamed El Khalifi b, Richard Gerlach c, Roger Pettersson d a

Dpto. Ecuaciones Diferenciales y Análisis Numérico, Facultad de Matemáticas, Universidad de Sevilla, Tarfia s/n, 41012, Sevilla, Spain Dept. of Mathematics, Faculty of Sciences, Ibn Tofail University, Kénitra BP 133, Morocco Discipline of Business Analytics, Business School, The University of Sydney, Australia d Linnaeus University, Department of Mathematics, Växjö 351 95, Sweden b c

a r t i c l e

i n f o

Article history: Received 15 December 2019 Revised 19 January 2020 Accepted 21 January 2020

Keywords: Epidemic model Distributed delay Extinction Persistence Stationary distribution

a b s t r a c t In this work, we investigate a stochastic epidemic model with relapse and distributed delay. First, we prove that our model possesses and unique global positive solution. Next, by means of the Lyapunov method, we determine some sufficient criteria for the extinction of the disease and its persistence. In addition, we establish the existence of a unique stationary distribution to our model. Finally, we provide some numerical simulations for the stochastic model to assist and show the applicability and efficiency of our results.

1. Introduction Statistics reported by the World Health Organisation shows the threats of infectious diseases to global health. It has estimated that tuberculosis caused 1.18 millions of deaths and 2.56 millions of peoples died from pneumonia in 2017. At the end of 2018, a total of 32 millions of peoples were killed by HIV since the beginning of the epidemic. Biological literature has been enriched by contribution of mathematicians to understand, predict behaviors and control the spread of detrimental epidemics. Since the basic SIR model introduced in works of Kermack and McKendrick [12– 14], many important extensions have been developed including different characteristics of the epidemics [1–4,7,8,11,17,19,23,25,27]. Since many species should reveal time delay, the introduction of time delays into compartmental epidemic models presents a fascinating improvement as it models sojourn times in a specified state, e.g. the infective state, see [5,20,22,26]. In [18], authors studied the following stochastic SIR epidemic model with distributed time delay

∗

Corresponding author. E-mail address: [email protected] (T. Caraballo).

https://doi.org/10.1016/j.chaos.2020.109643 0960-0779/© 2020 Elsevier Ltd. All rights reserved.

© 2020 Elsevier Ltd. All rights reserved.

⎧   t D(t − s )I (s )ds dt + σ SdB(t ), ⎨ dS(t ) = μ − μS(t ) − β S(t ) −∞  t dI (t ) = β S(t ) −∞ D(t − s )I (s )ds − (μ + λ + δ )I (t ) dt ⎩ dR(t ) = (λI (t ) − μR(t ) )dt, (1) where S(t), I(t) and R(t) represent the densities of susceptible, infected and recovered individuals at time t respectively. The parameter μ is the recruitment rate, β the infection transmission rate, λ the recovery rate and δ the death rate revealed by the t disease. The quantity −∞ D(t − s )I (s )ds is the force of infection at time t. The function D: [0, ∞) → [0, ∞), called delay kernel, ∞ is an L1 −function verifying 0 D(s )ds = 1. The standard Brownian motion B(t) with intensity σ is used to model the environmental noise on susceptible class. Connected with the above statements, we suggest to study the dynamics of a stochastic SIR epidemic model incorporating relapse, possessing distributed time delay with environmental noise proportional to S(t), I(t) and R(t) as

  ⎧ t dS(t ) = μ − μS(t ) − β S(t ) −∞ D(t − s )I (s )ds dt + σ1 S(t )dB1 (t ), ⎪   ⎨ t dI (t ) = β S(t ) −∞ D(t − s )I (s )ds − (μ + λ + δ )I (t ) + γ R(t ) dt ⎪ +σ2 I (t )dB2 (t ) ⎩ dR(t ) = (λI (t ) − (μ + γ )R(t ) )dt + σ3 R(t )dB3 (t ), (2)

2

T. Caraballo, M.E. Fatini and M.E. Khalifi et al. / Chaos, Solitons and Fractals 133 (2020) 109643

where γ is the rate of reinfection, and Bi (t ), i = 1, 2, 3 are independent Brownian motions with intensities σi , i = 1, 2, 3. In practice, it is convenient to use the kernel delay D with Gamma distribution [21], that is

D (s ) =

sn α n+1 e−α s , n!

s > 0,

where the positive number α is the rate of decay of effect of past memories. In this paper, we consider the weak kernel (i.e. n = 1) and for simplicity we let

Z (t ) =



t

1. There exists a positive number κ such that

i, j=1 ai j ξi ξ j ≥ κ|ξ |2 , ξ ∈ D, where D ⊂ Rn is a bounded open set with compact

closure. 2. There exists a nonnegative C 2 −function V : Dc → R such that LV is negative for any x ∈ Dc = Rn \ D. Then, the Markov process X(t) has a unique ergodic stationary distribution μ( · ) with density in Rn such that



P lim

T →∞

α e−α (t−s) I (s )ds,

n

1 T



T

0

Rn



h(X (t ))dt =

Rn

h(x )μ(dx ) = 1,

then it follows that

for any x ∈ measure μ.

dZ (t ) = α (I (t ) − Z (t ) )dt.

2. Existence and uniqueness of global positive solution

We also allow the component Z(t) to be affected by the environmental variability as it depends of the perturbed density of infected individuals. Hence, system (2) can be rewritten as

In the sequel, we show that system (3) admits a unique global nonnegative solution for any initial condition on R4+ .

−∞

⎧ dS(t ) = (μ − μS(t ) − β S(t )Z (t ) )dt + σ1 S(t )dB1 (t ), ⎪ ⎪ ⎨ dI (t ) = (β S(t )Z (t ) − (μ + λ + δ )I (t ) + γ R(t ) )dt +σ2 I (t )dB2 (t ), ⎪ ⎪ ⎩dR(t ) = (λI (t ) − (μ + γ )R(t ) )dt + σ3 R(t )dB3 (t ), dZ (t ) = α (I (t ) − Z (t ) )dt + σ4 Z (t )dB4 (t ),

(3)

where B4 (t) stands for Brownian motion with intensity σ 4 . Despite the realism presented by epidemic models with distributed delay, their stochastic analysis is still quite difficult and limited. To the best of our knowledge, only sufficient conditions are obtained for the extinction and persistence of the diseases with such kind of delay (see, for instance, (1)). Actually, more general results are proved for models with maturation and bounded delays [9,24]. The content of this paper is structured as follows. Section 2 is devoted to establish the existence of a unique global positive solution for the stochastic epidemic model (3). In Section 3, by means of Lyapunov functions, we establish a sufficient criterion for the exponential extinction of the disease. The persistence in mean result under some appropriate conditions is proved in Section 4. In Section 5, we analyze the existence of an ergodic stationary distribution to system (3). Finally, in Section 6, some numerical examples are presented to provide a practical exhibition and explanation to our analytical results.  Throughout this paper, let , F, {Ft }t≥0 , P be a complete probability space with a filtration {Ft }t≥0 satisfying usual conditions. Define the space Rn+ as follows

Rn+ = {x ∈ Rn : xi > 0, i = 1, · · · , n}.

dXt = f (X )dt +

gr (X )dBr (t ),

X (0 ) = X0 ,

(4)

r=1

with the diffusion matrix defined by





A ( x ) = ai j ( x ) ,

ai j ( x ) =

k

Theorem 2.1. For any given initial value (S(0 ), I (0 ), R(0 ), Z (0 ) ) ∈ R4+ , there exists a unique local solution (S(t), I(t), R(t), Z(t)) to system (3) on t ≥ 0 and the solution remains in R4+ almost surely. Proof. Let (S(0 ), I (0 ), R(0 ), Z (0 ) ) ∈ R4+ . Since the considered system has a locally Lipschitz coefficients of linear growth, there exists a unique solution (S(t), I(t), R(t)) on t ∈ [0, τ e ), where τ e is the explosion time. Now, we shall prove that τe = ∞ a.s. To this end, we need to construct a C 2 −function V : R4+ → R+ ∪ {0} such that

lim inf

k→∞,(S,I,R,Z )∈R4+ \Dk

V (S, I, R, Z ) = ∞



and LV (S, I, R, Z ) ≤ C,

4

where Dk = 1k , k for k ≥ k0 ≥ 1, k0 is a sufficiently large integer such that (S(0 ), I (0 ), R(0 ), Z (0 ) ) ∈ Dk0 and C is a positive constant. For any k ≥ k0 , we define τk = inf{t ∈ [0, τe ) : (S, I, R, Z ) ∈ / D k }. Obviously, (τk )k≥k is an increasing sequence with limit τ ∞ veri0 fying τ ∞ ≤ τ e a.s. It is sufficient to prove that τ∞ = ∞ a.s. If this is false, then there is a pair of constants T > 0 and ε ∈ (0, 1) such that P{τ∞ ≤ T } > ε . Hence, there is an integer k1 ≥ k0 such that P{τk ≤ T } ≥ ε for all k ≥ k1 . Let a and b be two positive constants and define



V (S, I, R, Z ) = S − a − a ln



+b Z − b − b ln

 Z

S a



+ (I − 1 − ln I ) + (R − 1 − ln R )

.

b

Since − ln x −−−→ ∞ and −x − m − m ln x→0+

Let Xt be a regular homogeneous Markov process in Rn described by the SDE k

and h( · ) is an integrable function with respect to the

x −−→ m − x→∞

∞ for any m > 0,

one can easily obtain that

lim inf

k→∞,(S,I,R,Z )∈R4+ \Dk

V (S, I, R, Z ) = ∞.

Now, using Itô’s formula on V, we obtain that



   1 σ1 SdB1 + 1 − σ2 IdB2 I     1 b2 σ3 RdB3 + b − + 1− σ4 ZdB4 ,

dV (S, I, R, Z ) =LV (S, I, R, Z )dt + 1 −

gir (x )grj (x ).

a S

R

Z

r=1

We introduce the differential operator L associated to (4) and acts on any twice continuously differentiable function V as follows

LV (x ) =

n

i=1

f i (x )

n ∂ V (x ) 1

∂ 2V (x ) + . ∂ xi 2 ∂ xi ∂ x j i, j=1

The following lemma is a worthwhile criterion for the existence of ergodic stationary distributions [15]. Lemma 1. Assume the following assumptions hold,

where







a 1 (μ − μS − β SZ ) + 1 − S I (β SZ − (μ + λ + δ )I + γ R )

LV (S, I, R, Z ) = 1 −





+ 1−

1 (λI − (μ + γ )R ) R

+

b2 Z



b−



(α I − α Z ) +

aσ12 2



T. Caraballo, M.E. Fatini and M.E. Khalifi et al. / Chaos, Solitons and Fractals 133 (2020) 109643

σ22

+

+

2

σ32 2

+

bσ42 2

p(x ) = Qx aσ12 2

=3μ + aμ + λ + δ + γ +

σ32

σ22

+

I

γR



Let q(x ) = exp −2

 q(x ) = exp −2

LV (S, I, R, Z ) ≤ C,

(5)

aσ12 σ 2 σ 2 bσ42 where C = 3μ + aμ + λ + δ + γ + + 2 + 3 + . 2 2 2 2 The remaining part of the proof is similar to Theorem 1.2 of [6]. However, for the convenience of the reader, we will include it here. By (5), we have

   1 σ1 SdB1 + 1 − σ2 IdB2 I     2 1 b σ3 RdB3 + b − + 1− σ4 ZdB4 . a S

R



∞ 1

∞ 0



σ1

∧

k



k



− a + a ln(ka ) ∧ (k − 1 − ln k )

− 1 + ln k ∧ b k − b − b ln

k b



1

∧b

k



2μ

σ12

∞



x−2 q(x )−1 dx < ∞.

,

In this section, we establish sufficient conditions for the extinction of the disease in our model (3). We set the parameter

∞

x

− 2μ2 −1 − 22μ σ1

e

σ1 x

dx

0

− 2μ2   σ1 2μ  2 σ1 σ12

2μ



∞

x

− 2μ2 σ1

e

− 22μ

σ1 x

dx

0



− 2μ2 +1   σ1 2μ =Q  −1 σ12 σ12 2μ = . 2μ − σ12

− b + b ln(kb) .

2μ

Then,



∞ 0

(x − 1 )2 p(x )dx =



∞

0

(x2 − 2x + 1 ) p(x )dx =

σ12 . 2μ − σ12

Next, we consider the nonnegative C 2 −function V1 (I, R, Z ) = α1 I + α2 R + α3 Z, where α1 = μ+ωλ1+δ , α2 = μω+2γ , and α3 = ωα3 . ωi , i = 1, · · · , 3 are positive constants to be determined. Using Itô’s formula for ln V1 , we have

β λγ + . μ + λ + δ (μ + λ + δ )(μ + γ )

Then, we have the following theorem.

d ln V1 =

Theorem 3.1. Let (S(t), I(t), R(t), Z(t)) be the solution to system (3) for an initial value (S(0), I(0), R(0), Z(0)) in R4+ . Assume that

−

σ12

Then, the disease dies out exponentially with probability one. Moreover, the distribution of S(t) converges weakly to the measure which has the density

1 {α1 (β SZ − (μ + λ + δ )I + γ R ) + α2 (λI − (μ + γ )R ) V1 +α3 α (I − Z )}dt

≤ 2μ and √ βασ1 R0 min{μ + λ + δ, μ + γ , α}( R0 − 1 ) + < 0.  (μ + λ + δ ) 2μ − σ12 R0 < 1,

0

1

σ12

1

x2 p(x )dx = Q



3. Extinction of the disease



∞

and

 2μ2 +1   −1 σ1  2σμ2 + 1 . By the comparison theo-



0

Letting k → ∞ leads to the contradiction ∞> V (S(0 ), I (0 ), R(0 ), Z (0 )) + CT = ∞. Hence τ∞ = ∞ a.s. This concludes the proof. 

R0 =

(6)

= 1,





q(x )dx = ∞

σ1 x

and

a

1 0

=Q



≥ εθk ,

1

X ( 0 ) = S ( 0 ) > 0.

μ − μu du . We have (σ1 u )2

x

1



e

xp(x )dx = Q

≥ E 1{ τk ≤ T }V (S(τk ∧ T ), I (τk ∧ T ), R(τk ∧ T ), Z (τk ∧ T ) )

∧



rem, we have S(t) ≤ X(t), t ≥ 0 a.s. By direct computations, one can obtain



V (S(0 ), I (0 ), R(0 ), Z (0 )) + CT



− 2μ2 −2 − 22μ

which gives that

θk = k − a − a ln

1

Therefore, by Theorem 1.16 in [16], (6) has the ergodic property and the density of its invariant law is given by

≤ V (S(0 ), I (0 ), R(0 ), Z (0 )) + CT ,

  k 1

 2μ2 +1   −1 σ1  2σμ2 + 1 . σ2

2μ

 − 2μ2 2μ2 22μ 1 μ 1 − ln x − = e σ1 x σ1 e σ1 x . 2 x σ1

q(x )dx = ∞,

where Q =

Z

EV (S(τk ∧ T ), I (τk ∧ T ), R(τk ∧ T ), Z (τk ∧ T ) )





Hence,

p(x ) = Qx

Integrating both sides of the last inequality over [0, τ k ∧T] yields

where

x > 0,

dXt = (μ − μX )dt + σ1 X dB1 (t ),

Choosing a = μβ+δ and b = μα+δ leads to



,

Proof. First, we consider Xt solution to the following SDE

b2 α I + + + . I R Z

dV (S, I, R, Z ) ≤Cdt + 1 −

σ1 x

e

1



λI

σ1

where the normalisation constant Q =

bσ42 + + + + (bα − (μ + δ ) )I 2 2 2    a + (aβ − bα )Z − μ S + R + S

β SZ

− 2μ2 −2 − 22μ

3

 1  (α1 σ2 I )2 + (α2 σ3 R )2 + (α3 σ4 Z )2 dt 2 2V1

1 (σ2 IdB2 + σ3 RdB3 + σ4 ZdB4 ) V1 1 =L ln V1 dt + (σ2 IdB2 + σ3 RdB3 + σ4 ZdB4 ), V +

where

4

T. Caraballo, M.E. Fatini and M.E. Khalifi et al. / Chaos, Solitons and Fractals 133 (2020) 109643

L ln V1 =

1 {α1 (β Z − (μ + λ + δ )I + γ R ) + α2 (λI − (μ + γ )R ) V1 +α3 α (I − Z )} +

≤

α1 β V1

+

 1  (α1 σ2 I )2 + (α2 σ3 R )2 + (α3 σ4 Z )2 2V 2

( S − 1 )Z −

<0

 ω α1 β 1 1 |X − 1| + ( β Z − ( μ + λ + δ )I + γ R ) α3 V1 μ + λ + δ

ω2 + (λI − (μ + γ )R ) + ω3 α (I − Z ) . μ+γ

Let M0 be the matrix defined by

0

λ M 0 = ⎝ μ+ γ

γ μ+λ+δ

β μ+λ+δ

0 0

0 0

1

⎞



√ λγ R0 = μ+βλ+δ + (μ+λ+δ )(μ+γ ) . Choosing the √  λ , 1 leads to vector (ω1 , ω2 , ω3 ) = R 0 , μ+ γ

(ω1 , ω2 , ω3 )M0 =

R0 (ω1 , ω2 , ω3 ).

  α1 β 1 L ln V1 ≤ |X − 1| + (ω1 ω2 ω3 ) M0 (I, R, Z )T − (I, R, Z )T α3 V1   α1 β 1  = |X − 1| + R0 − 1 ( ω1 I + ω2 R + ω3 Z ) α3 V1   α1 β ≤ |X − 1| + min{μ + λ + δ, μ + γ , α} R0 − 1 . α3

ln I (t ) < 0, t

lim sup t→∞

R0s =

μ+λ+δ+



t

σ3 R



d B3 +

lim inf



lim

dV˜2 = −

0

|X (s ) − 1|ds =

0

|x − 1| p(x )dx ≤



=  Consequently

2

βK σ2 + β K 2α4 −

λγ μ+ γ +

, σ32 2

. We then have the following result.

2





t

I (s )ds ≥ Λ 1 −

0

1 R0s



a.s.,

(7)

μα  . In other words, the epidemic will be σ42 β α+ 

∞ 0

(x − 1 )2 p(x )dx

σ1

2μ − σ12

.

β1

(μ − μS − β SZ ) + β1

−

 12

β2 R

− β3

(λI − (μ + γ )R ) + β2

α

σ12

σ32 2

+ β4 α (I − Z ) 2 − (β1 σ1 dB1 + σ2 dB2 + β2 σ3 dB3 + σ4 (β3 − β4 Z )dB4 ) Z

(I − Z ) + β3

σ42

=LV˜2 dt − (β1 σ1 dB1 + σ2 dB2 + β2 σ3 dB3 ), where

(8)

constants to be specified later. goes to ∞ as (S, I, R, Z) goes to has a lower bound M. Then we V˜2 = V2 (S, I, R, Z ) − M. Applying

S 2 σ2 1 − (β SZ − (μ + λ + δ )I + γ R ) + 2 I 2

0

1 t

1 t

t

then

t→∞

σ12

σ22

where βi , i = 1, · · · , 4 are positive The function V2 is continuous and the boundary of R4+ . Thus it must define the nonnegative function Itô’s formula to V˜2 , we obtain

the strong law of large numbers for local martingales. ∞ As the solution to system (6) is ergodic and xp(x )dx < ∞,

∞

a.s.

V2 (S, I, R, Z ) = −β1 ln S − ln I − β2 ln R − β3 ln Z + β4 Z,



t→∞



ln Z (t ) <0 t

Proof. Let (S(0 ), I (0 ), R(0 ), Z (0 )) ∈ R4+ and consider the function

σ4 Z dB4 is a local marV1 0 V1 0 V1 tingale with finite quadratic variation. Hence, lim Ntt = 0 a.s. by 0

d B2 +

t→∞

2

+ min{μ + λ + δ, μ + γ , α} R0 − 1 t α1 β 1 N + |X (s ) − 1|ds + t , α3 t 0 t

σ2 I

lim sup

Theorem 4.1. If R0s > 1, then for any given initial value (S(0 ), I (0 ), R(0 ), Z (0 )) ∈ R4+ , the corresponding solution to (3) verifies

permanent.

ln V1 (I (t ), R(t ), Z (t )) ln V1 (I (0 ), R(0 ), Z (0 )) ≤ t t

t

μ+

1 + (σ2 IdB2 + σ3 RdB3 + +σ4 ZdB4 ). V Integrating both sides and dividing by t gives



μ

where K =

where Λ =

t

a.s.,

The study of the persistence of diseases is an important approach to know more on epidemics behaviors, since it provides the conditions under which the diseases are prevalent. For reasons of simplification, we define the quantity

  α1 β |X − 1| + min{μ + λ + δ, μ + γ , α} R0 − 1 α3





This implies that I(t) goes to 0 exponentially with probability 1. The proof is then complete. 

t→∞

Hence, we obtain the inequality

where Nt =

R0 − 1

which is the desired result. Moreover,

Therefore,

d ln V1 ≤



4. Persistence in mean of the disease

⎠,

with eigenvalue



α1 β σ  1 α3 2 μ − σ 2 1

= min{μ + λ + δ, μ + γ , α} √ βασ1 R0 +  (μ + λ + δ ) 2μ − σ12

1 {α1 (β Z − (μ + λ + δ )I + γ R ) V1

+α2 (λI − (μ + γ )R ) + α3 α (I − Z )}

⎛

t→∞

  ln V1 (I (t ), R(t ), Z (t )) ≤ min{μ + λ + δ, μ + γ , α} R0 − 1 t

α1 β |X − 1| α3 +

≤

lim sup

T. Caraballo, M.E. Fatini and M.E. Khalifi et al. / Chaos, Solitons and Fractals 133 (2020) 109643

5

Fig. 1. Trajectories of S(t), I(t) and R(t) components of (S(t), I(t), R(t), Z(t)) the solution of system (3) with parameters given in Table 1.

LV˜2 = −

β SZ I

−

β1 μ

−

S

β3 α I

−

Z

γR I

−

β2 λI

We further obtain



R

  σ2 + (β1 β − β4 α )Z + β1 μ + 1 2   2 σ2 σ + μ + λ + δ + 2 + β2 μ + γ + 3 2 2   2 σ + β3 α + 4 + β4 α I,

dV˜2 ≤

−

I γR − − I

β1 μ S

β2 λI R

we obtain



LV˜2 ≤ − 3 3

−

≤ −3 3

Z

≤ −2



βμαβ1 β3 and

λγ β2 ,

t 1 Mt = − (β1 σ1 B1 (t ) + σ2 B2 (t ) + β2 σ3 B3 (t ) + σ4 (β3 − t 0 β4 Z )dB4 ) is a continuous local martingale. Using the large numbers theorem for martingales, we obtain lim Mt = 0. Therefore where

σ22

lim inf t→∞

2

    σ2 σ2 μ + γ + 3 + β3 α + 4 + β4 α I. 2

LV˜2 ≤ −β K + μ + λ + δ +



1 t



t 0

I (s )ds ≥

K



β1

This completes the proof.

1−

1 R0s



t→∞

a.s.



2

  σ2 2 For β2 = λγ / μ + γ + 23 ,   σ2 β3 = β K/ α + 24 , we get

1 = −β K 1 − 0 Rs

+ β1 β I dt

V˜2 (S(t ), I (t ), R(t ), Z (t ) ) t   1 V˜2 (S(0 ), I (0 ), R(0 ), Z (0 ) ) ≤ − βK 1 − 0 t Rs 1 t + β1 β I (s )ds + Mt , t 0

2





0≤

   σ2 βμαβ1 β3 − 2 λγ β2 + β1 μ + 1

+μ+λ+δ+ + β2



β3 α I



Integrating both sides of the last inequality and dividing by t, leads to

Choosing β4 = β1 β /α and using the facts that

β SZ

1 −β K 1 − 0 Rs

−(β1 σ1 dB1 + σ2 dB2 + β2 σ3 dB3 + σ4 (β3 − β4 Z )dB4 ).

2

−



σ22 2

β1 =

+ βK

+ β1 β I.

βK2

 α+

σ42 2

σ42 λγ − 2α μ+γ +

5. Stationary distribution and positive recurrence



/(μα )

σ32 2

and

+ β1 β I

In the previous sections, we have investigated both the extinction and the permanence of the disease under some conditions. In this section, we will show that the solution to system (3) possesses a unique stationary distribution by means of the method given in [15] and used in [28]. Theorem 5.1. For any (S(0 ), I (0 ), R(0 ), Z (0 ) ) ∈ R4+ , the corresponding solution to system (3) is positive recurrent and has a unique ergodic stationary distribution if the condition R0s > 1 holds.

6

T. Caraballo, M.E. Fatini and M.E. Khalifi et al. / Chaos, Solitons and Fractals 133 (2020) 109643

Fig. 2. Simulation of paths S(t), I(t) and R(t) where (S(t), I(t), R(t), Z(t)) is the solution to system (3) using data of Table 2.



Table 1 Parameters values used in Fig. 1.

μ

β

λ

δ

γ

α

σ1

σ2

σ3

σ4

0.12

0.27

0.48

0.1

0.113

0.11

0.2

0.23

0.15

0.07

m   μ+δ μ+δ μ+δ Z I − μR − Z LV4 = S + I + R + μ − μS − 2α 2 2  m−1  2 2  m μ+δ + S+I+R+ Z σ1 S + σ22 I2 + σ32 R2 + σ42 Z 2 2 2α + cλI − c (μ + γ )R.

Proof. It is easy to verify the first condition of Lemma 1. We only have to prove that the second condition is also verified. Let us consider the C 2 -function V6 defined by

We choose m sufficiently small so that the following inequality holds,

V5 (S, I, R, Z ) =MV2 + V3 + V4 ,

m

where the

function

V2

is given in the



previous section,

m+1

V3 = − ln S − α1 ln Z and V4 = m1+1 S + I + R + μ2+αδ Z + cR. M, m and c are positive constants satisfying some conditions as we shall see later. Since V6 is continuous and tends to ∞ as (S, I, R, Z) tends to the boundary of R4+ , it admits a lower bound M∗ and achieves this lower bound at a point in R4+ . Therefore, we construct the nonnegative function V˜5 by

V˜5 (S, I, R, Z ) = V5 (S, I, R, Z ) − M∗ . By simple computations, we can get that

LV2 ≤ − + ββ1 I, where

=

K



1 1− 0 Rs

β1

 2  μ + δ ( μ + δ )2 , σ1 ∨ σ22 ∨ σ32 ∨ σ42 < min{μ, }. 2 4α

Then, we have

LV4 ≤ −μmSm+1 − −



1 2α m

where

ϒ=

μ+δ

1



2α m

> 0.

Moreover,

σ2 σ2 I LV3 = − − + β Z + 1 + μ + 1 + 4 . S Z 2 2α μ

Similarly, we can deduce that

Rm+1 − c (μ + γ )R

Z m+1 + ϒ ,

μ+δ 4

Im+1 −

μ 2

Rm+1

m+1 Z m+1 + cλI

2



 m 2 μ+δ σ ∨ σ22 ∨ σ32 ∨ σ42 S + I + R + Z 2 1 2α 



2

m+1

−μ(1 − m )Sm+1 −

μ+δ

μ

Im+1 −



(S,I,R,Z )∈R4+

+

4

2

sup

−

μ+δ

+μ S + I + R +

μ+δ Z 2α

m−1



S 2 + I 2 + R2

m

.

Hence,

LV˜5 ≤ − M + ββ1 I − c (μ + γ )R − −

μ+δ 4

Im+1 −

μ 2

Rm+1 −

1 4α m

μ S



−

I − μmSm+1 Z

μ+δ 2

m+1

Z m+1 +  ,



T. Caraballo, M.E. Fatini and M.E. Khalifi et al. / Chaos, Solitons and Fractals 133 (2020) 109643

7

Fig. 3. The kernel density functions of susceptible, infected and removed compartments of (3) at t = 10 0 0 (red), t = 20 0 0 (blue) and t = 30 0 0 (green).(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

where

  = sup − Z>0

1 4α m



μ+δ

m+1 Z m+1 + β Z + ϒ + 1 + μ +

2

σ12 2



+

σ42 . 2α

D=

1

(S, I, R, Z ) ∈ R4+ ; ε1 ≤ S ≤ and

ε3 ≤ Z ≤

ε1

 1 ε3

, ε2 ≤ I ≤

1

ε2

, ε22 ≤ S ≤

1

ε

2 2

,



1



,

ε1  D3 = (S, I, R, Z ) ∈ R4+ ; Z < ε3 ,   1 D4 = (S, I, R, Z ) ∈ R4+ ; Z > , ε3   D5 = (S, I, R, Z ) ∈ R4+ ; S ≥ ε1 , I < ε2 , Z ≥ ε3 ,

1 4 D6 = (S, I, R, Z ) ∈ R+ ; R > 2 , ε2   D7 = (S, I, R, Z ) ∈ R4+ ; I > ε2 , R < ε22 ,   1 D8 = (S, I, R, Z ) ∈ R4+ ; I > . ε2 

μ+δ 4

μ + 1 ε1

Im+1 + 



Z>0

Case 2. Let (S, I, R, Z) ∈ D2 . Since −μmSm+1 < − μmm +1 , we choose ε 1 sufficiently small such that

LV˜5 ≤ − M − μmSm+1 − mu

D1 = (S, I, R, Z ) ∈ R4+ ; S < ε1 , D2 = (S, I, R, Z ) ∈ R4+ ; S >

−



We divide the domaine R4+ \ D into eight domains



S

where 1 = sup ββ1 I − μ4+δ Im+1 +  .

.



≤ − M −

μ

<0,

Let εi , i = 1, · · · , 3 be positive reals and define the closed set



LV˜5 ≤ − M + ββ1 I −

Let us show that LV˜5 is negative at any vector in R4+ \ D. Case 1. If (S, I, R, Z) ∈ D1 , we use the fact that − μS < − εμ and 1 choose ε 1 sufficiently small such that

≤ − M −

μm + 1 ε1m+1

ε1

m m+1 S + 1 2

<0. Case 3. If (S, I, R, Z) ∈ D3 , then, there exists a sufficiently large M such that

LV˜5 ≤ −M + 1 < 0. Case 4. If (S, I, R, Z) ∈ D4 . With the fact that −Z m+1 < − choose ε 3 sufficiently small such that

LV˜5 ≤ − M − ≤ − M −

1 4α m 1 4α m

 

μ+δ 2

μ+δ 2ε3

1

ε3m+1

, we

m+1 Z m+1 + 1

m+1 + 1

<0. Case 5. If (S, I, R, Z) ∈ D5 , there is a sufficiently large number M such that

8

T. Caraballo, M.E. Fatini and M.E. Khalifi et al. / Chaos, Solitons and Fractals 133 (2020) 109643

Table 2 Parameters values used in Figs. 2 and 3.

Declaration of Competing Interest

μ

β

λ

δ

γ

α

σ1

σ2

σ3

σ4

0.09

0.37

0.328

0.1

0.11

0.07

0.05

0.04

0.035

0.04

LV˜5 ≤ −M + ββ1 ε2 + 1 < 0. Case 6. Similarly to the third case, if (S, I, R, Z) ∈ D5 , with a large value of M one can get

LV˜5 ≤ −M + 1 < 0. Case 7. If (S, I, R, Z) ∈ D7 ∪ D8 , there is a sufficiently large number M such that

LV˜5 ≤ −M + 1 < 0. With the above choices of M, ε 1 and ε 3 and for a given ε 2 < 1, we get

LV˜5 (S, I, R, Z ) < 0,

for any

This concludes the proof.

(S, I, R, Z ) ∈ R4+ \ D.



6. Numerical simulations and conclusion Adopting Milstein’s scheme for SDE discretization [10], we simulated system (3) for many parameter sets with the aim of illustrating the results exhibited in this paper. We choose the initial value (S(0 ), I (0 ), R(0 )) = (0.5, 0.15, 0.1 ) and use R software to perform simulations. Theorem 3.1 presents conditions under which the disease dies out exponentially. This is illustrated by Fig. 1 where sufficient conditions are satisfied.

σ12 = 0.04 ≤ 2μ = 0.24 and √  βασ1 R0 min{μ + λ + δ, μ + γ , α}( R0 − 1 ) +  (μ + λ + δ ) 2μ − σ12

R0 = 0.7183 < 1,

= −7 · 10−4 < 0. When R0s exceeds 1, the disease should be prevalent with probability one. This confirms the behaviour of the epidemic in Fig. 2 with parameters given in Table 2 and a value of R0s equals to 1.0636. By sampling 10, 0 0 0 stochastic paths of the solution to system (3) and keeping the same parameter set presented in Table 2, the stationary probability density functions of infected and removed individuals at different times are shown in Fig. 3 respectively. We see that they are close to each other at t = 20 0 0 and t = 30 0 0, which is confirmed by Theorem 5.1. To conclude, a stochastic SIRI epidemic model with distributed delay is analyzed in this paper. The existence and uniqueness of a global positive solution is the first result we established. Some sufficient conditions are obtained to promise the extinction of the disease. R0s > 0 is an adequate condition to the permanence of the disease as well as the existence of unique stationary distribution to our model. As an important question, it will be noteworthy to study the behaviour of the considered model in the case where R0 ≤ 1 ≤ R0s . We leave the response to this question to further works. Authors’ statement We attest here that all the authors of the paper have contributed jointly to the research and writing of the paper approximately in the same amount. All the authors of the paper.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments The research of T. Caraballo has been partially supported by Ministerio de Ciencia, Innovación y Universidades (Spain), FEDER (European Community) under grant PGC2018-096540-B-I00. References [1] Berrhazi B-E, Fatini ME, Caraballo T, Pettersson R. A stochastic SIRI epidemic model with lévy noise. Discrete Contin Dyn Syst Ser B 2018;23(6):2415–31. [2] Berrhazi B-E, Fatini ME, Pettersson R, Laaribi A. Media effects on the dynamics of a stochastic SIRI epidemic model with relapse and lévy noise perturbation. Int J Biomath 2019;12(3):1950037. 21 [3] Caraballo T, Colucci R. A comparison between random and stochastic modeling for a SIR model. Commun Pure Appl Anal 2017;16(1):151–62. [4] Caraballo T, Fatini ME, Pettersson R, Taki R. A stochastic SIRI epidemic model with relapse and media coverage. Discrete Contin Dyn Syst Ser B 2018;23(8):3483–501. [5] Fatini ME, Boukanjime B. Stochastic analysis of a two delayed epidemic model incorporating lévy processes with a general non-linear transmission. Stoch Anal Appl 2019. doi:10.1080/07362994.2019.1680295. [6] Fatini ME, Khalifi ME, Gerlach R, Laaribi A, Taki R. Stationary distribution and threshold dynamics of a stochastics SIRS model with a general incidence. Physica A 2019;534:120696. [7] Fatini ME, Laaribi A, Pettersson R, Taki R. Lévy noise perturbation for an epidemic model with impact of media coverage. Stochastics 2019;91(7):998–1019. [8] Fatini ME, Lahrouz A, Pettersson R, Settati A, Taki R. Stochastic stability and instability of an epidemic model with relapse. Appl Math Comput 2018;316:326–41. [9] Fatini ME, Sekkak I, Laaribi A. A threshold of a delayed stochastic epidemic model with crowly-martin functional response and vaccination. Physica A 2019;520:151–60. [10] Higham DJ. An algorithmic introduction to numerical simulation of stochastic differential equations. SIAM Rev 2001;43:525–46. [11] Ji C, Jiang D, Yang Q, Shi N. Dynamics of a multigroup SIR epidemic model with stochastic perturbation. Automatica 2012;48(1):121–31. [12] Kermack WO, Kendrick AGM. A contribution to the mathematical theory of epidemics. Proc R Soc Lond A 1927;115:700–21. [13] Kermack WO, Kendrick AGM. Contributions to the mathematical theory of epidemics, part. II. Proc R Soc Lond A 1932;138:55–83. [14] Kermack WO, Kendrick AGM. Contributions to the mathematical theory of epidemics, part. III. Proc R Soc Lond A 1933;141:94–112. [15] Khasminskii R. Stochastic stability of differential equations. Springer; 2012. [16] Kutoyants AY. Statistical inference for ergodic diffusion processes. London: Springer; 2003. [17] Lahrouz A, Settati A, Fatini ME, Pettersson R, Taki R. Probability analysis of a perturbed epidemic system with relapse and cure. Int J Comput Methods 2018;16. doi:10.1142/s0219876218501402. [18] Liu Q, Jiang D, Hayat T, Alsaedi A. Dynamics of a stochastic SIR epidemic model with distributed delay and degenerate diffusion. J Franklin Inst 2019;356(13):7347–70. [19] Lu Q. Stability of sirs system with random perturbations. Physica A 2009;388(18):3677–86. [20] Ma W, Takeuchi Y, Hara T, Beretta E. Permanence of an SIR epidemic model with distributed time delays. Tohoku Math J 2002;54:581–91. [21] Macdonald N. Time lags in biological models. Lecture notes in biomathematics. Heidelberg: Springer-Verlag; 1978. [22] Shu H, Fan D, Wei J. Global stability of multi-group SEIR epidemic models with distributed delays and nonlinear transmission. Nonlinear Anal: RWA 2012;13:1581–92. [23] Tornatore E, Buccellato S, Vetro P. Stability of a stochastic sir system. Physica A 2005;354(C):111–26. [24] Xu C, Li X. The threshold of a stochastic delayed SIRS epidemic model with temporary immunity and vaccination. Chaos Solitons Fract 2018;111:227–34. [25] Yu J, Jiang D, Shi N. Global stability of two-group sir model with random perturbation. J Math Anal Appl 2009;360:235–44. [26] Zhang T, Meng X, Zhang T, Song Y. Global dynamics for a new high-dimensional SIR model with distributed delay. Appl Math Comput 2012;218:11806–19. [27] Zhao Y, Jiang D. The threshold of a stochastic sirs epidemic model with saturated incidence. Appl Math Lett 2014;34:90–3. [28] Zhu C, Yin G. Asymptotic properties of hybrid diffusion systems. SIAM J Control Optim 2007;46(4):1155–79.