\documentclass{mc} %%%%% journal info -- DO NOT CHANGE %%%%%%%%% \setcounter{page}{1} \renewcommand\thisnumber{x} \renewcommand\thisyear {2025} \renewcommand\thismonth{xxx} \renewcommand\thisvolume{30} \renewcommand\datereceived{July 14, 2023} \renewcommand\dateaccepted{July 6, 2024} %%%%%%%% end journal info %%%%%%%%%%%%% %Rad zaprimljen dana: 14.7.2023. %Rad prihvacen dana: 6.7.2024. %%%%% author macros %%%%%%%%% % place your own macros HERE %%%%% end %%%%%%%%% %%%%%%%%% already loaded packages -- DO NOT RELOAD %%%%% % \usepackage{amssymb} % \usepackage{amsmath,amsthm} % \usepackage{latexsym} % \usepackage{amsfonts} % \usepackage{amsbsyt} % \usepackage{graphicx} %%%%%% end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{document} %\begin{linenumbers} %%%%% title and author(s): % \markboth{Author(s)}{Short Title} % \title{Title} \markboth{B.\,Jovanovi\'c}{Stochastic delayed epidemic model with L\'evy jump} \title[Stochastic delayed epidemic model with L\'evy jump]{Dynamical analysis and stationary distribution of a stochastic delayed epidemic model with L\'evy jump\thanks{This work was supported by the Ministry of Science, Technological Development and Innovation of the Republic of Serbia, agreement no. 451-03-66/2024-03/200124.}} % single author: % \author[AUTHOR]{AUTHOR} % \address{address of AUTHOR} % \email{{\tt email address of AUTHOR} (AUTHOR)} \author[B.\,Jovanovi\'c]{Bojana Jovanovi\'c\corrauth} \address{Faculty of Sciences and Mathematics, University of Ni\v s, Vi\v{s}egradska 33, 18\,000 Ni\v{s}, Serbia} \email{{\tt bojana.jovanovic@pmf.edu.rs} (B.\,Jovanovi\'c)} % multiple authors: % Please mark \corrauth after the name of the corresponding author. % different addresses: %\author[AUTHOR1 and AUTHOR2]{AUTHOR1\affil{1}\comma\corrauth and AUTHOR2\affil{2}} %\address{\affilnum{1}\ address of AUTHOR1\\ \affilnum{2}\ address of AUTHOR2} % %same address: %\author[AUTHOR1, AUTHOR2 and AUTHOR3]{AUTHOR1, AUTHOR2\corrauth and AUTHOR3} %\address{address of AUTHOR1, AUTHOR2 and AUTHOR3} % %\emails{{\tt email of AUTHOR1} (AUTHOR1), {\tt email of AUTHOR2} (AUTHOR2), {\tt email of AUTHOR3} (AUTHOR3)} % %%%%% Begin Abstract %%%%%%%%%%% \begin{abstract} In this paper, a stochastic delayed SVIR (susceptible-vaccinated-infected-\linebreak recovered) model with logistic growth of population, a saturated incidence function and distributed delay is analyzed. Sufficient conditions for the extinction and persistence in mean of the disease and existence of a stationary distribution are obtained. The theoretical results are illustrated via numerical simulations. \end{abstract} %%%%% end %%%%%%%%%%% %%%%% Keywords %%%%%%%%%%% \keywords{stochastic epidemic model, distributed delay, L\'evy jump, extinction, persistence in mean, stationary distribution} %%%% AMS subject classifications %%%% \ams{34F05, 60H10, 92D30} %%%% maketitle %%%%% \maketitle %%%% Start %%%%%% \section{Introduction} \label{intro} Infectious diseases, such as AIDS, Ebola, Zika, influenza, hepatitis B virus, COVID-19 and so on, are all public health concerns of the world. They have an important impact on public health, social and economic development and even human life (see \cite{Binder_1999,Blower_1991,Din_2021}). In order to better understand the transmission mechanism of infectious diseases, and obtain good strategies for preventing and controlling the spread of infectious diseases, mathematicians have proposed many mathematical models to discuss the dynamic behavior of infectious disease transmission (see, for example, \cite{Blower_1991,HHethcote_2000,HHethcote_Driessche_1991}). Various environmental variations in nature, such as people-to-environment interactions, person-to-person contacts and vaccination efficiency, are inherently subject to a continuous spectrum of disturbances (see \cite{Aksendal_2010,Allen_2008}, and references therein). In \cite{Mao_Marion_2002} and \cite{Li_Mao_2009} Mao et al. further demonstrated that white noise, which reflects the stochasticity, has a large destabilizing effect on epidemic transmission and species evolutions, suggesting that stochastic models can provide an additional degree of realism in comparison with their deterministic counterparts. There are some sudden environmental random perturbations which can produce sudden reduction in population size, for example, epidemics, earthquakes, and hurricanes. These sudden environmental random perturbations cannot be modeled by white noises. Some authors have suggested that one can use a Lévy jump process to model these sudden environmental random perturbations. Incorporation of L\'evy noise in the epidemic model is very significant in the investigation of model dynamics (see, for instance, \cite{Liu,Zhang1,Fatini17,Zhang16}). In many cases the constant recruitment rate is not a rational assumption for describing the epidemic model with considerable variability in the population size. For that reason, when describing dynamics of the population with fast mobility, for example due to fast transportation by metros, trains and airplanes within a period of time, authors assumed that the density of the susceptible in the municipal cities on holidays or weekends is described by logistic growth (see \cite{Li_Mao_2022,Liu_Jiang_2017,Cao_Feng_2019}). Additionally, the incidence rate of infectious diseases plays an important role in the investigation of epidemic models. A common incidence rate considered in the literature is the proportionate mixing standard incidence rate $\beta \frac{SI}{N}$ (see \cite{Cooke_Driessche_1996,WWang_2002}), which is based on the homogeneous mixing assumption, where it is assumed that a susceptible person (from the $S$ class) becomes infected when interacting with an infected person (from the $I$ class), $N$ is the total population and the coefficient $\beta >0$ denotes the transmission coefficient. There are a number of reasons why this incidence rate may require modification. The first is that the underlying assumption of homogeneous mixing may be invalid. However, incidence rates that increase more gradually than linearly in I and S may also result from saturation effects. Saturation incidence rates $\frac{\beta SI}{1+\alpha S}$ (see \cite{May_Anderson_1978}) or $\frac{\beta SI^l}{1+\alpha I^l}$ (see \cite{Ruan_Wang_2003}), where $\alpha$ is a saturation constant, are important because the number of effective contacts between infective and susceptible individuals may saturate due to the protection measures taken by susceptible individuals or due to awareness among individuals oft the infection. In addition, actual incidence rates are probably not strictly linear in each variable over the entire range of I and S, and that is described by $l$ in \cite{Ruan_Wang_2003}. Recently, considerable attention has been paid to studying the dynamics of epidemic models with time delays. Constant time delay was investigated in \cite{Miljana_Marija} and \cite{Marija2018}, while distributed time delay was considered by Ma et al. in \cite{Ma_Hara_2002}. It has been presented that time delays play a crucial role in dynamical properties of epidemic systems. Motivated by the above discussion, the main aim of this paper is to investigate the long-time behavior of a stochastic delayed SVIR model with logistic growth of population, a saturated incidence function and distributed delay. The delay, presented in this paper, reflects the time that it takes for an infection to develop and a person to become infectious after contact. To the best of the author's knowledge, this paper is the first attempt to study dynamic behavior and the existence of ergodic stationary distribution to a stochastic epidemic model with logistic growth, the saturated incidence rate and distributed delay perturbed by both white and L\'evy noise. The outline of this paper is as follows. In Section 2, a deterministic and stochastic SVIR epidemic model with logistic growth of population and the saturated incidence rate are described. On the basis of epidemic models described in Section 2, a stochastic epidemic model that contains L\'evy jumps and distributed time delay is constructed in Section 3. Preliminary results used in the paper are listed in Section 4. Sections 5 and 6 are devoted to giving sufficient conditions of the extinction and persistence in mean of the disease. In Section 7, a sufficient condition for the existence of an ergodic stationary distribution is derived. In Section 8, some numerical simulations are presented in order to demonstrate the main theoretical results. Section 9 contains concluding remarks. \section{Epidemic model } In \cite{Li_Mao_2022}, Li et al. proposed a stochastic SVIR epidemic system with saturation incidence rates and logistic growth. In this section, their result is briefly presented. In \cite{Li_Mao_2022}, Li et al. consider long-term properties of a stochastic epidemic model where human population is subdivided into four time-dependant classes, namely, susceptible $S $, vaccinated $V $, infected $I $ and recovered $R $. They assumed that the virus spreads when a susceptible person comes into contact with an infected person who is vaccinated or not. The dynamics of the population is given by the system \begin{align} \label{Li SVIR deterministic} \begin{split} dS(t)= &\left( \gamma S(t) \left( 1- \frac{S(t)}{K} \right) - \frac{\beta S(t) I(t)}{1+a_1 I(t)} - \zeta S(t) +v V(t) \right) dt , \\ dV(t)=& \left( \zeta S(t) - \frac{\beta V(t) I(t) }{1+a_2 I(t)} -v V(t) -\mu V(t) \right) dt , \\ dI(t)= &\left( \frac{\beta S(t) I(t)}{1+a_1 I(t)} + \frac{\beta V(t) I(t)}{1+a_2 I(t)} -\left( \tau+\mu+\delta \right) I(t) \right) dt , \\ dR(t)=& \left( \tau I(t) -\left( r+\mu \right) R(t) \right) dt. \end{split} \end{align} Here the intrinsic rate $\gamma$ equals the difference of the birth rate $b$ and the natural death rate $\mu$ ($\gamma = b- \mu >0$), and $K$ is the carrying capacity of a local population. Logistic growth of susceptible population $\gamma S(t) \left( 1- \frac{S(t)}{K} \right), \, t\geq 0,$ describes fast mobility of a local population that occurs due to fast transportation by metros, trains and airplanes within a period of time, as it is assumed that the density of the susceptible in the municipal cities on holidays or weekends is described by logistic growth. $\frac{\beta S(t) I(t)}{1+a_1 I(t)}$ and $\frac{\beta V(t) I(t)}{1+a_2 I(t)}$ are saturation incidence rates, where $\beta$ is the transmission rate measuring the average number of adequate contacts of a person per unit time and $a_1$ and $a_2$ denote half-saturation constants that satisfy the condition $a_1 < a_2$, that is, the probability that the vaccinated are infected by the infective individual is less than the probability that the susceptible are infected by the infective. In terms of these incidence rates, the overall growth of the infected population is less compared to a bilinear incidence rate due to awareness among individuals of the disease. The other parameters in model \eqref{Li SVIR deterministic} are described as follows: $\zeta$ is the proportion of the susceptibles who take the vaccine; $v$ is the rate of immunity loss of vaccines to the vaccinated; $\tau$ is the recovered rate of the infected; and $\delta$ means the mortality rate caused by infectious diseases to the infected. %measure the crowdedness of the infected when infectious diseases invade a local population In \cite{Li_Mao_2022}, based on the deterministic model \eqref{Li SVIR deterministic}, the corresponding stochastic model is constructed by assuming that the environmental noises are proportional to the variables $S, V, I , R$, where $B_1 (t), B_2 (t), B_3 (t)$ and $B_4(t)$ are four independent standard Brownian motions (or Wiener processes), while $\sigma_1, \sigma_2, \sigma_3$ and $\sigma_4$ are the intensities of white noises \begin{align} \label{Li SVIR stochastic} \begin{split} dS(t)= &\left( \!\gamma S(t) \left(\! 1- \frac{S(t)}{K} \!\right) - \frac{\beta S(t) I(t)}{1+a_1 I(t)} - \zeta S(t) +v V(t) \!\right) dt + \sigma_1 S(t) dB_1(t), \\ dV(t)=& \left( \zeta S(t) - \frac{\beta V(t) I(t) }{1+a_2 I(t)} -v V(t) -\mu V(t) \right) dt + \sigma_2 V(t) dB_2(t), \\ dI(t)= &\left( \frac{\beta S(t) I(t)}{1+a_1 I(t)} + \frac{\beta V(t) I(t)}{1+a_2 I(t)} -\left( \tau+\mu+\delta \right) I(t) \right) dt + \sigma_3 I(t) dB_3(t), \\ dR(t)= &\left( \tau I(t) -\left( r+\mu \right) R(t) \right) dt + \sigma_4 R(t) dB_4(t), \end{split} \end{align} respectively. In \cite{Li_Mao_2022}, the authors considered long-term properties of a stochastic SVIR model \eqref{Li SVIR stochastic} and established sufficient conditions for the extinction or persistence in mean of the disease. \section{Construction of the stochastic epidemic model with distributed delay and Lévy jumps} In this section, the deterministic epidemic model introduced in \cite{Li_Mao_2022} is modified in the following way. First, saturation incidence rates $\frac{\beta S(t) I(t)}{1+a_1 I(t)}$ and $\frac{\beta V(t) I(t)}{1+a_2 I(t)}$ describing the crowdedness of the infected are replaced with saturation incidence rates $\frac{\beta S(t) I(t)}{1+m_1 S(t)}$ and $\frac{\beta V(t) I(t)}{1+m_2 V(t)}$ describing appropriate preventions taken by susceptible or vaccinated individuals in order to control the spread of the epidemic. Second, following the approach presented in \cite{Caraballo_Fatini_2020}, a delay kernel $g (s)$ is introduced into the deterministic SVIR epidemic model. The constructed model is a system of four first order ordinary differential equations shown as below: \begin{align} \label{SVIR deterministic} \begin{split} \frac{dS(t)}{dt} =& \gamma S(t) \left(\! 1- \frac{S(t)}{K} \!\right) - \frac{\beta S(t)}{1+m_1 S(t)} \int_{- \infty}^{t}\!\! g(t- \tau) I( \tau) d\tau -v S(t) +v_1 V(t) , \\ \frac{dV(t)}{dt} = & v S(t) - \frac{\beta V(t) }{1+m_2 V(t)} \int_{- \infty}^{t} g(t- \tau) I( \tau) d\tau -v V(t) -\mu V(t) , \\ \frac{dI(t)}{dt} = & \frac{\beta S(t)}{1+m_1 S(t)}\int_{- \infty}^{t} g(t- \tau) I( \tau) d\tau + \frac{\beta V(t)}{1+m_2 V(t)} \int_{- \infty}^{t} g(t- \tau) I( \tau) d\tau \\ & -\left( \rho+\mu+\delta \right) I(t), \\ \frac{dR(t)}{dt} = & \rho I(t) -\left( r+\mu \right) R(t), \end{split} \end{align} where all parameters are positive numbers. The delay kernel $g : [0, \infty ) \rightarrow [0, \infty )$ is an $L^1$ -function, normalized so that $\int_0^{\infty} g(s) ds = 1$. The quality $\varrho = \int_0^{\infty} s g(s) ds $ is known as the average delay for the kernel. In biological systems, the delay kernel with Gamma distribution is often used $g(s) = \frac{s^n \xi^{n+1} e^{-\xi s}}{n!}, \,\,\, s \in \left( 0, \infty \right)$, where $\xi > 0$ is a parameter denoting exponentially fading memory also known as the retrogradation of the past memories effect, and $n$ is a nonnegative integer (see \cite{Macdonald_1978}). For this type of delay kernel, it can be proved that the average delay is $\varrho = \frac{n+1}{\xi}$. In the literature, the kernel with $n = 0$, i.e., $ g(s) = \xi e^{- \xi s}$ is called the weak kernel and it is considered in this paper. By substituting function $g(s)$, system \eqref{SVIR deterministic} becomes: \begin{align} \label{SVIR deterministic1} \begin{split} \frac{dS(t)}{dt} = & \gamma S(t) \left( \!1\!-\! \frac{S(t)}{K}\! \right) - \frac{\beta S(t)}{1+m_1 S(t)}\! \int_{- \infty}^{t}\!\! \xi e^{-\xi (t-\tau)} I( \tau) d\tau -v S(t) +v_1 V(t) , \\ \frac{dV(t)}{dt} = & v S(t) - \frac{\beta V(t) }{1+m_2 V(t)} \int_{- \infty}^{t} \xi e^{-\xi (t-\tau)} I( \tau) d\tau -v V(t) -\mu V(t) , \\ \frac{dI(t)}{dt} = & \frac{\beta S(t)}{1+m_1 S(t)}\int_{- \infty}^{t} \xi e^{-\xi (t-\tau)} I( \tau) d\tau + \frac{\beta V(t)}{1+m_2 V(t)} \int_{- \infty}^{t} \xi e^{-\xi (t-\tau)} I( \tau) d\tau \\ & -\left( \rho+\mu+\delta \right) I(t), \\ \frac{dR(t)}{dt} = & \rho I(t) -\left( r+\mu \right) R(t). \end{split} \end{align} Setting $X(t) = \int_{- \infty}^{t} \xi e^{-\xi (t-\tau)} I( \tau) d\tau,$ it follows from the linear chain technique that system \eqref{SVIR deterministic1} can be transformed into the following equivalent system: \begin{align} \label{SVIRX deterministic} \begin{split} dS(t)= &\left( \gamma S(t) \left( 1- \frac{S(t)}{K} \right) - \frac{\beta S(t) X(t)}{1+m_1 S(t)} -v S(t) +v_1 V(t) \right) dt , \\ dV(t)= &\left( v S(t) - \frac{\beta V(t) X(t)}{1+m_2 V(t)} -v_1 V(t) -\mu V(t) \right) dt , \\ dI(t)= &\left( \frac{\beta S(t) X(t)}{1+m_1 S(t)} + \frac{\beta V(t) X(t)}{1+m_2 V(t)} -\left( \rho+\mu+\delta \right) I(t) \right) dt , \\ dR(t)= &\left( \rho I(t) -\left( r+\mu \right) R(t) \right) dt , \\ dX(t)= &\xi \left( I(t) - X(t) \right) dt . \end{split} \end{align} Based on the deterministic epidemic model \eqref{SVIRX deterministic}, appropriate stochastic is constructed by assuming that environmental fluctuations are white noise and Lévy type, which are directly proportional to $S, V, I, R$ and $X$ class, in order to describe the continuous environmental variability and sudden and significant impacts that do not occur often. Throughout this paper, let $(\Omega, \mathcal{F}, \left\{ \mathcal{F}_{t} \right\}_{t\geq0},\textbf{P})$ be a complete probability space with a natural filtration $\left\{ \mathcal{F}_{t} \right\}_{t\geq0}$, which is right continuous and such that $\mathcal{F}$ contains all $\textbf{P}$-null sets. Let $B_i(t), i=1,2,3,4,5,$ be five independent standard Brownian motions and $\bar{N}$ a Poisson counting measure, independent of $B_i, i=1,2,3,4,5,$ and let them be defined on the given probability space $(\Omega, \mathcal{F}, \left\{ \mathcal{F}_{t} \right\}_{t\geq0},\textbf{P})$. $\tilde{N}(dt,du)= \bar{N}(dt,du) -\lambda(du)dt$ denotes the compensator of $\bar{N}$ and the characteristic measure $\lambda$ is defined on a measurable subset $Y$ of $(0, + \infty)$ satisfying $\lambda(Y)< + \infty.$ A stochastic epidemic model takes the following form: \begin{align*} \begin{split} dS(t) = & \left( \gamma S(t) \left( 1- \frac{S(t)}{K} \right) - \frac{\beta S(t) X(t)}{1+m_1 S(t)} -v S(t) +v_1 V(t) +r_1 R(t) \right) dt \\ & + \sigma_1 S(t) dB_1(t) +\int_{Y} \eta_1(u) S(t^-) \tilde{N}(dt,du), \end{split}\end{align*} \begin{align} \label{SVIR stochastic}\begin{split} dV(t) = & \left( v S(t) - \frac{\beta V(t) X(t)}{1+m_2 V(t)} -v_1 V(t) -\mu V(t) \right) dt + \sigma_2 V(t) dB_2(t) \\ & +\int_{Y} \eta_2(u) V(t^-) \tilde{N}(dt,du), \\ dI(t) = & \left( \frac{\beta S(t) X(t)}{1+m_1 S(t)} + \frac{\beta V(t) X(t)}{1+m_2 V(t)} -\left( \rho+\mu+\delta \right) I(t) \right) dt + \sigma_3 I(t) dB_3(t) \\ & +\int_{Y} \eta_3(u) I(t^-) \tilde{N}(dt,du), \\ dR(t) = & \left( \rho I(t) -\left( r_1 +\mu \right) R(t) \right) dt + \sigma_4 R(t) dB_4(t) +\int_{Y} \eta_4(u) R(t^-) \tilde{N}(dt,du), \\ dX(t) = & \xi \left( I(t) - X(t) \right) dt + \sigma_5 X(t) dB_5(t) +\int_{Y} \eta_5(u) X(t^-) \tilde{N}(dt,du), \end{split} \end{align} where notation $S(t^-), V(t^-), I(t^-), R(t^-), X(t^-)$ is used for the left limits of $S(t),$ $V(t) ,I (t), R(t), X (t)$, respectively, and $\eta_i : Y \times \Omega \rightarrow \mathbb{R}, i=1,2,3,4,5,$ are bounded and continuous with respect to $\lambda$ and $\mathcal{B}(Y) \times \mathcal{F}_{t}$-measurable, where $\mathcal{B}(Y)$ is a $\sigma$-algebra with respect to the set $Y$. Let $N(t)= S(t)+ V(t) +I (t)+ R(t), \,\, t\geq 0$, be the total population. To study long-term properties of infectious disease systems, the primary concern is to show that the solution is unique, positive, and global in time. In order to prove that there is a unique global positive solution to system \eqref{SVIR stochastic}, it is necessary to make the following hypotheses: \begin{itemize} \item $(\mathcal{H}1)$: Jump intensities $\eta_i (u)$ satisfy $\int_{Y} \eta^2_i (u) \lambda(du) < \infty \,\,\, (i = 1, 2, 3, 4, 5).$ \item ($\mathcal{H}2$): $1 + \eta_i(u) > 0$ and $\int_{Y} \left( \eta_i (u) - ln \left( 1 + \eta_i(u)\right) \right) \lambda(du) < \infty \,\,\, (i = 1, 2, 3, 4, 5).$ \item $(\mathcal{H}3)$: $p>2$, $b= \min \{ \mu, \mu+\delta-\xi, \xi \} -\frac{p-1}{2} \max_{i \in \{1,2,3,4,5\}} \{ \sigma_i^2 \} -\frac{1}{p} l\geq 0,$ where $l= \int_Y \max \{ [1+ \check{\eta}(u)]^p-1-p \check{\eta}(u), [1+ \hat{\eta}(u)]^p-1-p \hat{\eta}(u) \} \} \lambda(du),$ and $\check{\eta}(u)=\max_{ i \in \{1,2,3,4,5\}} \{ \eta_i(u) \}$, $ \hat{\eta}(u)=\min_{ i \in \{1,2,3,4,5\} } \{ \eta_i(u) \}$. \end{itemize} \begin{theorem} \label{theorem1.A} Let the assumptions $(\mathcal{H}1)$ and $(\mathcal{H}2)$ hold. For any given initial value $(S(0),V(0),I(0),R(0),X(0) )\in \mathbb{R}^5_+, $ there is a unique strong solution to SDEs\eqref{SVIR stochastic}, given by $(S(t),V(t),I(t),$ $R(t), X(t)), \, t \geq 0 $, and the solution will remain in $\mathbb{R}^5_+$ with probability one. \end{theorem} \begin{proof} Since the assumption $(\mathcal{H}1)$ ensures that the coefficients of system \eqref{SVIR stochastic} are locally Lipschitz continuous, the proof of the theorem is rather standard (see, for instance, \cite{Zhou_Yuan_2016}). The difference is that here the nonnegative $C^2$ function $F: \mathbb{R}^5_+ \rightarrow \mathbb{R}_+$ is defined as \begin{equation*} \begin{array}{lcl} F(S,V,I,R,X) & = & \left( S- a-a \ln \frac{S}{a} \right)+\left( V- a- \ln \frac{V}{a} \right)+\left( I- 1- \ln I \right) \\ & & +\left( R- 1- \ln R \right)+\frac{\mu+\delta}{\xi} \left( X- 1- \ln X \right), \end{array} \end{equation*} where $a>0$ is a positive constant to be determined later. Applying the operator $L$ \eqref{L ito} to the function $F$, we can get the following upper boundary: \begin{align*} LF &(S,V,I,R,X) \\ & \leq \gamma S \left( 1- \frac{S}{K}\right) - \mu V - (\mu + \delta) I - \mu R + \frac{\mu+\delta}{\xi} \xi \left( I-X \right) \\ & -a \gamma \left( 1- \frac{S}{K}\right) +a \frac{\beta X}{1+m_1 S} +av + a \frac{\beta X}{1+m_2 V} + a(v_1 + \mu) + r_1 + 3 \mu + \rho+ 2\delta \\ & + \frac{a \sigma_1^2}{2} + \frac{a \sigma_2^2}{2} + \frac{\sigma_3^2}{2} + \frac{\sigma_4^2}{2} + \frac{\mu+\delta}{\xi} \frac{\sigma_5^2}{2} + \sum_{i=1}^2 \int_Y a ( \eta_i (u) - \ln (1+\eta_i(u)) \lambda (du) \\ & + \sum_{i=3}^4 \int_Y ( \eta_i (u) - \ln (1+\eta_i(u)) \lambda (du) + \frac{\mu+\delta}{\xi} \int_Y ( \eta_i (u) - \ln (1+\eta_i(u)) \lambda (du) \\ \leq & \gamma (S-a) \left( 1- \frac{S}{K} \right) +(2a \beta - (\mu+\delta) ) X +C_1 \leq C_1+C_2, \end{align*} where we choose $a=\frac{\mu + \delta}{2 \beta}$ and as $ \max_{S \in \mathbb{R}_+} \{ \gamma (S-a) \left( 1- \frac{S}{K} \right) \} \leq C_2$. \end{proof} \section{Preliminaries} In this section, we shall introduce some necessary results which will be used in the following sections. For convenience, we recall some results on the stochastic differential equations with Lévy jumps, as we will need them in the sequel of the paper. For more details, we refer to \cite{Situ} and \cite{Oksendal_2005}. Let $(x(t), t \geq 0)$ be the stochastic hybrid system with jumps described by the following stochastic differential equations with Lévy jumps: \begin{equation} \label{p} \left\{ \begin{array}{l} dx(t)= b(x(t),t)dt + \sigma(x(t),t)dB(t)+\int_{Y} H(x(t),t,u) d\tilde{N}(dt,du), \\[1ex] x(0)=x_0, \, \, \, a(0)=a. \end{array} \right. \end{equation} Let $B(\cdot)$ be $d$-dimensional standard Brownian motions and $\bar{N}$ a Poisson counting measure, independent of $B$, and let them be defined on the given probability space $(\Omega, \mathcal{F}, \left\{ \mathcal{F}_{t} \right\}_{t\geq0},\textbf{P})$. $\tilde{N}(dt,du)= \bar{N}(dt,du) -\lambda(du)dt$ denotes the compensator of $\bar{N}$ and the characteristic measure $\lambda$ is defined on a measurable subset $Y$ of $(0, + \infty)$ satisfying $\lambda(Y)< + \infty$. $H(\cdot, \cdot, \cdot) : \mathbb{R}^n \times \mathbb{R}_+ \times Y \rightarrow \mathbb{R}^{n}$, and it is a right continuous function whose left limit exists. Furthermore, drift and diffusion coefficients are such that $b(\cdot, \cdot) : \mathbb{R}^n \times \mathbb{R}_+ \rightarrow \mathbb{R}^n, \sigma(\cdot, \cdot) : \mathbb{R}^n \times \mathbb{R}_+ \rightarrow \mathbb{R}^{n \times d} $ and $\sigma(x, t) \sigma^T (x, t) =: (d_{ij}(x, t)), \, i,j \in \left\lbrace 1,2 , \ldots , n \right\rbrace $. Let $F(x, t)$ be any twice continuously differentiable function with respect to $x$ and a continuously differentiable function with respect to $t$. The operator $L$ of system \eqref{p} is defined by \begin{align} \label{L ito} \begin{split} LF(x,t) = & \frac{\partial F(x,t)}{\partial t} + \sum\limits_{i=1}^{n} b_i(x,t) \frac{\partial F(x,t)}{\partial x_i} + \frac{1}{2} \sum\limits_{i=1}^{n} d_{ij}(x,t) \frac{\partial^2 F(x,t)}{\partial x_i \partial x_j} \\ & + \int\limits_{Y} (F(x+H(x,t,u),t)-F(x,t)-F_{x}(x,t)H(x,t,u)) \lambda (du) . \end{split} \end{align} By virtue of the generalized Itô’s formula - Theorem 1.16 in \cite{Oksendal_2005} (p.\,8), if $x(t) \in \mathbb{R}^n$, then \begin{align} \begin{split} dF(x,t) = & LF(x,t) dt + F_x(x,t) \sigma(x,t) dB(t) \\ & + \int_Y (F(x+H(x,t,u),t)-F(x,t)) \tilde{N} (dt,du). \end{split}\label{ito} \end{align} The following lemma will be useful in the analysis. \begin{definition} A stationary distribution for a Markov process is a probability measure $\vartheta$ on the measurable space $\left( \mathbb{R}^n, \mathcal{B} (\mathbb{R}^n) \right)$ that satisfies \[\int_{\mathbb{R}^n} P \left( s, x, B\right) \vartheta (dx)= \vartheta (B) , \qquad \text{for all} \qquad x \in \mathbb{R}^n. \] \end{definition} \begin{lemma}[Mutually exclusive possibilities, see \cite{Stettner,Meyn1993}] \label{lemaA1} Let $X(t) \in \mathbb{R}^d$ be a Feller process; then either an invariant probability measure $\vartheta $ exists, or \[ \lim_{t \rightarrow \infty} \sup_{\nu} \frac{1}{t} \int_0^t \int P(s, x, C) \nu (dx) ds =0,\] for any compact set $C \in \mathbb{R}^d$, where the supremum is taken over all initial distributions $\nu$ on $\mathbb{R}^d$, and $P(t, x, C)$ is the probability for $X(t) \in C$ with $X(0)=x$. \end{lemma} \section{Extinction} This section is devoted to establishing sufficient conditions so that the disease disappears from the population. Let us introduce notations \[ \left\langle A(t) \right\rangle := \frac{1}{t} \int_0^t A(s) \, ds , \, \, \, t \geq 0, \,\,\, \text{for an arbitrary function} \,\, A(t),\] \begin{equation} \label{R*aa} \begin{array}{lcl} R^*= \frac{ (1+ \sigma_3^2) \left( \frac{\beta }{m_1 }+\frac{\beta }{m_2 }\right) }{ \rho+\mu+\delta } . \end{array} \end{equation} \begin{lemma} \label{lim SVIRX} Assume that is assumptions $(\mathcal{H}1)$, $(\mathcal{H}2)$ and $(\mathcal{H}3)$ hold. Let $(S(t),V(t),I(t), $ $R(t),X(t))$ be the solution to model (7) with any given initial value $(S(0),V(0),I(0), $ $R(0), X(0)) \in \mathbb{R}^5_+$. Then \begin{align*} &\lim_{t \rightarrow \infty} \frac{S(t)}{t} =0, \quad \lim_{t \rightarrow \infty} \frac{V(t)}{t} =0, \quad \lim_{t \rightarrow \infty} \frac{I(t)}{t} =0, \quad \lim_{t \rightarrow \infty} \frac{R(t)}{t} =0, \\ &\lim_{t \rightarrow \infty} \frac{X(t)}{t} =0, \quad \lim_{t \rightarrow \infty} \frac{1}{t} \int_0^t S(s) dB_1(s) =0, \quad \lim_{t \rightarrow \infty} \frac{1}{t} \int_0^t V(s) dB_2(s) =0, \\ &\lim_{t \rightarrow \infty} \frac{1}{t} \int_0^t I(s) dB_3(s) =0, \quad \lim_{t \rightarrow \infty} \frac{1}{t} \int_0^t R(s) dB_4(s) =0, \quad \\ &\lim_{t \rightarrow \infty} \frac{1}{t} \int_0^t X(s) dB_5(s) =0, \quad \lim_{t \rightarrow \infty} \frac{1}{t} \int_0^t \int_Y \eta_1(u) S(s) \tilde{N}(ds,du) =0, \quad \\ &\lim_{t \rightarrow \infty} \frac{1}{t} \int_0^t \int_Y \eta_2(u) V(s) \tilde{N}(ds,du) =0, \quad \lim_{t \rightarrow \infty} \frac{1}{t} \int_0^t \int_Y \eta_3(u) I(s) \tilde{N}(ds,du) =0, \quad \\ &\lim_{t \rightarrow \infty} \frac{1}{t} \int_0^t \int_Y \eta_4(u) R(s) \tilde{N}(ds,du) =0, \quad \lim_{t \rightarrow \infty} \frac{1}{t} \int_0^t \int_Y \eta_5(u) X(s) \tilde{N}(ds,du) =0. \end{align*} \end{lemma} The proof can be obtained by the results from \cite{Zhou_Yuan_2016}, following proofs of lemmas 3.3 and 3.4, so it is omitted here. \begin{theorem} \label{theorem2.A} Let the assumptions $(\mathcal{H}1)$, $(\mathcal{H}2)$ and $(\mathcal{H}3)$ hold and let $(S(t),V(t),I(t),$ $R(t), X(t) ), \, t \geq 0 $ be the solution to system \eqref{SVIR stochastic} with any initial value $(S(0),V(0),$ $ I(0),R(0), X(0) ) \in \mathbb{R}^5_+ $. If $R^* < 1$, then the epidemic disappears from the population with probability one. That is, \[ \lim\limits_{t \rightarrow \infty} I(t) = 0, \, \, \, a.s. \] In addition, we have $ \lim\limits_{t \rightarrow \infty} R(t) = 0, \qquad \text{and} \qquad \lim\limits_{t \rightarrow \infty} X(t) = 0, \, \, a.s. $ \end{theorem} \begin{proof} Applying generalized Ito's formula to the function $I(t)+\sigma_3^2 \ln (I(t)+1)$, we obtain \begin{align*} L &\left( I(t) +\sigma_3^2 \ln (I(t)+1) \right) \leq \frac{\beta S(t) X(t)}{1+m_1 S(t)}+\frac{\beta V(t) X(t)}{1+m_2 V(t)}\\ & -\left( \rho+\mu+\delta \right) I(t) +\frac{\sigma_3^2}{I(t)+1} \left( \frac{\beta S(t) X(t)}{1+m_1 S(t)}+\frac{\beta V(t) X(t)}{1+m_2 V(t)}\right)\\ &+\int\limits_{Y} \left( \ln \left( 1+ \eta_3(u) \frac{I(s)}{I(s)+1} \right) - \eta_3(u) \frac{I(s)}{I(s)+1} \right) \lambda(du) \\ \leq &(1+\sigma_3^2) \left( \frac{\beta }{m_1 }+\frac{\beta }{m_2 }\right) X(t) -\left( \rho+\mu+\delta \right) I(t) . \\ % & + \sigma_3 dB(t) + \int\limits_{Y} \ln \left( 1+ \eta_3(u)\frac{I(s)}{I(s)+1} \right) \tilde{N}(dt,du) \\ \end{align*} Integrating both sides of the previous relation and dividing by $t$, we get \begin{align} \label{d ln Ia} \begin{split} & \hspace*{-1cm}\frac{I(t) -I(0)+\sigma_3^2(\ln (I(t)+1)- \ln (I(0)+1))}{t} \\ \leq & \frac{1}{t} (1+\sigma_3^2) \left( \frac{\beta }{m_1 }+\frac{\beta }{m_2 }\right) \int\limits_0^t X(s) ds \\ & - \frac{1}{t} \left( \rho+\mu+\delta \right) \int\limits_0^t I(s) ds + \frac{1}{t} \int\limits_0^t \sigma_3 (1+I(s)) dB(s) \\ & + \frac{1}{t} \int\limits_0^t \int\limits_{Y} \eta_3(u) I(s) + \ln \left( 1+ \eta_3(u)\frac{I(s)}{I(s)+1} \right) \tilde{N}(ds,du). \end{split} \end{align} According to the strong law of large numbers for local martingales and Lemma \ref{lim SVIRX}, one can easily obtain \begin{align} \label{noiseI} \begin{split} & \frac{1}{t} \int\limits_0^t \sigma_3 dB(s) \rightarrow 0, \,\,\,\, \,\,\,\,\, \frac{1}{t} \int\limits_0^t \int\limits_{Y} \ln \left( 1+ \eta_3(u) \frac{I(s)}{I(s)+1} \right) \tilde{N}(ds,du) \rightarrow 0, \,\,\,\, \,\,\,\,\, \\ & \frac{1}{t} \int\limits_0^t \sigma_3 I(s) dB(s) \rightarrow 0 \,\,\,\, \text{and} \,\,\,\,\, \frac{1}{t} \int\limits_0^t \int\limits_{Y} \eta_3(u) I(s) \tilde{N}(ds,du) \rightarrow 0, \,\,\,\,\,\,\, \text{as} \,\,\,\,\,\,\,\, t \rightarrow \infty. \end{split} \end{align} Integrating the fifth equation of model \eqref{SVIR stochastic} from $0$ to $t$, and then dividing by $t$ implies \begin{align} \label{e5} \begin{split} \frac{1}{t} \left( X(t)-X(0)\right) = & \xi \left\langle I(t) \right\rangle - \xi \left\langle X(t) \right\rangle + \frac{1}{t} \int\limits_0^t \sigma_5 X(s) dB(s) \\ & + \frac{1}{t} \int\limits_0^t \int\limits_{Y} \eta_5(u) X(s) \tilde{N}(ds,du). \end{split} \end{align} Let $\phi_1(t)=\frac{1}{t} \int\limits_0^t \sigma_5 X(s) dB(s) + \frac{1}{t} \int\limits_0^t \int\limits_{Y} \eta_5(u) X(s) \tilde{N}(ds,du) - \frac{1}{t} \left( X(t)-X(0)\right).$ By Lemma \ref{lim SVIRX}, we get $ \lim_{t \rightarrow \infty} \phi_1(t) =0.$ Taking the limit on both sides of \eqref{e5}, we have \begin{equation} \label{limXI} \lim_{t \rightarrow \infty} \left\langle X(t) \right\rangle = \lim_{t \rightarrow \infty} \left\langle I(t) \right\rangle. \end{equation} Similarly, integrating the sum of the first four equations of model \eqref{SVIR stochastic} from $0$ to $t$, and then dividing by $t$ gives \begin{align} \label{N A1} \begin{split} &\hspace*{-1cm}\frac{1}{t} \left( S(t)-S(0)\right) + \frac{1}{t} \left( V(t)-V(0)\right)+\frac{1}{t} \left( I(t)-I(0)\right)+\frac{1}{t} \left( R(t)-R(0)\right) \\ = & \left\langle \gamma S(t)\left( 1 - \frac{S(t)}{K}\right) \right\rangle - \mu \left\langle V(t) \right\rangle - (\mu+\delta) \left\langle I(t) \right\rangle - \mu \left\langle R(t) \right\rangle \\ & + \frac{1}{t} \int\limits_0^t \sigma_1 S(s) dB(s) + \frac{1}{t} \int\limits_0^t \int\limits_{Y} \eta_1(u) S(s) \tilde{N}(ds,du) \\ & + \frac{1}{t} \int\limits_0^t \sigma_2 V(s) dB(s) + \frac{1}{t} \int\limits_0^t \int\limits_{Y} \eta_2(u) V(s) \tilde{N}(ds,du) \\ & + \frac{1}{t} \int\limits_0^t \sigma_3 I(s) dB(s) + \frac{1}{t} \int\limits_0^t \int\limits_{Y} \eta_3(u) I(s) \tilde{N}(ds,du) \\ & + \frac{1}{t} \int\limits_0^t \sigma_4 R(s) dB(s) + \frac{1}{t} \int\limits_0^t \int\limits_{Y} \eta_4(u) R(s) \tilde{N}(ds,du) . \end{split} \end{align} Noticing that \begin{equation} \label{max Recruitment} \max_{S(t) \in \mathbb{R}_+} \left\lbrace \gamma S(t)\left( 1 - \frac{S(t)}{K}\right) \right\rbrace \leq \frac{\gamma K}{4}, \end{equation} the last two relations imply \begin{equation} \label{boundA} \begin{array}{lcl} \left\langle I(t) \right\rangle & \leq & \frac{1}{\mu + \delta} \frac{\gamma K}{4} + \frac{1}{\mu + \delta} \phi_2(t), \\ \end{array} \end{equation} where \begin{align} \label{fi2A} \begin{split} \phi_2(t) = & - \frac{1}{t} \left( S(t)-S(0)\right) + \frac{1}{t} \int\limits_0^t \sigma_1 S(s) dB(s) + \frac{1}{t} \int\limits_0^t \int\limits_{Y} \eta_1(u) S(s) \tilde{N}(ds,du) \\ & -\frac{1}{t} \left( V(t)-V(0)\right) + \frac{1}{t} \int\limits_0^t \sigma_2 V(s) dB(s) + \frac{1}{t} \int\limits_0^t \int\limits_{Y} \eta_2(u) V(s) \tilde{N}(ds,du) \\ & -\frac{1}{t} \left( I(t)-I(0)\right) + \frac{1}{t} \int\limits_0^t \sigma_3 I(s) dB(s) + \frac{1}{t} \int\limits_0^t \int\limits_{Y} \eta_3(u) I(s) \tilde{N}(ds,du) \\ & -\frac{1}{t} \left( R(t)-R(0)\right) + \frac{1}{t} \int\limits_0^t \sigma_4 R(s) dB(s) + \frac{1}{t} \int\limits_0^t \int\limits_{Y} \eta_4(u) R(s) \tilde{N}(ds,du). \end{split} \end{align} Again, by Lemma \ref{lim SVIRX}, it follows that $\limsup_{t \rightarrow \infty} \phi_2(t) =0.$ Taking the limit superior on both sides of \eqref{boundA}, we have \begin{equation} \label{limsupA} \begin{array}{lcl} \limsup_{t \rightarrow \infty} \left\langle I(t) \right\rangle & \leq & \frac{1}{\mu + \delta} \frac{\gamma K}{4}. \end{array} \end{equation} From \eqref{limXI} and \eqref{limsupA} it follows that \begin{equation} \label{limsup} \limsup_{t \rightarrow \infty} \left\langle X(t) \right\rangle \leq \frac{1}{\mu + \delta} \frac{\gamma K}{4}. \end{equation} Taking the limit superior on both sides of expression \eqref{d ln Ia} and using estimates \eqref{noiseI} and \eqref{limsup}, we obtain \begin{align*} \limsup\limits_{t \rightarrow \infty} \frac{I(t) + \sigma_3^2 \ln I(t)}{t} \leq & \left( (1+\sigma_3^2) \left( \frac{\beta }{m_1 }+\frac{\beta }{m_2 }\right) - \left( \rho+\mu+\delta \right) \right) \frac{\gamma K}{4(\mu + \delta)} \\ \leq & \left( \rho+\mu+\delta \right) \left( R^* -1 \right) \frac{\gamma K}{4(\mu + \delta)} < 0, \,\,\, a.s,\\ \limsup\limits_{t \rightarrow \infty} \frac{\sigma_3^2 \ln I(t)}{t} \leq & \limsup\limits_{t \rightarrow \infty} \frac{I(t) + \sigma_3^2 \ln I(t)}{t} \\ \leq &\left( \rho+\mu+\delta \right) \left( R^* -1 \right) \frac{\gamma K}{4(\mu + \delta)} < 0, \,\,\, a.s, \end{align*} by the theorem condition and the fact that $(S(t),V(t),I(t),R(t),X(t)) \in \mathbb{R}^5_+$. Then, \begin{equation*} \lim\limits_{t \rightarrow \infty} I(t) = 0, \,\,\, a.s. \end{equation*} Furthermore, this implies that $\lim\limits_{t \rightarrow \infty} R(t) = 0,$ and $\lim\limits_{t \rightarrow \infty} X(t) = 0$, a.s. The proof is completed. \end{proof} \section{Persistence in mean} In this section, sufficient conditions for strong persistence in mean of the disease will be established. \begin{definition} \label{Per.A} System \eqref{SVIR stochastic} is said to be strongly persistent in mean if $$ \liminf_{t \rightarrow \infty} \left\langle I(t) \right\rangle = \liminf_{t \rightarrow \infty} \frac{1}{t} \int\limits_0^t I(s) ds >0, \,\, \, a.s.\ $$ \end{definition} Let us define \begin{equation} \label{R0S A} R_0^S = \frac{\gamma + 2 \sqrt{v_1 v} -\frac{1}{2} \sigma_1^2 - \frac{1}{2} \sigma_2^2 - \sum_{i=1}^2 \int\limits_{Y} \left( \eta_i(u)-\ln \left( 1+ \eta_i(u) \right) \right) \lambda(du) }{\left( \frac{\beta }{m_2} \frac{1 }{\mu+ \delta} + \frac{v_1+\mu}{\mu} \right) \frac{\gamma^2}{4 v} +v + \mu +v_1 }. \end{equation} \begin{theorem} \label{theorem3.A} Let the assumptions $(\mathcal{H}1), (\mathcal{H}2)$, and $(\mathcal{H}3)$ hold and let $(S(t),V(t),I(t), $ $ R(t), X(t))$, $t \geq 0$ be the solution to system \eqref{SVIR stochastic} with any given initial value $(S(0),V(0), $ $I(0),R(0), X(0) ) \in \mathbb{R}_+^5 $. If $R_0^S >1,$ then \begin{align*} \label{I limit A} \liminf\limits_{t \rightarrow \infty} \frac{1}{t} \int\limits_{0}^t I(s) ds \geq & \frac{ \left( \left( \frac{\beta }{m_2} \frac{1 }{\mu+ \delta} + \frac{v_1+\mu}{\mu} \right) \frac{\gamma^2}{4 v} +v + \mu +v_1 \right) \left( R_0^S -1\right) }{2\beta } >0 \,\,\,\,\,\, a.s.,\\ \liminf\limits_{t \rightarrow \infty} \frac{1}{t} \int\limits_{0}^t\!\! R(s) ds \geq & \frac{\rho}{r_1+\mu} \frac{ \left( \left( \frac{\beta }{m_2} \frac{1 }{\mu+ \delta} + \frac{v_1+\mu}{\mu} \right) \frac{\gamma^2}{4 v} \! +\! v\! +\! \mu\! +\! v_1 \right) \left( R_0^S -1\right) }{2\beta } >0, \, \,\, a.s,\\ \liminf\limits_{t \rightarrow \infty} \frac{1}{t} \int\limits_{0}^t X(s) ds \geq & \frac{ \left( \left( \frac{\beta }{m_2} \frac{1 }{\mu+ \delta} + \frac{v_1+\mu}{\mu} \right) \frac{\gamma^2}{4 v} +v + \mu +v_1 \right) \left( R_0^S -1\right) }{2\beta } >0 \,\,\,\,\,\, a.s., \end{align*} that is, the disease is strongly persistent in mean. \end{theorem} \begin{proof} From the second equation of model \eqref{SVIR stochastic}, we obtain \begin{equation*} \hspace{-0.1cm} d(-V(t)) %& = & -v S(t) + \frac{\beta V(t) X(t)}{ 1+m_2 V(t) } + (v_1+\mu)V(t) - \sigma_2 V(t) dB(t) - \int\limits_{Y} \eta_2(u,t) V(t) d\tilde{N}(dt,du) \\ \leq \! (-v S(t) + \frac{\beta }{ m_2 }\! X(t) + (v_1+\mu)V(t)) dt - \sigma_2 V(t) dB(t) - \! \int\limits_{Y} \!\!\eta_2(u) V(t) \tilde{N}(dt,du) . \end{equation*} Integrating the previous equation from $0$ to $t$, and then dividing by $t$ implies \begin{align*} \frac{1}{t} \left( -V(t)+V(0)\right) \leq & - v \left\langle S(t) \right\rangle + \frac{\beta }{m_2} \left\langle X(t)\right\rangle + (v_1+\mu) \left\langle V(t)\right\rangle - \frac{1}{t} \int\limits_0^t \sigma_2 V(s) dB(s) \\ & - \frac{1}{t} \int\limits_0^t \int\limits_{Y} \eta_2(u) V(s) \tilde{N}(ds,du). \end{align*} Let $\phi_3(t)= \frac{1}{t} \left( V(t)-V(0)\right) - \frac{1}{t} \int\limits_0^t \sigma_2 V(s) dB(s) - \frac{1}{t} \int\limits_0^t \int\limits_{Y} \eta_2(u) V(s) \tilde{N}(ds,du) .$ Then \begin{equation} \label{e2A} \left\langle S(t) \right\rangle \leq \frac{\beta }{m_2 v} \left\langle X(t)\right\rangle + \frac{v_1+\mu}{v} \left\langle V(t)\right\rangle + \frac{1}{v} \phi_3. \end{equation} By Lemma \ref{lim SVIRX}, it follows that $\lim_{t \rightarrow \infty} \phi_3(t) =0.$ Taking the limit superior on both sides of \eqref{e2A}, we have \begin{equation} \label{limsup A} \limsup_{t \rightarrow \infty} \left\langle S(t) \right\rangle \leq \frac{\beta }{m_2 v} \limsup_{t \rightarrow \infty} \left\langle X(t) \right\rangle + \frac{v_1+\mu}{v} \limsup_{t \rightarrow \infty} \left\langle V(t) \right\rangle. \end{equation} Having in mind that relations \eqref{max Recruitment} and \eqref{fi2A} are satisfied, taking the limit superior on both sides of \eqref{N A1}, we have \begin{equation} \label{limsup A} \limsup_{t \rightarrow \infty} \left\langle V(t) \right\rangle \leq \frac{1 }{\mu} \frac{\gamma K}{4 }. \end{equation} Substituting \eqref{limsup} and \eqref{limsup A} in \eqref{limsup A}, we obtain \begin{equation} \label{limsup } \limsup_{t \rightarrow \infty} \left\langle S(t) \right\rangle \leq \left( \frac{\beta }{m_2 v} \frac{1 }{\mu+ \delta} + \frac{(v_1+\mu)}{v} \frac{1 }{\mu} \right) \frac{\gamma K}{4 }. \end{equation} Applying the operator $L$ to the function $- \ln S(t)- \ln V(t)$, we have \begin{align} \label{-lnS -lnV} \begin{split} &\hspace*{-1cm} L ( - \ln S(t)- \ln V(t)) %& = & - \gamma + \frac{\gamma S(t)}{ K } + \frac{\beta X(t)}{ 1+m_1 S(t) } +v - v_1\frac{V(t)}{S(t)} - r_1\frac{R(t)}{S(t)} + \frac{1}{2} \sigma_1^2 - v \frac{S(t)}{V(t)} +\frac{\beta X(t)}{ 1+m_2 V(t)} + \mu \\ %& & +v_1 + \frac{1}{2} \sigma_2^2 +\int\limits_{Y} \left( \eta_1-\ln \left( 1+ \eta_1(u) \right) \right) \lambda(du) +\int\limits_{Y} \left( \eta_2-\ln \left( 1+ \eta_2(u) \right) \right) \lambda(du) \\ %& & - \sigma_1 dB(t) + \int\limits_{Y} \ln \left( 1+ \eta_1(u,t) \right) d\tilde{N}(dt,du) \\ %& & - \sigma_2 dB(t) + \int\limits_{Y} \ln \left( 1+ \eta_2(u,t) \right) d\tilde{N}(dt,du) \\ %& \leq & - \gamma + \frac{\gamma }{ K } S(t) + 2 \beta X(t) +v - v_1\frac{V(t)}{S(t)} - v \frac{S(t)}{V(t)} + \mu +v_1 + \frac{1}{2} \sigma_1^2 + \frac{1}{2} \sigma_2^2 \\ %& & +\int\limits_{Y} \left( \eta_1-\ln \left( 1+ \eta_1(u) \right) \right) \lambda(du) +\int\limits_{Y} \left( \eta_2-\ln \left( 1+ \eta_2(u) \right) \right) \lambda(du) \\ %& & - \sigma_1 dB(t) + \int\limits_{Y} \ln \left( 1+ \eta_1(u,t) \right) d\tilde{N}(dt,du) \\ %& & - \sigma_2 dB(t) + \int\limits_{Y} \ln \left( 1+ \eta_2(u,t) \right) d\tilde{N}(dt,du) \\ \\ \leq &- \gamma + \frac{\gamma }{ K } S(t) + 2 \beta X(t) +v - 2 \sqrt{v_1 v} + \mu +v_1 \\ & + \frac{1}{2} \sigma_1^2 + \frac{1}{2} \sigma_2^2 + \sum_{i=1}^2 \int\limits_{Y} \left( \eta_i(u)-\ln \left( 1+ \eta_i(u) \right) \right) \lambda(du). \\ %& & - \sigma_1 dB(t) + \int\limits_{Y} \ln \left( 1+ \eta_1(u) \right) \tilde{N}(dt,du) \\ %& & - \sigma_2 dB(t) + \int\limits_{Y} \ln \left( 1+ \eta_2(u) \right) \tilde{N}(dt,du) . \end{split} \end{align} Integrating both sides of the previous inequality and dividing by $t$, we obtain \begin{align*} &\hspace*{-1cm} \frac {- \ln S(t)+\ln S(0)- \ln V(t)+\ln V(0)}{t} \\ \leq& - \gamma + \frac{\gamma }{ K } \left\langle S(t) \right\rangle + 2 \beta \left\langle X(t) \right\rangle +v - 2 \sqrt{v_1 v} + \mu +v_1 \\ & + \frac{1}{2} \sigma_1^2 + \frac{1}{2} \sigma_2^2 + \sum_{i=1}^2 \int\limits_{Y} \left( \eta_i(u)-\ln \left( 1+ \eta_i(u) \right) \right) \lambda(du) - \frac{1}{t} \int_{0}^{t} \sigma_1 dB(s) \\ & - \frac{1}{t} \int_{0}^{t} \sigma_2 dB(s) + \frac{1}{t} \sum_{i=1}^2 \int_{0}^{t} \int\limits_{Y} \ln \left( 1+ \eta_i(u) \right) \tilde{N}(ds,du) . \end{align*} Setting $\phi_4 = \frac{\ln S(t)-\ln S(0)+ \ln V(t)-\ln V(0)}{t} - \frac{1}{t} \int_{0}^{t} \sigma_1 dB(s) + \frac{1}{t} \int_{0}^{t} \int\limits_{Y} \ln \left( 1+ \eta_1(u) \right) \tilde{N}(ds,du)\\ - \frac{1}{t} \int_{0}^{t} \sigma_2 dB(s) + \frac{1}{t} \int_{0}^{t} \int\limits_{Y} \ln \left( 1+ \eta_2(u) \right) \tilde{N}(ds,du)$, we have \begin{align} \label{} \begin{split} 2 \beta \left\langle X(t) \right\rangle \geq & \gamma - \frac{\gamma }{ K } \left\langle S(t) \right\rangle -v + 2 \sqrt{v_1 v} - \int\limits_{Y} \left( \eta_1(u)-\ln \left( 1+ \eta_1(u) \right) \right) \lambda(du) \\ & -\mu -v_1 - \int\limits_{Y} \left( \eta_2(u)-\ln \left( 1+ \eta_2(u) \right) \right) \lambda(du) -\frac{1}{2} \sigma_1^2 - \frac{1}{2} \sigma_2^2 + \phi_4 .\\ \end{split} \end{align} By Lemma \ref{lim SVIRX}, it follows that $\lim_{t \rightarrow \infty} \phi_4(t) =0,$ taking the limit inferior on both sides of \eqref{}, one can obtain \begin{align*} 2 \beta \liminf_{t \rightarrow \infty} \left\langle X(t) \right\rangle % \geq & \gamma - \frac{\gamma }{ K } \limsup_{t \rightarrow \infty} \left\langle S(t) \right\rangle -v + 2 \sqrt{v_1 v} - \mu -v_1 - \frac{1}{2} \sigma_1^2 - \frac{1}{2} \sigma_2^2 \\ % & - \int\limits_{Y} \left( \eta_1(u)-\ln \left( 1+ \eta_1(u) \right) \right) \lambda(du) - \int\limits_{Y} \left( \eta_2(u)-\ln \left( 1+ \eta_2(u) \right) \right) \lambda(du) \\ \geq & \gamma - \left( \frac{\beta }{m_2} \frac{1 }{\mu+ \delta} + \frac{v_1+\mu}{\mu} \right) \frac{\gamma^2}{4 v} -v + 2 \sqrt{v_1 v} - \mu -v_1 - \frac{1}{2} \sigma_1^2 \\ & - \frac{1}{2} \sigma_2^2 - \sum_{i=1}^2 \int\limits_{Y} \left( \eta_i(u)-\ln \left( 1+ \eta_i(u) \right) \right) \lambda(du) \\ = & \left( \left( \frac{\beta }{m_2} \frac{1 }{\mu+ \delta} + \frac{v_1+\mu}{\mu} \right) \frac{\gamma^2}{4 v} +v + \mu +v_1 \right) \left( R_0^S -1\right) , \end{align*} where \eqref{limsup } is used and $R_0^S$ is defined by \eqref{R0S A}. Thus \begin{align*} \liminf\limits_{t \rightarrow \infty } \frac{1}{t } \int\limits_0^t X(s)ds & \geq & \frac{ \left( \left( \frac{\beta }{m_2} \frac{1 }{\mu+ \delta} + \frac{v_1+\mu}{\mu} \right) \frac{\gamma^2}{4 v} +v + \mu +v_1 \right) \left( R_0^S -1\right) }{2\beta } >0 \,\,\,\,\,\, a.s., \end{align*} by the theorem condition. From the last relation and \eqref{limXI}, one can conclude \begin{align*} \liminf\limits_{t \rightarrow \infty } \frac{1}{t } \int\limits_0^t I(s)ds \geq \frac{ \left( \left( \frac{\beta }{m_2} \frac{1 }{\mu+ \delta} + \frac{v_1+\mu}{\mu} \right) \frac{\gamma^2}{4 v} +v + \mu +v_1 \right) \left( R_0^S -1\right) }{2\beta } >0 \,\,\,\,\,\, a.s. \end{align*} Therefore, the infected class has a lower positive boundary as $R_0^S > 1$ holds, which implies that the infection will prevail for a long time. In addition, integrating the fourth equation from system \eqref{SVIR stochastic} from $0$ to $t$, dividing by $t$ and taking the limit inferior yields \begin{align*} \liminf\limits_{t \rightarrow \infty } \frac{1}{t } \int\limits_0^t \!\! R(s)ds \geq \frac{\rho}{r_1+\mu} \frac{ \left( \left( \frac{\beta }{m_2} \frac{1 }{\mu+ \delta} + \frac{v_1+\mu}{\mu} \right) \frac{\gamma^2}{4 v} +v + \mu +v_1 \right) \left( R_0^S -1\right) }{2\beta } >0 \,\,\,\,\,\, a.s. \end{align*} The proof is completed. \end{proof} \section{Stationary distribution} The aim of the following theorem is to give sufficient conditions under which stochastic model \eqref{SVIR stochastic} has a unique stationary distribution. \begin{theorem} \label{theorem4.A} Let the assumptions $(\mathcal{H}1), (\mathcal{H}2)$ and $(\mathcal{H}3)$ hold and let $(S(t),V(t),I(t), $ $R(t), X(t) ), \, t \geq 0$ be the solution to \eqref{SVIR stochastic} with any initial value $(S(0),V(0),I(0),R(0),$ $ X(0) ) \in \mathbb{R}_+^5 $. If $R_0^S >1$, then stochastic model \eqref{SVIR stochastic} admits a unique stationary distribution which has the ergodic property. \end{theorem} \begin{proof} In order to prove this theorem, we only need to prove that conditions in Lemma \ref{lemaA1} are true for system \eqref{SVIR stochastic}. According to Theorem 4.1 in \cite{Xi_Zhu_2017}, it follows that $(S(t), V(t), I(t), R(t), X(t)), t \geq 0,$ has the strong Feller property. Let us define a $C^2$-function $\tilde{W}(S,V,I,R,X) : \mathbb{R}_+^5 \rightarrow \mathbb{R}$ as follows: \begin{align*} \tilde{W}(S,V,I,R,X) = & M \left( -\ln S-\ln V + \frac{2 \beta}{\rho+\mu +\delta} \left( V+I+ \frac{\rho+\mu +\delta}{\xi} X\right) \right) \\ & + \frac{1}{2} \left(S+V+I+R \right)^2 +\frac{1}{2 \xi}X^2 - \ln V - \ln I- \ln R. \end{align*} Here $M>0$ is a constant that will be explicitely determined later. Function $\tilde{W}$ is continuous in $(S,V,I,R,X)$, which takes the minimum at the point $\tilde{W}(S(t_0),V(t_0),I(t_0),$ $R(t_0),X(t_0))$. Define a nonnegative $C^2$-function $W(S,V,I,R,X) : \mathbb{R}_+^5 \rightarrow \mathbb{R}$ as \begin{align*} &\hspace*{-0.3cm} W(S,V,I,R,X)\\ \quad =& M \left( -\ln S-\ln V + \frac{2 \beta}{\rho+\mu +\delta} \left( V+I+ \frac{\rho+\mu +\delta}{\xi} X\right) \right) +\frac{1}{2 \xi}X^2 \\ \quad & +\frac{1}{2} \left(S+V+I+R\right)^2 \! -\! \ln V\! -\! \ln I \! -\! \ln R \! -\! \tilde{W}(S(t_0),V(t_0),I(t_0),R(t_0),X(t_0)). \end{align*} Let us denote \begin{align*} W_1(S,V,I,R,X)= &-\ln S-\ln V + \frac{2 \beta}{\rho+\mu +\delta} \left( V+I+ \frac{\rho+\mu +\delta}{\xi} X\right), \\ W_2(S,V,I,R,X)=&\frac{1}{2} \left(S+V+I+R \right)^2 +\frac{1}{2 \xi}X^2 , \\ W_3(S,V,I,R,X)= & - \ln V - \ln I - \ln R. \end{align*} Applying the operator $L$ to the function $W_1(S,V,I,R,X)$, we obtain \begin{align} \label{LW1} \begin{split} &\hspace*{-1cm}LW_1(S,V,I,R,X) \\ \leq & 2 \beta X + \frac{\gamma }{ K } S -\left( \frac{\beta }{m_2} \frac{1 }{\mu+ \delta} + \frac{v_1+\mu}{\mu} \right) \frac{\gamma^2}{4 v} \\ & - \left( \left( \frac{\beta }{m_2} \frac{1 }{\mu+ \delta} + \frac{v_1+\mu}{\mu} \right) \frac{\gamma^2}{4 v} +v + \mu +v_1 \right) \left( R_0^S -1\right) \\ & + \frac{2 \beta}{\rho+\mu +\delta} \left( vS + \beta \frac{SX}{1+m_1S} - \left( \rho+\mu +\delta\right) X\right) \\ = : & \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) S + \frac{2 \beta^2 SX}{(\rho+\mu +\delta)(1+m_1S)} - \psi_1 , \end{split} \end{align} where we used \eqref{-lnS -lnV} and $R_0^S$ is defined by \eqref{R0S A}. Applying the operator $L$ to the function $W_2(S,V,I,R,X)$, we obtain \begin{align} \label{LW2} \begin{split} & \hspace{-0.3cm} LW_2 (S,V,I,R,X) \\ =& \left(S+V+I+R \right) \left( \gamma S \left(1-\frac{S}{K} \right)-\mu V - (\mu+\delta)I -\mu R \right) \\ & +\frac{1}{2} \left( \sigma_1^2 S^2 + \sigma_2^2 V^2+ \sigma_3^2 I^2+ \sigma_4^2 R^2 \right) + \frac{1}{2} \int_{Y} S(t) \left( \eta_1(u) - \ln (1+\eta_1(u)) \right) \lambda(du) \\ & + \frac{1}{2} \int_{Y} V(t) \left( \eta_2(u) - \ln (1+\eta_2(u)) \right) \lambda(du) \\ &+ \frac{1}{2} \int_{Y} I(t) \left( \eta_3(u) - \ln (1+\eta_3(u)) \right) \lambda(du) \\ & + \frac{1}{2} \int_{Y} R(t) \left( \eta_4(u) - \ln (1+\eta_4(u)) \right) \lambda(du) \\ & + X \left(I-X \right) + \frac{1}{2 \xi} \sigma_5^2 X^2 + \frac{1}{2 \xi} \int_{Y} X(t) \left( \eta_5(u) - \ln (1+\eta_5(u)) \right) \lambda(du) \\ %& \leq & \gamma S \left(S+V+I+R \right) -\frac{\gamma}{K} S^3 -\mu V^2 - (\mu+\delta)I^2 -\mu R^2 \\ %& & +\frac{1}{2} \max \left\lbrace \sigma_1^2 ,\sigma_2^2 , \sigma_3^2 , \sigma_4^2 \right\rbrace \left( S^2 + V^2+ I^2+ R^2 \right) + \frac{1}{2} \int_{Y} S(t) \left( \eta_1(u) - \ln (1+\eta_1(u)) \right) \lambda(du) \\ %& & + \frac{1}{2} \int_{Y} V(t) \left( \eta_2(u) - \ln (1+\eta_2(u)) \right) \lambda(du) + \frac{1}{2} \int_{Y} I(t) \left( \eta_3(u) - \ln (1+\eta_3(u)) \right) \lambda(du) \\ %& & + \frac{1}{2} \int_{Y} R(t) \left( \eta_4(u) - \ln (1+\eta_4(u)) \right) \lambda(du) \\ %& & - X^2+IX + \frac{1}{2 \xi} \sigma_5^2 X^2 + \frac{1}{2 \xi} \int_{Y} X(t) \left( \eta_5(u) - \ln (1+\eta_5(u)) \right) \lambda(du) \\ \leq & -\frac{\gamma}{2 K} S^3 -\frac{1}{2}\mu V^2 - \frac{1}{2}(\mu+\delta)I^2 -\frac{1}{2}\mu R^2 -\frac{1}{2} X^2 + C, \end{split} \end{align} where \begin{align*} C =&\max_{(S,V,I,R,X) \in \mathbb{R}_+^5 }\left\lbrace \gamma S \left(S+V+I+R \right) -\frac{\gamma}{2 K} S^3 -\frac{1}{2}\mu V^2 -\frac{1}{2} (\mu+\delta)I^2 \right. \\ & -\frac{1}{2}\mu R^2 +\frac{1}{2} \max \left\lbrace \sigma_1^2 ,\sigma_2^2 , \sigma_3^2 , \sigma_4^2, \sigma_5^2 \right\rbrace \left( S^2 + V^2+ I^2+ R^2+ \frac{X^2}{\xi} \right) \\ & + \frac{1}{2} \int_{Y} S(t) \left( \eta_1(u) - \ln (1+\eta_1(u)) \right) \lambda(du)\\ & + \frac{1}{2} \int_{Y} V(t) \left( \eta_2(u) - \ln (1+\eta_2(u)) \right) \lambda(du) \\ & + \frac{1}{2} \int_{Y} I(t) \left( \eta_3(u) - \ln (1+\eta_3(u)) \right) \lambda(du)\\ & + \frac{1}{2} \int_{Y} R(t) \left( \eta_4(u) - \ln (1+\eta_4(u)) \right) \lambda(du) \\ & \left. - \frac{1}{2}X^2+IX + \frac{1}{2 \xi} \int_{Y} X(t) \left( \eta_5(u) - \ln (1+\eta_5(u)) \right) \lambda(du) \right\rbrace. \end{align*} Applying the operator $L$ to $W_3(S,V,I,R,X)$ , we obtain \begin{align} \label{LW3} LW_3(S,V,I,R,X) %& = & -\frac{1}{V} \left( v S - \frac{\beta V X}{1+m_2 V} -v V -\mu V \right) +\frac{1}{2 } \sigma_2^2 + \int_{Y} \left( \eta_2(u) - \ln (1+\eta_2(u)) \right) \lambda(du) \\ %& & -\frac{1}{R} \left( \rho I -\left( r_1 +\mu \right) R \right) +\frac{1}{2 } \sigma_4^2 + \int_{Y} \left( \eta_4(u) - \ln (1+\eta_4(u)) \right) \lambda(du) \\ % & \leq & -\frac{v S}{V} + \frac{\beta X}{1+m_2 V} +v +\mu +\frac{1}{2 } \sigma_2^2 + \int_{Y} \left( \eta_2(u) - \ln (1+\eta_2(u)) \right) \lambda(du) \\ % & & -\frac{\rho I}{R} + r_1 +\mu +\frac{1}{2 } \sigma_4^2 + \int_{Y} \left( \eta_4(u) - \ln (1+\eta_4(u)) \right) \lambda(du) \\ % & \leq & -\frac{v S}{V} + \frac{\beta X}{1+m_2 V} +v +\mu +\frac{1}{2 } \sigma_2^2 + \int_{Y} \left( \eta_2(u) - \ln (1+\eta_2(u)) \right) \lambda(du) \\ % & & -\frac{\rho I}{R} + r_1 +\mu +\frac{1}{2 } \sigma_4^2 + \int_{Y} \left( \eta_4(u) - \ln (1+\eta_4(u)) \right) \lambda(du) \\ % & \leq & -\frac{v S}{V} - \frac{\beta XS}{(1+m_1 S)I}-\frac{\rho I}{R} +v +3 \mu +\delta+\rho +\frac{1}{2 } \sigma_2^2 +\frac{1}{2 } \sigma_3^2 + r_1 +\frac{1}{2 } \sigma_4^2 + 3 \alpha, \\ \leq -\frac{v S}{V} - \frac{\beta XS}{(1+m_1 S)I} -\frac{\rho I}{R} +q, \end{align} where $q=+v +3 \mu +\delta+\rho +\frac{1}{2 } \sigma_2^2 +\frac{1}{2 } \sigma_3^2 + r_1 +\frac{1}{2 } \sigma_4^2 + 3 \alpha$ and the condition $(\mathcal{H}2)$ is used, that is, there is $\alpha >0$, so that $ \int_{Y} \left( \eta_i(u) - \ln (1+\eta_i(u)) \right) \lambda(du) \leq \alpha, \,\,\, i=1,2,3,4,5.$ Then, from \eqref{LW1}, \eqref{LW2} and \eqref{LW3} it follows that \begin{align} \label{LWaa} \begin{split} &\hspace*{-0.3cm}LW (S,V,I,R,X)\\ \leq& M \left( \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) S + \frac{2 \beta^2 SX}{(\rho+\mu +\delta)(1+m_1S)} - \psi_1\right) -\frac{\gamma}{2 K} S^3 \\ & -\frac{1}{2}\mu V^2 - \frac{1}{2}(\mu+\delta)I^2 -\frac{1}{2}\mu R^2 -\frac{1}{2} X^2 + C -\frac{v S}{V} - \frac{\beta SX}{(1+m_1 S) I} -\frac{\rho I}{R} +q. \end{split} \end{align} Let $M>0$ be a sufficiently large positive constant such that \begin{equation} \label{cond MA1a} -M \psi_1 + A_1 \leq -2, \end{equation} where $A_1$ is defined by \eqref{A1a}. Define a bounded closed set as follows: \begin{align*} \hspace{-0.3cm} D= \left\lbrace \left( S, V,I,R,X \right) \in \mathbb{R}_+^5\right. :& \epsilon \leq S \leq \frac{1}{\epsilon},\epsilon^2 \leq V \leq \frac{1}{\epsilon^2},\\ &\left. \epsilon \leq I \leq \frac{1}{\epsilon},\epsilon^2 \leq R \leq \frac{1}{\epsilon^2},\epsilon \leq X \leq \frac{1}{\epsilon} \right\rbrace , \end{align*} where $0 < \epsilon < 1$ is a sufficiently small constant such that the following conditions hold: \begin{align} \label{uD1} M \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) \epsilon +& \ \sup_{X \in \mathbb{R}_+} \left\lbrace M \frac{2 \beta^2 \epsilon X}{\rho+\mu +\delta} -\frac{1}{4} X^2 \right\rbrace - M \psi_1 + A_1 \leq -1 ,\\ \label{uD2} -\frac{v }{\epsilon} + A_2 \leq &-1 , \\ \label{uD5} \frac{2 M \beta^2 }{\rho+\mu +\delta} \epsilon +&\sup_{S \in \mathbb{R}_+} \left\lbrace \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) M S -\frac{\gamma}{4 K} S^3 \right\rbrace - M \psi_1 + A_3 \leq -1 , \\ \label{uD3} -\frac{\beta }{(1+m_1 \epsilon)\epsilon} + A_2 \leq& -1 ,\\\label{uD4} -\frac{\rho }{\epsilon} + A_2 \leq& -1 , \\ \label{uD6} -\frac{\gamma}{4 K} \frac{1}{\epsilon^9} + A_4 \leq& -1 , \\ \label{uD7} -\frac{1}{4}\mu \frac{1}{\epsilon^8} + A_5 \leq& -1 , \\ \label{uD8} - \frac{1}{4}(\mu+\delta)\frac{1}{\epsilon^2} + A_6 \leq& -1, \\ \label{uD9} -\frac{1}{4}\mu \frac{1}{\epsilon^4} + A_7 \leq &-1 , \\ \label{uD10} -\frac{1}{4} \frac{1}{\epsilon^2} + A_8 \leq& -1 , \end{align} where $A_2 - A_8$ are constants which can be found from expressions \eqref{A2a}, \eqref{A3a}, \eqref{A4a}, \eqref{A5a}, \eqref{A6a}, \eqref{A7a} and \eqref{A8a}, respectively. Note that for a sufficiently small $ \epsilon $, conditions \eqref{uD1} and \eqref{uD5} hold due to \eqref{cond MA1a}. For convenience, $\mathbb{R}_+^5 \setminus D$ will be divided into the following ten domains: \begin{align*} D_1=& \left\lbrace \left( S, V,I,R,X \right) \in \mathbb{R}_+^5: S < \epsilon \right\rbrace, \ D_2= \left\lbrace \left( S, V,I,R,X \right) \in \mathbb{R}_+^5 : S \geq \epsilon, \, V < \epsilon^2 \right\rbrace, \\ D_3=& \left\lbrace \left( S, V,I,R,X \right) \in \mathbb{R}_+^5: X < \epsilon \right\rbrace, \ D_4= \left\lbrace \! \left( S, V,I,R,X \right) \in \mathbb{R}_+^5 : S,X \geq \epsilon, I\!< \epsilon^3 \!\right\rbrace, \\ D_5=&\left\lbrace \left( S, V,I,R,X \right) \in \mathbb{R}_+^5: I \geq \epsilon^3, \, R < \epsilon^4 \right\rbrace, \ D_6= \left\lbrace \left( S, V,I,R,X \right) \in \mathbb{R}_+^5 :\! S\!>\!\frac{1}{\epsilon} \!\right\rbrace, \\ D_7= &\left\lbrace \left( S, V,I,R,X \right) \in \mathbb{R}_+^5: V>\frac{1}{\epsilon^2}\right\rbrace , \ D_8= \left\lbrace \left( S, V,I,R,X \right) \in \mathbb{R}_+^5 : I>\frac{1}{\epsilon} \right\rbrace , \\ D_9= &\left\lbrace \left( S, V,I,R,X \right) \in \mathbb{R}_+^5: R>\frac{1}{\epsilon^2}\right\rbrace , \ D_{10}= \left\lbrace \left( S, V,I,R,X \right) \in \mathbb{R}_+^5 : X>\frac{1}{\epsilon}\right\rbrace . \end{align*} Next, we will prove that $LW (S, V,I,R,X) \leq -1$ for any $(S, V,I,R,X) \in \mathbb{R}_+^5 \setminus D$, which is equivalent to proving it on the above ten domains. Case 1. Having in mind estimate \eqref{LWaa}, for any $(S,V, I, R,X) \in D_1 $, we have \begin{align*} &\hspace{-0.3cm}LW(S,V,I,R,X)\\ %& \leq & M \left( \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) S + \frac{2 \beta^2 SX}{(\rho+\mu +\delta)(1+m_1S)} - \psi_1\right) -\frac{\gamma}{2 K} S^3 -\frac{1}{2}\mu V^2 \\ %& & - \frac{1}{2}(\mu+\delta)I^2 -\frac{1}{2}\mu R^2 -\frac{1}{2} X^2 + C -\frac{v S}{V} + \frac{\beta X}{1+m_2 V} +v +2 \mu +\frac{1}{2 } \sigma_2^2 \\ % & & -\frac{\rho I}{R} + r_1 +\frac{1}{2 } \sigma_4^2 + 2 \alpha \\ %& \leq & M \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) S + M \frac{2 \beta^2 SX}{\rho+\mu +\delta} - M \psi_1 -\frac{1}{4} X^2 + A_1 \\ %& \leq & M \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) \epsilon + M \frac{2 \beta^2 \epsilon X}{\rho+\mu +\delta} - M \psi_1 -\frac{1}{4} X^2 + A_1 \\ & \leq M \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) \epsilon + \sup_{X \in \mathbb{R}_+} \left\lbrace M \frac{2 \beta^2 \epsilon X}{\rho+\mu +\delta} -\frac{1}{4} X^2 \right\rbrace - M \psi_1 + A_1, \end{align*} where \begin{equation} \label{A1a} A_1= \sup_{(S,V, I, R,X) \in \mathbb{R}_+^5} \left\lbrace -\frac{\gamma}{2 K} S^3 -\frac{\mu}{2} V^2 - \frac{\mu+\delta}{2} I^2 -\frac{\mu}{2} R^2 -\frac{1}{4} X^2 + C +q \right\rbrace . \end{equation} By inequality \eqref{uD1}, we obtain $L W(S,V,I,R,X) \leq -1$ for all $(S,V,I,R,X) \in D_1$. Case 2. For any $(S,V, I, R,X) \in D_2 $, we have \begin{align*} LW(S,V,I,R,X) %& \leq & M \left( \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) S + \frac{2 \beta^2 SX}{(\rho+\mu +\delta)(1+m_1S)} - \psi_1\right) -\frac{\gamma}{2 K} S^3 -\frac{1}{2}\mu V^2 - \frac{1}{2}(\mu+\delta)I^2 \\ %& & -\frac{1}{2}\mu R^2 -\frac{1}{2} X^2 + C -\frac{v S}{V} + \frac{\beta X}{1+m_2 V} +v +2 \mu +\frac{1}{2 } \sigma_2^2 -\frac{\rho I}{R} + r_1 +\frac{1}{2 } \sigma_4^2 + 2 \alpha \\ \leq -\frac{v S}{V} + A_2 \leq -\frac{v }{\epsilon} + A_2, \end{align*} where \begin{align} \label{A2a} \begin{split} A_2 = & \sup_{(S,V, I, R,X) \in \mathbb{R}_+^5} \left\lbrace M \left( \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) S + \frac{2 \beta^2 SX}{(\rho+\mu +\delta)(1+m_1S)} - \psi_1\right) \right.\\ & \left. -\frac{\gamma}{2 K} S^3 -\frac{1}{2}\mu V^2 - \frac{1}{2}(\mu+\delta)I^2 -\frac{1}{2}\mu R^2 -\frac{1}{2} X^2 + C +q \right\rbrace . \end{split} \end{align} By inequality \eqref{uD2}, we obtain $L W(S,V,I,R,X) \leq -1$ for all $(S,V,I,R,X) \in D_2$. Case 3. For any $(S,V, I, R,X) \in D_3 $, we have \begin{align*} &\hspace*{-0.3cm}LW(S,V,I,R,X) \\%& \leq & M \left( \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) S + \frac{2 \beta^2 SX}{(\rho+\mu +\delta)(1+m_1S)} - \psi_1\right) -\frac{\gamma}{2 K} S^3 -\frac{1}{2}\mu V^2 - \frac{1}{2}(\mu+\delta)I^2 \\ %& & -\frac{1}{2}\mu R^2 -\frac{1}{2} X^2 + C -\frac{v S}{V} + \frac{\beta X}{1+m_2 V} +v +2 \mu +\frac{1}{2 } \sigma_2^2 -\frac{\rho I}{R} + r_1 +\frac{1}{2 } \sigma_4^2 + 2 \alpha \\ \leq & M \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) S +M \frac{2 \beta^2 X}{(\rho+\mu +\delta)} - M \psi_1 -\frac{\gamma}{4 K} S^3 + A_3 \\ \leq & \frac{2 M \beta^2 }{\rho+\mu +\delta} \epsilon +\sup_{S \in \mathbb{R}_+} \left\lbrace \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) M S -\frac{\gamma}{4 K} S^3 \right\rbrace - M \psi_1 + A_3, \end{align*} where \begin{align} \label{A3a} A_3= \sup_{(S,V, I, R,X) \in \mathbb{R}_+^5} \left\lbrace -\frac{\gamma}{4 K} S^3 -\frac{\mu}{2} V^2 - \frac{\mu+\delta}{2} I^2 -\frac{\mu}{2} R^2 -\frac{1}{2} X^2 + C +q \right\rbrace . \end{align} By inequality \eqref{uD5}, we obtain $L W(S,V,I,R,X) \leq -1$ for all $(S,V,I,R,X) \in D_3$. Case 4. For any $(S,V, I, R,X) \in D_4 $, we have \begin{align*} LW(S,V,I,R,X) %& \leq & M \left( \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) S + \frac{2 \beta^2 SX}{(\rho+\mu +\delta)(1+m_1S)} - \psi_1\right) -\frac{\gamma}{2 K} S^3 -\frac{1}{2}\mu V^2 - \frac{1}{2}(\mu+\delta)I^2 \\ %& & -\frac{1}{2}\mu R^2 -\frac{1}{2} X^2 + C -\frac{v S}{V} + \frac{\beta X}{1+m_2 V} +v +2 \mu +\frac{1}{2 } \sigma_2^2 -\frac{\rho I}{R} + r_1 +\frac{1}{2 } \sigma_4^2 + 2 \alpha \\ \leq -\frac{\beta SX}{(1+ m_1 S)I} + A_2 \leq -\frac{\beta }{(1+m_1 \epsilon)\epsilon} + A_2 . \end{align*} By inequality \eqref{uD3}, we obtain $L W(S,V,I,R,X) \leq -1$ for all $(S,V,I,R,X) \in D_4$. Case 5. For any $(S,V, I, R,X) \in D_5 $, we have \begin{align*} LW(S,V,I,R,X) %& \leq & M \left( \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) S + \frac{2 \beta^2 SX}{(\rho+\mu +\delta)(1+m_1S)} - \psi_1\right) -\frac{\gamma}{2 K} S^3 -\frac{1}{2}\mu V^2 - \frac{1}{2}(\mu+\delta)I^2 \\ %& & -\frac{1}{2}\mu R^2 -\frac{1}{2} X^2 + C -\frac{v S}{V} + \frac{\beta X}{1+m_2 V} +v +2 \mu +\frac{1}{2 } \sigma_2^2 -\frac{\rho I}{R} + r_1 +\frac{1}{2 } \sigma_4^2 + 2 \alpha \\ \leq -\frac{\rho I}{R} + A_2 \leq -\frac{\rho }{\epsilon} + A_2 . \end{align*} By inequality \eqref{uD4}, we obtain $L W(S,V,I,R,X) \leq -1$ for all $(S,V,I,R,X) \in D_5$. Case 6. For any $(S,V, I, R,X) \in D_6 $, we have \begin{align*} LW(S,V,I,R,X) %& \leq & M \left( \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) S + \frac{2 \beta^2 SX}{(\rho+\mu +\delta)(1+m_1S)} - \psi_1\right) -\frac{\gamma}{2 K} S^3 -\frac{1}{2}\mu V^2 - \frac{1}{2}(\mu+\delta)I^2 \\ %& & -\frac{1}{2}\mu R^2 -\frac{1}{2} X^2 + C -\frac{v S}{V} + \frac{\beta X}{1+m_2 V} +v +2\mu +\frac{1}{2 } \sigma_2^2 -\frac{\rho I}{R} + r_1 +\frac{1}{2 } \sigma_4^2 + 2 \alpha \\ \leq -\frac{\gamma}{4 K} S^3 + A_4 \leq -\frac{\gamma}{4 K} \frac{1}{\epsilon^9} + A_4, \end{align*} where \begin{align} \label{A4a} \begin{split} A_4 =& \sup_{(S,V, I, R,X) \in \mathbb{R}_+^5} \left\lbrace M \left( \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) S + \frac{2 \beta^2 SX}{(\rho+\mu +\delta)(1+m_1S)} - \psi_1\right) \right.\\ & \left. -\frac{\gamma}{4 K} S^3 -\frac{1}{2}\mu V^2 - \frac{1}{2}(\mu+\delta)I^2 -\frac{1}{2}\mu R^2 -\frac{1}{2} X^2 + C + q \right\rbrace . \end{split} \end{align} In view of \eqref{uD6}, it follows that $L W(S,V,I,R,X) \leq -1$ for all $ (S,V,I,R,X) \in D_6$. Case 7. For any $(S,V, I, R,X) \in D_7 $, we have \begin{align*} LW(S,V,I,R,X) %& \leq & M \left( \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) S + \frac{2 \beta^2 SX}{(\rho+\mu +\delta)(1+m_1S)} - \psi_1\right) -\frac{\gamma}{2 K} S^3 -\frac{1}{2}\mu V^2 - \frac{1}{2}(\mu+\delta)I^2 \\ %& & -\frac{1}{2}\mu R^2 -\frac{1}{2} X^2 + C -\frac{v S}{V} + \frac{\beta X}{1+m_2 V} +v +2 \mu +\frac{1}{2 } \sigma_2^2 -\frac{\rho I}{R} + r_1 +\frac{1}{2 } \sigma_4^2 + 2 \alpha \\ \leq -\frac{1}{4}\mu V^2 + A_5 \leq -\frac{\mu}{4} \frac{1}{\epsilon^8} + A_5, \end{align*} where \begin{align} \label{A5a} \begin{split} A_5 =& \sup_{(S,V, I, R,X) \in \mathbb{R}_+^5} \left\lbrace M \left( \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) S + \frac{2 \beta^2 SX}{(\rho+\mu +\delta)(1+m_1S)} - \psi_1\right) \right.\\ & \left. -\frac{\gamma}{2 K} S^3 -\frac{1}{4}\mu V^2 - \frac{1}{2}(\mu+\delta)I^2 -\frac{1}{2}\mu R^2 -\frac{1}{2} X^2 + C +q \right\rbrace . \end{split} \end{align} In view of \eqref{uD7}, it follows that $L W(S,V,I,R,X) \leq -1$ for all $ (S,V,I,R,X) \in D_7$. Case 8. For any $(S,V, I, R,X) \in D_8 $, we have \begin{align*} LW(S,V,I,R,X) %& \leq & M \left( \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) S + \frac{2 \beta^2 SX}{(\rho+\mu +\delta)(1+m_1S)} - \psi_1\right) -\frac{\gamma}{2 K} S^3 -\frac{1}{2}\mu V^2 - \frac{1}{2}(\mu+\delta)I^2 \\ %& & -\frac{1}{2}\mu R^2 -\frac{1}{2} X^2 + C -\frac{v S}{V} + \frac{\beta X}{1+m_2 V} +v +2 \mu +\frac{1}{2 } \sigma_2^2 -\frac{\rho I}{R} + r_1 +\frac{1}{2 } \sigma_4^2 + 2 \alpha \\ \leq - \frac{\mu+\delta}{4} I^2 + A_6 \leq -\frac{\mu+\delta}{4} \frac{1}{\epsilon^2} + A_6, \end{align*} where \begin{equation} \label{A6a} \begin{split} A_6 =& \sup_{(S,V, I, R,X) \in \mathbb{R}_+^5} \left\lbrace M \left( \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) S + \frac{2 \beta^2 SX}{(\rho+\mu +\delta)(1+m_1S)} - \psi_1\right) \right.\\ & \left. -\frac{\gamma}{2 K} S^3 -\frac{1}{2}\mu V^2 - \frac{1}{4}(\mu+\delta)I^2 -\frac{1}{2}\mu R^2 -\frac{1}{2} X^2 + C + q \right\rbrace . \end{split} \end{equation} In view of \eqref{uD8}, it follows that $L W(S,V,I,R,X) \leq -1$ for $(S,V,I,R,X) \in D_8$. Case 9. For any $(S,V, I, R,X) \in D_9 $, we have \begin{align*} LW(S,V,I,R,X) %& \leq & M \left( \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) S + \frac{2 \beta^2 SX}{(\rho+\mu +\delta)(1+m_1S)} - \psi_1\right) -\frac{\gamma}{2 K} S^3 -\frac{1}{2}\mu V^2 - \frac{1}{2}(\mu+\delta)I^2 \\ %& & -\frac{1}{2}\mu R^2 -\frac{1}{2} X^2 + C -\frac{v S}{V} + \frac{\beta X}{1+m_2 V} +v +2 \mu +\frac{1}{2 } \sigma_2^2 -\frac{\rho I}{R} + r_1 +\frac{1}{2 } \sigma_4^2 + 2 \alpha \\ \leq -\frac{1}{4}\mu R^2 + A_7 \leq -\frac{\mu}{4} \frac{1}{\epsilon^4} + A_7, \end{align*} where \begin{equation} \label{A7a} \begin{split} A_7 =& \sup_{(S,V, I, R,X) \in \mathbb{R}_+^5} \left\lbrace M \left( \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) S + \frac{2 \beta^2 SX}{(\rho+\mu +\delta)(1+m_1S)} - \psi_1\right) \right.\\ & \left. -\frac{\gamma}{2 K} S^3 -\frac{1}{2}\mu V^2 - \frac{1}{2}(\mu+\delta)I^2 -\frac{1}{4}\mu R^2 -\frac{1}{2} X^2 + C + q \right\rbrace . \end{split} \end{equation} In view of \eqref{uD9}, it follows that $L W(S,V,I,R,X) \leq -1$ for $ (S,V,I,R,X) \in D_9$. Case 10. For any $(S,V, I, R,X) \in D_{10} $, we have \begin{align*} LW(S,V,I,R,X) %& \leq & M \left( \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) S + \frac{2 \beta^2 SX}{(\rho+\mu +\delta)(1+m_1S)} - \psi_1\right) -\frac{\gamma}{2 K} S^3 -\frac{1}{2}\mu V^2 - \frac{1}{2}(\mu+\delta)I^2 \\ %& & -\frac{1}{2}\mu R^2 -\frac{1}{2} X^2 + C -\frac{v S}{V} + \frac{\beta X}{1+m_2 V} +v +2 \mu +\frac{1}{2 } \sigma_2^2 -\frac{\rho I}{R} + r_1 +\frac{1}{2 } \sigma_4^2 + 2 \alpha \\ \leq -\frac{1}{4} X^2 + A_8 \leq -\frac{1}{4} \frac{1}{\epsilon^2} + A_8, \end{align*} where \begin{equation} \label{A8a} \begin{split} A_8 = & \sup_{(S,V, I, R,X) \in \mathbb{R}_+^5} \left\lbrace M \left( \left( \frac{\gamma }{ K } + \frac{2 \beta v}{\rho+\mu +\delta}\right) S + \frac{2 \beta^2 SX}{(\rho+\mu +\delta)(1+m_1S)} - \psi_1\right) \right.\\ & \left. -\frac{\gamma}{2 K} S^3 -\frac{1}{2}\mu V^2 - \frac{1}{2}(\mu+\delta)I^2 -\frac{1}{2}\mu R^2 -\frac{1}{4} X^2 + C +q \right\rbrace . \end{split} \end{equation} In view of \eqref{uD10}, it follows that $L W(S,V,I,R,X) \leq -1$ for $ (S,V,I,R,X) \in D_{10}$. On the basis of the above discussion, it follows that \begin{align*} LW(S,V,I,R,X) \leq -1 , \quad \text{for all} \quad (S, V,I,R,X) \in \mathbb{R}_+^5 \setminus D. \end{align*} In addition, it is satisfied \begin{align*} \hspace{-0.2cm} 0 \leq & \frac{1}{t} E {W}(S(t), V(t),I(t),R(t),X(t)) \\ = & \frac{1}{t} E {W}(S(0), V(0),I(0),R(0),X(0)) + \!\frac{1}{t} \!\int_{0}^{t} \!\!E \left( L {W}(S(s), V(s),I(s),R(s),X(s)) \right) ds \\ = & \liminf_{t \rightarrow \infty} \frac{1}{t} \int_{0}^{t} E \left( L {W}(S(s),V(s),I(s),R(s),X(s)) \right) ds. \\ = & \liminf\limits_{t \rightarrow \infty} \frac{1}{t} \int_{0}^{t} \!\!E \left(\! L {W}(S(s),\! V(s)\!,\!I(s)\!,R(s)\!,\!X(s)) \!\right) 1_{\left\lbrace (S(s),V(s),I(s),R(s),X(s)) \in D\right\rbrace } ds \\ & + \!\liminf\limits_{t \rightarrow \infty}\! \frac{1}{t}\! \int_{0}^{t} \!\!\!\! E \left(\! L {W}(S(s), \!V(s)\!,\!I(s),\!R(s)\!,\!X(s)) \!\right) 1_{\left\lbrace (S(s), V(s),I(s),R(s),X(s)) \in \mathbb{R}_+^5 \setminus D\right\rbrace } ds \\ \leq & K \liminf_{t \rightarrow \infty} \frac{1}{t} \int_{0}^{t} 1_{\left\lbrace (S(s), V(s),I(s),R(s),X(s)) \in D\right\rbrace } ds \\ & - \!\liminf_{t \rightarrow \infty} \! \!\frac{1}{t}\! \int_{0}^{t}\! 1_{\left\lbrace (S(s), V(s),I(s),R(s),X(s)) \in \mathbb{R}_+^5 \setminus D\right\rbrace } ds \\ \leq & \left( K+1\right) \liminf_{t \rightarrow \infty} \frac{1}{t} \int_{0}^{t} 1_{\left\lbrace (S(s),V(s),I(s),R(s),X(s)) \in D\right\rbrace } ds -1, \end{align*} where $K= \sup_{\left( S(s), V(s),I(s),R(s),X(s) \right) \in \mathbb{R}_+^5} L {W}(S(s), V(s),I(s),R(s),X(s)) $. Then \begin{align*} \liminf_{t \rightarrow \infty} \frac{1}{t} \int_{0}^{t} 1_{\left\lbrace (S(s),V(s),I(s),R(s),X(s)) \in D\right\rbrace } ds \geq \frac{1}{K+1} >0, \,\,\,\,\, \text{a.s.} \end{align*} Let ${P}(t, (S(t), V(t), I(t), R(t), X(t)), \Omega)$ be the transition probability that \\ $(S(t), V(t), I(t), R(t), X(t)) \in \Omega$. Making use of Fatou’s lemma, we have \begin{align*} \liminf_{t \rightarrow \infty} \frac{1}{t} \int_{0}^{t} {P}(s, (S(s), V(s), I(s), R(s), X(s)), D) ds \geq \frac{1}{K+1} >0, \,\,\,\,\, \text{a.s.} \end{align*} Regarding Lemma \ref{lemaA1}, we can conclude that model \eqref{SVIR stochastic} has a unique stationary distribution and ergodicity holds. The proof is completed. \end{proof} \section{Numerical simulations} In this section, some numerical simulations on the example of Ebola virus disease are presented in order to illustrate the results of theorems \ref{theorem2.A} and \ref{theorem3.A}. Nowadays, Ebola infection is still a healthcare problem in many parts of the world, namely in impoverished regions where it has devastating consequences. Its control is the goal of many studies in the last years. Table \ref{Model-estpar1a} gives realistic parameter values for the spread of Ebola virus disease used in simulations, taken from previous studies in which they were estimated using different statistical methods, derived on the basis of World Health Organization data (see \cite{WHO,d.Torres_ebola_2016,Din_Khan_Sabbar_2022,Mbah_2023} and references therein) or rationally assumed (a). \begin{table}[!ht] \centering %\ra{1.3} \begin{footnotesize} \begin{tabular}{@{}lrrrrrrrrr@{}} & $\gamma$ & $K$ & $\beta$ & $v_1$ & $r_1$ &$m_1$ & $m_2$ & $\rho$ \\ \hline \textbf{Extinction} & 0.0007 & 10000 & 0.0015 & 0.0011 & 0.0011 & 1 & 1.4 & 0.03 \\ \hline \textbf{Persistence}& 0.35 & 120 & 0.01 & 0.005345 & 0.0011 & 0.01 &0.1 & 0.1 \\ \hline Reference & \cite{WHO} & (a) & \cite{WHO} - \cite{Mbah_2023} & \cite{WHO} & \cite{WHO} & (a) & (a) & \cite{WHO} - \cite{Mbah_2023} \\ \hline \hline & $v$ & $\delta$ & $\mu$ & $\xi$ & & & & & \\ \hline \textbf{Extinction} & 0.0005 & 0.25 & 0.05 & 0.03 & & & & & \\ \hline \textbf{Persistence} & 0.07 & 0.45 & 0.1& 0.1 & & & & & \\ \hline Reference & (a) & \cite{WHO} - \cite{Mbah_2023} & \cite{WHO} & (a) & & & & & \\ \hline \hline \end{tabular} \end{footnotesize} \caption{\rm \em Model parameter values taken in extinction and persistence simulations.} \label{Model-estpar1a} \end{table} Figures \ref{Figure1A} and \ref{Figure2A} illustrate the trajectories of stochastic model \eqref{SVIR stochastic} and its deterministic counterpart (without noises), obtained with the parameters given in Table \ref{Model-estpar1a}. Parameters used in simulations and presented in Figure \ref{Figure1A} satisfy the condition of Theorem \ref{theorem2.A}. In Figure \ref{Figure1A}, we can clearly seed that the disease is eradicated a.s. which supports the theorem result. Here, it is assumed that noise intensities are $\sigma_1=0.001, \sigma_2=0.1, \sigma_3=0.2, \sigma_4=0.12, \sigma_5= 0.1, \eta_1=0.00025, \eta_2=-0.1, \eta_3=0.08, \eta_4=-0.15, \eta_5=0.05$ and initial conditions $S(0)=8000 , V(0)=200, I(0)=100, R(0)=500$. \begin{figure}[htb!] \includegraphics[height=1.9cm]{ext-S1a.eps} \includegraphics[height=1.9cm]{ext-V1a.eps} \includegraphics[height=1.9cm]{ext-I1a.eps} \includegraphics[height=1.9cm]{ext-R1a.eps} \caption{\rm \em Extinction of the disease. Deterministic and stochastic trajectories for susceptible, vaccinated, infected and recovered classes, respectively, of model \eqref{SVIR stochastic}. } \label{Figure1A} \end{figure} Figure \ref{Figure2A} illustrates the trajectories of the model \eqref{SVIR stochastic} obtained with the parameters that satisfy the condition of Theorems \ref{theorem3.A} and \ref{theorem4.A}. Therefore, by Theorems \ref{theorem3.A} and \ref{theorem4.A} we obtain that the disease is persistent with probability one and there exists a stationary distribution of system \eqref{SVIR stochastic}. The numerical simulations given in Figure \ref{Figure2A} confirm the validity of theoretical results. Here, it is assumed that noise intensities are $\sigma_1=0.03, \sigma_2=0.08, \sigma_3=0.05, \sigma_4=0.025, \sigma_5= 0.02, \eta_1=0.05, \eta_2=-0.1, \eta_3=0.03, \eta_4=-0.015, \eta_5=0.025$ and $S(0)=100 , V(0)=3, I(0)=3, R(0)=2$. \begin{figure}[htb!] \includegraphics[height=1.9cm]{pers-S1a.eps} \includegraphics[height=1.9cm]{pers-V1a.eps} \includegraphics[height=1.9cm]{pers-I1a.eps} \includegraphics[height=1.9cm]{pers-R1a.eps} \caption{\rm \em Persistence of the disease. Deterministic and stochastic trajectories for susceptible, vaccinated, infected and recovered classes, respectively, of model \eqref{SVIR stochastic}. } \label{Figure2A} \end{figure} \section{Conclusion} In this paper, a stochastic SVIR epidemic model is established, in which growth of susceptible individuals follows the logistic function, and saturated incidence rates and distributed delay denoting the incubation period are considered. It is assumed that the epidemic system is stochastically affected by mixture of environmental perturbations modeled by white noise and Poisson jumps. For the obtained model, sufficient conditions for the extinction and persistence in mean of the disease are obtained and the existence of ergodic stationary distribution is studied. Particularly, analysis showed that the disease is eradicated if the coefficient $R^* <1$ and the disease persists in the population and there exists ergodic stationary distribution if $R^S_0 > 1$. From the expression of $R^S_0$, we can conclude that white noise and Lévy jump have a positive effect in controlling the progression of the disease. That is, lower intensities of white noise and Lévy jump guarantee persistence of the disease. Some interesting topics deserve further consideration. One may consider more realistic but more complex models, for example, by adding the effects of telegraph noise to the proposed model or by introducing time-varying parameters. %%%% Acknowledgment %%%%%%%% \section*{Acknowledgement} The author would like to thank the editor and anonymous reviewers for their helpful and valuable comments and suggestions that significantly improved the paper. %%%% Bibliography %%%%%%%%%% \begin{thebibliography}{99} \bibitem{Allen_2008} {\sc L.\,J.\,S.\,Allen}, {\em An Introduction to Stochastic Epidemic Models}, in: {\em Mathematical Epidemiology}, (F.\,Brauer, P.\,van den Driessche, J.\,Wu (Eds.)), Lecture Notes in Mathematics, Vol {\bf 1945}, Springer, 2008. %Allen, L.J.S. (2008). An Introduction to Stochastic Epidemic Models. In: Brauer, F., van den Driessche, P., Wu, J. (eds) Mathematical Epidemiology. Lecture Notes in Mathematics, vol 1945. Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-540-78911-6_3 \bibitem{Binder_1999} {\sc S.\,Binder, A.\,Levitt, J.\,Sacks, J.\,Hughes}, {\em Emerging infectious diseases: Public health issues for the 21st century}, Science. {\bf 284}(1999), 1311--1313. \bibitem{Blower_1991} {\sc S.\,Blower, A.\,McLean}, {\em Mixing ecology and epidemiology}, Proc. R. Soc. Lond. B. {\bf 245}(1991), 187--192. \bibitem{Cao_Feng_2019} {\sc Z.\,Cao, W.\,Feng, X.\,Wen}, {\em Dynamical behavior of a stochastic SEI epidemic model with saturation incidence and logistic growth}, Physica A. {\bf 523}(2019), 894--907. \bibitem{Caraballo_Fatini_2020} {\sc T.\,Caraballo, M.\,E.\,Fatini, M.\,E.\,Khalifi}, {\em Analysis of a stochastic distributed delay epidemic model with relapse and gamma distribution kernel}, Chaos, Solitons and Fractals. {\bf 133}(2020), 109643. \bibitem{Cooke_Driessche_1996} {\sc K.\,L.\,Cooke, P.\,van Den Driessche}, {\em Analysis of an SEIRS epidemic model with two delays}, J. Math. Biol. {\bf 35}(1996), 240--260. \bibitem{Din_Khan_Sabbar_2022} {\sc A.\,Din, A.\,Khan, Y.\,Sabbar}, {\em Long-Term Bifurcation and Stochastic Optimal Control of a Triple-Delayed Ebola Virus Model with Vaccination and Quarantine Strategies}, Fractal and Fractional {\bf 6}(2022), 578. \bibitem{Din_2021} {\sc A.\,Din, Y.\,Li, T.\,Khan, K.\,Anwar, G.\,Zaman}, {\em Stochastic dynamics of hepatitis B epidemics}, Results Phys. {\bf 20}(2021), 103730. \bibitem{Zhang1} {\sc K.\,Fan, Y\,Zhang, S.\,Gao, S.\,Chen}, {\em A delayed vaccinated epidemic model with nonlinear incidence rate and Lévy jumps}, Physica A. {\bf 544}(2019), 123379, doi: doi.org/10.1016/j.physa.2019.123379 \bibitem{Ma_Hara_2002} {\sc T. Hara, E. Beretta, W. Ma, Y. Takeuchi}, {\em Permanence of an SIR epidemic model with distributed time delays}, Tohoku Math. J. {\bf 54}(2002), 581--591. \bibitem{Liu} {\sc T.\,Hayat, B.\,Ahmad, Q.\,Liu, D.\,Jiang}, {\em Analysis of a delayed vaccinated SIR epidemic model with temporary immunity and Lévy jumps}, Nonlinear Analysis: Hybrid Systems. {\bf 27}(2018), 29--43. \bibitem{HHethcote_2000} {\sc H.\,Hethcote}, {\em The mathematics of infectious diseases}, SIAM Rev. {\bf 42}(2000), 599--653. \bibitem{HHethcote_Driessche_1991} {\sc H.\,Hethcote, P.\,Van den Driessche}, {\em Some epidemiological models with nonlinear incidence}, J. Math. Biol. {\bf 29}(1991), 271--287. \bibitem{Miljana_Marija} {\sc M.\,Jovanovi\'{c}, M.\,Krsti\'{c}}, {\em Stochastically perturbed vector-borne disease models with direct transmission}, Appl. Math. Mod. {\bf 36}(2012), 5214--5228. \bibitem{Marija2018} {\sc M.\,Krsti\'{c}}, {\em On stability of stochastic delay model for tumor-immune interaction}, Filomat. {\bf 32}(2018), 1273--1283. \bibitem{Li_Mao_2009} {\sc X.\,Li, X.\,Mao}, {\em Population dynamical behavior of non-autonomous Lotka-Volterra competitive system with random perturbation}, Dis. Con. Dyn. Sys. A. {\bf 24}(2009), 523. \bibitem{Fatini17} {\sc X.\,Li, X.\,Mao}, {\em A stochastic SIRS epidemic model incorporating media coverage and driven by Lévy noise}, Chaos, Solitons and Fractals. {\bf 105}(2017), 60--68. \bibitem{Li_Mao_2022} {\sc D.\,Li, F.\,Wei, X.\,Mao}, {\em Stationary distribution and density function of a stochastic SVIR epidemic model}, Journal of the Franklin Institute. {\bf 359}(2022), 9422--9449. \bibitem{Liu_Jiang_2017} {\sc Q.\,Liu, D.\,Jiang, N.\,Shi, T.\,Hayat, A.\,Alsaedi}, {\em Dynamical behavior of a stochastic HBV infection model with logistic hepatocyte growth}, Acta Math. Sci. {\bf 37}(2017), 927--940. \bibitem{Macdonald_1978} {\sc N.\,MacDonald}, {\em Time Lags in Biological Models, in: Lecture Notes in Biomathematics}, Springer-Verlag, Heidelberg, 1978. \bibitem{Mao} {\sc X.\,Mao}, {\em Stochastic Differential Equations and Applications}, Woodhead Publishing Limited, 2011. \bibitem{Mao_Marion_2002} {\sc X.\,Mao, G.\,Marion, E.\,Renshaw}, {\em Environmental Brownian noise suppresses explosions in population dynamics}, Stoch. Pro. Appl. {\bf 97}(2002), 95--110. \bibitem{May_Anderson_1978} {\sc R.\,M.\,May, R.\,M.\,Anderson}, {\em Regulation and stability of host-parasite population interactions II: Destabilizing process}, J. Anim. Ecol. {\bf 47}(1978), 219--267. \bibitem{Mbah_2023} {\sc G.\,C.\,E.\,Mbah, I.\,S.\,Onah, Q.\,O.\,Ahman, O.\,C.\,Collins, C.\,C.\,Asogwa, C.\,Okoye}, {\em Mathematical modeling approach of the study of Ebola virus disease transmission dynamics in a developing country}, Afr. J. Infect. Dis. {\bf 17}(2023), 10--26. \bibitem{Meyn1993} {\sc S.\,P.\,Meyn, R.\,L.\,Tweedie}, {\em Stability of Markovian processes III: Foster-Lyapunov criteria for continuous-time processes}, Adv. Appl. Probab. {\bf 25}(1993), 518--548. \bibitem{Aksendal_2010} {\sc B.\,{\O}ksendal}, {\em Stochastic differential equations: An introduction with applications}, Springer, 2010. \bibitem{Oksendal_2005} {\sc B.\,Øksendal, A.\,Sulem}, {\em Applied Stochastic Control of Jump Diffusions}, Heidelberg: Springer-Verlag, 2005. \bibitem{d.Torres_ebola_2016} {\sc A.\,Rachah, D.\,F.\,M.\,Torres}, {\em Mathematical modelling, simulation, and optimal control of the 2014 Ebola outbreak in west Africa}, Discrete Dynamics in Nature and Society. {\bf 842792}(2015), 9 pages \bibitem{Ruan_Wang_2003} {\sc S. Ruan, W. Wang}, {\em Dynamical behavior of an epidemic model with nonlinear incidence rate}, J. Differential Equations. {\bf 188}(2003), 135--163. \bibitem{Situ} {\sc R.\,Situ}, {\em Theory of stochastic differential equations with jumps and applications}, Springer, 2005. \bibitem{Stettner} {\sc L.\,Stettner}, {\em On the existence and uniqueness of invariant measure for continuous time Markov processes} Technical Report, LCDS 86-18, Brown University, Providence, RI, 1986. \bibitem{WWang_2002} {\sc W.\,Wang}, {\em Global behavior of an SEIRS epidemic model with time delays}, Appl. Math. Lett. {\bf 15}(2002), 423--428. \bibitem{WHO} {\sc World Health Organization}, https://www.who.int/health-topics/ebola, 2023. \bibitem{Xi_Zhu_2017} {\sc F.\,Xi, C.\,Zhu}, {\em On Feller and strong Feller properties and exponential ergodicity of regime-switching jump diffusion processes with countable regimes}, SIAM J. Control Optim. {\bf 55}(2017), 1789–1818. \bibitem{Zhang16} {\sc X.\,Zhang, F.\,Chen, K.\,Wang, H.\,Du}, {\em Stochastic SIRS model driven by Lévy noise}, Acta Math. Sci. {\bf 36B}(2016), 740--752. \bibitem{Zhou_Yuan_2016} {\sc Y.\,Zhou, S.\,Yuan, D.\,Zhao}, {\em Threshold behavior of a stochastic SIS model with L\'evy jumps}, Applied Mathematics and Computation {\bf 275}(2016), 255–267. \end{thebibliography} %\end{linenumbers} \end{document}