A New Stochastic Split-Step θ-Nonstandard Finite Difference Method for the Developed SVIR Epidemic Model with Temporary Immunities and General Incidence Rates

In this paper, an SVIR epidemic model with temporary immunities and general incidence rates is constructed and analyzed. By utilizing Lyapunov functions, we prove the existence and uniqueness of the positive global solution of the constructed model, as well as the sufficient conditions of extinction and persistence of disease, are provided. Due to the difficulty of obtaining the analytical solution to our model, we construct two numerical schemes to generate an approximate solution to the model. The first one is called the split-step θ-Milstein (SSTM) method, and the second one is called the stochastic split-step θ-nonstandard finite difference (SSSNSFD) method, which is designed by merging split-step θ method with stochastic nonstandard finite difference method for the first time in this paper. Further, we prove the positivity, boundedness, and stability of the SSSTNSFD method. By employing the two mentioned methods, we support the validity of the studied theoretical results, as well, the effect of the length of immunity periods, parameters values of the incidence rates, and noise on the dynamics of the model are discussed and simulated. The increase in the size of time step size plays a vital role in revealing the method that preserves positivity, boundedness, and stability. To this end, a comparison between the proposed numerical methods is carried out graphically.


Introduction
Stochastic modeling is considered as one of the widely used strategies in the modeling of infectious diseases for the purpose of studying the dynamics of the diseases. Moreover, it is observed that stochastic models are usually more informative than deterministic models, where a deterministic model can predicts only a single result based on a given set of circumstances. In contrast, a stochastic model predicts a set of possible outcomes. In recent years, many researchers have proposed many mathematical models by using stochastic differential equations to describe the dynamics of epidemics (see, e.g., Refs. [1][2][3][4]). However, to obtain more realistic systems of population interactions, the authors included the temporal delays in such models and investigated their dynamics properties (see, e.g., Refs. [5][6][7]).
The vaccination can play an important role in controlling the diseases because it lowers the reproduction number and possibly decreases the number of infected individuals in the endemic area. As is known and confirmed that some vaccines confer a lifetime immunity against infection, while others provide only temporary immunity. So, Infection-induced or vaccination-induced immunity period is considered as one of the delay factors used in constructing the epidemic models, which were used by authors in many published papers (see, e.g., Refs. [8][9][10]). Due to the effective strategy of vaccines for controlling diseases, the authors in [11] established the stochastic SVIR epidemic model based on the corresponding deterministic model, which was constructed and analyzed in [12]. In particular, the authors in [11] proved the existence and uniqueness of the positive global solution of the following SVIR epidemic model: (1) Moreover, they provided sufficient conditions for the extinction and persistence of the disease. The studied stochastic model (1) was constructed with bilinear incidence rates and without considering the temporary immunity.
On the other hand, the duration of immunity is one of the most critical aspects of disease and vaccine that effectively affects the impact of vaccines on public health for population communities. It was observed that the duration of immunity acquired by the individual against infectious diseases ranges from almost non-existent to lifelong [20]. For instance, both of vaccine of varicella [21] and the vaccine of pertussis [22] provide only temporary immunity for the vaccinated individual against the infectious. For many infectious diseases, the immunity (whether infection-induced immunity or vaccinationinduced immunity) wanes either due to the loss of immune memory or evolution of the disease itself [23].
Motivated by the facts mentioned above, we develop in this article the model (1) by modeling the disease incidence rates via general functional responses and incorporating the temporary immunities. Regarding temporary immunities, we assume that, due to loss of immune memory, the evolution of the disease itself, or any other reasons, a fraction of recovered infected individuals may lose their infection-induced immunity and returns to the susceptible compartment. Also, a fraction of recovered vaccinated individuals may lose their vaccination-induced immunity and then moves to the vaccinees compartment to have booster or additional doses for enhancing or restoring the protection, which faded over time after the primary series vaccination was taken.
To distinguish between the booster and additional doses, the authors in [24] mentioned that Booster doses are given for those individuals who responded adequately to primary series vaccination in order to restore protection after it would have waned, while additional doses are given for those immunocompromised individuals who did not respond adequately to the primary series vaccination.
Specifically, the developed stochastic SVIR model is described by the following stochastic itô equations: where the letters S, I and R represent the densities of susceptible, infected, and recovered individuals, respectively, whilst the letter V represents the density of the individuals who have begun the vaccination process. The total size of population will be represented by N (i.e., N = S + V + I + R). The Biological meanings of all parameters in model (2) are listed in Table 1. In our model, we have two terms to denote the temporary immunities. The first term e −µτ 1 I(t − τ 1 ) indicates those individuals who survive natural death after they have become infected and then become susceptible due to the loss of infection-induced immunity for a specific time τ 1 > 0. The second term e −µτ 2 V (t − τ 2 ) indicates to those individuals who survive from natural death after they have completed their primary vaccine series and then move to the vaccinees compartment (V ) to have booster or additional doses due to the loss of vaccination-induced immunity for a specific time τ 2 > 0. However, due to the possibility that the vaccinees individuals have some partial immunity during the vaccination process, it is assumed that β 2 less than β 1 .

Parameter
Biological Meaning µ The recruitment rate and natural rate of population α The rate at which susceptible individuals are moved into the vaccination process β 1 The transmission coefficient between the two compartments S and I β 2 The disease transmission rate when the vaccinees contact with infected individuals before obtaining the immunity against the disease γ 1 The recovery rate of infected individuals The average rate for the vaccinees to obtain immunity and move into recovery compartment τ 1 The length of the immunity period of the recovered infected individuals τ 2 The length of the immunity period of recovered vaccinated individuals For λ i ≥ 0, η i ≥ 0, ∀i ∈ N 3 1 , the specific nonlinear functions F 1 , F 2 represent the incidence rates, in which they have the following forms: and satisfying the below properties: In addition, for i ∈ N 4 1 , W i are independent standard Brownian motions defined on a complete probability space with a filtration (Ω, f , P, {f t } t≥0 ) satisfying the usual conditions (i.e., it is right continuous and f 0 contains all P−null sets), and σ i are their intensities.
In looking at the structure of the model (2), we observe that the term R is absent in the first three equations. Therefore, the fourth equation can be ignored from the model without loss of generality. Thus, we here only debate the following model: where R = N − S − V − I. Obviously, the disease free equilibrium of model (3) is Discussing the effect of the length of immunity periods, parameter values of the incidence rates and noise on the dynamics of the model.
The strategy of this paper is broken down as follows: Section 2 is devoted to the stochastic analysis of the model (3) and divided into three sections: in Section 2.1, the existence and uniqueness of the positive global solution of the model (3) are investigated (see Theorem 1). According to Section 2.2, the stochastic reproduction number is defined and used to provide sufficient conditions for the extinction of disease (see Theorem 2). At the same time, sufficient conditions for the persistence of disease are provided in Section 2.3 (see Theorem 3). Section 3 is allocated to the numerical analysis of the developed model and involves three sections: In Section 3.1, we construct SSTM scheme for the model (3). In Section 3.2, we design and analyze a new stochastic method, namely, SSSTNSFD scheme for (3) is constructed and analyzed. With regard to the validity and effectiveness of the obtained results, both of the two constructed methods are employed to support those results graphically in Section 3.3. Finally, the conclusion of the study is discussed in Section 4.

Stochastic Model Analysis
This section is dedicated to show that model (3) has a positive global solution. Furthermore, by establishing a stochastic reproduction number (R s 0 ), we provide sufficient conditions for extinction and persistence of the disease. For some upcoming proofs, we need to reformulate the two functionsF 1 andF 2 as follows: Now, from (4), we deduceF Also, from (5), we concludeF

Existence and Uniqueness of Positive Global Solution
Let τ = max{τ 1 , τ 2 }, and we define R 3 be the Banach space of continuous functions mapping from the interval[−τ, 0] into R 3 + and is equipped by the norm ||Ψ|| = sup −τ≤ξ≤0 |Ψ(ξ)|. Biologically, we assume the initial conditions of model (3) to be : According to the following theorem, we can deduce that there is a unique positive global solution of model (3) with the described conditions (10). (3) admits a unique positive global solution given by (S(t), V (t), I(t)), ∀t ≥ 0, and with probability one, the solution will remain in R 3 + almost surely (a.s.).

Extinction of Disease
The proofs of Lemmas 1 and 2 are almost the same to those in Lemma 3.1 [10] and Lemma 2 [25]. Thus, herein, they can be omitted. Now, we define the reproduction number for model (3) as: , which is arguably the most significant quantity in infectious disease epidemiology that used for estimating the average number of new infections caused by an infectious individual in a population where a fraction of the susceptible individuals is vaccinated. Nonetheless, when noise does not exist, we get the reproduction number of the corresponding deterministic model as follows: For simplicity, the below notation is introduced: By utilizing the previous lemmas, sufficient conditions of extinction of disease can be obtained through the theorem below, which is one of the main results of the current article.
Proof. First, we define the function G 2 (S(t), I(t)) = S + I + γ 1 e −µτ 1 t t−τ 1 I(s)ds. Then, applying Itô's formula to get Taking the integral for (15) over [0, t], dividing by t and utilizing (6), we have From (16), we get where Second, we define another Lyapunov function as follows: Then, by utilizing from Itâ formula, we have Taking the integral for (18) over [0, t] and dividing by t, we get From (19), we obtain where t t 0 I(s)dW 3 . It follows from Lemmas 1 and 2 that lim t→∞ M 1 (t) = lim t→∞ M 2 (t) = 0 a.s. Third, we define G 4 (t) = ln I(t), then employing Itô's formula and using (6), (8) will lead us to the following: Integrating the inequality (21) over [0, t] and dividing by t, we get According to the large number theorem for martingale (see Theorem 3.4 in [26]), we have lim t→∞ The obtained result in (23) leads to lim t→∞ I(t) = 0, and then lim t→∞ I(t) = 0 a.s. Therefore, from (17) and (20), we obtain

Persistence of Disease in Mean
Definition 1. We say that model (3) is persistent in the mean, if In the forthcoming theorem, we determine the sufficient conditions for the persistence of disease for the model (3) based on Lemma 3.
if R s 0 > 1, then the solution is said to be persistent in the mean and the following properties are satisfied: Proof. Applying Itô's formula again on G 4 (t) and using (7), (9), we have Integrating inequality (25) over [0, t], and dividing by t, one can obtain From (26), we get where Since, lim t→∞ M 1 (t) = lim t→∞ M 2 (t) = lim t→∞ 2 )(R s 0 − 1) Also, from (20), we get The proof is completed.

Numerical Model Analysis
In this section, by constructing two effective methods that give dynamically consistent solutions with the continuous-time model, we intend to demonstrate the validity and effectiveness of the studied results. The first method is SSTM, and we are considering it here because it is computationally simplified. The second method will be modern, designed by incorporating the split-step θ method with a nonstandard finite difference method for the model (3) and called the SSSTNSFD method. Our modern method is constructed based on Mickens'framework, where it is free of any numerical instabilities regardless of the size of the step-size used in the numerical simulation. In what follows, h ξ i,n where ξ i,n is independent Gaussian random variable N(0, 1).

Split-Step θ-Milstein Scheme
The proposed SSTM method is easy to construct to get an approximate solution for the model (3), which was designed and used for the first time in [28]. Therefore, it is directly constructed from the model (3) as follows:

Stochastic Split-Step θ-Nonstandard Finite Difference Method
Specifically, this section aims to design and analyze a dynamically consistent SSSTNSFD method to obtain an approximate solution of the model (3), where this newly constructed method enjoys the properties of elementary stability and preservation of the positivity of solution regardless of the size of "h" used in the numerical simulations. Therefore, in order to construct our new scheme in the sense of Mickens (see [29,30]), we use non-local approximations (i.e., the terms on the right hand side of model (3) must be contain terms with the form S n+1 , S n , V n+1 , V n , I n+1 and I n ), and the terms which contain time delay are approximated by split-step θ method (i.e., the term I(t − τ 1 ) is approximated by θI n−m 1 +1 + (1 − θI n−m 1 ), and V (t − τ 2 ) is approximated by θV n−m 2 +1 + (1 − θV n−m 2 )). Moreover, for any g(t) ∈ C 1 (R), we choose an equivalent derivative which can be defined is a real-valued nonnegative function called the denominator function in which satisfies the condition v(h) = h + O(h 2 ). Therefore, for any h > 0, a common function such v(h) = 1 − e −h can be used. Consequently, the SSSTNSFD scheme is constructed as: Thence, re-arranging equations in (31) yields the following explicit SSSTNSFD scheme for model (3): It should be noted that one of the essential features of the SSSTNSFD method (32) is represented in its ease of implementation since the numerical computation of the discrete solutions of the model (3) is carried out explicitly by the following sequential process, which is mainly similar to Gauss-Seidel-method, where, we first compute S n+1 , then V n+1 , and then I n+1 . It is worth mentioning that, in the implementation of the above method (32), once a variable (e.g., S n+1 ) is computed, it is instantly used for the computations of subsequent variables (e.g., S n+1 is used for the computation of V n+1 , and then, both are used for the computation of I n+1 ). In fact, this asserts the Gauss-Seidel natural of implementing the SSSTNSFD method.
Proof. The proof of the current theorem can be obtained straightforwardly due to the fact that the constraint of biological problems is always positive.
Proof. According to Theorem 4, the numerical solution of (32) is nonnegative, and based on Theorem 5, the solution is bounded. Therefore, for all n ≥ 0, we can find positive constants The proof is completed.

Illustration and Discussion
In this section, after constructing the explicit SSTM and SSSTNSFD schemes for (3), we have two aims. The first aim is to use these schemes to show the validity of the results obtained in Section 2, as well as to discuss the effects of some parameters on the dynamics of the model (3). The second aim is to compare the two schemes in terms of dynamic properties, such as positivity, boundedness, and consistency, to show the new method's efficiency, especially when applied to larger time step.
First, if β 1 = 10.5, β 2 = 5.5, we get R s 0 = 0.9561 < 1. In addition, Figure 1a displays the results of Theorem 2 through the generated simulation by SSTM method (30), whilst Figure 1b displays the results of the same theorem through the generated simulation by SSSTNSFD method (32). It is noted that, in both simulations, the disease extincts as long as R s 0 < 1. Moreover, to show the effect of temporary immunities on dynamics of model (3), we only change the length of immunity periods as: τ 1 = 0.8, τ 2 = 1.3, and other parameter values kept unaltered. In this case, we find that R s 0 = 1.0177 > 1, and this indicates that the disease will persist if the temporary immunities are sufficiently small as simulated in Figure 2. Further, in order to show the vital role of the parameters of incidence functions on the dynamics of model (3), we only change the values of the parameters of the incidence functions as: λ 1 = 0.02, λ 2 = 0.03, λ 3 = 0.04, η 1 = 0.04, η 2 = 0.02, η 3 = 0.01, and keeping other parameters fixed, where we get R s 0 = 1.0994 > 1. It follows from this that disease persists when the parameter values of incidence functions are sufficiently small. The simulations in Figure 3 clarify that through using the two proposed methods.  Second, if β 1 = 14, β 2 = 7, then the reproduction number R s 0 = 1.2396 > 1, and the condition c 1 α+µ − β 2 αµ η 1 αµ+c 2 = 0.0977 > 0 holds. Therefore, it follows that I * = 0.125, S * = 0.0898, and V * = 0.2168. In addition, Figure 4a displays the results of Theorem 3 through the generated simulation by SSTM method (30), whilst Figure 4b displays the results of the same theorem through the generated simulation by SSSTNSFD method (32). It is noted that, in both simulations, the disease persists as long as R s 0 > 1. Moreover, in order to examine the effectiveness of noise on the dynamics of the model (3), we choose σ 3 to be sufficiently large, e.g., σ 3 = 1.5 with keeping other parameters unchanged. In this case, we get R s 0 = 0.8553 < 1 > R d 0 = 1.2402, and this reveals that the stochastic model (3) has an extinct disease, while the corresponding deterministic model has an endemic with probability one. This emphasizes that noise can suppress the disease outbreak. The simulations in Figure 5 display this fact.
Finally, in order to compare the two schemes, we increase the size of h, e.g., h = 0.5, and considering β 1 = 14, β 2 = 7, with keeping other parameters as mentioned in the set (36). As a result, we observe Figure 6a shows that the SSTM method is extremely sensitive to the time step size and fails to preserve the properties of elementary stability and preservation of positivity of the solution of the model (3) unlike the SSSTNSFD method that adheres to those properties regardless of the size of h as demonstrated in Figure 6b. Hence, this reveals the axial role that the time step size plays in detecting the most efficient numerical method.  Figure 6. Numerical simulations to compare between the two methods when the size of h is increased, e.g., h = 0.5 on [0, 10 2 ]. (a) The simulation of S, V and I paths by SSTM method. (b) The simulation of S, V and I paths by SSSTNSFD method.

Conclusions
Without a doubt, the length of the immunity period plays a significant role in the extinction or persistence of disease, where the short temporary immunity helps in the persistence of the disease, whilst the long-life immunity helps in the extinction of the disease. To support this fact, a new SVIR model has been developed and studied theoretically and numerically. In the theoretical aspect, Lyapunov functions were utilized to prove the existence and uniqueness of the global solution and to establish sufficient conditions for the extinction and persistence of disease. Additionally, the stochastic reproduction number R s 0 was established, and then we proved that if R s 0 is less than unity, the disease will die out; if R s 0 is greater than unity, it will persist in the mean. In the numerical aspect, a simulation and discussion were conducted by employing both of the two constructed methods to assess the validity of the theoretical results. However, the generated simulations in Figures 2, 3 and 5 have shown the effect of the length of immunity periods, parameters of the incidence rates, and noise on the dynamics of the model, respectively. Based on the comparison between the two used methods, we have observed through Figure 6a that the SSTM method exhibited unexpected results regarding positivity and boundedness of solutions. It was noticed that this method converges only for a small time step size, but when we increased the time step size, the method failed to restore the desired properties. On the other hand, Figure 6b showed that the SSSTNSFD method preserves all the desired properties of the model regardless of the size of the time step size used in the numerical simulation. This supports the idea that the SSSTNSFD scheme is dynamically consistent and more proper for studying the asymptotic dynamics of our model.