Numerical Analysis of Bacterial Meningitis Stochastic Delayed Epidemic Model through Computational Methods

1  Introduction

Meningitis caused by bacteria results in inflammation of the membranes surrounding the brain and spinal cord. The most common bacteria that cause meningitis are Streptococcus pneumoniae, Neisseria meningitidis, and Haemophilus influenzae type b. The signs and symptoms of bacterial meningitis include a sharply rising temperature, headache, stiff neck, and light sensitivity. Serious complications such as brain damage or death are caused by bacterial meningitis, so swift medical attention is crucial. Specific bacterial strains that cause meningitis can be controlled by effective vaccination [1]. Medical officers usually prescribe medicines such as ceftriaxone or cefotaxime to target invading bacteria and cure bacterial meningitis. Corticosteroids such as dexamethasone may be dispensed to decrease inflammation. Shin et al. [2] studied a mathematical model incorporating hospitalization and treatment factors. Madaki et al. [3] designed a mathematical model of meningococcal meningitis in Nigeria. In [4,5], the study was based on acute bacterial meningitis in adults. In [6,7], the study was designed on reviewed reports from the WHO and a model that describes a bacterial disease. Mohanty et al. [8] explained the different types of vaccines and preventive measures to control bacterial meningitis in Sweden. Van et al. [9] investigated in great detail all the stages of bacterial meningitis, from its initiation to its outbreak, and formulated various preventive measures to control it. Rauti et al. [10] developed a mathematical model for preventing bacterial meningitis in animals that can be employed to prevent infections of this deadly disease. Hou et al. [11] presented a dynamic model for preventing bacterial meningitis in young children using metagenomic sequencing, in which this approach can largely control the disease. Signing et al. [12] created a mathematical model that can visualize the changes caused by bacterial meningitis disease and their effect and can be treated accordingly to control the disease. Afridi et al. [13] identified different preventive measures and vaccines while reviewing bacterial meningitis cases brought to a civil hospital in Pakistan. Buonomo et al. [14] reported that the meningitis outbreak in Nigeria is analyzed using a behavior-focused vaccination model. This strategy considers behavioral elements while evaluating immunization strategies. It contributes to developing more effective preventative methods tailored to Nigeria’s socioeconomic environment and sheds light on the dynamics of meningitis transmission. Babatunde et al. [15] employed an online tool using logistic regression to predict cerebrospinal meningitis incidence is currently being developed. The goal of this initiative is to enhance illness management and early diagnosis. The initiative uses statistical modeling to evaluate relevant data, which can give useful information to policymakers and healthcare practitioners. Bradley et al. [16] reviewed the English medical literature and introduced a method that can be utilized to control bacterial meningitis. Choi et al. [17] created a mathematical model using AI to control bacterial meningitis in the future. Chen et al. [18] stated that it is essential to identify clinical and laboratory parameters to control bacterial meningitis disease. Olu et al. [19] created a mathematical model that can be very useful to control the growing epidemic of bacterial meningitis in the population. In [20], the authors studied a new analyzing technique for nonlinear time fractional Cauchy reaction-diffusion model equations. The stochastic epidemic models with delays are constructed given the need to incorporate the proper transmission process of the epidemic in any populace. These models make use of such lags and lags to mimic the period within which a person gets exposed and becomes infectious and lags in the uptake and efficacy of the implemented interventions. These models help to explain why so many epidemics exhibit cyclic behavior, this is in terms of the waves and seasonality that characterize many epidemics, and delay, noise, and heterogeneity are shown to play crucial roles. Also, the stochastic with delay effect results close to the actual statistics of the problem. So, in the present study, we extend the deterministic model in the sense of stochastic with delay. Then our new system was based on stochastic delay differential equations (SDDEs), due to complexity, we used standard and nonstandard methods to understand the dynamics of bacterial meningitis.

The structure of the paper is as follows: a brief overview and thorough analysis of bacterial meningitis-like disorders found in the literature are provided in Section 1. Establishing the delayed model and the ensuing mathematical analysis are the focus of Section 2. Section 3 presents the stochastic formulation of the model and the extinction and persistence of the model. The nonstandard computational method with delay and stochastic is introduced in Section 4, along with its global and local stability assessment. Numerical simulations and the presentation of the results are the explicit focus of Section 5. The final opinions provide a comprehensive conclusion of the work in Section 6.

2  Formulation of Model

The dynamics of population theory is a foundation for the model’s development. The classes susceptible S(t), carrier C(t), infected I(t), and recovered R(t) add up to the population N(t), and they use the law of mass action to represent the condition of bacterial meningitis (Fig. 1).

Figure 1: Flow of bacterial meningitis [3]

The parameter values of the constants are shown as follows in Table 1.

The law of mass action specifies the continuous model based on the presumptions. The transmission flow of bacterial meningitis-type disease is determined using the nonlinear delayed differential equations (DDEs) as follows:

S′(t)=Λ−βSIe−µτN−(µ+ρ)S(t)+ωR(t).(1)

C′(t)=βSIe−µτN−(γ+µ)C(t).(2)

I′(t)=γC(t)−(τ1+δ+μ)I(t).(3)

R′(t)=τ1I(t)+ρS(t)−(ω+µ)R(t).(4)

with the initial assumptions S(0)≥0,C(0)≥0,I(0)≥0,R(0)≥0, as well as the constraint that t≥0 and τ≤t.

2.1 Model Properties

For the critical analysis of the model to continue, all of the variables S(t),C(t),I(t), and R(t) must be non-negative. As a result, every time t that is greater than or equal to 0 and less than or equal to τ inside a feasible zone is covered by the conclusions drawn from the model’s analysis.

ℋ={(S,C,I,R)ϵR+4:N(t)≤Λµ,S≥0,C≥0,I≥0,R≥0}.

2.2 Analysis of Model

This section presents a concise analysis of the equilibria in the delayed model for bacterial meningitis, i.e., Meningitis-Free Equilibrium (MFE−D0), and Meningitis Existing Equilibrium (MEE−D∗). Hence,

D0=(S0,C0,I0,R0)=(Λ(ω+µ)[(ω+µ)(µ+ρ)−ωρ],0,0,Λρ[(ω+µ)(µ+ρ)−ωρ]), and D∗=(S∗,C∗,I∗,R∗).

S∗=N[(Λ+τ1I∗)βγe−µτ+ρN(µ+γ)(τ1+δ+μ)]βγe−μτ[βγe−μτ+N(µ+ρ)], C∗=(τ1+δ+μ)I∗γ,I∗=Λβγe−μτ−Nμ(µ+γ)(τ1+δ+μ)(µ+γ)(τ1+δ+μ)−τ1βγe−μτ,

R∗=τ1βγe−μτ+ρN(µ+γ)(τ1+δ+μ)βγe−μτ.

2.3 Reproduction Number

The systems (1)–(4) incorporate the advanced matrix approach to compute the threshold number R0. Neglected S′, the following can be obtained:

[C′I′R′]=[0βSe−μτN0000000][CIR]−[(µ+γ)00−γ(τ1+δ+μ)00−τ1(ω+µ)][CIR]

F=[0βSe−μτN0000000],G=[(µ+γ)00−γ(τ1+δ+μ)00−τ1(ω+µ)].

FG−1|D0=[βγS0e−μτN(µ+γ)(τ1+δ+μ)βS0e−μτN(τ1+δ+μ)0000000].

The spectral radius of FG−1|D0, which is also referred to as the threshold number, is defined as follows:

ℛ0=Λ(ω+µ)βγe−μτ(µ+γ)(τ1+δ+μ)[(ω+µ)(µ+ρ)−ωρ].(5)

3  Transition Probabilities of the Delayed Model

Let us examine the vector U(t)=[S(t),C(t),I(t),R(t)]T, along with the corresponding event probabilities as outlined in Table 2.

Here,

E∗[ΔU]=∑i=110PiMi=[Λ−βS(t)I(t)e−µτ−(µ+ρ)S(t)+ωR(t)βS(t)I(t)e−µτ−(γ+µ)C(t)γC(t)−(τ1+δ+μ)I(t)τ1I(t)+ρS(t)−(ω+µ)R(t)]Δt.

Var=E∗[ΔUΔUT]=∑i=110Pi[Mi][Mi]T.

=[P1+P2+P3+P4+P5−P20−P3−P5−P2P2+P6+P7−P600−P6P6+P8+P9−P8−P3−P50−P8P3+P5+P8+P10]Δt.

drift =H(U(t),t)=E∗[ΔU]Δt, diffusion =K(U(t),t)=E∗[ΔUΔUT]Δt, so

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

Eq. (6) is the stochastic differential equation, where B(t) is the Brownian motion.

3.1 Euler Maruyama Method

This section discusses a conventional numerical technique for approximating a stochastic delayed model’s result of (6). In this regard, Iq={0,1,2,3,…,q} is admitted for each q ∈ N. Let N∈N, and with the effect of time ∆t divides the partition into equal intervals [0, T] with constant delay. As follows:

0=to<t1<t2<…<tN=T,

for each n ∈ IN. Needless to mention tn=τn, for each n ∈ IN. Moreover, Un=U(tn) is accepted whenever n ∈ IN and U = S, E, I, R. Also, it can be set

ΔBn =B(tn+1) −B(tn),∀n∈ IN−1.

The mean of each ∆ Bn follows a normal distribution with a variance of one and an average of zero.

3.2 Stochastic Delayed Model

The system of stochastic delayed differential equations (SDDEs) is a mathematical model that describes the evolution of a set of variables over time, where the equations involve both deterministic time delays and stochastic (random) components. Where the stochastic term σi:(i=1,2,3,4)d(B(t)) introduces randomness into the system of differential equations as follows:

dS(t)=[Λ−βS(t)I(t)e−µτN−(µ+ρ)S(t)+ωR(t)]dt+σ1S(t)d(B(t)).(7)

dC(t)=[βS(t)I(t)e−µτN−(γ+µ)C(t)]dt+σ2C(t)d(B(t)).(8)

dI(t)=[γC(t)−(τ1+δ+μ)I(t)]dt+σ3I(t)d(B(t)).(9)

dR(t)=[τ1I(t)+ρS(t)−(ω+µ)R(t)]dt+σ4R(t)d(B(t)).(10)

Subject to non-negative initial conditions S(t)=S0≥0, C(t)=C0≥0,I(t)=I0≥0,R(t)=R0≥0.

3.2.1 Positivity and Boundedness

Let us examine a probability space represented as (Ω,F,P) and filtered by {Ft}tϵR. All P-null sets that fulfill the criteria of being both right-continuous and growing are encompassed inside Fo. Represent

U(t)=(S(t),C(t),I(t),R(t)).(11)

The norm |U(t)|=S2(t)+C2(t)+I2(t)+R2(t), and denote O3,1(R4×(0,∞);R+).

This study analyzes the collection of all non-negative functions V (U, t) defined on R4×(0,∞), which possess continuous second-order differentiability in U and first-order differentiability in t. Within this specific framework, the differential operator Q is precisely defined, which is associated with a three-dimensional stochastic differential equation.

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

As, Q=∂∂t+∑i=14Hi(U,t)∂∂ui+12∑i,j=14(KT(U,t)K(U,t))i,j×∂2∂Ui∂Uj

If Q acts on a function V∈O3,1(R4×(0,∞);R+), then

QV(U,t)=Vt(U,t)+VU(U,t)H(U,t)+12Trace(KT(U,t)VUU(U,t)K(U,t)),

where T means transportation.

Theorem 1. For a given model (7)–(10) and initial values (S(0),C(0),I(0),R(0))∈R+4, there exists a single solution (S(t),C(t),I(t),R(t)) on t≥0,τ≤t, and this solution will always remain in R+4 with a probability of one.

Proof: Based on Ito’s formula, the model (7)–(10) has a unique local positive solution in the interval [0,τe], and the time at which it explodes is represented by τe. Since all the coefficients in the model above satisfy the local Lipschitz criterion.

Next, this study demonstrates that the provided model (7)–(10) possesses a solution in the global sense, specifically when τe approaches infinity with practically certain probability.

Assume that mo=0 is a suitably high value such that S(0),C(0),I(0), and R(0) fall within the range [1mo,mo]. Let’s define a sequence of interrupting times for each integer m≥mo.

τm=inf{tϵ[0,τe]:S(t)∉(1m,m)orC(t)∉(1m,m)orI(t)∉(1m,m)orR(t)∉(1m,m)}.(13)

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

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

Thereafter, the limit τ∞ is less than or equal to τe nearly surely. Now, this study must demonstrate that τ∞ approaches infinity with almost certain probability.

If this condition is not met, then there are positive values T and ε, where T is more than zero, and ε is between 0 and 1, such that

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

there is an integer m1>m0 such that

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

Define a O4− function f:R+4→R+ by

f(S,C,I,R)=(S−1−ln⁡S)+(C−1−ln⁡C)+(I−1−ln⁡I)+(R−1−ln⁡R).(17)

The following can be calculated using Ito’s formula:

df(S,C,I,R)=(1−1S)dS+(1−1C)dC+(1−1I)dI+(1−1R)dR+σ12+σ22+σ32+σ422dt.

df(S,C,I,R)=(1−1S)((Λ−βS(t)I(t)e−µτN−(µ+ρ)S(t)+ωR(t))dt+σ1S(t)d(B(t)))+(1−1C)((βS(t)I(t)e−µτN−(γ+µ)C(t))dt+σ2C(t)d(B(t)))+(1−1I)((γC(t)−(τ1+δ+μ)I(t))dt+σ3I(t)d(B(t)))+(1−1R)(τ1I(t)+ρS(t)−(ω+µ)R(t))dt+σ4R(t)d(B(t))

df(S,C,I,R)≤[Λ+4μ+ρ+τ1+δ+γ+ω+σ12+σ22+σ32+σ422]dt+[σ1S(t)+σ2C(t)+σ3I(t)σ4R(t)]d(B(t)).(18)

For simplification, this study assumes

N1=Λ+4μ+ρ+τ1+δ+γ+ω+σ12+σ22+σ32+σ422, Then Eq. (18) could be written as

df(S,C,I,R)≤N1dt+[σ1S(t)+σ2C(t)+σ3I(t)σ4R(t)]d(B(t)).(19)

Assume N1 be a positive constant. Integrating from 0 to τm∧τ yields

∫0τm∧τdf(S(t),C(t),I(t),R(t))≤∫0τm∧τN1ds+∫0τm∧τ[σ1S(t)+σ2C(t)+σ3I(t)σ4R(t)]d(B(t)).(20)

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

Ef(S(τm∧τ),C(τm∧τ),I(τm∧τ),R(τm∧τ))≤f(S(0),C(0),I(0),R(0))+N1T.(21)

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

Hence, f(S(τm,l),C(τm,l),I(τm,l),R(τm,l)) is less than min{m−1−ln⁡m,1m−1−ln⁡1m}.

Thereafter, the following can be obtained:

f(S(0),C(0),I(0),R(0))+N1T≥E(IΩm(l)V(S(τm),C(τm),I(τm),R(τm)))≥{min(m−1−ln⁡m,1m−1−ln⁡1m)}.(22)

The indicator function is represented by IΩm(l) of Ωm. Letting m→∞ leads to the contradiction

∞=f(S(0),C(0),I(0),R(0))+N1T<∞.

as desired.

3.2.2 Extinction and Persistence

Definition: Consider B(t) as a Brownian motion and I(t) as an Ito drift-diffusion process that obeys the stochastic differential equation:

dI(t)=[γC(t)−(τ1+δ+μ)I(t)]dt+σ3I(t)d(B(t)).

If f(I,t)∈C2(ℛ2,ℛ) then f(I(t),t) is also an Ito drift-diffusion process, which satisfies as follows:

d(f(I(t),t))=∂f∂t(I(t),t)dt+f′((I(t),t))dB(t)+12f′′((I(t),t))dB(t)2.

Let us introduce ℛ0S=ℛ0d−σ322(τ1+δ+μ).

Theorem 2. If ℛ0S < 1 and σ32<Λ(ω+µ)βγSe−μτ(µ+γ)(τ1+δ+μ)[(ω+µ)(µ+ρ)−ωρ], then the number of infected individuals in the system will exponentially approach zero (limt→∞I(t)=0).

Proof: Given the initial data (S(0),C(0),I(0),R(0))∈R+4, this study determines if the system (7)–(10) fulfills the stochastic differential equation and allows for a solution in the form of (S(t),C(t),I(t),R(t)).

dI(t)=[γC(t)−(τ1+δ+μ)I(t)]dt+σ3cI(t)d(B(t)).

Let, f(I)=ln⁡(I), then

dln(I)=f′(I)dI+12f′′(I)I2σ2dt.

dln⁡(I)=1IdI+12(−1I2)I2σ2dt.

dln⁡(I)=(γC(t)I(t)−(τ1+δ+μ)−12σ32)dt+σ3cd(B(t)).

ln⁡(I)=ln⁡I(0)+∫0t(γC(t)I(t)−(τ1+δ+μ)−12σ32)dt+∫0tσ3cd(B(t)),

Notice that, W(t)=∫0tσ3cd(B(t)) with W(0)=0.

Ifσ32>Λ(ω+µ)βγSe−μτ(µ+γ)(τ1+δ+μ)[(ω+µ)(µ+ρ)−ωρ],

ln⁡(I)>(Λ(ω+µ)βγSe−μτ(µ+γ)(τ1+δ+μ)[(ω+µ)(µ+ρ)−ωρ]−(τ1+δ+μ)−12Λ(ω+µ)βγSe−μτ(µ+γ)(τ1+δ+μ)[(ω+µ)(µ+ρ)−ωρ])t+W(t)+ln⁡I(0),

ln⁡(I)t>(Λ(ω+µ)βγSe−μτ2(µ+γ)(τ1+δ+μ)[(ω+µ)(µ+ρ)−ωρ]−(τ1+δ+μ))+W(t)t+ln⁡I(0)t,

limt→∞ln⁡(I)t>(Λ(ω+µ)βγSe−μτ2(µ+γ)(τ1+δ+μ)[(ω+µ)(µ+ρ)−ωρ]−(τ1+δ+μ))>0,withlimt→∞W(t)t=0,

Ifσ32<Λ(ω+µ)βγSe−μτ(µ+γ)(τ1+δ+μ)[(ω+µ)(µ+ρ)−ωρ], then

ln⁡(I(t))<(Λ(ω+µ)βγSe−μτ(µ+γ)(τ1+δ+μ)[(ω+µ)(µ+ρ)−ωρ]−(τ1+δ+μ)−12σ32)t+W(t)+ln⁡I(0),

ln⁡(I)t<(τ1+δ+μ)(Λ(ω+µ)βγSe−μτ(µ+γ)(τ1+δ+μ)2[(ω+µ)(µ+ρ)−ωρ]−1)+W(t)t+ln⁡I(0)t,

limt→∞supln⁡(I)t<(τ1+δ+μ)(ℛ0S−1), when ℛ0S<1, then limt→∞supln⁡(I)t≤0 can be obtained,

limt→∞I(t)=0, as desired.

4  Nonstandard Computational Method

The proposed non-parametric perturbation model’s first Eq. (10) can be stated using an unconventional computing technique; specifically, Eqs. (10)–(13) could be solved using a stochastic nonstandard finite difference technique.

dS(t)=[Λ−βS(t)I(t)e−µτN−(µ+ρ)S(t)+ωR(t)]dt+σ1S(t)d(B(t)).(23)

For the stochastic NSFD approach, the structure of the equation is as follows:

Sn+1−Snh=[Λ−βSnIne−µτN−(µ+ρ)Sn+ωRn+σ1SnΔBn](24)

The system (10)–(13) can be broken down using the stochastic NSFD process, and the complete system can then be written as follows:

Sn+1=Sn+h[Λ+ωRn+σ1SnΔBn]1+h(βIne−µτ+(µ+ρ))(25)

Cn+1=Cn+h[βSnIne−µτ+σ2CnΔBn]1+h(γ+μ)(26)

In+1=In+h[γCn+σ3InΔBn]1+h(τ1+δ+μ)(27)

Rn+1=Rn+h[τ1In+σ4RnΔBn]1+h(ω+μ)(28)

where n=0,1,2,…, and the discretization gap is denoted by “h.”

4.1 Convergence Analysis

Theorem 3. There is only one positive solution (S,C,I,R)∈R4,∀n>0 for any initial value (S(0),C(0),I(0),R(0))∈R4 for Eq. (25) through Eq.(28) with ΔBn=0.

Proof: The proof is easily verifiable due to the non-negative nature of the restriction of biological problems.

Theorem 4. For the region ℋ={(Sn,Cn,In,Rn)∈ℛ+4:Sn+Cn+In+Rn=N≤Λµ,Sn≥0,Cn≥0,In≥0,Rn≥0}. For every, n≥0 is an area of equations that is feasible and positive invariant (25) to (28) with ΔBn=0.

Proof: The system (25) to (28) can be deconstructed as follows:

Sn+1−Snh=[Λ−βSnIne−µτN−(µ+ρ)Sn+ωRn],

Cn+1−Cnh=[βSnIne−µτN−(γ+µ)Cn],

In+1−Inh=[γCn−(τ1+δ+μ)In],

Rn+1−Rnh=[τ1In+ρSn−(ω+µ)Rn],

After adding the above system of equations, the following can be obtained:

(Sn+1+Cn+1+In+1+Rn+1)−(Sn+Cn+In+Rn)h≤Λ−µ(Sn+Cn+In+Rn).

(Sn+1+Cn+1+In+1+Rn+1)−(Sn+Cn+In+Rn)≤hΛ−hµ(Sn+Cn+In+Rn).

(Sn+1+Cn+1+In+1+Rn+1)≤Λµ, as desired.

4.2 Local Stability Analysis of Nonstandard Computational Method

Consider the right-hand side of the Eqs. (25)–(28) functions F, G, H, and P and ΔBn=0.

F=S+h[Λ+ωR]1+h(βIe−µτ+(µ+ρ)),G=C+h[βSIe−µτ]1+h(γ+μ),H=I+h[γC]1+h(τ1+δ+μ),P=R+h[τ1I+ρS]1+h(ω+μ).

It is commonly understood that a system of the forms (25)–(28) converges to the model’s optimal state if and only if the Jacobian’s spectral radius, J,

JD=[∂F∂S∂F∂E∂F∂I∂F∂R∂G∂S∂G∂E∂G∂I∂G∂R∂H∂S∂P∂S∂H∂E∂P∂E∂H∂I∂P∂I∂H∂R∂P∂R](29)

∂F∂S=11+h(βIe−µτ+(µ+ρ)),∂F∂C=0,∂F∂I=−hβe−µτ(S+h[Λ+ωR])(1+h(βIe−µτ+(µ+ρ)))2,∂F∂R=hω1+h(βIe−µτ+(µ+ρ)),∂G∂S=0,∂G∂C=11+h(γ+μ),∂G∂I=hβSe−µτ1+h(γ+μ),∂G∂R=0,∂H∂S=0,∂H∂C=hγ1+h(τ1+δ+μ),∂H∂I=11+h(τ1+δ+μ),∂H∂R=0,∂P∂S=hρ1+h(ω+μ),∂P∂C=0,∂P∂I=hτ11+h(ω+μ),∂P∂R=11+h(ω+μ).