A Stochastic Model to Assess the Epidemiological Impact of Vaccine Booster Doses on COVID-19 and Viral Hepatitis B Co-Dynamics with 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

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.