Numerical simulation of a SIS epidemic model based on a nonlinear Volterra integral equation

We consider a SIS epidemic model based on a Volterra integral equation and we compare the dynamical behavior of the analytical solution and its numerical approximation obtained by direct quadrature methods. We prove that, under suitable assumptions, the numerical scheme preserves the qualitative properties of the continuous equation and we show that, as the stepsize tends to zero, the numerical bifurcation points tend to the continuous ones.


1.
Introduction. Volterra Integral Equations (VIEs) appear in fields as automatic control theory, network theory, dynamics of nuclear reactors and in a variety of applied problems [1,6,14,19]. A major source of Volterra equations is population dynamics [4, 7-9, 15, 17, 26]. The numerical simulation of these problems provides an useful support to analytical studies and an essential tool to analyze situations that are difficult to treat analytically.
Here, we are interested in finding approximate solutions to the following nonlinear VIE which represents a disease transmission model, where I(t) is the fraction of infective. The nonlinearity is defined as f (I) = 1 λ0 λ(I)(1 − I), where λ(I) is the contact rate, λ 0 = λ(0) if λ(0) > 0 and λ 0 = 1 if λ(0) > 0, and R 0 is the reproduction number. This model, due to van den Driessche-Watmough [26], has received great attention because of the presence of the nonlinear contact rate λ(I) which provides a very rich dynamical behavior (multiple stable equilibria, backward bifurcation and hysteresis). It is indeed this dynamical behavior that we want to reproduce in the numerical simulation.
In order to guarantee the reliability of the numerical approximation one has to carry out a comparative analysis of the qualitative properties and the long-time behavior of the numerical and analytical solution, when the stepsize h of the discretization is a realistic fixed constant different from zero. In the literature (see [6] and the bibliography therein) we find a well-consolidated theory on the numerical stability for VIEs which studies the behavior of the numerical solution on very simple test problems like and y(t) = y 0 + t 0 (λ + µ(t − s))y(s)ds, t ≥ 0.
More recently, classes of more significative linear convolution equations (see, for example [2] and [21]) have been considered as test equations for numerical purposes, whereas the stability theory on nonlinear problems is much less developed and much more tuned to the characteristics of the problem itself [11,13,22,23]. It turns out that, for the specific equation (1), the numerical stability theory on nonlinear problems developed in recent years do not apply and, furthermore, classical stability results on linearized equations may not take into account important changes in the behavior of the solution. These considerations motivated us to focus our investigations on the performances of numerical methods applied to solve the model (1). This work has been inspired by [10], where the preservation of qualitative properties of the analytical solution in the numerical approximation is analyzed for an integro-differential equation with fading memory.The discrete model obtained upon the application of the Direct Quadrature (DQ) methods with convolution weights to equation (1) reads where P n = P (nh), I n ≈ I(nh), h = const., w n,j and ω n are the weights and are bounded. Here we assume that I 0 (t 0 ), I 1 , . . . , I m−1 are given starting values. For the convergence and the order of the method we refer to [6,21]. In Section 2 we will outline the very well known results on the continuous model and, in Section 3 we carry out a corresponding analysis for the discrete scheme. Of course we expect that all the dynamics of the numerical solution depend on the stepsize h of the discretization and tend to the behavior of the continuous solution as h → 0. The numerical experiments in Section 4 will show this behavior. Some concluding remarks are reported in Section 5.
2. Preliminaries. As described in [26], in equation (1) the convolution kernel P (t) = 1 τ p(t)e −bt is related to the fraction of individuals p(t) infected at time zero and that survive at time t, here b > 0 is the constant birth and death rate and τ = Hence, P (t) is a non-negative and non-increasing function for all t ≥ 0 and +∞ 0 P (s)ds = 1.
In [26] it is shown that, in these assumptions, equation (1) has a unique solution for each I 0 (t) and this solution remains on the interval [0, 1] for all t ≥ 0. Furthermore, equation (1) depends on the value of the single parameter R 0 , its equilibria are the solutions of and displays the following dynamical behavior: i: always has the disease free equilibrium ii: has no endemic equilibrium if R 0 < R c 0 iii: has at least one endemic equilibrium if R 0 > R c As for the local stability, the following result is proved in [26,Th.3]: 3. Numerical equilibria and asymptotic stability. From now on we assume that hypotheses h 0 , h 1 , h 2 , listed in Section 2, hold. These assumptions are the same of assumptions (A1), (A2), (A3) reported in [26,Sec.2] for the continuous model (1) and represent the starting point for our numerical investigation. Since in these assumptions the solution I n of the numerical scheme (3) converges to I(t) as h → 0, we are allowed to choose the stepsize h small enough to have 0 ≤ I n ≤ 1, for all n ≥ 0. The search for equilibrium solutions in the numerical model (3) passes through the following considerations: • since ω n are the bounded weights of a convolution quadrature formula and P (t) is a non-increasing function, we are allowed to write where ϕ(h) is the quadrature error and ϕ(h) → 0 as h → 0; • since lim t→+∞ P (t) = 0, lim t→+∞ I 0 (t) = 0 and the starting weignts w n,j are bounded, then lim n→+∞ I 0n = 0, where I 0n is defined in (4), • since 0 ≤ I n ≤ 1 and λ(I n ) is bounded, then I n f (I n ) is bounded for all n ≥ 0. In this situation, we have that the asymptotic behavior of I n is described by that is the discrete equivalent of the result proved in [20, Th.1] for a nonlinear VIE. Hence, if I n is the solution of (3) with lim n→+∞ I n =Ī h , and 0 ≤ I n ≤ 1, for all n ≥ 0, then, taking into account that lim n→+∞ I 0n = 0,Ī h must satisfy the equation Thus,Ī h is an equilibrium solution of the model (3) if and only ifĪ h is a solution of (8). Let I h = 0 be the disease free equilibrium andĪ h =Ī e h > 0 the endemic equilibrium, the fraction of infective individuals at an endemic equilibriumĪ e h is given implicitly as a function of R 0 by the relation R 0 f (Ī e h ) (1 − ϕ(h)) = 1. We are now in the position to state the following theorem which is the discrete analogue of [26,Th.2] outlined at the beginning of Section 2. ϕ(h) and tend to the continuous ones as h tends to zero.
In order to study the local stability of these bifurcation points, we linearize equation (3) around the equilibrium solutionĪ h with Let r n be the resolvent of the kernel RĪ h hω n P n (for the theory on the resolvent kernel of discrete Volterra equations we refer to [16], [27] and [28]). By the variation of constant formula [11,18,25], it turns out that equation (3) is equivalent to and the asymptotic behavior of I n −Ī h is described by the solution x n of (9) given by x n = I 0 (nh) + n j=0 r n−j I 0 (jh), with I 0 (nh) → 0 as n → ∞.
From the final value theorem for the Z−transform, since (z − 1)x(z) has no poles outside the circle |z| > 1: lim If hP (1) = RĪ h (1 − ϕ(h)) > 1, then (11) has a pole for |z| > 1 and the series x n does not converge.
Assertions i and ii are consequence of the fact that, for the disease free equilibrium, RĪ h = R 0 if λ(0) > 0 and RĪ h = 0 if λ(0) = 0. In order to obtain assertion iii, observe that Taking into account equation (7), which describes the asymptotic behavior of the solutions of the nonlinear equation (3), it is possible to prove the following theorem on the global stability of equilibria, which is the discrete analogue of theorems 5 and 6 in [26]. • the disease free equilibrium is globally asymptotically stable for R 0 (1 − ϕ(h)) < R c 0 , • if λ(0) > 0, for R 0 (1 − ϕ(h)) > R m 0 the unique endemic equilibrium is globally asymptotically stable.

Numerical experiments.
As an example we consider equation (1) and its discrete counterpart (3) with f (I) = (1 − I)(1 + 5I). (12) In this case the threshold parameters are R c 0 = 5 9 and R m 0 = 1 and the equilibria situation is shown in Figure 2 where the number of the equilibria, both for the continuous and the discrete model, are given by the intersections of the lines with the graph of f (I). The nature of the numerical equilibria is shown in Figure 3. By using the results in theorems 3.1-3.3, it is clear from the figure that: for R 0 (1 − ϕ(h)) < R c 0 , the numerical model (3), with f given in (12), has only the disease free equilibrium which is also globally asymptotically stable; for R c 0 < R 0 (1 − ϕ(h)) < R m 0 , the disease free equilibrium is locally stable and two endemic equilibria appear, one stable, corresponding to the branch where f (I e h ) (1 − ϕ(h)) < 0, and the other one unstable; for R 0 (1 − ϕ(h)) > R m 0 , the disease free equilibrium becomes unstable and the unique endemic equilibrium is globally asymptotically stable. The bifurcation diagram shown in Figure 3 as a function of R 0 (1 − ϕ(h)) is the same of the one in [26, Fig.2c] for the continuous model as a function of R 0 . Here, all the dynamic is shifted by a factor 1 − ϕ(h) and tends to that of the continuous model as the stepsize h of the discretization tends to zero. The incidence of the factor 1 − ϕ(h) on the dynamic depends on the numerical method chosen within the class of DQ methods with convolution weights. For this example we use the Trapezoidal DQ method of order 2 and we integrate, with stepsize h = 0.1, equation (1) with P (t) = 1 2 e − 1 2 t , f (I) given in (12) and I 0 (t) → 0. Figure 4 shows the results of the numerical integration for three values of R 0 in different positions with respect to the bifurcation points. It turns out that the numerical solution behaves accordingly to the analytical one showing the same asymptotic behavior. The good performance of the numerical method was to be expected since R 0 , and hence R 0 (1 − ϕ(h)), are far from the bifurcation points. In Figure 5a the Trapezoidal DQ method is used to solve the same problem with R 0 smaller than, but very close to R c 0 (R 0 = R c 0 − , with = 10 −4 ). According to Theorem 5 in [26], with such a value of R 0 the only equilibrium is the disease free one and this is globally asymptotically stable. However, the numerical solution tends to the endemic stable equilibrium that arises when R 0 > R c 0 . This is clearly an example where the quadrature error ϕ(h) plays a crucial role since it affects the position of the numerical parameter with respect to the bifurcation point R c 0 . We expect that reducing h or integrating by an higher order method ϕ(h) becomes smaller and the numerical solution displays the correct behavior. As a matter of fact, Figure  5b shows the result of the numerical integration with h = 0.01 and in Figure 5c an order 4 Gregory DQ method is used. In both cases the numerical solution behaves, once again, with the expected dynamic.
5. Concluding remarks. In this paper we have performed the numerical simulation of the SIS epidemic model with nonlinear contact rate described in [26]. The qualitative behavior of the numerical solution has been analyzed related to the corresponding features of the continuous model. The investigations we have carried out here show that all the dynamics of the numerical model depend on the stepsize h of the discretization through the quadrature error ϕ(h). As we expected, the bifurcation points and the dynamical behavior of the discrete scheme approach the ones of the continuous problem as h → 0. However, when h is positive and constant, the quadrature error ϕ(h) may influence the numerical behavior when the bifurcation parameter R 0 is close to the bifurcation points. In particular, the experiments in Section 4 show that the trapezoidal rule, which is A-stable (that is unconditionally stable with respect to the linear basic test equation (2) (see [3])), may show instabilities when more complicated dynamics are taken into account.
The results obtained here motivates us to extend our studies to other classes of methods and to approach numerically more complicated epidemic models.      (1), with f (I) given in (12).