A stochastic differential equation SIS epidemic model with two correlated Brownian motions

In this paper, we introduce two perturbations in the classical deterministic susceptible-infected-susceptible epidemic model with two correlated Brownian motions. We consider two perturbations in the deterministic SIS model and formulate the original model as a stochastic differential equation with two correlated Brownian motions for the number of infected population, based on the previous work from Gray et al. (SIAM J Appl Math 71(3):876–902, 2011) and Hening and Nguyen (J Math Biol 77:135–163, 2017. https://doi.org/10.1007/s00285-017-1192-8). Conditions for the solution to become extinction and persistence are then stated, followed by computer simulation to illustrate the results.

infected people after recovery, while the classical susceptible-infected-susceptible (SIS) models (1) can be used to explain disease transmission.
with initial value I 0 ∈ (0, N ). Specifically, S(t) and I (t) are the numbers of susceptibles and infected at time t. N is the total size of the population where the disease is found. μ is the per capita death rate, and γ represents the rate of infected individuals becoming cured, while β is the per capita disease contact rate. Moreover, environmental noises, such as white noise and telegraph noise, are taken into consideration in deterministic models to help us understand dynamic behaviours in epidemic models. There are many examples studying the behaviour of both deterministic [1,2] and stochastic [3][4][5][6] SIS epidemic models. Different medical means on controlling the disease are also mathematically applied in SIS model such as [7][8][9]. Also, Gray et al. [10] apply a perturbation on the disease transmission coefficient in SIS model.
However, to the best of our knowledge, there is not enough work on incorporating white noise with μ + γ in the SIS epidemic model (1). Here, we suppose that the variance of estimating μ + γ is proportional to the number of susceptible population. Consequently, we then add another perturbation on per capita death rate and infectious period μ + γ .

(t)
− σ 2 I (t) N − I (t) dB 4 (t) (2) with initial value I (0) = I 0 ∈ (0, N ), B 3 and B 4 are two independent Brownian motions. Moreover, it is necessary to consider if there is a relationship between these two perturbation. And if we use the same data in real world to construct these two Brownian motions, they are very likely to be correlated [11]. And there is the previous work focusing on correlation of Brownian motions in dynamic systems. Hu et al. [12] consider two correlative stochastic disturbances in the form of Gaussian white noise in an epidemic deterministic model constructed by Roberts and Jowett [13]. Also, Hening and Nguyen [14] construct a stochastic Lotka-Volterra food chain system by introducing n number of correlated Brownian motions into the deterministic food chain model. n is the total species number in the food chain and they use a coefficient matrix to convert the vector of correlated Brownian motions to a vector of independent standard Brownian motions. Inspired by Emmerich [11], Hu et al. [12] and Hening and Nguyen [14], we are going to replace B 3 and B 4 by two correlated Brownian motions to introduce correlation of noises in SIS epidemic model. Considering two correlated Brownian motions, one with linear diffusion coefficient and the other with Hölder continuous diffusion coefficient, is clearly different from other work on stochastic SIS model. Though Hölder continuous diffusion coefficient and correlations of white noises are often involved in stochastic financial and biological models [15], there is no related work based on deterministic SIS model. As a result, this paper aims to fill this gap.
We now consider B 3 and B 4 in our model (2) to be correlated. Replace B 3 and B 4 by correlated Brownian motions E 1 and E 2 .
Note that E 1 and E 2 can be written as is a vector of independent Brownian motions and A is the coefficient matrix Also which gives the correlation of E 1 and E 2 ρ = a 1 a 2 , 0 < |ρ| < 1 Note that when ρ = 0, B 1 and B 2 are independent Brownian motion.
Substituting (4) into (3), we have with initial value I (0) = I 0 ∈ (0, N ) and this is our new model. Throughout this paper, unless otherwise specified, we let ( , {F t } t 0 , P) be a complete probability space with a filtration {F t } t 0 satisfying the usual conditions. We also define F ∞ = σ ( t 0 F t ), i.e. the σ -algebra generated by t 0 F t . Let B(t) = (B 1 (t), B 2 (t) T be a two-dimensional Brownian motion defined on this probability space. We denote by R 2 + the positive cone in R 2 , that is R 2 + = {x ∈ R 2 : x 1 > 0, x 2 > 0}. We also set inf ∅ = ∞. If A is a vector or matrix, its transpose is denoted by A T . If A is a matrix, its trace norm is |A| = trace(A T A) while its operator norm is ||A|| = sup{|Ax| : |x| = 1}. If A is a symmetric matrix, its smallest and largest eigenvalues are denoted by λ min (A) and λ max (A). In the following sections, we will focus on the long-time properties of the solution to model (5).

Existence of unique positive solution
We firstly want to know if the solution of our model (5) has a unique solution. Also, we need this solution to be positive and bounded within (0, N ) because it is meaningless for the number of infected population to exceed the number of whole population. So here we give Theorem 2.1.

Theorem 2.1
If μ + γ ≥ 1 2 (a 2 2 + a 2 3 )σ 2 2 N , then for any given initial value I (0) = I 0 ∈ (0, N ), the SDE has a unique global positive solution I (t) ∈ (0, N ) for all t ≥ 0 with probability one, namely, Proof By the local Lipschitz condition, there must be a unique solution for our SDE (5) for any given initial value. So there is a unique maximal local solution I (t) on t ∈ [0, τ e ), where τ e is the explosion time [15]. Let k 0 ≥ 0 be sufficient large to satisfy 1 Clearly, τ k is increasing when k → ∞. And we set τ ∞ = lim k→∞ τ k . It is obvious that τ ∞ ≤ τ e almost sure. So if we can show that τ ∞ = ∞ a.s., then τ e = ∞ a.s. and I (t) ∈ (0, N ) a.s. for all t ≥ 0.
Assume that τ ∞ = ∞ is not true. Then, we can find a pair of constants T > 0 and ∈ (0, 1) such that P{τ ∞ ≤ T } > So we can find an integer k 1 ≥ k 0 large enough, such that and . By Ito formula [15], we have, for any t ∈ [0, T ] and any k Then, we have which yields that Set k = {τ k ≤ T } for k ≥ k 1 and we have P( k ) ≥ . For every ω ∈ k , I (τ k , ω) equals either 1/k or N − 1/k and we have So the assumption is not reasonable and we must have τ ∞ = ∞ almost sure, whence the proof is now complete. Compared to the result from Gray et al. [10], the condition is now related to (a 2 2 + a 2 3 ). The square root terms are the reasons for us to give initial conditions in this section: when N − I (t) → 0, √ N − I (t) changes rapidly. This can also be an explanation to the readers that the condition is dependent on a 2 and a 3 instead of ρ = a 1 a 2 .

Extinction
The previous section has already provided us with enough evidence that our model has a unique positive bounded solution. However, we do not know under what circumstances the disease will die out or persist and they are of great importance in study of epidemic models. In this section, we will discuss the conditions for the disease to become extinction in our SDE model (5). Here, we give Theorem 3.1 and we will discuss persistence in the next section.

Theorem and proof
Theorem 3.1 Given that the stochastic reproduction number of our model < 1, then for any given initial value I (0) = I 0 ∈ (0, N ), the solution of SDE obeys lim sup if one of the following conditions is satisfied 16 a 2 2 σ 2 2 and 0 < ρ < 1 namely, I(t) will tend to zero exponentially a.s. And the disease will die out with probability one.
Proof Here, we use Ito formula LṼ (I (s)) ds and according to the large number theorem for martingales [15], we must have lim sup So if we want to prove lim sup t→∞ 1 t log I (t) < 0 almost sure, we need to find the conditions for LṼ (x) to be strictly negative in (0, N ). LṼ is defined by And it is clear that is ensured by R S 0 < 1. However, we do not know the behaviour of LṼ in (0, N ) and it is no longer quadratic as [10], which is very easy to analyse. As a result, we derive the first derivative of LṼ .
√ N In this case, the value of D(z) is always positive within z ∈ (0, √ N ), which leads to the similar result in Case 1. So we have We have 1 2 (a 2 2 + a 2 3 )σ 2 2 ≥ β + 9 16 a 2 2 σ 2 2 . In this case, D(z) will be positive in (0, √ N ), so LṼ increases when x increases. Similarly, extinction still maintains in this case.
In the deterministic SIS model, we have the result that if R D 0 < 1, the disease will die out. However, from our results in this section, we can see that our stochastic reproduction number R S 0 is always less that the deterministic reproduction number R D 0 = β N μ+γ , which indicates that the noise and correlation in our model help expand the conditions of extinction. For those parameters that will not cause the dying out of disease in the deterministic SIS model as well as Gray's model [10], extinction will become possible if we consider the new perturbation and the correlation.

Simulation
In this section, we use Euler-Maruyama Method [10,16] implemented in R to simulate the solutions in those 4 cases. For each case, we initially give a complete set of parameter to satisfy not only the extinction conditions, but also μ + γ ≥ 1 2 (a 2 2 + a 2 3 )σ 2 2 N to make sure the uniqueness and boundedness of solutions. Also, both large and small initial values are used in all 4 cases for better illustration. We then choose the step size = 0.001 and plot the solutions with 500 iterations. Throughout the paper, we shall assume that the unit of time is one day and the population sizes are measured in units of 1 million, unless otherwise stated (Figs. 1 2, 3, 4). The simulation results are clearly supporting our theorem and illustrating the extinction of the disease. Note that these conditions are not all the conditions for extinction. We only consider the situation that D(z) is either strictly positive or strictly negative. Otherwise, there will be much more complicated cases when LṼ is not monotonic in (0, N ).

Persistence
In this section, we firstly define persistence in this paper as there are many definitions in stochastic dynamic models to define persistence [3,4,10,15,[17][18][19][20]. However, our model (5) is based on [10]. As a result, we want to give a similar definition of persistence in our model (5). So here we give Theorem 4.1 to give a condition for the solution of (5) oscillating around a positive level.   (0, N ). I (t) will be above or below the level ξ infinitely often with probability one.
Proof When R S 0 > 1, recall that 3 2 We have LṼ (0) > 0 which is guaranteed by R S 0 > 1, LṼ (N ) = −(μ + γ ) < 0. As LṼ (x) is a continuous function in (0, N ), there must be a positive root of LṼ (x) = 0 in (0, N ). Moreover, from the behaviour of D(z), it is clear that LṼ will either increase to max value and then decrease to minimum, or increase to maximum, decrease to minimum and then increase to LṼ (N ) < 0. So In both cases, LṼ (x) = 0 will only have one unique positive root ξ in (0, N ).

Simulation
In this section, we choose the values of our parameter as follows: In order to prove the generality of our result, we use two different ρ, one positive and the other negative.
In both cases, we firstly use Newton-Raphson method [21] in R to find a approximation to the roots ξ of both LṼ , which are 7.092595 and 9.680572, respectively. Then, we use Euler-Maruyama method [10,16] implemented in R to plot the solutions of our SDE with both large and small initial values, following by using red lines to indicate the level ξ . The step size is 0.001, and 20000 iterations are used in each example. In the following figures, Theorem 4.1. is clearly illustrated and supported (Figs. 5, 6).

Stationary distribution
To find a stationary distribution of our SDE model (5) is of great important. From the simulation of the last section, we can clearly see that the results strongly indicate the existence of a stationary distribution. In order to complete our proof, we need to initially use a wellknown result from Khasminskii as a lemma. [22] Now we give the following Theorem 5.1 and the proof by using Lemma 6. Proof Firstly, we need to fix any (a, b) such that, 0 < a < ξ < b < N recall the discussion ofLV in last section, we can see that

Theorem and proof
also, recall (11) log and define Case 1 For all t ≥ 0 and any I 0 ∈ (0, a], use (20) in (11), we have From definition of τ , we know that Hence, we have Case 2 For all t ≥ 0 and any I 0 ∈ (b, N ), use (21) in (11), we have From definition of τ , we know that Hence, we have Combine the results from both Case 1 and Case 2, and we complete the proof of Theorem 5.1. Now we need to give the mean and variance of the stationary distribution.
Proof For any I 0 ∈ (0, N ), we firstly recall (5) in the integral form Dividing both sides by t and when t → ∞, applying the ergodic property of the stationary distribution [22] and also the large number theorem of martingales, we have the result that where m, m 2 are the mean and second moment of the stationary distribution. So we have We have tried to get other equations of higher-order moment of I (t) in order to solve m and v but fail to get a result. Though we do not have an explicit formulae of mean and variance of stationary distribution like [10], simulations can still be effective to prove (22).

Simulation
In this section, we use the same values of our parameters in the last section Now we simulate the path of I (t) for a long run of 200,000 iterations with step size 0.001 by using the Euler-Maruyama method [10,16] in R. And we only reserve the last 10,000 iterations for the calculations. These 10,000 iterations can be considered as stationary in the long term, so the mean and variance of this sample are the mean and variance of the stationary distribution of our solution. In both cases, the results of left side and right side of the equation (22) are 10.84587 and 10.68102, 1.743753 and 1.848107, respectively, so we can conclude that the mean and variance of the stationary distribution satisfy Eq. (22) . Figures 7 and 8 are the histograms and empirical cumulative distribution plots for each case of last 10000 iterations.

Conclusion
In this paper, we replaced independent Brownian motions in our previous model by correlated Brownian motions which leads to not only the increasing number of noises compared to Greenhalgh's work [4,10], but also turning the drifting coefficient into a nonlinear term. Then, we prove that the stochastic reproduction number R S 0 is the Keynes to define the extinction and persistence of the solution. Similar to the deterministic SIS model, with R S 0 < 1 and extra conditions, the disease will die out. When R S 0 > 1, we prove that the solution will oscillate around a certain positive level. Compared to [10], our LṼ is not linear, which results in more general and complicated conditions to both extinction and persistence sections. Moreover, compared to [23], this paper assumes that the Brownian motions are correlated, and hence, the effects of the correlations on the behaviours of our SIS system are studied. The analytical results including the form of R S 0 and the additional restrictions indicate that the correlations between the Brownian motions do make a significant difference. Though we do not know the explicit expression of that level, numerical method are then used to find the exact value under certain circumstance. Moreover, when R S 0 > 1, there is a unique stationary distribution of the solution. On the other hand, we have tried to get the explicit expression of the mean and variance histogram of I by deducing higher moments of I (t), but we seems not able to get results at this time. Consequently, we will leave this as a future work.