Introduction

Recently, global social economics and human health are greatly affected by the infectious disease. The awareness has been increased for preventing and to control the world-wide spreading of COVID-19. Mathematical model is one of the most important tools to describe the behavior of the epidemic in epidemiology. Kermack and Mckendric (1927) investigated a susceptible-infected-susceptible (SIS) model and discussed its dynamical properties. The transmissions of various epidemics (Liu et al. 2008; Li et al. 2017, 2001; Jerubet and Kimathi 2019; Hove-Musekwa and Nyabadza 2009; Iwami et al. 2007; Cai and Wu 2009; Vincenzo and Gabriella 1978; Carter et al. 2020; Mahato et al. 2021; Das et al. 2021) were developed with realistic ordinary differential equation. (Kuniya and Wang 2018) established a susceptible-infected-recovered (SIR) epidemic model. They studied the existence of the global stability of disease-free equilibrium for basic reproduction number \({R}_{0}<1\) and investigated the uniform persistence of the system under \({R}_{0}>1.\) (Naik et al. 2020a) formulated a SIR epidemic model with Crowley–Martin incidence rate and Holling type–II treatment. With the help of La-Salle invariance principle and the Lyapunov function, they examined the existence and stability of both equilibrium points. A nonlinear fractional infectious disease model in HIV transmission was discussed by Naik et al. (2020b). A SIRS epidemic model with mixed vaccination strategies was established by Gao et al. (2011). They added the seasonal variability on infection in their study. A SVIR (Susceptible-Vaccinated-Infected-Recovered) epidemic model was studied by Liu et al. (2008). They discussed the vaccine strategies for controlling the disease. Xing and Li (2021) established a SVIR epidemic model with relapse and investigated the persistence of the epidemic. Lee and Lao (2018) described an epidemic model with bilinear incidence rate and investigated the transmission dynamics of the disease.

Environmental variations (Liu et al. 2018; Cai and Kang 2015) may perturb the travel of population, design of control strategies. Many researchers formulated the stochastic differential equation for epidemic and analyzed their dynamical behavior. A stochastic SIS epidemic model with vaccination was studied by Zhao and Jiang (2014). The theory of extinction and persist in mean of the epidemic was investigated in their work. In Caraballo et al. (2020), Carballo et al. formulated a stochastic SIRS epidemic model and discussed the condition of stationary distribution. The existence of ergodic stationary distribution and the probability density function of the SVIS epidemic model were studied by Zhou et al. (2020). Zhou et al. (2021) solved the general three-dimensional Fokker–Planck equation. The existence of stationary distribution and ergodicity of the system has been discussed. They have studied the impact of random noises on the disease extinction.

Concentration in the vaccination on infected individuals and random oscillation, the main theme of our study is to improve SIVS epidemic model in stochastic nature with vaccination strategies. The disease which will be persistent depends on the corresponding basic reproduction number of the system. There is no positive equilibrium exists in the stochastic system for the environmental fluctuation. So, the stochastic permanence of the epidemic can be greatly affected by the existence of a stationary distribution and the properties of ergodicity. For controlling the outbreak of epidemic, we must need for statistical data of the disease in our real life. For difficulties to solve the higher-order Fokker–Planck equations, we have analyzed the probability density function. In this work, some studies of probability density function of stationary distribution are discussed. For this, we focus for three points.

  1. (i)

    Comprise stochastic threshold \({R}_{S}^{*}\) with respect to the basic reproduction number \({R}^{*}.\)

  2. (ii)

    Look into the persistence of the disease of stochastic SVIS system under the condition of \({R}_{S}^{*}>1.\)

  3. (iii)

    Discuss the numerical simulation and the sensitivities of the ecological parameters of the system get a clear view of our study.

This work is represented as follows: Sect. "Model calibration and dynamical behavior:" presents model formulation and necessary notation for the model. Sect. "Persistence and extinction of the system" introduces the persistence and the extinction of the stochastic system. The ergodic stationary distribution under \({R}_{S}^{*}>1\) is investigated in the subsection of 3. Sect. "Numerical results:" shows some numerical results and Sect. "Sensitivities of the parameters:" represents the sensitivities of the parameters. Finally, results are discussed and conclusions are drawn in Sect. "Discussion of results."

Model calibration and dynamical behavior

In this section, the deterministic and stochastic SVIS epidemic models are developed after considering some suitable conditions.

Deterministic SVIS epidemic system

Suppose, \(N(t)\) is the total investigated population. This population is divided into susceptible \(S\left(t\right),\) vaccinated \(V\left(t\right),\) infected \(I(t)\) populations at any instant \(t.\) In this study, the susceptible individuals obey the rule of logistic growth model which is the growth process of species in natural way (Jiang et al. 2007; Arino et al. 2006; Xu et al. 2015). We formulate a deterministic SVIS epidemic system (Zhou et al. 2020, 2021) with saturated incidence and vaccination strategies, which is given as

$$\frac{dS\left( t \right)}{{dt}} = pS\left( t \right)\left( {1 - \frac{S\left( t \right)}{q}} \right) - \frac{\lambda S\left( t \right)I\left( t \right)}{{1 + k_{1} I\left( t \right)}} - \psi S\left( t \right) - \alpha_{1} S\left( t \right) + \beta_{2} I\left( t \right) + \beta_{1} V\left( t \right)$$
$$\frac{dV\left( t \right)}{{dt}} = \alpha_{1} S\left( t \right) - \left( {\beta_{1} + \psi } \right)V\left( t \right)$$
$$\frac{dI\left( t \right)}{{dt}} = \frac{\lambda S\left( t \right)I\left( t \right)}{{1 + k_{1} I\left( t \right)}} - \left( {\psi + \alpha_{2} + \beta_{2} } \right)I\left( t \right)$$
(1)

where \(p,q,\lambda , {k}_{1}, \psi , {\alpha }_{1},{\beta }_{2},{\beta }_{1}\) and \({\alpha }_{2}\) are all positive constants. The parameter \(p\) is the intrinsic growth rate of susceptible individuals, and \(q\) depicts the carrying capacity of \(S\left(t\right).\) The term \(\frac{\lambda S\left(t\right)I\left(t\right)}{1+{k}_{1}I\left(t\right)}\) is more sensitive than the bilinear incidence rate in epidemiological study (Zhu et al. 2020; Xu et al. 2016; Batabyal and Batabyal 2021; Chong et al. 2014). Here, \(\lambda\) denotes the disease transmission rate between susceptible and infected individuals, and \({k}_{1}\) is the half-saturation constant. The parameters \(\psi\) and \({\alpha }_{1}\) represent the natural mortality rate of all individuals and vaccination rate of susceptible individuals, respectively. The disease-related death rate of the infected individuals is represented by \({\alpha }_{2}\). The immunity loss coefficient of the vaccinated individuals and the recovery rate of the infected individuals are depicted by the parameters \({\beta }_{1}\) and \({\beta }_{2}\), respectively. Ecological description of the parameters and its units are shown in Appendix 1.

The disease-free equilibrium point of the system is given by

$$E^{\# } = \left( {S^{\# } , V^{\# } , I^{\# } } \right) = \left[ {\frac{q}{p}\left\{ {\frac{{\left( {p - \psi - \alpha_{1} } \right)\left( {\beta_{1} + \psi } \right) + \beta_{1} \alpha_{1} }}{{\beta_{1} + \psi }}} \right\},\frac{{\alpha_{1} }}{{(\beta_{1} + \psi )}}\frac{q}{p}\left\{ {\frac{{\left( {p - \psi - \alpha_{1} } \right)\left( {\beta_{1} + \psi } \right) + \beta_{1} \alpha_{1} }}{{\beta_{1} + \psi }}} \right\},0} \right].$$

Now, we compute the basic reproduction number (Driessche and Watmough 2002) of the deterministic system.

\(R^{*} = \frac{{\lambda q\left\{ {\left( {p - \psi - \alpha_{1} } \right)\left( {\beta_{1} + \psi } \right) + \beta_{1} \alpha_{1} } \right\}}}{{p(\beta_{1} + \psi )\left( {\psi + \alpha_{2} + \beta_{2} } \right)}}\) (see Appendix 2). Actually, the basic reproduction number \(R^{*}\) can be defined as the number of new infections started from infective individuals at disease-free equilibrium. \(R^{*} < 1\) indicates infected population creates less than one new infected population at time of its infective situation and the disease becomes extinct. In other way, \(R^{*} > 1\) represents that each infected populations create more than one new infection, and then, the infection can spread over the population.

For the endemic equilibrium point of the system (1), we obtain.

\(E^{@} = \left( {S^{@} , V^{@} ,I^{@} } \right)\) and \(S^{@} = \frac{{(\beta_{1} + \psi )V^{@} }}{{\alpha_{1} }}\), \(I^{@} = \frac{1}{{\psi + \alpha_{2} }}\left[ {\frac{{(\beta_{1} + \psi )}}{{\alpha_{1} }}\left\{ {p\left( {1 - \frac{{k_{1} }}{q}} \right) - \left( {\psi + \alpha_{1} } \right)} \right\} + \beta_{1} } \right]V^{@}\). where \(V^{@}\) is the root of the following equation

\(E_{1} \left( {V^{@} } \right)^{2} + E_{2} V^{@} + E_{3} = 0\) and \(E_{1} = \frac{{\left( {\psi + \alpha_{2} + \beta_{2} } \right)p}}{{\alpha_{1} \left( {\psi + \alpha_{2} } \right)q}}\), \(E_{2} = \frac{\lambda }{{\alpha (\beta_{1} + \psi )}} - \frac{{p\left( {\psi + \alpha_{2} + \beta_{2} } \right)}}{{\alpha_{1} \left( {\psi + \alpha_{2} } \right)\left( {\beta_{1} + \psi } \right)}}\),

$$E_{3} = \frac{{\psi + \alpha_{2} + \beta_{2} }}{{(\beta_{1} + \psi )^{2} }}.$$

Stochastic SVIS epidemic system

In real life, the dynamical properties of the maximum epidemiological model are greatly influenced by random perturbation in the environments. With the help of relevant study (Cai and Kang 2015; Zhao and Jiang 2014; Khan and Khan 2018; Zhang 2017; Caraballo et al. 2020; Wang and Jiang 2019; Wang and Wang 2018; Liu et al. 2019a; Zhou et al. 2020), we consider the stochastic perturbations are directly proportional to \(S\left( t \right),V\left( t \right)\), \(I\left( t \right)\). Stochastic perturbations are influenced by multiplicative noises. These multiplicative noises are considered for describing the non-equilibrium systems to better understand the fluctuations which are not self-originating. So, the corresponding stochastic SVIS epidemic system with vaccination strategies is formulated by

$$dS\left( t \right) = \left[ {pS\left( t \right)\left( {1 - \frac{S\left( t \right)}{q}} \right) - \frac{\lambda S\left( t \right)I\left( t \right)}{{1 + k_{1} I\left( t \right)}} - \psi S\left( t \right) - \alpha_{1} S\left( t \right) + \beta_{2} I\left( t \right) + \beta_{1} V\left( t \right)} \right]dt +\tau_{1} SdG_{1} \left( t \right)$$
$$dV\left( t \right) = \left[ {\alpha_{1} S\left( t \right) - \left( {\beta_{1} + \psi } \right)V\left( t \right)} \right]dt +\tau_{2} VdG_{2} \left( t \right)$$
$$dI\left( t \right) = \left[ {\frac{\lambda S\left( t \right)I\left( t \right)}{{1 + k_{1} I\left( t \right)}} - \left( {\psi + \alpha_{2} + \beta_{2} } \right)I\left( t \right)} \right]dt +\tau_{3} IdG_{3} \left( t \right)$$
(2)

Here, \(G_{i} \left( t \right)\left( {i = 1,2,3} \right)\) are the independent standard Brownian motions and \(\tau_{i} \left( {i = 1,2,3} \right)\) are their intensities.

Existence and uniqueness

Now, we have investigated the existence and uniqueness of the solutions.

Theorem 1

The solutions \(\left\{ {S\left( t \right), V\left( t \right), I\left( t \right)} \right\} \in R_{3}^{ + }\) of the stochastic system (2) are unique and exist in the region \(R_{3}^{ + }\) for the initial values \(\left\{ {S\left( 0 \right), I\left( 0 \right), V\left( 0 \right)} \right\} \in R_{3}^{ + } .\)

Proof

The coefficients of the system (2) are locally Lipschitz continuous as reported by of Mao et al. (Mao 1997). \(\left\{ {S\left( t \right), V\left( t \right), I\left( t \right)} \right\} \in R_{3}^{ + } , t \in \left( {0, \sigma_{0} } \right)\) is the solution of the system (2) for the initial values \(\left\{ {S\left( 0 \right), I\left( 0 \right), V\left( 0 \right)} \right\} \in R_{3}^{ + }\). \(\sigma_{0}\) depicts the explosion time. For proving the global solution, we must prove \(\sigma_{0} = \infty\) a.s. The stopping time is satisfied the following conditions.

\(\sigma_{n} = {\text{inf}}\{ 0 < t <\sigma_{0} :{\text{min}}\left\{ {S\left( t \right), V\left( t \right), I\left( t \right)} \right\} \le \frac{1}{\delta }\) or max \(\left\{ {S\left( t \right), V\left( t \right), I\left( t \right)} \right\} \ge \delta \}\).where, \(\left\{ {S\left( 0 \right), I\left( 0 \right), V\left( 0 \right)} \right\} \in \left[ {\frac{1}{{\delta_{0} }},\delta_{0} } \right]\) and \(\delta \ge \delta_{0} .\)

We also assume the set inf{Φ} = \(\infty\), when \(n \to \infty\) then \(\sigma_{n}\) is increasing. Thus, we get \(\sigma_{\infty } = \mathop {\lim }\limits_{n \to \infty}\sigma_{n} .\) If we show that \(\sigma_{\infty } = \infty\) almost surely, then \(\sigma_{0} = \infty\) almost surely. This implies that \(\left\{ {S\left( t \right), V\left( t \right), I\left( t \right)} \right\} \in R_{3}^{ + }\) almost surely for \(t \ge 0.\) When \(\sigma_{\infty } < \infty\), then there exist two positive constants \(M\) and \(\varepsilon\) such that \(P\left\{ {\sigma_{\infty } \le M} \right\} > \varepsilon\). Again, we define.

\(P\left\{ {\sigma_{n} \le M} \right\} \ge \varepsilon\) for any integer \(\delta_{1} \ge \delta_{0} .\)

Let us assume a fundamental function \(\widehat{F:}R_{3}^{ + } \to R\) such that \(\widehat{F }\left( {S, V, I} \right) = \left( {S - \text{In}S - 1} \right) + \left( {V - \text{InV} - 1} \right) + \left( {I - 1 - \text{In}I} \right).\) As \(w > 0\), then \(\left( {w - 1 - \text{In}w} \right) \ge 0\). Thus \(\widehat{F }\left( {S, V, I} \right)\) is non-negative \(C^{2}\) function.

Using It \({\hat{\text{o}}}\)’s formula, we obtain

$${\text{d}}\widehat{F }\left( {S, V, I} \right) = {\mathbb{L}}\widehat{F }\left( {S, V, I} \right){\text{d}}t +\tau_{1} \left( {S - 1} \right){\text{d}}G_{1} \left( t \right) +\tau_{2} \left( {V - 1} \right){\text{d}}G_{2} \left( t \right) +\tau_{3} \left( {I - 1} \right){\text{d}}G_{3} \left( t \right)$$

where \({\mathbb{L}}\widehat{{F{ }}}\left( {S,{ }V,{ }I} \right) = \left( {1 - \frac{1}{S}} \right)\left[ {pS\left( {1 - \frac{S}{q}} \right) - \frac{\lambda SI}{{1 + k_{1} I}} - \psi S - \alpha_{1} S + \beta_{2} I + \beta_{1} V} \right] + \frac{{\tau_{1}^{2} }}{2} + \left( {1 - \frac{1}{V}} \right)\left[ {\alpha_{1} S - \left( {\beta_{1} + \psi } \right)V} \right] + \frac{{\tau_{2}^{2} }}{2} + \left( {1 - \frac{1}{I}} \right)\left[ {\frac{\lambda SI}{{1 + k_{1} I}} - \left( {\psi + \alpha_{2} + \beta_{2} } \right)I} \right] + \frac{{\tau_{3}^{2} }}{2}\)

$$\le p + 3\psi + \alpha_{1} + \alpha_{2} + \beta_{1} + \beta_{2} + \frac{{\tau_{1}^{2} +\tau_{2}^{2} +\tau_{3}^{2} }}{2} - p\left( {1 - \frac{S}{q}} \right) + \frac{\lambda I}{{1 + k_{1} I}}$$
$$\le p + 3\psi + \alpha_{1} + \alpha_{2} + \beta_{1} + \beta_{2} + \frac{{\tau_{1}^{2} +\tau_{2}^{2} +\tau_{3}^{2} }}{2} \le \hat{H}.$$

Therefore, we get. \({\text{d}}\widehat{F }\left( {S, V, I} \right) \le \hat{H}{\text{d}}t + \left[ {\tau_{1} \left( {S - 1} \right){\text{d}}G_{1} \left( t \right) +\tau_{2} \left( {V - 1} \right){\text{d}}G_{2} \left( t \right) +\tau_{3} \left( {I - 1} \right){\text{d}}G_{3} \left( t \right)} \right]\).

We integrate both sides with respect to limit runs from \(0\) to \(\sigma_{n}\wedge{\text{ M }}\) and take expectation,

$${\mathbb{E}}\hat{F}\left( {S\left( {\sigma_{n}\wedge M} \right), V\left( {\sigma_{n}\wedge M} \right), I\left( {\sigma_{n}\wedge {\text{M}}} \right)} \right) \le \hat{H}{\mathbb{E}}\left( {\sigma_{n}\wedge {\text{ M }}} \right) + \hat{F}\left( {S\left( 0 \right),{ }V\left( 0 \right),{ }I\left( 0 \right)} \right) \le \hat{H}M + \hat{F}\left( {S\left( 0 \right),{ }V\left( 0 \right),{ }I\left( 0 \right)} \right).$$

Let us consider \(\Omega_{m} = \left\{ {\sigma_{n} \le M} \right\}\). for \(m \ge m_{1}\)., then \(P\left( {\Omega_{m} } \right) \ge \varepsilon .\).

For every \(\eta\in \Omega_{m} ,\) such that \(S\left( {\sigma_{n} ,\eta} \right), V\left( {\sigma_{n} ,\eta} \right), I\left( {\sigma_{n} ,\eta} \right)\) is equal to \(\frac{1}{m}\) or \(m.\)

So, \(\hat{F}\left( {S\left( {\sigma_{m} ,\eta} \right), V\left( {\sigma_{m} ,\eta} \right), I\left( {\sigma_{m} ,\eta} \right)} \right)\) is not less than either \(m - 1 - {\text{ln}}m\) or \(\frac{1}{m} - 1 + {\text{ln}}m.\)

$$\hat{H}M + \hat{F}\left( {S\left( 0 \right), V\left( 0 \right), I\left( 0 \right)} \right) \ge {\mathbb{E}}[T_{{\Omega_{m} }} (\eta)\hat{F}\left( {S\left( {\sigma_{m} ,\eta} \right), V\left( {\sigma_{m} ,\eta} \right), I\left( {\sigma_{m} ,\eta} \right)} \right)]$$
$$\ge \varepsilon \left\{ {m - 1 - {\text{ln}}}m \right\} \wedge\left\{ {\frac{1}{m} - 1 + {\text{ln}}}m \right\}$$

Here, \(T_{{\Omega_{m} }}\) represents the indicator function of \(\Omega_{m} .\) Applying \(m \to \infty\) both sides, that gives a contradiction

$$\infty = \hat{H}M + \hat{F}\left( {S\left( 0 \right), V\left( 0 \right), I\left( 0 \right)} \right) < \infty .$$

This gives us \(\sigma_{\infty } = \infty\) a.s. Hence, the proof is completed.

Persistence and extinction of the system

In this section, we have discussed the persistence and the extinction of the system. For this purpose, we define that the stochastic reproductive ratio of the system (2) is \(R_{S}^{*} = \frac{{\lambda q\left\{ {\alpha_{1} \beta_{1} + \left( {\beta_{1} + \psi } \right)\left( {p - \psi - \alpha_{1} } \right)} \right\}}}{{\left( {p + \frac{{\tau_{1}^{2} }}{2}} \right)\left( {\beta_{1} + \psi + \frac{{\tau_{2}^{2} }}{2}} \right)\left( {\psi + \alpha_{2} + \beta_{2} + \frac{{\tau_{3}^{2} }}{2}} \right)}}\).

Stationary Distribution and ergodic property of the Stochastic System

Since the stochastic system (2) has no endemic equilibrium point, we examine the ergodic stationary distribution that represents the persistence of the disease. Now, we study some important lemma of the theory of Khas’minskii (Khas’miniskii RZ 1980).

Let us assume a stochastic differential equation

$${\text{d}}Y\left( t \right) = g\left( y \right){\text{d}}t + \mathop \sum \limits_{k = 1}^{n} h_{k} \left( y \right){\text{d}}G_{k} \left( t \right)$$

where \(Y\left( t \right)\) is a homogeneous Markov process in d-dimensional Euclidean space \(\gamma_{d} .\) The diffusion matrix is \(B\left( Y \right) = \left( {b_{ij} \left( y \right)} \right);\) and \(b_{ij} \left( y \right) = \mathop \sum \limits_{l = 1}^{n} h_{l}^{\left( i \right)} \left( y \right)h_{l}^{\left( j \right)} \left( y \right).\)

Lemma 1

Zhou et al. 2021; Khas’miniskii RZ 1980): The Markov process \(Y\left( t \right)\) has a unique ergodic stationary distribution \(\pi \left( . \right)\) for any bounded region \(R\) with boundary \(\pi\) and.

  1. (i)

    There is a non-negative integer \(N\) such that \(\mathop \sum \limits_{i,j = 1}^{d} b_{ij} \left( y \right)\kappa_{i} \kappa_{j} \ge N\left| k \right|^{2}\), \(y \in R, \kappa \in R^{d} .\)

  2. (ii)

    There is a non-negative \(C^{2}\) function \(F\left( {Y\left( t \right)} \right)\) such that \(LF\left( {Y\left( t \right)} \right)\) is negative for \(y \in \Omega_{R} \backslash R\). Then, for all \(y \in R^{d}\) and integral function \(\chi \left( . \right)\) with respect to the measure \(\chi \left( . \right)\), it follows that

\(P\left\{ {\mathop {\lim }\limits_{q \to \infty } \frac{1}{q}\mathop \smallint \limits_{0}^{q} \chi (y\left( t \right)) = \smallint \chi \left( y \right)\pi \left( {{\text{d}}y} \right)} \right\} = 1\).

Theorem 2

The stochastic system (2) has ergodic property and a unique stationary distribution \(\pi \left( . \right)\) for the initial value \(\left\{ {S\left( 0 \right), V\left( 0 \right), I\left( 0 \right)} \right\} \in R_{3}^{ + }\) and \(R_{S}^{*} > 1.\)

Proof

In the previous theorem, we prove the unique global positive solution for the initial values \(\left\{ {S\left( 0 \right), V\left( 0 \right), I\left( 0 \right)} \right\} \in R_{3}^{ + } .\) For proving this theorem, we only verify Lemma 1.

We choose a \(C^{2}\) function \(Z\) such that

$$Z\left( {S, V, I} \right) = N_{0} \left( {S + V + I - b_{1} {\text{ln}}S - b_{1}b_{2} {\text{ln}}V - b_{3} \text{ln}I} \right) - {\text{ln}}S - {\text{ln}}V + \left( {S + V + I} \right)$$

where \(b_{1} , b_{2} ,b_{3}\) all are positive constants and \(N_{0} = \frac{{p\left( {1 - \frac{S}{q}} \right) + 2\psi + \beta_{1} + \alpha_{1} + \frac{{\tau_{1}^{2} +\tau_{2}^{2} }}{2} + 2}}{{3p\left( {\sqrt[3]{{R_{S}^{*} }} - 1} \right)\left( {1 - \frac{S}{q}} \right)}} > 0.\)

From above expression, we have

$$- 3N_{0} p\left( {1 - \frac{S}{q}} \right)\left( {\sqrt[3]{{R_{S}^{*} }} - 1} \right) + p\left( {1 - \frac{S}{q}} \right) + 2\psi + \beta_{1} + \alpha_{1} + \frac{{\tau_{1}^{2} +\tau_{2}^{2} }}{2} = - 2$$

For simplification,

$$Z_{1} = \left( {S + V + I - b_{1} {\text{ln}}S - b_{1} b_{2} {\text{ln}}V - b_{3} {\text{ln}}I} \right)$$
$$Z_{2} = - {\text{ln}}S - {\text{ln}}V; Z_{3} = S + V + I$$

Using It \(\hat{o}^{\prime}\) s formula to \(Z_{1}\) (see Appendix 3) we derive

$$LZ_{1} = p\left( {1 - \frac{S}{q}} \right) - \psi \left( {S + V + I} \right) - \alpha_{2} I - b_{1} \left[ {p\left( {1 - \frac{S}{q}} \right) - \frac{\lambda I}{{1 + k_{1} I}} + \beta_{2} \frac{I}{S} + \beta_{1} \frac{V}{S} - \left( {\psi + \alpha_{1} + \frac{{\tau_{1}^{2} }}{2}} \right)} \right] - b_{1} b_{2} \left[ {\alpha_{1} \frac{S}{V} - \left( {\beta_{1} + \psi + \frac{{\tau_{2}^{2} }}{2}} \right)} \right] - b_{3} \left[ {\frac{\lambda S}{{1 + k_{1} I}} - \left( {\psi + \alpha_{2} + \beta_{2} + \frac{{\tau_{3}^{2} }}{2}} \right)} \right]$$
$$\le p\left( {1 - \frac{S}{q}} \right) - \left[ {\psi N + b_{1} p\left( {1 - \frac{S}{q}} \right) + \frac{{b_{3} \lambda S}}{{1 + k_{1} I}}} \right] + b_{1} \left( {\psi + \alpha_{1} + \frac{{\tau_{1}^{2} }}{2}} \right) + b_{3} \left( {\psi + \alpha_{2} + \beta_{2} + \frac{{\tau_{3}^{2} }}{2}} \right) - \left( {\frac{{b_{1} \beta_{2} I}}{S} + \frac{{b_{1} b_{2} }}{V}\alpha_{1} S} \right) + b_{1} b_{2} \left( {\beta_{1} + \psi + \frac{{\tau_{2}^{2} }}{2}} \right) + \frac{{b_{1} \lambda I}}{{1 + k_{1} I}}$$
$$\le p\left( {1 - \frac{S}{q}} \right) - 3\left\{ {\lambda b_{1} b_{3} p\left( {1 - \frac{S}{q}} \right)} \right\}^{\frac{1}{3}} - 2(b_{1}^{2} b_{2} \beta_{1} \psi )^{1/2} + b_{1} \left( {\psi + \alpha_{1} + \frac{{\tau_{1}^{2} }}{2}} \right) + b_{3} \left( {\psi + \alpha_{2} + \beta_{2} + \frac{{\tau_{3}^{2} }}{2}} \right) + \frac{{b_{1} \lambda I}}{{1 + k_{1} I}}$$

where \(b_{1} , b_{2} ,b_{3}\) represents

$$b_{2} \left( {\psi + \alpha_{1} + \frac{{\tau_{1}^{2} }}{2}} \right) = \lambda \alpha_{1} ,$$
$$b_{1} \left[ {\psi + \alpha_{2} + \beta_{2} + \frac{{\tau_{3}^{2} }}{2} - \frac{{\lambda \alpha_{1} }}{{\psi + \alpha_{1} + \frac{{\tau_{1}^{2} }}{2}}}} \right] = b_{3} \left( {\psi + \alpha_{2} + \beta_{2} + \frac{{\tau_{3}^{2} }}{2}} \right) = p\left( {1 - \frac{S}{q}} \right)$$

Above expressions imply the following values.\(b_{1} = \frac{{\lambda \alpha_{1} }}{{\left( {\psi + \alpha_{2} + \frac{{\tau_{2}^{2} }}{2}} \right)^{2} }}, b_{2} = \frac{{p\left( {1 - \frac{S}{q}} \right)}}{{\psi + \alpha_{1} + \frac{{\tau_{1}^{2} }}{2} - \frac{{\lambda \alpha_{1} }}{{\psi + \alpha_{1} + \frac{{\tau_{2}^{2} }}{2}}}}}\) and \(b_{3} = \frac{{p\left( {1 - \frac{S}{q}} \right)}}{{\psi + \alpha_{2} + \beta_{2} + \frac{{\tau_{3}^{2} }}{2}}}\)

Therefore,

$$LZ_{1} \le 3p\left( {1 - \frac{S}{q}} \right) - 3p\left( {1 - \frac{S}{q}} \right)\left( {R_{S}^{*} } \right)^{1/3} + \frac{{b_{1} \lambda I}}{{1 + k_{1} I}}$$
(3)

Next, we use It \(\hat{o}^{\prime}s\) formula on \(Z_{2}\) and \(Z_{3}\). We obtain

$$LZ_{2} = \left[ { - p\left( {1 - \frac{S}{q}} \right) + \left( {\psi + \alpha_{1} + \frac{{\tau_{1}^{2} }}{2}} \right) + \frac{\lambda I}{{1 + k_{1} I}} - \frac{{\beta_{1} V}}{S} - \frac{{\beta_{2} I}}{S}} \right] - \frac{{\alpha_{1} S}}{V} + \left( {\psi + \alpha_{1} + \frac{{\tau_{1}^{2} }}{2}} \right)$$
$$\le - p\left( {1 - \frac{S}{q}} \right) - \frac{{\beta_{1} V}}{S} - \frac{{\beta_{2} I}}{S} + 2\psi + \alpha_{1} + \frac{{\tau_{1}^{2} +\tau_{2}^{2} }}{2}$$
(4)
$$LZ_{3} = p\left( {1 - \frac{S}{q}} \right) - \psi \left( {S + I} \right) - \left( {\psi + \alpha_{2} } \right) \le p\left( {1 - \frac{S}{q}} \right) - \psi N - \frac{{\alpha_{1} S}}{V}$$
(5)

Since \(Z\left( {S, V,I} \right)\) is continuous function,

$$\mathop {\lim }\limits_{t \to + \infty } {\text{inf}} Z\left( {S, V, I} \right) = \infty , \left( {S, V, I} \right) \in R_{ + }^{3} \backslash E_{\mu } .$$

Now, we obtain.

\(Q\left( {S, V, I} \right) = Z\left( {S, V, I} \right) - Z_{0}\), where \(Z_{0}\) is minimum value of \(Z_{1}\) and \(Q\left( {S, V, I} \right) \in R_{ + }^{3}\).

Using (2, 3, 4, 5), we get

$$LQ \le - 3N_{0} p\left( {1 - \frac{S}{q}} \right)\left\{ {\left( {R_{S}^{*} } \right)^{\frac{1}{3}} - 1} \right\} + p\left( {1 - \frac{S}{q}} \right) - \frac{{\alpha_{1} S}}{V} + \frac{\lambda I}{{1 + k_{1} I}} + \left\{ {p\left( {1 - \frac{S}{q}} \right) + 2\psi + \alpha_{1} + \frac{{\tau_{1}^{2} +\tau_{2}^{2} }}{2}} \right\} - \psi N$$
$$= - 2 + \frac{{\left( {b_{1} N_{0} + 1} \right)\lambda I}}{N} - p\left( {1 - \frac{S}{q}} \right) - \frac{{\alpha_{1} S}}{V} - \psi N$$
(6)

We define a compact subset \(E_{\mu } = \left\{ {\left( {S, V, I} \right) \in R_{ + }^{3} :S \ge \mu ;V \ge \mu^{2} ;I \ge \mu^{3} ;S + V + I \le \frac{1}{\mu }} \right\}\).

Here, μ is a small positive constant which satisfy the condition

$$- 2 + \left( {b_{1} N_{0} + 1} \right)\lambda - \frac{{\min \left( {\psi , \alpha_{1} ,p} \right)}}{\mu } \le - 1$$
(7)
$$- 2 + \left( {b_{1} N_{0} + 1} \right)\lambda \mu \le - 1$$
(8)

Using the following subsets of \(R_{ + }^{3} \backslash E_{\mu }\)

$$E_{1,\mu }^{c} = \left\{ {\left( {S, V,I} \right) \in R_{ + }^{3} :S < \mu } \right\};E_{2,\mu }^{c} = \left\{ {\left( {S, V,I} \right) \in R_{ + }^{3} :V < \mu^{2} ,S \ge \mu } \right\}$$
$$E_{3,\mu }^{c} = \{ \left( {S, V,I} \right) \in R_{ + }^{3} :I < \mu^{3} , V \ge \mu^{2} \} ;\;E_{4,\mu }^{c} = \{ \left( {S, V,I} \right) \in R_{ + }^{3} :S + V + I > \frac{1}{\mu }\}$$

Therefore, \(LQ \le - 1,\) where \(\left( {S, V,I} \right) \in E_{i,\mu } \left( {i = 1,2,3,4} \right)\).

  1. (i)

    When \(\left( {S, V,I} \right) \in E_{1,\mu }\) then we get from (6) & (7)

    $$LQ \le - 2 + \left( {b_{1} N_{0} + 1} \right)\lambda - \frac{{\alpha_{1} S}}{V} \le - 2 + \left( {b_{1} N_{0} + 1} \right)\lambda - p\left( {1 - \frac{S}{q}} \right) \le - 1$$
  2. (ii)

    When \(\left( {S, V,I} \right) \in E_{2,\mu }\) then we get from (6) & (7)

    $$LQ \le - 2 + \left( {b_{1} N_{0} + 1} \right)\lambda - \frac{{{\text{min}}\left( {\psi , \alpha_{1} ,p} \right)}}{\mu } \le - 1$$
  3. (iii)

    When \(\left( {S, V,I} \right) \in E_{3,\mu }\) then we get from (6) & (7)

    $$LQ \le - 2 + \left( {b_{1} N_{0} + 1} \right)\frac{\lambda I}{V} \le - 2 + \left( {b_{1} N_{0} + 1} \right)\lambda \mu \le - 1$$
  4. (iv)

    When \(\left( {S, V,I} \right) \in E_{4,\mu }\) then we get from (6) & (8)

    $$LQ \le - 2 + \left( {b_{1} N_{0} + 1} \right)\lambda - \psi N \le - 2 + \left( {b_{1} N_{0} + 1} \right)\lambda - \frac{{{\text{min}}\left( {\psi , \alpha_{1} ,p} \right)}}{\mu } \le - 1$$

So, \(LQ \le - 1\), for all \(\left( {S, V, I} \right) \in R_{ + }^{3} \backslash E_{\mu } = E_{1,\mu } \cup E_{2,\mu } \cup E_{3,\mu } \cup E_{4,\mu }\).

Therefore, the condition (i) of Lemma 1 holds. The corresponding diffusion matrix is given by.

\(L = \left( {\begin{array}{*{20}c} {\tau_{1}^{2} S^{2} } & 0 & 0 \\ 0 & {\tau_{2}^{2} V^{2} } & 0 \\ 0 & 0 & {\tau_{3}^{2} I^{2} } \\ \end{array} } \right)\).

Here, \(\tau_{1}^{2} S^{2} > 0,\ \tau_{2}^{2} V^{2} > 0,\ \tau_{3}^{2} I^{2} > 0\), so \(L\) is positive definite matrix. Then, the condition (ii) of Lemma 1 also holds. Therefore, the global positive solution \(\left( {S\left( t \right), V\left( t \right),I\left( t \right)} \right)\) of the system (2) satisfies a unique ergodic stationary distribution \(\pi \left( \cdot \right)\). This completes the proof of Theorem (2).

Density function analysis

In this portion, we have to obtain the probability density function of the stochastic model (2). For this purpose, we apply logarithmic transformation and equilibrium offset transformation.

Let us assume that \(g_{1} = \text{ln}S,\ g_{2} = \text{ln}V,\ g_{3} = \text{ln} I.\) Using It \(\hat{o}\) s formula in the system (2), we derive

$$dg_{1} = p\left( {1 - \frac{{e^{{g_{1} }} }}{q}} \right) - \frac{{\lambda e^{{g_{3} }} }}{{1 + k_{1} e^{{g_{3} }} }} - \left( {\psi + \alpha_{1} + \frac{{\tau_{1}^{2} }}{2}} \right) + \beta_{2} \frac{{e^{{g_{3} }} }}{{e^{{g_{1} }} }}$$
$$dg_{2} = \alpha_{1} \frac{{e^{{g_{1} }} }}{{e^{{g_{2} }} }} - \left( {\beta_{1} + \psi + \frac{{\tau_{2}^{2} }}{2}} \right)$$
$$dg_{3} = \frac{{\lambda e^{{g_{1} }} }}{{1 + k_{1} e^{{g_{3} }} }} - \left( {\psi + \alpha_{2} + \beta_{2} + \frac{{\tau_{3}^{2} }}{2}} \right)$$
(9)

For considering random effect, we have a critical value

$$R_{0}^{c} = \frac{{\lambda q\left\{ {\alpha_{1} \beta_{1} + \left( {\beta_{1} + \psi + \frac{{\tau_{2}^{2} }}{2}} \right)\left( {p - \psi - \alpha_{1} - \frac{{\tau_{3}^{2} }}{2}} \right)} \right\}}}{{p\left( {\beta_{1} + \psi + \frac{{\tau_{2}^{2} }}{2}} \right)\left( {\psi + \alpha_{2} + \beta_{2} + \frac{{\tau_{3}^{2} }}{2}} \right)}}$$

If \(R_{0}^{c} > 1,\) the equation.\(f\left( x \right) = E_{1} V^{2} + E_{2} V + E_{3}\) has unique positive root \(V_{ + }^{*}\). where \(E_{1} = \frac{{\left( {\psi + \alpha_{2} + \beta_{2} } \right)p}}{{\alpha_{2} \left( {\psi + \alpha_{2} } \right)q}} > 0;E_{2} = \frac{\lambda }{{\alpha_{1} \left( {\beta_{1} + \psi } \right)}} - \frac{{p\left( {\psi + \alpha_{2} + \beta_{2} } \right)}}{{\alpha_{1} \left( {\beta_{1} + \psi } \right)\left( {\psi + \alpha_{2} } \right)}} > 0;E_{3} = \frac{{\left( {\psi + \alpha_{2} + \beta_{2} } \right)}}{{\left( {\beta_{1} + \psi } \right)^{2} }} > 0\)

Clearly, \(\left( {S_{ + }^{*} , V_{ + }^{*} , I_{ + }^{*} } \right) = \left( {e^{{g_{1}^{*} }} , e^{{g_{2}^{*} }} , e^{{g_{3}^{*} }} } \right)\) are same with the endemic equilibrium point \(E^{@} = \left( {S^{@} , V^{@} ,I^{@} } \right)\) of the system (2).

Now, we choose \(\left( {h_{1} ,h_{2} ,h_{3} } \right) = \left( {g_{1} - g_{1}^{*} , g_{2} - g_{2}^{*} , g_{3} - g_{3}^{*} } \right)\) such that \(g_{1}^{*} = {\text{ln}}S_{ + }^{*} , g_{2}^{*} = {\text{ln}}V_{ + }^{*} ,g_{3}^{*} = {\text{ln}}I_{ + }^{*} .\) The system (9) can be represented as linear form,

$${\text{d}}g_{1} = \left( { - b_{11} h_{1} - b_{12} h_{2} - b_{13} h_{3} } \right){\text{d}}t +\tau_{1} {\text{d}}G_{1} \left( t \right)$$
$${\text{d}}g_{2} = \left( {b_{21} h_{1} - b_{21} h_{2} } \right){\text{d}}t +\tau_{2} {\text{d}}G_{2} \left( t \right)$$
$${\text{d}}g_{3} = \left( {b_{31} h_{1} - b_{32} h_{2} } \right){\text{d}}t +\tau_{3} {\text{d}}G_{3} \left( t \right)$$
(10)

where \(b_{11} = \frac{p}{q}e^{{g_{1}^{*} }}\), \(b_{12} = \frac{{\lambda e^{{g_{3}^{*} }} }}{{\left( {1 + k_{1} e^{{g_{3}^{*} }} } \right)^{2} }}, b_{13} = \psi + \alpha_{1} , b_{21} = \beta_{1} + \psi\), \(b_{31} = \frac{{\lambda e^{{g_{1}^{*} }} }}{{1 + k_{1} e^{{g_{3}^{*} }} }}\), \(b_{32} = \psi + \alpha_{2} + \beta_{2} .\) It is obvious that \(b_{11} > 0\), \(b_{12} > 0,b_{13} > 0,b_{21} > 0,b_{31} > 0\), \(b_{32} > 0\).

Theorem 3

When \(R_{0}^{C} > 1,\) then there exists a local normal density function \({\mathcal{T}}\left( {h_{1} ,h_{2} ,h_{3} } \right)\) at the quasi-stationary state \(\left( {S_{ + }^{*} , V_{ + }^{*} , I_{ + }^{*} } \right)\) with the initial value \(\left( {h_{1} \left( 0 \right), h_{2} \left( 0 \right), h_{3} \left( 0 \right)} \right)\) which satisfy the following condition.

$${\mathcal{T}}\left( {h_{1} ,h_{2} ,h_{3} } \right) = \left( {2\pi } \right)^{{ - \frac{3}{2}}} \left| \rho \right|^{{ - \frac{1}{2}}} e^{{ - \frac{1}{2}\left( {h_{1} ,h_{2} ,h_{3} } \right)\rho^{ - 1} \left( {h_{1} ,h_{2} ,h_{3} } \right)^{T} }}$$

where \(\rho\) can be written as following way, and \(\rho\) is a positive definite matrix. If \(\Lambda \ne V_{ + }^{*}\), then

$$\rho = r_{1}^{2} \left( {N_{1} L_{1} } \right)^{ - 1} \rho_{0} \left[ {\left( {N_{1} L_{1} } \right)^{ - 1} } \right]^{T} + r_{2}^{2} \left[ {\left( {N_{2} L_{2} } \right)^{ - 1} \rho_{0} \left( {N_{2} L_{2} } \right)^{ - 1} } \right] +$$
$$r_{3}^{2} \left( {N_{3} Q_{3} L_{3} } \right)^{ - 1} \rho_{0} \left[ {\left( {N_{3} Q_{3} L_{3} } \right)^{ - 1} } \right]^{T}$$

If \(\Lambda = V_{ + }^{*}\), then

$$\rho = r_{1}^{2} \left( {N_{1} L_{1} } \right)^{ - 1} \rho_{0} \left[ {\left( {N_{1} L_{1} } \right)^{ - 1} } \right]^{T} + r_{2}^{2} \left[ {\left( {N_{2} L_{2} } \right)^{ - 1} \rho_{0} \left( {N_{2} L_{2} } \right)^{ - 1} } \right] + r_{3u}^{2} \left( {N_{3u} Q_{3} L_{3} } \right)^{ - 1} \xi_{0} \left[ {N_{3u} Q_{3} L_{3} } \right)^{ - 1} ]^{T}$$

where \(\Lambda = \frac{{\left( {\psi + \frac{{\tau_{2}^{2} }}{2}} \right)k_{1} - \left( {\psi + \frac{{\tau_{1}^{2} }}{2}} \right)k_{1} }}{{\psi + \frac{{\tau_{1}^{2} }}{2} + k_{1} - \psi - \frac{{\tau_{2}^{2} }}{2}}},u = \frac{{b_{21} \left( {b_{11} - b_{21} } \right)}}{{b_{13} }},\ \)\(r_{1} = - b_{21} b_{31}\tau_1\), \(r_{2} = - b_{13} b_{32}\tau_2 ,r_{3} = - b_{21} \left( {b_{11} - b_{22} } \right)\tau_{3}\), \(r_{3u} = - b_{13}\tau_3\),

$$\rho_{0} = \left( {\begin{array}{*{20}c} {\frac{{B_{2} }}{{2\left( {B_{1} B_{2} - B_{3} } \right)}}} & 0 & { - \frac{1}{{2\left( {B_{1} B_{2} - B_{3} } \right)}}} \\ 0 & {\frac{1}{{2\left( {B_{1} B_{2} - B_{3} } \right)}}} & 0 \\ { - \frac{1}{{2\left( {B_{1} B_{2} - B_{3} } \right)}}} & 0 & {\frac{{B_{1} }}{{2B_{3} \left( {B_{1} B_{2} - B_{3} } \right)}}} \\ \end{array} } \right),\;\xi_{0} = \left( {\begin{array}{*{20}c} {\frac{1}{{2d_{1} }}} & 0 & 0 \\ 0 & {\frac{1}{{2d_{1} d_{2} }}} & 0 \\ 0 & 0 & 0 \\ \end{array} } \right)$$

\(B_{1} = b_{11} + b_{21} + b_{32}\), \(B_{2} = b_{21} \left( {b_{11} - b_{12} + b_{32} } \right) + \left[ {b_{11} b_{13} \left( {b_{32} + b_{32} } \right)} \right]\), \(B_{3} = b_{21} b_{32} (b_{11} - b_{12} - b_{32} )\),\(d_{1} = b_{11} + b_{21} ,d_{2} = b_{21} \left( {b_{11} - b_{12} - b_{32} } \right)\)

$$N_{3} = \left( {\begin{array}{*{20}c} { - b_{32} u} & { - \left( {b_{11} + b_{21} } \right)u} & {b_{21}^{2} } \\ 0 & u & { - b_{21} } \\ 0 & 0 & 1 \\ \end{array} } \right),\;N_{3u} = \left( {\begin{array}{*{20}c} { - b_{13} } & { - b_{11} } & 0 \\ 0 & u & 0 \\ 0 & 0 & 1 \\ \end{array} } \right),\,Q_{3} = \left( {\begin{array}{*{20}c} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & { - \frac{{b_{21} }}{{b_{31} }}} & 1 \\ \end{array} } \right)$$
$$L_{1} = \left( {\begin{array}{*{20}c} 0 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{array} } \right),\;L_{3} = \left( {\begin{array}{*{20}c} 0 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{array} } \right),\;N_{1} = \left( {\begin{array}{*{20}c} { - b_{21} b_{31} } & {b_{21} \left( {b_{22} + b_{31} } \right)} & {b_{21}^{2} + b_{21} b_{32} } \\ 0 & { - b_{23} } & { - b_{21} } \\ 0 & 0 & 1 \\ \end{array} } \right)$$
$$N_{2} = \left( {\begin{array}{*{20}c} { - b_{13} b_{32} } & {b_{13} \left( {b_{11} + b_{32} } \right)} & {b_{11}^{2} + b_{13} b_{31} } \\ 0 & { - b_{13} } & { - b_{11} } \\ 0 & 0 & 1 \\ \end{array} } \right)$$

Lemma 2

(Zhou et al. 2020, 2021) Let us consider the equation \(H_{0}^{2} + B_{0} \rho_{0} + \rho_{0} B_{0}^{T} = 0\), with \(H_{0} =\) diag \(\left( {1, 0, 0} \right)\), \(B_{0} = \left( {\begin{array}{*{20}c} { - B_{1} } & { - B_{2} } & { - B_{3} } \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{array} } \right).\) If \(B_{1} > 0,B_{2} > 0\) and \(B_{1} B_{2} - B_{3} > 0\) then matrix \(\rho_{0}\) have all positive eigen values. Here, \(B_{0}\) is the standard \(R_{1}\) matrix and \(B_{1} , B_{2} , B_{3}\) are the coefficients of the characteristic polynomial \(x^{3} + B_{1} x^{2} + B_{2} x + B_{3}\) of \(B_{0} .\)

Lemma 3

(Zhou et al. 2020, 2021) Let us consider the equation \(H_{0}^{2} + C_{0} \xi_{0} + \xi_{0} C_{0}^{T} = 0\). Here, \(H_{0} =\) diag \(\left( {1, 0, 0} \right)\),

$$C_{0} = \left( {\begin{array}{*{20}c} { - d_{1} } & { - d_{2} } & { - d_{3} } \\ 1 & 0 & 0 \\ 0 & 0 & {d_{33} } \\ \end{array} } \right)$$

If \(d_{1} > 0\) and \(d_{2} > 0\) then \(\xi_{0}\) is semi-positive. Here, \(d_{1} - d_{33} ,d_{2} - d_{1} d_{33}\) and \(- d_{2} d_{33}\) are the coefficient of the polynomial \(x^{3} + \left( {d_{1} - d_{33} } \right)x^{2} + (d_{2} - d_{1} d_{33} )x + - d_{2} d_{33}\) of \(C_{0} .\) Here, \(C_{0}\) represents standard \(R_{2}\) matrix.

Proof of Theorem

From the system of equation of (10), we construct.

\(X = \left( {h_{1} ,h_{2} ,h_{3} } \right)^{T}\), \(Y =\) diag \(\left( {\tau_{1} , \tau_{2} , \tau_{3} } \right).\) The system becomes.\(dX = BXdt + YdB\left( t \right)\) where \(B = \left( {\begin{array}{*{20}c} { - b_{11} } & 0 & { - b_{13} } \\ 0 & { - b_{22} } & { - b_{23} } \\ {b_{31} } & {b_{32} } & { - b_{33} } \\ \end{array} } \right).\)

The Fokker–Plank equation is given by

$$- \mathop \sum \limits_{j = 1}^{3} \frac{{\tau_{j}^{2} }}{2}\frac{{\partial^{2} }}{{\partial h_{j}^{2} }}T + \frac{\partial }{{\partial h_{1} }}\left[ {\left( { - b_{11} h_{1} - b_{13} h_{3} } \right)T} \right] + \frac{\partial }{{\partial h_{2} }}\left[ {\left( { - b_{21} h_{2} - b_{23} h_{3} } \right)T} \right] + \frac{\partial }{{\partial h_{3} }}\left[ {(b_{31} h_{1} + b_{32} h_{2} - b_{32} h_{3} } \right)T] = 0.$$
(11)

The above three dimension Fokker–Plank equation approximates the density function \(T\left( X \right) = T\left( {h_{1} ,h_{2} ,h_{3} } \right)\) with help of (Roozen 1989) research work.

Eq. (11) can be approximated with the help of Gaussian distribution

$$T\left( X \right) = c {\text{exp}}\left\{ { - \frac{1}{2}\left( {X - X^{*} } \right)\theta \left( {X - X^{*} } \right)^{T} } \right\}$$
(12)

Here, \(X = \left( {0, 0, 0} \right)\) and the real symmetric matrix \(\theta\) satisfies the equation.

\(\theta Y^{2} \theta + B^{T} \theta + \theta B = 0.\) If all the eigenvalues of \(\theta\) are positive and \(\theta^{ - 1} = \rho\), we obtain

$$Y^{2} + B\rho + \rho B^{T} = 0$$
(13)

Equation (13) is equivalent to the equation

$$Y_{j}^{2} + B\rho_{j} + \rho_{j} B^{T} = 0,\left( {j = 1,2,3} \right)$$

where \(Y_{1} =\) diag \(\left( {\tau_{1} ,0,0} \right)\), \(Y_{2} =\) diag \(\left( {0,\tau_{2} ,0} \right)\), \(Y_{3} =\) diag \(\left( {0,0,\tau_{3} } \right)\)

$$\rho_{j} = \rho_{1} + \rho_{2} + \rho_{3}$$

\(Y_{j}^{2} = Y_{1}^{2} + Y_{2}^{2} + Y_{3}^{2}.\)

Now, we choose a polynomial of \(B\) such that

$$\psi_{B} \left( x \right) = x^{3} + b_{1} x^{2} + b_{2} x + b_{3}$$

where \(b_{1} = b_{11} + b_{21} + b_{32} > 0\), \(b_{2} = b_{11} \left( {b_{21} + b_{32} } \right) + b_{21} b_{32} + b_{12} b_{32} + b_{13} b_{31} > 0\), \(b_{3} = b_{21} b_{32} \left( {b_{11} - b_{12} - b_{13} } \right) > 0.\)

Applying the theory of matrix similar transformation \(\psi_{B} \left( x \right)\) is similarly invariant.

For solving Eq. (13), we have followed some steps.

For step 1

We choose the algebraic equation

$$Y_{1}^{2} + B\rho_{1} + \rho_{1} B^{T} = 0$$
(14)

where \(Y_{1} =\) diag \(\left( {\tau_{1} , 0, 0} \right)\)

Let, \(B_{1} = J_{1} BJ_{1}^{T},\) \(J_{1} = \left( {\begin{array}{*{20}c} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \\ \end{array} } \right)\) then, \(B_{1} = \left( {\begin{array}{*{20}c} { - b_{11} } & { - b_{13} } & 0 \\ {b_{31} } & { - b_{32} } & 0 \\ 0 & { - b_{21} } & { - b_{21} } \\ \end{array} } \right).\)

\(C_{1} = N_{1} B_{1} N_{1}^{ - 1},\) where \(N_{1} = \left( {\begin{array}{*{20}c} { - b_{21} b_{31} } & {b_{21} \left( {b_{21} + b_{33} } \right)} & {b_{21}^{2} + b_{21} b_{32} } \\ 0 & { - b_{21} } & { - b_{32} } \\ 0 & 0 & 1 \\ \end{array} } \right).\)

With the help of uniqueness of standard \(R_{1}\) matrix of \(B\), we derive.

\(B = \left( {\begin{array}{*{20}c} { - b_{1} } & { - b_{2} } & { - b_{3} } \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{array} } \right)\).

Then, we obtain \(b_{1} b_{2} - b_{3} = \left( {b_{21} + b_{32} } \right)\left[ {b_{11} \left( {b_{21} + b_{32} } \right) + b_{21} b_{32} } \right] + \left( {b_{11} + b_{32} } \right)b_{13} b_{31} ] > 0\).

From Eq. (14), it can be seen that

$$(N_{1} J_{1} )Y_{1}^{2} (N_{1} J_{1} )^{T} + (N_{1} J_{1} )B(N_{1} J_{1} )^{ - 1} (N_{1} J_{1} )\rho_{1} (N_{1} J_{1} )^{T} + (N_{1} J_{1} )\rho_{1} (N_{1} J_{1} )^{T} \left[ {(N_{1} J_{1} } \right)B(N_{1} J_{1} )^{ - 1} ]^{T} = 0.$$

This gives \(Y_{0}^{2} + C_{1} \rho_{0} + \rho_{0} C_{1}^{T} = 0\). where \(\rho_{0} = \frac{{(N_{1} J_{1} )\rho_{1} (N_{1} J_{1} )^{T} }}{{r_{1}^{2} }}\), \(r_{1} = - b_{32} b_{31}\tau_1\) and \(\rho_{0} = \left( {\begin{array}{*{20}c} {\frac{{b_{2} }}{{2\left( {b_{1} b_{2} - b_{3} } \right)}}} & 0 & { - \frac{1}{{2\left( {b_{1} b_{2} - b_{3} } \right)}}} \\ 0 & {\frac{1}{{2\left( {b_{1} b_{2} - b_{3} } \right)}}} & 0 \\ { - \frac{1}{{2\left( {b_{1} b_{2} - b_{3} } \right)}}} & 0 & {\frac{{b_{1} }}{{2b_{3} \left( {b_{1} b_{2} - b_{3} } \right)}}} \\ \end{array} } \right)\).

Hence, \(\rho_{1} = r_{1}^{2} (N_{1} J_{1} )^{ - 1} \rho_{0} \left[ {(N_{1} J_{1} )^{ - 1} } \right]^{T}\) has positive eigen value, i.e., positive definite.

For step 2

$${\text{Let us consider the equation be}}\;Y_{2}^{2} + B\rho_{2} + \rho_{2} B^{T} = 0$$
(15)

where \(Y_{2} =\) diag \(\left( {0, \tau_{2} , 0} \right)\), \(B_{2} = \left( {J_{2} BJ_{2}^{ - 1} } \right)\) and \(J_{2} = \left( {\begin{array}{*{20}c} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \\ \end{array} } \right)\).

Thus, the matrix \(B_{2}\) can be represented into the following form.

\(B_{2} = \left( {\begin{array}{*{20}c} { - b_{21} } & { - b_{32} } & 0 \\ {b_{32} } & { - b_{32} } & {b_{31} } \\ 0 & { - b_{13} } & { - b_{11} } \\ \end{array} } \right)\).

Again, \(C_{2} = N_{2} B_{2} N_{2}^{ - 1}\) and \(N_{2} = \left( {\begin{array}{*{20}c} { - b_{13} b_{32} } & {b_{13} \left( {b_{11} + b_{32} } \right)} & {b_{11}^{2} + b_{13} b_{31} } \\ 0 & { - b_{13} } & { - b_{11} } \\ 0 & 0 & 1 \\ \end{array} } \right)\).

Here, \(C_{2}\) is standard \(R_{1}\) matrix and \(C_{2} = C_{1}\).

Equation (15) is transformed into the following form

$$\left( {N_{2} J_{2} } \right)Y_{2}^{2} \left( {N_{2} J_{2} } \right)^{T} + \left( {N_{2} J_{2} } \right)B(N_{2} J_{2} )^{ - 1} (N_{2} J_{2} )\rho_{2} (N_{2} J_{2} )^{T} + (N_{2} J_{2} )\rho_{2} (N_{2} J_{2} )^{T} \left( {\left( {N_{2} J_{2} } \right)B(N_{2} J_{2} )^{ - 1} } \right)^{T} = 0$$

This gives us \(Y_{0}^{2} + C_{2} \rho_{0} + \rho_{0} C_{2}^{T} = 0.\)

All the eigen values of the matrix \(\rho_{0}\) are positive and \(\rho_{2} = r_{2}^{2} (N_{2} J_{2} )^{ - 1} \rho_{0} [(N_{2} J_{2} )^{ - 1} ]^{T}\) is positive definite.

For step 3

Consider the equation

$$Y_{3}^{2} + B\rho_{3} + \rho_{3} B^{T} = 0$$
(16)

where \(Y_{3} =\) diag \(\left( {0,0,\tau_{3} } \right)\), \(B_{3} = J_{3} BJ_{3}^{T}\) and \(J_{3} = \left( {\begin{array}{*{20}c} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{array} } \right)\).

Therefore, we obtain \(Y_{3} = \left( {\begin{array}{*{20}c} { - b_{32} } & {b_{31} } & {b_{32} } \\ { - b_{31} } & { - b_{11} } & 0 \\ { - b_{21} } & 0 & { - b_{21} } \\ \end{array} } \right)\).

Using the transformation \(C_{3} = N_{3} B_{3} N_{3}^{ - 1}\), where \(N_{3} = \left( {\begin{array}{*{20}c} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & { - \frac{{b_{21} }}{{b_{13} }}} & 1 \\ \end{array} } \right)\), We derive \(C_{3} = \left( {\begin{array}{*{20}c} { - b_{32} } & {\frac{{b_{13} + b_{11} + b_{21} b_{32} }}{{b_{13} }}} & {b_{32} } \\ { - b_{13} } & { - b_{11} } & 0 \\ 0 & {\frac{{b_{21} \left( {b_{11} - b_{21} } \right)}}{{b_{13} }}} & { - b_{21} } \\ \end{array} } \right)\).

Case-1

When \(U = \frac{{b_{21} \left( {b_{11} - b_{21} } \right)}}{{b_{13} }} \ne 0\) and \(\Lambda \ne 0\), from step (1) and step (2),we obtain \(D_{3} = N_{3} C_{3} N_{3}^{ - 1} , N_{3} = \left( {\begin{array}{*{20}c} { - b_{13} U} & { - \left( {b_{11} + b_{21} } \right)U} & {b_{21}^{2} } \\ 0 & U & { - b_{21} } \\ 0 & 0 & 1 \\ \end{array} } \right)\).

Transforming Eq. (15), it can be seen that.

\(Y_{0}^{2} + D_{3} \rho_{0} + \rho_{0} D_{3}^{T} = 0\),where \(\rho_{0} = \frac{{(N_{3} D_{3} J_{3} )\rho_{3} (N_{3} D_{3} J_{3} )^{T} }}{{r_{3}^{2} }}\) and \(r_{3} = - b_{21} \left( {b_{11} - b_{21} } \right)\tau_{3} .\)

Therefore, the matrix \(\rho_{3}\) has positive eigenvalues.

Case-2: If \(U = \frac{{b_{21} \left( {b_{11} - b_{21} } \right)}}{{b_{13} }} = 0\), \(\Lambda = 0\) and \(D_{3U} = N_{3U} C_{3} N_{3U}^{ - 1}\), then the transformation matrix is given by \(N_{3U} = \left( {\begin{array}{*{20}c} { - b_{13} } & { - b_{11} } & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{array} } \right)\).

The matrix \(D_{3U}\) is represented by \(\left( {\begin{array}{*{20}c} { - C_{1} } & { - C_{2} } & { - C_{3} } \\ 1 & 0 & 0 \\ 0 & 0 & { - b_{11} } \\ \end{array} } \right).\)

With the help of similarity invariant of the characteristic polynomial of \(B\) is

$$\Phi_{B} \left( x \right) = x^{3} + b_{1} x^{2} + b_{2} x + b_{3} = \Phi_{{D_{3U} }} \left( x \right) = x^{3} + \left( {C_{1} + b_{11} } \right)x^{2} + \left( {C_{2} + C_{1} b_{11} } \right)x + C_{2} b_{11} .$$

Above expression gives us,

$$C_{1} = b_{1} - b_{11} = b_{21} + b_{32} > 0$$
$$C_{2} = b_{2} - C_{1} b_{11} = b_{21} b_{32} + b_{21} b_{12} + b_{13} b_{31} > 0.$$

Transforming Eq. (16), we get

$$Y_{0}^{2} + D_{3U} W_{0} + W_{0} D_{3U}^{T} = 0$$

where \(W_{0} = \frac{{(N_{3U} D_{3} J_{3} )\rho_{3} (N_{3U} D_{3} J_{3} )^{T} }}{{r_{3u}^{2} }}\), \(N_{3U} = - b_{13} \tau_{3}\).

Hence, \(W_{0}\) is semi-positive definite and \(W_{0} = \left( {\begin{array}{*{20}c} {\frac{1}{{2C_{1} }}} & 0 & 0 \\ 0 & {\frac{1}{{2C_{1} C_{2} }}} & 0 \\ 0 & 0 & 0 \\ \end{array} } \right)\).

Therefore, \(\rho_{3} = r_{3U}^{2} (N_{3U} D_{3} J_{3} )^{ - 1} W_{0} \left( {(N_{3U} D_{3} J_{3} )^{ - 1} } \right)^{T}\).

Now, \(\rho = \rho_{1} + \rho_{2} + \rho_{3}\) becomes real symmetric matrix, and hence, Eq. (13) represents the positive definite matrix.

The proof is completed.

Remark

Theorem 3 shows that there is local and exact probability density function around the quasi-stationary state \(\left( {S_{ + }^{*} , V_{ + }^{*} , I_{ + }^{*} } \right)\) if \(R_{0}^{c} > 1\), we easily obtain \(R_{0}^{c} < R_{S}^{*} .\) In addition \(R_{S}^{*} = R_{0}^{c} = R^{*}\) while \(\tau_{i} = 0 \left( {i = 1, 2, 3} \right)\). This means that the disease persistence is critically affected by the random fluctuation of the susceptible and infected individuals.

Extinction of the stochastic system

In this section, we discuss the suitable condition for the extinction of the system (2).

Theorem 4

Let, the solutions \(\left( {S\left( t \right), V\left( t \right), I\left( t \right)} \right)\) of the system (2) with initial values \(\left( {S\left( 0 \right), V\left( 0 \right), I\left( 0 \right)} \right) \in R_{ + }^{3} ,\) the epidemic will be eradicated for the long run if \(R_{0}^{*} = \frac{\lambda }{{\psi + \alpha_{2} + \beta_{2} + \frac{{\tau_{3}^{2} }}{2}}} < 1\) and \(\mathop {\lim }\limits_{t \to \infty } \sup \frac{{{\text{ln}}I\left( t \right)}}{t} \le \left( {\psi + \alpha_{2} + \beta_{2} + \frac{{\tau_{3}^{2} }}{2}} \right)\left( {R_{0}^{*} - 1} \right) < 0\) a.s.

Proof

Applying It \({\hat{\text{o}}}\)’s formula in \(\text{In}I\left( t \right)\), we obtain

$${\text{d ln}}I\left( t \right) = \left[ {\frac{\lambda S\left( t \right)}{{1 + k_{1} I\left( t \right)}} - \left( {\psi + \alpha_{2} + \beta_{2} + \frac{{\tau_{3}^{2} }}{2}} \right)} \right]{\text{d}}t + \tau_{3} {\text{d}}G_{3} \left( t \right)$$
(17)

By integrating both sides and limiting runs from \(0\) to \(t_{1}\), we can get

$$\frac{{{\text{ln}}I\left( t \right)}}{t} \le \frac{{{\text{ln}}I\left( 0 \right)}}{t} + \frac{1}{t}\mathop \smallint \limits_{0}^{{t_{1} }} \left[ {\frac{\lambda S\left( n \right)}{{1 + k_{1} I\left( n \right)}} - \left( {\psi + \alpha_{2} + \beta_{2} + \frac{{\tau_{3}^{2} }}{2}} \right)} \right]{\text{d}}n + \frac{{\mathop \smallint \nolimits_{0}^{{t_{1} }} \tau_{3} dG_{3} \left( t \right)}}{t}$$
$$\le \frac{{{\text{ln}}I\left( 0 \right)}}{t} + \frac{1}{t}\left[ {\lambda - \left( {\psi + \alpha_{2} + \beta_{2} + \frac{{\tau_{3}^{2} }}{2}} \right)} \right]{\text{d}}n + \frac{{\mathop \smallint \nolimits_{0}^{{t_{1} }} \tau_{3} {\text{d}}G_{3} \left( t \right)}}{t}$$
$$= \frac{{{\text{ln}}I\left( 0 \right)}}{t} + \left( {\psi + \alpha_{2} + \beta_{2} + \frac{{\tau_{3}^{2} }}{2}} \right)\left( {R_{0}^{*} - 1} \right) + \frac{{\mathop \smallint \nolimits_{0}^{{t_{1} }} \tau_{3} {\text{d}}G_{3} \left( t \right)}}{t}$$
(18)

With the help of Strong law of large number (Lipster 1980), we derive

$$\mathop {\lim }\limits_{t \to \infty } \frac{{\mathop \smallint \nolimits_{0}^{{t_{1} }} \tau_{3} dG_{3} \left( t \right)}}{t} = 0\;a.s.$$
(19)

Using (19) and applying the superior limit of \(t \to \infty\), we get \(R_{0}^{*} \le R_{S}^{*}\).

Hence, the proof is completed.

Numerical results

In this section, applying the higher-order method improved by Milstein (Higham 2001) we have the following discretization equation of the system (2)

$$S^{\omega + 1} = S^{\omega } + \left[ {pS^{\omega } \left( {1 - \frac{{S^{\omega } }}{q}} \right) - \frac{{\lambda S^{\omega } I^{\omega } }}{{1 + k_{1} I^{\omega } }} - \psi S^{\omega } - \alpha_{1} S^{\omega } + \beta_{2} I^{\omega } + \beta_{1} V^{\omega } } \right]\Delta t + \frac{{\tau_{1}^{2} }}{2}S^{\omega } \left( {\eta_{1\omega }^{2} - 1} \right)\Delta t +\tau_{1} S^{\omega } \sqrt {\Delta t} \eta_{1\omega }$$
$$V^{\omega + 1} =V^{\omega }+ \left[ {\alpha_{1} S^{\omega } - \left( {\beta_{1} + \psi } \right)V^{\omega } } \right]\Delta t + \frac{{\tau_{2}^{2} }}{2}V^{\omega } \left( {\eta_{2\omega }^{2} - 1} \right)\Delta t +\tau_{2} V^{\omega } \sqrt {\Delta t} \eta_{2\omega }$$
$$I^{\omega + 1} = I^{\omega }+\left[ {\frac{{\lambda S^{\omega } I^{\omega } }}{{1 + k_{1} I^{\omega } }} - \left( {\psi + \alpha_{2} + \beta_{2} } \right)I^{\omega } } \right]\Delta t + \frac{{\tau_{3}^{2} }}{2}I^{\omega } \left( {\eta_{3\omega }^{2} - 1} \right)\Delta t +\tau_{3} I^{\omega } \sqrt {\Delta t} \eta_{3\omega }$$
(20)

where \(\Delta t > 0\) represents the time increment, and the independent Gaussian random variables are indicated by the variables \(\eta_{1\omega } , \eta_{2\omega }\) and \(\eta_{3\omega }\). Here, all the random variables follow the distribution \(N\left( {0, 1} \right)\) for \(\omega = 1,2, \ldots ,n.\)

Result 1

We choose the environmental noise intensities \(\left( {\tau_{1} ,\tau_{2} ,\tau_{3} } \right) = \left( {0.08, 0.02,0.01} \right)\) and the other parameters are \(\left( {p, q, \lambda , k_{1} , \psi , \alpha_{1} ,\alpha_{2} , \beta_{1} ,\beta_{2} } \right) = \left( {2, 50, 0.9, 0.65, 0.02, 0.9, 0.03, 0.2, 0.05} \right)\). We can obtain.

\(R^{*} = \frac{{\lambda q\left\{ {\left( {p - \psi - \alpha_{1} } \right)\left( {\beta_{1} + \psi } \right) + \beta_{1} \alpha_{1} } \right\}}}{{p(\beta_{1} + \psi )\left( {\psi + \alpha_{2} + \beta_{2} } \right)}} > 1\) and \(R_{S}^{*} = \frac{{\lambda q\left\{ {\alpha_{1} \beta_{1} + \left( {\beta_{1} + \psi } \right)\left( {p - \psi - \alpha_{1} } \right)} \right\}}}{{\left( {p + \frac{{\tau_{1}^{2} }}{2}} \right)\left( {\beta_{1} + \psi + \frac{{\tau_{2}^{2} }}{2}} \right)\left( {\psi + \alpha_{2} + \beta_{2} + \frac{{\tau_{3}^{2} }}{2}} \right)}} > 1.\)

If \(R^{*} > 1\) and \(R_{S}^{*} > 1\) satisfies then Theorem 2 represents unique ergodic stationary distribution \(\pi \left( . \right).\) From Theorem 3, we have a unique log-normal density function around the quasi-endemic equilibrium \(E^{@} .\)

Result 2

We choose the environmental noise intensities \(\left( {\tau_{1} ,\tau_{2} ,\tau_{3} } \right) = \left( {0.08, 0.02,0.01} \right)\) and the other parameters are.

\(\left( {p, q, \lambda , k_{1} , \psi , \alpha_{1} ,\alpha_{2} , \beta_{1} ,\beta_{2} } \right) = \left( {2, 50, 0.00009, 0.65, 0.02, 0.9, 0.03, 0.2, 0.05} \right)\). We derive.

\(R_{0}^{*} = \frac{\lambda }{{\psi + \alpha_{2} + \beta_{2} + \frac{{\tau_{3}^{2} }}{2}}} = 0.00899 < 1\) and \(\mathop {\lim }\limits_{t \to \infty } \sup \frac{\text{ln}I\left( t \right)}{t} \le \left( {\psi + \alpha_{2} + \beta_{2} + \frac{{\tau_{3}^{2} }}{2}} \right)\left( {R_{0}^{*} - 1} \right) = - 0.0995 < 0\). It follows that Theorem 4 is verified and the epidemic will be eradicated for long run (Fig. 1). Figure 2 depicts the extinction of the system (2) for the long-term behavior.

Fig. 1
figure 1

Population trajectories in the system of (i) Deterministic system (ii) Stochastic system under the noise intensities \(\tau_{1} = 0.08, \tau_{2} = 0.02, \tau_{3} = 0.01.\)

Full size image
Fig. 2
figure 2

Extinction of the compartment of \(I\left( t \right)\) under the noise \(\tau_{1} = 0.08, \tau_{2} = 0.02, \tau_{3} = 0.01\) and the other parameters \(\left( {p, q, \lambda , k_{1} , \psi , \alpha_{1} ,\alpha_{2} , \beta_{1} ,\beta_{2} } \right) = \left( {2, 50, 0.00009, 0.65, 0.02, 0.9, 0.03, 0.2, 0.05} \right)\)

Full size image

Figure 3 depicts the impact of random noises \(\tau_{i} \left( {i = 1,2,3} \right)\) on all compartment of \(S\left( t \right),V\left( t \right), I\left( t \right)\) and Fig. 4 and Fig. 5 represent the effect of the noise intensities of \(\tau_{1} ,\tau_{2} ,\tau_{3}\) on the compartment of \(V\left( t \right), I\left( t \right)\). The stochastic perturbations are.

  1. (i)

    \(\tau_{1} = 0.08, \tau_{2} = 0.02, \tau_{3} = 0.01\) ii)\(\tau_{1} = 0.008, \tau_{2} = 0.02, \tau_{3} = 0.01\)

  2. (iii)

    \({ }\tau_{1} = 0.08, \tau_{2} = 0.002, \tau_{3} = 0.01\) iv) \(\tau_{1} = 0.08, \tau_{2} = 0.02, \tau_{3} = 0.001\), respectively.

Fig. 3
figure 3

Profile of the compartment \(S\left( t \right), V\left( t \right), I\left( t \right)\) in deterministic system and stochastic system with intensities \(\tau_{1} = 0.08, \tau_{2} = 0.02\), \(\tau_{3} = 0.01\) and other parameters are \(\left( {p, q,\lambda ,k_{1} ,\psi , \alpha_{1} , \alpha_{2} , \beta_{1} ,\beta_{2} } \right) = \left( {2, 50, 0.9, 0.65, 0.02, 0.9, 0.03, 0.2, 0.05} \right).\)

Full size image
Fig. 4
figure 4

Effect of the noise intensities \(\tau_{1} ,\tau_{2} ,\tau_{3}\) on the compartment of \(V\left( t \right)\)

Full size image
Fig. 5
figure 5

Effect of the noise intensities \(\tau_{1} ,\tau_{2} ,\tau_{3}\) on the compartment of \(I\left( t \right)\)

Full size image

For all the above cases, Fig. 3 and Fig. 4 represent the existence of a stationary distribution. The stationary distribution has an ergodicity property. But we mainly concentrate on the compartment of \(V\left( t \right)\) and \(I\left( t \right).\) When the perturbation intensities of vaccinated population or the infected population (\(\tau_{2}\) or \(\tau_{3}\)) increase, then the disease infection will be controlled.

Sensitivities of the parameters

In this section, we have discussed the sensitivities of the biological parameters which influence greatly for all compartment. We consider that the epidemiological parameters are \(\left( {p, q,\lambda ,k_{1} ,\psi , \alpha_{1} , \alpha_{2} , \beta_{1} ,\beta_{2} } \right) = \left( {2, 50, 0.9, 0.65, 0.02, 0.9, 0.03, 0.2, 0.05} \right)\) and the random noises are \(\tau_{1} = 0.08, \tau_{2} = 0.02\), and \(\tau_{3} = 0.01\). For the corresponding carrying capacity of the susceptible individuals \(q = 50, q = 75,q = 95\), the solutions \(S\left( t \right), V\left( t \right)\) and \(I\left( t \right)\) of the system (2) are presented in Fig. 6. Vaccinated population increases as carrying capacity increases. That implies disease infection can be controlled.

Fig. 6
figure 6

Effect of all the compartment \(\left( {S\left( t \right), V\left( t \right), I\left( t \right)} \right)\) for the parameter \(q\)

Full size image

Figure 7 depicts the effect of all compartment for the different values of disease transmission rate. For the stochastic perturbation \(\left( {\tau_{1} ,\tau_{2} ,\tau_{3} } \right) = \left( {0.08, 0.02,0.01} \right)\), the ecological parameters \(\left( {p, q,k_{1} ,\psi , \alpha_{1} , \alpha_{2} , \beta_{1} ,\beta_{2} } \right) = \left( {2, 50, 0.65, 0.02, 0.9, 0.03, 0.2, 0.05} \right)\) and the disease transmission rate \(\lambda = 0.9, 0.09,0.009\) the corresponding population intensities of susceptible, vaccinated and infected are presented in Fig. 7. The infection will decrease as the disease transmission rate increases.

Fig. 7
figure 7

Impact of all the compartment for the parameter \(\lambda\)

Full size image

Choose the epidemiological parameter \(\left( {p, q,\lambda ,\psi , \alpha_{1} , \alpha_{2} , \beta_{1} ,\beta_{2} } \right) = \left( {2, 50, 0.9, 0.02, 0.9, 0.03, 0.2, 0.05} \right)\) and the random noises \(\left( {\tau_{1} ,\tau_{2} ,\tau_{3} } \right) = \left( {0.08, 0.02,0.01} \right)\). Consider the saturation constant \(k_{1} = 0.65, 0.85, 1.05\), the corresponding population number of susceptible, vaccinated and infected individuals is represented in Fig. 8. The saturation rate is inversely proportional to infected individuals, i.e., the saturation increases, and then, the epidemic infection decreases.

Fig. 8
figure 8

Impact of all the compartment for the parameter \(k_{1}\)

Full size image

Discussion of results

In this study, to the best of our knowledge, we establish the extinction of disease and disease persistence of the SVIS epidemic model, which adds the existence of ergodic stationary distribution and analysis of probability density function. The main study will be delivered in the following ways.

  • With the linear random fluctuation (Cai and Kang 2015; Zhao and Jiang 2014; Khan and Khan 2018; Zhang 2017; Caraballo et al. 2020; Wang and Jiang 2019; Wang and Wang 2018; Liu et al. 2019a; Zhou et al. 2020, 2021; Qi and Jiang 2020; Li et al. 2010), we emphasize on SVIS epidemic model in stochastic nature with saturated incidence and vaccination strategies. With the help of Lyapunov functions, we derive the stochastic critical value \(R_{S}^{*}\). Under \(R_{S}^{*} > 1\), a unique ergodic stationary distribution is obtained of the system (2) by the Khas’minskii theory. As the same expression of \(R_{S}^{*}\) and basic reproduction number \(R^{*}\), the dynamical properties of susceptible, vaccinated, infected individuals determine the stochastic positive equilibrium state. In the present study, \(\tau_{1} , \tau_{2} , \tau_{3}\) are their corresponding random fluctuations. By using numerical results and sensitivities of parameters, we discuss some important measure to control the infectious disease.

  • We cannot determine the clear view of statistical nature of disease persistence from the existence of ergodic stationary distribution. With the help Zhou et al. (Zhou et al. 2020), some algebraic equations are improved by the three-dimensional probability density function. The existence of stationary distribution with ergodic properties can obtain the corresponding persistence in the studies (Ma et al. 2015; Liu et al. 2019b). We obtain the expression of log-normal three-dimensional density function \({\mathcal{T}}\left( {h_{1} ,h_{2} ,h_{3} } \right).\) Using the algebraic equation \(Y^{2} + B\rho + \rho B^{T} = 0\), we solve the covariance matrix \(\rho .\) It is more difficult to obtain \(\rho\). Then, some standard matrices \(R_{1} ,R_{2} ,R_{3}\) are taken. We can investigate the matrix \(\rho\) is positive definite and the diffusion matrix \(W_{0}\) is semi-positive definite.

For our future scope, we should add another population for recovered and analyze the disease persistence and the existence of ergodic stationary distribution. By considering the important effect of telegraph noises (Zhang et al. 2016) into the system and regime switching can be studied. It is expected that this task can be solved in later research work.