Stability of stochastic SIS model with disease deaths and variable diffusion rates

The SIS model is a fundamental model that helps to understand the spread of an infectious disease, in which infected individuals recover without immunity. Because of the random nature of infectious diseases, we can estimate the spread of a disease in population by stochastic models. In this article, we present a class of stochastic SIS model with births and deaths, obtained by superimposing Wiener processes (white noises) on contact and recovery rates and allowing variable diffusion rates. We prove existence of the unique, positive and bounded solution of this nonlinear system of stochastic differential equations (SDEs) and examine stochastic asymptotic stability of equilibria. In addition, we simulate the model by considering a numerical approximation based on a balanced implicit method (BIM) on an appropriately bounded domain D ⊂ R2.


Introduction
Each year infectious diseases kill almost 9 million people, many of them children under five, and they also cause enormous burdens through life-long disability [37]. Because of ethical concerns and the nature of diseases, it is difficult to do experiments searching an effective strategy for the management of diseases. Mathematical models may be needed. Mathematical epidemiology studies the spread of diseases in populations by using tools from mathematics, statistics, and computer science. SIS and SIR (Susceptible, Infected, Removed) are building blocks of modern mathematical epidemiology. If the infectives recover from the disease without immunity, then they immediately become susceptible. Such a model is called an SIS model. SIS models are appropriate for most diseases transmitted by bacterial or helminth agents [7].
Since we repeatedly use some terms from epidemiology, we briefly define them here. An outbreak of a disease that spreads rapidly and widely is called epidemic, and a disease that exists permanently in a particular location is called endemic. One of the most important quantities in epidemiology is a basic reproduction number, the expected number of secondary infections produced when one infected individual entered a fully susceptible population [14]. It determines whether there is an epidemic or not.
Generally, epidemic models admit two types of equilibria; disease-free and endemic. If the disease-free equilibrium is globally asymptotically stable then the disease dies out. If the endemic equilibrium is globally asymptotically stable then the disease persists in the population at the equilibrium level.
Lyapunov's direct method [26] is a useful tool to establish the global stability of equilibria of an epidemic model and has been used in [5, 10, 13, 15, 19-25, 35, 36]. A historic function V(x 1 , . . . , ) is considered as a good candidate for the Lyapunov function for epidemic models [10].
We consider a SIS model with births and deaths of the form S (t) = − βS(t)I(t) + µ K − S(t) + αI(t) (1.1) In this model, the non-negative constant µ represents the per capita birth rate. So µK is the number of births (immigrants) where K is a carrying capacity (i.e. maximum total population size). Here we consider all the newborns are susceptible. The model (1.1) assumes that the birth rate µ is equal to the natural death rate. The contact rate β is the average number of contacts per infective. βSI represents the number of new infections in unit time. 1/α, the reciprocal of the recovery rate α, is the mean infected period. Infected individuals spread the disease to susceptibles, and remain in the infected class in that time, before recovered from the disease without immunity. That is, individuals become susceptible immediately once they have recovered. Non-negative constant γ is the disease related death rate such that 1/γ is a mean of the disease related death period.
Since N(t), the affected population size at time t, is N(t) = S(t) + I(t), we obtain by adding the above equations (1.1). Therefore the affected population size is not constant and may vary in time. In what follows, we shall introduce a stochastic version of model (1.1) preserving the same equation for affected population size N(t), but now under the presence of martingale-type noises.

Stochastic model with martingale-type noise
In this paper we present a family of stochastic SIS model with births and deaths, obtained by superimposing white noises on contact and recovery rates. This is only one of several ways of introducing random noise into model [2,4,16,33]. We consider that contact and recovery rates are subject to random disturbances. We perturb the deterministic system (1.1) by white noises and dW 2 (t) dt ; and obtain a stochastic model by replacing β and α by β + F 1 S(t), I(t) dW 1 (t) dt and α + F 2 S(t), I(t) dW 2 (t) dt respectively, where the perturbation functions F 1 and F 2 are locally Lipschitz-continuous (almost surely w.r.t. 2D Lebesgue measure λ 2 , i.e. except for a λ 2 -null subset of D) on W 1 and W 2 are i.i.d. Wiener processes defined on a complete filtered probability space (Ω, F , {F t } t≥0 , P) for all t ≥ 0. Therefore we obtain a stochastic model (interpreted in Itô sense) where the parameters α, β, γ and µ are non-negative constants. We ignore t in the above SDE and express the stochastic SIS model with disease deaths in the form (2.2) Similar stochastic models with specific diffusion terms are discussed in [12,16]. Since the perturbation functions F 1 and F 2 are arbitrary, our model represents a rather non-parametric approach to the class of stochastic SIS models with respect to the diffusion terms, in contrast to the models in literature.
Gray et al. [12] investigated the properties of a stochastic SIS model for constant populations in the form Obviously, this represents a subclass of models (2.2) with constant diffusion rates F 2 = 0 and F 1 = σ. In [12] they showed the positivity of solutions and established conditions for extinction and persistence of infective population. Imhof and Walcher [16] studied stochastic single-substrate chemostat model by considering biomass concentrations of two microbe species. They showed that, under certain conditions, the stochastic model leads to extinction even though the deterministic counterpart predicts persistence.
In this paper, we first prove the existence of unique strong solution of the stochastic SIS model (2.2) on D. Then we discuss stochastic asymptotic stability of disease free and endemic equilibria. We show that the disease free equilibrium (S 1 , I 1 ) = (K, 0) is globally stochastically asymptotically stable under the condition of βK ≤ α. Furthermore, we prove that the endemic equilibrium (S 2 , Generally speaking, if F i (S 2 , I 2 ) = 0 for the diffusion rates F i then the related stochastic system will not have a positive steady state solution (endemic equilibrium). In here, we consider random fluctuations around endemic equilibrium by the assumption F i (S 2 , I 2 ) = 0, i.e. no perturbation at the equilibrium point. Stochastic perturbations which are proportional to the deviation of the system state from the endemic equilibrium have been first considered by Beretta et al. [6], Carletti [9] and Lahrouz et al. [25]. Finally, we demonstrate the applicability of the mathematical approach with simulations and show parametric dependence of asymptotic stability of related equilibria in view of expectations and variances. An appendix explains the mean square convergence of the class of balanced implicit numerical methods which we use to generate the results of positive and D-invariant simulations for our examples.

Preliminary
Consider a d-dimensional Itô stochastic differential equation of the form as the minimally generated σ-algebra generated by W and X 0 is an R d -valued random variable.
The infinitesimal generator L associated with Itô SDE (3.1) is given by V(x, t) → ∞ as n → ∞ then, for any X 0 which is independent of σ-algebra σ(W) such that the initial condition P(X 0 ∈ D) = 1 is met, there is a unique Markovian, global, continuous time solution X of (3.1) with X(t 0 ) = X 0 , and X(t) ∈ D for all t ∈ [t 0 , T] (a.s.).

Existence of a unique global solution of stochastic SIS model
In the statements below and its proof, consider the function where we refer to the essential supremum here (according to Lebesgue's measure theory) and [·] + denotes the positive part of inscribed expression. Then, the stochastic SIS model (2.1) has a unique, continuous time, Markovian global solution process S(t), I(t) t≥t 0 on t ≥ t 0 and this solution is almost surely invariant with respect to D, i.e. for all initial data (S 0 , I 0 ) ∈ D independent of σ(W) P(∀t ≥ t 0 : (S(t), I(t)) ∈ D) = 1.
Proof. We use stochastic invariance theorem stated by Khas'minskiȋ [18] and follow the ideas in [29]. For arbitrary terminal time T > t 0 , consider the stochastic process X defined by Since the coefficients of the system (2.1) are a.s. locally Lipschitz-continuous and satisfy linear growth condition on D, for any initial value (S 0 , I 0 ) ∈ D, there is a unique continuous time, Markovian local solution on t ∈ t 0 , τ T (D) , where τ T (D) is the random time of first exit of stochastic process S(t), I(t) t 0 ≤t≤t≤T from the interior of domain D before or at T, started in S(t 0 ), I(t 0 )) ∈ D at the initial time t 0 ∈ R. We define τ T (D) = ∞ if X does not hit the boundary of D before T or at T. To make this solution X global, we shall prove that P τ T (D) = ∞ = 1 (a.s.) for all finite terminal times T. For fixed initial values (S 0 , I 0 ) inside D, we consider the fully extended system (which is equivalent to original SIS model (2.2)) for n ∈ N. Obviously, for any initial value (S 0 , I 0 ) inside D n , the system (2.1) has a unique solution up to stopping time τ T (D n ). Define where 0 < S + I < K on D and suppose that EV(S 0 , I 0 ) < ∞. Note that 0 < V(S, I) < +∞ for (S, I) ∈ D. Furthermore, let the essential supremum be where [·] + is the positive part of inscribed expression. Obviously, under hypothesis (3.5) with almost surely continuous F i , the constant c is positive and finite. Then, the infinitesimal generator of process X = (S(t), I(t)) t 0 ≤t≤T is of the form Recall that all parameters α, β, γ, and µ are non-negative. Evaluating the aforementioned expression for LV yields that ∀ (S, where the finite c ≥ 0 is defined as in (3.7) and [·] + is the positive part of inscribed expression. Therefore LV(S, I) ≤ c on D.
Note that, by elementary calculus, we find that inf (S,I)∈∂D n V(S, I) > n + 2 ln (2) since ∂D n consists of the boundary of triangle bounded by the lines I = Ke −n , S = Ke −n , and . Now, define τ n := min{t, T, τ T (D n )} and apply Dynkin's formula to get to Therefore, since D n ⊂ D for all n ∈ N, we have for all (S 0 , I 0 ) ∈ D n (for large n), and for all fixed t ∈ [t 0 , T).
This proves the invariance property and the existence of the strong solution S(t), I(t) t≥t 0 on D, provided that (S 0 , I 0 ) ∈ D. Remark 3.3. It remains to discuss the boundary cases for initial data. Recall that I = 0 and S = 0 are not inside of the domain D . We study these cases separately.
is left with positive values of I(t) ≤ K. Note that, under (3.5) and due to continuity of F 2 (0, ·), we find that the term IF 2 (0, Moreover, if the initial condition is given by Basically, for validity of (3.5), it would also suffice to take any a.s. continuous functions F 2 on D which satisfy the limit condition lim sup A further remarkable fact is that, for existence of unique, global, continuous time Markovian solutions (S, I), it basicly suffices to impose no significant assumption on the choice of diffusion rate F 1 other than almost sure Lipschitz-continuity of Carathéodory function F 1 on compact set D (i.e. in 2D Lebesgue's almost sure sense). However, the choice of diffusion rates F k 's plays an important role for asymptotic stability of equilibria as seen below.

Stability of equilibria 4.1 Preliminary
Consider a d-dimensional SDE Assume that the functions f and g satisfy, in addition to the existence and uniqueness assumptions, f (x * , t) = 0 and g(x * , t) = 0 for equilibrium solution x * for t ≥ t 0 . Let D ⊆ R d be nonempty and simple-connected.
where X s,x 0 (t) denotes the solution of (4.1) at time t ≥ s, satisfying X(s) = x 0 ∈ D at initial time s.

Definition 4.2.
The equilibrium solution x * ∈ D of the SDE (4.1) is said to be stochastically asymptotically stable iff it is stochastically stable and at every s we have Remark 4.5. Since the equations of our stochastic SIS model do not have time-dependent coefficients, the requirement of "for every s" in above definitions can be dropped.
Theorem 4.6. Assume that f and g satisfy the existence and uniqueness assumptions and they have continuous coefficients with respect to t. Let x 0 be constant with probability 1 and P(∀t ≥ t 0 : Then, the equilibrium solution x * of (4.1) is stochastically stable.
ii) If, in addition, V is decrescent (there exists a positive definite function V 1 such that V(x, t) ≤ V 1 (x) for all x ∈ U h ) and LV(x, t) is negative definite, then the equilibrium solution x * is stochastically asymptotically stable.
iii) If assumptions ii) hold for all x ∈ D for a radially unbounded function V ∈ C 2,1 R d × [t 0 , ∞) defined everywhere then the equilibrium solution x * is globally stochastically asymptotically stable.
Remark 4.7. This Theorem 4.6 is a slight modification of the famous stability theorem on R d due to L. Arnold [3], and its proof is very similar to that in [3], hence we may omit its proof here.

Global stochastic asymptotic stability of disease free equilibrium
Recall that, for our stochastic SIS model, in previous section we have justified to confine our related domains D ⊂ R 2 + satisfying where its entire dynamic takes place (a.s.). For the rest of this paper, we assume that D satisfies that compact triangular form in R 2 .  This theorem also says that an outbreak can be controlled by decreasing the number β of contacts and infection period 1/α. It is remarkable that the disease related death rate γ, birth rate µ (= natural death rate) and the perturbation functions F 1 and F 2 do not play a role in the sufficient condition for the stochastic stability of the disease free equilibrium. Remark 4.9. By writing the stability condition, βK ≤ α, in terms of the basic reproduction number R 0 = βK α+γ+µ of the deterministic model, we obtain So our sufficient condition for stability of the disease free equilibrium does not contradict with the well-known R 0 < 1 rule.
then the disease free equilibrium (K, 0) is globally stochastically asymptotically stable if βK α ≤ 1. Example 4.12. In order to illustrate an application of Theorem 4.8, we simulate the solution of the stochastic SIS model (2.1) by considering a numerical approximation of process (S, I) based on Balanced Implicit Methods (BIMs) [27,31]. We use the following scheme (n ∈ N) (as proposed in [34] based on [29]) S n+1 = S n + − βS n I n + µ(K − S n ) + αI n ∆ n − S n I n F 1 (S n , I n )∆W 1 n + I n F 2 (S n , I n )∆W 2 n + A n (S n − S n+1 ) I n+1 = I n + βS n I n − (α + γ + µ)I n ∆ n + S n I n F 1 (S n , I n )∆W 1 n − I n F 2 (S n , I n )∆W 2 n + A n (I n − I n+1 ) (4.10) for a discretization of (2.1) along partitions 0 = t 0 < t 1 < · · · < t n < t n+1 < · · · with step sizes ∆ n = t n+1 − t n , where (for S n > 0) and W j 's are standard Wiener processes defined on a (complete) filtered probability space (Ω, F , {F t } t≥0 , P) and which are independent of the initial value (S 0 , I 0 ) ∈ D with finite second moments E (S 0 , I 0 ) 2 < ∞. Here, as in all of our simulations, the independent Gaussian N (0, ∆ n )-distributed increments ∆W k n are generated by the Polar Marsaglia method. The algorithm for our simulations is programmed in C++. The weights A n are carefully chosen such that we achieve positivity and invariance of BIMs (4.10) on bounded domain D whenever we start at (S 0 , I 0 ) ∈ D. For simplicity, all our figures below are generated by equidistant step sizes ∆.
Based on the pioneering works of [27], [29] and [31], a more detailed study on the convergence of this numerical method (4.10) along with positivity and stability of the numerical solution has been carried out in [34] (cf. remarks in appendix). Thus, in this article, we may leave out further quantitative analysis of related numerical algorithm because we focus more on qualitative analysis of our SIS model.
We consider the stochastic SIS model    Simulations agree with analytic study βK α = 0.31 < 1 for different perturbation functions since trajectories of solutions of stochastic system (4.11) stay nearby the disease free equilibrium when time goes to infinity.
Since R 0 = 1.92, the R 0 > 1 condition satisfies. The inequality (4.13) also satisfies because is negative definite for all (S, I) ∈ D = (S, I) : S > 0, I > 0, S + I ≤ K . Hence Theorem 4.14 guaranties the stochastic asymptotic stability of the endemic equilibrium (S 2 , I 2 ) = (260.13, 0.24). In Figure 4.4, we fix all the parameters other than the recovery rate α and plot the graphs of expected values of Susceptible and Infected versus time versus α. These pictures show the effects of the recovery rate on the asymptotic stability of equilibria. If α gets large then we lose existence of an endemic equilibrium and verify the stochastic asymptotic stability of the disease free equilibrium of the system. Similarly in Figure 4.5, we fix all the parameters other than the contact rate β and plot the graphs of expected values of Susceptible and Infected versus time versus β. These pictures show the effects of the contact rate on the asymptotic stability of equilibria. If β gets small then we lose existence of an endemic equilibrium and verify the stochastic asymptotic stability of the disease free equilibrium.    In order to investigate stochastic effects for the classical model (4.15), we put noises into the contact and recovery rates. We replace β by β + S−S 2 is not negative for all (S, I) ∈ D (for instance LV(20, 1) = 38 > 0). In that case the stochastic noise destabilizes a stable system.

Conclusion and summary
We obtained a family of stochastic SIS models with generally state-dependent diffusion coefficients F i by introducing random fluctuations on contact and recovery rates in well-known deterministic model. We extended some major results of paper [29] on stochastic logistic equations in R 1 to R 2 here. The existence of unique strong solution of IVPs for non-linear SDEs has been shown by a mathematically rigorous proof of Theorem 3.2. This is a non-trivial task since we deal with non-linear SDEs with non-globally Lipschitz-continuous coefficients. This investigation revealed a biologically relevant bounded domain D ⊂ R 2 + on which random dynamics of Susceptible and Infected with non-globally Lipschitz-continuous coefficients takes places. This fact is non-trivial under the presence of erratic unbounded martingale-type noises. Furthermore we discussed stochastic asymptotic stability of disease free and endemic equilibria. As common, stochastic asymptotic stability of equilibria is connected to the basic reproduction number R 0 . Stochastic asymptotic stability of the disease free equilibrium depending on the recovery rate α and the contact rate β is shown by Theorem 4.8 for almost sure locally Lipschitz-continuous perturbation functions F 1 and F 2 of Carathéodory-type under the condition βK ≤ α. This condition does not contradict with the well-known R 0 < 1 rule. This is verified by the help of invariance principle and Lyapunov's second method. Finally the stochastic asymptotic stability of the endemic equilibrium is investigated by Theorem 4.14.
A sufficient condition for stochastic asymptotic stability is found in terms of parameters and the perturbation functions (4.13). A remarkable fact of the criteria (4.13) is that a sufficient condition for stability can be found even for arbitrary locally Lipschitz-continuous functions F 1 and F 2 which are Lipschitz-continuous on D. Some graphical illustrations demonstrate the applicability of the mathematical approach. The simulations show parametric dependence of asymptotic stability of related equilibria in view of expectations and variances. Moreover, we have exploited our main idea of studying effects of variable diffusion rates F k on the qualitative behavior of stochastic models -an idea which is originally introduced by the authors in [28] (but there for the case of stochastic SIR model).

Appendix A Mean square convergence of well-posed BIMs
Mean square convergence of numerical sequences Y n := (S n , I n ) to exact solution X(t) := (S(t), I(t)) of our stochastic SIS models (2.2) can be verified along any nonrandom partitions (η = η ∆ ([0, T])) ∆>0 of fixed, finite time-intervals [0, T] when maximum step size ∆ = max ∆ i = max {|t i+1 − t i | : t i , t i+1 ∈ η} tends to zero. The criterion of mean square convergence is given by where γ is said to be the (least) order (rate) of mean square convergence of numerical sequence Y = (Y i ) i∈N and K = K(T) the leading error constant which depends on several system-parameters (here · is the Euclidean vector norm in R 2 ). Then, the numerical approximation (Y n ) n∈N with Y n = (S n , I n ) governed by S n+1 = S n + − βS n I n + µ(K − S n ) + αI n ∆ n − S n I n F 1 (S n , I n )∆W 1 n + I n F 2 (S n , I n )∆W 2 n + A n (S n − S n+1 ) I n+1 = I n + βS n I n − (α + γ + µ)I n ∆ n + S n I n F 1 (S n , I n )∆W 1 n − I n F 2 (S n , I n )∆W 2 n + A n (I n − I n+1 ) with weights (∆W k n are Gaussian N (0, ∆ n )-distributed) A n = (α + γ + µ + βI n ) ∆ n + K |F 1 (S n , I n ) ∆W 1 n | + K S n |F 2 (S n , I n ) ∆W 2 n | is mean square convergent with order γ = 0.5 towards the exact solution X = (S, I) of (2.2) on D, whenever started at Y 0 = X(0) = (S(0), I(0)) ∈ D (a.s.).
For detailed proof of Theorem A.1, one may combine the results of papers [31] and [32] with Lyapunov function V(S, I) = 1 + S 2 + I 2 . Note that the BIMs with specific weights A n are positive and D-invariant almost surely (for detailed verification of those facts, see also [34]). For example, every polynomial F k (S, I) in the variables S and I guarantees the Lipschitz-continuity of terms SIF 1 (S, I) and IF 2 (S, I) on compact domain D.
In passing, we note that the numerical scheme (A.2) can be equivalently rewritten to its increment form ∆S n = − βS n I n + µ(K − S n ) + αI n ∆ n − S n I n F 1 (S n , I n )∆W 1 n + I n F 2 (S n , I n )∆W 2 n 1 + A n ∆I n = βS n I n − (α + γ + µ)I n ∆ n + S n I n F 1 (S n , I n )∆W 1 n − I n F 2 (S n , I n )∆W 2 n 1 + A n in terms of their increments ∆S n = S n+1 − S n and ∆I n = I n+1 − I n . This is obviously possible due to the linear-implicit character of BIMs (A.2). Moreover, the scheme values (S n , I n ) and A n cannot explode to infinity whenever S 0 > 0 due to its nonnegativity and Lipschitz-continuous coefficients F k . In fact, when S n tends to 0 then A n with any coefficient |F 2 (0, I n )| > 0 would tend to +∞ and ∆S n = S n+1 − S n to 0 at the same time. So, we may define the increment ∆S n+1 = 0 whenever S n = 0. However, although of less importance, this case is physically relevant (from the point of view of any practically relevant biological modelling) since any disease cannot spread out more if there are no susceptable species. Furthermore, the positivity of both S and I, and the boundedness S + I ≤ K of all numerical values (S n , I n ) of BIMs (A.2) are guaranteed by the choice of any step sizes ∆ n satisfying (α + γ + µ)∆ n ≤ 1, βK∆ n ≤ 1 for all n ∈ N. Consequently, the numerical scheme of BIMs (A.2) with those restrictions on step sizes ∆ n is well-posed for any initial data (S 0 , I 0 ) ∈ D.