# An unconditional boundary and dynamics preserving scheme for the stochastic
epidemic model ††thanks: This work was supported by Natural Science Foundation
of China (12071488, 11971488) and Natural Science Foundation of Hunan Province
(2020JJ2040).
Emails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
Ruishu Liu, Xiaojie Wang, Lei Dai
School of Mathematics and Statistics, HNP-LAMA, Central South University,
Changsha, China

###### Abstract

By combining a logarithm transformation with a corrected Milstein-type method,
the present article proposes an explicit, unconditional boundary and dynamics
preserving scheme for the stochastic susceptible-infected-susceptible (SIS)
epidemic model that takes value in $(0,N)$. The scheme applied to the model is
first proved to have a strong convergence rate of order one. Further, the
dynamic behaviors are analyzed for the numerical approximations and it is
shown that the scheme can unconditionally preserve both the domain and the
dynamics of the model. More precisely, the proposed scheme gives numerical
approximations living in the domain $(0,N)$ and reproducing the extinction and
persistence properties of the original model for any time discretization step-
size $h>0$, without any additional requirements on the model parameters.
Numerical experiments are presented to verify our theoretical results.

Keywords:  Stochastic SIS epidemic model; Boundary preserving; Extinction;
Persistence; Explicit scheme; Order 1 strong convergence

##  1 Introduction

As the divergent genotypes of the Corona Virus Disease 2019 (COVID-19) virus,
the swine flu and other influenza viruses still bring out massive infections
and even deaths this year, it is of great importance to shed light on the
pattern of the diseases spread and figure out efficient countermeasures
against contagions. In the past few decades, mathematic models have been
constructed to describe the evolution of different kinds of epidemics. The
classical Kermack-McKendrick model [15] is 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. This is also called the susceptible-infected-removed (SIR) model in
similar situations as it literally means. However, there are kinds of viruses
that possess extremely strong mutabilities, especially the single-stranded RNA
viruses, which are also the prime culprits that have made the world suffer in
the past three years. Due to their instability, cured people can rarely get
long-time immunity and turns to be susceptible to some mutations of the same
original virus. For these cases, the susceptible-infected-susceptible (SIS)
model is constructed in [20]. However, the fact that the world we live in is
hardly predictable is overlooked in the deterministic models. Moreover, the
incompleteness of knowledge for the systems is also an impetus for adding
randomness to the epidemic models. The stochastic version of the SIS model
proposed by Gray et al. in [9] is presented as the following Itô SDEs:

$\displaystyle{\rm d}S_{t}$ $\displaystyle=\left[\mu N-\beta S_{t}I_{t}+\gamma
I_{t}-\mu S_{t}\right]\,{\rm d}t-\sigma S_{t}I_{t}\,{\rm d}W_{t}\,,$ (1.1)
$\displaystyle{\rm d}I_{t}$ $\displaystyle=\left[\beta
S_{t}I_{t}-(\mu+\gamma)I_{t}\right]\,{\rm d}t+\sigma S_{t}I_{t}\,{\rm
d}W_{t}\,,$ (1.2)

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 that
$S_{0}+I_{0}$ equals the total population $N$. In this model, the total
population is divided into susceptibles and infecteds, with the assumption
that cured people can be infected again. Here $\mu$, $\gamma$, $\beta$ and
$\sigma$ are positive constants representing the per capita death rate, the
cure rate among infecteds, the disease transmission coefficient and the
variance of the number of potentially infectious contacts that a single
infected individual makes with another individual in the time period
$[t,t+{\rm d}t)$. $\\{W_{t}\\}_{t\in[0,T]}$, where $T$ is an arbitrary
positive constant, stands for a standard Brownian motion defined on a complete
probability space $(\Omega,\mathcal{F},\mathbb{P})$. Observed that ${\rm
d}\,[S_{t}+I_{t}]=[\mu N-\mu(S_{t}+I_{t})]\,{\rm d}t$, it can be easily
deduced that $S_{t}+I_{t}=N$, which reduces the original SDE system
(1.1)-(1.2) to the following scalar SDE for $I_{t}$:

${\rm d}I_{t}=I_{t}(\beta N-\mu-\gamma-\beta I_{t})\,{\rm d}t+\sigma
I_{t}(N-I_{t})\,{\rm d}W_{t},\quad t\in(0,T]$ (1.3)

with initial value $I_{0}\in(0,N)$. It was asserted in [9] that SDE (1.3)
admits a unique global solution in $(0,N)$, whose dynamic behaviors were also
discussed. On the one hand, the extinction of SDE (1.3) 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

$\displaystyle R_{0}^{D}:=\tfrac{\beta N}{\mu+\gamma}\quad\text{and}\quad
R_{0}^{S}:=R_{0}^{D}-\tfrac{\sigma^{2}N^{2}}{2(\mu+\gamma)}$ (1.4)

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 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 turns out to be a bad candidate since it cannot preserve
the domain of the original model and the super-linear drift and diffusion
coefficients of (1.3) make the resulting approximations divergent [13]. The
past few years have also witnessed a prosper growth of implicit [3, 4, 10, 26,
2] and explicit [14, 22, 25, 19, 24, 12, 5, 23, 6, 31] schemes designed for
SDEs with super-linearly growing coefficients. However, boundary and dynamics
preserving schemes are much less studied. In [1, 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 [16] proposed a positivity-preserving logarithmic transformed
truncated Euler-Maruyama method for scalar SDEs and proved its strong and weak
convergence rates. For the particular model (1.3), some authors [7, 28, 27,
17] 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 [29, 30], 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 dynamic behaviors of the original model.

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 [8], where one is
required to use many simulations with large discretization time step $h>0$. It
is therefore natural to look for more advanced numerical schemes that capture
such a property unconditionally. The present work aims to give a positive
answer, by introducing an explicit, strongly convergent scheme that preserves
both the dynamic behavior and the boundaries of the analytic solutions for any
given step-size $h>0$. Indeed, 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), to which a corrected Milstein method (2.10) is applied. Transforming
back by $I_{k}=\text{e}^{Y_{k}}$ gives the desired approximation (2.12) that
is boundary preserving and strongly convergent with order one (see Theorem
3.2, 3.3):

$\mathbb{E}\Big{[}\sup_{k=0,...,M}\big{|}I_{t_{k}}-I_{k}\big{|}^{p}\Big{]}\leq
Ch^{p},$ (1.5)

where $h=\tfrac{T}{M}$ is the uniform step-size. 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 [7, 17, 28, 27] to design
explicit boundary preserving schemes. Different from [7, 17, 28, 27], we turn
to a direct logarithm transformation $X_{t}=\log(I_{t})$, which combined with
a corrected Milstein time-stepping scheme (2.10) helps preserve the both the
boundary and dynamic behaviors of the original model for any step-size $h>0$.
Under exactly the same conditions as required in [9], 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 dynamic behaviors 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.

##  2 The stochastic epidemic model and the proposed scheme

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 [9, Theorem 3.1].

######  Lemma 2.1.

For any initial value $I_{0}\in(0,N)$, the SDE (1.3) has a unique global
positive solution $I_{t}\in(0,N)$ for all $t\geq 0$ with probability one,
namely,

$\mathbb{P}(I_{t}\in(0,N),\forall t\geq 0)=1.$ (2.1)

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

$\left\\{\begin{array}[]{l}X_{t}-X_{0}=\int_{0}^{t}f(X_{s}){\rm
d}s+\int_{0}^{t}g(X_{s}){\rm d}W_{s},\ t\in[0,T],\\\ X_{0}=\log
I_{0},\end{array}\right.$ (2.2)

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

$\displaystyle f(x)$
$\displaystyle:=-\tfrac{1}{2}\sigma^{2}\text{e}^{2x}+(\sigma^{2}N-\beta)\text{e}^{x}+(\beta
N-\mu-\gamma-\tfrac{1}{2}\sigma^{2}N^{2}),$ (2.3) $\displaystyle g(x)$
$\displaystyle:=\sigma(N-\text{e}^{x}).$ (2.4)

One can readily obtain from Lemma 2.1 that

$\mathbb{P}\big{(}X_{t}\in(-\infty,\log N),\,\forall t\geq 0\big{)}=1.$ (2.5)

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

######  Lemma 2.2.

For any $x,y\in(-\infty,\log N)$, there exists a non-negative constant $L$
such that

$|f(x)-f(y)|\vee|g(x)-g(y)|\leq L|x-y|$ (2.6)

and

$|g(x)g^{\prime}(x)-g(y)g^{\prime}(y)|\leq L|x-y|.$ (2.7)

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 use $X_{k}$ to denote $X_{t_{k}}$ for any $k=0,...,M$ and
define

$\displaystyle\Delta W_{k}:=W_{t_{k+1}}-W_{t_{k}},\quad\forall k=0,...,M-1.$
(2.8)

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

$\displaystyle\kappa(s):=\sup\limits_{k=0,...,M}\\{t_{k}:t_{k}\leq s\\}.$
(2.9)

Letting $\alpha\in(0,1]$, $\theta\geq\frac{3}{2}$ be fixed, we propose a
corrected Milstein scheme for the transformed SDE (2.2), given by

$\left\\{\begin{array}[]{l}\bar{Y}_{k+1}=Y_{k}+f(Y_{k})h+g(Y_{k})\Delta
W_{k}+g(Y_{k})g^{\prime}(Y_{k})\Delta\zeta_{k},\\\
{Y}_{k+1}=\bar{Y}_{k+1}\wedge(\log N-\alpha h^{\theta}),\quad
k=0,1,...,M-1,\end{array}\right.$ (2.10)

where for brevity we denote

$\Delta\zeta_{k}:=\int_{t_{k}}^{t_{k+1}}\int_{t_{k}}^{s}\,{\rm d}W_{r}{\rm
d}W_{s}=\tfrac{1}{2}\Delta W_{k}^{2}-\tfrac{1}{2}h.$ (2.11)

Then transforming the numerical approximation $Y_{k}$ back gives numerical
approximations $I_{k}$ for the original model (1.3), given by

$I_{k}=e^{Y_{k}},\quad k=0,...,M.$ (2.12)

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

$\mathbb{P}\big{(}Y_{k}\in(-\infty,\log N),\,k=0,1,2,...\big{)}=1$ (2.13)

and thus

$\mathbb{P}\big{(}I_{k}\in(0,N),\,k=0,1,2,...\big{)}=1.$ (2.14)

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).

######  Lemma 2.3 (Exponential moment bounds).

Let $T>0$ and $X_{t}$, $Y_{k}$ be defined as (2.2), (2.10). For any $q\geq 0$,
it holds that

$\sup_{t\in[0,T]}\mathbb{E}\big{[}{\rm
e}^{qX_{t}}\big{]}\bigvee\sup_{k=0,...,M}\mathbb{E}\big{[}{\rm
e}^{qY_{k}}\big{]}<\infty.$ (2.15)

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

######  Lemma 2.4 (Moment bounds).

Let $T>0$ and $X_{t}$, $Y_{k}$ be defined as (2.2), (2.10). For any $q\geq 0$,
it holds that

$\sup_{t\in[0,T]}\mathbb{E}\big{[}|X_{t}|^{q}\big{]}\bigvee\sup_{k=0,...,M}\mathbb{E}\big{[}|Y_{k}|^{q}\big{]}<\infty.$
(2.16)

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.

##  3 The first-order strong convergence of the scheme

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 [11, Lemma 4.1].

######  Lemma 3.1.

Let $M\in\mathbb{N}$ and $\xi_{1},...,\xi_{M}:\Omega\rightarrow\mathbb{R}$ be
$\mathcal{F}/\mathcal{B}(\mathbb{R})$-measurable mappings with
$\mathbb{E}\big{[}|\xi_{n}|^{2}\big{]}<+\infty$ for all $n\in\\{1,...,M\\}$
and with $\mathbb{E}\big{[}\xi_{n+1}|\xi_{1},...,\xi_{n}\big{]}=0$ for all
$n\in\\{1,...,M-1\\}$. Then

$\displaystyle\|\xi_{1}+\cdots+\xi_{M}\|_{\mathcal{L}^{p}}\leq
C_{p}\cdot\Big{(}\|\xi_{1}\|_{\mathcal{L}^{p}}^{2}+\cdots+\|\xi_{M}\|_{\mathcal{L}^{p}}^{2}\Big{)}^{\tfrac{1}{2}}$
(3.1)

for every $p\in[2,+\infty)$, where $C_{p},p\in[2,+\infty)$ are universal
constants.

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.

######  Theorem 3.2.

Let $T>0$ and $X_{t}$, $Y_{k}$ be defined as (2.2), (2.10). For any $q\geq 1$,
it holds that

$\mathbb{E}\Big{[}\sup\limits_{k=0,...,M}|X_{k}-Y_{k}|^{q}\Big{]}\leq Ch^{q}.$
(3.2)

Proof: It is easy to deduce that

$\displaystyle|X_{k+1}-Y_{k+1}|^{2}$
$\displaystyle=\mathbbm{1}_{\\{\bar{Y}_{k+1}\leq\log N-\alpha
h^{\theta}\\}}|X_{k+1}-\bar{Y}_{k+1}|^{2}$ (3.3)
$\displaystyle\quad+\mathbbm{1}_{\\{\bar{Y}_{k+1}>\log N-\alpha
h^{\theta}\\}}\cdot\mathbbm{1}_{\\{X_{k+1}\leq\log N-\alpha
h^{\theta}\\}}|X_{k+1}-(\log N-\alpha h^{\theta})|^{2}$
$\displaystyle\quad+\mathbbm{1}_{\\{\bar{Y}_{k+1}>\log N-\alpha
h^{\theta}\\}}\cdot\mathbbm{1}_{\\{\log N>X_{k+1}>\log N-\alpha
h^{\theta}\\}}|X_{k+1}-(\log N-\alpha h^{\theta})|^{2}$
$\displaystyle\leq|X_{k+1}-\bar{Y}_{k+1}|^{2}+\alpha^{2}h^{2\theta}.$

Recalling (2.2) and (2.10) we have

$\displaystyle X_{k+1}-\bar{Y}_{k+1}$
$\displaystyle=X_{k}-Y_{k}+\Big{(}f(X_{k})-f(Y_{k})\Big{)}h+\Big{(}g(X_{k})-g(Y_{k})\Big{)}\Delta
W_{k}$ (3.4)
$\displaystyle\quad+\Big{(}g(X_{k})g^{\prime}(X_{k})-g(Y_{k})g^{\prime}(Y_{k})\Big{)}\Delta\zeta_{k}+R_{k},$

where we denote

$\displaystyle R_{k}$
$\displaystyle:=\int_{t_{k}}^{t_{k+1}}f(X_{s})-f(X_{k})\,{\rm
d}s+\int_{t_{k}}^{t_{k+1}}g(X_{s})-g(X_{k})\,{\rm
d}W_{s}-\int_{t_{k}}^{t_{k+1}}\int_{t_{k}}^{s}g(X_{k})g^{\prime}(X_{k})\,{\rm
d}W_{r}{\rm d}W_{s}.$ (3.5)

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

$\displaystyle R_{k}$
$\displaystyle=\int_{t_{k}}^{t_{k+1}}\int_{t_{k}}^{s}f^{\prime}(X_{r})f(X_{r})+\tfrac{1}{2}f^{\prime\prime}(X_{r})g^{2}(X_{r})\,{\rm
d}r{\rm
d}s+\int_{t_{k}}^{t_{k+1}}\int_{t_{k}}^{s}f^{\prime}(X_{r})g(X_{r})\,{\rm
d}W_{r}{\rm d}s$ (3.6)
$\displaystyle\quad+\int_{t_{k}}^{t_{k+1}}\int_{t_{k}}^{s}g^{\prime}(X_{r})f(X_{r})+\tfrac{1}{2}g^{\prime\prime}(X_{r})g^{2}(X_{r})\,{\rm
d}r{\rm d}W_{s}$
$\displaystyle\quad+\int_{t_{k}}^{t_{k+1}}\int_{t_{k}}^{s}g(X_{r})g^{\prime}(X_{r})-g(X_{k})g^{\prime}(X_{k})\,{\rm
d}W_{r}{\rm d}W_{s}$ $\displaystyle=R_{k}^{(1)}+R_{k}^{(2)},$

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

$\displaystyle R_{k}^{(1)}$
$\displaystyle:=\int_{t_{k}}^{t_{k+1}}\int_{t_{k}}^{s}f^{\prime}(X_{r})f(X_{r})+\tfrac{1}{2}f^{\prime\prime}(X_{r})g^{2}(X_{r})\,{\rm
d}r{\rm d}s$ (3.7) $\displaystyle R_{k}^{(2)}$
$\displaystyle:=\int_{t_{k}}^{t_{k+1}}\int_{t_{k}}^{s}f^{\prime}(X_{r})g(X_{r})\,{\rm
d}W_{r}{\rm
d}s+\int_{t_{k}}^{t_{k+1}}\int_{t_{k}}^{s}g^{\prime}(X_{r})f(X_{r})+\tfrac{1}{2}g^{\prime\prime}(X_{r})g^{2}(X_{r})\,{\rm
d}r{\rm d}W_{s}$
$\displaystyle\quad+\int_{t_{k}}^{t_{k+1}}\int_{t_{k}}^{s}g(X_{r})g^{\prime}(X_{r})-g(X_{k})g^{\prime}(X_{k})\,{\rm
d}W_{r}{\rm d}W_{s}.$ (3.8)

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

$\displaystyle\mathbb{E}\Bigg{[}\bigg{|}\int_{t_{k}}^{t_{k+1}}\int_{t_{k}}^{s}g(X_{r})g^{\prime}(X_{r})-g(X_{k})g^{\prime}(X_{k})\,{\rm
d}W_{r}{\rm d}W_{s}\bigg{|}^{2q}\Bigg{]}$ (3.9) $\displaystyle\leq
h^{2q-2}\cdot\int_{t_{k}}^{t_{k+1}}\int_{t_{k}}^{s}\mathbb{E}\bigg{[}\Big{|}g(X_{r})g^{\prime}(X_{r})-g(X_{k})g^{\prime}(X_{k})\Big{|}^{2q}\bigg{]}\,{\rm
d}r{\rm d}s$ $\displaystyle\leq
Lh^{2q-2}\cdot\mathbb{E}\Bigg{[}\int_{t_{k}}^{t_{k+1}}\int_{t_{k}}^{s}\Big{|}X_{r}-X_{k}\Big{|}^{2q}\,{\rm
d}r{\rm d}s\Bigg{]}$ $\displaystyle\leq
CLh^{4q-3}\int_{t_{k}}^{t_{k+1}}\int_{t_{k}}^{s}\int_{t_{k}}^{r}\mathbb{E}\bigg{[}\Big{|}f(X_{u})\Big{|}^{2q}\bigg{]}\,{\rm
d}u{\rm d}r{\rm d}s$
$\displaystyle\quad+CLh^{3q-3}\int_{t_{k}}^{t_{k+1}}\int_{t_{k}}^{s}\int_{t_{k}}^{r}\mathbb{E}\bigg{[}\Big{|}g(X_{u})\Big{|}^{2q}\bigg{]}\,{\rm
d}u{\rm d}r{\rm d}s,$

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$

$\displaystyle\mathbb{E}\Big{[}\big{|}R_{k}^{(1)}\big{|}^{2q}\Big{]}\leq
Ch^{4q},\quad\mathbb{E}\Big{[}\big{|}R_{k}^{(2)}\big{|}^{2q}\Big{]}\leq
Ch^{3q}$ (3.10)

and thus

$\mathbb{E}\Big{[}\big{|}R_{k}\big{|}^{2q}\Big{]}\leq Ch^{3q}.$ (3.11)

Denote for brevity

$\displaystyle e_{k}$ $\displaystyle:=X_{k}-Y_{k},$ $\displaystyle\quad A_{k}$
$\displaystyle:=\Big{(}f(X_{k})-f(Y_{k})\Big{)}h,$ $\displaystyle B_{k}$
$\displaystyle:=\Big{(}g(X_{k})-g(Y_{k})\Big{)}\Delta W_{k},$
$\displaystyle\quad D_{k}$
$\displaystyle:=\Big{(}g(X_{k})g^{\prime}(X_{k})-g(Y_{k})g^{\prime}(Y_{k})\Big{)}\Delta\zeta_{k}.$

Thanks to Lemma 2.2, we have

$\displaystyle|A_{k}|\leq Lh|e_{k}|,\quad|B_{k}|\leq L|e_{k}||\Delta
W_{k}|,\quad|D_{k}|\leq L|e_{k}||\Delta\zeta_{k}|.$ (3.12)

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

$\displaystyle e^{2}_{k+1}$ $\displaystyle\leq
e^{2}_{k}+A^{2}_{k}+B^{2}_{k}+D^{2}_{k}+R^{2}_{k}+2e_{k}A_{k}+2e_{k}B_{k}+2e_{k}D_{k}+2e_{k}R_{k}$
(3.13)
$\displaystyle\quad+2A_{k}B_{k}+2A_{k}R_{k}+2A_{k}D_{k}+2B_{k}D_{k}+2B_{k}R_{k}+2D_{k}R_{k}+\alpha^{2}h^{2\theta}$
$\displaystyle\leq(1+3h)e^{2}_{k}+\big{(}4+h^{-1}\big{)}(A^{2}_{k}+D^{2}_{k})+4\big{(}B^{2}_{k}+R^{2}_{k}\big{)}+2e_{k}B_{k}$
$\displaystyle\quad+h^{-1}|R^{(1)}_{k}|^{2}+2e_{k}R^{(2)}_{k}+\alpha^{2}h^{2\theta}.$

By iteration, we arrive at

$\displaystyle e^{2}_{k+1}$
$\displaystyle\leq(1+3h)^{k}\sum_{i=0}^{k}\tfrac{1}{(1+3h)^{i}}\Big{[}(4+h^{-1})(A^{2}_{i}+D^{2}_{i})+4\big{(}B^{2}_{i}+R^{2}_{i}\big{)}+h^{-1}|R^{(1)}_{i}|^{2}\Big{]}$
(3.14)
$\displaystyle\quad+(1+3h)^{k}\sum_{i=0}^{k}\tfrac{1}{(1+3h)^{i}}\Big{[}2e_{i}B_{i}+2e_{i}R^{(2)}_{i}\Big{]}+Ch^{2\theta-1}.$

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

$\displaystyle\mathbb{E}\Big{[}\sup_{k=0,...,n}|e_{k+1}|^{2p}\Big{]}$
$\displaystyle\leq
Ch\sum_{i=0}^{n}\mathbb{E}\Big{[}|e_{i}|^{2p}\Big{]}+C(1+3h)^{np}\cdot\mathbb{E}\Bigg{[}\sup_{k=0,...,n}\bigg{|}\sum_{i=0}^{k}\tfrac{1}{(1+3h)^{i}}2e_{i}B_{i}\bigg{|}^{p}\Bigg{]}$
(3.15)
$\displaystyle\quad+C(1+3h)^{np}\cdot\mathbb{E}\Bigg{[}\sup_{k=0,...,n}\bigg{|}\sum_{i=0}^{k}\tfrac{1}{(1+3h)^{i}}2e_{i}R^{(2)}_{i}\bigg{|}^{p}\Bigg{]}+Ch^{2p}.$

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

$\displaystyle\mathbb{E}\Bigg{[}\bigg{|}\sup_{k=0,...,n}\sum_{i=0}^{k}\tfrac{1}{(1+3h)^{i}}2e_{i}B_{i}\bigg{|}^{p}\Bigg{]}$
$\displaystyle\leq
C\cdot(\tfrac{1}{h})^{\frac{p}{2}-1}\sum_{i=0}^{n}\mathbb{E}\bigg{[}\Big{|}2e_{i}B_{i}\Big{|}^{p}\bigg{]}\leq
Ch\sum_{i=0}^{n}\mathbb{E}\Big{[}|e_{i}|^{2p}\Big{]}.$ (3.16)

With the help of the Young inequality additionally, we have

$\displaystyle\mathbb{E}\Bigg{[}\bigg{|}\sup_{k=0,...,n}\sum_{i=0}^{k}\tfrac{1}{(1+3h)^{i}}2e_{i}R_{i}^{(2)}\bigg{|}^{p}\Bigg{]}$
$\displaystyle\leq
C\cdot(\tfrac{1}{h})^{\frac{p}{2}-1}\sum_{i=0}^{n}\mathbb{E}\bigg{[}\Big{|}2e_{i}R_{i}^{(2)}\Big{|}^{p}\bigg{]}\leq
Ch\sum_{i=0}^{n}\mathbb{E}\Big{[}|e_{i}|^{2p}\Big{]}+Ch^{2p}.$ (3.17)

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

$\displaystyle\mathbb{E}\Big{[}\sup_{k=0,...,n}|e_{k+1}|^{2p}\Big{]}\leq
Ch\sum_{i=1}^{n}\mathbb{E}\Big{[}|e_{i}|^{2p}\Big{]}+Ch^{2p}<\infty.$ (3.18)

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.

######  Theorem 3.3.

Let $T>0$ and $I_{t}$, $I_{k}$ be defined as (1.3), (2.12). For any $q\geq 1$,
there exists some positive constant $C$ independent of $h$ such that

$\mathbb{E}\Big{[}\sup_{k=0,...,M}\big{|}I_{t_{k}}-I_{k}\big{|}^{q}\Big{]}\leq
Ch^{q}.$ (3.19)

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

$\displaystyle\mathbb{E}\Big{[}\sup_{k=0,...,M}\big{|}e^{X_{k}}-e^{Y_{k}}\big{|}^{p}\Big{]}\leq
C\mathbb{E}\Big{[}\sup_{k=0,...,M}\big{|}X_{k}-Y_{k}\big{|}^{p}\Big{]}\leq
Ch^{p}.$ (3.20)

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

##  4 Dynamic behaviors of numerical approximations

In this section we investigate the ability of the proposed scheme to reproduce
two significant dynamic behaviors 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
[9], the proposed scheme is able to preserve these two dynamic behaviors of
the original model for any step-size $h>0$.

###  4.1 Unconditional extinction preserving

As the first dynamic behavior of the stochastic model, we focus on the
extinction property that is defined as follows.

######  Definition 4.1 (Extinction).

If the unique positive solution $I_{t}$, $t\geq 0$ of the stochastic SIS model
tends to zero in almost sure sense as $t\rightarrow\infty$, that is,

$\lim\limits_{t\rightarrow\infty}I_{t}=0\quad a.s.,$ (4.1)

then we say that the infected population system has the extinction property.
Similarly, the numerical approximation $\\{I_{k}\\}_{k\geq 0}$ is said to have
the extinction property if

$\lim\limits_{kh\rightarrow\infty}I_{k}=0\quad a.s..$ (4.2)

The next theorem presents some sufficient conditions to ensure the extinction
property of the original stochastic SIS model (1.3), which is indeed quote
from [9, Theorems 4.1, 4.3].

######  Theorem 4.2.

Let the basic reproduction number of the stochastic SIS model (1.3) be less
than one, i.e.,

$\displaystyle R_{0}^{S}=\tfrac{\beta
N}{\mu+\gamma}-\tfrac{\sigma^{2}N^{2}}{2(\mu+\gamma)}<1.$ (4.3)

For any given initial value $I_{0}\in(0,N)$, we have the following extinction
properties.

  1. (i)

If $\sigma^{2}\leq\frac{\beta}{N}$, then

$\displaystyle\limsup\limits_{t\rightarrow\infty}\frac{\log{I_{t}}}{t}\leq\beta
N-\tfrac{1}{2}\sigma^{2}N^{2}-\mu-\gamma<0\quad a.s..$ (4.4)

  2. (ii)

If
$\sigma^{2}>\max\Big{\\{}\frac{\beta}{N},\,\frac{\beta^{2}}{2(\mu+\gamma)}\Big{\\}}$,
then

$\displaystyle\limsup\limits_{t\rightarrow\infty}\frac{\log{I_{t}}}{t}\leq\tfrac{\beta^{2}}{2\sigma^{2}}-\mu-\gamma<0\quad
a.s..$ (4.5)

We mention that the assertions (4.4) and (4.5) imply that the solution $I_{t}$
of the stochastic SIS model tends to zero exponentially almost surely as time
grows. In other words, the disease dies out with probability one, as time
tends to infinity. In what follows, we investigate under which conditions the
numerical approximation will admit the extinction property in the sense that
$I_{k}$ tends to zero exponentially almost surely, which is to be stated in
the next theorem.

######  Theorem 4.3.

Let the condition (4.3) hold, i.e., $R_{0}^{S}<1$. For any given initial value
$I_{0}\in(0,N)$ and for any step-size $h>0$, the numerical approximations
$I_{k}$ possess the following extinction properties.

  1. (i)

If $\sigma^{2}\leq\frac{\beta}{N}$, then

$\displaystyle\limsup\limits_{kh\rightarrow\infty}\frac{\log{I_{k}}}{kh}\leq\beta
N-\tfrac{1}{2}\sigma^{2}N^{2}-\mu-\gamma<0\quad a.s..$ (4.6)

  2. (ii)

If
$\sigma^{2}>\max\Big{\\{}\frac{\beta}{N},\,\frac{\beta^{2}}{2(\mu+\gamma)}\Big{\\}}$,
then

$\displaystyle\limsup\limits_{kh\rightarrow\infty}\frac{\log{I_{k}}}{kh}\leq\tfrac{\beta^{2}}{2\sigma^{2}}-\mu-\gamma<0\quad
a.s..$ (4.7)

In other words, the numerical approximations produced by (2.10) tends to zero
exponentially almost surely as time evolves.

Proof: Recall first that

$\displaystyle f(x)$
$\displaystyle=-\tfrac{1}{2}\sigma^{2}\text{e}^{2x}+(\sigma^{2}N-\beta)\text{e}^{x}+(\beta
N-\mu-\gamma-\tfrac{1}{2}\sigma^{2}N^{2}),\quad x\in(-\infty,\log N).$

In the first case that $\sigma^{2}\leq\tfrac{\beta}{N}$, clearly
$f(x)\leq\beta N-\mu-\gamma-\frac{1}{2}\sigma^{2}N^{2}<0$. In the second case
that
$\sigma^{2}>\max\Big{\\{}\frac{\beta}{N},\,\frac{\beta^{2}}{2(\mu+\gamma)}\Big{\\}}$,
$f(x)$ reaches its maximum value $\tfrac{\beta^{2}}{2\sigma^{2}}-\mu-\gamma<0$
when $e^{x}=\tfrac{\sigma^{2}N-\beta}{\sigma^{2}}$. By denoting

$\displaystyle f_{max}^{\sigma}:=\begin{cases}\beta
N-\mu-\gamma-\frac{1}{2}\sigma^{2}N^{2},\quad\sigma^{2}\leq\tfrac{\beta}{N};\\\
\tfrac{\beta^{2}}{2\sigma^{2}}-\mu-\gamma,\quad\sigma^{2}>\max\Big{\\{}\frac{\beta}{N},\,\frac{\beta^{2}}{2(\mu+\gamma)}\Big{\\}},\end{cases}$
(4.8)

in any case we have

$\displaystyle f(x)\leq f_{max}^{\sigma}<0.$ (4.9)

Recall that

$\left\\{\begin{array}[]{l}\bar{Y}_{k+1}=Y_{k}+f(Y_{k})h+g(Y_{k})\Delta
W_{k}+g(Y_{k})g^{\prime}(Y_{k})\Delta\zeta_{k},\\\
{Y}_{k+1}=\bar{Y}_{k+1}\wedge(\log N-\alpha h^{\theta}).\end{array}\right.$

It follows from (4.9) and by induction that

$\displaystyle\bar{Y}_{k+1}$ $\displaystyle\leq
Y_{k}+hf_{max}^{\sigma}+g(Y_{k})\Delta
W_{k}+g(Y_{k})g^{\prime}(Y_{k})\Delta\zeta_{k}$ $\displaystyle\leq
Y_{0}+(k+1)hf_{max}^{\sigma}+\sum_{i=0}^{k}g(Y_{i})\Delta
W_{i}+\sum_{i=0}^{k}g(Y_{i})g^{\prime}(Y_{i})\Delta\zeta_{i},$ (4.10)

and thus

$\displaystyle Y_{k+1}\leq\bar{Y}_{k+1}\leq
Y_{0}+(k+1)hf_{max}^{\sigma}+\sum_{i=0}^{k}g(Y_{i})\Delta
W_{i}+\sum_{i=0}^{k}g(Y_{i})g^{\prime}(Y_{i})\Delta\zeta_{i}.$ (4.11)

Therefore,

$\displaystyle\limsup\limits_{kh\rightarrow\infty}\frac{Y_{k}}{kh}$
$\displaystyle\leq
f_{max}^{\sigma}+\limsup\limits_{kh\rightarrow\infty}\frac{1}{kh}\sum_{i=0}^{k-1}g(Y_{i})\Delta
W_{i}+\limsup\limits_{kh\rightarrow\infty}\frac{1}{kh}\sum_{i=0}^{k-1}g(Y_{i})g^{\prime}(Y_{i})\Delta\zeta_{i}\quad
a.s..$

Further, note that

$\Xi^{(3)}_{t}:=\int_{0}^{t}g(Y_{\kappa(s)}){\rm d}W_{s},\
\Xi^{(4)}_{t}:=\int_{0}^{t}g^{\prime}(Y_{\kappa(s)})g(Y_{\kappa(s)})(W_{s}-W_{\kappa(s)}){\rm
d}W_{s}$

are both square-integrable continuous time martingales. By the large number
theorem (see e.g., [18, Theorem 3.4.]) for martingales we have

$\lim\limits_{t\rightarrow\infty}\frac{\Xi^{(3)}_{t}}{t}=\lim\limits_{t\rightarrow\infty}\frac{\Xi^{(4)}_{t}}{t}=0,\
a.s.,$

which yields that

$\displaystyle\limsup\limits_{kh\rightarrow\infty}\frac{1}{kh}\sum_{i=0}^{k-1}g(Y_{i})\Delta
W_{i}=\limsup\limits_{kh\rightarrow\infty}\frac{1}{kh}\sum_{i=0}^{k-1}g(Y_{i})g^{\prime}(Y_{i})\Delta\zeta_{i}=0\quad
a.s..$ (4.12)

Therefore

$\displaystyle\limsup\limits_{kh\rightarrow\infty}\frac{Y_{k}}{kh}\leq
f_{max}^{\sigma}<0\quad a.s..$ (4.13)

Noting that $\log I_{k}=Y_{k},$ we finally arrive at

$\displaystyle\limsup\limits_{kh\rightarrow\infty}\frac{\log I_{k}}{kh}$
$\displaystyle=\limsup\limits_{kh\rightarrow\infty}\frac{Y_{k}}{kh}\leq
f_{max}^{\sigma}<0\quad a.s..$

The proof is thus completed. ∎

Theorem 4.3 shows that the numerical approximation produced by the LCM scheme
is able to reproduce the extinction property of the original model for any
step-size $h>0$.

###  4.2 Unconditional persistence preserving

In this subsection, let us turn to the persistence property of the model,
defined as follows.

######  Definition 4.4 (Persistence).

If the unique positive solution $I_{t}$, $t\geq 0$ of the stochastic SIS model
obeys

$\limsup\limits_{t\rightarrow\infty}I_{t}\geq\lambda\quad
a.s.\quad\text{and}\quad\liminf\limits_{t\rightarrow\infty}I_{t}\leq\lambda\quad
a.s.,$ (4.14)

where $\lambda$ is a positive constant, namely, $I_{t}$ will rise to or above
$\lambda$ infinitely often with probability one, we say that the infected
population system has the persistence property. Similarly, the numerical
approximation $\\{I_{k}\\}_{k\geq 0}$ is said to have the persistence property
if

$\limsup\limits_{kh\rightarrow\infty}I_{k}\geq\lambda\quad
a.s.\quad\text{and}\quad\liminf\limits_{kh\rightarrow\infty}I_{k}\leq\lambda\quad
a.s..$ (4.15)

We are now in a position to provide some sufficient conditions to ensure the
persistence property of the original stochastic SIS mode (1.3), which have
been established in [9, Theorem 5.1].

######  Theorem 4.5.

If

$\displaystyle R_{0}^{S}=\tfrac{\beta
N}{\mu+\gamma}-\tfrac{\sigma^{2}N^{2}}{2(\mu+\gamma)}>1,$ (4.16)

then for any given initial value $I_{0}\in(0,N)$, the solution of the
stochastic SIS model (1.3) obeys

$\limsup\limits_{t\rightarrow\infty}I_{t}\geq\lambda\quad
a.s.\quad\text{and}\quad\liminf\limits_{t\rightarrow\infty}I_{t}\leq\lambda\quad
a.s.,$ (4.17)

where $\lambda$ is the unique root in $(0,N)$ of

$\displaystyle\beta
N-\mu-\gamma-\beta\lambda-\tfrac{1}{2}\sigma^{2}(N-\lambda)^{2}=0.$ (4.18)

In other words, $I_{t}$ will rise to or above $\lambda$ infinitely often with
probability one.

Next, let us investigate the ability of the numerical approximation produced
by the LCM scheme to preserve the persistence property in the sense that
$I_{k}$ rise to or above a positive constant infinitely often with probability
one.

######  Theorem 4.6.

Let the condition (4.16) hold, i.e., $R_{0}^{S}>1$. For any step-size $h>0$,
if we chose a proper constant $\alpha\in(0,1]$ in the scheme (2.10) satisfying

$\alpha<h^{-\theta}\log\Big{(}\tfrac{\sigma^{2}N}{\sigma^{2}N-\beta+\sqrt{\beta^{2}-2\sigma^{2}(\mu+\gamma)}}\Big{)},$
(4.19)

then for any given initial value $I_{0}\in(0,N)$ the numerical approximation
obeys:

$\limsup\limits_{kh\rightarrow\infty}I_{k}\geq\lambda\quad
a.s.\quad\text{and}\quad\liminf\limits_{kh\rightarrow\infty}I_{k}\leq\lambda\quad
a.s.,$ (4.20)

where $\lambda$ is the unique root of

$\displaystyle\beta
N-\mu-\gamma-\beta\lambda-\tfrac{1}{2}\sigma^{2}(N-\lambda)^{2}=0.$ (4.21)

In other words, the numerical approximations produced by (2.10) will rise to
or above $\lambda$ infinitely often with probability one.

Proof: By the assumption (4.16), there exists a positive constant $\delta$
such that $\delta<R_{0}^{S}-1.$ Recall that

$f(x)=-\tfrac{1}{2}\sigma^{2}\text{e}^{2x}+(\sigma^{2}N-\beta)\text{e}^{x}+(\beta
N-\mu-\gamma-\tfrac{1}{2}\sigma^{2}N^{2}).$

Clearly we have

$\displaystyle\lim\limits_{x\rightarrow-\infty}f(x)$ $\displaystyle=\beta
N-\mu-\gamma-\tfrac{1}{2}\sigma^{2}N^{2}>\delta(\mu+\gamma)>0,$ (4.22)
$\displaystyle\lim\limits_{x\rightarrow(\log N-\alpha h^{\theta})^{-}}f(x)$
$\displaystyle=-\tfrac{1}{2}\sigma^{2}N^{2}\big{(}1-\text{e}^{-\alpha
h^{\theta}}\big{)}^{2}+\beta N\big{(}1-\text{e}^{-\alpha
h^{\theta}}\big{)}-\mu-\gamma.$ (4.23)

It is easy to derive that

$\lim\limits_{x\rightarrow(\log N-\alpha h^{\theta})^{-}}f(x)<0$ (4.24)

holds true, provided

$1-\text{e}^{-\alpha
h^{\theta}}<\tfrac{\beta-\sqrt{\beta^{2}-2\sigma^{2}(\mu+\gamma)}}{\sigma^{2}N}.$
(4.25)

This can be fulfilled under the assumption (4.19), that is,

$\alpha<h^{-\theta}\log\Big{(}\tfrac{\sigma^{2}N}{\sigma^{2}N-\beta+\sqrt{\beta^{2}-2\sigma^{2}(\mu+\gamma)}}\Big{)}.$

Noting

$f^{\prime}(x)=-\sigma^{2}\text{e}^{2x}+(\sigma^{2}N-\beta)\text{e}^{x},$

one can readily infer that $f^{\prime}(x)<0$ and $f(x)$ is decreasing in the
case $\sigma^{2}N-\beta\leq 0$. In the other case $\sigma^{2}N-\beta>0$,

  * •

$f(x)$ is strictly increasing as
$x<\log(\tfrac{\sigma^{2}N-\beta}{\sigma^{2}})$,

  * •

$f(x)$ is strictly decreasing as
$x>\log(\tfrac{\sigma^{2}N-\beta}{\sigma^{2}})$.

Thus the equation $f(x)=0$ has a unique root, denoted as $\eta$, satisfying
$\eta<\log N-\alpha h^{\theta}$.

Now we validate the two assertions in (4.20) in two steps. If the first
assertion of (4.20) is not true, then there exists a positive constant
$\epsilon$ such that

$\displaystyle P(\Omega_{1})>\epsilon,$ (4.26)

where
$\Omega_{1}:=\\{\omega:\limsup\limits_{kh\rightarrow\infty}Y_{k}(\omega)\leq\eta-2\epsilon\\}$.
Hence $\forall\omega\in\Omega_{1}$, there exists a $m=m(\omega)>0$ such that

$\displaystyle\bar{Y}_{k}(\omega)=Y_{k}(\omega)\leq\eta-\epsilon<\log N-\alpha
h^{\theta}$ (4.27)

whenever $k>m$. We can choose $\epsilon$ sufficiently small such that
$f(\eta-\epsilon)<\delta(\mu+\gamma)$, which, together with (4.22), implies
that

$\displaystyle f(Y_{k}(\omega))\geq f(\eta-\epsilon),$ (4.28)

whenever $k>m$. Moreover, by the large number theorem of martingales, there is
an $\Omega_{2}\subset\Omega$ with $P(\Omega_{2})=1$ such that
$\forall\omega\in\Omega_{2}$

$\displaystyle\lim\limits_{kh\rightarrow\infty}\frac{1}{kh}\sum_{i=0}^{k-1}g(Y_{i}(\omega))\Delta
W_{i}(\omega)+\lim\limits_{kh\rightarrow\infty}\frac{1}{kh}\sum_{i=0}^{k-1}g(Y_{i}(\omega))g^{\prime}(Y_{i}(\omega))\Delta\zeta_{i}(\omega)=0.$
(4.29)

Fix any $\omega\in\Omega_{1}\cap\Omega_{2}$ and let $k>m+1$. Then it follows
from (4.28) and by induction that

$\displaystyle\bar{Y}_{k+1}(\omega)$
$\displaystyle=\bar{Y}_{k}(\omega)+f(Y_{k}(\omega))h+g(Y_{k}(\omega))\Delta
W_{k}(\omega)+g(Y_{k}(\omega))g^{\prime}(Y_{k}(\omega))\Delta\zeta_{k}(\omega)$
$\displaystyle=Y_{m+1}(\omega)+\sum_{i=m+1}^{k}f(Y_{i}(\omega))h+\sum_{i=m+1}^{k}\Big{[}g(Y_{i}(\omega))\Delta
W_{i}(\omega)+g(Y_{i}(\omega))g^{\prime}(Y_{i}(\omega))\Delta\zeta_{i}(\omega)\Big{]}$
$\displaystyle\geq Y_{m+1}(\omega)+(k-m)f(\eta-\epsilon)h$
$\displaystyle\quad+\sum_{i=m+1}^{k}\Big{[}g(Y_{i}(\omega))\Delta
W_{i}(\omega)+g(Y_{i}(\omega))g^{\prime}(Y_{i}(\omega))\Delta\zeta_{i}(\omega)\Big{]}.$
(4.30)

Thus

$\displaystyle\liminf\limits_{kh\rightarrow\infty}\frac{\bar{Y}_{k}(\omega)}{kh}\geq
f(\eta-\epsilon)>0,$ (4.31)

and consequently

$\displaystyle\lim\limits_{kh\rightarrow\infty}\bar{Y}_{k}(\omega)=+\infty.$
(4.32)

Therefore

$\lim\limits_{kh\rightarrow\infty}Y_{k}(\omega)=\lim\limits_{kh\rightarrow\infty}\Big{[}\bar{Y}_{k}(\omega)\wedge(\log
N-\alpha h^{\theta})\Big{]}=\log N-\alpha h^{\theta},$ (4.33)

which contradicts (4.27). Thus

$\displaystyle\limsup\limits_{kh\rightarrow\infty}Y_{k}\geq\eta\quad a.s..$
(4.34)

Further, if the second assertion of (4.20) is not true, then there exists a
positive constant $\epsilon$ such that

$\displaystyle P(\Omega_{3})>\epsilon,$ (4.35)

where
$\Omega_{3}:=\\{\omega:\liminf\limits_{kh\rightarrow\infty}Y_{k}(\omega)\geq\eta+2\epsilon\\}$.
Hence $\forall\omega\in\Omega_{3}$, there exists a $m=m(\omega)>0$ such that

$\displaystyle Y_{k}(\omega)\geq\eta+\epsilon$ (4.36)

whenever $k>m$. Simple analysis on $f(x)$ implies that

$\displaystyle f(Y_{k}(\omega))\leq f(\eta+\epsilon)$ (4.37)

whenever $k>m$. Fix any $\omega\in\Omega_{3}\cap\Omega_{2}$. Then it follows
from (4.37) and by induction that

$\displaystyle\bar{Y}_{k+1}(\omega)$
$\displaystyle=Y_{k}(\omega)+f(Y_{k}(\omega))h+g(Y_{k}(\omega))\Delta
W_{k}(\omega)+g(Y_{k}(\omega))g^{\prime}(Y_{k}(\omega))\Delta\zeta_{k}(\omega)$
$\displaystyle\leq
Y_{0}(\omega)+\sum_{i=0}^{m}f(Y_{i}(\omega))h+(k-m)f(\eta+\epsilon)h$
$\displaystyle\quad+\sum_{i=0}^{k}\Big{[}g(Y_{i}(\omega))\Delta
W_{i}(\omega)+g(Y_{i}(\omega))g^{\prime}(Y_{i}(\omega))\Delta\zeta_{i}(\omega)\Big{]}.$
(4.38)

Thus by (4.29) we have

$\displaystyle\limsup\limits_{kh\rightarrow\infty}\frac{\bar{Y}_{k}(\omega)}{kh}\leq
f(\eta+\epsilon)<0,$ (4.39)

and consequently

$\displaystyle\lim\limits_{kh\rightarrow\infty}\bar{Y}_{k}(\omega)=-\infty.$
(4.40)

Therefore we have

$\displaystyle\lim\limits_{kh\rightarrow\infty}Y_{k}(\omega)=\lim\limits_{kh\rightarrow\infty}\Big{[}\bar{Y}_{k}(\omega)\wedge(\log
N-\alpha h^{\theta})\Big{]}=-\infty,$ (4.41)

which contradicts (4.36). Thus

$\displaystyle\liminf\limits_{kh\rightarrow\infty}Y_{k}\leq\eta\quad a.s..$
(4.42)

The assertion (4.20) therefore stands for $\lambda=\text{e}^{\eta}$ and the
proof is completed. ∎

######  Remark 4.7.

Theorem 4.6 illustrates that, when the scheme parameters $\alpha,\theta$ obey
the condition (4.19), the numerical approximation produced by the LCM scheme
with any step-size $h>0$ is able to reproduce the persistence property of the
original model without any addition restriction on the model parameters. It is
worthwhile to point out that, the condition (4.19) can be easily fulfilled, by
taking $\alpha$ to be sufficiently small.

##  5 Numerical experiments

In this section, we provide numerical experiments to illustrate the previous
theoretical findings. Consider the stochastic SIS epidemic model

${\rm d}I_{t}=I_{t}(\beta N-\mu-\gamma-\beta I_{t})\,{\rm d}t+\sigma
I_{t}(N-I_{t})\,{\rm d}W_{t},\quad 0\leq t\leq T$

with a fixed total population $N=10$ and initial value $I_{0}=1$. The
population sizes are measured in units of 1 billion and the unit of time is
assumed to be one day throughout this section.

The approximation errors will be calculated in terms of

$\bigg{(}\mathbb{E}\Big{[}\sup\limits_{k=0,...,M}\big{|}I_{t_{k}}-I_{k}\big{|}^{2}\Big{]}\bigg{)}^{\frac{1}{2}}.$

The newly proposed LCM scheme (2.10) and the Lamperti Euler-Maruyama (LEM)
scheme proposed in [28] will be both tested for comparison. Moreover, we use
numerical approximations with the step-size $h_{exact}=2^{-14}$ to represent
the analytic solutions of the original model and those with step-sizes
$h=2^{-i},i=4,5,6,7,8$ for numerical schemes, over $10^{4}$ Brownian paths.

The following two sets of parameters are taken for tests.

######  Example 5.1.

$T=1,\beta=0.5,\sigma=0.2,\mu+\gamma=4,\alpha=0.1,\theta=2$.

######  Example 5.2.

$T=1,\beta=0.7,\sigma=0.1,\mu+\gamma=2,\alpha=1,\theta=3$.

Figure 1:  Example 5.1

Figure 2:  Example 5.2

Table 1: Sample average errors |  $h$ |  $2^{-10}$ |  $2^{-9}$ |  $2^{-8}$ |  $2^{-7}$ |  $2^{-6}$
---|---|---|---|---|---|---
Ex 5.1 |  LCM |  0.0006 |  0.0011 |  0.0023 |  0.0045 |  0.0094
|  LEM |  0.0015 |  0.0030 |  0.0061 |  0.0123 |  0.0250
Ex 5.2 |  LCM |  0.0025 |  0.0052 |  0.0105 |  0.0211 |  0.0425
|  LEM |  0.0014 |  0.0028 |  0.0058 |  0.0116 |  0.0235
Table 2:  Least-squares fit for the convergence rate $q$ |  LCM scheme |  LEM scheme
---|---|---
Ex 5.1 |  $q=1.0042$, resid = 0.0359 |  $q=1.0240$, resid = 0.0132
Ex 5.2 |  $q=1.0209$, resid = 0.0154 |  $q=1.0240$, resid=0.0125

Figs 2 and 2 show the average sample errors against various step-sizes on a
log-log scale. Detailed data is presented in Table 1. The approximation errors
of the LCM and LEM scheme are plotted in the blue and red solid lines. The
black dashed lines are reference lines of slope 1. From the figures we can
easily identify convergence rates of order 1 for both the LCM and LEM scheme.
Moreover, by assuming that $error\approx Ch^{q}$ such that
$\log\\{error\\}\approx\log C+q\log h$, the convergence rate $q$ and the least
square residual is obtained with a least-squares fitting as presented in Table
2. These results confirm the expected convergence rate.

As for the dynamic behaviors of the approximations, we test three different
schemes for comparison: the newly proposed LCM scheme, the LEM scheme and the
usual Euler-Maruyama scheme directly applied to the original SDE without any
transformation. The following four sets of model parameters, carefully chosen
to meet all conditions in Theorems 4.3 or 4.6, with fixed $\alpha=0.1$ and
$\theta=2$, are taken as examples.

######  Example 5.3.

$\beta=0.2,\sigma=0.2,\mu+\gamma=2,h=2^{-3}$.

######  Example 5.4.

$\beta=1.2,\sigma=0.8,\mu+\gamma=6,h=2^{-3}$.

######  Example 5.5.

$\beta=0.8,\sigma=0.2,\mu+\gamma=2,h=2^{-3}$.

######  Example 5.6.

$\beta=3,\sigma=0.6,\mu+\gamma=6,h=2^{-i},i=5,4,3,2,1,0$.

We mention that parameters in Example 5.3 and 5.4 satisfy the conditions of
extinction in Theorem 4.2 and 4.3, while those in Example 5.5 and 5.6 meet the
requirements of persistence in Theorem 4.5 and 4.6. Moreover, we can conclude
by simple calculations that

$\liminf\limits_{kh\rightarrow\infty}I_{k}\leq
7.3205\leq\limsup\limits_{kh\rightarrow\infty}I_{k}\quad a.s.$

when $\beta=0.8$, $\sigma=0.2$ and $\mu+\gamma=2$, and

$\liminf\limits_{kh\rightarrow\infty}I_{k}\leq
7.6759\leq\limsup\limits_{kh\rightarrow\infty}I_{k}\quad a.s.$

when $\beta=3$, $\sigma=0.6$ and $\mu+\gamma=6$.

Figure 3: Example 5.3: $\beta=0.2,\sigma=0.2,\mu+\gamma=2,h=2^{-3}$ Figure 4:
Example 5.4: $\beta=1.2,\sigma=0.8,\mu+\gamma=6,h=2^{-3}$ Figure 5:  Example
5.5: $\beta=0.8,\sigma=0.2,\mu+\gamma=2,h=2^{-3}$,

Figure 6: Example 5.6: $\beta=3,\sigma=0.6,\mu+\gamma=6,h=2^{-5}$

Figure 7: Example 5.6: $\beta=3,\sigma=0.6,\mu+\gamma=6,h=2^{-4}$

Figure 8: Example 5.6: $\beta=3,\sigma=0.6,\mu+\gamma=6,h=2^{-3}$

Figure 9: Example 5.6: $\beta=3,\sigma=0.6,\mu+\gamma=6,h=2^{-2}$

Figure 10: Example 5.6: $\beta=3,\sigma=0.6,\mu+\gamma=6,h=2^{-1}$

Figure 11: Example 5.6: $\beta=3,\sigma=0.6,\mu+\gamma=6,h=1$

Figures 3-7 plot one-path simulations for the three different schemes. We can
clearly observe from Fig. 3 and 4 that both the LCM and LEM schemes using the
step-size $h=2^{-3}$ preserve the extinction property of the original SIS
model for these two examples. On the contrary, the Euler-Maruyama scheme turns
out to be not domain-preserving for Example 5.3, and even provide explosive
results for Example 5.4. From Fig. 5 one can observe that the three considered
schemes applied to Example 5.5 using the step-size $h=2^{-3}$ well preserve
the persistence property of the model. However, as indicated in Figures 7-11,
for Example 5.6 the newly proposed LCM scheme shows a persistent property for
all step-sizes $h=2^{-i},i=0,1,...,5$, even with large step-sizes $h=1/2,1$.
Instead, the LEM scheme for Example 5.6 can only preserve the persistence
property for sufficiently small step-sizes $h=2^{-5},2^{-4}$ but shows an
extinct property for step-sizes $h=2^{-i},i=3,2,1,0$. The Euler-Maruyama
scheme even explodes for all the chosen step-sizes. These results confirm the
best dynamics-preserving ability of the LCM scheme.

##  6 Conclusion

In this paper we propose a first-order strongly convergent scheme for the
stochastic SIS epidemic model which preserves the domain as well as the
dynamic behaviors of the model unconditionally. Based on a logarithm
transformation, the numerical scheme is designed to be explicit and easily
implementable. The strong convergence rate of the scheme is carefully analyzed
and proved to be order $1$. Moreover, the proposed scheme is able to reproduce
the extinction and persistence properties of the original model for any
stepsize $h>0$. Numerical experiments are provided to verify the theoretical
results. In the future we plan to work on the stationary distribution for the
numerical approximations, and design higher-order schemes that are also able
to reproduce the dynamic properties of the considered model as stated in [17].

## References

  * [1] A. Alfonsi.  Strong order one convergence of a drift implicit Euler scheme: Application to the CIR process.  Statistics & Probability Letters, 83(2):602–607, 2013.
  * [2] A. Andersson and R. Kruse.  Mean-square convergence of the BDF2-Maruyama and backward Euler schemes for SDE satisfying a global monotonicity condition.  BIT Numerical Mathematics, 57(1):21–53, 2017.
  * [3] W.-J. Beyn, E. Isaak, and R. Kruse.  Stochastic C-stability and B-consistency of explicit and implicit Euler-type schemes.  Journal of Scientific Computing, 67(3):955–987, 2016.
  * [4] W.-J. Beyn, E. Isaak, and R. Kruse.  Stochastic C-stability and B-consistency of explicit and implicit Milstein-type schemes.  Journal of Scientific Computing, 70(3):1042–1077, 2017.
  * [5] C.-E. Bréhier.  Approximation of the invariant distribution for a class of ergodic SDEs with one-sided Lipschitz continuous drift coefficient using an explicit tamed Euler scheme.  arXiv preprint arXiv:2010.00512, 2020.
  * [6] Y. Cai, J. Hu, and X. Mao.  Positivity and boundedness preserving numerical scheme for the stochastic epidemic model with square-root diffusion term.  Applied Numerical Mathematics, 182:100–116, 2022.
  * [7] L. Chen, S. Gan, and X. Wang.  First order strong convergence of an explicit scheme for the stochastic SIS epidemic model.  Journal of Computational and Applied Mathematics, 392:113482, 2021\.
  * [8] M. B. Giles.  Multilevel monte carlo path simulation.  Operations research, 56(3):607–617, 2008.
  * [9] A. Gray, D. Greenhalgh, L. Hu, X. Mao, and J. Pan.  A stochastic differential equation SIS epidemic model.  SIAM Journal on Applied Mathematics, 71(3):876–902, 2011.
  * [10] D. J. Higham, X. Mao, and A. M. Stuart.  Strong convergence of Euler-type methods for nonlinear stochastic differential equations.  SIAM Journal on Numerical Analysis, 40(3):1041–1063, 2002.
  * [11] M. Hutzenthaler and A. Jentzen.  Convergence of the stochastic Euler scheme for locally Lipschitz coefficients.  Foundations of Computational Mathematics, 11(6):657–706, 2011.
  * [12] M. Hutzenthaler and A. Jentzen.  On a perturbation theory and on strong convergence rates for stochastic ordinary and partial differential equations with nonglobally monotone coefficients.  The Annals of Probability, 48(1):53–93, 2020.
  * [13] M. Hutzenthaler, A. Jentzen, and P. E. Kloeden.  Strong and weak divergence in finite time of Euler’s method for stochastic differential equations with non-globally Lipschitz continuous coefficients.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 467(2130):1563–1576, 2011.
  * [14] M. Hutzenthaler, A. Jentzen, and P. E. Kloeden.  Strong convergence of an explicit numerical method for SDEs with nonglobally Lipschitz continuous coefficients.  The Annals of Applied Probability, 22(4):1611–1641, 2012.
  * [15] W. O. Kermack and A. G. Mckendrick.  A contribution to the mathematical theory of epidemics.  Proceedings of the Royal Society of London Series A, 115, 1927.
  * [16] Z. Lei, S. Gan, and Z. Chen.  Strong and weak convergence rates of logarithmic transformed truncated EM methods for SDEs with positive solutions.  Journal of Computational and Applied Mathematics, 419:114758, 2023\.
  * [17] R. Liu and X. Wang.  A higher order positivity preserving scheme for the strong approximations of a stochastic epidemic model.  Communications in Nonlinear Science and Numerical Simulation, page 107258, 2023.
  * [18] X. Mao.  Stochastic Differential Equations and Applications.  Elsevier, 2007.
  * [19] X. Mao.  The truncated Euler–Maruyama method for stochastic differential equations.  Journal of Computational and Applied Mathematics, 290:370–384, 2015\.
  * [20] W. Meiske.  Hethcote, H. W., J. A. Yorke: Gonorrhea. Transmission Dynamics and Control. Springer‐Verlag, Berlin — Heidelberg — New York — Tokyo 1984 (Lecture Notes in Biomathematics, vol. 66), 105 S., DM 25.  Biometrical Journal, 28, 2010.
  * [21] A. Neuenkirch and L. Szpruch.  First order strong approximations of scalar SDEs defined in a domain.  Numerische Mathematik, 128(1):103–136, 2014.
  * [22] S. Sabanis.  A note on tamed Euler approximations.  Electronic Communications in Probability, 18:1–10, 2013.
  * [23] S. Sabanis and Y. Zhang.  On explicit order 1.5 approximations with varying coefficients: the case of super-linear diffusion coefficients.  Journal of Complexity, 50:84–115, 2019.
  * [24] M. V. Tretyakov and Z. Zhang.  A fundamental mean-square convergence theorem for SDEs with locally Lipschitz coefficients and its applications.  SIAM Journal on Numerical Analysis, 51(6):3135–3162, 2013.
  * [25] X. Wang and S. Gan.  The tamed Milstein method for commutative stochastic differential equations with non-globally Lipschitz continuous coefficients.  Journal of Difference Equations and Applications, 19(3):466–490, 2013.
  * [26] X. Wang, J. Wu, and B. Dong.  Mean-square convergence rates of stochastic theta methods for SDEs under a coupled monotonicity condition.  BIT Numerical Mathematics, 60(3):759–790, 2020.
  * [27] H. Yang.  First order strong convergence and extinction of positivity preserving logarithmic truncated Euler-Maruyama method for the stochastic SIS epidemic model.  arXiv preprint arXiv:2111.00771, 2021.
  * [28] H. Yang and J. Huang.  First order strong convergence of positivity preserving logarithmic Euler–Maruyama method for the stochastic SIS epidemic model.  Applied Mathematics Letters, 121:107451, 2021.
  * [29] H. Yang, Y. Pan, W. Liu, and Z. Mu.  Numerical analysis of split-step $\theta$ methods with truncated Wiener process for a stochastic SIS epidemic model.  Journal of Computational and Applied Mathematics, 415:114433, 2022\.
  * [30] X. Yang, M. Li, Z. Yang, and C. Zhang.  Numerical analysis of a linearly backward Euler method with truncated Wiener process for a stochastic SIS model.  Numerical Algorithms, 93(2):563–579, 2023.
  * [31] Y. Yi, Y. Hu, and J. Zhao.  Positivity preserving logarithmic Euler-Maruyama type scheme for stochastic differential equations.  Communications in Nonlinear Science and Numerical Simulation, 101:105895, 2021.