Global dynamics and density function in a class of stochastic SVI epidemic models with Lévy jumps and nonlinear incidence

1.   Introduction

As is well known, in real life, there are many infectious diseases that seriously threaten human health. Especially, in recent decades, the repeated epidemic of infectious diseases has brought great disasters to human survival and the national economy and people's livelihood.

It has been confirmed that vaccination is one of the most practical and efficient strategies to prevent and control the spread of many diseases, such as measles, pertussis, influenza, Hepatitis B virus (HBV) and human tuberculosis (TB) (see ^[1,2,3,4]). The spectacular successful cases were seen to be the eradication of small-pox in 1977 ^[5], the control of poliomyelitis and measles throughout Central and South America ^[6,7], and in the United Kingdom the vaccination campaign against measles in 1994 ^[8]. In order to analyze the dynamical properties of vaccination, in recent years various types of epidemic dynamical models with vaccination are established and investigated widely (see^{[9,10,11,12,13,14,15]} and the references cited therein).

On the other hand, to the best of our knowledge, in modeling the dynamics of epidemic systems the incidence rate is an important substance. In many practicalities, such as media reports, vaccination, quarantine, catch and kill, protection, and population density, which may directly or indirectly influence the incidence rate. At this time, the nonlinear incidence, such as the saturated incidence $\frac{\beta SI}{1+\alpha I}$, Beddington-DeAngelis incidence $\frac{\beta SI}{1+\omega S+\alpha I}$, nonlinear incidences $\beta Sg(I)$ and $\beta f(S, I)$, is more realistic and achieving more exact results (see ^{[15,16,17,18,19,20,21,22,23]}and the references cited therein).

Motivated by the previous works, this paper describes the effects of vaccination prevention strategies for the newly susceptible individuals and the vaccinated can also be affected under the nonlinear incidence of disease, we propose the following deterministic Susceptible-Vaccinated-Infected (SVI) epidemic model with nonlinear incidence:

where the definitions of all state variables, parameters and functions in model (1.1) are listed in Tables 1–3. It will be seen below that the basic reproduction number $R_0$ of the model (1.1) depends directly on the vaccination rate and nonlinear incidence, and that the main dynamical properties of the model, such as the stability of equilibria, the extinction and persistence of the disease, etc., are fully determined by $R_0$.

However, the results of studies over the past few years have shown that birth and death rates are more or less influenced by random white noise during disease transmission. As a result, a growing number of authors have studied the associated stochastic epidemic models (see ^{[13,22,23,24,25,26,27,28,29,30,31]} and the references cited therein). The main research subjects include the global existence of a positive solution with any positive initial value in probability, the persistence and extinction of the disease in probability, the asymptotical behaviors in probability around the disease-free and endemic equilibria of the corresponding deterministic models, the existence of stationary distribution as well as ergodicity, the expressions of probability density functions, etc. Especially, in articles ^[13,18,31] the stochastic epidemic models with vaccination are proposed and studied.

From the perspective of epidemiology, the existence and ergodicity of stationary distribution indicates that an infectious disease will prevail and persist for a long time. More importantly, the corresponding probability density function of the stationary distribution can reflect all statistical properties of different compartment individuals. It can be regarded as a great intersection of epidemiological dynamics and statistics. We see that recently the probability density functions for stochastic epidemic models are studied in articles ^{[25,29,30,31]} by solving the corresponding algebraic equations which are equivalent to the Fokker-Planck equation. It should be pointed out that until now there are still relatively few works devoted to the expressions of probability density functions due to the difficulty of solving the corresponding Fokker-Planck equation.

It is a pity that the stochastic model with only white noise cannot reasonably describe some random disturbance in the actual environment, such as the outbreak of bird flu and SARS, earthquakes, hurricanes, floods, discharge of toxic pollutants, etc., this is because these processes are discontinuous. In order to accurately describe these phenomena, it is a feasible and more realistic method to introduce the Lévy jumps noise to the original basic dynamical model. We see that a lot of research has been done to direct at the epidemic models with Lévy jumps noise (see ^{[32,33,34,35,36,37,38,39,40,41]} and the references cited therein). Particularly, in articles ^[34,41], the stochastic epidemic models with vaccination and Lévy jumps are proposed and studied. It was shown that the white noises and Lévy jumps could make the stationary distribution vanish as well as appear.

Motivated by the above works, in order to describe the effects of Lévy jumps in the transmission of disease under the vaccination prevention strategies, on the basis of model (1.1), in this paper we propose the following stochastic SVI models with white noise, Lévy jumps and nonlinear incidence:

In this paper, we always assume that model (1.2) is defined on a complete probability space $(\Omega, \{\mathcal{F}\}_{t\geq0}, P)$ with a filtration $\{\mathcal{F}\}_{t\geq0}$ satisfying the usual conditions (that is to say, it is increasing and right continuous while $\mathcal{F}_{0}$ contains all $P$-null sets). In model (1.2), $B_i(t)\; (i = 1, 2, 3)$ are standard Brownian motion, $\sigma_{i}\; (i = 1, 2, 3)$ are the intensities of $B_i(t)$; $\widetilde{W}(t, u)$ is the compensated Poisson random measure with characteristic measure $\nu$ on $Z$, where $Z$ is a measurable subset of $(0, +\infty)$ with the measure $\nu(Z) < \infty$. $\widetilde{W}(t, u) = W(t, u)-t\nu(u)$, and $W(t, u)$ is a Poisson counting measure with characteristic measure $\nu$ which is defined on $Z$ and often used to describe jumps process. $\eta_{i}(u)\; (i = 1, 2, 3)$ are the density of jumps process defined for all $u\in Z$. $S(t-)$, $V(t-)$ and $I(t-)$ are the left limit of $S(t)$, $V(t)$ and $I(t)$, respectively.

For the jumps process in model (1.2), we always require the following assumptions (see ^[39]):

(H1) $\eta_{i}(u)$ is continuous function for $u\in Z$, and $\int_{Z}\eta^2_i(u)v(du) < \infty, \; i = 1, 2, 3.$

(H2) $1+\eta_{i}(u) > 0$ for all $u\in Z$, and $\int_{Z}[\eta_{i}(u)-\ln(1+\eta_{i}(u))]v(du) < \infty, \; i = 1, 2, 3.$

The main purpose in this paper is to investigate the stochastic dynamics of model (1.2) including the stochastic extinction and persistence of disease, the existence of ergodic stationary distribution and expressions of probability density functions, which are corresponding to the global stability of disease-free and endemic equilibria and the uniform persistence of positive solutions for corresponding deterministic model (1.1). The main contribution and innovations are summarized as follows:

(1) The basic reproduction number $R_0$ of deterministic model (1.1) is calculated, which depends directly on the vaccination and nonlinear incidence. The stability of disease-free and endemic equilibria is fully determined by $R_0$.

(2) The threshold value $R_0^s$ of stochastic model (1.2) is defined which depends on the white noise, Lévy jumps, vaccination and nonlinear incidence. The threshold criteria for the stochastic extinction and persistence in the mean of the disease are established.

(3) The existence and ergodicity of stationary distribution for model (1.2) are obtained by constructing a suitable Lyapunov function, and it is determined by threshold value $R_0^s$.

(4) The expression of a log-normal density function around the quasi-endemic equilibrium of stochastic model is calculated, and a new technique for the calculation of probability density function is proposed.

The rest of this article is organized as follows. In Section 2, the dynamical behavior for model (1.1) is discussed. In Sections 3 and 4, we investigate the stochastic extinction and persistence in the mean of the disease, and the existence of stationary distribution for model (1.2), respectively. In Section 5, the criterion on the existence of log-normal probability density function of model (1.2) is established, here we will adopt a new technique for the calculation of density function. Furthermore, in Section 6, we present the numerical simulations to support the main results obtained in this paper. Lastly, in Section 7, we give a brief conclusion.

2.   Stability in deterministic model (1.1)

The initial condition of any solution for model (1.1) is given by

Based on the fundamental theory of ordinary differential equations, it is easy to get that the unique nonnegative solution $(S(t), V(t), I(t))$ for any initial value $(\hat{S}_0, \hat{V}_0, \hat{I}_0)\in R_+^3$ model (1.1) is defined for all $t\geq 0$. Moreover, in the light of model (1.1) we have

This implies $\limsup_{t\rightarrow \infty}(S+V+I)\leq \frac{\Lambda}{\mu}$, and the set

is a positive invariant set of model (1.1). This shows that, without loss of generality, we only need to consider the solutions of model (1.1) in the region $\Pi$.

Model (1.1) always has a unique disease-free equilibrium $P_{0} = (S_{0}, V_{0}, 0)$ with $S_0 = \frac{(1-\pi)\Lambda}{\mu}$ and $V_0 = \frac{\pi\Lambda}{\mu}$. By using the next generation matrix method we can obtain that the basic reproduction number of model (1.1) is

When $R_{0} > 1$, we can easily obtain that model (1.1) has a unique endemic equilibrium $P^* = (S^*, V^*, I^*)$ with $S^* = \frac{(1-\pi)\Lambda}{\mu+\beta_{b}f(I^*)}$, $V^* = \frac{\pi\Lambda}{\mu+\beta_{v}g(I^*)}$ and $I^*$ is the unique positive solution of the equation

Furthermore, we can build the following result.

Theorem 2.1. For model (1.1), the following conclusions hold.

$(i)$ If $R_{0}\leq 1$, then disease-free equilibrium $P_{0}$ is globally asymptotically stable.

$(ii)$ If $R_{0} > 1$, then $P_{0}$ is unstable and endemic equilibrium $P^*$ is locally asymptotically stable.

Proof. $(i)$ In view of model (1.1), we have $dS\leq [(1-\pi) \Lambda-\mu S]dt$ and $dV\leq [\pi \Lambda-\mu V]dt$, which implies that $\limsup_{t\rightarrow \infty}S\leq\frac{(1-\pi)\Lambda}{\mu}: = S_{0}$ and $\limsup_{t\rightarrow \infty}V\leq\frac{\pi\Lambda}{\mu}: = V_{0}$. Choose a Lyapunov function $U(t) = \frac{1}{2}I^{2}(t)$, then when $R_{0}\leq1$ we have

for any $(S, V, I)\in \Pi$. Furthermore, we easily prove that the maximal invariant set in $\{(S, V, I)\in \Pi: \frac{dU(t)}{dt} = 0\}$ is equilibrium $P_0$. Therefore, by the LaSalle's invariant principle, $P_0$ is globally asymptotically stable.

(ii) Calculating the Jacobi matrix of model (1.1) at equilibrium $P_0$ implies the characteristic equation: $(\lambda+\mu)^2(\lambda-(\mu+\delta)(R_{0}-1)) = 0.$ When $R_{0} > 1$, it is clear that $J(P_{0})$ has an eigenvalue $\lambda = (\mu+\delta)(R_{0}-1) > 0$. Hence, $P_{0}$ is unstable.

The Jacobi matrix of model (1.1) at equilibrium $P^*$ is

where $\hat{l}_{11} = \mu+\beta_{b}f(I^*)$, $\hat{l}_{13} = \beta_{b}S^*f'(I^*)$, $\hat{l}_{22} = \mu+\beta_{v}g(I^*)$, $\hat{l}_{23} = \beta_{v}V^*g'(I^*)$, $\hat{l}_{31} = \beta_{b}f(I^*)$, $\hat{l}_{32} = \beta_{v}g(I^*)$ and $\hat{l}_{33} = -\beta_{b}S^*f'(I^*)-(\beta_{v}V^*g'(I^*)+(\mu+\delta) > -\frac{\beta_{b}S^*f(I^*)}{I^*}-\frac{\beta_{v}V^*g(I^*)}{I^*}+(\mu+\delta) = 0$. By directly calculating, we can obtain the characteristic equation of $J(P^*)$

where $\hat{l}_{1} = \hat{l}_{11}+\hat{l}_{22}+\hat{l}_{33} > 0$, $\hat{l}_{2} = \hat{l}_{11}(\hat{l}_{22}+\hat{l}_{33})+\hat{l}_{22}\hat{l}_{33}+\hat{l}_{23}\hat{l}_{32}+\hat{l}_{13}\hat{l}_{31} > 0$ and $\hat{l}_{3} = \hat{l}_{11}(\hat{l}_{22}\hat{l}_{33}+\hat{l}_{23}\hat{l}_{32})+\hat{l}_{13}\hat{l}_{22}\hat{l}_{31} > 0$. Since

owing to the Routh-Hurwitz criterion, we can obtain that $P^*$ is locally asymptotically stable. This completes the proof.

Remark 2.1. The conclusion $(i)$ of Theorem 2.1 shows that the disease in model (1.1) is extinct. Furthermore, from conclusion $(ii)$ of Theorem 2.1 and by using the persistence theory of dynamical systems we can easily prove that when $R_0 > 1$ the disease in model (1.1) is uniformly persistent, that is, there is a constant $m > 0$ such that for any positive solution $(S(t), V(t), I(t))$ of model (1.1), one has $\liminf_{t\to\infty}(S(t), V(t), I(t))\geq (m, m, m)$. However, it is regrettable that when $R_0 > 1$ we can not get the global stability of $P^*$. Therefore, this will be an interesting open problem.

3.   Extinction and persistence in the mean

Firstly, as the based properties of solutions for model (1.2), we introduce the following lemmas.

Lemma 3.1. For any initial value $(S(0), V(0), I(0))\in\mathbb{R}_{+}^{3}$, model (1.2) has a unique solution $(S(t), V(t), I(t))$ defined for all $t\geq0$, and this solution remains in $\mathbb{R}_{+}^{3}$ with probability one.

The proof of this lemma is similar to that given in ^[37], we here omit it.

Lemma 3.2. For any initial value $(S(0), V(0), I(0)) \in R_{+}^{3}$, then solution $(S(t), V(t), I(t))$ of model (1.2) satisfies

Moreover, if $\mu > \frac{1}{2}\left(\sigma_{1}^{2} \vee \sigma_{2}^{2} \vee \sigma_{3}^{2}\right)$, we obtain

where $x_1(t) = S(t)$, $x_2(t) = V(t)$ and $x_3(t) = I(t)$.

The proof of Lemma 3.2 is similar to that given in ^[37], we hence omit it here.

In this section, we mainly discuss the extinction and persistence in mean of the disease in model (1.2). Firstly, define a threshold value as follows:

where $\mu_3 = \mu+\frac{\sigma_{3}^{2}}{2}+\int_{Z}(\eta_{3}(u)-\ln[1+\eta_{3}(u)])v(du)$.

For the following stochastic system:

with the initial values $\bar{S}(0) = S(0) > 0$ and $\bar{V}(0) = V(0) > 0$, owing to the stochastic comparison theorem, we have

where $(S(t), V(t), I(t))$ is the solution of model (1.2) with initial value $(S(0), V(0), I(0))$. Moreover, by integrating (3.4) we easily obtain that

Now, we can establish the following main conclusions in this section.

Theorem 3.1. Let $(S(t), V(t), I(t))$ be any positive solution of system (1.2) with initial value $(S(0), V(0), I(0))\in R_{+}^{3}$. Then the following conclusions hold.

$(i)$ If $R_{0}^{s} < 1$, then $\limsup _{t \rightarrow \infty} \frac{\ln I(t)}{t} \leq(\mu_{3}+\delta)\left(R_{0}^{s}-1\right) < 0\; a.s.$ That is, $\lim_{t \rightarrow \infty} I(t) = 0\; a.s.$ Furthermore, $\lim_{t\rightarrow \infty}\frac{1}{t}\int_{0}^{t}S(s)ds = \frac{(1-\pi)\Lambda}{\mu}$ and $\lim_{t\rightarrow \infty}\frac{1}{t}\int_{0}^{t}V(s)ds = \frac{\pi\Lambda}{\mu}\; a.s.$

$(ii)$ If $R_{0}^{s} > 1$, then $\liminf_{t \rightarrow \infty} \frac{1}{t} \int_{0}^{t} I(s)ds \geq \frac{1}{\gamma}(\mu_{3}+\delta)(R_{0}^{s}-1)\; a.s.$, where $\gamma: = \frac{\mu+\delta}{\mu}\left[\{\mu+\beta_{b}\}f'(0)+\{\mu+\beta_{v}\}g'(0)\right] > 0$. Furthermore, when $0\leq\pi < 1$, then $\liminf_{t\to\infty}\frac{1}{t} \int_{0}^{t} S(s)ds > 0$, and when $0 < \pi\leq 1$, then $\liminf_{t\to\infty}\frac{1}{t} \int_{0}^{t}V(s)ds > 0$.

Proof. $(i)$ Using Ito formula to $\ln I(t)$, we obtain

Hence, we have

In the light of the strong law of large numbers ^[42], we have $\lim_{t\rightarrow \infty}\frac{1}{t}\int_{0}^{t}dB_{3}(s) = 0\; a.s.$ and $\lim_{t\rightarrow \infty}\frac{1}{t}\int_{0}^{t}\int_{Z}[\ln(1+\eta_{3})I(s)-\ln I(s)]\widetilde{W}(ds, du) = 0\; a.s.$ Divided by $t$ the both sides of (3.7), then taking the superior limit and combining (3.6), we can obtain

It shows that $\lim_{t\rightarrow \infty}I(t) = 0\; a.s.$ Namely, the disease will be eliminated in the future. Furthermore, when $\lim_{t\rightarrow \infty}I(t) = 0\; a.s.$, we further obtain the limit system as follows:

Clearly, by integrating, as in (3.6), we can directly obtain

$(ii)$ Define function $U_{1}$ as follows:

Using Ito formula, we can obtain

Therefore,

Integrating both sides of (3.8) on interval $[0, t]$, and then divide by $t$, in the light of Lemma 3.2, we can obtain $\liminf_{t\rightarrow \infty}\frac{1}{t}\int_{0}^{t}I(s)ds\geq\frac{1}{\gamma}(\mu_{3}+\delta)(R_{0}^{s}-1) > 0\; a.s.$ It shows that if $R_{0}^{s} > 1$, then the disease will persistence in the mean.

Let $N(t) = S(t)+V(t)+I(t)$, then from model (1.2) and Lemma 3.1 we obtain

where $H(t) = \sum_{i = 1}^3[\sigma_ix_i(t)dB_i(t)+\int_{Z}\eta_i(u)x_i(t-)\widetilde{W}(dt, du)]$, here we denote $x_1(t) = S(t)$, $x_2(t) = V(t)$ and $x_3(t) = I(t)$ for the convenience. By the comparison principle of stochastic differential equations, we further obtain

where

Define $X(t) = N(0)+A(t)-U(t)+E(t)$, where $A(t) = \frac{\Lambda}{\mu}(1-e^{-\mu t})$ and $U(t) = N(0)(1-e^{-\mu t})$. From (3.9) we have $N(t)\leq X(t)$ for all $t\geq 0$. Obviously, $X(t)\geq 0\; a.s.$ for all $t\geq 0$, and $A(t)$ and $U(t)$ are continuous adapted increasing processes on $t\geq 0$ with $A(0) = U(0) = 0$. Therefore, by Theorem 3.9 in ^[42], we can obtain that $\lim_{t\to\infty}X(t) < \infty\; a.s.$ Consequently, $\limsup_{t\to\infty}N(t) < \infty\; a.s.$ Thus, there is a constant $M_0 > 0$, which is dependent on the solution $(S(t), V(t), I(t))$, such that for all $t\geq 0$,

When $0\leq\pi < 1$, from the first equation of model (1.2) and (3.10), we have

From Lemma 3.2 we finally obtain

When $0 < \pi\leq 1$, a similar argument as in above, we can obtain

This completes the proof.

Remark 3.1. When $\sigma_i = 0$ and $\eta_i(u)\equiv 0\; (i = 1, 2, 3)$, we see that model (1.2) becomes into model (1.1), and threshold value $R_0^s$ reduces to the basic reproduction number $R_0$ of model (1.1). Therefore, a comparison of Theorem 3.1 with Theorem 2.1 and Remark 2.1 shows that the conclusions on the extinction and persistence of the disease for model (1.1) is extended to the conclusions on the extinction and persistence in the mean with probability one for model (1.2).

Remark 3.2. Regrettably, in conclusion $(ii)$ of Theorem 3.1 we do not obtain a more strong conclusion for $S(t)$ and $V(t)$ just as for $I(t)$. That is, there is a common constant $m > 0$ such that for any positive solution $(S(t), V(t), I(t))$ of model (1.2) one has $\liminf_{t\rightarrow \infty}\frac{1}{t}\int_{0}^{t}S(s)ds\geq m\; a.s.$ or $\liminf_{t\rightarrow \infty}\frac{1}{t}\int_{0}^{t}V(s)ds\geq m\; a.s.$ Here, we will leave this problem in the future study.

4.   Stationary distribution

In this section, we discuss the ergodicity and the existence of stationary distribution in model (1.2). We can directly establish the following result.

Theorem 4.1. Assume $R_{0}^{s} > 1$. Then any positive solution $(S(t), V(t), I(t))$ of model (1.2) with initial value $(S(0), V(0), I(0))\in R_{+}^{3}$ is ergodic and has a unique stationary distribution $\pi(\cdot)$.

Proof. Define function $W: \mathbb{R}_{+}^{3}\rightarrow \mathbb{R}_+$ as follows:

where $U_{1}$ is defined in the proof of conclusion $(ii)$ of Theorem 3.1, $Q > 0$ is an enough large constant which is determined below, $\theta\in (0, 1)$ is an enough small constant satisfying $\rho = \mu-\frac{\theta}{2}(\sigma_{1}^{2}\vee\sigma_{2}^{2}\vee\sigma_{3}^{2}) > 0.$ Additionally, $U_{2} = -\ln S$, $U_{3} = -\ln V$ and $U_{4} = \frac{1}{\theta+1}(S+V+I)^{\theta+1}$. For any integer $k > 0$, let set $H_k = (\frac{1}{k}, k)\times(\frac{1}{k}, k)\times(\frac{1}{k}, k)\subset R_+^3$. It is clear that

Hence, function $W(S, V, I)$ has a minimum point $(S_{0}, V_{0}, I_{0})$ in the interior of $\mathbb{R}_{+}^{3}$. Then, we can define the nonnegative function $U: \mathbb{R}_{+}^{3}\rightarrow \mathbb{R}_{+}$ as follows:

In view of Ito formula, we have

From the proof of conclusion $(ii)$ of Theorem 3.1, we have obtained

Calculating $\mathcal{L}U_2$, $\mathcal{L}U_3$ and $\mathcal{L}U_4$ respectively, we further obtain

and

where

Choose constant $Q > 0$ satisfying the following inequality:

Then, from (4.1)–(4.6) we can obtain

Now, we define the set $H = \{(S, V, I): \xi < S < \frac{1}{\xi}, \xi < V < \frac{1}{\xi}, \xi < I < \frac{1}{\xi}\},$ where $\xi > 0$ is an enough small constant satisfying the following inequalities:

where $B = \sup_{I\in \mathbb{R}_{+}^{1}}\{[Q\gamma+\beta_{b}f'(0)+\beta_{v}g'(0)] I-\frac{\rho}{4}I^{\theta+1}\} < \infty.$

Divide the set $H^c$ into the following six domains:

Now, we deduce $\mathcal{L}U(S, V, I)\leq-1$ for any $(S, V, I)\in H^{c}$ in the following four cases.

Case 1. If $(S, V, I)\in H_1$, by (4.8), we can get

Case 2. If $(S, V, I)\in H_2$, by (4.9), we can get

Case 3. If $(S, V, I)\in H_3$, by (4.10), we can get

Case 4. If $(S, V, I)\in H_4\cup H_5\cup H_6$, by (4.11), we can get

In addition, for model (1.2) the diffusion matrix is $A = diag(\sigma_{1}^{2}S^{2}, \sigma_{2}^{2}V^{2}, \sigma_{3}^{2}I^{2})$. Choose $M_{0} = \min_{(S, V, I)\in\overline{H}}\{\sigma_{1}^{2}S^{2}, \sigma_{2}^{2}V^{2}, \sigma_{3}^{2}I^{2}\},$ we can obtain

for any $(S, V, I)\in\overline{H}$ and $\tau = (\tau_{1}, \tau_{2}, \tau_{3})\in\mathbb{R}_{+}^{3}$.

Thus, by Rayleigh's principle in ^[43] and ^[44], the conditions $(i)$ and $(ii)$ in [45,Lemma 4.4] are verified, respectively. It follows from the conclusions in [45,Lemma 4.4] that model (1.2) is ergodic and has a unique stationary distribution $\pi(\cdot)$. This completes the proof.

Remark 4.1. From Theorems 3.1 and 4.1, we see that when threshold value $R_0^s > 1$, for model (1.2), we not only establish the persistence of positive solutions in the mean, but also the ergodicity and the existence of stationary distribution of positive solutions. This shows that model (1.2) has more rich dynamical properties than the corresponding deterministic model (1.1).

Remark 4.2. From conclusion $(i)$ in Theorem 3.1, we can obtain that if any positive solution $(S(t), V(t), I(t))$ of model (1.2) is persistent in the mean then threshold value $R_0^s\geq 1$. If we further can get $R_0^s > 1$, then from Theorem 4.1 it will show that the solution $(S(t), V(t), I(t))$ has a unique stationary distribution. This seems to indicate that the persistence in the mean for model (1.2) may imply the existence of stationary distribution. Here, we will leave this issue in the future study.

5.   Density function

In this section, we investigate the existence and calculation of log-normal probability density function for model (1.2). Firstly, in order to facilitate calculation and demonstration, we introduce the logarithmic transformation $x_1 = \ln S$, $x_2 = \ln V$ and $x_3 = \ln I$ to model (1.2), then by Ito formula, we have

where $\mu_{i} = \mu+\frac{\sigma_{i}^{2}}{2}+\int_{Z}(\eta_{i}(u)-\ln[1+\eta_{i}(u)])v(du), i = 1, 2, 3.$ Define a constant as follows:

If $\tilde{R}_{0}^{s} > 1$, then function equation:

has a unique positive root $I_{*}^{+}.$ Define $E_{*}^{+} = (S_{*}^{+}, V_{*}^{+}, I_{*}^{+}) = (e^{x_{1}^{*}}, e^{x_{2}^{*}}, e^{x_{3}^{*}})$, where

It can be seen that $E_{*}^{+} = (S_{*}^{+}, V_{*}^{+}, I_{*}^{+})$ coincides with endemic equilibrium $P^* = (S^*, V^*, I^*)$ of model (1.1) when $\sigma_{1} = \sigma_{2} = \sigma_{3} = 0.$ In the general, $E_*^+$ is called the quasi-stationary state of model (1.2).

Let $(y_{1}, y_{2}, y_{3}) = (x_{1}-x_{1}^{*}, x_{2}-x_{2}^{*}, x_{3}-x_{3}^{*})$, where $x_{1}^{*} = \ln S_{*}^{+}$, $x_{2}^{*} = \ln V_{*}^{+}$ and $x_{3}^{*} = \ln I_{*}^{+},$ the linearization of system (5.1) at $E_*^+$ is

where

To prove the existence of a log-normal probability density function for model (1.2), we firstly introduce the following auxiliary lemmas.

Lemma 5.1. ^[31] For the algebraic equation $G_{0}^{2}+A_{0}\Sigma_{0}+\Sigma_{0}A_{0}^{T} = 0$, where $G_{0} = diag(1, 0, 0)$,

If $l_1 > 0$, $l_3 > 0$ and $l_1l_2 -l_3 > 0$, then $\Sigma_{0}$ is a positive definite matrix and has the expression:

Lemma 5.2. ^[31] For the algebraic equation $G_{0}^{2}+A_{0}\Gamma_{0}+\Gamma_{0}A_{0}^{T} = 0$, where $G_{0} = diag(1, 0, 0)$,

If $b_1 > 0$ and $b_2 > 0$, then $\Gamma_{0}$ is a positive semidefinite matrix and has the expression:

Next, by introducing a new calculation technique for the density function, the conclusion on the existence of log-normal probability density function is given as follows.

Theorem 5.1. Let $(y_{1}, y_{2}, y_{3})$ be any solution of system (5.2) with initial value $(y_{1}(0), y_{2}(0), y_{3}(0)) \in R^3$, If $\tilde{R}_{0}^{s} > 1$, then there is a log-normal probability density function $\Phi(y_{1}, y_{2}, y_{3})$ around quasi-stationary state $E_{*}^{+}$, which has the following expression:

where $\Sigma$ is a positive definite matrix, and the specific form of $\Sigma$ is given as follows:

where

and if $l_{11}\neq l_{22}$, then

if $l_{11} = l_{22}$, then

where $\Delta = \frac{l_{13}l_{23}(l_{11}-l_{22})}{l_{11}l_{22}}$ and $\rho_{i}^{2} = \sigma_{i}^{2}+2\int_{Z}(\eta_{i}(u)-\ln[1+\eta_{i}(u)])v(du)\; (i = 1, 2, 3).$

Proof. In view of (5.2), let $Y = (y_{1}, y_{2}, y_{3})^{T}$ and the coefficients matrix

Similar to the method in ^[46], the density function $\Phi(Y) = \Phi(y_{1}, y_{2}, y_{3})$ of the quasi-stationary distribution of system (5.1) around the origin point can be obtained by solving the following three-dimensional Fokker-Plank equation

its form may be expressed approximately as a Gaussian distribution

where $Y^* = (0, 0, 0)$, and $\tilde{Q}$ is a real symmetric matrix which satisfies the following equation:

where $G^2 = diag(\rho_{1}^{2}, \rho_{2}^{2}, \rho_{3}^{2})$ with $\rho_{i}^{2} = \sigma_{i}^{2}+2\int_{Z}(\eta_{i}(u)-\ln[1+\eta_{i}(u)])v(du), i = 1, 2, 3.$

If $\tilde{Q}$ is positive definite matrix, let $\tilde{Q}^{-1} = \Sigma$, then

Therefore, if a positive definite matrix $\Sigma$ is calculated, then positive definite matrix $\tilde{Q}$ will be obtained. Thus, density function $\Phi(Y)$ will be concretely acquired. According to ^[47], Eq (5.8) can be formed from the sum of the following three equations:

where $\Sigma = \Sigma_{1}+\Sigma_{2}+\Sigma_{3}$ and $G^{2} = G_{1}^{2}+G_{2}^{2}+G_{3}^{2}$ with

Next, it will be shown that $A$ is a stable matrix. In fact, the characteristic equation of the matrix $A$ is

where $l_{1} = l_{11}+l_{22}+l_{33} > 0$, $l_{2} = l_{11}(l_{22}+l_{33})+l_{22}l_{33}+l_{23}l_{32}+l_{13}l_{31} > 0$ and $l_{3} = l_{11}(l_{22}l_{33}+l_{23}l_{32})+l_{13}l_{22}l_{31} > 0$. By calculation, we can obtain

which means that $A$ is a stable matrix by the Hurwitz criterion.

Now, the special expression of $\Sigma$ can be found in three steps as follows: $\Sigma = \Sigma_{1}+\Sigma_{2}+\Sigma_{3}$.

Step 1. We consider the following algebraic equation:

We will choose a reversible matrix $J_{1}$ with the expression

such that Eq (5.11) changes to the following form:

which satisfies $J_1G_1^2J_1^T = G_0^2 = diag(1, 0, 0)$ and $J_1AJ_1^{-1} = A_0.$ Let $\Sigma_{0} = J_{1}\Sigma_{1}J_{1}^{T}$, then Eq (5.12) is equivalently rewritten by $G_0^2+A_0\Sigma_0+\Sigma_0A_0^{T} = 0$.

According to $G_{0}^{2} = J_{1}G_{1}^{2}J_{1}^{T}$, we have

it implies that

In view of $A_{0} = J_{1}AJ_{1}^{-1}$ and (1.2), namely, $J_{1}A = A_{0}J_{1}$, we have

Hence, we have $l_{31}\mu_{33} = 0$, $-l_{22}\mu_{32}+l_{32}\mu_{33} = \mu_{22}$, $l_{31}\mu_{23} = \mu_{11}$, $-l_{23}\mu_{32}-l_{33}\mu_{33} = \mu_{23}$, $-l_{22}\mu_{22}+l_{32}\mu_{23} = \mu_{12}$ and $-l_{23}\mu_{22}-l_{33}\mu_{23} = \mu_{13}$. Solving these equations we further can obtain

In addition, by carefully calculating we can verify that $-l_{11}\mu_{11}+l_{31}\mu_{13} = -l_1\mu_{11}$, $-l_{22}\mu_{12}+l_{32}\mu_{13} = -(l_1\mu_{12}+l_2\mu_{22}+l_3\mu_{32})$ and $-(l_{13}\mu_{11}+l_{23}\mu_{12}+l_{33}\mu_{13}) = -(l_1\mu_{13}+l_2\mu_{23}).$ Thus, by (5.13) and (5.14), the specific expression of $J_1$ is calculated as

Clearly, $J_1$ is reversible. From Lemma 5.1, we have known that $\Sigma_{0}$ is positive definite and

Therefore, from $\Sigma_{0} = J_{1}\Sigma_{1}J_{1}^{T}$, we finally obtain a positive definite matrix

Step 2. We consider the following algebraic equation:

Similarly to Step 1, we will choose a reversible matrix $J_{2}$ with the expression

such that Eq (5.16) changes to the following form:

which satisfies $J_2G_2^2J_2^T = G_0^2 = diag(1, 0, 0)$ and $J_2AJ_2^{-1} = A_0.$ According to $G_{0}^{2} = J_{2}G_{2}^{2}J_{2}^{T}$, we have

which implies

In view of $A_{0} = J_{2}AJ_{2}^{-1}$ and (5.17), namely, $J_{2}A = A_{0}J_{2}$, we have

Hence, we obtain $l_{32}m_{33} = 0$, $l_{32}m_{23} = m_{12}$, $-l_{13}m_{31}-l_{33}m_{33} = m_{23}$, $-l_{11}m_{31}+l_{31}m_{33} = m_{21}$, $-l_{11}m_{21}+l_{31}m_{23} = m_{11}$, $-l_{13}m_{21}-l_{33}m_{23} = m_{13}$. Solving these equations, we further obtain

In addition, by carefully calculating we can verify that $-l_{11}m_{11}+l_{31}m_{13} = -(l_{1}m_{11}+l_{2}m_{21}+l_{3}m_{31})$, $-l_{22}m_{12}+l_{32}m_{13} = -l_{1}m_{12}$ and $-(l_{13}m_{11}+l_{23}m_{12}+l_{33}m_{13}) = -(l_{1}m_{13}+l_{2}m_{23})$. Thus, by (5.17) and (5.18), we can also obtain

Clearly, $J_2$ is reversible. From Lemma 5.1, we finally obtain a positive definite matrix

Step 3. We consider the algebraic equation

Likewise, we will choose a reversible matrix $J_{3}$ such that Eq (5.20) changes to the following form:

which satisfies $G_{0}^{2} = J_{3}G_{3}^{2}J_{3}^{T}$ and $A_0 = J_{3}AJ_{3}^{-1}$, where

In the light of $G_{0}^{2} = J_{3}G_{3}^{2}J_{3}^{T}$, we have

Then, we can obtain

In view of $A_{0} = J_{3}AJ_{3}^{-1}$ and (5.21), namely, $J_{3}A = A_{0}J_{33}$, we have

Thus, we have

which implies

Assume $l_{11}\neq l_{22}$. By solving Eq (5.23), we have

where $\Delta = \frac{l_{13}l_{23}(l_{11}-l_{22})}{l_{11}l_{22}}$. According to (5.22) and (5.24), one can easily obtain $n_{11} = -l_{11}n_{21} = \frac{l_{11}l_{23}}{l_{22}\rho_{3}\Delta}$, $n_{12} = -l_{22}n_{22} = -\frac{l_{13}l_{22}}{l_{11}\rho_{3}\Delta}$, $n_{31} = -\frac{n_{21}}{l_{11}} = \frac{l_{23}}{l_{11}l_{22}\rho_{3}\Delta}$ and $n_{32} = -\frac{n_{22}}{l_{22}} = -\frac{l_{13}}{l_{11}l_{22}\rho_{3}\Delta}$. In addition, by carefully calculating we can verify that $-l_{11}n_{11}+l_{31}n_{13} = -(l_{1}n_{11}+l_{2}n_{21}+l_{3}n_{31})$, $-l_{22}n_{12}+l_{32}n_{13} = -(l_{1}n_{12}+l_{2}n_{22}+l_{3}n_{32})$ and $-(l_{13}n_{11}+l_{23}n_{12}+l_{33}n_{13}) = -l_{1}n_{13}$. Thus, we can obtain

Clearly, $J_3$ is reversible. From Lemma 5.1, we finally obtain a positive definite matrix

When $l_{11} = l_{22}$, we will use Lemma 5.2. Let

where $b_{1} = l_{22}+l_{33}$, $b_{2} = l_{22}l_{33}+l_{23}l_{32}+l_{13}l_{31}$ and $b_3$ is given below. From $A_{0} = J_{3}AJ_{3}^{-1}$ and (5.21), namely, $J_{3}A = A_{0}J_{33}$, we have

Thus, we have

We further have

and then

Choose $n_{32} = 1$ and $b_3 = -l_{32}n_{13} = -\frac{l_{32}}{\rho_3}$, then from (5.26) and (5.27) we easily obtain $n_{12} = 0$, $n_{31} = -\frac{l_{23}}{l_{13}}$, $n_{11} = \frac{l_{11}}{l_{13}\rho_{3}}$, $n_{21} = -\frac{1}{l_{13}\rho_{3}}$ and $n_{22} = 0$. Thus, we finally have

Clearly, $J_3$ is reversible. From Lemma 5.2, we can choose a semipositive definite $\Gamma_{0}$ as follows:

By $\Gamma_{0} = J_{3}\Sigma_{3}J_{3}^{T}$, we finally obtain a semipositive definite matrix

Summarizing the above calculations, we finally conclude that there exists a real symmetric positive definite matrix $\Sigma = \Sigma_{1}+\Sigma_{2}+\Sigma_{3}$ satisfying (5.8). As a result, there is a locally approximate log-normal probability density function

near the quasi-stationary state $E_{*}^{+}$. This completes the proof.

Remark 5.1. From the proof of Theorem 5.1 it is shown that a new calculation technique for matrix $\Sigma$ is proposed. Obviously, this technique is different from the calculation method given in ^[31].

6.   Numerical examples

In this section, we present the simulation results to give the reader a clear understanding of our results were achieved by using the method mentioned in ^[48]. Throughout the following numerical simulations, we choose the nonlinear incidence functions as follows:

Example 1. In model (1.2), we choose the parameters $\mu = 0.5$, $\lambda = 2.5$, $\beta_{b} = 0.4$, $\beta_{\nu} = 0.2$, $p = 0.4$, $\delta = 0.85$, $H_{b} = H_{\nu} = 1$, $\alpha = 1$, $(\sigma_{1}, \sigma_{2}, \sigma_{3}) = (0.2, 0.2, 0.65)$, $(\eta_{1}, \eta_{2}, \eta_{3}) = (0.01, 0.01, 0.02)$ and $\nu(Z) = 1$. By calculating, from (3.3) we obtain $R_{0}^{s} = 0.98 < 1$, which means that disease $I(t)$ will disappear with probability one by conclusion $(i)$ of Theorem 3.1. However, model (1.1) has an endemic equilibrium $P^* = (S^*, V^*, I^*)$, which is local asymptotically stable because the basic reproduction number $R_{0} = 1.0963 > 1$ by Theorem 2.1.

The numerical simulations are presented in Figure 1 in allusion to the deterministic, white noise and Lévy jumps, respectively. We easily see from Figure 1 that the solution $(S(t), V(t), I(t))$ of deterministic model (1.1) converges to its endemic equilibrium as $t\to\infty$, and the solution $(S(t), V(t), I(t))$ for stochastic model (1.2) satisfies that $I(t)$ is extinct with probability one, and $S(t)$, $V(t)$ are persistence in the mean. Therefore, conclusion $(ii)$ in Theorem 2.1 and conclusion $(i)$ in Theorem 3.1 are verified by the numerical simulations. This also demonstrates that the jump noise have a positive impact on control the diseases. Hence, the impact of the noise cannot be overlooked in modeling process.

Example 2. In model (1.2), we take the parameters $(\sigma_{1}, \sigma_{2}, \sigma_{3}) = (0.02, 0.02, 0.06)$ and other parameters are given as in Example 1. By calculation, we obtain $R_{0}^{s} = 2.949 > 1$, which shows that the diseases will persistence in the mean and any positive solution of model (1.2) is ergodic and has a unique stationary distribution by Theorems 3.1 and 4.1, respectively.

The numerical simulations are presented in Figure 2. We easily see from Figure 2 that the solution $(S(t), V(t), I(t))$ of model (1.1) converges to its endemic equilibrium as $t\to\infty$, and the solution $(S(t), V(t), I(t))$ for stochastic model (1.2) is persistence in the mean and has a unique stationary distribution. Therefore, conclusion $(ii)$ in Theorem 2.1, conclusion $(ii)$ in Theorem 3.1 and Theorem 4.1 are verified by the numerical simulations.

In addition, by calculating we also have $\tilde{R}_{0}^{s} = 2.9478 > 1$. Hence, there is a log-normal probability density function $\Phi(y_{1}, y_{2}, y_{3})$ in the quasi-stationary state $E_{*}^{+}$ by Theorem 5.1. Moreover, it is calculated that $\Delta = \frac{l_{13}l_{23}(l_{11}-l_{22})}{l_{11}l_{22}} = 0.00122\neq0$, which implies

By simple calculation, one can get that $E_{+}^{*} = (S_{+}^{*}, V_{+}^{*}, I_{+}^{*}) = (5.074, 3.952, 1.154)$. Thus, the log-normal probability density function $\Phi(S, V, I)$ of system (1.2) is derived as

The numerical simulations are given in Figure 3, which present the visual expressions of marginal density functions of solution $(S(t), V(t), I(t))$ for model (1.2).

7.   Conclusions

In this paper, we investigated the dynamical behavior for a stochastic SVI epidemic model with white noise, Lévy jumps and nonlinear incidence. In order to observe the influence of randomness, deterministic model (1.1) is first discussed. The basic reproduction number $R_{0}$ is calculated, and then it is proved the disease-free equilibrium is globally asymptotically stable if $R_{0} > 1$, otherwise, the endemic equilibrium is local asymptotically stable if $R_{0} > 1$. For stochastic model (1.2), a new threshold value $R_0^s$ is defined, and when $R_0^s < 1$ then the extinction with probability one of the disease is proved, and when $R_0^s > 1$ then the persistence in the mean and the existence of stationary distribution for any positive solution are established. Furthermore, we also established the existence criterion for the log-normal probability density function by solving the corresponding Fokker-Planck equation. Particularly, a new technique for the calculation of probability density function is introduced. The results show that Lévy noise can effectively alter the dynamical behaviour of the disease and that it can also contribute to its extinction.

Theorems 3.1 and 4.1 imply that the disease dies out as threshold value $R_{0}^{s} < 1$, and otherwise the disease persists and possesses a unique stationary distribution as $R_{0}^{s} > 1$. This shows that $R_{0}^{s}$ plays a similar role to the basic reproduction number $R_{0}$ of model (1.1). Comparing $R_{0}$ and $R_{0}^{s}$ given in (2.1) and (3.3), we have $R_{0}^{s} < R_{0}$. This means that Lévy jumps can also inhibit the outbreak of the disease. These important results demonstrate that the Lévy jumps process may have a greater impact on the dynamical properties of model (1.2).

There is a few interesting topics to worthy of further study. It is possible to come up with more realistic and complex stochastic models, such as considering the impact of vaccination of susceptible individuals, vaccine effectiveness on model (1.2). In addition, it is important to note that the approach used in this paper can also be applied to the study of other interesting models, such as COVID-19 spread model, SVEIS model, SVIRS model and so on. We will study these problems in the future.

Furthermore, in ^[39,40] we see that the stochastic SIC epidemic system with quadratic white noise and Lévy jumps and an application to COVID-19 in Morocco, and the stochastic and fractal-fractional Atangana-Baleanu order hepatitis B model with Lévy noise are proposed and investigated. Particularly, in ^[40] the probability density function is discussed for quadratic white noise intensity by the numerical simulations. Therefore, an interesting and challenging issue is to establish the theoretical results in allusion to probability density function {for the} above two kinds of models.

Acknowledgments

This research was supported by Program for Tianshan Innovative Research Team of Xinjiang Uygur Autonomous Region, China (2020D14020) and National Natural Science Foundation of China (Grant No. 72174175, 72064036) and Chinese Foundation for Hepatitis Prevention and Control(YGFK20200059)

Conflict of interest

We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with this work.