DEMOGRAPHIC STOCHASTICITY IN THE SDE SIS EPIDEMIC MODEL

. In this paper we discuss the stochastic diﬀerential equation (SDE) susceptible-infected-susceptible (SIS) epidemic model with demographic stoch- asticity. First we prove that the SDE has a unique nonnegative solution which is bounded above. Then we give conditions needed for the solution to become extinct. Next we use the Feller test to calculate the respective probabilities of the solution ﬁrst hitting zero or the upper limit. We conﬁrm our theoreti- cal results with numerical simulations and then give simulations with realistic parameter values for two example diseases: gonorrhea and pneumococcus.


1.
Introduction. Epidemics of infectious diseases have been a constant threat towards our society. In the past, Europe suffered from 25 million deaths out of a population of 100 million due to the Black Death [7]; Russia suffered from about 25 million cases of typhus with a death rate of about 10 percent, whilst smallpox wiped out half of the population of the Aztecs of three and a half million in 1520 [5]. Although in the 21st century, diseases such as smallpox no longer pose a threat towards mankind, there is still a high proportion of the population that is under threat of diseases such as malaria and HIV (e.g. [7,13]).
As a result mathematical models have been constructed in order to predict the behaviour of a disease and help to control a particular epidemic. Epidemics can be modelled by compartmental models such as Susceptible-Infected-Susceptible (SIS) and Susceptible-Infected-Removed (SIR) models where each individual has been assigned to a different subgroup representing a specific stage of disease. In 1927, Kermack and McKendrick [23] constructed the SIR model to describe the behaviour of diseases such as chickenpox and measles [18]. A typical individual starts off susceptible, at some stage catches the disease and after a short infectious period becomes completely immune. However, there are many diseases such as tuberculosis and meningitis [18] where a recovered individual will immediately become susceptible again, in other words individuals do not become immune to these diseases. For such a disease, an SIS model would be more suitable. There has been previous work done on the SIS epidemic model, for example Hethcote [18] studied the SIS epidemic model involving different factors such as disease mortality and migration. The SIS epidemic model is the simplest possible epidemic model and has been widely studied.
An epidemic of an infectious disease can be modelled by using either the deterministic model or the stochastic model. The deterministic model is often formulated as a system of differential equations where its solution is uniquely dependent on the initial value. On the other hand a stochastic model is a stochastic process with a collection of random variables where its solution is a probability distribution for each of the random variables. There has been much work done on deterministic models already, however there are some limitations in using these in analysing infectious diseases. A deterministic model is more appropriate when we are dealing with a large population. However, if we consider an epidemic outbreak in a small community such as school, a stochastic model would be more appropriate as the element of variability would become significant [7,8,10]. In addition, the real world is not deterministic, and there are many factors that can influence the behaviour of a disease and thus it is not always possible to predict with certainty what would happen. Consequently, a stochastic model is introduced to compensate for this problem. There are also many properties that are unique to the stochastic epidemic model which could enhance our understanding towards the behaviour of a particular disease. For example, the probability that an endemic will not occur, the final size distribution of an epidemic and the expected duration of an epidemic [2]. Clearly, we can see that introducing stochasticity into an epidemic model will provide some additional information that will improve the realism of our results compared to the deterministic approach.
There are three different types of stochastic models commonly used in population biology, namely the discrete time Markov chain (DTMC), continuous time Markov chain (CTMC) and stochastic differential equations (SDEs) [2,3]. In a DTMC model, the time and the state space variables are discrete. In a CTMC model, the time is continuous but the state variables are discrete, while the SDE is based on a diffusion process where both the time and the state variables are continuous [2]. These three stochastic models all consider the random behaviour occurring within the birth and death process of an individual [2,3]. Allen and Allen [3], made a thorough comparison between these three epidemic models with respect to the persistence time. They found that, when consistently formulated, the three stochastic models produced similar results for the mean and variance in persistence time.
In this paper, we shall focus on working with the SDE SIS epidemic model with demographic stochasticity.
The stochastic aspects of the SIS model for infectious diseases have been studied by many authors. In [11], Cavender considered the SIS model as an example of a birth and death process, which is a stochastic population model used to model demographic stochasticity [9]. Norden [32] described the stochastic SIS model as a stochastic logistic population model and aimed to investigate the distribution of the extinction times both numerically and theoretically. Kryscio and Lefévre [24] also looked at the stochastic SIS logistic model (also known as the stochastic SIS model). They extended and combined the results mentioned by Norden and Cavender. Kryscio and Lefévre obtained the approximations of the quasi-stationary distribution by studying two approximations of the process as well as the approximation to the mean time to extinction for the stochastic SIS logistic model. Furthermore, Clancy and Pollett [12] also considered the SIS logistic model as a birth and death process with a different death/recovery rate, µ i = µi, than the one mentioned in [24], namely µ i = µ(i − 1). By using Theorem 1 mentioned in Clancy and Pollett's paper, they have managed to prove one of the conjectures mentioned by Kryscio and Lefévre which Kryscio and Lefévre did not prove in their paper [24], namely q ≺ ST m, where q and m represent the quasi-stationary distribution and the stationary distribution of the process respectively. The notation ≺ ST represents the concept of stochastic ordering where it denotes the idea of one random variable bigger than another random variable. In other words, if A ≺ ST B where A and B are random variables then it means that P(A > x) ≤ P(B > x) for all x and P(A > x) < P(B > x) for some x, and that A is stochastically strictly less than B (e.g. [34]). Ovaskainen [33] looked at the other aspect of the quasi-stationary distribution of the stochastic SIS model with a different infection rate of susceptible individuals, λ i = λi(N − i), than the conventional one mentioned in most papers, namely λ i = (λ/N )i(N − i) [2,12,24], where N denotes the total population size. This is because by scaling λ by N , it implies that the number of contacts per person is independent of the population size while the one used by Ovaskainen [33] will take into account that there are more contacts per person in a large population than in a small population. Ovaskainen improved on previous approximation formulae and obtained a rigorous mathematical formula for the quasi-stationary distribution as N → ∞. Nasell [28] showed that for the stochastic SIS model with no demography, both the quasi-stationary distribution and the expected time to extinction from quasi-stationarity have three qualitatively different behaviours as a function of N and R 0 where R 0 is the basic reproduction number which determines whether or not a disease will die out or persist. Other than the papers we have mentioned above, there are still many other papers that dealt with the stochastic SIS model [4,29,30,37].
There are several ways that stochasticity can be introduced into the SIS epidemic model, one of which is introducing the random effects of the environment into how the disease spreads. This is known as environmental stochasticity. Dalal et al. [13] introduced environmental stochasticity into the disease transmission term of a deterministic model and looked at the effect of condom use on the spread of HIV amongst a homosexual population. In a similar way, Gray et al. [15] introduced stochasticity via perturbation of the disease transmission parameter in the standard SIS epidemic model. Moreover, inspired by Takeuchi et al. [14], Gray et al. [16] introduced the effect of telegraph noise, which is an example of environmental noise, into the SIS epidemic model by using a Markov switching process where the parameters switch between two or more regimes of environment and the switching is memoryless.
However, in this paper we shall introduce stochasticity into our SIS epidemic model by demographic stochasticity. This is when we introduce births and deaths into the population and derive a stochastic model. This is described as a set of differential equations for p i (t), the probability that there are exactly i susceptible individuals at time t. From these equations we can then derive the stochastic differential equations satisfied by the number of susceptible and infected individuals. Nasell [27] mentioned that an SIS model without demography is based on the assumption that an individual will live forever and is clearly unrealistic. So adding demography makes the SIS model more realistic and stochasticity makes the SIS model more realistic in a different way by adding random fluctuation in the population size, as the birth and death of an individual is a discrete and probabilistic event [26]. Some of the literature has dealt with the stochastic SIS model with demographic stochasticity, for example Andersson and Lindenstrand [6] looked at a two dimensional Markov process and analysed the behaviour of the model close to quasi-stationarity and the time it took for the system to become extinct with the help of a diffusion approximation. On the other hand, Nasell [27] focused on finding approximations of the quasi-stationary distribution and the time to extinction for his SIS model of the form of a bivariate Markov population process with appropriate transition rates. In addition, Nasell also derived an approximation for the expected time to extinction in the stochastic SIR model with demography, where he looked at two SIR models which each vary with a different demographic force.
The SDE SIS model with demographic stochasticity that we shall look at in this paper has been derived fully in detail in [1] with the model given in section 2. For this model, we have assumed that an infected individual who dies is immediately replaced by a susceptible individual and thus the population size is kept constant. This SDE SIS model is an approximation to the continuous time Markov Chain model (CTMC). The SDE model is preferred over the CTMC model when it comes to computing numerical simulations to illustrate the behaviour of the model [3], which we will do later on in this paper. In order to get a good estimation of the probability distribution for the CTMC model, computational costs can be very high. The SDE model is especially useful in situations where there are several random variables and several interacting populations [3]. As far as we know there have not been any detailed analyses on this model and we hope to fill this gap by providing a thorough analytical and numerical investigation of the model in this paper. Furthermore, other useful stochastic properties such as extinction will also be proven for this SDE SIS model with demographic stochasticity, in the hope to further enhance our knowledge on the impact of demographic variability on an epidemic.
The paper is organised as follows. In the next chapter we shall describe the basic model. In the following section we shall show the existence of a unique nonnegative solution. In section 4, we shall look at conditions for extinction and in section 5 we shall look at the Feller test which gives probabilities of hitting the top and bottom limits. In section 6 we perform some simulations with theoretical parameter values to verify the results and simulations with realistic parameter values for gonorrhea and pneumococcus.
2. Demographic stochasticity for the SDE SIS epidemic model. Throughout this paper, we let (Ω, F, {F t } t≥0 , P) be a complete probability space with filtration {F t } t≥0 satisfying the usual conditions (i.e. it is increasing and right continuous while F 0 contains all P-null sets). Let us consider the following deterministic SIS model with two populations S(t) and I(t) Here, S and I denote two populations representing respectively the number of susceptible and infected individuals in the population. N is the total size of the population, β is the disease transmission coefficient and β = λ/N where λ is the disease contact rate for each individual, that is the rate at which susceptible individuals come into contact with infected individuals. µ is the per capita death rate and γ is the rate at which an infected individual becomes cured. By looking at the interaction occurring between the two populations, Allen [1] constructed a list of possible changes with their corresponding probabilities p i (t), i = 1, 2, 3 . . . . In order to introduce demographic stochasticity into the deterministic SIS model (1), the mean change E(∆x) and the covariance matrix V for the time interval ∆t are calculated [1]. The stochastic SIS model has the form: In other words: If we write B = (W 1 −W 2 )/ √ 2 then B is a Wiener process so the SDE SIS model with demographic stochasticity becomes: In fact if in equations (5.8) and (5.9) on p.147 of [1] we replace γ by (µ + γ) and α by βN then equations (5.10) and (5.11) on p.148 of [1] are our equations (4). By using S(t) + I(t) = N , we can combine the two SDEs shown in (4) into one SDE for I(t), namely: The corresponding deterministic SIS model to the SDE SIS model (5) is given by: An alternative derivation of equation (5) based on (6) is given by Allen [2] who applies the procedure outlined above to (6). Equations (4) then follow from S + I = N . Note that the diffusion coefficient of the SDE SIS model (5) vanishes when I(t) = N + µ+γ β . It is appropriate to take an initial value I(0) = I 0 ∈ (0, N ). For the rest of the paper, we shall focus on analysing the SDE SIS model with demographic stochasticity (5). Throughout this paper, unless stated otherwise, we shall assume that the unit of time is one day, and the population sizes are measured in units of one million.
3. Existence of unique nonnegative solution. In order to prove the existence and uniqueness of the solution to the SDE SIS model (5), let us denote where λ(x) and σ(x) are the drift and diffusion coefficients of the SDE SIS model (5) respectively and x ∈ [0, N + µ+γ β ]. We shall now extend the domain of our SDE SIS model (5) into the whole domain, i.e. λ(x), σ(x): R → R by considering the following equations: and As this SDE is a special case, the standard existence and uniqueness theorems on SDEs are not applicable here (e.g. [25]) and other methods [15] do not adapt here as well. We are now ready to prove the existence and uniqueness of the solution for the SDE SIS model (5) by using the following lemma which is mentioned in [21] (Theorem 3.2 of Chapter IV).
, the SDE (9) with its coefficients defined by (7) and (8) has a strong pathwise unique solution.
Proof. We shall split the proof into two sections by showing condition (i) is satisfied first.
(i) The first derivative of equation (7) is defined as: For x, y ∈ 0, N + µ+γ β , by the Mean Value Theorem we have that for some It is easy to see that the result follows for x, y in (−∞, ∞). Therefore, condition (i) is satisfied with κ(u) = M u for some constant M for all x, y ∈ R and that λ(x) is Lipschitz continuous.
We will now complete the proof by looking at the second condition: (ii) From equation (8), it is clear that the Mean Value Theorem does not apply in this case. In addition, if we are able to show condition (ii) is satisfied for x, y ∈ [0, N + µ+γ β ] then the rest will follow. In other words, if we could find a constant L such that for x, y ∈ [0, N + µ+γ β ], then condition (ii) is proved. By choosing ε = 1 4 (N + µ+γ β ) and considering separately the regions ε ≤ and 0 < x < ε, it is straightforward to show that (10) holds. As a result, condition (ii) is satisfied with ρ(u) = L √ u for some constant L for all x, y ∈ R. In other words, σ(x) is Hölder continuous with exponent 1/2. Moreover, by definitions (7) and (8) both λ(x) and σ(x) are bounded so the theorem follows from Lemma 3.1.
Note that by Theorem 2.4 of Chapter IV of [21], since σ(I(t)) and λ(I(t)) are bounded, the solution to the SDE SIS model (5) will not explode. Hence, we have shown that a unique strong pathwise non-explosive solution does in fact exist for our SDE SIS model (5).
All we have left to show now is the non-negativity of our solution. Note that in order for this SDE SIS model (5) to make sense, the term inside the square root has to be nonnegative. We consider the SDE (9). We show that provided that I 0 ∈ (0, N ) then I(t) ∈ [0, N + µ+γ β ]. Theorem 3.3. For any given initial value I(0) = I 0 ∈ (0, N ), the probability that the SDE (9) has a unique and nonnegative solution almost surely for all t ≥ 0.
Note that this result is slightly unusual as based on biological considerations we would have expected I(t) ∈ (0, N ). This is the result of introducing stochasticity into the model using the technique mentioned in [1] and this SDE SIS model (5) is a well-established model. In section 6, we shall show that although it is theoretically possible for I(t) to exceed N for simulations carried out with realistic parameter values, in practice I(t) never went over N .
Proof. The following proof for Theorem 3.3 is established based on the framework of the "Square Root Process" illustrated in [25]. We shall divide this proof into two parts. First of all, we shall prove the left hand side inequality, I(t) ≥ 0 and by using a similar strategy we shall prove the right hand side inequality and thus complete the proof. Let a 0 = 1 and a k = e −k(k+1)/2 for every integer k ≥ 1. Note that Let Ψ k (u) be a continuous function such that its support is contained in the interval (a k , a k−1 ) where It can be shown that such a function exists. Define ϕ k (x) = 0 for x ≥ 0 and It is easy to see that ϕ k ∈ C 2 (R, R). Furthermore, , respectively. As in [25]: where we define x − = −x if x < 0 or otherwise x − = 0. Now, by Itô's formula, we get that for any t ≥ 0: where λ, σ : R → R are defined as before.
Taking the expectation yields: Now for all t, I − (t) ≥ 0, so EI − (t) ≥ 0, hence from our result (13), we must have: Furthermore, by using equation (14) and proof by contradiction, it is straightforward to show that for all t > 0, Therefore, I(t) ≥ 0 almost surely and this completes the left hand side of the proof for equation (11).
To complete the proof we shall now show that I(t) ≤ N + µ+γ β . Let us define Then, from Itô's formula on equation (15), we get: for J ≥ N + µ+γ β . By Itô's formula, we derive that: where P, Q : R → R are defined by: So P (x) ≤ 0 and Q(x) = 0 for all x. Thus Now take the expectations to get Eϕ k (J(t)) ≤ 0. Hence, Similarly to the argument we used for proving the left hand side of equation (11), it is clear that for all t > 0, P(J(t) < 0) = 0, ⇒ P(J(t) ≥ 0) = 1, i.e., I(t) ≤ N + µ+γ β almost surely ∀t ≥ 0, which completes the entire proof.
Hence we have proven the nonnegative property of the solution of the SDE SIS model with demographic stochasticity (9) and provided that I(t) ∈ [0, N + µ+γ β ] we could express equation (9) as equation (5). This has completed our proof on the existence of a unique nonnegative solution for the SDE SIS model (5). 4. Extinction of our solution. In this section, we shall focus on the extinction aspect of the nonnegative solution I(t) ∈ [0, N + µ+γ β ] to the SDE SIS model (5). Let us define the basic reproduction number R 0 as: where the parameters as denoted as before.
Theorem 4.1. For any given initial value I(0) = I 0 ∈ (0, N ), if R 0 ≤ 1 or if R 0 > 1 and N < 1 4 + µ+γ β , then I(t) will hit zero with probability one in finite time. In other words, the disease will certainly die out in finite time.
Proof. Let us define the stopping time τ n = inf{t : I(t) ≤ n} for 0 ≤ n < I 0 , where we set inf ∅ = ∞. We need to show that We will show this by using proof by contradiction. If (16) were false, then P(τ 0 = ∞) > 0. Noting that lim n→0 τ n = τ 0 , we could find an n sufficiently small so that δ := P(τ n = ∞) > 0.
By Itô's formula, we have that d( I(t)) = I(0) + q(I(t))dt for 0 ≤ t ≤ τ n , where q : R → R is defined by: In order to prove this theorem, we will now show that there exists a negative upper bound for q(x) when x ∈ [n, N + µ+γ β ]. • Case 1. R 0 = βN µ+γ ≤ 1 . By using Theorem 3.3, it is clear that the second term in equation (19) is negative. In addition, due to the fact that βN µ+γ ≤ 1, then the first term for (19) is negative, which makes q(x) in (19) negative. Therefore, we have that for x ∈ [n, N + µ+γ β ], As a result from (20), we could conclude that for x ∈ n, N + µ+γ β , q(x) is negative and thus there must exist some ε > 0 such that q(x) < −ε < 0 for t ≥ 0 for x ∈ n, N + µ+γ β . Now by substituting the negative upper bound of q(x) into equation (18) and integrating, we get that: By taking the expectation of equation (21) and using the result given by (17), we obtain that: Now letting t → ∞, we have that 0 ≤ −∞, which clearly is a contradiction. Therefore the result given by (16) must be true, in other words the disease will die out almost surely for the case where R 0 ≤ 1. Similarly, we shall apply the same argument to the second case: • Case 2. R 0 = βN µ+γ > 1 where N < 1 4 + µ+γ β . First of all we shall rewrite equation (19) as: where U (x) is defined as: for x ∈ n, N + µ+γ β . Clearly, U (x) is a quadratic function so therefore it must have at most two real roots. From (24), We can choose n sufficiently small, thus making U (n) negative. Additionally Also, U (x) has a maximum turning point at x * = 1 suppose that x * ∈ n, N + µ+γ β , then by substituting x * into (24) we could see that the second term of (24) is negative. Furthermore, the first term of (24) becomes: which is negative as N < 1 4 + µ+γ β . So U (x) < 0 when x ∈ n, N + µ+γ β . Now consider the case where x * > N + µ+γ β . By using the fact that U (n) < 0 and U (N + µ+γ β ) < 0 and that U (x) has one unique turning point, at x = x * , U (x) is negative for x ∈ n, N + µ+γ β . By combining both results we could conclude from equation (4.10) that, for R 0 > 1, ∃ε > 0 such that q(x) < −ε < 0 for Arguing as in case 1, we deduce that P(τ 0 < ∞) = 1. This completes the proof.

5.
Probabilities of hitting the top and bottom limits. Now that we know under certain situations, the number of infected individuals will die out, it is also useful to know the probability of it hitting zero and the probability of it hitting N + µ+γ β . For the rest of this paper, we shall work on the SDE SIS model (5) Proof. The proof of Theorem 5.1 is established based on the framework of the "Mean Reverting Square Root Process" illustrated in [25]. Let us define the drift and diffusion coefficients of our SDE SIS model (5) as and respectively.
As mentioned in [25], we know that for any given pair of nonnegative constants a and b with a < I 0 < b there is a unique solution, say M (x) satisfying the equation with boundary conditions M (a) = M (b) = 0. This equation for M (x) is solved in [22] and we shall outline the important aspects of the working for the purpose of completeness. Let us introduce the speed measure and the Green function where the scale function p(x) is defined in [22] as where c ∈ R is a fixed number. This scale function p(x) is a monotonic increasing function of x. By using equations (31), (32) and applying the boundary conditions, we could obtain the explicit solution M (x) as illustrated in [22] that satisfies the equation (30), namely: Since G a,b (x, y) is a non-negative function, it is clear that M (x) ≡ M a,b (x) is also a non-negative function. Now, let us define the stopping times: τ a = inf{t ≥ 0 : I(t) ≤ a}, τ b = inf{t ≥ 0 : I(t) ≥ b}, where a < I 0 < b. By the Itô formula we get that: Taking the expectations yields the following results which are similar to the ones that have been illustrated in [25]: This indicates that I(t) exits from every compact subinterval of 0, N + µ+γ β in finite expected time, which means that we must have P(τ a ∧ τ b < ∞) = 1. In addition, by referring to the boundary conditions we get that lim t→∞ EM (I(t ∧ τ a ∧ τ b )) = 0, and so E(τ a ∧ τ b ) = M (I 0 ). Let us now define where x ∈ 0, N + µ+γ β and we define x 0 = 1 2 (N + µ+γ β ). This function has continuous first and second derivatives V (x) and V (x) in (0, N + µ+γ β ) with strictly nonnegative V (x), and V (x) satisfies where v(x) and w(x) are defined as equation (28) and (29) respectively. By the Itô formula, we could derive that: Taking the expectations and letting t → ∞ yields that: By using the fact that the two probabilities must add up to one, we obtain from equation (34) that: and where equation (35) represents the probability of hitting a before it reaches b and vice versa for equation (36). Now by substituting v(z) and w(z) into equation (33), we get that: For the case where 4(µ+γ) β ≥ 1, we have that V (x) tends to a finite (strictly negative) limit as x → 0 + and that V (x) tends to infinity as x → N + µ+γ β − , namely: Recall that:

DEMOGRAPHIC STOCHASTICITY SIS EPIDEMIC 2873
Define τ [a,b] = τ a ∧ τ b . From equations (35) and (36) and the above notations, we can work out the probability of I(t) hitting the bottom limit is as follows: For a ∈ 0, N + µ+γ .
But since this holds for any a > 0 we must therefore have Similarly, the probability of I(t) hitting the top limit N + µ+γ β is: By letting a ↓ 0 we get that But since this holds for any b < N + µ+γ β , then letting b ↑ N + µ+γ β , As a result, for the case 4(µ+γ) β ≥ 1, I(t) will reach 0 first before it reaches N + µ+γ β almost surely. By applying a similar argument for the case where 4(µ+γ) β < 1, we get that as x → 0 + , V (x) tends to a finite strictly negative limit, whereas as x → N + µ+γ β − , V (x) tends to a finite strictly nonnegative limit. In other words: Furthermore, the probability of I(t) reaching 0 before it reaches N + µ+γ β is given as follows: Letting b ↑ N + µ+γ β in equation (37), we get that: and since this holds for any a > 0, we have that Similarly, the probability that I(t) reaches N + µ+γ β before it reaches 0 is given as: >0.
If P(τ < ∞) = 1, then and thus the inequalities (38) and (39) are actually equalities. This indicates that wherever I(t) starts, there is a nonnegative probability that I(t) will first hit each of zero and N + µ+γ β . If I(t) starts exactly halfway between zero and N + µ+γ β , then there is a higher probability that I(t) will hit zero before it hits N + µ+γ β . However, if P(τ = ∞) > 0 then all that we can say is as described in inequalities (38) and (39), namely Hence in this section we have used the Feller test to calculate the probabilities that I(t) will hit zero before it hits N + µ+γ β and vice versa. In the next section we shall look at some of our analytical results using computer simulations. 6. Simulations. In this section we shall use the Milstein numerical simulation method for SDEs (e.g. [35]) to numerically illustrate Theorem 4.1 and Theorem 5.1. The Milstein method is superior to the simpler Euler-Maruyama method, for example used in [15], because as the integration time-step goes to zero the Milstein method is strongly convergent with order 1 as opposed to 0.5 for the Euler-Maruyama method [20]. Our numerical integration program was written in R and comprehensively verified using a large number of runs. For Theorem 4.1, we first show that the disease will die out in finite time if R 0 ≤ 1, or R 0 > 1 and N < 1 4 + µ+γ β , and explore numerically the situation where R 0 > 1 and N ≥ 1 4 + µ+γ β . 6.1. Simulations on extinction. In this section, we shall focus on highlighting the results shown in Theorem 4.1.
Moreover, by substituting the parameters (40) into the corresponding SIS deterministic model (6), we have: By applying the Milstein method on the SDE SIS model (41) and its corresponding SIS deterministic model (42), we have managed to construct the computer simulations illustrated in Figure 1 for parameters given by (40). Figure 1 illustrates two different simulations constructed with different initial values. The simulation on the left hand side represents the behaviour of the model when I(0) = 90, while the one on the right hand side represents the behaviour of the model when I(0) = 1. For both cases we could see that no matter what we choose our initial value to be, I(t) will eventually die out and hit zero and thus the disease will go extinct. The simulation was repeated with many different parameter values satisfying the condition R 0 ≤ 1 and different initial conditions and in each case the disease died out in finite time. This supports the results of Theorem 4.1 on extinction.
The simulation was repeated with many different parameter values satisfying R 0 > 1 and N < 1 4 + µ+γ β and in each case the disease died out in finite time as predicted by Theorem 4.1. One such simulation is shown in Figure 2 with parameter values as in (43). and It is easy to see that 4(µ+γ) β = 440 > 1 and thus from Theorem 5.1, we conclude that the disease hits zero before N + µ+γ β . The numerical simulations support these results as expected. Two typical simulations are shown in Figure 3. The numerical simulations were repeated with a variety of parameter values and initial conditions.  and its corresponding SIS deterministic model (6) becomes: For this example, Theorem 5.1 says that it is possible for I(t) to hit either zero or N + µ+γ β first. Figure 4 illustrates simulations which clearly shows that this is the case.  In sections 6.1 and 6.2, we have been focusing on using arbitrary parameters to support our theories proved in Theorems 4.1 and 5.1 respectively. However, it would be better to use parameters for reallife diseases. In this section, we shall look at two different diseases for which an SIS model is suitable: gonorrhea amongst homosexuals and pneumococcus amongst very young children in Scotland. We shall first look at gonorrhea amongst homosexuals. Throughout the section the unit of time is still one day but the population sizes are not scaled as previously. Note that the demographic SDE SIS epidemic model (5) is a well established model. This model approximates the system of ordinary differential equations describing the probabilities that there are exactly i infected individuals at time t by a single stochastic differential equation. In the system of ordinary differential equations i never exceeds N but we have shown that in the stochastic differential equation approximation I may possibly exceed N . However, it is important to note that when carrying out the simulations with realistic parameter values, in practice I(t) never exceeded N , although the theoretical possibility remains that it could do so.
For this case, as R 0 > 1 and N > 1 4 + µ+γ β , Theorem 4.1 is inconclusive. For this case Theorem 5.1 predicts that I(t) will almost surely hit zero in finite time. However, as the time realistically looks likely to be very high it is not feasible to run the simulations for that long. For the simulations shown and the other simulations not shown with both different starting values, and different realistic parameter values with R 0 > 1 and N ≥ 1 4 + µ+γ β , after an initial transient stage the stochastic simulations oscillated about the deterministic level.
(51) Clearly in this case, R 0 = 0.98 < 1 when we could conclude from Theorem 4.1 that for any given initial value I(0) ∈ [0, N + µ+γ β ], the solution I(t) of the SDE SIS model (5) will die out almost surely with probability one. Furthermore 4(µ+γ) β > 1 whence, from Theorem 5.1 we could also conclude that I(t) will hit zero before N + µ+γ β with probability one. The simulation produced by the Milstein method for SDE SIS model (4) and the corresponding SIS deterministic model (5) with parameters given by (51) supports both Theorems 4.1 and 5.1. In other words the disease almost surely hits zero before the upper bound and hits zero in finite time almost surely.
The numerical simulations were repeated with different values of N where R 0 < 1, and similar results were obtained each time.
Next, we shall look at pneumococcus, especially focussing on children under two years old in Scotland mentioned in Greenhalgh, Lamb and Robertson [17].
Example 6.7. (Pneumococcus Model) In Greenhalgh, Lamb and Robertson's paper [17], they have chosen N = 150, 000, µ = 1/104 per week = 1.37363×10 −3 /day. In Weir's thesis [36], she chose γ = 1/7.1 per week = 0.02011/day and in Zhang et al. [39] they chose β = 2 × 10 −6 /week = 2.857 × 10 −7 /day. It is easy to see that in this case R 0 = 2 > 1 and N > 1 4 + µ+γ β and thus from Theorem 4.1 we are unable to conclude anything. For these parameter values 4(µ+γ) β > 1 and so ultimately the disease goes extinct, but the time taken for this to happen again seems very large. The numerical simulation produced by these parameters is shown in Figure  6. Again for other simulations not shown with different initial values and different parameter values with R 0 > 1 and N ≥ 1 4 + µ+γ β , after an initial transient stage the stochastic simulations oscillated about the deterministic level.
For illustrative purposes, we change N to 68, 000 so that R 0 = 0.907 < 1 and that 4(µ+γ) β > 1. The numerical simulation produced for this case support both our results in Theorems 4.1 and 5.1 and thus the disease almost surely hits zero before the upper bound and hits zero in finite time almost surely. Again the numerical simulations were repeated with different values of N where R 0 < 1, and similar results were obtained each time.
As we mentioned at the beginning of this section the theoretical results show that if we use the stochastic differential equation approximation suggested by Allen [1] to incorporate demographic stochasticity it becomes theoretically possible for the number of infected individuals to exceed the population size and the number of susceptibles to become negative. This may make us question whether the model is practically useful. However extensive simulations with realistic parameter values for real diseases were performed (some examples have been illustrated above) and in these simulations we never once actually observed the number of susceptibles become negative, although it remains a theoretical possibility. Thus this approximate model may still be useful to illustrate the effect of inherent stochasticity in population dynamics.  Figure 6. Computer simulations of the path I(t) for the Pneumococcus Model with parameters N = 150, 000, µ = 1.37363 × 10 −3 /day, γ = 0.02011/day and β = 2.8650 × 10 −7 /day, using the Milstein method with step size ∆ = 0.001 and with initial value I(0) = 70, 000 7. Conclusion and discussion. The use of epidemic models to control infectious diseases is becoming increasingly common. The SIS epidemic model is one of the simplest epidemic models possible and has been widely used practically to predict the spread of infectious diseases such as gonorrhea and pneumococcus and examine the effect of control strategies. However it ignores random variability in the population. One way to include random variation into the SIS epidemic model is to model the transitions as Markov processes with the appropriate rates and then either perform Monte-Carlo simulations, or derive the differential equations satisfied by p i (t), the probability that there are exactly i individuals infected by the disease at time t. The latter approach is illustrated by Bailey [7]. However for realistically large population sizes these approaches rapidly become very cumbersome and use a lot of computational power. Allen [1] suggested to use a stochastic differential equation approximation to simplify the analytical stochastic model so that one had essentially a single stochastic differential equation instead of a very large set of ordinary differential equations. This model has previously been formulated but never analysed. In this paper we have filled this gap. We showed that this SDE SIS epidemic model has a unique nonnegative bounded solution. Then we derived sufficient conditions for the disease to go extinct in a finite time. This behaviour is different than the behaviour for the SIS model with environmental stochasticity studied in [15] where environmental noise altered the threshold value R 0 from the deterministic model. If the stochastic threshold value R s 0 exceeded one then the disease would persist and oscillate about a non-zero level. In our model, the demographic noise does not alter the threshold value. Next we used the Feller test to establish the probabilities of the number of infectious individuals hitting the lower and upper boundaries. Finally we used numerical simulations to confirm our analytical results and examine the behaviour of the model for realistic parameter values for gonorrhea and pneumococcus.
The analytical results show that it is theoretically possible for the number of susceptibles to become negative in the solution to the stochastic differential equation model. However in many simulation runs with realistic parameter values this was never actually observed so the stochastic differential equation model remains a useful approximation to illustrate the possible effects of demographic stochasticity on population dynamics.