Global Stability of an SIR Epidemic Model with Delay and General Nonlinear Incidence

An SIR model with distributed delay and a general incidence function is studied. Conditions are given under which the system exhibits threshold behaviour: the disease-free equilibrium is globally asymptotically stable if R0 < 1 and globally attracting if R0 = 1; if R0 > 1, then the unique endemic equilibrium is globally asymptotically stable. The global stability proofs use a Lyapunov functional and do not require uniform persistence to be shown a priori. It is shown that the given conditions are satisfied by several common forms of the incidence function.


1.
Introduction. The prevalence of disease in a population is often described by an SIR model where the population is subdivided into three classes: susceptibles, infecteds and recovereds (or removeds). The simplest forms of these models are ordinary differential equations (ODEs) [10,11]. In [4], a discrete delay model is given to account for transmission by vectors (e.g. mosquitoes), where the delay τ is used to account for a latent period in the vector. Allowing the vectors' latency periods to vary according to some distribution gives a model with a distributed delay [23].
The delay appears in the incidence term which is typically the only nonlinearity, and is therefore the "cause" of all "interesting behaviour". Various forms have been used for the incidence term, both for ODEs and for delay equations. Common forms include mass action βSI [2,18,23], saturating incidence βS I 1+cI [3,24], and standard (or proportional) incidence β SI N [10]. Changing the form of the incidence function can potentially change the behaviour of the system.
In [14] a system of ODEs with a general incidence term f (S, I) is studied. Conditions are found on f under which the standard threshold behaviour occurs: the disease-free equilibrium is globally asymptotically stable for R 0 < 1 and the endemic equilibrium is globally asymptotically stable for R 0 > 1.
The goal of this paper is to present a similar analysis for equations with a bounded distributed delay and a general nonlinear incidence function. The conditions given here are equivalent to those given in [14] for the ODE case. In Section 7, the conditions are shown to apply to mass action, saturating incidence and, for an SI model, standard incidence.
The case of separable incidence, where f (S, I) = a(S)b(I) is studied in [13] for discrete delay. The authors study an SIR system where the delay is included in b(I), modelling vector transmission, and also an SEIR system where the delay appears in both a(S) and b(I), modelling a fixed duration of latency. The work in this paper extends their SIR analysis to distributed delay, while also allowing for a more general class of incidence functions. Example 3 in Section 7 provides the distributed delay version of their SIR analysis.
The approach here is to use a Lyapunov functional of the type used by Goh for ODE models in ecology [7]. Earlier work [17,19,20,21,22] which used a similar Lyapunov functional relied on knowing a priori that the system was uniformly persistent. This is important because the Lyapunov functional is not defined if any of the state variables are zero.
In this paper, we note that if the delay is bounded, then solutions for which the disease is initially present will move to the interior of the state space and so the Lyapunov function becomes defined. This approach greatly simplifies the work (given that the delay is bounded) since it does not require uniform persistence to be shown; rather, uniform persistence follows from the global stability. Issues related to infinite delay are mentioned briefly in a remark at the end of Section 6.
The paper is organized as follows. The model is described in Section 2. In Section 3, the basic reproduction number R 0 is determined and the equilibria are found. Local stability of the equilibria is studied in Section 4. The global dynamics are resolved for R 0 ≤ 1 in Section 5 and for R 0 > 1 in Section 6. In Section 7, examples are given of incidence functions that satisfy the assumptions that are used throughout the paper. 2. The model. A population is divided into susceptible, infectious and recovered (or removed) classes with sizes S, I and R, respectively. Recruitment of new individuals is into the susceptible class, at constant rate Λ. The death rates for the classes are µ S , µ I and µ R , respectively. The average time spent in class I before recovery (or removal) is 1/γ. Thus, the total exit rate for infectives is µ I + γ, which, for biological reasons we assume is at least as large as µ S ; that is, µ I + γ ≥ µ S . Transmission of the disease is through vectors which undergo fast dynamics. Following [4] and [23], the vectors can be omitted from the equations by including a distributed delay τ in the incidence term up to a maximum delay h > 0. The incidence at time t is β h 0 k(τ )f (S(t), I(t − τ ))dτ where k is a Lebesgue integrable function which gives the relative infectivity of vectors of different infection ages. We choose β so that h 0 k(τ )dτ = 1. It is assumed that the support of k has positive measure in any open interval having supremum h so that the interval of integration is not artificially extended by concluding with an interval for which the integral is automatically zero.
The form of the function f is of fundamental importance. In this paper we want to work with a function that is as general as possible, but still possesses the properties necessary for conclusions to be made through mathematical analysis. Because of this, we will introduce conditions on f which may appear technical. However, as shown in Section 7, many commonly used incidence functions satisfy these conditions. For now, we assume only the following.
(H1) f is a non-negative differentiable function on the non-negative quadrant. Furthermore, f is positive if and only if both arguments are positive.
The partial derivatives of f are denoted by f 1 and f 2 . In Sections 4, 5 and 6, it will be shown how extra conditions on f imply the local and global stability of either a disease-free equilibrium or an endemic equilibrium.
In order to avoid excessive use of parentheses in some of the later calcuations, we use the notation S = S(t), I = I(t) and I τ = I(t − τ ). The model equations are and Since R does not appear in the equations for dS dt and dI dt , it is sufficient to analyze the behaviour of solutions to (1). The initial condition for (1) is . Standard theory of functional differential equations [8] can be used to show that solutions of (1) exist and are differentiable for all t > 0. Furthermore, the state space R ≥0 ×C is positively invariant.
Note that d dt (S + I) = Λ − µ S S − (µ I + γ)I ≤ Λ − µ S (S + I), and so lim sup It follows that the system is point dissipative. Without loss of generality, we assume that S(t) + I(t) ≤ 2Λ µS for all t ≥ −h. A consequence of this is that we may assume I is bounded above, which in turn implies dS dt is positive for small S, and so S is positive for t > 0.
We say that disease is initially present if the initial condition satisfies I(θ 0 ) > 0 for some θ 0 ∈ [−h, 0]; since the initial condition is continuous, this means that I is positive on some interval about θ 0 . Then, either I(t) is positive for some t ∈ [0, θ 0 + h] or dI dt (θ 0 + h) > 0, and so I becomes positive. In either case, there exists t 1 ≥ 0 such that I(t 1 ) > 0. Then for t ∈ [t 1 , t 1 +h], we have dI dt ≥ −(µ I +γ)I(t) and so I is strictly positive on this interval. Furthermore, I(t) will remain bounded below by the exponential, Therefore, without loss of generality, we assume that any initial condition for which the disease is initially present satisfies I(θ) > 0 for all θ ∈ [−h, 0]. Furthermore, we assume that for any t ≥ 0, we have I(t + θ) bounded away from 0 for all θ ∈ [−h, 0]. This is useful in Section 6 when evaluating a Lyapunov functional along solutions.
3. Equilibria and R 0 . For any values of the parameters, the disease-free equilibrium is given by The basic reproduction number [5] for the model is The presence and number of endemic equilibria depend on the form of the nonlinearity f , as well as the values of the parameters. In searching for equilibria, we note that the equilibria of Equation (1) are the same as the equilibria of the corresponding ordinary differential equation system. Sufficient conditions for the existence of an endemic equilibrium are given in [6,14,15]. Here, we give the following result. Proof. We look for solutions (S * , I * ) of the equations dS dt = 0, dI dt = 0. We first note that dS dt + dI dt = 0 implies Λ − µ S S * − (µ I + γ)I * = 0, and so S * = Λ−(µI +γ)I *
The function H is continuous and so a sufficient condition for H to have a zero in 0, Λ µI +γ , is that H is increasing at 0. Thus, there is an endemic equilibrium if Since f (S, 0) = 0 for all S, it follows that f 1 (E 0 ) = 0 and so (3) is equivalent to R 0 > 1.

4.
Local stability of the equilibria.
Proof. We begin by linearizing Equation (1) at E 0 . In doing so, we note that f (S, 0) = 0 for all S and so f 1 (E 0 ) = 0. The linearization is Substituting the Ansatz (s(t), i(t)) = e λt (s 0 , i 0 ) into (4) gives Cancelling e λt from each term and rearranging gives the homogeneous linear equation .
There exist non-zero solutions if and only if det(A 0 ) is zero. Thus, the characteristic equation is We show that all solutions λ have negative real part. Suppose λ has non-negative real part. Then, λ + µ S = 0. Also, and so λ cannot be a solution of (5). Hence, all characteristic roots have negative real part and therefore E 0 is locally asymptotically stable [16, Chapter 2, Theorem 4.2].
We now give conditions on f that are used here to show that an endemic equilibrium is locally asymptotically stable, and in Section 6 to show that it is globally asymptotically stable. As a precondition, we assume that R 0 > 1, guaranteeing the existence of an endemic equilibrium E * = (S * , I * ) (see Theorem 3.1). (H4) Either f 1 (S * , I * ) > 0 or f 2 (S * , I * ) < f (S * ,I * ) In order to appreciate that the hypothesis (H4) is not very restrictive, we consider (H2) in a neighbourhood of S = S * , deducing Similarly, considering (H3) at S = S * and in a neighbourhood of I = I * , we deduce We can now see that (H4) is merely requiring that at least one of (H2) and (H3) leads to a strict inequality. On the other hand if (H4) fails to be satisfied, then the endemic equilibrium is still globally attractive (see Section 6) but is not locally asymptotically stable, as the characteristic equation of the linearization at E * will have λ = 0 as a root. Proof. The linearization of Equation (1) at an endemic equilibrium E * = (S * , I * ) is We demonstrate that all zeros of the characteristic equation have negative real part. In order to find the characteristic equation, we substitute the Ansatz (s(t), i(t)) = e λt (s 0 , i 0 ) into (8) to get Cancelling e λt from each term and rearranging gives the homogeneous linear equation .

There exist non-zero solutions if and only if det(A) is zero. Thus, the characteristic equation is
Since (H2) and (H3) hold, the inequalities (6) and (7) are satisfied. Using the equilibrium equation to replace f (E * ) in (7), it follows that Suppose λ is a solution of (10) with non-negative real part. Then, using (6) and (11), we deduce and so the characteristic equation (10) has solutions λ with non-negative real part only if all of the inequalities in (12) are in fact equality. The final inequality is strict unless f 1 (E * ) = 0 (and λ = 0). The second last inequality is strict unless f 2 (E * ) = f (E * )/I * . Assumption (H4) implies at least one is strict and therefore solutions to (10) must have negative real part. Thus, the endemic equilibrium E * is locally asymptotically stable.

5.
Global asymptotic stability for R 0 ≤ 1. The expression for R 0 given in Equation (2) depends on the behaviour of f near the disease-free equilibrium E 0 , which is locally asymptotically stable for R 0 < 1. Results on the global dynamics for R 0 less than one will necessarily require further assumptions on the form of f . Theorem 5.1. Suppose (H5) and (H6) hold. If R 0 < 1, then the disease-free equilibrium E 0 is globally asymptotically stable. If R 0 = 1, and one of (H7.1) and (H7.2) holds, then E 0 is globally attracting.
Proof. We begin by defining and differentiating the function U + (t), which will be one of the terms involved in a Lyapunov functional U . Let Note that ν(τ ) > 0 for 0 ≤ τ < h since the support of k has positive measure near h, and therefore, I ≥ 0 implies U + (t) ≥ 0 with equality if and only if I is identically zero on the interval [t − h, t].
We now find the time-derivative of U + .
Using integration by parts, we obtain From (13), it follows that ν(h) = 0 and dν dτ = −βk(τ ). Using these, as well as the expression for ν(0), we find Next, recall that (H1) requires that f (S, I) be positive if S and I are both positive. Combined with (H6), this implies f 2 (S, 0) > 0 for S > 0, which allows us to make the following definition without fear of division by zero. Let Then dG dx = 1 − f2(S0,0) f2(x,0) which (H5) implies changes from non-positive to nonnegative as x increases through S 0 . Thus, G is minimized at S 0 with G(S 0 ) = 0. Thus, G(x) ≥ 0 for all x > 0. Let Then U (t) is non-negative for S > 0 and I ≥ 0. Using (14) and Λ = µ S S 0 , we obtain Recalling that lim sup(S + I) ≤ Λ µS , we let A 0 = {(S, I) ∈ (0, Λ µS ] × C : dU dt = 0} and let M 0 be the largest invariant subset of A 0 . If dU dt is non-positive, then the Lyapunov-LaSalle Theorem [8, Theorem 5.3.1] implies every omega limit point is contained in M 0 . Case 1: R 0 < 1. By (H5) and (H6), we deduce with equality only if I(t) = 0. So, I is zero at each point in A 0 . Within the set A 0 , we have dS dt = Λ − µ S S, and so M 0 consists of just the point (S 0 , 0). Thus, E 0 is globally attracting. By Theorem 4.1, E 0 is locally asymptotically stable, and so we may conclude that it is globally asymptotically stable.
Recalling that the support of k has positive measure near h, that solutions are continuous and that f is positive when its arguments are positive, we must have I(t − τ ) = 0 for τ near h. Since dS dt = 0 must hold for all t, this implies I is identically zero in M 0 . Thus, M 0 consists of only the point (S 0 , 0). Thus, E 0 is globally attracting. 6. Global asymptotic stability for R 0 > 1. In this section, we resolve the global dynamics for R 0 > 1, given that certain assumptions on f are satisfied. We recall that Theorem 3.1 implies an endemic equilibrium E * exists if R 0 > 1.
Theorem 6.1. Suppose R 0 > 1, and (H2) and (H3) hold. Then the endemic equilibrium E * is unique and all solutions for which the disease is initially present tend to E * . If (H4) also holds, then E * is globally asymptotically stable.
Proof. The uniqueness of E * will follow from the fact that it is globally attracting. We now work towards demonstrating the attractivity of E * . Evaluating both sides of (1) at E * gives Λ = µ S S * + β h 0 k(τ )f (S * , I * )dτ (15) and which will be used as substitutions in the calculations below. Let where We note that ( This determines I to be a constant, and in fact gives I = I * for all t. Thus, each element of M satisfies S(t) = S * and I(t) = I * for all t. We may now conclude that lim t→∞ (S(t), I(t)) = (S * , I * ) = E * . By Theorem 4.2, if (H4) also holds, then E * is locally asymptotically stable, and so it now follows that it is globally asymptotically stable. (1) is modified to have infinite delay, then the basic Lyapunov calculation still works, as long as the delay kernal k is bounded above by a decaying exponential function and the phase space is chosen to be an appropriate fading memory space [1,9,12]. However, it becomes necessary to prove uniform persistence. Even then, since initial conditions could involve I(·) being zero on sets of positive measure, the Lyapunov functional would not be defined for these initial conditions. Since the delay is infinite, the problem would persist for all time. Thus, it becomes necessary to do the Lyapunov calculation for solutions lying in the omega limit sets (or attractor), which by uniform persistence are bounded away from zero. This would show that solutions in the attractor limit to the endemic equilibrium E * . Then, one argues that other solutions must also limit to E * . See [20] for an example of this approach. 7. Examples. We now give examples of incidence functions for which the required hypotheses are satisfied. Example 1: Mass Action Let f (S, I) = SI. Then hypotheses (H1)-(H6) and (H7.1) are satisfied and so the global dynamics are determined by the magnitude of R 0 . The global behaviour of this model was previously studied in [2,18,23] and was fully resolved in [21].

Remark 1. If Equation
Example 2: Saturating Incidence Let f (S, I) = S I 1+cI for some constant c > 0. Then hypotheses (H1)-(H6), as well as both of (H7.1) and (H7.2), are satisfied and so the global dynamics are determined by the magnitude of R 0 . The discrete delay version of this model was previously studied in [24] with the global dynamics being resolved in [22]. The discrete delay version of separable incidence is studied in [13].
Example 4: Standard Incidence For this example, it is necessary to interpret the R class as removed. Then, the total population is S + I. In this case, standard incidence is given by using f (S, I) = SI S+I . Hypothesis (H1) is not satisfied since f is not defined at (0, 0). However, since {S = 0} is repelling, we may relax (H1), requiring only that f ∈ C 1 (R >0 × R ≥0 → R ≥0 ), with f (S, I) = 0 if and only if I = 0. This condition is satisfied, as are hypotheses (H2)-(H6) and (H7.2) and so the global dynamics are determined by the magnitude of R 0 .