Persistence-time estimation for some stochastic SIS epidemic models

In this paper, we study two stochastic SIS epidemic models: the first one with a constant population size, and the second one with a death factor. We analyze persistence and extinction behaviors for these models. The persistence time depends on the initial population size and satisfies a stationary backward Kolmogorov differential equation, which is a linear second-order partial differential equation with variable degenerate coefficients. We solve this equation numerically using a classical finite element method. We give computational evidence that the importance of understanding the dynamics of both the deterministic and the stochastic epidemic models is due to the numerical approximations to the mean persistence time. This can give more information about the model and may perhaps explain strange behaviors, such as the differences between the deterministic model and the stochastic one for long times.


(Communicated by Thomas Wanner)
Abstract. In this paper, we study two stochastic SIS epidemic models: the first one with a constant population size, and the second one with a death factor. We analyze persistence and extinction behaviors for these models. The persistence time depends on the initial population size and satisfies a stationary backward Kolmogorov differential equation, which is a linear second-order partial differential equation with variable degenerate coefficients. We solve this equation numerically using a classical finite element method. We give computational evidence that the importance of understanding the dynamics of both the deterministic and the stochastic epidemic models is due to the numerical approximations to the mean persistence time. This can give more information about the model and may perhaps explain strange behaviors, such as the differences between the deterministic model and the stochastic one for long times.
1. Introduction. A general stochastic differential system for the dynamics of n interacting populations has the following form: dX(t) = µ(t, X)dt + B(t, X) dW(t), t > 0, where X(t) = (X 1 (t), ..., X n (t)) T is an n-dimensional random variable, X 0 is the initial population size, and W(t) = (W 1 (t), . . . , W m (t)) T are m independent Wiener processes. The vectorial function µ(t, X) is called the drift, and the matrix B(t, X) is the diffusion matrix. Note that if B = 0, then this system is reduced to a standard deterministic model for the population dynamics. The stochastic model (1) can be simulated computationally, for example by the Euler-Maruyama method, which is simple and straightforward to implement, and has strong order of convergence equal to 1/2, and weak order equal to 1. The traditional convergence theory for numerical methods applied to stochastic differential equations (SDEs) requires a global Lipschitz assumption on the drift and diffusion coefficients. However, in practice, many important SDE models satisfy only a local Lipschitz property and, since Brownian paths can make arbitrarily large excursions, the global Lipschitz based theory is not directly relevant; more details, along with further references, can be found in [13]. We refer to [15] for a review on numerical solutions to SDEs.
One of the major goals of studying stochastic population dynamics is to predict finite extinction time, i.e., the time when one or more components of the population go to zero. In this paper, T is the random variable that indicates the extinction time: T := inf{t ≥ 0 : X i (t) = 0, for any i = 1, . . . , n}.
(2) Note that there are other definitions of persistence or extinction: weak persistence, weak persistence in the mean, strong persistence in the mean, etc. (see for example [17, p. 444]). Obviously, T depends on the initial population size X 0 .
On the other hand, the expected time until extinction or mean extinction time is given by where E(T ) is the expectation of T , and the p−th moment for p > 1 is τ (p) := E(T p ). It is known that the mean extinction time τ and the p−th moments τ (p) can be computed by solving a suitable stationary backward Kolmogorov equation (see [1, p.150] or [4, p.443]; we refer to [3, Ch. 8] for a clear review on these equations).
The backward Kolmogorov equations are linear second-order partial differential equations (PDEs) with variable degenerate coefficients. It it well known that it is possible to find an analytical solution only in a few simple cases (see for example [4], for a population model with linear birth and death rates, or [10], where the authors give exact expressions for the extinction time of a class of birth-death processes). In general, we can only compute numerical approximations, as for example in [2], where central difference approximations are used. Recently, in [9], a finite element method (FEM) has been chosen. There is an ample bibliography on finite element methods (see for instance [14] for a good introduction, and the classics [7] and [8] for details), and their implementation is quite standard (see for instance [11]). In [9], after transforming the problem into a well-suited variational form, the PDE is discretized by means of rectangular elements with cubic basic functions; finally, a direct solver is applied to the resulting sparse linear system, which can be formulated very conveniently in Matlab c by means of sparse matrices. In this paper, we have followed two approaches: one similar to that in [9], and a new implementation using FreeFem++ [12].
The paper is organized as follows. In Section 2, we consider a classic SIS epidemic model with constant population size. In Section 3, we introduce a death factor over the infected population. Finally, in Section 4, we draw the main conclusions.
Our numerical methods have been implemented in Matlab c and FreeFem++. The experiments have been carried out in an Intel(R) Core(TM)2 Duo CPU E6850 @ 3.00GHz. The codes for the numerical tests are available on request.
2. A SIS model with constant population size. In this section, we consider a SIS epidemic model for a single species, such that the susceptible individuals become infected, recover, and become susceptible again; this means that there is no immunity to the disease. We also describe the equilibrium for the deterministic part of the model, and present its stability properties depending on the so-called basic reproduction number. Following [1, p. 148], the stochastic SIS model has the following form: where S = S(t) is the susceptible population size at time t; I = I(t) is the infected population size at time t; s and y are the initial population sizes; W 1 (t) and W 2 (t) are two independent Wiener processes; and the drift and diffusion functions are given respectively by µ(S, I) = − αSI S + I + γI, ρ(S, I) = 1 2 The parameter α > 0 is the contact rate for transmitting the infection, i.e., the average number of individuals with whom an infected individual makes sufficient contact to pass an infection, and γ > 0 is the probability that an infected individual is removed from the infection process during a unit time interval (relative removal rate). Note that the sum S(t) + I(t) = N is constant for t ≥ 0; this means that it is a model without demography. In fact, the system (4) can be reduced to a single stochastic differential equation for I: nevertheless, this reduction will be impossible for the model with deaths that we will consider later in Section 3.

2.1.
Extinction time for the one-dimensional problem. The deterministic part of (6) is There are two fixed points: I = 0 and I = N (R − 1)/R, where R = α/γ is called the basic reproductive number. This number is defined as the average number of secondary infections caused by one infective individual during his/her infective period in an entirely susceptible population, and determines the asymptotic behavior with two different possibilities: (i) If R ≤ 1, solutions to (7) approach the disease-free equilibrium, because lim t→∞ I(t) = 0.
(ii) If R > 1, solutions to (7) approach a unique positive endemic equilibrium, because lim t→∞ The time to disease extinction T = inf{t ≥ 0 : I(t) = 0} is the extinction time for the SDE (6), and its expected value τ = τ (y) := E(T ) is the solution to the following boundary value problem for ordinary differential equations (ODEs): We have solved numerically this problem in Matlab c using Chebyshev differential matrices, as is explained in [21,Ch. 7]. In Figure 1, we have plotted, in black, the numerical solutions for N = 500; in blue, for N = 1000; and, in red, for N = 2000.
On the left-hand side of Figure 1, α = 0.25, γ = 0.75, so R = 1/3 < 1; on the right-hand side, α = 0.50, γ = 0.49, so R > 1. The difference in the type of the solutions is evident: whereas in the first case (left), the values of τ depend basically on y = I(0); in the second case (right), the solution changes depending on how large N is. Moreover, additional tests show that this situation is exacerbated when the basic reproductive number R increases.

2.2.
Extinction time for the two-dimensional problem. The deterministic version of (4) leads to the following system of ordinary differential equations: Theorem 1 from [5] proves that the asymptotic behavior of (9) has two possibilities: (i) If R ≤ 1, solutions to (9) approach the disease-free equilibrium, because lim The Euler-Maruyama algorithm applied to (4) is straightforward and may be described as follows: Algorithm 2.1. Given α and γ, let t be the time-step, nrun, the number of simulations, and T ol 1 , T ol 2 , the tolerances. Then, we have: Generate two random numbers w 1 , w 2 normally distributed (0, 1) Compute µ n and ρ n using (5) Compute the mean value and standard deviation of extinction_time In Table 1, we present the results of 10000 trials of Algorithm 2.1, for α = 0.25, γ = 0.75 (R < 1). We can observe that I n always hits zero, and, even if S(0) varies widely, the results change only slightly. In Table 2, we give the results of 2000 trials of Algorithm 2.1, for α = 0.50, γ = 0.49 (R > 1). Algorithm 2.1 obtains the mean and the standard deviation of the stopping time T * := inf{t ≥ 0 : S(t) = 0 or I(t) = 0}, although the biological interest is in the time to disease extinction T E := inf{t ≥ 0 : I(t) = 0}. In order to obtain T E , it is almost trivial to write an algorithm applying the Euler-Maruyama method to the SDE (6); but all the numerical experiments indicate that there is no difference between both algorithms, because, in Algorithm 2.1, I(t) is always the variable that hits zero, while S(t) never does.
As we have already said at the beginning, in order to compute the extinction time, we will use other techniques, related to the backward Kolmogorov equations. More precisely, the stationary backward Kolmogorov equation for (4) is the partial differential equation M S and M I being positive constants. Note that this procedure can be extended to higher-order moments τ (p) (s, y), which are solution to the following PDE: defined over Ω = [0, M S ] × [0, M I ], and with the same boundary conditions as in (11). In order to apply the FEM technique, we need to write a well-suited variational form for (10)- (11), and make formal computations. Let us multiply (10)   Integrating by parts the second order terms, and taking into account the Neumann boundary conditions from (11), we obtain the following weak formulation: where the last boundary term cannot be removed or replaced by the boundary conditions in (11).
Using similar ideas to those in [9], i.e., discretizing by means of rectangular elements with quadratic polynomials as basic functions, and a nodal distribution of the form  Table 1, which contains the results of 10000 trials of Algorithm 2.1. I n always hits zero, and, although S(0) varies widely, the results change only slightly. In Figure 3, we have plotted, on the left-hand side, three cuts of the numerical solution to (10)-(11), at s = 889 (blue), s = 121 (green), and s = 15 (magenta). The blue points represent the means of 1000 trials of Algorithm 2.1, with T ol 1 = T ol 2 = 1, for S(0) = s = 899; and they give, in our opinion, an acceptable approximation. On the right-hand side, we have also plotted, with black dotted stroke, the numerical solution to (8) with N = 1000, in order to check that the differences between both solutions are small enough.
On the other hand, in this example (α = 0.25, γ = 0.75), the basic reproduction number is R = 1/3 < 1, so the solution to the corresponding deterministic system (9) approaches the disease-free equilibrium, i.e., lim t→∞ I(t) = 0; this same behavior  11), is compared with the blue points, which represent the means of 1000 trials of Algorithm 2.1. Additionally, we have plotted cuts at s = 121 (green), and s = 15 (magenta). Right: the numerical solution to (8), with N = 1000, is also plotted, using black dotted stroke.
is observed in the stochastic model, with Algorithm 2.1. In all the numerical tests with R < 1, a disease-free equilibrium has been achieved in finite time, and the graphics are very similar to those shown in Figure 2.
The case R > 1 is different and more complicated. In  tests show that this situation is exacerbated when the parameter R increases. Observe that, in that case, the situation is very different, because the contour lines on the right-hand side of Figure 4 are roughly S(0) + I(0) = N , i.e., the extinction time now depends basically on the size of the population, N , whereas in the previous case, it depended on the initial number of infected individuals I(0) = y, conclusion already noted in [16] and [20]. This behavior is confirmed in Table 2, with the results of 2000 trials of Algorithm 2.1. N = 1000, and I n always hits zero. The elapsed time in all cases was greater than one hour.
In Figure 5, we have plotted, on the left-hand side, three cuts of the numerical solution to (10)-(11), at s = 1005 (blue), s = 500 (red), and s = 101 (green); the blue and green points represent the means of 1000 trials of Algorithm 2.1, for S(0) = s = 1005, and S(0) = s = 101, respectively, and give, in our opinion, an acceptable approximation. On the right-hand side, we have replotted the left-hand side, but adding the numerical solution to (8), with N = 2000 (dotted black), and N = 1000 (dotted cyan). It is important to observe that the blue and black curves intersect at about I(0) = y = 1000, and the red and cyan ones, at I(0) = y = 500, as expected.
Although our codes perform well, the previous remarks suggest that a more specific software should be convenient. More precisely, we have chosen FreeFem++ [12], which is freely available. FreeFem++ greatly simplifies the implementation of FEM, and is particularly efficient. Therefore, we are able to consider much larger domains, with much more nodes.
In Figure 6, we plot, on the left-hand side, the numerical solution with our Matlab c implementation (blue), and the numerical solution with FreeFem++ (magenta), for S(0) = s = 1000, and 0 ≤ I(0) = y ≤ 300; and, on the right-hand side, a zoom for 0 ≤ I(0) = y ≤ 70. We can clearly see that the solutions with both Matlab c and FreeFem++ are virtually identical, and that the differences with respect to the results of Algorithm 2.1, depicted with blue points, become bigger with I(0). Moreover, the red line is the standard deviation obtained with FreeFem++, which is also compared with the estimates of Algorithm 2.1, depicted with red points.
It is important to note two facts: • The difference in computation time: while the first implementation requires more than 100 seconds, the FreeFem++ program obtains estimates for both the mean τ and the standard deviation τ (2) − (τ ) 2 in just 7 seconds. • Figure 6 highlights that, for I(0) small, τ (2) − (τ ) 2 is greater than τ , i.e., predicting accurately the disappearing time of the infection is very difficult for those values of I(0).
In Figure 8, we show the numerical solution to (21)-(22), using FreeFem++. As commented in the previous section, this procedure can be extended to higherorder moments τ (p) (s, y) in the same way as in (12). In Figure 9, we plot their standard deviations. For each case, we have used third-order FEM, with a mesh discretization having 3400 triangles; the CPU times are between 7 and 8 seconds. Figure 9 highlights the standard deviation sizes for a very small number of initially infected individuals, in the three cases.
By analyzing the three cases in Figure 8, from left to right, we can observe that, when R takes its smallest value (left case), the stochastic model approaches the disease-free state. Later, as R continues growing (middle and right cases), a ridge near I(0) = 0 develops. This ridge is a hindrance for the infection to disappear; therefore, the population collapses, because N → 0, and an endemic state is not reached.  Figure 9. Numerical estimate of the standard deviation for the three previous cases, with R = 5/6, 1 and 5/4 respectively. 4. Conclusions. In this paper, we have considered two different types of stochastic SIS epidemic models: the first one, with a constant population size, and the second one, with deaths. In both cases, there is a basic reproduction number R, such that (i) If R ≤ 1, the two models, deterministic and stochastic, approach the diseasefree equilibrium. Moreover, we can estimate the mean and standard deviation times, in order that the infection disappears. Moreover, from Figures 2 and  3, the time to disappearance of the infection depends essentially on I(0). (ii) If R > 1, the asymptotic behavior of both models is quite different. While the deterministic model approaches a unique positive endemic equilibrium, the stochastic model always tends to a disease-free state, although the time needed for the infection to disappear is much greater than in the previous case, and increases with the coefficient R. Moreover, from Figures 4 and 5, we can deduce that the time to disappearance of the infection depends essentially on the size of the population N . (iii) After the introduction of the death of infected individuals into the model, we have found that, although the differential equation is more complicated than in the previous case, the model has nonetheless a very similar asymptotic behavior from a qualitative point of view. Moreover, from Figures 8 and 9, it seems that, for R large, the population collapses. A minor issue of these techniques is that we have to choose a large domain, because the Neumann boundary conditions do not seem completely natural. Indeed, as we can see in the right-hand side of Figure 4, the stationary backward Kolmogorov equation works well for small I(0), whereas for large I(0), the enforcement of homogeneous boundary conditions makes the numerical solution diverge from the results of Algorithm 2.1. In general, this type of problems will arise when we truncate an unbounded domain to a bounded one; in order to avoid them, a promising idea may be the use of transparent boundary conditions (TBCs). This quite recent technique consists in obtaining a solution over a truncated domain that coincides with the general solution, by imposing boundary conditions specifically designed for each individual problem. A very complete description of this technique can be found, for instance, in [6].
Trying to apply TBCs to our problems goes beyond the purpose of this paper. Furthermore, finding accurate transparent boundary conditions can be difficult. However, in our problems, the numerical solutions appear to converge quickly to a plane, as shown in the left-hand side of Figure 4. Therefore, a good approximation of the transparent boundary conditions seems to be where c 1 and c 2 are some constants that need to be chosen carefully. Obviously, the non-homogeneous Neumann boundary conditions need to be introduced into the weak formulation (14). In Figure 10 we replot Figure 5, incorporating the solutions for non-homogeneous Neumann boundary conditions, taking c 1 = c 2 = 0.1 (black) and c 1 = c 2 = 0.25 (magenta). The improvement with respect to c 1 = c 2 = 0 (blue) is clearly visible.
To finish this paper, let us recall an important conclusion of [9]: the numerical solution for the backward Kolmogorov equation can provide valuable information for the epidemic model. Indeed, with a single simulation that takes just a few seconds, we are able to obtain a valuable estimation about the mean persistence time of the infected individuals in a population.