Analysis of a Stochastic SICR Epidemic Model Associated with the Lévy Jump

Abstract

We propose and study a Susceptible-Infected-Confined-Recovered (SICR) epidemic model. For the proposed model, the driving forces include (for example) the Brownian motion processes and the jump Lévy noise. Usually, in the existing literature involving epidemiology models, the Lévy noise perturbations are ignored. However, in view of the presence of strong fluctuations in the SICR dynamics, it is worth including these perturbations in SICR epidemic models. Quite frequently, this results in several discontinuities in the processes under investigation. In our present study, we consider our SICR model after justifying its used form, namely, the component related to the Lévy noise. The existence and uniqueness of a global positive solution is established. Under some assumptions, we show the extinction and the persistence of the infection. In order to give some numerical simulations, we illustrate a new numerical method to validate our theoretical findings.

1. Introduction

In a large variety of ecosystems and in humans, which are the most affected species, epidemics are known to pose a significant threat to living organisms. Every year, millions of humans die from epidemics, which has motivated researchers in mathematical and biological sciences to develop control and mitigation strategies in the fight against epidemic diseases. Using various mathematical modeling techniques, epidemiologists have gained deeper insights into the dynamics by using this mathematical and cost-effective analytical setup, whereas traditional experimental techniques are time-consuming and more expensive. This has encouraged many mathematicians and biologists to develop competent epidemic models that can present a vivid picture of the reality (see, for details, [1,2,3,4,5,6]).

The first SIR model describing the dynamics of the three principal populations: the susceptible $\left(S\right)$, the infected $\left(I\right)$ and the recovered $\left(R\right)$, was studied by Kermack and McKendricks in 1927 (see [7]). Their proposed model played a cruciel role in initianting different research works in the field of disease dynamics. Based upon the natural history of hepatitis C, a recent study proposed a Susceptible–Infectious–Chronic–Recovered (SICR) type model (see [8]).

A natural phenomenon is always affected by environmental factors that can aggravate or mitigate the spread of the epidemic. The stochastic quantification of several real-life phenomena has been immensely helpful in understanding the random nature of their incidence or occurrence. It has also helped to find solutions to those problems that arise from it, either in the form of minimizing their undesirability or maximizing their rewards. In addition, infectious diseases are subject to chance and uncertainty in terms of the normal progression of the infection. Therefore, stochastic models are more suitable when compared to deterministic models, keeping in view the fact that stochastic systems consider not only the moving mean, but also the standard deviation behavior around it. On the other hand, deterministic systems generate similar results for fixed initial values, but stochastic systems may give different predicted results. Several stochastic infection models describing the effect of the Brownian motion on viral dynamics have been developed in the literature (see, for example, [9,10,11,12,13]; see also [14]).

Motivated essentially by the above-mentioned and other recent developments (see [15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32]), we propose here a mathematical model for the transmission of an SICR model. Our goal is to show the effect of the Lévy jump in population dynamics, in which the Lévy noise is used to describe contingency and outburst. It will, therefore, be interesting to consider the following stochastic model jointly driven by white noise and Lévy noise:

$$\left\{\begin{array}{cc}\hfill dS\left(t\right)=& \left[\lambda -\beta I\left(t\right)S\left(t\right)-\mu S\left(t\right)\right]dt-\sigma I\left(t\right)S\left(t\right)d{W}_t\hfill \\ \hfill & -{\int }_{U}Q\left(u\right)I(t-)S(t-)\check{N}(dt,du),\hfill \\ \hfill dI\left(t\right)=& \left[\beta I\left(t\right)S\left(t\right)-({\gamma }_1+{\gamma }_2+\mu )I\left(t\right)+ϵR\left(t\right)+\delta C\left(t\right)\right]dt+\sigma I\left(t\right)S\left(t\right)d{W}_t\hfill \\ \hfill & +{\int }_{U}Q\left(u\right)I(t-)S(t-)\check{N}(dt,du),\hfill \\ \hfill dC\left(t\right)=& \left[{\gamma }_{2}I\left(t\right)-(\delta +\mu +\mathfrak{d})C\left(t\right)\right]dt,\hfill \\ \hfill dR\left(t\right)=& \left[{\gamma }_{1}I\left(t\right)-(ϵ+\mu )R\left(t\right)\right]dt,\hfill \end{array}\right.$$

(1)

where $S\left(t\right)$ denotes the number of susceptible individuals at time t. The infected individuals with no clinical symptoms (that is, the virus is living or developing in the individuals, but without producing symptoms or only mild ones) are able to transmit the infection to other individuals $I\left(t\right)$. The infected individuals under treatment (that is, those in the so-called chronic stage) with a viral load remaining low $C\left(t\right)$, and $R\left(t\right)$ denotes the number of the recovered individuals at time t. The parameters of the SICR model (1) are given in

Table 1. Parameters and their meanings in the suggested SICR model.

Parameter	Meaning
$\lambda$	The recruitment rate
$\mu$	Natural death rate
$\beta$	The transmission rate
${\gamma }_2$	The average to the infected I individuals becoming chronically infected
${\gamma }_1$	Default treatment rate for I individuals
$ϵ$	The average of the recovered individuals returning infected
$\delta$	The average the chronically infected individuals returning infected
$\mathfrak{d}$	Death rate due to the infection in chronic stage

and the two-strain SEIR diagram is illustrated in Figure 1. Moreover, $W_t$ is a standard Brownian motion with intensity $\sigma$ defined on a complete filtered probability space:

$$\left(\mathsf{\Omega },\mathcal{F},{\left({\mathcal{F}}_t\right)}_{t\geqq 0},\mathbb{P}\right)$$

with the filtration ${\left({\mathcal{F}}_t\right)}_{t\geqq 0}$ satisfying the usual conditions. We denote by $S(t-)$, $I(t-)$, $C(t-)$ and $R(t-)$ the left limits of $S\left(t\right)$, $I\left(t\right)$, $C\left(t\right)$ and $R\left(t\right)$, respectively. $N(dt,du)$ is a Poisson counting measure with the stationary compensator $\nu \left(du\right)dt$ given by

$$\tilde{N}(dt,du)=N(dt,du)-\nu \left(du\right)dt,$$

where $\nu$ is defined on a measurable subset U of the non-negative half-line with $\nu \left(U\right)<\infty$. $Q\left(u\right)$ represents the jump intensity.

2. Existence and Uniqueness of the Global Positive Solution

In this section, we first define $\Delta$ as follows:

$$\Delta :=\left\{(S;I;C;R)\in {\mathbb{R}}_{+}^4:\frac{\lambda }{\mu +\mathfrak{d}}\leqq S+I+C+R\leqq \frac{\lambda }{\mu }\right\}$$

and we consider the following assumption:

$$\left(\mathrm{A}\right):0<Q\left(u\right)\leqq \frac{\mu }{\lambda }(u\in U).$$

Theorem 1.

With the initial condition $\left(S\left(0\right);I\left(0\right);C\left(0\right);R\left(0\right)\right)\in \Delta ,$ the system (1) has a unique global solution $\left(S\left(t\right);I\left(t\right);C\left(t\right);R\left(t\right)\right)\in \Delta$ for any $t\geqq 0$ $\mathrm{a}.\mathrm{s}.$

Moreover, each of the following inequalities holds true:

$${\displaystyle \underset{t\to \infty }{lim\; sup}S\left(t\right)\leqq \frac{\lambda }{\mu }\mathrm{a}.\mathrm{s}.and\underset{t\to \infty }{lim\; inf}S\left(t\right)\geqq \frac{\lambda }{\mu +\mathfrak{d}}\mathrm{a}.\mathrm{s}.,}$$

$${\displaystyle \underset{t\to \infty }{lim\; sup}I\left(t\right)\leqq \frac{\lambda }{\mu }\mathrm{a}.\mathrm{s}.and\underset{t\to \infty }{lim\; inf}I\left(t\right)\geqq \frac{\lambda }{\mu +\mathfrak{d}}\mathrm{a}.\mathrm{s}.,}$$

$${\displaystyle \underset{t\to \infty }{lim\; sup}C\left(t\right)\leqq \frac{\lambda }{\mu }\mathrm{a}.\mathrm{s}.and\underset{t\to \infty }{lim\; inf}C\left(t\right)\geqq \frac{\lambda }{\mu +\mathfrak{d}}\mathrm{a}.\mathrm{s}.,}$$

and

$${\displaystyle \underset{t\to \infty }{lim\; sup}R\left(t\right)\leqq \frac{\lambda }{\mu }\mathrm{a}.\mathrm{s}.and\underset{t\to \infty }{lim\; inf}R\left(t\right)\geqq \frac{\lambda }{\mu +\mathfrak{d}}\mathrm{a}.\mathrm{s}.}$$

Proof. 

Since the drift and the diffusion are local Lipschitzian, and in view of the initial condition $\left(S\left(0\right);I\left(0\right);C\left(0\right);R\left(0\right)\right)\in \Delta ,$ we prove the existence and uniqueness of the local solution $\left(S\left(t\right);I\left(t\right);C\left(t\right);R\left(t\right)\right)\in \Delta$ for $t\in [0,t_e)$, where $t_e$ is the explosion time.

In order to prove that this solution is global, we define the stopping time $\tau$ given by

$$\tau =\{t\in [0,t_e):S\left(t\right)\leqq 0,I\left(t\right)\leqq 0,C\left(t\right)\leqq 0\mathrm{and}R\left(t\right)\leqq 0\}.$$

If we suppose that $t_e<\infty$, then $\tau \leqq t_e$, so there exist $T>0$ and $ϵ>0$ such that $P(\tau \leqq T)>ϵ$.

Consider the function V on ${\mathbb{R}}_{+}^4$ defined by

$$V(S,I,C,R)=log\left(SICR\right).$$

Using the Itô’s formula, we get

$$\begin{array}{cc}\hfill dV& =LVdt+\left(\frac{\lambda }{S}-\beta I+\beta S-{\gamma }_1-{\gamma }_2-2\mu +ϵ\frac{R}{I}+\delta \frac{C}{I}\right)d{W}_t\hfill \\ \hfill & +{\int }_{U}log(1-QI)(1+QS)\check{N}(dt,du),\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill LV=& (\lambda -\beta IS-\mu S)\frac{1}{S}+(\beta IS-({\gamma }_1+{\gamma }_2+\mu )I+ϵR+\delta C)\frac{1}{I}-\frac{{\sigma }^2{I}^2}{2}\hfill \\ \hfill & -\frac{{\sigma }^2{S}^2}{2}+{\int }_U[log(1-QI)(1+QS)-QI+QS]\nu \left(du\right).\hfill \end{array}$$

Noting that $1-QI>0$ from Assumption (A), this fact implies that

$$\begin{array}{cc}\hfill LV& \geqq -\frac{\beta \lambda }{\mu }-2\mu -{\gamma }_1-{\gamma }_2-\frac{{\sigma }^2{\lambda }^2}{{\mu }^2}+{\int }_U[log(1-QI)+QI]\nu \left(du\right)\hfill \\ \hfill & +{\int }_U[log(1+QS)-QS]\nu \left(du\right)\hfill \\ \hfill & :=K.\hfill \end{array}$$

Since $x⟼log(1-x)+x$ and $x⟼log(1+x)-x$ are nonpositive functions, we find that

$$\begin{array}{cc}\hfill dV& \geqq Kdt+\left(\frac{\lambda }{S}-\beta I+\beta S-{\gamma }_1-{\gamma }_2-2\mu +ϵ\frac{R}{I}+\delta \frac{C}{I}\right)d{W}_t\hfill \\ \hfill & +{\int }_{U}log(1-QI)(1+QS)\check{N}(dt,du),\hfill \end{array}$$

which, upon integrating from 0 to t, yields

$$\begin{array}{cc}\hfill & V\left(S\left(t\right),I\left(t\right),C\left(t\right),R\left(t\right)\right)\geqq Kt+{\int }_0^t{\int }_U\left(\frac{\lambda }{S}-\beta I+\beta S-{\gamma }_1-{\gamma }_2-2\mu +ϵ\frac{R}{I}+\delta \frac{C}{I}\right)d{W}_s\hfill \\ \hfill & +V\left(S\left(0\right),I\left(0\right),C\left(0\right),R\left(0\right)\right)+{\int }_0^t{\int }_{U}log(1-QI)(1+QS)\check{N}(ds,du).\hfill \end{array}$$

Now, by using the continuity property, some components of $\left(S\left(\tau \right),I\left(\tau \right),C\left(\tau \right),R\left(\tau \right)\right)$ will be equal to 0, so we obtain

$$\underset{t\to \tau }{lim}V\left(\tau \right)=-\infty .$$

For $t\to \tau$, we get

$$\begin{array}{cc}\hfill -\infty \geqq & V\left(S\left(0\right),I\left(0\right),C\left(0\right),R\left(0\right)\right)+Kt\hfill \\ \hfill & +{\int }_0^t{\int }_U\left(\frac{\lambda }{S}-\beta I+\beta S-{\gamma }_1-{\gamma }_2-2\mu +ϵ\frac{R}{I}+\delta \frac{C}{I}\right)d{W}_s\hfill \\ \hfill & +{\int }_0^t{\int }_{U}log(1-QI)(1+QS)\check{N}(ds,du)>\infty ,\hfill \end{array}$$

which contradicts our hypothesis. Thus, clearly, we have shown that the solution of (1) is global.

Next, we will prove the boundedness of the solution. The equation of the system (1) implies that

$${\displaystyle \frac{dN\left(t\right)}{dt}=\lambda -\mu N\left(t\right)-dC\left(t\right),}$$

where

$$N\left(t\right)=S\left(t\right)+I\left(t\right)+C\left(t\right)+R\left(t\right).$$

Therefore, we have

$${\displaystyle \lambda -(\mu +\mathfrak{d})N\left(t\right)\leqq \frac{dN\left(t\right)}{dt}\leqq \lambda -\mu N\left(t\right),}$$

so that

$$\begin{array}{c}\hfill {\displaystyle e^{\mu t}\frac{dN\left(t\right)}{dt}\leqq e^{\mu t}\left(\lambda -\mu N\left(t\right)\right),}\end{array}$$

$${\displaystyle {\int }_0^t{e}^{\mu s}\frac{dN\left(s\right)}{ds}ds\leqq {\int }_0^t{e}^{\mu s}\left(\lambda -\mu N\left(s\right)\right)ds,}$$

$${\displaystyle e^{\mu t}N\left(t\right)\leqq \frac{\lambda }{\mu }\left(e^{\mu t}-1\right)+N\left(0\right),}$$

$${\displaystyle N\left(t\right)\leqq \frac{\lambda }{\mu }(1-e^{-\mu t})+N\left(0\right)e^{-\mu t}}$$

and

$$\begin{array}{c}\hfill {\displaystyle \underset{t\to \infty }{lim\; sup}N\left(t\right)\leqq \frac{\lambda }{\mu }\mathrm{a}.\mathrm{s}.}\end{array}$$

Using the same technique as above, we get

$${\displaystyle \underset{t\to \infty }{lim\; inf}N\left(t\right)\geqq \frac{\lambda }{\mu +\mathfrak{d}}\mathrm{a}.\mathrm{s}.}$$

This proves the boundedness. □

The above result shows the biologically well-posedness of our model and also that the solution is ultimately bounded such that

$${\displaystyle \underset{t\to \infty }{lim\; sup}N\left(t\right)\leqq \frac{\lambda }{\mu }\mathrm{a}.\mathrm{s}.}$$

and

$${\displaystyle \underset{t\to \infty }{lim\; inf}N\left(t\right)\geqq \frac{\lambda }{\mu +\mathfrak{d}}\mathrm{a}.\mathrm{s}.}$$

3. Extinction of the Infection

In this section, we will prove the extinction of the infection under some sufficient condition. We first define $M_1$ by

$$M_1:=({\gamma }_1+{\gamma }_2+\mu )+(ϵ+\delta )\frac{\lambda }{\mu }-\frac{{\beta }^2}{2{\sigma }^2}.$$

Theorem 2.

If $M_1>0,$ then $I\left(t\right)\to 0$ when $t\to \infty$ a.s.

Proof. 

Let

$$V\left(I\right)=log\left(I\right).$$

Then, by using Itô’s formula, we obtain

$$\begin{array}{c}\hfill dV\left(t,X\left(t\right)\right)=LV\left(t,X\left(t\right)\right)dt+\sigma S\left(t\right)d{W}_t+{\int }_U\left(log(1+QS)\right)\check{N}(dt,du),\end{array}$$

with

$$LV=(\beta IS-({\gamma }_1+{\gamma }_2+\mu )I+ϵA+\delta C)\frac{1}{I}-\frac{{\sigma }^2{S}^2}{2}+{\int }_{U}log(1+QS)-QS]\nu \left(du\right),$$

so we have

$$LV\leqq \frac{{\beta }^2}{2{\sigma }^2}-({\gamma }_1+{\gamma }_2+\mu )+(ϵ+\delta )\frac{\lambda }{\mu },$$

so that

$$dV\leqq \left(\frac{{\beta }^2}{2{\sigma }^2}-({\gamma }_1+{\gamma }_2+\mu )+(ϵ+\delta )\frac{\lambda }{\mu }\right)dt+\sigma S\left(t\right)d{W}_t+{\int }_{U}log(1+QS)\check{N}(dt,du).$$

Therefore, we get

$$\begin{array}{cc}\hfill {\displaystyle \frac{logI\left(t\right)}{t}}& {\displaystyle \leqq \frac{log\left(I_0\right)}{t}+\frac{1}{t}{\int }_0^t\left(\frac{{\beta }^2}{2{\sigma }^2}-({\gamma }_1+{\gamma }_2+\mu )+(ϵ+\delta )\frac{\lambda }{\mu }\right)ds+\frac{1}{t}{\int }_0^t\sigma S\left(s\right)dW\left(s\right)}\hfill \\ \hfill & {\displaystyle +\frac{1}{t}{\int }_0^t{\int }_{U}log(1+QS)\check{N}(ds,du)}\hfill \\ \hfill & {\displaystyle \leqq \frac{log\left(I_0\right)}{t}+\frac{{\beta }^2}{2{\sigma }^2}-({\gamma }_1+{\gamma }_2+\mu )+(ϵ+\delta )\frac{\lambda }{\mu }+\frac{\sigma }{t}{\int }_0^{t}S\left(s\right)dW\left(s\right)}\hfill \\ \hfill & {\displaystyle +\frac{1}{t}{\int }_0^t{\int }_{U}log(1+QS)\check{N}(ds,du).}\hfill \end{array}$$

If we now define $M_t$ by

$$M_t:={\int }_0^t\sigma S\left(s\right)d{W}_s,$$

then

$$\begin{array}{c}\hfill {\displaystyle \underset{t\to \infty }{lim\; sup}\frac{\langle M_t,M_t\rangle }{t}=\underset{t\to \infty }{lim\; sup}\frac{{\sigma }^2}{t}{\int }_0^t{u}^2\left(s\right)ds\leqq {\sigma }^2{\left({\displaystyle \frac{\lambda }{\mu }}\right)}^2<\infty .}\end{array}$$

Thus, by using the fact that the solution of the system (1) is bounded and also the strong law of large numbers theorem for martingales, we have (see [33])

$$\begin{array}{c}\hfill {\displaystyle \underset{t\to \infty }{lim\; sup}\frac{M_t}{t}=0,}\end{array}$$

so that

$$\begin{array}{cc}\hfill {\displaystyle \underset{t\to \infty }{lim\; sup}\frac{logI\left(t\right)}{t}}& \leqq \frac{{\beta }^2}{2{\sigma }^2}-({\gamma }_1+{\gamma }_2+\mu )+(ϵ+\delta )\frac{\lambda }{\mu }.\hfill \end{array}$$

So, if

$${\displaystyle \frac{{\beta }^2}{2{\sigma }^2}<({\gamma }_1+{\gamma }_2+\mu )+(ϵ+\delta )\frac{\lambda }{\mu }}$$

when $M_1>0$, then $I\left(t\right)\to 0$ when $t\to \infty$ a.s. □

This result shows that the disease dies out when we have the critical threshold of a possibly greater amplitude of the volatility $\sigma$, which is described as follows:

$$\frac{{\beta }^2}{2{\sigma }^2}<({\gamma }_1+{\gamma }_2+\mu )+(ϵ+\delta )\frac{\lambda }{\mu }.$$

4. Persistence in the Mean

Firstly, we define

$${\displaystyle \langle x\left(t\right)\rangle :=\frac{1}{t}{\int }_0^{t}x\left(s\right)ds,}$$

then the persistence in the mean of $x\left(t\right)$ is defined as follows:

$$\underset{t\to \infty }{lim\; inf}\langle x\left(t\right)\rangle >0.$$

Now, we will give some condition to prove the persistence of $S\left(t\right)$, $I\left(t\right)$, $C\left(t\right)$ and $R\left(t\right)$ in the mean. For this purpose we define $M_2$ by

$$M_2:=\beta \frac{\lambda }{\mu +\mathfrak{d}}-({\gamma }_1+{\gamma }_2+\mu )-\frac{{\sigma }^2{\lambda }^2}{2{\mu }^2}.$$

Theorem 3.

If, for any solution $\left(S\left(t\right),I\left(t\right),C\left(t\right),R\left(t\right)\right)$ of the system (1), $M_2>0,$ then

$$\begin{array}{c}\hfill \underset{t\to \infty }{lim\; inf}\langle S\left(t\right)\rangle \geqq \frac{\lambda \mu }{\lambda \beta +{\mu }^2}>0,\end{array}$$

$$\underset{t\to \infty }{lim\; inf}\langle I\left(t\right)\rangle \geqq \frac{M_2}{({\gamma }_1+{\gamma }_2+\mu )}>0,$$

$$\underset{t\to \infty }{lim\; inf}\langle C\left(t\right)\rangle \geqq \frac{{\gamma }_2{M}_2}{({\gamma }_1+{\gamma }_2+\mu )(\delta +\mu +\mathfrak{d})}>0$$

and

$$\underset{t\to \infty }{lim\; inf}\langle R\left(t\right)\rangle \geqq \frac{{\gamma }_1{M}_2}{({\gamma }_1+{\gamma }_2+\mu )(ϵ+\mu )}>0.$$

Proof. 

Using the first equation of the system (1), we have

$$\begin{array}{cc}\hfill {\displaystyle \frac{S\left(t\right)-S\left(0\right)}{t}}& =\frac{1}{t}{\int }_0^t(\lambda -\beta IS-\mu S)ds-\frac{1}{t}{\int }_0^t\sigma ISd{W}_s-\frac{1}{t}{\int }_0^t{\int }_{U}Q\left(u\right)IS\check{N}(ds,du)\hfill \\ \hfill & \geqq \frac{1}{t}{\int }_0^t(\lambda -\beta \frac{\lambda }{\mu }S-\mu S)ds-\frac{1}{t}{\int }_0^t\sigma ISd{W}_s-\frac{1}{t}{\int }_0^t{\int }_{U}Q\left(u\right)IS\check{N}(ds,du)\hfill \end{array}$$

and

$$\begin{array}{c}\hfill {\displaystyle \left(\frac{\lambda \beta }{\mu }+\mu \right)\langle S\left(t\right)\rangle \geqq \lambda -\frac{S\left(t\right)-S\left(0\right)}{t}-\frac{1}{t}{\int }_0^t\sigma ISd{W}_s-\frac{1}{t}{\int }_0^t{\int }_{U}Q\left(u\right)IS\check{N}(ds,du).}\end{array}$$

Thus, by the boundedness of solution and the strong law of large numbers for the martingales, we obtain

$$\underset{t\to \infty }{lim\; inf}\langle S\left(t\right)\rangle \geqq \frac{\lambda \mu }{\lambda \beta +{\mu }^2}>0.$$

Next, by integrating the second equation of the system (1) from 0 to t and dividing both sides by t, we obtain

$$\begin{array}{cc}\hfill {\displaystyle \frac{I\left(t\right)-I\left(0\right)}{t}}& =\frac{1}{t}{\int }_0^t\left[\beta I\left(s\right)S\left(s\right)-({\gamma }_1+{\gamma }_2+\mu )I\left(s\right)+ϵR\left(s\right)+\delta C\left(s\right)\right]ds+\frac{1}{t}{\int }_0^t\sigma ISd{W}_s\hfill \\ \hfill & +\frac{1}{t}{\int }_0^t{\int }_{U}Q\left(u\right)IS\check{N}(ds,du)\hfill \\ \hfill & \geqq -({\gamma }_1+{\gamma }_2+\mu )\langle I\left(t\right)\rangle +\frac{1}{t}{\int }_0^t\sigma ISd{W}_s+\frac{1}{t}{\int }_0^t{\int }_{U}Q\left(u\right)IS\check{N}(ds,du).\hfill \end{array}$$

Using the Itô’s formula on $V\left(I\right)=log\left(I\right)$, we have

$$\begin{array}{cc}\hfill dV& =(\beta IS-({\gamma }_1+{\gamma }_2+\mu )I+ϵR+\delta C)\frac{1}{I}-\frac{{\sigma }^2{S}^2}{2}\hfill \\ \hfill & +{\int }_{U}log(1+QS)-QS]\nu \left(du\right)dt\hfill \\ \hfill & +\sigma S\left(t\right)d{W}_t+{\int }_{U}log(1+QS)\check{N}(dt,du)\hfill \\ \hfill & \geqq \left(\beta \frac{\lambda }{\mu +\mathfrak{d}}-({\gamma }_1+{\gamma }_2+\mu )+ϵ\frac{\mu }{\mu +\mathfrak{d}}+\delta \frac{\mu }{\mu +\mathfrak{d}}\right)-\frac{{\sigma }^2{S}^2}{2}\hfill \\ \hfill & +{\int }_{U}log(1+QS)-QS]\nu \left(du\right)dt\hfill \\ \hfill & +\sigma S\left(t\right)d{W}_t\hfill \\ \hfill & +{\int }_{U}log(1+QS)\check{N}(dt,du)\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill \frac{logI\left(t\right)-logI\left(0\right)}{t}& \geqq (\beta \frac{\lambda }{\mu +\mathfrak{d}}-({\gamma }_1+{\gamma }_2+\mu )+ϵ\frac{\mu }{\mu +\mathfrak{d}}+\delta \frac{\mu }{\mu +\mathfrak{d}})-\frac{{\sigma }^2{\lambda }^2}{2{\mu }^2}\hfill \\ \hfill & +\frac{1}{t}{\int }_{U}log(1+QS)-QS]\nu \left(du\right)ds+\frac{1}{t}{\int }_0^t\sigma S\left(s\right)d{W}_s\hfill \\ \hfill & +\frac{1}{t}{\int }_0^t{\int }_{U}log(1+QS)\check{N}(ds,du).\hfill \end{array}$$

Now, if we sum

$${\displaystyle \frac{logI\left(t\right)-logI\left(0\right)}{t}\mathrm{and}\frac{I\left(t\right)-I\left(0\right)}{t}}$$

by using the positivity and the boundedness of the solution as well as the strong law of large numbers for the martingales, we find that

$$\begin{array}{c}\hfill \underset{t\to \infty }{lim\; inf}\langle I\left(t\right)\rangle \geqq \frac{1}{{\gamma }_1+{\gamma }_2+\mu }\left(\frac{\beta \lambda }{\mu +\mathfrak{d}}-({\gamma }_1+{\gamma }_2+\mu )-\frac{{\sigma }^2{\lambda }^2}{2{\mu }^2}\right)>0.\end{array}$$

Hence, we have

$$\begin{array}{c}\hfill \underset{t\to \infty }{lim\; inf}\langle I\left(t\right)\rangle \geqq \frac{M_2}{{\gamma }_1+{\gamma }_2+\mu }>0.\end{array}$$

Consequently, the persistence in the mean of $I\left(t\right)$ is proved.

Using the third equation of the model (1), we have

$$\begin{array}{ccc}\hfill {\displaystyle \frac{C\left(t\right)-C\left(0\right)}{t}}& =\frac{1}{t}{\int }_0^t\left({\gamma }_{2}I\left(s\right)-(\delta +\mu +\mathfrak{d})C\left(s\right)\right)ds\hfill & \hfill ={\gamma }_2\langle I\left(t\right)\rangle +(\delta +\mu +\mathfrak{d})\langle C\left(t\right)\rangle ,\end{array}$$

so that

$$\begin{array}{c}\hfill \underset{t\to \infty }{lim\; inf}\langle C\left(t\right)\rangle \geqq \frac{{\gamma }_2{M}_2}{({\gamma }_1+{\gamma }_2+\mu )(\delta +\mu +\mathfrak{d})}>0.\end{array}$$

This last equation of the model (1) implies that

$$\begin{array}{cc}\hfill {\displaystyle \frac{R\left(t\right)-R\left(0\right)}{t}}& =\frac{1}{t}{\int }_0^t\left({\gamma }_{1}I\left(s\right)-(ϵ+\mu )R\left(s\right)\right)ds,\hfill \\ \hfill & ={\gamma }_1\langle I\left(t\right)\rangle +(ϵ+\mu )\langle R\left(t\right)\rangle .\hfill \end{array}$$

Therefore, we finally obtain

$$\begin{array}{c}\hfill {\displaystyle \underset{t\to \infty }{lim\; inf}\langle R\left(t\right)\rangle \geqq \frac{{\gamma }_1{M}_2}{({\gamma }_1+{\gamma }_2+\mu )(ϵ+\mu )}>0.}\end{array}$$

□

The theorem involving the persistence in the mean has been proved under the following sufficient condition:

$$({\gamma }_1+{\gamma }_2+\mu )+\frac{{\sigma }^2{\lambda }^2}{2{\mu }^2}<\beta \frac{\lambda }{\mu +\mathfrak{d}}.$$

This means that, with an adopted smallest magnitude of volatility $\sigma$, the model is persistent in the mean.

5. Numerical Results

In order to illustrate the numerical simulations of the model (1), we consider the following problem:

$$dX\left(t\right)=f\left(t,X\left(t\right)\right)dt+\sigma \left(t,X\left(t\right)\right)d{W}_t+{\int }_{U}Q(t,u)g\left(t-,X(t-)\right)\check{N}(dt,du).$$

(2)

The solution of the system (2) will now be given by

$$\begin{array}{cc}\hfill X\left(t\right)& =X\left(0\right)+\underset{\mathrm{Part}1}{\underbrace{{\int }_0^{t}f\left(s,X\left(s\right)\right)ds+{\int }_0^t\sigma \left(s,X\left(s\right)\right)d{W}_s}}\hfill \\ \hfill & +\underset{\mathrm{Part}2}{\underbrace{g\left(t-,X(t-)\right){\int }_0^t{\int }_{U}Q(s,u)\check{N}(ds,du)}},\hfill \end{array}$$

(3)

The Milsteins Higher Order Method consists of approximating Part 1 in the system (2) in the following form:

$$\begin{array}{cc}\hfill X_{k+1}& =X_k+f(t_k,X_k)\Delta t+\sigma (t_k,X_k)(W_{k+1}-W_k)\hfill \\ \hfill & {\displaystyle +\frac{1}{2}{\sigma }^2(t_k,X_k)\left[{(W_{k+1}-W_k)}^2-(t_{k+1}-t_k)\right].}\hfill \end{array}$$

About the approximation of Part 2, we have two cases. Consider any infinitesimal interval $[T_i,T_{i+1})\subset ({\tau }_i,{\tau }_{i+1})$.

There is no jump point in this interval:

$$g\left(t-,X(t-)\right){\int }_{T_i}^{T_{i+1}}{\int }_{U}Q(s,u)\check{N}(ds,du)=0.$$

If there is only one jump point ${\tau }_i\in [T_i,T_{i+1}),$ then

$$g\left(t-,X(t-)\right){\int }_{T_i}^{T_{i+1}}{\int }_{U}Q\left(u\right)\check{N}(ds,du)=g\left({\tau }_i-,X({\tau }_i-)\right)Q\left({\tau }_i,\zeta \left({\tau }_i\right)\right),$$

so that

$$g\left(t-,X(t-)\right){\int }_0^t{\int }_{U}Q(s,u)\check{N}(ds,du)=\sum _{i=0}^{i=n}g\left({\tau }_i-,X({\tau }_i-)\right)Q\left({\tau }_i,\zeta \left({\tau }_i\right)\right).$$

Therefore, the Milsteins Higher Order Method of the system (2) will be given by

$$\begin{array}{cc}\hfill X_{k+1}& =X_k+f\left(t_k,X_k)\right)\Delta t+\sigma (t_k,X_k)(W_{k+1}-W_k)\hfill \\ \hfill & {\displaystyle +\frac{1}{2}{\sigma }^2(t_k,X_k)\left({(W_{k+1}-W_k)}^2-\Delta t\right)+g\left(t_k-,X\left(t_k\right)\right)Q\left(t_k,\zeta \left(t_k\right)\right)}\hfill \end{array}$$

We now use the previous method to solve the system (1).

Figure 2 and Figure 3 show the dynamics of the susceptible population for the case of the extinction of the infection. From these figures, we observe clearly that the infected compartment converges to 0. It should be noted that, in this case, the susceptible class increases to attain their maximum, which means the the disease dies out. This is consistent with our theoretical results for the extinction of the infection.

The behavior of the infection for both the infected and susceptible population with the Lévy jump process is graphically illustrated in Figure 4 and Figure 5 in the case of the persistence of the infection. In this epidemic scenario, we note from all of the four SICR compartments, that is, the infected, the susceptible, the chronic sub-population and the recovered, that the infection persists. This shows that our theoretical and numerical results are compatible.

6. Conclusions

In our present investigation, we have considered the stochastic epidemic SICR model which is driven by the Brownian motion and the Lévy noise jointly so as to better describe the sudden social fluctuations. We have also performed some adequate numerical simulations not only to support our theoretical results, but also for predicting the asymptotic behavior of the solutions of the corresponding system in this study. More precisely, first of all, by using the Lyapunov analysis method, we have demonstrated the existence and uniqueness of the derived solutions. In addition, we have shown the positivity and the boundedness:

$${\displaystyle \underset{t\to \infty }{lim\; sup}N\left(t\right)\leqq \frac{\lambda }{\mu }\mathrm{a}.\mathrm{s}.}$$

and

$${\displaystyle \underset{t\to \infty }{lim\; inf}N\left(t\right)\geqq \frac{\lambda }{\mu +\mathfrak{d}}\mathrm{a}.\mathrm{s}.}$$

The extinction of the infection is established with a critical threshold of a possible greater amplitude of the volatility $\sigma$, which is described as follows:

$$M_1=({\gamma }_1+{\gamma }_2+\mu )+(ϵ+\delta )\frac{\lambda }{\mu }-\frac{{\beta }^2}{2{\sigma }^2}>0.$$

The persistence of the infection has been proved under the following sufficient condition:

$$M_2=\beta \frac{\lambda }{\mu +\mathfrak{d}}-({\gamma }_1+{\gamma }_2+\mu )-\frac{{\sigma }^2{\lambda }^2}{2{\mu }^2}>0.$$

Finally, we have constructed a modified Milsteins Higher Order Method, which is implemented in order to validate our theoretical findings in relation to the model considered, with the aim to provide additional information to help the decision makers choose a good disease control strategy. On the one hand, we have decreased or increased the intensity of fluctuations. On the other hand, we have considered the effect of the Lévy noise in the development of the system variables. In possible sequels to our investigation in this paper, one can extend the model (1) to various fractional-order models, which are based upon several operators of fractional derivatives such as those known as the Riemann–Liouville, Liouville–Caputo and other fractional derivatives (see, for example, [34,35]) and also extend our model to an impulsive model such as that considered in [36].