Stochastic Analysis for the Dynamics of a Poliovirus Epidemic Model

1  Introduction

Jenkins et al. in 2006 formulated a model in which he concluded that the use of monovalent is better than other vaccines. It provides outstanding outbreak control [1]. Haldar et al. [2] introduced the poliovirus vaccine in India. Kalkowska et al. in 2020 represent a differential equation-based stochastic model for poliovirus transmission. The model shows the poliovirus transmission for 2019 to 2023 with a strategic eradication plan [3]. Minor studied the types of polioviruses, vaccination, and eradication of the virus worldwide [4]. Thompson [5] investigated the transmission dynamics of the poliovirus in Nigeria. Duque-Marin et al. [6] studied two types of vaccines in the mathematical model. Denes et al. [7] presented a model which describes polio transmission in tropical regions. Cheng et al. [8] discussed a polio vaccination model in two different age classes. Alba et al. [9] addressed the correlation between climate and vaccination through a mathematical model. Shaghaghi et al. in 2018, studied that the OPV and IVPPvs vaccine was helpful for the eradication of the virus last few years [10]. Shimizu in 2014 [11] explained IPV is very effective against the poliovirus, and the author reviewed the introduction, development, and characterization of the OPV vaccine. In addition, his place in the world was told. Rafique et al. in 2020 presented a mathematical model in which they discovered the dynamics of poliovirus transmission using standard methods with vaccination [12]. Nidia et al. [13] in 2007, examined the effects of the poliovirus on human life and the steps taken to eradicate the virus and discussed what steps we could take in the future to get rid of it. Thompson et al. [14] presented polio outbreaks in the USA. Kalkowska et al. introduced a model to identify poliovirus and opportunities to increase population immunity [15]. Kim et al. [16] presented a model to examine the transmission of virulent circulating vaccine-derived polioviruses. Hillis [17] formulated a model in different regions before using artificial poliovirus vaccination. Mendrazitsky et al. [18] explained a disease of epidemic development model. The model analyzed other properties of polio and its non-equilibrium outbreak dynamics. Debanne et al. [19] presented a mathematical model of poliovirus in America. Naik et al. [20,21] studied the fractional modeling of cancer and HIV infection with the well-known results of stabilities.

The strategy of the paper is as follows. The first section goes to literature, and Section 2 goes to stochastic modeling of poliovirus and its fundamental properties. Section 3 goes to the proposed numerical method and its simulation with current approaches in the literature. Section 4 goes to the paper’s conclusion and remarks.

2  Poliovirus Model

For any time t, S: represents the class that is influenced by infection, E: represents the class that is disclosed by infection, I: represents an infective class, V: represents immunization class, A: represents the constant immigration rate of the human population. β: is the per unit time probability of infection transmission by the infective population. r: is the reduction in the exposed class due to transmission of infection. v: represents the proportion of recruits in the susceptible class moving to the vaccinated class, v1: is the number of vaccinated exposed populations, b: number of exposed populations moving to the infection class. µ: natural death of the human population, α: disease death rate. The first order, nonlinear, and coupled ordinary differential equations of the poliovirus epidemic model are assumed as follows:

dS=(A−βSI−rβSE−(μ+v)S)dt+σ1SdB(t).(1)

dE=(βSI+rβSE−(b+μ+v1)E)dt+σ2EdB(t).(2)

dI=((b+v1)E−(μ+α)I)dt+σ3IdB(t).(3)

dV=(vS−μV)dt+σ4VdB(t).(4)

with initial condition S (0) ≥ 0; E (0) ≥ 0; I (0) ≥ 0; V (0) ≥ 0, and (σi:i=1,2,3,4) is the peutrbation term with B(t) is the Brownian motion [22,23].

2.1 Properties [24]

This section studies the positivity and boundedness of the system ((1)–(4)). Let us consider the vector as follows:

U(t)=(S(t),E(t),I(t),V(t)).(5)

And the norm |U(t)|=S2(t)+E2(t)+I2(t)+V2(t).

dU(t)=H(U,t)dt+K(U,t)dW(t).(6)

As, L=∂∂t+∑i=14Hi(U,t)∂∂ui+12∑i,j=14(KT(U,t)K(U,t))i,j×∂2∂Ui∂Uj, where L is differential operator.

If L acts on a function V1∈C2,1(R4×(0,∞);R+) then we denote

LV(U,t)=V1t(U,t)+V1U(U,t)H(U,t)+12Trace(KT(U,t)VUU(U,t)K(U,t)).

where transportation is denoted by T.

Theorem 1: For model ((1)–(4)) and any given initial value (S(0),E(0),I(0),V(0))∈R+4, there is a unique solution (S(t),E(t),I(t),V(t)) on t≥0 and will remain in R+4 with probability one.

Proof: By Ito’s formula, the model ((1)–(4)) admits a positive solution in the unique local on [0,τe], and explosion time is denoted by τe. Because the local Lipschitz condition is satisfied by all the coefficients of the model as mentioned earlier.

Next, let us show that the given model ((1)–(4)) admits this solution in the global sense; that is, τe=∞ almost sure.

Let mo=0 be sufficiently large for S(0),E(0),I(0), and V(0) lying with the interval [1mo,mo]. For each integer m≥mo, define a sequence that is so-called stopping times as

τm=inf{tϵ[0,τe]:S(t)∉(1m,m)orE(t)∉(1m,m)orI(t)∉(1m,m)orV(t)∉(1m,m)}.(7)

where we set infϕ=∞(ϕ represents the empty set). Since τm is non-decreasing as m→∞,

τ∞=limm→∞τm.(8)

Then τ∞≤τe. To prove, τ∞=∞.

In case of violation of statement, then T>0 and εϵ(0,1) such that

P{τ∞≤T}>ε.(9)

this, there is an integer m1>m0 such that

P{τm≤T}≥ε,∀m≥m1.(10)

Define a C4− function V:R+4→R+ by

V1(S,E,I,V)=(S−1−ln⁡S)+(E−1−ln⁡E)+(I−1−ln⁡I)+(V−1−ln⁡V).(11)

By using Ito’s formula, we calculate

dV1(S,E,I,V)=(1−1S)dS+(1−1E)dE+(1−1I)dI+(1−1V)dV+σ22dt.

dV1(S,I,R)=(1−1S)((A−βSI−rβSE−(μ+v)S)dt+σ1SdB(t))+(1−1E)×((βSI+rβSE−(b+μ+v1)E)dt+σ2EdB(t))+(1−1I)×(((b+v1)E−(μ+α)I)dt+σ3IdB(t))+(1−1R)[(vS−μV)dt+σ4VdB(t)]dt+σ22dt.

dV1(S,E,I,V)≤[A+μ+v+α+σ22]dt+σ1SdB(t)++σ2EdB(t)+σ3IdB(t)++σ4VdB(t).(12)

Let, N=A+μ+v+α+σ22, Then Eq. (12) could be written as

dV1(S,E,I,V)≤Ndt+σ1SdB(t)++σ2EdB(t)+σ3IdB(t)++σ4VdB(t).(13)

By integrating from 0 to τm∧τ, we get

∫0τm∧τdV1(S(s),E(s),I(s),V(s))≤∫0τm∧τNds+∫0τm∧τ(σ1S+σ2E+σ3I+σ4V)dB(s).(14)

where τm∧τ=min(τm,T), the taking the expectations to lead to

EV1(S(τm∧τ),I(τm∧τ),E(τm∧τ),V(τm∧τ))≤V(S(0),E(0),I(0),V(0))+NT(15)

Set Ωm={τm≤T} for m>m1, and from (15), we have P(Ωm≥ε). For every v∈Ωm there are some i such that ui(τm,v) equals either m or 1m for i=1,2,3;

Hence, V1(S(τm,v),E(τm,v),I(τm,v),V(τm,v)) is less than min{m−1−ln⁡m,1m−1−ln⁡1m}.

Then we obtain

V1(S(0),E(0),I(0),V(0))+NT≥E(IΩm(v)V1(S(τm),E(τm),I(τm),V(τm)))≥{min(m−1−ln⁡m,1m−1−ln⁡1m)}.(16)

IΩm(v) of Ωm represents the indicator function. Letting m→∞ leads to the contradiction ∞=V1(S(0),E(0),I(0),V(0))+NT<∞.

As desired.

2.2 Equilibria of Model

The disease-free equilibrium of the model is K 0=(Aμ+v,0,0,vμS).

The endemic equilibrium of the model is denoted by K1 = (S*, E*, I*, V*).

S∗=(b+μ+v1)(μ+α)β(b+v1)+rβ(μ+α), E∗=A−(μ+v)S∗rβS∗+βS∗(b+v1μ+α), I∗=b+v1μ+αE∗, and V∗=vμS∗.

2.3 Local Stability of Model

Theorem 2: The disease-free equilibrium K 0=(Aμ+v,0,0,vμS) is locally asymptotically stable if R0<1; otherwise, unstable if R0>1.

Proof: Considering the function from the system ((1)–(4)) as follows:

F1=A−βSI−rβSE−(μ+v)S.

F2=βSI+rβSE−(b+μ+v1)E.

F3=(b+v1)E−(μ+α)I.

F4=vS−μV.

The elements of the given Jacobean Matrix at KO is as follows:

J(KO)=[−(μ+v)−rβAμ+V−βAμ+V00rβAμ+V−(b+μ+v1)βAμ+V00b+v1−(μ+α)0V00−μ].

J(KO−λI)=|−(μ+v)−λ−rβAμ+V−βAμ+V00rβAμ+V−(b+μ+v1)−λβAμ+V00b+v1−(μ+α)−λ0V00−μ−λ|=0.

−(−μ−λ)|−(μ+v)−λ−rβAμ+V−βAμ+V0rβAμ+V−(b+μ+v1)−λβAμ+V0b+v1−(μ+α)−λ|=0.

λ1=−μ<0,

|−(μ+v)−λ−rβAμ+V−βAμ+V0rβAμ+V−(b+μ+v1)−λβAμ+V0b+v1−(μ+α)−λ|=0.

−(μ+v)−λ|rβAμ+V−(b+μ+v1)−λβAμ+Vb+v1−(μ+α)−λ|=0.

λ2=−(μ+v)<0.

λ2+λ[(μ+α)−rβAμ+v+λ2]−rβAμ+v(μ+α)−βAμ+v(b+v1)=0.

λ2+λ[(μ+α)+(b+μ+v1)(1−R2)]+(μ+α)(b+μ+v1)(1−R0)=0.

We obtain the following results by applying Routh Hurwitz criteria for 2nd order.

A1>0,A2>0, if R0,R2<1, where A1=[(μ+α)+(b+μ+v1)(1−R2)], A2=(μ+α)(b+μ+v1)(1−R0).

Theorem 3: The endemic equilibrium K1 = (S*, E*, I*, V*) is locally asymptotically stable if Ro>1.

Proof: The Jacobean matrix at K1=(S∗,E∗,I∗,V∗) is as follows:

J(K1)=[−βI∗−rβE∗−(μ+v)−rβS∗−βS∗0βI∗+rβE∗rβS∗−(b+μ+v1)βS∗00b+v1−(μ+α)0V00−μ].

|J(K1−λI)|=|−βI∗−rβE∗−(μ+v)−λ−rβS∗−βS∗0βI∗+rβE∗rβS∗−(b+μ+v1)−λβS∗00b+v1−(μ+α)−λ0V00−μ−λ|=0.

−(−μ−λ)|−βI∗−rβE∗−(μ+v)−λ−rβS∗−βS∗βI∗+rβE∗rβS∗−(b+μ+v1)−λβS∗0b+v1−(μ+α)−λ|=0.

λ1=−μ<0,

|−βI∗−rβE∗−(μ+v)−λ−rβS∗−βS∗βI∗+rβE∗rβS∗−(b+μ+v1)−λβS∗0b+v1−(μ+α)−λ|=0

−βI∗−rβE∗−(μ+v)−λ|rβS∗−(b+μ+v1)−λβS∗b+v1−(μ+α)−λ|−(βI∗+rβE∗)|−rβS∗−βS∗b+v1−(μ+α)−λ|=0.

−βI∗−rβE∗−(μ+v)−λ[(rβS∗−(b+μ+v1)−λ)(−(μ+α)−λ)−(βS∗)(b+v1)]−(βI∗+rβE∗)[(−rβS∗)(−(μ+α)−λ)−(−βS∗)(b+v1)]=0.

−βI∗−rβE∗−(μ+v)−λ[−rβS∗(μ+α)−rβS∗λ+(b+μ+v1)(μ+α)+(b+μ+v1)λ+(μ+α)λ+λ2−(βS∗)(b+v1)]−βI∗−rβE∗[rβS∗(μ+α)+rβS∗λ+βS∗b+v1]=0.

rβ2S∗I∗(μ+α)+rβ2S∗I∗λ−βI∗(b+μ+v1)(μ+α)−βI∗(b+μ+v1)λ−βI∗(μ+α)λ−βI∗λ2+β2S∗I∗(b+v1)+r2β2E∗S∗(μ+α)+r2β2E∗S∗λ−rβE∗(b+μ+v1)(μ+α)−rβE∗(b+μ+v1)λ−rβE∗(μ+α)λ−rβE∗λ2+rβE∗S∗(b+v1)+rβS∗(μ+α)(μ+v)+rβS∗(μ+v)λ−(b+μ+v1)(μ+α)(μ+v)−(μ+v)(b+μ+v1)λ−λ(μ+α)(μ+v)−(μ+v)λ2+(μ+v)βS∗(b+v1)+rβS∗(μ+α)λ+rβS∗λ2−(b+μ+v1)(μ+α)λ−(b+μ+v1)λ2−(μ+α)λ2−λ3+(βS∗)(b+v1)λ−rβ2S∗I∗(μ+α)−rβ2S∗I∗λ−β2S∗I∗(b+v1)−r2β2E∗S∗(μ+α)−r2β2E∗S∗λ−rβ2S∗E∗(b+v1)=0.

−λ3−λ2[βI∗+rβE∗+(μ+v)−rβS∗+(b+μ+v1)+(μ+α)]−λ[−rβ2S∗I∗+βI∗(b+μ+v1)+βI∗(μ+α)+rβE∗(b+μ+v1)+rβE∗(μ+α)−rβS∗(μ+v)+(μ+v)(b+μ+v1)+(μ+α)(μ+v)−rβS∗(μ+α)+(b+μ+v1)(μ+α)−(βS∗)(b+v1)+rβ2S∗I∗]−[−rβ2S∗I∗(μ+α)+βI∗(b+μ+v1)(μ+α)−r2β2E∗S∗(μ+α)+rβE∗(b+μ+v1)(μ+α)−rβE∗S∗(b+v1)−rβS∗(μ+α)(μ+v)+(b+μ+v1)(μ+α)(μ+v)−(μ+v)βS∗(b+v1)+rβ2S∗I∗(μ+α)+r2β2E∗S∗(μ+α)+rβ2S∗E∗(b+v1)]=0.

−[λ3+λ2A+λC+D]=0.

λ3+λ2A+λC+D=0.

where A= βI∗+rβE∗+(μ+v)−rβS∗+(b+μ+v1)+(μ+α),

B=−rβ2S∗I∗+βI∗(b+μ+v1)+βI∗(μ+α)+rβE∗(b+μ+v1)+rβE∗(μ+α)−rβS∗(μ+v)+(μ+v)(b+μ+v1)+(μ+α)(μ+v)−rβS∗(μ+α)+(b+μ+v1)(μ+α)−(βS∗)(b+v1)+rβ2S∗I∗,

C=−rβ2S∗I∗(μ+α)+βI∗(b+μ+v1)(μ+α)−r2β2E∗S∗(μ+α)+rβE∗(b+μ+v1)(μ+α)−rβE∗S∗(b+v1)−rβS∗(μ+α)(μ+v)+(b+μ+v1)(μ+α)(μ+v)−(μ+v)βS∗(b+v1)+rβ2S∗I∗(μ+α)+r2β2E∗S∗(μ+α)+rβ2S∗E∗(b+v1),

D=−rβ2S∗I∗(μ+α)+βI∗(b+μ+v1)(μ+α)−r2β2E∗S∗(μ+α)+rβE∗(b+μ+v1)(μ+α)−rβE∗S∗(b+v1)−rβS∗(μ+α)(μ+v)+(b+μ+v1)(μ+α)(μ+v)−(μ+v)βS∗(b+v1)+rβ2S∗I∗(μ+α)+r2β2E∗S∗(μ+α)+rβ2S∗E∗(b+v1).

Applying Routh-Hurwitz Criterion for 3rd order, A>0,D>0, and AC>D, if R0>1.

Hence the given system is locally asymptotically stable.

2.4 Reproduction Number

The idea of reproduction number is presented in [25] by considering Eqs. (3) and (4), we get the following matrices:

[E′I′V′]=[rSββS0000000][EIV]−[b+μ+v100−b−v1μ+α000μ][EIV].

Here, F =[rSββS0000000] and G =[b+μ+v100−b−v1μ+α000μ]

FG−1=[rβA(μ+v)(b+μ+v1)+βA(b+v1)(μ+v)(μ+α)(b+μ+v1)βA(μ+v)(μ+α)0000000].

The spectral radius of the FG−1 is called the reproduction number is as follows:

R0=R1+R2

R0=βA(b+v1)(μ+v)(μ+α)(b+μ+v1)+rβA(μ+v)(b+μ+v1)

where R1=βA(b+v1)(μ+v)(μ+α)(b+μ+v1) , R2=rβA(μ+v)(b+μ+v1).

3  Stochastic Poliovirus Epidemic Model

Let us consider the vector C=[S,E,I,V]T of stochastic differential equations (SDEs) of the poliovirus epidemic model ((1)–(4)). We want to calculate the expectation and variance (see Table 1).

Drift=G(C,t)=[A−βSI−rβSE−μS−VSβSI+rβSE−bE−μE−v1EbE+v1E−μI+αIA−μV]Δt.

Diffusion=[A+βSI+rβSE+μS+vS−βSI−rβSE0−vS−βSI−rβSEβSI+rβSE+bE+μE+v1E−bE−v1E00−bE−v1EbE+v1E+μI+αI0−vS00vS+μV]

The equation of poliovirus epidemic model ((1)–(4)) can be written as

d[SEIV]=[A−βSI−rβSE−μS−VSβSI+rβSE−bE−μE−v1EbE+v1E−μI+αIA−μV]dt+[A+βSI+rβSE+μS+vS−βSI−rβSE0−vS−βSI−rβSEβSI+rβSE+bE+μE+v1E−bE−v1E00−bE−v1EbE+v1E+μI+αI0−vS00vS+μV]dB.(17)

The Euler Maruyama approach is cast-off to determine the numerical result of the Eq. (17) by using the values of the parameters given in Table 2 as follows:

Cn+1=Cn+f(Cn,t)Δt+L(Cn,t)dB.

[Sn+1En+1In+1Vn+1]=[SnEnInVn]+[[A−βSnIn−rβSnEn−μSn−vSnβSnIn+rβSnEn−bEn−μEn−v1EnbEn+v1En−μIn+αInA−μVn]]Δt+[A+βSnIn+rβSnEn+μSn+vSn−βSnIn−rβSnEn0−vSn−βSnIn−rβSnEnβSnIn+rβSnEn+bEn+μEn+v1En−bEn−v1En00−bEn−v1EnbEn+v1En+μIn+αIn0−vSn00vSn+μVn]ΔBn(18)

where C(0)=Co=[0.5,0.3,0.2,0.1,]T,0≤t≤C and  Brownian motion is denoted as B.

3.1 Stochastic Nonstandard Finite Difference Method

The stochastic NSFD can be developed for the system ((1)–(4)) as

dSdt=A−βSI−rβSE−(μ+v1)S.

The breakdown of the proposed method for the above equation.

Sn+1=Sn+h[A−βSn+1In−rβSn+1En−(μ+v1)Sn+1+σ1SnΔB1].

Sn+1=Sn+hA+hσ1SnΔB11+hβIn+hrβEn+h(μ+v1).(19)

Similarly, we break the remaining system into a proposed method like (19), as follows:

En+1=En+hβSnIn+hrβSnEn+hσ2EnΔB21+h(b+μ+v1).(20)