images
Computer Modeling in Engineering & Sciences
images

DOI: 10.32604/cmes.2023.029681

ARTICLE

A Stochastic Model to Assess the Epidemiological Impact of Vaccine Booster Doses on COVID-19 and Viral Hepatitis B Co-Dynamics with Real Data

Andrew Omame1,2,*, Mujahid Abbas3,6 and Dumitru Baleanu4,5,6

1Department of Mathematics, Federal University of Technology, Owerri, 460114, Nigeria
2Abdus Salam School of Mathematical Sciences, Government College University, Lahore, 54000, Pakistan
3Department of Mathematics, Government College University, Lahore, 54000, Pakistan
4Department of Computer Science and Mathematics, Labanese American University, Beirut, Lebanon
5Institute of Space Sciences, Magurele-Bucharest, Romania
6Department of Medical Research, China Medical University Hospital, China Medical University, Taichung, 40402, Taiwan
*Corresponding Author: Andrew Omame. Email: andrew.omame@futo.edu.ng
Received: 02 March 2023; Accepted: 04 August 2023; Published: 15 December 2023

Abstract: A patient co-infected with COVID-19 and viral hepatitis B can be at more risk of severe complications than the one infected with a single infection. This study develops a comprehensive stochastic model to assess the epidemiological impact of vaccine booster doses on the co-dynamics of viral hepatitis B and COVID-19. The model is fitted to real COVID-19 data from Pakistan. The proposed model incorporates logistic growth and saturated incidence functions. Rigorous analyses using the tools of stochastic calculus, are performed to study appropriate conditions for the existence of unique global solutions, stationary distribution in the sense of ergodicity and disease extinction. The stochastic threshold estimated from the data fitting is given by: ℛ0S=3.0651. Numerical assessments are implemented to illustrate the impact of double-dose vaccination and saturated incidence functions on the dynamics of both diseases. The effects of stochastic white noise intensities are also highlighted.

Keywords: Viral hepatitis B; COVID-19; stochastic model; extinction; ergodicity; real data

1  Introduction

According to the report issued on April 24, 2023, by Johns Hopkins University Coronavirus Center: The “Coronavirus Disease 2019” (COVID-19) caused by the “Severe Acute Respiratory Syndrome Coronavirus-2” (SARS-CoV-2) has infected 676,609,955 individuals which resulted in 6,881,955 deaths globally, and 13,338,833,198 COVID-19 vaccine doses have been administered [1]. Viral hepatitis B virus (HBV) poses a great threat to public health globally [2]. The “Centers for Disease Control and Prevention” has estimated that over 290 million individuals are infected with HBV with nearly 0.9 million deaths globally [3]. HBV is among the high-risk factors for chronic liver infections: cirrhosis, liver fibrosis and hepatocellular carcinoma. It is worth mentioning that HBV is responsible for more than 40% of hepatocellular carcinoma cases and more than 25% of liver cirrhosis cases [4]. The “Global Burden of Disease” has reported that HBV infection is one of the leading causes of adult mortality worldwide with more than 750,000 deaths annually associated with it [5]. The prevalence of HBV cases differs across different countries and regions [6]. Most HBV cases have been reported in the Middle East countries, Asia, and Africa [7]. In many of these countries, most of the reported cases are transmitted via mother-child [7]. Particularly, in Pakistan, the prevalence of viral hepatitis B is around 5 million [8]. Recent studies have indicated that between 2% and 11% of individuals infected with COVID-19 were already suffering from the liver infection and about 25% of liver problems are linked with COVID-19 [9]. SARS-CoV-2 infection can be a big risk factor for severe illness among under-diagnosed patients with viral hepatitis B [10].

According to the report issued by the “McGill COVID-19 vaccine tracker” system on December 02, 2022; 50 out of 242 developed vaccines were approved and are now available in over 200 countries [11]. Some of the approved against COVID-19 include: “BNT162b2 (Pfizer/BioNTech), AZD-1222/ChAdOx1-nCoV (Oxford/AstraZeneca), NVX-CoV2373 (Novavax), CoronaVac (Sinovac), mRNA-1273 (Moderna), rAd26-S+rAd5-S (Sputnik V), and Ad26.COV2.S (Janssen)” [12]. The effectiveness of many recommended vaccines ranges between 60% and above 90% [12]. Although antiviral agents can treat HBV, medicines alone cannot eradicate the virus from a host. Thus, vaccination has become an important resource to eliminate HBV infections [13]. That is why, different countries like China and the USA have made infant and adult vaccination programme a government priority [13,14]. This measure has significantly reduced the prevalence of HBV, as reported by Zheng et al. [14] and Zhao et al. [13]. However, developing nations in Asia and sub-Saharan Africa are still far away from achieving the target of mass vaccination, especially among the adult population.

On the other hand, mathematical modeling play a significant role not only in understanding the dynamics of the disease but also in suggesting the cost effective measures to eliminate the disease. Different epidemiological models have been proposed and analyzed to understand the dynamics of COVID-19 [15–21], Hepatitis B virus [22,23] as well as the co-infection both diseases [24–26]. Specifically, the authors [24] considered an HBV-COVID-19 model in resource limitation settings. Omame et al. [25] investigated the impact of incident co-infection in a model for HBV and COVID-19, and showed that this phenomenon could trigger backward bifurcation. Din et al. [26] investigated a co-dynamical model for HBV and COVID-19 with a bilinear incidence rate. Mathematical modelling can mainly be classified into deterministic and stochastic modelling. Though, deterministic models have their own merits and advantages but stochastic models have a potential to incorporate uncertainties and randomness and to explain the dynamics of epidemics more effectively. This is why, stochastic models have been applied successfully in understanding the patterns of diseases in recent years [27–32].

As a matter of fact, the nature of epidemic growth and spread is random due to the unpredictability in human-to-human contacts [33]. Therefore, handling the variability and randomness in the different states of the disease dynamical system is a big challenge [34] (see also, [35]). In many such cases, stochastic models could best describe the randomness of infectious contacts which takes place at different infectious stages [36]. It has been shown that stochastic models provide a higher level of realistic outputs when compared with their deterministic counterparts [37]. It was shown in [37] that an endemic equilibrium appearing in a deterministic model could disappear in its corresponding stochastic model because of stochastic fluctuations. Also, the authors in [38] explained that stochastic models give better interpretation to the question of disease extinction than their deterministic equivalents. Nasell [39] pointed out that stochastic models present better approaches in describing epidemics given large range of realistic parameter values when compared with the corresponding deterministic models.

Vaccination has become an important measure in managing the co-spread of both COVID-19 and viral hepatitis B; a feature that has not been considered by previous studies. This is the main motivation behind this work which aims to fill the gap in the existing studies by proposing a robust stochastic model for the two viral diseases co-dynamics incorporating vaccine booster doses, logistic growth and saturated incidence rates. The logistic growth deals with the population growth over time taking the carrying capacity into consideration when there is limited resources whereas, in exponential growth assumption, the population of interest grows over time and the carrying capacity is not considered which is not realistic and hence the logistic growth is most suitable [40]. Also, effective contacts between infected and uninfected humans may saturate at high infection peaks as a result of the crowding effect of infected individuals or as a result of control measures by the uninfected individuals [41]. The saturated incidence function best explains this, and has been adopted successfully in many epidemic models [42–44]. In biological models which deal with high co-endemicity of two diseases, this serves as the best incidence rate.

This study develops and analyzes the model using stochastic calculus tools. It seeks to suggest comprehensive intervention measures against the two viral diseases with the inclusion of vaccination programs for both infections. The study’s findings contribute significantly to aid in the battle against COVID-19 and viral hepatitis.

This paper contributes in the following:

(i)   A novel comprehensive stochastic model incorporating the logistic growth and saturated incidence rates is designed to assess the epidemiological effect of vaccine booster doses on the co-dynamics of viral hepatitis B and COVID-19.

(ii)   The perturbed system is fitted to the real data from Pakistan and stochastic thresholds for both diseases are estimated.

(iii)   Rigorous analysis using the tools of stochastic calculus are performed to establish appropriate conditions needed for existence of unique global solution, stationary distribution in the sense of ergodicity and disease extinction.

(iv)   The optimal levels for COVID-19 and viral hepatitis B primary and booster vaccination rates that eliminate both infections are determined.

(v)   Numerical assessments are carried out which highlight the impact of saturated incidence functions and stochastic white noise intensities on the dynamics of both diseases.

1.1 Preliminaries

We shall now recall some basic tools of stochastic calculus required in the sequel.

Theorem 1.1. [45] Let ε(t) be a “one-dimensional Ito process” on t≥0 with “stochastic differential” given by:

dε(t)=p(t,ε(t))dt+g(t,ε(t))dBt,ε(0)=ε0,(1)

where, p∈ℒ1(R+×R;R), g∈ℒ2(R+×R;R), and Bt is a “one-dimensional Brownian motion”. If ℋ∈𝒞2(R+×R;R+). Then ℋ(t,ε(t)) is another Ito process whose “stochastic differential” is given by:

dℋ(t,ε(t))=[∂ℋ(t,ε(t))∂t+∂ℋ(t,ε(t))∂εp(t,ε(t))+12∂2ℋ(t,ε(t))∂ε2g2(t,ε(t))]dt+∂ℋ(t,ε(t))∂εg(t,ε(t))dBt.(2)

The “infinitesimal generator” ℒ associated with system (1) [45] is given by:

ℒ=∂∂t+∂∂εp(t,ε(t))+12∂2∂ε2g2(t,ε(t)).(3)

Lemma 1.1. [45] If ℒ applies to a function ℋ(t,ε(t))∈C2(R+×R;R+), then

ℒℋ(t,ε(t))=∂ℋ(t,ε(t))∂t+∂ℋ(t,ε(t))∂εp(t,ε(t))+12∂2ℋ(t,ε(t))∂ε2g2(t,ε(t)).(4)

Also, the one dimensional Ito’s lemma [44] can be re-written as:

dℋ(t,ε(t))=ℒℋ(t,ε(t))dt+∂ℋ(t,ε(t))∂εg(t,ε(t))dBt.(5)

2  Model Formulation

The compartments for the formulation of the proposed model are defined as follows: S(t): susceptibles, Va(t): those vaccinated against COVID-19, Vb(t): those vaccinated against viral hepatitis B, Ia(t): infected with COVID-19, Ib(t): infected with viral hepatitis B, Icb(t): infected with the dual infections, Ra(t): those who have recored from COVID-19, Rb(t): those who have recovered from viral hepatitis B, Rab(t): those who have recovered from the dual infections, and at any time t, the total population is given by N(t)=S(t)+Va(t)+Vb(t)+Ic(t)+Ia(t)+Iab(t)+Ra(t)+Rb(t)+Rab(t). As the constant recruitment rate which has been used in many epidemic models is unrealistic, this study assumes the logistic growth in the proposed model, where r denotes per capita birth rate, K stands for the carrying capacity of the environment. The COVID-19 and viral hepatitis B transmission rates are defined as: βa and βb, respectively. Unlike several existing models dealing with COVID-19 and hepatitis B which assume the bilinear or standard incidence rates, the saturated incidence functions are used in this model, where the saturation effects associated with COVID-19 and viral hepatitis B are denoted by α1 and α2, respectively. This form of incidence has been considered most favorable when dealing with disease models involving large number of infectives. COVID-19 and viral hepatitis B primary vaccination rates are defined by ψ and ρ. The parameters: θc and θh denote the COVID-19 and viral hepatitis B booster dose vaccination rates. Immunity due to COVID-19 and viral hepatitis B vaccines are not lifelong, and wanes at the rates δa and δb, respectively. It is assumed that Individuals who recovered from COVID-19 are immune to re-infection. Similar assumption is made for individuals who have recovered from viral hepatitis B. However, infection with the other disease is possible. Due to the imperfect nature of the COVID-19 and viral hepatitis B vaccines, vulnerable or unvaccinated individuals have reduced rates of getting infected with either of the viral diseases with σ(0<σ<1) and γ(0<γ<1) denoting the COVID-19 and viral hepatitis B vaccine inefficacies. The model flow chart is given in Fig. 1, and parameters description is given in Table 1. The proposed model (both unperturbed and equivalent perturbed form) is described by the nonlinear systems defined by Eqs. (6) and (7), respectively.

Figure 1: Model’s schematic diagram

dS(t)dt=rN(1−NK)−βaIa1+α1IaS−βbIb1+α2IbS+δaVa+δbVb−(μ+ψ+ρ)SdVa(t)dt=ψS−σβaIa1+α1IaVa−βbIb1+α2IbVa−(δa+θa+μ)VadVb(t)dt=ρS−βaIa1+α1IaVb−γβbIb1+α2IbVb−(δb+θb+μ)VbdIa(t)dt=βaIa1+α1Ia(S+Rb+σVa+Vb)−(ξa+ηa+μ)Ia−φ1βbIb1+α2IbIadIb(t)dt=βbIb1+α2Ib(S+Ra+Va+γVb)−(ξb+ηb+μ)Ib−φ2βaIa1+α1IaIbdIab(t)dt=φ1βbIb1+α2IbIa+φ2βaIa1+α1IaIb−(ξab+ηab+μ)IabdRa(t)dt=θaVa+ξaIa−μRa−βbIb1+α2IbRadRb(t)dt=θbVb+ξbIb−μRb−βaIa1+α1IaRbdRab(t)dt=ξabIab−μRab(6)

The stochastic model analogue of the above deterministic model is given by:

dS(t)=(rN(1−NK)−βaIa1+α1IaS−βbIb1+α2IbS+δaVa+δbVb−(μ+ψ+ρ)S)dt+ζ1SdB1(t)dVa(t)=(ψS−σβaIa1+α1IaVa−βbIb1+α2IbVa−(δa+θa+μ)Va)dt+ζ2VadB2(t)dVb(t)=(ρS−βaIa1+α1IaVb−γβbIb1+α2IbVb−(δb+θb+μ)Vb)dt+ζ3VbdB3(t)dIa(t)=(βaIa1+α1Ia(S+Rb+σVa+Vb)−(ξa+ηa+μ)Ia−φ1βbIb1+α2IbIa)dt+ζ4IadB4(t)dIb(t)=(βbIb1+α2Ib(S+Ra+Va+γVb)−(ξb+ηb+μ)Ib−φ2βaIa1+α1IaIb)dt+ζ5IbdB5(t)dIab(t)=(φ1βbIb1+α2IbIa+φ2βaIa1+α1IaIb−(ξab+ηab+μ)Iab)dt+ζ6IabdB6(t)dRa(t)=(θaVa+ξaIa−μRa−βbIb1+α2IbRa)dt+ζ7RadB7(t)dRb(t)=(θbVb+ξbIb−μRb−βaIa1+α1IaRb)dt+ζ8RadB8(t)dRab(t)=(ξabIab−μRab)dt+ζ9RadB9(t)(7)

where, Bi(t), i=1,…,9 denote “independent Brownian motions” with Bi(0)=0, and ζi, i=1,…,9 are the intensities of white noise which reflects all random components that could impact the disease dynamics in each compartment.

3  Analysis of the Unperturbed System (6)

In this section, we now present the analysis of the unperturbed system (6).

3.1 Unperturbed System’s Reproduction Number

The disease free equilibrium (DFE) for the deterministic system is:

Q0=(S0,Va0,Vb0,Ia0,Ib0,Iab0,Ra0,Rb0,Rab0),

with,

S0=μK(r−μ)(δa+θa+μ)(δb+θb+μ)r[ψ(θa+μ)(δb+θb+μ)+ρ(θb+μ)(δa+θa+μ)+μ(δa+θa+μ)(δb+θb+μ)],Va0=ψμK(r−μ)(δb+θb+μ)r[ψ(θa+μ)(δb+θb+μ)+ρ(θb+μ)(δa+θa+μ)+μ(δa+θa+μ)(δb+θb+μ)],Vb0=ρμK(r−μ)(δa+θa+μ)r[ψ(θa+μ)(δb+θb+μ)+ρ(θb+μ)(δa+θa+μ)+μ(δa+θa+μ)(δb+θb+μ)],Ra0=θaψK(r−μ)(δb+θb+μ)r[ψ(θa+μ)(δb+θb+μ)+ρ(θb+μ)(δa+θa+μ)+μ(δa+θa+μ)(δb+θb+μ)],Rb0=θbρK(r−μ)(δa+θa+μ)r[ψ(θa+μ)(δb+θb+μ)+ρ(θb+μ)(δa+θa+μ)+μ(δa+θa+μ)(δb+θb+μ)].(8)

Associated transfer matrices are:

F=(βa(S0+σVa0+Vb0+Rb0)000βa(S0+Va0+γVb0+Ra0)0000), andV=(ξa+ηa+μ000ξb+ηb+μ000ξab+ηab+μ).(9)

The basic reproduction employing the approach in [51] is: ℛ0=ρ(FV−1)=max{ℛ0a,ℛ0b}, where ℛ0b and ℛ0a denote the deterministic reproduction numbers for viral hepatitis B and COVID-19, and are given by:

ℛ0a=βa[μ(δa+θa+μ)(δb+θb+μ)+μσψ(δa+θa+μ)+ρ(μ+θb)(δa+θa+μ)]S0μ(δa+θa+μ)(δb+θb+μ)(ξa+ηa+μ),ℛ0b=βb[μ(δa+θa+μ)(δb+θb+μ)+μγρ(δb+θb+μ)+ψ(μ+θa)(δb+θb+μ)]S0μ(δa+θa+μ)(δb+θb+μ)(ξb+ηb+μ).

3.2 Local Asymptotic Stability the Unperturbed System’s DFE (6)

Theorem 3.1. The unperturbed system’s DFE, Q0, “is locally asymptotically stable” whenever ℛ0<1, and unstable for ℛ0>1.

Proof:

Note that, the stability within neighbourhood of DFE for unperturbed system (6) is analyzed with the help of its Jacobian matrix evaluated at DFE given by:

J=(−(r−μ+ψ+ρ)Δ+δaΔ+δbΔ−βaS0Δ−βbS0Δ−(βa+βb)S0ΔΔΔ−βaS0ψ−Υ10−σβaVa0−βbVa0−(σβa+βb)Va0000ρ0−Υ20−βaVb0−γβbVa0−(βa+γβb)Vb000000βaA0−K10βaA00000000βbB0−K2βbB000000000−K30000θa0ξa00−μ0000θb0ξb00−μ000000−ξab00−μ),(10)

where,

Δ=2μ−r,A0=(S0+σVa0+Vb0+Rb0),B0=(S0+Va0+γVb0+Ra0),Υ1=(δa+θa+μ),Υ2=(δb+θb+μ),K1=ξa+ηa+μ,K2=ξb+ηb+μ,K3=ξab+ηab+μ(δa+θa+μ).

The first seven eigenvalue are given by:

Φ1=−μ(with multiplicity of three),Φ2=−(r−μ+ψ+ρ),Φ3=−(δa+θa+μ)Φ4=−(δb+θb+μ),Φ5=−(ξab+ηab+μ),

and the solutions of these equations:

(Φ+(ξa+ηa+μ)(1−ℛ0a))=0,(Φ+(ξb+ηb+μ)(1−ℛ0b))=0.(11)

As all parameters of the model are positive, it is concluded that the equations in (11) will both possess roots with negative real parts whenever ℛ0=max{ℛ0a,ℛ0b}<1. Thus, the unperturbed system’s DFE, Q0 is locally asymptotically stable if ℛ0=max{ℛ0a,ℛ0b}<1.

4  Analysis of the Perturbed System (6)

4.1 Existence of Solution

In this section, we study appropriate conditions for the existence of a unique global solution to the stochastic system with the help of well defined stochastic Lyapunov functions. We now present the following result:

Theorem 4.1. Given the initial conditions

𝒢0=(S(0),Va(0),Vb(0),Ia(0),Ib(0),Iab(0),Ra(0),Rb(0),Rab(0))∈ℳ,

the perturbed system (7) has a unique global solution

𝒢t=𝒢(t)=(S(t),Va(t),Vb(t),Ia(t),Ib(t),Iab(t),Ra(t),Rb(t),Rab(t))

for t≥0, which is invariant in ℳ with unit probability, where,

ℳ={(S(t),Va(t),Vb(t),Ia(t),Ib(t),Iab(t),Ra(t),Rb(t),Rab(t))∈R9:S(t)>0,Va(t)>0,Vb(t)>0,Ia(t)>0,Ib(t)>0,Iab(t)>0,Ra(t),Rb(t)>0,Rab(t)>0,S(t)+Va(t)+Vb(t)+Ia(t)+Ib(t)+Iab(t)+Ra(t)+Rb(t)+Rab(t)≤μK(r−μ)r}.

Proof:

For given initial states 𝒢0=(S(0),Va(0),Vb(0),Ia(0),Ib(0),Iab(0),Ra(0),Rb(0),Rab(0))∈ℳ, unique local solution 𝒢t=𝒢(t)=(S(t),Va(t),Vb(t),Ia(t),Ib(t),Iab(t),Ra(t),Rb(t),Rab(t)) exists for t∈[0,τe) with τe denoting explosion time [45]. Suppose ϕ0>0 is such that

S(0),Va(0),Vb(0),Ia(0),Ib(0),Iab(0),Ra(0),Rb(0),Rab(0) stay within [1ϕ0,ϕ0]. Then, for every ϕ≥ϕ0, we define the stopping time [29]:

τϕ=inf{t∈[0,τe):min{S(t),Va(t),Vb(t),Ia(t),Ib(t),Iab(t),Ra(t),Rb(t),Rab(t)}≤1ϕormax{S(t),Va(t),Vb(t),Ia(t),Ib(t),Iab(t),Ra(t),Rb(t),Rab(t)}≥ϕ}(12)

Note that, τϕ increases as ϕ→∞. Therefore, τ∞=ϕ→∞limτϕ. It is to be shown that τ∞=∞ a.s., so that τe=∞ and

(S(t),Va(t),Vb(t),Ia(t),Ib(t),Iab(t),Ra(t),Rb(t),Rab(t))∈ℳ a.s for all t≥0.

If τ∞<∞, then ∃ T>0 and ϑ∈(0,1) such that P{τ∞≤T}>ϑ. Thus, there is an integer ϕ1≥ϕ0 such that

P{τϕ≤T}≥ϑ∀ϕ≥ϕ1.(13)

Define a 𝒞2 stochastic Lyapunov function ℋ1:R+9→R+ as:

ℋ1(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)=(S−1−ln⁡S)+(Va−1−ln⁡Va)+(Vb−1−ln⁡Vb)+(Ia−1−ln⁡Ia)+(Ib−1−ln⁡Ib)+(Iab−1−ln⁡Iab)+(Ra−1−ln⁡Ra)+(Rb−1−ln⁡Rb)+(Rab−1−ln⁡Rab),(14)

Applying Ito’s lemma [45] to (14), we have

dℋ1=(1−1S)dS+121S2(dS)2+(1−1Va)dVa+121Va2(dVa)2+(1−1Vb)dVb+121Vb2(dVb)2+(1−1Ia)dIa+121Ia2(dIa)2+(1−1Ib)dIb+121Ib2(dIb)2+(1−1Iab)dIab+121Iab2(dIab)2+(1−1Ra)dRa+121Ra2(dRa)2+(1−1Rb)dRb+121Rb2(dRb)2+(1−1Rb)dRb+121Rb2(dRb)2=(1−1S)(rN(1−NK)−βaIa1+α1IaS−βbIb1+α2IbS+δaVa+δbVb−(μ+ψ+ρ)S)dt+(1−1S)ζ1SdB1+12ζ12dt+(1−1Va)(ψS−σβaIa1+α1IaVa−βbIb1+α2IbVa−(δa+θa+μ)Va)dt+(1−1Va)ζ2VadB2+12ζ22dt+(1−1Vb)(ρS−βaIa1+α1IaVb−γβbIb1+α2IbVb−(δb+θb+μ)Vb)dt+(1−1Vb)ζ3VbdB3+12ζ32dt+(1−1Ia)(βaIa1+α1Ia(S+Rb+σVa+Vb)−(ξa+ηa+μ)Ia−φ1βbIb1+α2IbIa)dt+(1−1Ia)ζ4IadB4+12ζ42dt+(1−1Ib)(βbIb1+α2Ib(S+Ra+Va+γVb)−(ξb+ηb+μ)Ib−φ2βaIa1+α1IaIb)dt+(1−1Ib)ζ5IbdB5+12ζ52dt+(1−1Iab)(φ1βbIb1+α2IbIa+φ2βaIa1+α1IaIb−(ξab+ηab+μ)Iab)dt+(1−1Ib)ζ6IabdB6+12ζ62dt+(1−1Ra)(θaVa+ξaIa−μRa−βbIb1+α2IbRa)dt+(1−1Ra)ζ7RadB7+12ζ72dt+(1−1Rb)(θbVb+ξbIb−μRb−βaIa1+α1IaRb)dt+(1−1Rb)ζ8RbdB8+12ζ82dt+(1−1Rab)(ξabIab−μRab)dt+(1−1Rab)ζ9RabdB9+12ζ92dt.

On simplification, we obtain that

dℋ1=[(1−1S)(rN(1−NK)−βaIa1+α1IaS−βbIb1+α2IbS+δaVa+δbVb−(μ+ψ+ρ)S)+12ζ12+(1−1Va)(ψS−σβaIa1+α1IaVa−βbIb1+α2IbVa−(δa+θa+μ)Va)+12ζ22+(1−1Vb)(ρS−βaIa1+α1IaVb−γβbIb1+α2IbVb−(δb+θb+μ)Vb)+12ζ32+(1−1Ia)(βaIa1+α1Ia(S+Rb+σVa+Vb)−(ξa+ηa+μ)Ia−φ1βbIb1+α2IbIa)+12ζ42+(1−1Ib)(βbIb1+α2Ib(S+Ra+Va+γVb)−(ξb+ηb+μ)Ib−φ2βaIa1+α1IaIb)+12ζ52+(1−1Iab)(φ1βbIb1+α2IbIa+φ2βaIa1+α1IaIb−(ξab+ηab+μ)Iab)+12ζ62+(1−1Ra)(θaVa+ξaIa−μRa−βbIb1+α2IbRa)+12ζ72+(1−1Rb)(θbVb+ξbIb−μRb−βaIa1+α1IaRb)+12ζ82+(1−1Rab)(ξabIab−μRab)+12ζ92]dt+ζ1(S−a1)dB1+ζ2(Va−a2)dB2+ζ3(Vb−a3)dB3+ζ4(Ia−1)dB4+ζ5(Ib−1)dB5+ζ6(Iab−1)dB6+ζ7(Ra−a4)dB7+ζ8(Rb−a5)dB8+ζ9(Rab−a6)dB9.

Hence,

dℋ1=ℒℋ1dt+ζ1(S−a1)dB1+ζ2(Va−a2)dB2+ζ3(Vb−a3)dB3+ζ4(Ia−1)dB4+ζ5(Ib−1)dB5+ζ6(Iab−1)dB6+ζ7(Ra−a4)dB7+ζ8(Rb−a5)dB8+ζ9(Rab−a6)dB9

where,

ℒℋ1=(1−1S)(rN(1−NK)−βaIa1+α1IaS−βbIb1+α2IbS+δaVa+δbVb−(μ+ψ+ρ)S)+12ζ12+(1−1Va)(ψS−σβaIa1+α1IaVa−βbIb1+α2IbVa−(δa+θa+μ)Va)+12ζ22+(1−1Vb)(ρS−βaIa1+α1IaVb−γβbIb1+α2IbVb−(δb+θb+μ)Vb)+12ζ32+(1−1Ia)(βaIa1+α1Ia(S+Rb+σVa+Vb)−(ξa+ηa+μ)Ia−φ1βbIb1+α2IbIa)+12ζ42+(1−1Ib)(βbIb1+α2Ib(S+Ra+Va+γVb)−(ξb+ηb+μ)Ib−φ2βaIa1+α1IaIb)+12ζ52+(1−1Iab)(φ1βbIb1+α2IbIa+φ2βaIa1+α1IaIb−(ξab+ηab+μ)Iab)+12ζ62+(1−1Ra)(θaVa+ξaIa−μRa−βbIb1+α2IbRa)+12ζ72+(1−1Rb)(θbVb+ξbIb−μRb−βaIa1+α1IaRb)+12ζ82+(1−1Rab)(ξabIab−μRab)+12ζ92.

ℒℋ1≤μK(r−μ)−μK(r−μ)S+βaIa+βbIb−δaVbS+(μ+ψ+ρ)−ψSVa+σβaIa+βbIb+(δa+θa+μ)−φSVb+βaIa+γβbIb+(δb+θb+μ)−θaVaRa−ξaIaRa+μ+βbIb−θbVbRb−ξbIbRb+μ+βaIa−ηa−ηb−ηab−βaIaIa+(ξa+ηa+μ)+φ1βbIb+φ2βaIa−βbIbIb+(ξb+ηb+μ)−φ1βbIbIaIab−φ2βaIaIbIab+(ξab+ηab+μ)−ξabIabRab+μ+ζ122+ζ222+ζ322+ζ422+ζ522+ζ622+ζ722+ζ822+ζ922.

ℒℋ1≤μK(r−μ)+βaIa+βbIb+(μ+ψ+ρ)+σβaIa+βbIb+(δa+θa+μ)+βaIa+γβbIb+(δb+θb+μ)+μ+βbIb+μ+βaIa+(ξa+ηa+μ)+φ1βbIb+φ2βaIa+(ξb+ηb+μ)+(ξab+ηab+μ)+μ+ζ122+ζ222+ζ322+ζ422+ζ522+ζ622+ζ722+ζ822+ζ922.(15)

ℒℋ1≤μK(r−μ)+(3+σ)βaK(r−μ)r+(3+γ)βbK(r−μ)r+(μ+ψ+ρ)+(δa+θa+μ)+(δb+θb+μ)+μ+μ+(ξa+ηa+μ)+(ξb+ηb+μ)+(ξab+ηab+μ)+μ+ζ122+ζ222+ζ322+ζ422+ζ522+ζ622+ζ722+ζ822+ζ922=:Ψ

Thus, we have

dℋ1=Ψdt+(ζ1(S−1)dB1+ζ2(Va−1)dB2+ζ3(Vb−1)dB3+ζ4(Ia−1)dB4+ζ5(Ib−1)dB5+ζ6(Iab−1)dB6+ζ7(Ra−1)dB7+ζ8(Rb−1)dB8+ζ9(Rab−1)dB9).(16)

If both sides of Eq. (16) are integrated from 0 to τλ∧T, then we have

∫0τϕ∧Tdℋ1(S(ν),Va(ν),Vb(ν),Ia(ν),Ib(ν),Iab(ν),Ra(ν),Rb(ν),Rab(ν))≤∫0τϕ∧TΨdν+∫0τϕ∧T(ζ1(S(ν)−1)dB1+ζ2(Va(ν)−1)dB2+ζ3(Vb(ν)−1)dB3+ζ4(Ia(ν)−1)dB4+ζ5(Ib(ν)−1)dB5+ζ6(Iab(ν)−1)dB6+ζ7(Ra(ν)−1)dB7+ζ8(Ra(ν)−1)dB8+ζ9(Rb(ν)−1)dB9).(17)

On taking expectation on both sides of the above inequality, we obtain that

Eℋ1(S(τϕ∧T),Va(τϕ∧T),Va(τϕ∧T),Ia(τϕ∧T),Ib(τϕ∧T),Iab(τϕ∧T),Ra(τϕ∧T),Rb(τϕ∧T),Rab(τϕ∧T))≤ℋ1(S(0),Va(0),Vb(0),Ia(0),Ib(0),Iab(0),Ra(0),Rb(0),Rab(0))+E∫0τϕ∧TΨdν,≤ℋ1(S(0),Va(0),Vb(0),Ia(0),Ib(0),Iab(0),Ra(0),Rb(0),Rab(0))+ΨT.

Let Ωϕ={τϕ≤T}∀ϕ≥ϕ1. Then, by (13), we have P(Ωϕ)≥ϑ. Note that, for every ω∈Ω, at the least one of S(τϕ,ω) or Va(τϕ,ω) or Vb(τϕ,ω) or Ia(τϕ,ω) or Ib(τϕ,ω) or Iab(τϕ,ω) or Ra(τϕ,ω) or Rb(τϕ,ω) or Rab(τϕ,ω) is equivalent to ϕ or 1ϕ. Thus,

ℋ1(S(τϕ,ω),Va(τϕ,ω),Vb(τϕ,ω),Ia(τϕ,ω),Ib(τϕ,ω),Iab(τϕ,ω),Ra(τϕ,ω),Rb(τϕ,ω),Rab(τϕ,ω),)

is not less than

ϕ−1−ln⁡ϕor1ϕ−1−ln⁡ϕ(which equals:1ϕ−1+ln⁡ϕ).

Therefore,

ℋ1(S(τϕ,ω),Va(τϕ,ω),Vb(τϕ,ω),Ia(τϕ,ω),Ib(τϕ,ω),Iab(τϕ,ω),Ra(τϕ,ω),Rb(τϕ,ω),Rab(τϕ,ω),)≥min{(ϕ−1−ln⁡ϕ),(1ϕ−1+ln⁡ϕ)}.

Finally, we have

ℋ1(S(0),Va(0),Vb(0),Ia(0),Ib(0),Iab(0),Ra(0),Rb(0),Rab(0))+ΨT≥E(1Ωϕℋ1(S(τϕ∧T),Va(τϕ∧T),Vb(τϕ∧T),Ia(τϕ∧T),Ib(τϕ∧T),Iab(τϕ∧T),Ra(τϕ∧T),Rb(τϕ∧T),Rab(τϕ∧T)))=E(1Ωϕ(ω)ℋ1(S(τϕ,ω),Va(τϕ,ω),Vb(τϕ,ω),Ia(τϕ,ω),Ib(τϕ,ω),Iab(τϕ,ω),Ra(τϕ,ω),Rb(τϕ,ω),Rab(τϕ,ω))),≥E(min{1Ωϕ(ω)(ϕ−1−ln⁡ϕ),(1ϕ−1+ln⁡ϕ)})=min{(ϕ−1−ln⁡ϕ),(1ϕ−1+ln⁡ϕ)}E(1Ωϕ(ω))≥ϑmin{(ϕ−1−ln⁡ϕ),(1ϕ−1+ln⁡ϕ)},(18)

where, 1Ωϕ(ω) is “indicator function” of Ωϕ(ω). If ϕ→∞, then

∞>ℋ1(S(0),Va(0),Vb(0),Ia(0),Ib(0),Iab(0),Ra(0),Rb(0),Rab(0))+ΨT=∞, which is a contradiction and hence τ∞=∞.

4.2 Extinction of the Disease

The following notation and results are stated first:

⟨𝒦(t)⟩=1t∫0t𝒦(s)ds.

Theorem 4.2 (“Strong Law of Large Numbers”). [45] Let G={Gt}t≥0 be continuous and real valued local martingale vanishing at t=0 and ⟨G,G⟩t be its quadratic variation. Then

limt→∞⟨G,G⟩t=∞, a.s. ⇒ limt→∞Gt⟨G,G⟩t=0, a.s. and,limt→∞sup⟨G,G⟩tt<0, a.s. ⇒ limt→∞Gtt=0, a.s.(19)

Lemma 4.1. Let (S(t),Va(t),Vb(t),Ia(t),Ib(t),Iab(t),Ra(t),Rb(t),Rab(t)) be a solution of perturbed system (7) subject to

(S(0),Va(0),Vb(0),Ia(0),Ib(0),Iab(0),Ra(0),Rb(0),Rab(0))∈(R)+9, then

limt→∞S(t)+Va(t)+Vb(t)+Ia(t)+Ib(t)+Iab(t)+Ra(t)+Rb(t)+Rab(t)t=0,a.s.(20)

Moreover, if μ>(ζ12∨ζ22∨ζ32∨ζ42∨ζ52∨ζ62∨ζ72∨ζ82∨ζ92)2, then

limt→∞∫0tS(ν)dB1(ν)t=0,limt→∞∫0tVa(ν)dB2(ν)t=0,limt→∞∫0tVb(ν)dB3(ν)t=0,limt→∞∫0tIa(ν)dB4(ν)t=0,limt→∞∫0tIb(ν)dB5(ν)t=0,limt→∞∫0tIab(ν)dB6(ν)t=0,limt→∞∫0tRa(ν)dB7(ν)t=0,limt→∞∫0tRb(ν)dB8(ν)t=0,limt→∞∫0tRab(ν)dB9(ν)t=0.(21)

This Proof follows using the arguments similar to those given in Lemmas (2.1) and (2.2) in [32] and is omitted.

The threshold parameter ℛ0S for the perturbed system (7) is given as:

ℛ0S=max{ℛ0aS,ℛ0bS},(22)

where,

ℛ0aS=βa[μ(δa+θa+μ)(δb+θb+μ)+μσψ(δa+θa+μ)+ρ(μ+θb)(δa+θa+μ)]S0μ(δa+θa+μ)(δb+θb+μ)(ξa+ηa+μ)−ζ422(ξa+ηa+μ)=ℛ0a−ζ422(ξa+ηa+μ),ℛ0bS=βb[μ(δa+θa+μ)(δb+θb+μ)+μγρ(δb+θb+μ)+ψ(μ+θa)(δb+θb+μ)]S0μ(δa+θa+μ)(δb+θb+μ)(ξb+ηb+μ)−ζ522(ξb+ηb+μ)=ℛ0b−ζ522(ξb+ηb+μ).

The theorem below provides necessary requirement for disease extinction from a population.

Theorem 4.3. Given the initial states (S(0),Va(0),Vb(0),Ia(0),Ib(0),Iab(0),Ra(0),Rb(0),Rab(0))∈R+9, the solution

(S(t),Va(t),Vb(t),Ia(t),Ib(t),Iab(t),Ra(t),Rb(t),Rab(t)) of the stochastic system (7) has the properties defined below: If

(a) (i)   ζ42>βa2[μ(δa+θa+μ)(δb+θb+μ)+μσψ(δa+θa+μ)+ρ(μ+θb)(δa+θa+μ)]2(S0)22(δa+θa+μ)μ2(δb+θb+μ)2, then COVID-19 goes extinct “almost surely” (a.s).

      (ii)   ζ52>βb2[μ(δb+θb+μ)(δa+θa+μ)+μργ(δa+θa+μ)+ψ(μ+θa)(δb+θb+μ)]2(S0)22μ2(δb+θb+μ)(δa+θa+μ)2, viral hepatitis B goes extinct a.s.

(b) (i)   ℛ0aS<1, then COVID-19 is eliminated from the population with unit probability.

      (ii)   ℛ0bS<1, then viral hepatitis B is eliminated from the population with unit probability.

Thus, if the requirements (a) and (b) above are satisfied, then

limt→∞⟨log⁡Ia(t)⟩t<0,andlimt→∞⟨log⁡Ib(t)⟩t<0,andlimt→∞⟨log⁡Iab(t)⟩t<0,a.s.

That is, both viral diseases will be eliminated with probability one.

More-over,

limt→∞⟨S(t)⟩=μK(r−μ)(δa+θa+μ)(δb+θb+μ)r[ψ(θa+μ)(δb+θb+μ)+ρ(θb+μ)(δa+θa+μ)+μ(δa+θa+μ)(δb+θb+μ)],limt→∞⟨Va(t)⟩=ψμK(r−μ)(δb+θb+μ)r[ψ(θa+μ)(δb+θb+μ)+ρ(θb+μ)(δa+θa+μ)+μ(δa+θa+μ)(δb+θb+μ)],limt→∞⟨Va(t)⟩=ρμK(r−μ)(δa+θa+μ)r[ψ(θa+μ)(δb+θb+μ)+ρ(θb+μ)(δa+θa+μ)+μ(δa+θa+μ)(δb+θb+μ)],limt→∞⟨Ia(t)⟩=0,limt→∞⟨Ib(t)⟩=0,limt→∞⟨Iab(t)⟩=0,limt→∞⟨Ra(t)⟩=θaψK(r−μ)(δb+θb+μ)r[ψ(θa+μ)(δb+θb+μ)+ρ(θb+μ)(δa+θa+μ)+μ(δa+θa+μ)(δb+θb+μ)],limt→∞⟨Rb(t)⟩=θbρK(r−μ)(δa+θa+μ)r[ψ(θa+μ)(δb+θb+μ)+ρ(θb+μ)(δa+θa+μ)+μ(δa+θa+μ)(δb+θb+μ)],limt→∞⟨Rab(t)⟩=0.(23)

Proof:

(a) If Ito′s lemma is applied to the 2nd equation of (7), then we have

dlogIa(t)=[βaIa1+α1Ia(S+Rb+σVa+Vb)−(ξa+ηa+μ)Ia−φ1βbIb1+α2IbIa]1Iadt−ζ422dt+ζ4dB4(t)≤[βaIa1+α1Ia(S+Rb+σVa+Vb)−(ξa+ηa+μ)Ia]1Iadt−ζ422dt+ζ4dB4(t)≤[βaIa(S+Rb+σVa+Vb)−(ξa+ηa+μ)Ia]1Iadt−ζ422dt+ζ4dB4(t)≤[βaIa(S0+Rb0+σVa0+Vb0)−(ξa+ηa+μ)Ia]1Iadt−ζ422dt+ζ4dB4(t)=[βa(S0+Rb0+σVa0+Vb0)−(ξa+ηa+μ)]dt−ζ422dt+ζ4dB4(t)=[βa[μ(δa+θa+μ)(δb+θb+μ)+μσψ(δa+θa+μ)+ρ(μ+θb)(δa+θa+μ)]S0μ(δa+θa+μ)(δb+θb+μ)−(ξa+ηa+μ)]dt−ζ422dt+ζ4dB4(t)(24)

Re-writing (24) in the form of integral gives that

logIa(t)≤∫0t(βa[μ(δa+θa+μ)(δb+θb+μ)+μσψ(δa+θa+μ)+ρ(μ+θb)(δa+θa+μ)]S0μ(δa+θa+μ)(δb+θb+μ)−ζ422)dν(ξa+ηa+μ)t+∫0tζ4dB4(ν)+logIa(0)

This can also be written in the form given by:

logIa(t)≤−ζ422∫0t(1−βa[μ(δa+θa+μ)(δb+θb+μ)+μσψ(δa+θa+μ)+ρ(μ+θb)(δa+θa+μ)]S0ζ42μ(δa+θa+μ)(δb+θb+μ))2dr+∫0t(βa2[μ(δa+θa+μ)(δb+θb+μ)+μσψ(δa+θa+μ)+ρ(μ+θb)(δa+θa+μ)]2(S0)22ζ42μ2(δa+θa+μ)2(δb+θb+μ)2)dr−(δa+θa+μ)t+∫0tζ4dB4(ν)+logIa(0)≤−((ξa+ηa+μ)t−βa2[μ(δa+θa+μ)(δb+θb+μ)+μσψ(δa+θa+μ)+ρ(μ+θb)(δa+θa+μ)]2(S0)22ζ42μ2(δa+θa+μ)2(δb+θb+μ)2)t+∫0tζ4dB4(ν)+logIa(0)

Dividing by t, we get that

logIa(t)t≤−((ξa+ηa+μ)t−βa2[μ(δa+θa+μ)(δb+θb+μ)+μσψ(δa+θa+μ)+ρ(μ+θb)(δa+θa+μ)]2(S0)22ζ42μ2(δa+θa+μ)2(δb+θb+μ)2)+1t∫0tζ4dB4(ν)+logIa(0)t,(25)

Using the “Strong Law of Large numbers”, limt→∞1t∫0tζ4dB4(ν)=0.a.s.,

ζ42>βa2[μ(δa+θa+μ)(δb+θb+μ)+μσψ(δa+θa+μ)+ρ(μ+θb)(δa+θa+μ)]2(S0)22μ2(δa+θa+μ)2(δb+θb+μ)2, and on taking “limit superior on both sides of” (25), we have

limt→∞suplogIa(t)t≤−((ξa+ηa+μ)t−βa2[μ(δa+θa+μ)(δb+θb+μ)+μσψ(δa+θa+μ)+ρ(μ+θb)(δa+θa+μ)]2(S0)22ζ42μ2(δa+θa+μ)2(δb+θb+μ)2)<0.

That is, limt→∞Ia(0)=0.

Similarly, it can be shown that

limt→∞suplogIb(t)t≤−((ξb+ηb+μ)−βb2[μ(δb+θb+μ)(δa+θa+μ)+μργ(δa+θa+μ)+ψ(μ+θa)(δb+θb+μ)]2(S0)22ζ52μ2(δa+θa+μ)2(δb+θb+μ)2)<0,

which implies that limt→∞Ib(0)=0.

(b) Also, if Eq. (24) is integrated over the interval [0,t] and divided by t, we obtain

logIa(t)−logIa(0)t=βaIa1+α1Ia(S+Rb+σVa+Vb)Ia−(ξa+ηa+μ)−φ1βbIb1+α2Ib−ζ422+ζ4t∫0tdB4(ν),≤βa[μ(δa+θa+μ)(δb+θb+μ)+μσψ(δa+θa+μ)+ρ(μ+θb)(δa+θa+μ)]S0μ(δa+θa+μ)(δb+θb+μ)−(ξa+ηa+μ)−ζ422+ζ4t∫0tdB4(ν)=(ξa+ηa+μ)(βa[μ(δa+θa+μ)(δb+θb+μ)+μσψ(δa+θa+μ)+ρ(μ+θb)(δa+θa+μ)]S0μ(δa+θa+μ)(δb+θb+μ)(ξa+ηa+μ)−ζ422(ξa+ηa+μ)−1)+ζ4t∫0tdB4(ν)=(ξa+ηa+μ)(ℛ0aS−1)+ζ4t∫0tdB4(ν).(26)

Moreover, G(t)=ζ4t∫0tdB4(ν) is “locally continuous and G(0)=0”. Employing the Lemma (4.1) on taking t→∞, we get that

limt→∞supG(t)t=0(27)

If ℛ0aS<1, then Eq. (26) becomes

limt→∞suplogIa(t)t≤(ξa+ηa+μ)(ℛ0aS−1)<0a.s.(28)

Eq. (28) implies

limt→∞Ia(t)=0a.s.(29)

Applying Ito^ Lemma to the 5th equation of system (7), we have

logIb(t)−logIb(0)t≤(ξb+ηb+μ)(βb[μ(δa+θa+μ)(δb+θb+μ)+μγρ(δb+θb+μ)+ψ(μ+θa)(δb+θb+μ)]S0μ(δa+θa+μ)(δb+θb+μ)(ξb+ηb+μ)−ζ522(ξb+ηb+μ)−1)+ζ4t∫0tdB4(ν)=(ξb+ηb+μ)(ℛ0bS−1)+ζ4t∫0tdB4(ν).(30)

Note that G(t)=ζ4t∫0tdB4(ν), is “locally continuous martingale” and G(0)=0. By applying Lemma (4.1) with t→∞, we obtain

limt→∞supG(t)t=0.(31)

If ℛ0bS<1, then Eq. (30) results in

limt→∞suplogIb(t)t≤(ξb+ηb+μ)(ℛ0bS−1)<0a.s.(32)

Eq. (32) gives

limt→∞Ib(t)=0a.s.(33)

Using the last equation of the stochastic system, we have

limt→∞⟨Rab(t)⟩=0,a.s.

From the first equations of system (7), we get that

S(t)−S(0)=∫0t(rN(1−NK)−βaIa1+α1IaS−βbIb1+α2IbS+δaVa+δbVb)dν−(μ+ψ+ρ)∫0tS(ν)dν+ζ1∫0tS(ν)B1(ν),(34)

Considering the bound for N, and recalling (29) and (33), gives

S(t)−S(0)=μK(r−μ)rt+δa∫0tVa(ν)dν+δb∫0tVb(ν)dν−(μ+ψ+ρ)∫0tS(ν)dν+ζ1∫0tS(ν)B1(ν),(35)

Dividing by t and taking limt→∞, we obtain that

limt→∞S(t)−S(0)t=μK(r−μ)r+δalimt→∞1t∫0tVa(ν)dν+δblimt→∞1t∫0tVb(ν)dν−(μ+ψ+ρ)limt→∞1t∫0tS(ν)dν+ζ1limt→∞1t∫0tS(ν)B1(ν),=μK(r−μ)r+δalimt→∞⟨Va(t)⟩+δblimt→∞⟨Vb(t)⟩−(μ+ψ+ρ)limt→∞⟨S(t)⟩+ζ1limt→∞1t∫0tS(ν)B1(ν),(36)

which can be re-written as:

(μ+ψ+ρ)limt→∞⟨S(t)⟩−δalimt→∞⟨Va(t)⟩−δblimt→∞⟨Vb(t)⟩=μK(r−μ)r−limt→∞S(t)−S(0)t+ζ1limt→∞1t∫0tS(ν)B1(ν)(37)

On taking t→∞, the above becomes

(μ+ψ+ρ)limt→∞⟨S(t)⟩−δalimt→∞⟨Va(t)⟩−δblimt→∞⟨Vb(t)⟩=K(r−μ)μr(38)

Following the similar arguments to those given above, the following expressions can also be obtained from the stochastic system (7):

ψlimt→∞⟨S(t)⟩=(δa+θa+μ)limt→∞⟨Va(t)⟩,ρlimt→∞⟨S(t)⟩=(δb+θb+μ)limt→∞⟨Vb(t)⟩,θalimt→∞⟨Va(t)=μlimt→∞⟨Ra(t)⟩,θblimt→∞⟨Vb(t)=limt→∞⟨Rb(t)⟩(39)

Solving the Eqs. (38) and (39) simultaneously, the following bounds are obtained

limt→∞⟨S(t)⟩=μK(r−μ)(δa+θa+μ)(δb+θb+μ)r[ψ(θa+μ)(δb+θb+μ)+ρ(θb+μ)(δa+θa+μ)+μ(δa+θa+μ)(δb+θb+μ)],limt→∞⟨Va(t)⟩=ψμK(r−μ)(δb+θb+μ)r[ψ(θa+μ)(δb+θb+μ)+ρ(θb+μ)(δa+θa+μ)+μ(δa+θa+μ)(δb+θb+μ)],limt→∞⟨Va(t)⟩=ρμK(r−μ)(δa+θa+μ)r[ψ(θa+μ)(δb+θb+μ)+ρ(θb+μ)(δa+θa+μ)+μ(δa+θa+μ)(δb+θb+μ)],limt→∞⟨Ra(t)⟩=θaψK(r−μ)(δb+θb+μ)r[ψ(θa+μ)(δb+θb+μ)+ρ(θb+μ)(δa+θa+μ)+μ(δa+θa+μ)(δb+θb+μ)],limt→∞⟨Rb(t)⟩=θbρK(r−μ)(δa+θa+μ)r[ψ(θa+μ)(δb+θb+μ)+ρ(θb+μ)(δa+θa+μ)+μ(δa+θa+μ)(δb+θb+μ)].(40)

4.3 Existence of Ergodic Stationary Distribution

Definition 4.1. [52] “The transition probability function P(s,ε,t,A) is said to be time-homogeneous (and the corresponding Markov process is called time-homogeneous) if the function P(s,ε,t+s,A) is independent of s, where 0≤s≤t, t∈Rq, A∈ℬ and ℬ denotes the σ-algebra of Borel sets in Rq”.

Let Y(t) represent a “time-homogeneous regular Markov” process in R+9 defined by the SDE below:

dY(t)=p(Y)dt+∑r=1khr(Y)dBr(t).(41)

The corresponding diffusion matrix is defined as:

A(ε)=(aij(ε)),aij(ε)=∑r=1khri(ε)hrj(ε).

Lemma 4.2. [52] “Suppose there exists a bounded open domain ℳ⊂R+9 having regular boundary Γ, with the following properties”:

A1: “there exists a positive number G such that ∑i,j=1qaij(ε)ρiρj≥G|ρ|2,ε∈ℳ,ρ∈Rq”.

A2: “there exists a non-negative 𝒞2-function ℋ such that ℒℋ is negative for any Rq∖ℳ”.

“Then the Markov process Y(t) have a unique ergodic stationary distribution μ()”. That is

Pε{limT⟶∞1T∫0Tp(Y)dt=∫R+9p(ε)μ(dε)}=1,∀ε∈R+9,

with p(ε) denotes an integrable function relative to measure μ.

Theorem 4.4. Define the threshold

ℛ¯0S=μΔ1Δ2βaβb{(μ+ψ+ρ+ζ122)+(δa+θa+μ+ζ222)+(δb+θb+μ+ζ322)}(ξa+ηa+μ+ζ422)(ξb+ηb+μ+ζ522),(42)

where,

Δ1=βa[μ(δa+θa+μ)(δb+θb+μ)+μσψ(δa+θa+μ)+ρ(μ+θb)(δa+θa+μ)]S0μ(δa+θa+μ)(δb+θb+μ),Δ2=βb[μ(δa+θa+μ)(δb+θb+μ)+μγρ(δb+θb+μ)+ψ(μ+θa)(δb+θb+μ)]S0μ(δa+θa+μ)(δb+θb+μ),S0=μK(r−μ)(δa+θa+μ)(δb+θb+μ)r[ψ(θa+μ)(δb+θb+μ)+ρ(θb+μ)(δa+θa+μ)+μ(δa+θa+μ)(δb+θb+μ)].

Then the perturbed system (7) has a “unique ergodic stationary distribution” μ(⋅) whenever ℛ¯0S>1.

Proof:

For (S(0),Va(0),Vb(0),Ia(0),Ib(0),Iab(0),Ra(0),Rb(0),Rab(0))∈R+9, we have unique solution

(S(t),Va(t),Vb(t),Ia(t),Ib(t),Iab(t),Ra(t),Rb(t),Rab(t))∈R+9. Moreover, the diffusion matrix of (7) is given as:

A=[ζ12S2000000000ζ22Va2000000000ζ32Vb2000000000ζ42Ia2000000000ζ52Ib2000000000ζ62Iab2000000000ζ72Ra2000000000ζ82Rb2000000000ζ92Rab2](S(t),Va(t),Vb(t),Ia(t),Ib(t),Iab(t),Ra(t),Rb(t),Rab(t))∈R9

Suppose that G=min{ζ12S2,ζ22Va2,ζ32Vb2,ζ42Ia2,ζ53Ib2,ζ62Iab2,ζ72Ra2,ζ82Rb2,ζ92Rab2}.

where, (S(t),Va(t),Vb(t),Ia(t),Ib(t),Iab(t),Ra(t),Rb(t),Rab(t))∈D¯∈R+9

Then we have

∑i,j=19aij(S(t),Va(t),Vb(t),Ia(t),Ib(t),Iab(t),Ra(t),Rb(t),Rab(t))ρi¯ρj¯=ζ12S2ρ¯12+ζ22Va2ρ¯22+ζ32Vb2ρ¯32+ζ42Ia2ρ¯42+ζ52Ib2ρ¯52+ζ62Iab2ρ¯62+ζ72ρ¯72Ra2+ζ82ρ¯82Rb2+ζ92Rab2ρ¯92≥G|ρ¯|2,where,(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈D¯,andρ¯=(ρ1¯,ρ2¯,ρ3¯,ρ4¯,ρ5¯,ρ6¯,ρ7¯,ρ8¯,ρ9¯)∈R+9.

Therefore, the requirement A1 of Lemma (4.2) is fulfilled.

Now, consider a C2− function V:R+9→R+:

Let

ℋ3=−lnS− lnVa−lnVb−ω1lnIa−ω2lnIb+(S+Va+Vb+Ia+Ib+Iab+Ra+Rb+Rab),

where ω1>0 and ω2>0 are to be determined. Applying Ito’s Lemma, we have

ℒ(S+Va+Vb+Ia+Ib+Iab+Ra+Rb+Rab)=rN(1−NK)−μ(S+Va+Vb+Ia+Ib+Iab+Ra+Rb+Rab)−ηaIa−ηbIb−ηabIab,ℒ(−lnS)=−rNS(1−NK)+βaIa1+α1Ia+βbIb1+α2Ib−δaVaS−δbVbS+(μ+ψ+ρ)+ζ122ℒ(−lnVa)=−ψSVa+σβaIa1+α1Ia+βbIb1+α2Ib+(δa+θa+μ)+ζ222ℒ(−lnVb)=−ρSVb+βaIa1+α1Ia+γβbIb1+α2Ib+(δb+θb+μ)+ζ322ℒ(−lnIa)=−βaIa1+α1Ia(S+Rb+σVa+Vb)Ia+(ξa+ηa+μ)+φ1βbIb1+α2Ib+ζ422ℒ(−lnIb)=−βbIb1+α2Ib(S+Ra+Va+γVb)Ib+(ξb+ηb+μ)+φ2βaIa1+α1Ia+ζ522ℒ(−lnIab)=−φ1βbIb1+α2IbIaIab−φ2βaIa1+α1IaIbIab+(ξab+ηab+μ)+ζ622ℒ(−lnRa)=−θaVaRa−ξaIaRa+μ+βbIb1+α2Ib+ζ722,ℒ(−lnRb)=−θbVbRb−ξbIbRb+μ+βaIa1+α1Ia+ζ822ℒ(−lnRab)=−ξabIabRab+μ+ζ922.(43)

Thus, we get that

ℒℋ3=−rNS(1−NK)+βaIa1+α1Ia+βbIb1+α2Ib−δaVaS−δbVbS+(μ+ψ+ρ+ζ122)−ψSVa+σβaIa1+α1Ia+βbIb1+α2Ib+(δa+θa+μ+ζ222)−ρSVb+βaIa1+α1Ia+γβbIb1+α2Ib+(δb+θb+μ+ζ322)+ω1[−βaIa1+α1Ia(S+Rb+σVa+Vb)Ia+(ξa+ηa+μ+ζ422)+φ1βbIb1+α2Ib]+ω2[−βbIb1+α2Ib(S+Ra+Va+γVb)Ib+(ξb+ηb+μ+ζ522)+φ2βaIa1+α1Ia]+[rN(1−NK)−μ(S+Va+Vb+Ia+Ib+Iab+Ra+Rb+Rab)−ηaIa−ηbIb−ηabIab]≤−K(r−μ)μrS−ω1βa(S+Rb+σVa+Vb)1+α1Ia−ω2βb(S+Ra+Va+γVb)1+α2Ib+ω1(ξa+ηa+μ+ζ422)+ω2(ξb+ηb+μ+ζ522)+(μ+ψ+ρ+ζ122)+(δa+θa+μ+ζ222)+(δb+θb+μ+ζ322)+βa(2+σ)Ia1+α1Ia+βb(2+γ)Ib1+α2Ib+ω1φ1βbIb1+α2Ib+ω2φ2βaIa1+α1Ia+K(r−μ)μr−K(r−μ)μr−δaVaS−δbVbS−ψSVa−ρSVb≤−μ−ω1βa[μ(δa+θa+μ)(δb+θb+μ)+μσψ(δa+θa+μ)+ρ(μ+θb)(δa+θa+μ)]S0μ(δa+θa+μ)(δb+θb+μ)−ω2βb[μ(δa+θa+μ)(δb+θb+μ)+μγρ(δb+θb+μ)+ψ(μ+θa)(δb+θb+μ)]S0μ(δa+θa+μ)(δb+θb+μ)+ω1(ξa+ηa+μ+ζ422)+ω2(ξb+ηb+μ+ζ522)+(μ+ψ+ρ+ζ122)+(δa+θa+μ+ζ222)+(δb+θb+μ+ζ322)+βa(2+σ)Ia1+α1Ia+βb(2+γ)Ib1+α2Ib+ω1φ1βbIb1+α2Ib+ω2φ2βaIa1+α1Ia+K(r−μ)μr−K(r−μ)μr−δaVaS−δbVbS−ψSVa−ρSVb−ηaIa−ηbIb−ηabIab,

Since the arithmetic mean is always greater than or equal to the geometric mean, it implies that

ℒℋ3≤−3[μΔ1Δ2ω1ω2βaβb]13+ω1(ξa+ηa+μ+ζ422)+ω2(ξb+ηb+μ+ζ522)+K(r−μ)μr+(μ+ψ+ρ+ζ122)+(δa+θa+μ+ζ222)+(δb+θb+μ+ζ322)+βa(2+σ)Ia1+α1Ia+βb(2+γ)Ib1+α2Ib+ω1φ1βbIb1+α2Ib+ω2φ2βaIa1+α1Ia−δaVaS−δbVbS−ψSVa−ρSVb−ηaIa−ηbIb−ηabIab.

where,

Δ1=βa[μ(δa+θa+μ)(δb+θb+μ)+μσψ(δa+θa+μ)+ρ(μ+θb)(δa+θa+μ)]S0μ(δa+θa+μ)(δb+θb+μ)Δ2=βb[μ(δa+θa+μ)(δb+θb+μ)+μγρ(δb+θb+μ)+ψ(μ+θa)(δb+θb+μ)]S0μ(δa+θa+μ)(δb+θb+μ).

Let

ω1(ξa+ηa+μ+ζ422)=ω2(ξb+ηb+μ+ζ522)=μΔ1Δ2βaβb(ξa+ηa+μ+ζ422)(ξb+ηb+μ+ζ522),

with

ω1=μΔ1Δ2βaβb(ξa+ηa+μ+ζ422)2(ξb+ηb+μ+ζ522),ω2=μΔ1Δ2βaβb(ξa+ηa+μ+ζ422)(ξb+ηb+μ+ζ522)2.(44)

Consequently, we have

ℒℋ3≤−[(μ+ψ+ρ+ζ122)+(δa+θa+μ+ζ222)+(δb+θb+μ+ζ322)]−ηaIa−ηbIb−ηabIab×[μΔ1Δ2βaβb{(μ+ψ+ρ+ζ122)+(δa+θa+μ+ζ222)+(δb+θb+μ+ζ322)}(ξa+ηa+μ+ζ422)(ξb+ηb+μ+ζ522)−1]+βa(2+σ)Ia1+α1Ia+βb(2+γ)Ib1+α2Ib+ω1φ1βbIb1+α2Ib+ω2φ2βaIa1+α1Ia−δaVaS−δbVbS−ψSVa−ρSVb.(45)

It implies that

ℒℋ3≤−[(μ+ψ+ρ+ζ122)+(δa+θa+μ+ζ222)+(δb+θb+μ+ζ322)][(ℛ¯0S)−1]−ηaIa−ηbIb−ηabIab+βa(2+σ)Ia1+α1Ia+βb(2+γ)Ib1+α2Ib+ω1φ1βbIb1+α2Ib+ω2φ2βaIa1+α1Ia−δaVaS−δbVbS−ψSVa−ρSVb.

Moreover, define

ℋ4=ω3[−lnS− lnVa−lnVb−ω1lnIa−ω2lnIb+(S+Va+Vb+Ia+Ib+Iab+Ra+Rb+Rab)]−ln⁡S−ln⁡Va−ln⁡Vb−ln⁡Ra−ln⁡Rb−ln⁡Rab+S+Va+Vb+Ia+Ib+Iab+Ra+Rb+Rab.

where, ω3>0, shall be determined later.

Now, consider a C2−function ℋ:R+9→R+:

ℋ(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)=ℋ4(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)−ℋ4(S(0),Ia(0),Ib(0),Iab(0),R(0),V(0)).

Applying Ito’s Lemma, we have

ℒℋ≤ω3{−[(μ+ψ+ρ+ζ122)+(δa+θa+μ+ζ222)+(δb+θb+μ+ζ322)][(ℛ¯0S)−1]−ηaIa−ηbIb−ηabIab+βa(2+σ)Ia1+α1Ia+βb(2+γ)Ib1+α2Ib+ω1φ1βbIb1+α2Ib+ω2φ2βaIa1+α1Ia−δaVaS−δbVbS−ψSVa−ρSVb}+rN(1−NK)−μ(S+Va+Vb+Ia+Ib+Iab+Ra+Rb+Rab)−ηaIa−ηbIb−ηabIab−rNS(1−NK)+βaIa1+α1Ia+βbIb1+α2Ib−δaVaS−δbVbS+(μ+ψ+ρ)+ζ122−ψSVa+σβaIa1+α1Ia+βbIb1+α2Ib+(δa+θa+μ)+ζ222−ρSVb+βaIa1+α1Ia+γβbIb1+α2Ib+(δb+θb+μ)+ζ322−θaVaRa−ξaIaRa+μ+βbIb1+α2Ib+ζ722−θbVbRb−ξbIbRb+μ+βaIa1+α1Ia+ζ822−ξabIabRab+μ+ζ922,(46)

which results in

ℒℋ≤−ω3ω4+rN(1−NK)−μ(S+Va+Vb+Ia+Ib+Iab+Ra+Rb+Rab)−ηaIa−ηbIb−ηabIab−rNS(1−NK)+βaIa1+α1Ia+βbIb1+α2Ib−δaVaS−δbVbS+(μ+ψ+ρ)+ζ122−ψSVa+σβaIa1+α1Ia+βbIb1+α2Ib+(δa+θa+μ)+ζ222−ρSVb+βaIa1+α1Ia+γβbIb1+α2Ib+(δb+θb+μ)+ζ322−θaVaRa−ξaIaRa+μ+βbIb1+α2Ib+ζ722−θbVbRb−ξbIbRb+μ+βaIa1+α1Ia+ζ822−ξabIabRab+μ+ζ922,(47)

where,

ω4=[(μ+ψ+ρ+ζ122)+(δa+θa+μ+ζ222)+(δb+θb+μ+ζ322)][(ℛ¯0S)−1]>0.

Define the domain

D={ε1<S<1ε2,ε1<Va<1ε2,ε1<Vb<1ε2,ε1<Ia<1ε2,ε1<Ib<1ε2,ε1<Iab<1ε2,ε1<Ra<1ε2,ε1<Rb<1ε2,ε1<Rab<1ε2}.

where εi>0 for i=1,2 will be estimated later. First, we subdivide the domain R+9∖D as follows:

D1={(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈R+9,0<S≤ε1},D2={(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈R+9,0<Va≤ε1,S>ε2},D3={(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈R+9,0<Vb≤ε1,Va>ε2},D4={(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈R+9,0<Ia≤ε1,Vb>ε2},D5={(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈R+9,0<Ib≤ε1,Ia>ε2},D6={(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈R+9,0<Iab≤ε2,Ib>ε1},D7={(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈R+9,0<Ra<1ε2,Iab≥ε1},D8={(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈R+9,0<Rb<1ε2,Ra≥ε1},D9={(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈R+9,0<Rab<1ε2,Rb≥ε1},D10={(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈R+9,S≥1ε2},D11={(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈R+9,Va≥1ε2},D12={(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈R+9,Vb≥1ε2},D13={(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈R+9,Ia≥1ε2},D14={(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈R+9,Ib≥1ε2},D15={(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈R+9,Iab≥1ε2},D16={(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈R+9,Ra≥1ε2},D17={(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈R+9,Rb≥1ε2},D18={(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈R+9,Rab≥1ε2}.

Next, we show that ℒℋ<0 in all the above eighteen domains which then implies that

ℒℋ(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)<0 on R+9∖D.

Case 1. Suppose (S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈D1, then using Eq. (47), we have

ℒℋ≤−ω3ω4+rN(1−NK)−μ(S+Va+Vb+Ia+Ib+Iab+Ra+Rb+Rab)−ηaIa−ηbIb−ηabIab−rNS(1−NK)+βaIa1+α1Ia+βbIb1+α2Ib+(μ+ψ+ρ)+ζ122−ψSVa+σβaIa1+α1Ia+βbIb1+α2Ib+(δa+θa+μ)+ζ222−ρSVb+βaIa1+α1Ia+γβbIb1+α2Ib+(δb+θb+μ)+ζ322−θaVaRa−ξaIaRa+μ+βbIb1+α2Ib+ζ722−θbVbRb−ξbIbRb+μ+βaIa1+α1Ia+ζ822−ξabIabRab+μ+ζ922−δaVaS−δbVbS≤−ω3ω4+rN(1−NK)−μ(S+Va+Vb+Ia+Ib+Iab+Ra+Rb+Rab)−ηaIa−ηbIb−ηabIab−rNS(1−NK)+βaIa1+α1Ia+βbIb1+α2Ib+(μ+ψ+ρ)+ζ122−ψSVa+σβaIa1+α1Ia+βbIb1+α2Ib+(δa+θa+μ)+ζ222−ρSVb+βaIa1+α1Ia+γβbIb1+α2Ib+(δb+θb+μ)+ζ322−θaVaRa−ξaIaRa+μ+βbIb1+α2Ib+ζ722−θbVbRb−ξbIbRb+μ+βaIa1+α1Ia+ζ822−ξabIabRab+μ+ζ922−δaVaϵ1−δbVbϵ1.(48)

If ϵ1>0 is chosen to be sufficiently small such that the right hand sides (48) is not greater than zero, then ℒℋ<0 for (S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈D1.

Similarly, it can also be shown that ℒℋ<0 for (S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈D2,

(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈D3, (S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈D4, (S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈D5,

(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈D6, (S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈D7, (S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈D8,

(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈D9.

Case 2. If (S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈D10, then by Eq. (47), we get

ℒℋ≤−ω3ω4+rN(1−NK)−μ(S+Va+Vb+Ia+Ib+Iab+Ra+Rb+Rab)−ηaIa−ηbIb−ηabIab−rNS(1−NK)+βaIa1+α1Ia+βbIb1+α2Ib+(μ+ψ+ρ)+ζ122−ψSVa+σβaIa1+α1Ia+βbIb1+α2Ib+(δa+θa+μ)+ζ222−ρSVb+βaIa1+α1Ia+γβbIb1+α2Ib+(δb+θb+μ)+ζ322−θaVaRa−ξaIaRa+μ+βbIb1+α2Ib+ζ722−θbVbRb−ξbIbRb+μ+βaIa1+α1Ia+ζ822−ξabIabRab+μ+ζ922−δaVaS−δbVbS≤−ω3ω4+rN(1−NK)−μ(S+Va+Vb+Ia+Ib+Iab+Ra+Rb+Rab)−ηaIa−ηbIb−ηabIab−rNS(1−NK)+βaIa1+α1Ia+βbIb1+α2Ib+(μ+ψ+ρ)+ζ122−ψSVa+σβaIa1+α1Ia+βbIb1+α2Ib+(δa+θa+μ)+ζ222−ρSVb+βaIa1+α1Ia+γβbIb1+α2Ib+(δb+θb+μ)+ζ322−θaVaRa−ξaIaRa+μ+βbIb1+α2Ib+ζ722−θbVbRb−ξbIbRb+μ+βaIa1+α1Ia+ζ822−ξabIabRab+μ+ζ922−δaε2ε1−δbε3ε1(49)

If we take, ε1=ε22, for large value of ω3>0 and smallest value of ε2>0 such that right hand sides of (49) is not greater than zero, then ℒℋ<0 for every (S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈D7.

Similarly, we can obtain that (S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈D11,(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈D12, (S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈D13, (S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈D14, (S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈D15, (S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈D16.

Hence, there exist W>0 such that

ℒℋ(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)<−W<0for every(S,Va,Vb,Ia,Ib,Iab,Ra,Rb,Rab)∈R+9∖D.

Therefore, we have

dℋ(S,Ia,Ib,Iab,R)<−Wdt+[(ω3+1)S−(ω1+1)ζ1]dB1(t)+[(ω3+1)Ia−ω2ω3ζ2]dB2(t)+[(ω3+1)Ib−ω3ζ3]dB3(t)+[(ω3+1)Iab−ω3ζ4]dB4(t)+[(ω3+1)R−ζ5]dB5(t)+[(ω3+1)R−ζ5]dB6(t).(50)

Let (S(0),Va(0),Vb(0),Ia(0),Ib(0),Iab(0),Ra(0),Rb(0),Rab(0))=ε=(ε1,ε2,ε3,ε4,ε5,ε6,ε7,ε8,ε9)∈R+9∖D, and τε denotes the time in moving from ε to D

andτn=inf{t:z=|Y(t)|},andτ(z)(t)=min{τε,t,τz}.

Taking expectation, applying Dynkin’s formula [53] and integrating both sides of Eq. (50) over [0,τ(z)(t)], gives

Eℋ(S(τ(z)(t)),Va(τ(z)(t)),Vb(τ(z)(t)),Ia(τ(z)(t)),Ib(τ(z)(t)),Iab(τ(z)(t)),Ra(τ(z)(t)),Rb(τ(z)(t)),Rab(τ(z)(t)))−ℋ(0)=E∫0τ(z)(t)ℒℋ(S(ν),Va(ν),Vb(ν),Ia(ν),Ib(ν),Iab(ν),Ra(ν),Rb(ν),Rab(ν))du≤E∫0τ(z)(t)(−W)du=−WEτ(z)(t).(51)

Since ℋ is non-negative, we have

Eτ(z)(t)≤1Wℋ(ε).

Following arguments similar to those in the proof of the Theorem (4.4), we can show that P{τe=∞}=1. Hence, the system (7) is regular. But, if z→∞ and t→∞, then τ(z)(t)→τε a.s. Hence, using Fatou’s lemma [54], it is obtained that

Eτ(z)(t)≤1Wℋ(ε)<∞.

Now, supε∈KEτε<∞, where K denotes the compact subset of R+9; condition A2 in Lemma (4.2) is satisfied. Therefore, it is concluded that system (7) has “unique stationary distribution”.

5  Numerical Scheme/Simulations

The perturbed system (7) shall now be experimented numerically in this section. The higher order scheme by Milstein [55] shall be used and is defined below:

Si+1=Si+[rNi(1−NiK)−βaSiIc,i1+α1Ic,i+βbSiIh,i1+α2Ih,i+δaVc,i++δbVh,i−(μ+ψ+ρ)Si]△t+ζ1Si△tς1,i+ζ122Si(ϖ1,i2−1)△t,Vc,i+1=Vc,i+[δaSi−σβaVc,iIc,i1+α1Ic,i−βbVc,iIh,i1+α2Ih,i−(ψ+θa+μ)Vc,i]△t+ζ2Vc,i△tϖ2,i+ζ222Vc,i(ϖ2,i2−1)△t,Vhi+1=Vh,i+[δbSi−βaVh,iIc,i1+α1Ic,i−γβbVh,iIh,i1+α2Ih,i−(ρ+θb+μ)Vh,i]△t+ζ3Vh,i△tϖ3,i+ζ322Vh,i(ϖ3,i2−1)△t,Ic,i+1=Ic,i+[βa(Si+σVc,i+Vh,i+Rh,i)Ic,i1+α1Ic,i−(ξa+ηa+μ)Ic,i−φ1Ih,iIc,i1+α2Ih,i]△t+ζ4Ic,i△tϖ4,i+ζ422Ic,i(ϖ4,i2−1)△t,Ih,i+1=Ih,i+[βb(Si+Vc,i+γVh,i+Rc,i)Ih,i1+α2Ih,i+(ξb+ηb+μ)Ih,i−φ2Ic,iIh,i1+α1Ic,i]△t+ζ5Ih,i△tϖ5,i+ζ522Ih,i(ϖ5,i2−1)△t,Ich,i+1=Ich,i+[φ1Ih,iIc,i1+α2Ih,i+φ2Ic,iIh,i1+α1Ic,i−(ξa+ηa+μ)Ich,i]△t+ζ6Ich,i△tϖ6,i+ζ622Ich,i(ϖ6,i2−1)△t,Rc,i+1=Rc,i+[θaVc,i+ξaIc,i−μRc,i−Ih,iRc,i1+α2Ih,i]△t+ζ7Rc,i△tϖ7,i+ζ722Rc,i(ϖ7,i2−1)△t,Rhi+1=Rh,i+[θbVh,i+ξbIh,i−μRh,i−Ic,iRh,i1+α1Ic,i]△t+ζ8Rh,i△tϖ8,i+ζ822Rh,i(ϖ8,i2−1)△t,Rchi+1=Rch,i+[ξabIch,i−μRch,i]△t+ζ9Rch,i△tϖ9,i+ζ922Rch,i(ϖ9,i2−1)△t,(52)

with, ϖj,i(j=1,2,3,4,5,6,7,8,9) the “independent Gaussian random variables”, N(0,1), “normal distribution” and Δt represents step size. ζj>0,(i=1,2,3,4,5,6,7,8,9) represent white noise intensities. In the sequel, numerical experiments are implemented to substantiate qualitative results established in preceding sections. Demographic data related to Pakistan have been used for the simulations. The initial conditions for the state variables are assumed thus: S(0)=175,000,000,Va(0)=15,000,000,Vb(0)=15,000,000,Ia(0)=1,296,527,Ib(0)=250,895,Iab(0)=5,000,Ra(0)=0,Rb(0)=0,Rab(0)=0. For the fitting of the model to data, available records for daily reported COVID-19 cases in Pakistan [56] between January 01, 2022 and April 10, 2022 are used. The fitting shown in Fig. 2 reveals that our perturbed model (7) behaves very well with the data. Important parameters estimated from the fitting exercise are presented in Table 1. Other parameters are estimated or derived from relevant literature. As depicted by Figs. 3–5, assessments of the perturbed system (7) are carried out, when the white noise intensities are: ζi2=0.015, for i=1,2,…,9 and when α1=2.6891×10−6, α2=7.3×10−7, βa=5.1024×10−9, βb=1.0×10−9, so that ℛ0S=max{ℛ0aS,ℛ0bS}=max{3.0323,1.5789}=3.0323>1. These assessments show that both viruses will persist in the population with unit probability. This experiment is also compatible with the conclusion of Theorem 4.4. The associated frequency distributions showing the intensity of random fluctuations, for the different compartments under this scenario are also presented alongside their solution profiles.

Figure 2: Fitting of the proposed stochastic model to real data

Figure 3: Solution profiles for the various compartments when α1=2.6891×10−6, α2=7.3×10−7, βa=5.1024×10−9, βb=1.0×10−9, so that ℛ0S=max{ℛ0aS,ℛ0bS}=max{3.0323,1.5789}=3.0323>1

Figure 4: Solution profiles for the various compartments when α1=2.6891×10−6, α2=7.3×10−7, βa=5.1024×10−9, βb=1.0×10−9, so that ℛ0S=max{ℛ0aS,ℛ0bS}=max{3.0323,1.5789}=3.0323>1

Figure 5: Solution profiles for the various compartments when α1=2.6891×10−6, α2=7.3×10−7, βa=5.1024×10−9, βb=1.0×10−9, so that ℛ0S=max{ℛ0aS,ℛ0bS}=max{3.0323,1.5789}=3.0323>1

5.1 Impact of Primary and Booster Vaccination Rates

Numerical assessments of the epidemiological impact of COVID-19 and viral hepatitis B vaccination strategies are presented in Figs. 6 and 7, respectively. The solution profiles for infected components at different primary and booster vaccination rates for COVID-19 are shown in Fig. 6. It is observed that increasing primary and booster dose vaccination rates greatly caused reduction in infected classes with COVID-19 (Fig. 6a as expected). This measure also brought about reduction in the infected individuals with viral hepatitis B and the compartment co-infected with both diseases (as can be noted in Figs. 6b and 6c). It is interesting to observe that, stepping up the COVID-19 primary vaccination rate to, ψ=0.30 and booster dose vaccination rate to the level, θc=0.20, the least number of infections is recorded for all the infected components (including those infected with viral hepatitis B). This result also support the epidemiological records from the introduction [10] showing that severe COVID-19 infection can be an important risk factor for viral hepatitis B infection. The solution profiles for the infected components at different primary and booster vaccination rates for viral hepatitis B are shown in Fig. 7. It is observed from this experiment, that increasing primary and booster dose vaccination rates for viral hepatitis B not only caused reduction in infected classes with viral hepatitis B (as seen in Fig. 7a, which is expected), but also brought about noticeable reduction in the infected individuals with COVID-19 and the co-infection of both diseases (as can be seen in Figs. 7b and 7c). It is interesting to note that, keeping the viral hepatitis primary vaccination rate at a high level, ρ=0.20 per day and the booster dose vaccination rate at θh=0.15 per day, the least number of infections is recorded for all the infected components (including singly infected with only viral hepatitis B). This result also support the epidemiological evidences [9] from the introduction section that viral hepatitis B infection could enhance susceptibility to COVID-19 infection.

Figure 6: Solution profiles for the infected compartments at different primary and booster vaccination rates for COVID-19, where βa=5.1024×10−9, βb=1.0×10−9, so that ℛ0S=max{ℛ0aS,ℛ0bS}=max{3.0323,1.5789}=3.0323>1

Figure 7: Solution profiles for the infected compartments at different primary and booster vaccination rates for viral hepatitis, where βa=5.1024×10−9, βb=1.0×10−9, so that ℛ0S=max{ℛ0aS,ℛ0bS}=max{3.0323,1.5789}=3.0323>1

5.2 Impact of Saturation Effects

The numerical investigation of the epidemiological impact of saturation effect α1 on the various components of the model is presented in Figs. 8–11. It is observed that, the saturation effect greatly impacted the different compartments. In particular, while decreasing values for α1 caused some initial increase in random fluctuations for Ic and Ich compartments (as observed in Figs. 8d and 8f), marginal random fluctuations or impact is observed for the class of individuals infected with viral hepatitis B (as can be noted in Fig. 8e).

Figure 8: Solution profiles for the various compartments at different saturation effect, α1, when βa=5.1024×10−9, βb=1.0×10−9, so that ℛ0S=max{ℛ0aS,ℛ0bS}=max{3.0323,1.5789}=3.0323>1

Figure 9: Solution profiles for the various compartments at different saturation effect, α1, when βa=5.1024×10−9, βb=1.0×10−9, so that ℛ0S=max{ℛ0aS,ℛ0bS}=max{3.0323,1.5789}=3.0323>1

Figure 10: Solution profiles for the various compartments at different saturation effect, α2, when βa=5.1024×10−9, βb=1.0×10−9, so that ℛ0S=max{ℛ0aS,ℛ0bS}=max{3.0323,1.5789}=3.0323>1

Figure 11: Solution profiles for the various compartments at different saturation effect, α2, when βa=5.1024×10−9, βb=1.0×10−9, so that ℛ0S=max{ℛ0aS,ℛ0bS}=max{3.0323,1.5789}=3.0323>1

6  Conclusion

In this paper, a comprehensive stochastic model was developed to assess the epidemiological effect of vaccine booster doses on the co-dynamics of viral hepatitis B and COVID-19 using the real data from Pakistan. The proposed model incorporates logistic growth and saturated incidence functions. Rigorous analyses employing the tools of stochastic calculus have been carried out to find appropriate conditions required for the existence of unique global solutions, stationary distribution in the sense of ergodicity and disease extinction. The stochastic threshold estimated from the data fitting is given by: ℛ0S=3.0651. Numerical assessments are implemented to illustrate the impact of double dose vaccination and saturated incidence functions on the dynamics of both diseases. The effect of stochastic white noise intensities is also highlighted. Important highlights from the numerical assessment of the perturbed model are presented as follows:

(i) The perturbed system was fitted to the real COVID-19 data from Pakistan (depicted by Fig. 2) with stochastic threshold estimated at ℛ0S=3.0651.

(ii) Increasing the COVID-19 primary vaccination rate to ψ=0.30 and booster dose vaccination rate to the level θc=0.20, the least number of infections is recorded for all the infected components (including those infected with viral hepatitis B and co-infections, as observed in Figs. 6b and 6c).

(iii) It is noted from the experiments that increasing primary and booster dose vaccination rates for viral hepatitis B not only caused reduction in infected classes with viral hepatitis B (as seen in Fig. 7a, which is expected) but also brought about a noticeable reduction in the infected individuals with COVID-19 and the co-infection of both diseases (as can be seen in Figs. 7b and 7c).

However, further investigations to improve the present study with fewer limitations can lead to some new avenues of research. More efficient algorithms can be developed for the proposed model with some more realistic assumptions. The model can also consider time dependent contact rates and time delay. Asymptomatic classes can also be considered for both viruses, and data for both diseases to be used for more accurate model fittings. In the future, we shall consider the impact of quadratic Levy noise and variable diffusion rates on the dynamics of both diseases [57–59].

Acknowledgement: The authors wish to express their appreciation to the reviewers for their helpful suggestions which greatly improved the presentation of this paper.

Funding Statement: The authors received no specific funding for this study.

Author Contributions: The authors confirm contribution to the paper as follows: study conception and design: A. O., M. A., D. B.; data collection: A. O., M. A., D. B.; analysis and interpretation of results: A. O., M. A., D. B.; draft manuscript preparation: A. O., M. A., D. B. All authors reviewed the results and approved the final version of the manuscript.

Availability of Data and Materials: All data used for the fitting of the model are available at “Pakistan: Coronavirus Pandemic Country Profile. Available online: https://ourworldindata.org/coronavirus/country/pakistan (accessed on 19/02/2023)”.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

References

1. Johns Hopkins Coronavirus Resource Center. https://coronavirus.jhu.edu/map.html (accessed on 24/04/2023)

2. Wang, S., Tao, Y., Tao, Y., Jiang, J., Yan, L. et al. (2018). Epidemiological study of hepatitis B and hepatitis C infections in Northeastern China and the beneficial effect of the vaccination strategy for hepatitis B: A cross-sectional study. BMC-Biomedical Central Public Health, 18(3), 1088. https://doi.org/10.1186/s12889-018-5984-6

3. Centers for Disease Control and Prevention, About Global Hepatitis B. https://www.cdc.gov/globalhealth/immunization/diseases/hepatitis-b/about/index.html (accessed on 25/01/2023)

4. Vos, T., Flaxman, A. D., Naghavi, M., Lozano, R., Michaud, C. et al. (2012). Years lived with disability (YLDs) for 1160 sequelae of 289 diseases and injuries 1990–2010: A systematic analysis for the global burden of disease study 2010. Lancet, 380(9859), 2163–2196.

5. Lozano, R., Naghavi, M., Foreman, K., Lim, S., Shibuya, K. et al. (2012). Global and regional mortality from 235 causes of death for 20 age groups in 1990 and 2010: A systematic analysis for the global burden of disease study 2010. Lancet, 380(9859), 2095–2128.

6. Liang, X., Bi, S., Yang, W., Wang, L., Cui, G. et al. (2009). Epidemiological serosurvey of hepatitis B in China–declining HBV prevalence due to hepatitis B vaccination. Vaccine, 27(47), 6550–6557.

7. Mokhtari, A. M., Moghadami, M., Seif, M., Mirahmadizadeh, A. (2021). Association of routine hepatitis B vaccination and other effective factors with hepatitis B virus infection: 25 Years since the introduction of national hepatitis B vaccination in Iran. Iranian Journal of Medical Sciences, 46(2), 93–102. https://doi.org/10.30476/ijms.2019.83112.1199

8. Yousafzai, M. T., Qasim, R., Khalil, R., Kakakhel, M. F., Rehman, S. U. (2014). Hepatitis B vaccination among primary health care workers in Northwest Pakistan. International Journal of Health Sciences, 8(1), 67–76. https://doi.org/10.12816/0006073

9. Li, Y., Xiao, S. Y. (2020). Hepatic involvement in COVID-19 patients: Pathology, pathogenesis, and clinical implications. Journal of Medical Virology, 92(9), 1491–1494. https://doi.org/10.1002/jmv.25973

10. Zhang, C., Shi, L., Wang, F. S. (2020). Liver injury in COVID-19: Management and challenges. The Lancet Gastroenterology & Hepatology, 5(5), 428–430. https://doi.org/10.1016/s2468-1253(20)30057-1

11. McGill COVID-19 vaccine tracker. https://covid19.trackvaccines.org/ (accessed on 25/01/2023)

12. Mohammed, I., Nauman, A., Paul, P., Ganesan, S., Chen, K. H. et al. (2022). The efficacy and effectiveness of the COVID-19 vaccines in reducing infection, severity, hospitalization, and mortality: A systematic review. Human Vaccines & Immunotherapeutics, 18(1), 2027160. https://doi.org/10.1080/21645515.2022.2027160

13. Zhao, H., Zhou, X., Zhou, Y. H. (2020). Hepatitis B vaccine development and implementation. Human Vaccines & Immunotherapeutics, 16(7), 1533–1544. https://doi.org/10.1080/21645515.2020.1732166

14. Zheng, H., Wang, F. Z., Zhang, G. M., Cui, F. Q., Wu, Z. H. et al. (2015). An economic analysis of adult hepatitis B vaccination in China. Vaccines, 33(48), 6831–6839. https://doi.org/10.1016/j.vaccine.2015.09.011

15. Din, A., Li, Y., Khan, T., Zaman, G. (2020). Mathematical analysis of spread and control of the novel corona virus (COVID-19) in China. Chaos Solitons & Fractals, 141, 110286. https://doi.org/10.1016/j.chaos.2020.110286

16. Asamoah, J. K. K., Okyere, E., Abidemi, A., Moore, S. E., Sun, G. Q. et al. (2022). Optimal control and comprehensive cost-effectiveness analysis for COVID-19. Results in Physics, 33, 105177. https://doi.org/10.1016/j.rinp.2022.105177

17. Addai, E., Zhang, L., Preko, A. K., Asamoah, J. K. K. (2022). Fractional order epidemiological model of SARS-CoV-2 dynamism involving Alzheimer’s disease. Healthcare Analytics, 2, 100114. https://doi.org/10.1016/j.health.2022.100114

18. Aslefallah, M., Yuzbasi, S., Abbasbandy, S. A. (2023). A numerical investigation based on exponential collocation method for nonlinear SITR model of COVID-19. Computer Modeling in Engineering & Sciences, 136(2), 1687–1706. https://doi.org/10.32604/cmes.2023.025647

19. Sadek, L., Sadek, O., Alaoui, H. T., Abdo, M. S., Shah, K. et al. (2023). Fractional order modeling of predicting COVID-19 with isolation and vaccination strategies in Morocco. Computer Modeling in Engineering & Science, 136(2), 1931–1950. https://doi.org/10.32604/cmes.2023.025033

20. Alnahdi, A. S., Jeelani, M. B., Wahash, H. A., Abdulwasaa, M. A. (2023). A detailed mathematical analysis of the vaccination model for COVID-19. Computer Modeling in Engineering & Science, 135(2), 1315–1343. https://doi.org/10.32604/cmes.2022.023694

21. Iqbal, S., Baleanu, D., Ali, J., Younas, H. M., Riaz, M. B. (2021). Fractional analysis of dynamical novel COVID-19 by semi-analytical technique. Computer Modeling in Engineering & Science, 129(2), 705–727. https://doi.org/10.32604/cmes.2021.015375

22. Din, A., Li, Y., Khan, F. M., Khan, Z. U., Liu, P. (2022). On analysis of fractional order mathematical model of hepatitis B using Atangana-Baleanu Caputo (ABC) derivative. Fractals, 30(1), 2240017.

23. Din, A., Li, Y., Yusuf, A., Ali, A. Y. (2022). Caputo type fractional operator applied to hepatitis B system. Fractals, 30(1), 2240023.

24. Din, A., Li, Y., Omame, A. (2022). A stochastic stability analysis of a HBV-COVID-19 co-infection model in resource limitation settings. Waves in Random and Complex Media. https://doi.org/10.1080/17455030.2022.2147598

25. Omame, A., Abbas, M. (2023). Modeling SARS-CoV-2 and HBV co-dynamics with optimal control. Physica A, 615, 128607. https://doi.org/10.1016/j.physa.2023.128607

26. Din, A., Amine, S., Allali, A. A. (2023). stochastically perturbed co-infection epidemic model for COVID-19 and hepatitis B virus. Nonlinear Dynamics, 111, 1921–1945. https://doi.org/10.1007/s11071-022-07899-1

27. Din, A., Khan, A., Baleanu, D. (2020). Stationary distribution and extinction of stochastic coronavirus (COVID-19) epidemic model. Chaos Solitons & Fractals, 139, 110036. https://doi.org/10.1016/j.chaos.2020.110036

28. Rajasekar, S. P., Pitchaimani, M. (2019). Qualitative analysis of stochastically perturbed SIRS epidemic model with two viruses. Chaos Solitons & Fractals, 207–221. https://doi.org/10.1016/j.chaos.2018.11.023

29. Okuonghae, D. (2022). Analysis of a stochastic mathematical model for tuberculosis with case detection. International Journal of Dynamics and Control, 10, 734–747. https://doi.org/10.1007/s40435-021-00863-8

30. Okuonghae, D. (2022). Ergodic stationary distribution and disease eradication in a stochastic SIR model with telegraph noises and Levy jumps. International Journal of Dynamics and Control, 10, 1778–1793. https://doi.org/10.1007/s40435-022-00962-0

31. Omame, A., Abbas, M. (2023). Din A, global asymptotic stability, extinction and ergodic stationary distribution in a stochastic model for dual variants of SARS-CoV-2. Mathematics and Computers in Simulation, 204, 302–336. https://doi.org/10.1016/j.matcom.2022.08.012

32. Zhao, Y., Jiang, D. (2014). The threshold of a stochastic SIRS epidemic model with saturated incidence. Applied Mathematics Letters, 34, 90–93.

33. Spencer, S. (2008). Stochastic epidemic models for emerging diseases (Ph.D. Thesis). University of Nottingham, UK.

34. Truscott, J. E., Gilligan, C. A. (2003). Response of a deterministic epidemiological system to a stochastically varying environment. Proceedings of the National Academy of Sciences, 100, 9067–9072.

35. Liu, M., Wang, K. (2013). Dynamics of a two-prey one predator system in random environments. Journal of Nonlinear Science, 23, 751–775.

36. Britton, T., Lindenstrand, D. (2010). Epidemic modelling: Aspects where stochastic epidemic models: A survey. Mathematical Biosciences, 222, 109–116.

37. van Herwaarden, O. A., Grasman, J. (1995). Stochastic epidemics: Major outbreaks and the duration of the endemic period. Journal of Mathematical Biology, 33, 581–601.

38. Allen, L. J. S. (2008). An introduction to stochastic epidemic models. In: Mathematical epidemiology, pp. 81–130. Berlin, Germany: Springer.

39. Nasell, I. (2002). Stochastic models of some endemic infections. Mathematical Biosciences, 179, 1–19.

40. Rajasekar, S. P., Pitchaimani, M. (2020). Ergodic stationary distribution and extinction of a stochastic SIRS epidemic model with logistic growth and nonlinear incidence. Applied Mathematics and Computation, 377, 125143. https://doi.org/10.1016/j.amc.2020.125143

41. Omame, A., Abbas, M. (2023). The stability analysis of a co-circulation model for COVID-19, dengue, and zika with nonlinear incidence rates and vaccination strategies. Healthcare Analytics, 3, 100151. https://doi.org/10.1016/j.health.2023.100151

42. Olaniyi, S. (2018). Dynamics of Zika virus model with nonlinear incidence and optimal control strategies. Applied Mathematics and Information Science, 12(5), 969–982.

43. Abidemi, A., Owolabi, K. M., Pindza, E. (2022). Modelling the transmission dynamics of Lassa fever with nonlinear incidence rate and vertical transmission. Physica A: Statistical Mechanics and Its Applications, 597, 127259. https://doi.org/10.1016/j.physa.2022.127259

44. Opara, C. Z., Uche-Iwe, N., Inyama, S. C., Omame, A. A. (2020). Mathematical model and analysis of an SVEIR model for Streptococcus pneumonia with saturated incidence force of infection. Mathematical Modelling and Applications, 5(1), 16–38. https://doi.org/10.11648/j.mma.20200501.13

45. Mao, X. (1997). Stochastic differential equations and their applications. Chichester: Horwood.

46. Omame, A., Abbas, M., Onyenegecha, C. P. (2022). A fractional order model for the co-interaction of COVID-19 and hepatitis B virus. Results in Physics, 37, 105498. https://doi.org/10.1016/j.rinp.2022.105498

47. Din, A., Li, Y., Yusuf, A. (2021). Delayed hepatitis B epidemic model with stochastic analysis. Chaos Solitons & Fractals, 146, 110839.

48. Omame, A., Okuonghae, D., Nwajeri, U. K., Onyenegecha, C. P. (2022). A fractional-order multi-vaccination model for COVID-19 with non-singular kernel. Alexandria Engineering Journal, 61(8), 6089–6104. https://doi.org/10.1016/j.aej.2021.11.037

49. Nanduri, S., Pilishvili, T., Derado, G., Soe, M. M., Dollard, P. et al. (2021). Effectiveness of Pfizer-BioNTech and Moderna vaccines in preventing SARS-CoV-2 infection among nursing home residents before and during widespread circulation of the SARS-CoV-2 B.1.617.2 (Delta) variant—National healthcare safety network, March 1-August 1, 2021. Morbidity and Mortality Weekly Report, 70(34), 1163–1166.

50. https://www.indexmundi.com/pakistan/demographics_profile.html (accessed on 25/01/2023)

51. van den Driessche, P., Watmough, J. (2002). Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical Biosciences, 180, 29–48.

52. Khas’minskii, R. (2012). Stochastic stability of differential equations, stochastic modelling and applied probability 66, 2nd edition. Berlin Heidelberg: Springer-Verlag. https://doi.org/10.1007/978-3-642-23280-0_1

53. Dynkin, E. B., Fabius, J., Greenberg, V., Maitra, A., Majone, G. (1965). Markov processes. In: Grundlehren der mathematischen Wissenschaften, vol. 121. New York: Academic Press Inc.

54. Weir, A. J. (1973). The convergence theorems. In: Lebesgue integration and measure, pp. 93–118. Cambridge: Cambridge University Press.

55. Higham, D. (2001). An algorithmic introduction to numerical simulation of stochastic differential equations. Siam Review, 43(3), 525–546.

56. Pakistan: Coronavirus Pandemic Country Profile. https://ourworldindata.org/coronavirus/country/pakistan (accessed on 11/01/2023)

57. Rajasekar, S. P., Pitchaimani, M., Zhu, Q. (2022). Probing a stochastic epidemic hepatitis C virus model with a chronically infected treated population. Acta Mathematica Scientia, 42, 2087–2112. https://doi.org/10.1007/s10473-022-0521-1

58. Sabbar, Y., Kiouach, D., Rajasekar, S. P., El-idrissi, S. E. A. (2022). The influence of quadratic Lévy noise on the dynamic of an SIC contagious illness model: New framework, critical comparison and an application to COVID-19 (SARS-CoV-2) case. Chaos, Solitons & Fractals, 159, 112110. https://doi.org/10.1016/j.chaos.2022.112110

59. Pitchaimani, M., Rajasekar, S. P. (2018). Global analysis of stochastic SIR model with variable diffusion rates. Tamkang Journal of Mathematics, 49(2), 155–182.

images
This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.