An unconditional boundary and dynamics preserving scheme for the stochastic epidemic model

Over the past few decades, mathematical models have been constructed to describe the evolution of different kinds of epidemics. The classical Kermack-McKendrick model [1] was used for modelling common childhood diseases, where a typical individual starts off susceptible, at some stage catches the disease, and after a short infectious period becomes permanently immune. Nevertheless, in terms of many unstable sources of infection, particularly most viruses, cured people can rarely get long-time immunity and turn to be susceptible to some mutations of the same original virus. For these cases, the susceptible-infected-susceptible (SIS) model was constructed in [2]. By adding randomness to the model, the classical susceptible-infected-susceptible epidemic model was extended by Gray et al. in [3] from a deterministic framework to a stochastic one, which is formulated as the following stochastic differential equations (SDEs) of Itô type:

where $S_{t}$ denotes the number of susceptibles and $I_{t}$ the number of infecteds at time $t\in[0,T]$, with initial values satisfying $S_{0}+I_{0}=N$. In this model, the total population $N$ is divided into susceptibles $S_{t}$ and infecteds $I_{t}$, with the assumption that recovered people become susceptible again. Here and below, the model parameters $\mu$, $\beta$, $\gamma$, $\sigma>0$ are assumed to be positive constants. More accurately, $\mu$ stands for the per capita death rate, $\beta$ is the disease transmission coefficient, $\gamma$ means the cure rate and $\sigma$ represents the variance of the occurrence of potentially infectious contacts between infecteds and susceptibles. In addition, $\{W_{t}\}_{t\in[0,T]}$, $T>0$ is a standard Brownian motion defined on a complete probability space $(\Omega,\mathcal{F},\mathbb{P})$. Since ${\rm d}\,[S_{t}+I_{t}]=[\mu N-\mu(S_{t}+I_{t})]\,{\rm d}t$, it can be directly deduced that $S_{t}+I_{t}=N$. Thus the original SDE system (1.1)-(1.2) can be reduced to a scalar SDE for $I_{t}$:

It was asserted in [3] that SDE (1.3) admits a unique global solution in $(0,N)$, whose dynamics behavior was also discussed. On the one hand, the extinction property of SDE (1.3), which means that the disease will die out with probability one (see Definition 4.1), was confirmed under the conditions that

• •

  $R_{0}^{S}<1$,

• •

  $\sigma^{2}\leq\tfrac{\beta}{N}$ or $\sigma^{2}>\tfrac{\beta}{N}\vee\tfrac{\beta^{2}}{2(\mu+\gamma)}$,

where

are the basic reproduction numbers for the deterministic and stochastic SIS epidemic models respectively. On the other hand, SDE (1.3) was shown to possess the persistence property, that is, the disease will stay (see Definition 4.4), when $R_{0}^{S}>1$. Since the analytic solution of SDE (1.3) is not available, people resort to numerical solutions that are expected to be convergent and preserve as many properties of the original model as possible. We mention that the numerical scheme must be carefully designed to obtain good approximations. Indeed, the well-known Euler Maruyama (EM) scheme cannot preserve the domain of the original model and the super-linear drift and diffusion coefficients of (1.3) make the resulting approximations divergent [4]. In the past few years, many researchers have made a great of efforts to design implicit [5, 6, 7, 8, 9] and explicit [10, 11, 12, 13, 14, 15, 16, 17, 18, 19] schemes for SDEs with super-linearly growing coefficients. However, boundary and dynamics preserving schemes are much less studied. In [20, 21], the authors combined a Lamperti-type transformation with the backward Euler method to derive strongly convergent and domain-preserving schemes for scalar SDEs. More recently, the authors of [22] proposed a positivity-preserving logarithmic transformed truncated Euler-Maruyama method for scalar SDEs and proved its strong and weak convergence rates. For the model (1.3), some authors [23, 24, 25, 26] relied on the Lamperti-type transformation and applied some explicit schemes to the transformed SDE that has an additive noise. Transforming back yields the desired numerical approximations, which preserves the domain and possesses a strong convergence rate of order one in [23, 24, 25] and order 1.5 in [26]. In [27, 28], a truncated Wiener process was used to construct one-half order strong schemes for the model (1.3), where restrictive assumptions were imposed on the step size $h$ to preserve the domain and dynamics behavior of the original model. It is worthwhile to point out that the numerical approximations produced by the above domain-preserving schemes cover the entire domain $(0,N)$, except for that in [25], where a truncation strategy leads to approximations living in a sub-set of $(0,N)$. Also, we mention that the dynamics behavior of the numerical approximations was investigated in [25, 27, 26, 28], but not in [23, 24]. More precisely, the proposed schemes in [25, 27, 26, 28] were proven to reproduce extinction under certain conditions, and the preservation of the persistence property was achieved in [26, 27, 28]. A comparison between numerical schemes concerning dynamics behavior with ours is made in Table 1 presented below, to provide a more comprehensive overview of the properties of these domain-preserving schemes.

An interesting question thus arises as to whether one can construct any boundary and dynamics preserving scheme for the stochastic epidemic model, without any restriction on the step size $h>0$. This is desirable, particularly in the multi-level Monte Carlo (MLMC) setting [29], where one is required to use many simulations with large discretization time step size $h>0$. It is therefore natural to look for more advanced numerical schemes that capture such properties unconditionally. The present work aims to give a positive answer, by introducing an explicit, strongly convergent scheme that preserves both the dynamics behavior and the boundaries of the analytic solutions for any given step size $h>0$. To this end, for $M\in\mathbb{N}^{+}$ we construct a uniform mesh $\{t_{k}=kh\}_{k=0}^{M}$ over $[0,T]$ with a uniform step size $h=\tfrac{T}{M}$. Our strategy is to combine a logarithm transformation with a corrected Milstein method. More precisely, by a logarithm transformation $X_{t}=\log(I_{t})$, we obtain the transformed SDE (2.2), for which corrected Milstein approximations $\{Y_{k}\}_{k\geq 0}$, given by (2.10), are designed. The correction used here guarantees that the approximations $\{Y_{k}\}_{k\geq 0}$ live in $(-\infty,\log N)$. Transforming back by $\mathcal{I}_{k}:=\text{e}^{Y_{k}}$ gives the desired approximation of $I_{t_{k}}$, which is boundary preserving and strongly convergent with order one (see Theorem 3.2, 3.3):

Throughout this paper the notation $C$ might be slightly abused to denote a generic positive constant depending on $T$ but independent of $h$, and might vary from each time of appearance. Recall that a Lamperti-type transformation $X_{t}=\log(\tfrac{I_{t}}{N-I_{t}})$ was used in [23, 26, 24, 25] to design explicit boundary preserving schemes. Different from [23, 26, 24, 25], we turn to a direct logarithm transformation $X_{t}=\log(I_{t})$, which, combined with a corrected Milstein time-stepping scheme (2.10), helps covering the entire domain $(0,N)$ and preserving the dynamics behavior of original model for any step size $h>0$. Under exactly the same conditions as required in [3], the proposed scheme reproduces the extinction and persistence properties of the original model (see Theorem 4.3, 4.6). This is also confirmed in the numerical experiments, showing that the proposed scheme performs well in dynamics preserving, even with large step sizes.

The rest of this paper is organized as follows. In the next section, the basic properties of the stochastic SIS epidemic model are stated and the numerical scheme is proposed. In Section 3, the first order strong convergence is elaborated. Section 4 is devoted to the analysis of the dynamics behavior for the proposed scheme. Numerical experiments are presented in Section 5 to verify our theoretical findings. A short conclusion is made finally in Section 6.

The present section revisits the stochastic epidemic model (1.3) and attempts to design a desirable time-stepping scheme. For the well-posedness of (1.3), we recall the following lemma, which can be found in [3, Theorem 3.1].

We utilize a logarithmic transform $X_{t}=\log I_{t}$ to turn (1.3) into

where and throughout the article $f,g:\mathbb{R}\rightarrow\mathbb{R}$ are defined as

One can readily obtain from Lemma 2.1 that

Moreover, the coefficients $f$ and $g$ satisfy the Lipschitz conditions in the domain $(-\infty,\log N)$, which will be used in the error analysis.

Proof: The assertions are trivial due to the boundedness of $\text{e}^{x}$ in $(-\infty,\log N)$. ∎

For any chosen positive integer $M\in\mathbb{N}^{+}$, a uniform mesh $\{t_{n}=nh\}_{n=0}^{M}$ is constructed with the uniform step size $h=T/M$. For convenience, we define

In addition, for $s\in[0,T]$, denote

Letting $\alpha\in(0,1]$, $\theta\geq\frac{3}{2}$ be fixed, we propose a corrected Milstein scheme $\{Y_{k}\}_{k=0,...,M}$ starting from $Y_{0}=X_{0}$ for the transformed SDE (2.2), given by

where

and $\mathbbm{1}_{A}$ is the indicator function of some set $A$, i.e., $\mathbbm{1}_{A}(\omega)=1$ when $\omega\in A$ and $\mathbbm{1}_{A}(\omega)=0$ when $\omega\notin A$. Then transforming the $\{Y_{k}\}_{k=0}^{M}$ back gives the true numerical approximations $\{\mathcal{I}_{k}\}_{k=0}^{M}$ for the original model (1.3), given by

Such a scheme is termed as the logarithmic corrected Milstein (LCM) method. In light of (2.10), one knows

and thus

In other words, the LCM method is able to preserve the domain of the stochastic epidemic model.

In the next two lemmas we establish the moment and exponential moment bounds of the solution for the transformed SDE (2.2) and its numerical approximation (2.10).

Proof: The assertion follows due to (2.5) and (2.13). ∎

Proof: The assertion can be derived directly by definitions of $\{X_{t}\}_{t\in[0,T]}$, $\{Y_{k}\}_{k=0}^{M}$ and Lemma 2.3. ∎

In view of Lemmas 2.2, 2.3 and 2.4, moments of the analitical and numerical solutions we encounter later are all bounded, which will not be emphasized further.

In this section we attempt to present the analysis of the strong convergence rate for the proposed LCM scheme. To begin with, we recall first a discrete version of the Burkholder-Davis-Gundy inequality, quoted from [30, Lemma 4.1].

Here and below we denote $\|\xi\|_{\mathcal{L}^{q}}:=\big{(}\mathbb{E}[|\xi|^{q}]\big{)}^{\tfrac{1}{q}}\in[0,+\infty]$ for any $q\in[1,+\infty)$. Now we are ready to carry out error analysis on the mesh grids.

Proof: It is easy to deduce that

Recalling (2.2) and (2.10) we have

where we denote

Next we estimate $R_{k}$ first. For $k=0,...,M-1$, the Itô formula shows that

where we split $R_{k}$ into two parts as follows:

The moment inequality and the Jensen inequality together with Lemma 2.2 yield that, for any $q\geq 1$

where the Hölder inequality is also used in the last inequality. Therefore, the boundedness of moments ensure that, for any $q\geq 1$ and $k=0,1,...,M-1$

and thus

Denote for brevity

Thanks to Lemma 2.2, we have

Plugging (3.4) into (3.3) and by the Young inequality we get, for any $k=0,...,M-1$,

By iteration, we arrive at

For any $n=0,...,M-1$ and $p\geq 2$, raising both sides of (3.14) to the power of $p$, taking supremum and expectation as well as applying the Jensen inequality we get

By denoting $\xi_{i}^{(1)}:=(1+3h)^{-i}e_{i}B_{i}$ and $\xi_{i}^{(2)}:=(1+3h)^{-i}e_{i}R_{i}^{(2)},$ it is not difficult to check that $\mathbb{E}\big{[}\xi_{i+1}^{(1)}|\xi_{1}^{(1)},...,\xi_{i}^{(1)}\big{]}=0$ and $\mathbb{E}\big{[}\xi_{i+1}^{(2)}|\xi_{1}^{(2)},...,\xi_{i}^{(2)}\big{]}=0.$ Moreover, $\{\Xi_{k}^{(1)}\}_{k=1}^{M}:=\Big{\{}\sum\limits_{i=0}^{k-1}\xi_{i}^{(1)}\Big{\}}_{k=1}^{M}$ and $\{\Xi_{k}^{(2)}\}_{k=1}^{M}:=\Big{\{}\sum\limits_{i=0}^{k-1}\xi_{i}^{(2)}\Big{\}}_{k=1}^{M}$ are square-integrable discrete time martingales. Hence the Doob discrete martingale inequality, Lemma 3.1 and the Jensen inequality can be applied to infer that

With the help of the Young inequality additionally, we have

Inserting (3.16)-(3.17) into (3.15) we arrive at

Finally, the Gronwall inequality and the Lyapunov inequality finish the proof. ∎

Equipped with the above convergence result for the approximations of the transformed SDE, we therefore obtain a convergence rate of order one for the LCM scheme.

Proof: The Hölder inequality and boundedness of the exponential moment together imply that for $p\geq 2$,

Thus by the Lyapunov inequality the assertion holds for any $q\geq 1$. ∎

In this section we investigate the ability of the proposed scheme to reproduce significant dynamics behavior of the stochastic epidemic model. More precisely, we examine the extinction and persistence properties of the time discretization and show that, under exactly the same conditions as required in [3], the proposed scheme is able to preserve the dynamics behavior of the original model for any step size $h>0$.