key: cord-0031056-6lp20okz
authors: Din, Anwarud; Li, Yongjin
title: Mathematical analysis of a new nonlinear stochastic hepatitis B epidemic model with vaccination effect and a case study
date: 2022-05-06
journal: Eur Phys J Plus
DOI: 10.1140/epjp/s13360-022-02748-x
sha: b576aaae793cc6a1fb0164cff364070e3571d196
doc_id: 31056
cord_uid: 6lp20okz

This work present a detailed analysis of a stochastic delayed model which governs the transmission mechanism of the Hepatitis B virus (HBV) while considering the white noises and the effect of vaccinations. It is assumed that the perturbations are nonlinear and an individual may lose his/her immunity after the vaccination, that is, the vaccination can produce temporal immunity. Based on the characteristics of the disease and the underlying assumptions, we formulated the associated deterministic model for which the threshold parameter [Formula: see text] is calculated. The model was further extended to a stochastic model and it is well-justified that the model is both mathematically and biologically feasible by showing that the model solution exists globally, bounded stochastically and is positive. By utilizing the concepts of stochastic theory and by constructing appropriate Lyapunov functions, we developed the theory for the extinction and persistence of the disease. Further, it is shown that the model is ergodic and has a unique stationary distribution. The stochastic bifurcation theory is utilized and a detailed bifurcation analysis of the model is presented. By using the standard curve fitting tools, we fitted the model against the available HBV data in Pakistan from March 2018 to February 2019 and accordingly the parameters of the model were estimated. These estimated values were used in simulating the model, theoretical findings of the study are validated through simulations and predictions were drawn. Simulations suggest that for a complete understanding of HBV dynamics, one must include time delay into such studies, and improvements in every vaccination program are unavoidable.

COVID-19 1-14 days [7] Common cold 1-3 days [8] Cholera 0.5-4.5 days [9] HIV 2-3 weeks to months or longer [10, 11] SARS 1-10 days [12] MERS 2-14 days [13] Influenza 1-3 days [14] Ebola 1-21 days [15] in view such a large number, a realistic approach is to use the standard incidence rate instead of mass action law. Further, the time delay involved in the production of viruses is 1-2 days [5] which is negligible upon comparing it with the life expectancy of a normal hepatocyte that is 6-12 months or more than that [6] . The hepatitis B virus may persist for at least 7 days outside the body. During this period, if the virus enters into a body that is not vaccinated, it might still cause illness. The hepatitis B virus has an incubation period of 30-180 days. The virus can be diagnosed 30-60 days after infection and can persist and progress to chronic hepatitis B, especially if it is transmitted during childhood or infancy. Keeping in view the incubation period of HBV, it is argued that the time delay may affect the dynamics of HBV and hence must be considered while modeling HBV. Table 1 shows the period of incubation of a few common epidemics circulating around the world.

Many fundamental concerns about infectious diseases have been answered using mathematical simulation models. From a mathematical point of view, epidemic modeling produces some interesting results, particularly in the control and prevention of contagious diseases [6, 7, 10, 11] . In the literature, several models with various traits and characteristics of infectious illnesses have been established and examined. Mathematical models play an important role in better understanding the HBV transmission rules. HBV models could be formulated at both micro and macro levels. Majority of these articles are related to studying the dynamics and control theory of HBV [16, 17] . Zou and Ruan [18] established a mathematical model (comprises of six compartments) for understanding the dynamics and prevalence of HBV in China mainland. A rigorous stability analysis and uniform persistence of the models are investigated in [19, 20] . The authors in [16] formulated an MSLIR model for describing the role of treatment and vaccination in HBV transmission. An SEICRS epidemic model with vaccination and treatment is presented to describe the spreading behavior of HBV disease in [21] . To study the dynamics of HBV in New Zealand, Mann and Roberts [22] presented an SECIR epidemic model. However, it is well-known that the stochastic elements like the humidity, temperature and precipitation, etc., have a significant impact on the force of infection, particularly in viral infections related to human. We may add randomness into deterministic models by incorporating such and similar aspects in the modeling. In biological models, these randomness will reveal the variety of environmental impacts. The varying impacts might be changes in the parameters caused by the environmental factors or random noise caused by the system [23, 24] . Furthermore, compared to deterministic models, modeling with stochastic differential equations is close to reality and provides an extra degree of freedom. The literature on stochastic models affected by the Brownian motions (or white noises) is extensive, and a detailed study can be found in the work [25] . In recent decades, linear perturbation, the most basic and obvious hypothesis, has been employed to a variety of ODE models. Population expansion can be slow down by introducing the white noise into the environment, according to above-mentioned research.

The stochastic bifurcation theory deals with the qualitative changes in the parametrized families of stochastic dynamical systems. Literature has cited the two common approaches that are used for explaining the behavior of stochastic dynamical systems: dynamical approach and phenomenological approach. The former one is related to the sign changes of Lyapunov exponents of invariant measures, whereas the later one focuses on qualitative changes of stationary densities of the corresponding stochastic processes (it is also called p-bifurcation). The latter one has no concern with the stability measured by Lyapunov exponents. There are many difficulties associated with invariant measures, for instance, at parameter values where the Lyapunov exponent is equal to zero, invariant measures bifurcate from invariant measures. Therefore, in the study of stochastic Hopf bifurcation, the phenomenological approach is far more prevalent. Invariant measures are the random analogues for deterministic fixed points, hence form statistical point of view, stationary measures are appropriate choices to describe the long-term behaviors of stochastic systems [26, 27] .

One of the Hepatitis B virus's propagation properties is the incubation time within the humans. It's the period of time between infection and infectiousness. Because the symptoms are modest, many people may not know that they are sick. Because of the lack of public understanding on HBV prevention, they may contribute to an increase in HBV infection. As a result, when studying the HBV transmission law, the incubation period's impact cannot be overlooked.

In this paper, on the basis of the above statements, we will formulate and study a hepatitis B model by adding the effects of nonlinear perturbations, time delay and vaccination. In particular, we consider the white noises, the incubation period and the efficacy of vaccination program. We will use this as a starting model for understanding the dynamics of the infection keeping in view majority of the characteristics related to HBV. The importance of vaccination and the significance of the incubation period in case of HBV are the motivations of the study. Further, the rapid spread and difficulties in treatment of HBV in district Swat motivated the authors to conduct this research. The goal of this research is to develop a health care approach that relies on vaccines and on time delay for the effective HBV control in Pakistan using the tools of mathematical modeling. Particularly, we want to look into the role of vaccine and its coverage in avoiding the disease transmission across the country. This work is also a motivation of the studies [23] [24] [25] 28] which primarily worked on the comprehension and elimination of HBV using the tools of stochastic modeling. To be more specific, this work follows the idea of [17] where the authors used deterministic approach and here, we emphasize on stochastic modeling. Majority of these works assumed the division of total population into five subgroups and hence they formulated a SVACR (susceptible-vaccinated-acute-chronic-recovered) model. The studies emphasize on the effect of vaccination on the overall dynamics of the disease.

Here, we introduce the stochastic system with variable time delay, derive its corresponding conditions for the extinction and persistence of the disease and use it to determine the quantitative and qualitative effects of the model that lead to hepatitis B virus control. Our results indicate that early vaccination that mainly modulates and restores the immune responses against the virus is mandatory for the success of the therapy. The obtained results might help government and health authorities in devising strategic vaccination policies to address the immunization gaps and ultimately to discourage the fast spread of HBV outbreak.

Practically, vaccination and the utilization of time delay are considered to be the most cost-effective strategies in order to control the transmission rate of HBV. WHO specifically suggested that the effective HBV vaccination should be administered on a case-to-case basis in high-endemicity areas. As a result, numerous writers have developed several deterministic HBV models that include the vaccination strategy [6, 23, 28, 29] . According to the characteristic of HBV infection, we shall assume the nonlinear incidence rate for the transmission of the disease. To formulate the model, we will stratify the whole population in five different compartments namely; susceptible S(t), vaccinated individuals V(t), acute infectious A(t), chronic infectious C(t), and immunized R(t). In other words, at time t, the total individuals of the population is given by N(t) = S(t) + V(t) + A(t) + C(t) + R(t). Other main assumptions taken into consideration are:

A 1 : The model's state variables and parameters are all nonnegative. A 2 : In failure case of the treatment, the infected person will move to the chronic stage as patient stays in the acute compartment for a very short period of time. A 3 : As the treatment starts, the condition of the liver is going to improve day by day and thus the disease-related death rate is not considered in the present work. A 4 : The removed population is assumed to have immunity.

Assumptions A 1 -A 4 and other related characteristics of the disease can be modeled by the following system of equations:

The parameters used in the model are defined in Table 2 . In this work, the authors assumed that the vaccination program is imperfect. That is, the vaccines will not protect the population 100%. To be more specific, the people who are vaccinated may go to the acutely infected population whenever they make a contact with the infected individuals. If τ 1 denotes the vaccination rate, then 0 ≤ τ 1 ≤ 1 and τ 1 = 0 means no vaccination while τ 1 = 1 physically shows that the whole susceptible population is vaccinated. For the deterministic model (1), we have obtained the following theoretical results with the help of [7, 13, 23] . It is quite easy to show that system (1) always admits a disease-free fixed point E 0 (it can be obtained by setting the right-hand side of each equation in model (1) to zero) as follows: The recruitment rate β

The infection rate η

The natural death rate μ

The contact rate of sick people with vulnerable individuals γ 1 The rate at which chronic population recovers γ 2 Death rate due to chronic HBV α 1 Rate at which acute patients moved to the chronic compartment α 2 The recovery rate of severely acute infected individuals κ Proportion of newborns without effective vaccination ρ

Vaccination rate

Using the next-generation matrix technique, we computed the threshold parameter R D 0 (by using the infectious classes A and C classes) from model (1) as

It can be proved that system model (1) has the following dynamical properties:

• The DFE E 0 is locally as well as globally asymptotically stable whenever R D 0 < 1, and not stable for (1) which is locally and globally asymptotically stable. In such case, the infection will tend to persist in the population for the long run.

Nevertheless, for a better explanation of HBV transmission, system (1) may be gently criticized. The reasons include the effect of environmental fluctuations, time-delay and white noises on the transmission behavior of infections, treatment and human to human interactions. Such effects forced the researchers to include the notion of stochasticity in modeling HBV epidemics. Particularly, these factors motivate us to incorporate the randomness into the models which governs the dynamics of HBV due to its high destabilizing influences. Thus, system (1) can be easily converted into a stochastic model by considering the perturbation in the form of white noises and can be expressed as

Here the notions B i (t) = {1, 2, 3, 4, 5} stand for the standard independent Brownian motion, and, respectively, the symbols δ 1 , δ 2 , δ 3 , δ 4 , δ 5 denotes the intensities of the Gaussian standard white noises. Further, the terms δ 1 S(t)dB 1 (t), δ 2 V(t)dB 2 (t), δ 3 E(t)dB 3 (t), δ 4 I(t)dB 4 (t) and δ 5 R(t)dB 5 (t) biologically shows the individuals' interaction with the environment. Model (4) is accompanied by the following subsidiary conditions

where the symbol C stand for the Banach space C [−τ, 0] : R 5 + containing the continuous functions from [−τ, 0] to R 5 + . The space R 5

. . , 5 (Fig. 1) . Since we are dealing with the population, thus, it is vital to show that the model has a nonnegative solution. It is of great concern to show that such solution is unique. In Theorem 1, the authors intend to present the existence and uniqueness of a positive solution to the model. (4) for each t ≥ 0 and an appropriate initial data

Moreover, such solution must remain in the space R 5 + surely (with probability one), in other words, almost surely (a.s) the solution (S, V, A, C, R) ∈ R 5 + for all t ≥ 0.

Proof For a suitable initial data of the states S(t), V(t), A(t), C(t) and R(t), the system coefficients are Lipschitz locally as well as continuous. This implies that there exist a unique solution in the local sense of the stochastic system in the interval t ∈ [0, τ e ).

Here, the notation τ e stand for the explosion time and the interested readers are suggested to read the works [24, 25] . It is necessary to prove that actually this solution is naturally global and for this, we have to prove that a.s. τ e = ∞. For doing so, let us consider a large enough positive number g 0 such that the interval [ 1 g 0 , g 0 ] contains all of the initial populations of each compartment. Further, we assume the explosion time for each positive integer g ≥ g 0 of the form

Where throughout this paper, we set inf ∅ = ∞ (here ∅ denotes the empty set). The definition says that τ k is increasing as k → ∞. Set τ ∞ = lim k→∞ τ k , whence τ ∞ ≤ τ e a.s. In other word, we need to show that τ ∞ = ∞ a.s. If this assertion is false, then there exist a pair of constants T > 0 and ∈ (0, 1) such that

This implies the existence of positive real number g 1 ≥ g 0 which satisfies

Next, we will define a C 2 function H :

The function defined in (7) is positive and the constant c 1 need to be calculated at later stages. Also, for all positive y, we have 0 ≤ y − log y − 1. For arbitrary g 0 ≤ G and positive T, the application of Itô's formula to Eq. (7) will yield

Eur. Phys. J. Plus (2022) 137:558

The LH in Eq. (8) is from the space R 5 + into R + which is of the form

The remaining part is similar to the proof of Theorem 2.1 in Din et al. [23] and so we must omit it.

In addition to other notions, the most important concept in epidemiology is predicting the long-term effect of an infectious disease. In this part of the work, we will provide some sufficient criteria for the extinction of the disease by using model (4) and presenting the stochastic Lyapunov functions. Model (4) always has a constant solution in which the infected compartment may vanish and such a solution if called the disease-free equilibrium. If the disease-free equilibrium of the stochastic model (4) is globally asymptotically stable, it guarantees the eradication of the infection out of the population. Thus, it is of great concern to study the extinction theory of these stochastic models.

In this section, we will show the extinction of HBV from the population considering the proposed model and will prove certain lemmas and theorems in this regard. To proceed further, let

Then Lemma 1 [23] [24] [25] (Strong Law of Large Number) Let X = {X} t≥0 be a continuous real valued local martingale and vanishing at t = 0, then lim t→∞ X, X t = ∞, a.s. ⇒ lim t→∞ X t X, X t = 0, a.s., and also

Lemma 2 For any given initial value (S 0 , V 0 , A 0 , A 0 , R 0 ) ∈ R 5 + , the solution (S(t), V(t), A(t), C(t), R(t)) of model (4) will have the following properties:

Proof The proof of Lemma 2 is similar to Lemma 4.1 in [30] , thus we omit it here. Also if the following conditions holds then we can say system (4) has a stochastic disease-free equilibrium.

where

Theorem 2 Suppose that (C 1 ) and (C 1 ) holds, then the solution of (S(t), V(t), A(t), C(t), R(t)) of system (4) with initial condition

has the following properties:

this surely shows the extinction of the disease. Proof Upon integrating system (4), one can easily obtain the following relations

By the implication of formula due to Itô, the 3rd number equation in model (4) will take the form

Integrating Eq. (18) over the interval [0, t] and a division of both sides by t leads to

Using Lemmas 1 and 2 as well as conditions (C 1 ) and (C 1 ), one can easily obtain

Again by using Lemmas 1 and 2, and conditions (C 1 ) and (C 1 ), and Eq. (20) , we have

This finishes the proof.

When examining the dynamics of an epidemic, we are concerned not only with the possibility of the disease extinction, but also with the possibility that when the disease will persist in a community. It is well understood that the stochastic model under discussion has fixed equilibria (both disease-free and endemic). As a result, in this part, we will focus on the existence of a stationary distribution, which indicates whether the infection is spreading or not, utilizing Khasminskii's theory [31] . To do so, we need to present a few lemmas and definitions which must be used in proving our main result. Assume the dynamics of a regular Markov process X(t) which is homogenous in time from the space R n + of the form

The diffusion matrix is of the form

Lemma 3 [14, 15] . The Markov process G(t) has a unique ergodic stationary distribution m(·) if a bounded domain U ∈ R d with regular bounds λ which have following two properties:

1. In the neighborhood of U as well as in its domain, the eigenvalue of matrix A(t) with the smallest is bounded away from 0. 2. If x is from the space ∈ R d \U, the average-time τ which is required for covering the path from x to U is not infinite, and sup x∈K E x τ is also finite for each K ⊂ R n (being a compact subset). Besides, for every π-measurable function f (·), we have

Due to its importance, we must define the parameter

Here, in this section, we will discuss the existence of stationary distribution and will show that model (4) Proof First of all, for proving the assertion 2 in Lemma 3, we must formulate a C 2 -function from the space R 5 + into R + . The assumed function is of the form

here the notation c i 's stand for the positive real numbers and currently unknown. By the application of the formula due to Itô and using model (4), we have the following relations

Therefore, we have

We have

Let

Here, the constant are defined as follows:

Consequently

Upon simplifications, we get the following inequality

In addition, we obtain

where the constant c 4 > 0, to be determined at later stages. It is handy to show that lim inf

where , k) . The next step is to prove that V 2 (S, V, A, C, R) has one and only one minimum value

With respect to the state variables S, V, A, C and R, the partial derivative of the function V 2 (S, V, A, C, R) can be written as It can be notice that the function V 2 has one and only one point of stagnation

Moreover, the Hessian matrix of V 2 (S, V, A, C, R) at

Obviously, the Hessian matrix is positive-definite. Thus, V 2 (S, V, A, C, R) has a minimum value V 2 (S 0 , V 0 , A 0 , C 0 , R 0 ). According to Eq. (26) and from the continuity of V 2 (S, V, A, C, R), we can say that V 2 (S, V, A, C, R) has one and only one minimum value V 2 (S 0 , V 0 , A 0 , C 0 , R 0 ) inside R 5 + . Next, we will define a nonnegative C 2 -function V : By using the stochastic model as well as the I to s formula, we have

The above relation will leads to the following assertion

where

Next, we will define a set of the form

where i for (i = 1, 2, . . . , 10) are enough small positive real numbers and currently unknown. Due to simplicity purposes, we need to make a partition of the region R 5 + \D into the following sub-parts D 1 = (S, V, A, C, R) ∈ R 5 + , 0 < S ≤ 1 , 

We can choose sufficiently small 2 > 0, yields LV < 0 for each (S, V, A, C, R) ∈ D 6 .

In the similar way we can get LV < 0, for each (S, V, A, C, R) ∈ D 7 , D 8 , D 9 , D 10 .

Concluding the above, we can say that must exist a positive constant W which gives us

Hence dV(S, V, A, C, R) < −Wdt + [(c 4 + 1)S − (c 1 c 4 + 1)

Assume that (S 0 , V 0 , A 0 , C 0 , R 0 ) = (x 1 , x 2 , x 3 , x 4 , x 5 ) = x ∈ R 5 + \D. Further, if τ x stand for the time in which the curve starting from x to D is covered, then τ n = inf{t : |X(t)| = n} and τ (n) (t) = min{τ x , t, τ n }.

By taking the integral of relation (31) over the interval [0, τ (n) (t)], applying the formula due to Dynkin and taking the expectation, we have

As V(x) is nonnegative, thus

Reference to the proof of Theorem (3), we have shown that P{τ e = ∞} = 1. Thus, we must be sure about the regularity of model (4) . Further, by letting t, n → ∞, almost surely we have τ (n)(t) → τ x . In addition, by applying the Fatou's lemma, we have

is finite and this shows that sup x∈K Eτ x is finite, here K is a subset of R 5 + and is compact. Consequently, assertion 2 in Lemma 3 is proved. Besides, model (4) has the diffusion matrix

which shows the proof of assertion 1 in Lemma 3 where ξ = (ξ 1 , ξ 2 , ξ 3 , ξ 4 , ξ 5 ) belong to the space R 5 + . Accordingly, it shows the ergodicity of the system as well as it proves that the system has a unique stationary distribution. This completes the proof of the theorem.

For the values of the parameter β from the interval (0, 1), model (1) exhibits a bifurcation, and this is due to imperfect vaccination. As the reproduction number for model (1) is given by Table 3 , the bifurcation diagram is shown in Fig. 2 . It was observed that the bifurcation diagram, for other population sizes (arbitrary N ), has almost the same shape as that of Fig. 2 . The only difference between the two shapes are the range on the vertical axis which is approximately [0, N /4]. For instance, if one assume N = 1000, the values of the range set will lies approximately within the interval [0, 250]). For β < 0.62 the value of the basic reproduction number R D 0 is 7.432 which is clearly greater than one, whereas R D 0 > 1 for β > 0. 10 . In other words, the region of bistability extends from the node bifurcation to the transcritical bifurcation when the value of β varies from 0.62 to 0.10.

Proof For any initial value (S(0), V(0), A(0), C(0), R(0) ), we let

If the conclusion of the theorem is false, then there must exist a positive constant ε 3 from the interval (0, 1) such that P ( 4 ) = ε 3 where 4 = {ω : ρ(ω) < ∞}. 

∈ V 2 (S, V, A, C, R),

This complete the proof of the theorem. From Theorems 3-4, we concluded that the stationary distribution of system (4) is confined to the region V(S, V, A, C, R). While studying the bifurcation of a nonlinear dynamical system, calculating the extreme values of the invariant measures is one of the most popular and efficient method. Particularly, the invariant measure is an important characteristic value in the stochastic bifurcation. In the present study, the shape of stationary distribution related to model (4) depends on the critical parameter β and thus, this parameter will decide the nature of bifurcation. Taking into account the parameter values from Table 3 and assuming β = 0.05, we calculated one sample path with a final time 1000 and step size 0.001. As the average time spent by a solution path in each sample set is approximately equal to the measure of this set, as a result, we have stationary measures. By looking into Fig. 3 , one can verify that the stochastic bifurcation occurs when we change the shape of stationary distribution.

Among other hypothesis related to an epidemic model, the main question is whether the model satisfies the real data? Is it support the theory? The objective of the study is achievable through the model? For the question related to data supporting, it is crucial to estimate the parameters from the data if available. Once the parameter values were identified, the model may be fitted with the real data. In this study, we assumed the real HBV data from which the values of the parameters were estimated and consequently, the stochastic, the underlying deterministic and the real data were plotted from which accuracy of the model can be seen. Basically, the real HBV cases were taken into account from the district Swat KPK Pakistan for the year 2020 and the present population of this district is nearly 2,309,570 individuals.

According to the findings of the authors in [30] , the HBV vaccines provide immunity against the infection up to twenty years. The HBV data of district Swat is taken into consideration and by using the MATLAB routine scheme lsqcurvefit, we performed the parameters estimation as shown in Table 3 .

The main theme of the lsqcurvefit algorithm is as follows. Firstly, the proposed deterministic model (1) can be comprehensively expressed as

The function F depends on time t, the vectors of dependent or state variables Y and unknown parameters θ to be estimated. The purpose of using the least square technique is to estimate the best values of model parameter which is obtain by minimizing the error between the reported data pointsȳ t i and the solution of the model y t i associated with the model parameters θ . The objective function used in the minimization procedure is given as

Eur. Phys. J. Plus (2022) 137:558 Fig. 4 Fitting of both deterministic and stochastic models with the HBV statistical data by using the values of parameter from Table 3 where n denotes the available actual data points. To obtain the model parameters, we aimed to minimize the following objective function minθ, subject to Equation (33) .

The function lsqcurvefit uses the same algorithm as lsqnonlin, however, lsqcurvefit may be preferred over other available algorithm because its main purpose is to provide an interface designed specifically for data-fitting problems. Further, these techniques can effectively handle the medium and large-scale optimization problems [32] . For more detail about this technique please see [33] and reference therein. This algorithm will be coded in MATLAB and the parameters of the model will be estimated. The predictions of stochastic system (4) as well as its underlying deterministic model were fitted against the real HBV data from January 1st, 2020, till December 28th, 2020, in Fig. 4 . By looking into the figure, we can see that the proposed stochastic model (4) fitting the data to a great degree of accuracy. In order to measure the goodness of data fit, we calculated the mean relative fit' error through the relation 1 12 12 k=1

The results are further verified by using the sufficiently small relative error formula (1.5685e − 01) which suggest that this fitness is much better compared to any other fittings.

This section of the work emphasizes on the validation of the obtain analytical results through numerical simulations. To present an approximate solution of the stochastic system, first of all, we need to discretize model (4) as follows:

Here, the notions ζ i, j for i = 1, . . . , 5) denotes the independent Gaussian stochastic variables and strictly follows the normal distribution (i. e., N (0, 1) ). Further, the notations δ i (for i = 1, . . . , 5) stand for the intensity values of the white noises and must be positive. We assumed the uniform time step size and in the scheme, it is represented by t. To simulate the stochastic model, we shall implement the algorithm presented in (35) in MATLAB and will perform the simulations. Likewise, the standard RK4 In the subsequent section, we will present graphical illustrations which will depict the long-term behavior of the model. Furthermore, we will show how the dynamics of the disease could be affected if we increase/decrease the values of a few sensitive parameters of the model.

This subsection deals with simulating model (4) by utilizing the 1st order Milstein stochastic method. As the model contains five Brownian motions dB i (t), thus for an approximation of the double stochastic integral, it is proposed to use the tools of stochastic theory.

To show the extinction of the disease from the population, we shall take the parameter values = 5.5, β = 0.2, ρ = 0.2, τ 1 = 0.02, η = 0.05 , α 1 = 0.03, α 2 = 0.2, κ = 0.2, γ 1 = 0.3 and γ 2 = 0.1. Similarly, the initial data for the population is assumed as (S 0 , V 0 , A 0 , C 0 , R 0 ) = (30, 20, 15, 10, 5) and the noise intensities are given by (δ 1 , δ 2 , δ 3 , δ 4 , δ 5 ) = (0.2, 0.1, 0.06, 0.15, 0.11). It is well known that the basic reproduction number (R 0 ) is a very important parameter in the infectious disease model, which determines whether the disease could extinct or persist in the population. Therefore, initially, we calculated that R E 0 = 0.93 which is clearly less than one. Thus, hypothesis of Theorem 2 holds and hence the conclusion must be true as one can see from Fig. 5a . This figure shows that if we keep R E 0 < 1, then surely the disease tends to extinct out of the population. By using the same set of parameters, we calculated the basic reproduction number R D 0 for the underlying deterministic model (i.e., model (1)) and again it was observed that R D 0 = 0.95 which is clearly less than unity. Based on Theorem 2, we have a similar result for the deterministic model as well. That is, if R D 0 < 1, the disease-free equilibrium is locally as well as globally asymptotically stable. Thus, it is mandatory to show that the solution curves of the deterministic model tend to the disease-free equilibrium as time evolves. This scenario is well explained through Fig. 5 where the disease tends to die out of the population if one keeps R D 0 < 1.

In order to show that the proposed stochastic model has the property of existence of a unique stationary distribution, we shall consider another set of parameters of the form = 2.5,

Here, the intensities of white noises are (δ 1 , δ 2 , δ 3 , δ 4 , δ 5 ) = (0.025, 0.1, 0.058, 0.035, 03), where the initial data in respect of each compartment were kept the same as in Sect. 8.1. By using these values of the parameters, we calculated the threshold parameter R D 0 = 1.27 for the underlying deterministic model. Clearly, the value of R D 0 is greater than one and hence, the endemic equilibrium E 0 of the said model must be globally asymptotically stable as supported by Fig. 6a-e. To be more specific, we can say that the infection will persist in the population whenever the parameter satisfies the condition of R D 0 > 1. To check the behavior of stochastic model (4) for this set of parameters, initially we calculated R S 0 = 1.7. As this parameter crossed the threshold of one, therefore, the simulations must be inline with the conclusion of Theorem 3.

In such a case, the curves fluctuating around the respective endemic equilibrium points and this is because, the intensities of white noises are too much small. This means that the disease will persist in the population in the long run. Furthermore, Theorem 3 analytically shows the existence of an ergodic stationary distribution and numerically this is confirmed with the help of Fig. 7 . Moreover, we plot the frequency distribution histogram in Fig. 7 , by using the bootstrap estimate value for each compartment. This indicates that the disease is uniform persistence. Therefore, if no further effective prevention and control measures are taken, the disease will not vanish. Fig. 6 Trajectories of susceptible, vaccinated, acute infections, chronic carriers, and recovered individuals for system (4) and (1)

This section deals with the application of the present study to explain the HBV outbreak in district Swat, Pakistan by using the available data with effect from January 2020 to December 2020. In this district, the vaccination campaign covering almost 25,000 individuals above age 1 was initiated in the month of June 2020. Based on the statistical data, we estimated the values of parameter (shown in Table 3 ) by fitting model (4) against the reported cases. Taking into consideration the nonlinear perturbations used in the model and keeping in view the data of infected HBV individuals, we will provide the effect of various parameters on the infected classes A(t) and C(t). The relative impact of these parameters on the said classes is shown in the following subsections. 9 .1 The impact of α 1 on the stochastic system (4) Assume that the value of α 1 = (7.56, 7.06) with different stochastic noises (δ 1 , δ 2 , δ 3 , δ 4 , δ 5 ) = (0.045, 0.01, 0.068, 0.025, 01) and the rest of parameter' values are taken from Table 3 . The corresponding partial solutions A(t), C(t), and the average infected relative cure of stochastic system (4) with its corresponding deterministic model are shown in Fig. 8 . Clearly, as the stochastic fluctuation of acute individuals A(t) and C(t) increases, all of the infected individuals are tends to eradicate from the community in finite time. Fig. 7 Ergodic stationary distribution of system (4) This means that the HBV in district Swat, KPK, Pakistan can be controlled and prevented in the long run if we restrict the value of α 1 . Particularly, the extinction of the disease is possible if we assume sufficiently large values of the stochastic noises and decrease the transmission coefficient. 9.2 The impact of γ 1 on the stochastic system (4) Let us choose the value of γ 1 as γ 1 = (3.55, 1.55) and the stochastic noises (δ 1 , δ 2 , δ 3 , δ 4 , δ 5 ) = (0.035, 0.02, 0.08, 0.05, 01), whereas the rest of the parameter' values are as it is in Table 3 . Figure 9 shows the corresponding partial solutions of the functions A(t), C(t) and the average infected relative cure projected by the stochastic system (4) as well as by its corresponding deterministic version. Obviously, large stochastic perturbations, will lead to the vanishing of these populations. In contrast, they have a little influence on the dynamical behavior of infected individuals. According to the analysis of the simulations, small environmental noises have a great influence on the size of A(t), C(t) and the average relative infections. Therefore, it is crucial to introduce the nonlinear stochastic noises for studying epidemic models both from practical and theoretical perspectives. More importantly, the numerical simulations such as Fig. 9a -c, reveal that some effective measures (like quarantine strategies) should be imposed in order Fig. 8 The plot shows the impact of α 1 on the dynamics of A(t) and C(t) based on stochastic systems (4) and its corresponding deterministic system to restrict environment-to-human interactions and the human activities for decreasing the spread of HBV particularly in endemic areas.

9.3 The impact of μ on the stochastic system (4) If we choose the value of μ 1 from the interval (5.11e − 02, 5.11e − 02), the stochastic noises as (δ 1 , δ 2 , δ 3 , δ 4 , δ 5 ) = (0.035, 0.02, 0.08, 0.05, 01) and assume the rest of the parameters from Table 3 . Based on these values of the parameters, Fig. 9 shows the corresponding partial solutions A(t), C(t), and the average infected relative cure of stochastic system (4) with its corresponding deterministic version, respectively. For μ 1 = (5.11 − e02, 5.11 − e01), we can see from Fig. 10 that the disease is going to die out with probability one. On the other hand, by considering the corresponding deterministic model (1), one can notice that the disease free equilibrium E 0 is globally asymptotically stable. Further, it is worthy to notice that the disease goes for extinction and this fact is supported by Fig. 10a -c. Based on the simulation results, we notice that small environmental noises greatly affect the population of A(t), C(t) and the average relative infections. Here again, from Fig. 9a -c, we can say that in modeling diseases like HBV outbreak of Swat, one must include nonlinear stochastic noises into the models. 9.4 The impact of τ 1 on the stochastic system (4) Assume the values of the parameters as in Table 3 except the stochastic noises (δ 1 , δ 2 , δ 3 , δ 4 , δ 5 ) = (0.005, 0.002, 0.006, 0.005, 001) and the value of τ 1 is taken from the interval (0.77, 0.57). Both the stochastic and its deterministic counterpart were simulated by using these parameter' values and we plot solution curves of A(t), C(t). Figure 11a -c shows that τ 1 has an inverse relationship with these population, that is, by increasing the value of τ 1 the respective infected populations may be decreased remarkably. 9.5 The impact of β on the stochastic system (4) To show the effect of β on the dynamics of the infected classes, we assumed β 1 = (14. 25, 13.25) and the stochastic noises (δ 1 , δ 2 , δ 3 , δ 4 , δ 5 ) = (0.005, 0.002, 0.006, 0.005, 001) along with other parameters as shown in Table 3 . The behavior of A(t), C(t) and the average relative infections are shown graphically in Fig. 9 by using the proposed stochastic model and its underlying deterministic counterpart. Fig. 12a -c shows that if we decrease the value of β, the corresponding populations will converge fast to the disease extinction. Thus, for disease extinction, we have to decrease the values of β. Further, these simulations suggest that we Fig. 9 Trajectories of A(t), C(t) and the average infected relative cure predicted by the stochastic systems and its corresponding deterministic counterpart Fig. 10 Trajectories of A(t), C(t) and the average relative infections by using the proposed stochastic and deterministic models Eur. Phys. J. Plus (2022) 137:558 Fig. 11 The plot shows the relationship of parameter τ 1 with the infected HBV classes based on stochastic systems (4) and its corresponding deterministic model must incorporate the nonlinear stochastic noises for better comprehension of HBV dynamics like the outbreak happened in district Swat, Pakistan.

In this research article, we studied the dynamics of a novel stochastic HBV model keeping in view the effect of vaccination and timedelay. With the help of standard existence and uniqueness theory of stochastic differential equations, we proved that the model has a unique positive bounded solution. We obtained the threshold quantity R E 0 which was used for deriving the necessary conditions of extinction and persistence of the disease. Analytically, we proved that the disease will tend to extinct out of the population whenever R E 0 < 1. However, for R S 0 > 1, we showed that the proposed model has a unique stationary distribution as well as the property of ergodicity. Biologically, the properties of persistence and stationary distribution show the existence of the infection in the community in the long run. By using the bifurcation theory of stochastic differential equations, we showed that the model exhibit a bifurcation. Further, it is shown that the region of bistability extends from the node bifurcation to the transcritical bifurcation when the value of β varies from 0.62 to 0.10. We used the standard numerical algorithm for showing the analytical results through simulations. This research work was tested against the real data obtained from district Swat (Pakistan) between January 2020 and December 2020. Both stochastic and deterministic models were fitted against the data by using the nonlinear least squares fit routine in MATLAB. The graphical illustrations show a good agreement between the data and the model's curves. Further, simulations suggest that large values of the white noise can lead to the extinction of HBV, whereas the smaller values ensure the persistence of HBV. By simulating the model several times, the authors conclude that several effective measures like the quarantine strategy and decreasing the human interactions particularly in HBV endemic areas could be helpful in the elimination of HBV.

We fitted both the deterministic and stochastic models against the available data of HBV of district Swat (Pakistan) and the values of parameters are estimated by using the lsqcurvefit techniques. These values of the parameter are utilized and the reproduction number is calculated which is approximately 13. This shows that the HBV is very dangerous and could adversely affect this locality if proper control strategies were not implemented on time. It also indicates that if there are no effective measures, the number of HBV cases may rise in the next couple of years. To decrease the spreading rate further, the health officials and policy makers must start awareness campaigns, mass vaccinations particularly in children, and importantly the treatment and care of chronic HBV patients. It was also noticed that β is the most sensitive parameter to the reproduction number and hence, reducing the transmission coefficient Fig. 12 The plot shows the behavior of A(t), C(t) and the average relative infected populations by using both the stochastic and the associated deterministic models (by reducing the sexual and household contacts of susceptible with chronic) for acute infected individuals who become chronic is also an effective measure to decrease the further spread of HBV.

In future, our aim is to conduct studies on the fractional-order versions of this model and to solve such models by using the techniques of Weyl, Riesz, Riemann-Liouville, Caputo, Caputo-Fabrizio, Atangana-Baleanu, Atangana-Gomez, fractal-fractional and other operators which have the capability to capture the complex and anomalous behavior of related dynamical systems.

Plus (2022) 137:558 References 1. WHO

Virus dynamics: mathematical principles of immunology and virology: mathematical principles of immunology and virology

Hepatitis B Vaccines

The role of antigen-independent persistence of memory cytotoxic T lymphocytes

Viral dynamics in hepatitis B virus infection

Delayed hepatitis B epidemic model with stochastic analysis

Global stability of a class of virus dynamics models with general incidence rate and multitarget cells

Incubation periods of acute respiratory viral infections: a systematic review

The incubation period of cholera: a systematic review

Addressing Ebola-related stigma: lessons learned from HIV

Acute human immunodeficiency virus type 1 infection

Advisors of expert SARS group of hospital authority. Effectiveness of precautions against droplets and contact in prevention of nosocomial transmission of severe acute respiratory syndrome (SARS)

Epidemiology and predictors of survival of MERS-CoV infections in Riyadh region

A new fractional model for the dynamics of the hepatitis B virus using the Caputo-Fabrizio derivative

Ebola virus disease outbreak in Nigeria: transmission dynamics and rapid control

Mathematical model and analysis of hepatitis B virus transmission dynamics

Mathematical analysis and simulation of a hepatitis B model with time delay: a case study for Xinjiang

On the sexual transmission dynamics of hepatitis B virus in China

Dynamic analysis of a hepatitis B model with three-age-classes

A mathematical model for hepatitis B with infection-age structure

Stability analysis of the transmission dynamics of an HBV model

Modelling the epidemiology of hepatitis B in New Zealand

Stationary distribution extinction and optimal control for the stochastic hepatitis B epidemic model with partial immunity

Stability of SIRS system with random perturbations

Multigroup SIR epidemic model with stochastic perturbation

Stationary distribution and stochastic Hopf bifurcation for a predator-prey system with noises

Van den Driessche, Stochastic epidemic models with a backward bifurcation

Impact of information intervention on stochastic hepatitis B model and its variable-order fractional network

Two decades of universal hepatitis B vaccination in Taiwan: impact and implication for future strategies

Extinction and stationary distribution of a stochastic SIRS epidemic model with standard incidence rate and partial immunity

Stochastic stability of differential equations

An algorithm for least squares estimation of nonlinear parameters

Parameter estimation in nonlinear dynamical systems

Acknowledgements This work was supported by the National Natural Science Foundation of P. R. China (Nos. 11971493 and 12071491).

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

The authors declare that they have no competing interests.