- Email: [email protected]
On the Lyapunov stability for SIRS epidemic models with general nonlinear incidence rate
Applied Mathematics and Computation 217 (2010) 4010–4016
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
On the Lyapunov stability for SIRS epidemic models with general nonlinear incidence rate Bruno Buonomo ⇑, Salvatore Rionero Department of Mathematics and Applications, University of Naples Federico II, via Cintia, I-80126 Naples, Italy
a r t i c l e
i n f o
Keywords: Epidemic model Lyapunov direct method Lyapunov function Nonlinear incidence rate
a b s t r a c t We study an epidemic model for infections with non permanent acquired immunity (SIRS). The incidence rate is assumed to be a general nonlinear function of the susceptibles and the infectious classes. By using a peculiar Lyapunov function, we obtain necessary and sufficient conditions for the local nonlinear stability of equilibria. Conditions ensuring the global stability are also obtained. Unlike the recent literature on this subject, here no restrictions are required about the monotonicity and concavity of the incidence rate with respect to the infectious class. Among the applications, the noteworthy case of a convex incidence rate is provided. Ó 2010 Elsevier Inc. All rights reserved.
1. Introduction In the seventies, Capasso and his coworkers stressed the importance to consider nonlinear incidence rates for some disease spreads, [3,5,6]. The case study was the investigation on a cholera epidemic spread. Since then, many scholars have proposed various nonlinear forms of the incidence rate (see e.g., the brief surveys contained in [7,20]). However, as it has been underlined in [10], the details of transmission of infectious diseases are generally unknown, and may depend on several factors. For this reason there is a wide and growing interest in studying epidemic models with general incidence rates. Examples are given in [10,12,13]. In this paper, we study an epidemic model for infections with non permanent acquired immunity (SIRS). We consider a general nonlinear incidence rate depending on susceptibles (S) and infectious (I) classes, say f(S, I). In particular, the stability of the equilibria states is performed. We introduce a Lyapunov function, which is the ODEs ‘adaptation’ of a Lyapunov function introduced by Rionero in the context of L2-stability analysis for binary reaction–diffusion systems, see e.g. [14–18]. This function allows us to obtain the following results: (i) Exact coincidence between the linear and nonlinear stability of equilibria (at least locally). (ii) The stability result requires the function f(S, I) to be increasing respect to the variable S and does not require any restriction about the monotonicity and the concavity of the function f, with respect to the variable I. Recently, the same model has been considered by Korobeinikov, [9,10]. By using the Lyapunov direct method, he found sufficient conditions on f(S, I) ensuring that there exists a unique endemic equilibrium and it is globally asymptotically stable. Such conditions are satisfied by several well known incidence rates as, e.g., the ones reflecting the effect of saturation with respect to the infectious I. However, the sufficient conditions found by Korobeinikov rule out some other incidence rates, as the ones considered, e.g. in [19,20], which are convex (i.e. concave up) respect to the variable I. Here, among the applications,
⇑ Corresponding author. E-mail addresses: [email protected] (B. Buonomo), [email protected] (S. Rionero). 0096-3003/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2010.10.007
B. Buonomo, S. Rionero / Applied Mathematics and Computation 217 (2010) 4010–4016
4011
the case of a convex incidence rate is provided. Hence, in view of the results (i) and (ii) above, we enlarge the class of incidence rates ensuring the nonlinear stability. The paper is organized as follows. In Section 2 we introduce the SIRS model. In Section 3 we present our main result, i. e. the nonlinear stability/ instability thresholds. Some applications are given in Section 4. In particular, we consider three cases: the classical bilinear incidence rate; the incidence rate with ‘psichological effect’, which is described by a non monotone (respect to I) functional; a quadratic incidence rate, which is an example of convex (respect to I) incidence rate. For this last case, in Section 5, we enhance the result of the previous sections, by obtaining sufficient conditions for the global stability. Concluding remarks are given in Section 6. 2. The SIRS model We consider the following SIRS epidemic model (for the general theory about epidemic models, see e.g., [1,2,4,8]):
8 _ > < S ¼ l f ðS; IÞ lS þ cR; I_ ¼ f ðS; IÞ ðr þ lÞI; > :_ R ¼ rI ðl þ cÞR;
ð1Þ
where the upper dot denotes the time derivative, d /dt. The state variables are the fractions in which a host population is divided: the susceptibles, S, the infectious, I, and the recovered, R. As in [9], we assume the equality between birth rate and death rate. This assumption allows us a better comparison with the result obtained in [9], although this restriction is unnecessary for our analysis. The parameters (all positive constants) have the following meanings: l is the birth/death rate, r is the recovery rate and c is the rate at which recovered individuals lose immunity and return to the susceptible class. The nonlinear incidence, f, is required to satisfy the following general assumptions:
8 > < f : ðS; IÞ 2 Q ! f ðS; IÞ 2 Rþ ; f 2 C 1 ðQ Þ; > : f ðS; 0Þ ¼ f ð0; IÞ ¼ 0; for all S P 0;
I P 0;
R2þ .
where Q ¼ ½0; 1 ½0; 1 Denoting by N(t) the size of the total population, i.e. setting N = S + I + R, from (1) we get: N_ ¼ l lN, so that for any initial N(t0) P 0 it follows: limt?1N(t) = 1. Hence, the plane S + I + R = 1, is an invariant manifold of system (1), which is attracting in the first octant, [20]. Here, we assume that the population is in equilibrium and investigate the dynamics of system (1) on such plane. Equilibria E(S*, I*, R*) of system (1) satisfy the system:
8 f ðS ; I Þ þ lS ¼ l þ cR ; > > > < f ðS ; I Þ ¼ ðr þ lÞI ; > rI ¼ ðl þ cÞR ; > > : S ¼ 1 I R ;
ð2Þ
Because of f(S, 0) = 0, for all S P 0, there exists the disease-free equilibrium:
E0 ð1; 0; 0Þ: The existence of one or more endemic equilibria (i.e. equilibria with positive components) will depend on the form of the incidence rate, as we will see in Section 4. The existence of the invariant manifold allow us to refer to the reduced binary ODE:
(
I_ ¼ f ðSðR; IÞ; IÞ ðr þ lÞI; R_ ¼ rI ðl þ cÞR;
ð3Þ
where, S(R, I) = 1 I R, and (I(t), R(t)) 2 R, where,
R ¼ fðI; RÞ 2 R2 : I; R P 0; I þ R 6 1g: 3. Nonlinear stability/ instability Let us consider the perturbations to the generic equilibrium, u = S S*, the perturbation system:
_ ¼ ðr þ lÞw þ u½uðv ; wÞ; w; w v_ ¼ rw ðl þ cÞv ;
v = R R* and w = I I*. From (2) and (3) we get ð4Þ
where,
u½uðv ; wÞ; w ¼ f ðuðv ; wÞ þ S ; w þ I Þ f ðS ; I Þ;
ð5Þ
4012
B. Buonomo, S. Rionero / Applied Mathematics and Computation 217 (2010) 4010–4016
and u(v, w) = v w. We expand f into Mac Laurin series in powers of u and ditions for expansion) and obtain:
v (we assume that f satisfies the required con-
_ ¼ ða ðr þ lÞÞw þ bv þ Uðw; v Þ; w v_ ¼ rw ðl þ cÞv ;
ð6Þ
where U includes all the terms involving u and
a¼
@ u
b¼
@ u
@w u¼0;w¼0 @ v u¼0;w¼0
v with powers higher than unity, U(0, 0) = 0, and,
@f @S þ @f ; @I S¼S ;I¼I S¼S ;I¼I @f @S : S¼S ;I¼I
In the following, we will prove a stability/ instability theorem. In order to underline the wide generality of our result, we will refer to the following general system:
v_ ¼ Jv þ Gðv Þ;
ð7Þ
where J = [aij], whose entries are given by:
a11 ¼ a ðr þ lÞ;
a21 ¼ r;
a12 ¼ b;
a22 ¼ ðl þ cÞ;
ð8Þ
T
and G = (g1, g2) (U, 0) . Observe that (7) reduces to system (6) by setting: v (v1, v2) = (w, v)T. Denote the principal invariants of the Jacobian matrix J as: T = a11 + a22, and A = a11a22 a12a21. Now we prove the stability theorem for the general system (7). As corollary, the same result will be true for system (6). qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Theorem 1. Let g1 and g2 be infinitesimal of higher order than v 21 þ v 22 . If T < 0, and A > 0, then the null solution of (7) is nonlinearly (locally) asymptotically exponentially stable. If T > 0 and A > 0, or A < 0, then it is unstable. Proof. The result will be obtained following the discussion given in [18]. We introduce the following Lyapunov function:
V¼
i 1 h 2 T A v 1 þ v 22 þ ða11 v 2 a21 v 1 Þ2 þ ða12 v 2 a22 v 1 Þ2 ; 2
ð9Þ
= 1, for T > 0, and = 1, for T < 0. Along the solution of (7), it turns out:
where
2 _ Vj ð7Þ ¼ AT
v 21 þ v 22
þ Wðv 1 ; v 2 Þ;
ð10Þ
where,
W ¼ T½ða1 v 1 a3 v 2 Þg 1 ðv 1 ; v 2 Þ þ ða2 v 2 a3 v 1 Þg 2 ðv 1 ; v 2 Þ; and,
a1 ¼ A þ a221 þ a222 ; a2 ¼ A þ a211 þ a212 ; a3 ¼ a11 a21 þ a12 a22 :
ð11Þ
In view of the hypothesis, two positive constants m0 and g0 exist, such that:
jg 1 j þ jg 2 j 6 m0
v 21 þ v 22
g0 þ1
;
so that
gþ1 jWj 6 m v 21 þ v 22 ; where m and g are positive constants. Consider now the circle:
Bð0; rÞ ¼
8 < :
2
ðv 1 ; v 2 Þ 2 R : v þ v < 2 1
2 2
jAjT 2 m
!21g 9 = : ;
It follows that: (i) For T < 0, A > 0, = 1, in any circle B the function V is positive definite and V_ is negative definite. Therefore the zero solution of (7) is nonlinearly locally asymptotically stable under the same conditions of linear stability; (ii) For T > 0, A > 0, = 1, both V and V_ are positive definite in any circle B, and hence the zero solution of (7) is (nonlinearly) unstable; (iii) For A < 0, = 1, the function V takes positive values in any circle B and V_ is positive definite there. Then, the zero solution of (7) is (nonlinearly) unstable. h We now specialize the result obtained in Theorem 1 to the SIRS epidemic model case:
B. Buonomo, S. Rionero / Applied Mathematics and Computation 217 (2010) 4010–4016
4013
Corollary 1. If
ðl þ cÞ @f @I
r lþc þ 1
S¼S ;I¼I
@f @S
@f @S
S¼S ;I¼I
S¼S ;I¼I
@f @I
S¼S ;I¼I
þ r þ l > 0; ð12Þ
ðr þ lÞ ðl þ cÞ < 0;
then the generic equilibrium E = (S*, I*, R*) is nonlinearly locally asymptotically stable. If one of the two inequalities (12) holds reversed, then E is nonlinearly unstable.
l + c)[a Proof. The proof follows from Theorem 1 by observing that, from (8), A=p( ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi (r + l)] rb, and T = a r c 2l, and that the boundedness of the solutions of (1) guarantees that U ¼ o w2 þ v 2 (see system (6)). h Set now: R0 ¼
1 @f : ðr þ lÞ @I S¼1;I¼0
ð13Þ
Because of the hypothesis f(S, 0) = 0, it follows that (@ f/@S)S=1,I=0 = 0. Hence, in terms of the equilibrium E0, the inequalities (12) can be written:
ðr þ lÞð1 R0 Þ > 0;
ðr þ lÞðR0 1Þ ðl þ cÞ < 0:
ð14Þ
From inequalities (14) and Theorem 1, we get the following result concerning the stability/instability of the trivial equilibrium E0: if R0 < 1, then E0 is nonlinearly stable. If R0 > 1, then E0 is unstable. 4. Incidence rates In this section we will apply our result by considering several different incidence rates. We will consider the class:
f ðS; IÞ ¼ SgðIÞ:
ð15Þ
(i) Principle of mass action. In the classical SIR model with vital dynamics (see [4]), it is postulated that the spread of the disease among the population occurs according to the principle of mass action. This means that the incidence rate is bilinear respect to S and I. In other words,
gðIÞ ¼ kI;
ð16Þ 1
where k is a positive constant. From (13), (15) and (16) we have R0 = k(r + l) . Then, if R0 > 1, system (1), where the incidence rate f is given by (15) and (16), admits the unique endemic equilibrium E (S*, I*, R*), where:
S ¼ R1 0 ;
R ¼ rl1 I ;
I ¼ lg1 ðR0 1Þ;
and g = k[1 rc(l + r)1(l + c)1] > 0. It is easy to check that the inequalities (12) are satisfied, provided that R0 > 1. As a consequence, Theorem 1 ensures that the equilibrium E, when it exists, is nonlinearly stable. (ii) Incidence rates with ‘psychological’ effect. When the function g(I) in (15) is non monotone, in the sense that it is increasing when I is small and decreasing when I is large, then it may represent the so called ‘psychological’ effect. Large number of infectious may push the population to reduce the number of contacts per unit time, so that the infection force first increases and then decreases as I increases (inhibitory effect). An example is given by the following functional, which has been proposed in [20]:
gðIÞ ¼
kI 1 þ aI2
;
ð17Þ
where k and a are positive constants. In (17), kI measures the infection force of the disease and 1/(1 + a I2) describes the inhibitory effect from the behavioral change of the susceptible individuals when the number of infectious individuals is very large. The functional (17) is non monotone. From (13), (15) and (17) it follows: R0 = k(r + l)1. It is shown in [20] that if R0 > 1, system (1), where the incidence rate f is given by (15) and (17), admits the unique endemic equilibrium E (S*, I*, R*). In [20] it is also shown that such equilibrium is a stable node. In other word, the inequalities (12) are satisfied. Hence, Theorem 1 ensures that the equilibrium E is nonlinearly stable. (iii) A quadratic incidence rate. Consider:
gðIÞ ¼ kIð1 þ aIÞ;
ð18Þ
4014
B. Buonomo, S. Rionero / Applied Mathematics and Computation 217 (2010) 4010–4016
where k and a are positive constants. As stressed in [19], the use of this functional means the inclusion of an increased rate of infection due to two exposures over a short time period. The single contacts lead to infection at the rate kIS, whereas the new infectious individuals arise from double exposures at a rate akI2S. From (13), (15) and (18) we have R0 = k(r + l)1. The dynamics of this model is very rich. Following the same discussion given in [11], it can be seen that for R0 < 1, system (1), where the incidence rate f is given by (15)–(18), may admit none, one, or two endemic equilibria. Backward bifurcation occurs. For R0 > 1, there exists a unique endemic equilibrium E. Stability switches are also possible, so that the stability/instability of the endemic equilibria depends on the fulfillment of additional conditions on the parameters. For example, in the case R0 > 1, it may be shown that the determinant of the Jacobian matrix is always positive, whereas the sign of the trace depends on the parameters. The following result has been stated. Set:
r k 4ð1 þ qÞ2 q¼ ; B ¼ ; B0 ¼ 2 ; k ðl þ cÞ l q ðq þ 2Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð4q þ 4 þ 2 4q2 þ 8q þ 4 þ B2 q2 Þðq þ 1Þ ; p0 ¼ 2B2 q pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi B þ 2q þ 4 þ 2 ðB þ q þ 2Þðq þ 2Þ : p2 ¼ B2
p¼
aðl þ cÞ
;
Proposition 1 [11]. If B 6 B0 and p < p0, or B > B0 and p < p2, then there exist an endemic equilibrium E, which is locally stable. As a consequence, Theorem 1 ensures that the hypotheses of Proposition 1 guarantee also the nonlinear stability of E. 5. Global nonlinear asymptotic stability As stressed in Section 1, in [10] it has been shown that the SIRS model (1) admits a unique endemic equilibrium, which is globally asymptotically stable in R2þ , provided that R0 > 1 and the incidence rate satisfies some sufficient conditions. Precisely, the following inequalities must be satisfied, [10]:
f ðS; I Þ < f ðS ; I Þ for all 0 < S < S ; f ðS; I Þ > f ðS ; I Þ for all S > S ; and,
I=I 6 f ðS; IÞ=f ðS; I Þ 6 1 for 0 < I 6 I ; 1 6 f ðS; IÞ=f ðS; I Þ 6 I=I
ð19Þ
for I P I :
It can be seen that inequalities (19) rule out several incidence rates, as the ones considered, e.g. in [19,20], which are convex (i.e. concave up) respect to the variable I. Here, the use of the Lyapunov direct method, through the peculiar Lyapunov function (9), allows to show that the conditions for the linear stability are necessary and sufficient for the local nonlinear stability of hyperbolic equilibria, according to Theorem 1. This result implies that no restrictions are required about the monotonicity and concavity of the incidence rate with respect to the infectious class. Furthermore, if one is able to ‘control’ the functional W in (10), then sufficient conditions for the global stability may be obtained. In this section, we find such conditions in the case of quadratic incidence rate, (15)– (18), which is convex (concave up) with respect to the variable I. We may start from (4). Observe that (5) is now:
u ¼ kI ð1 þ aI Þu þ kS ð1 þ 2aI Þw þ ðk þ 2akI Þuw þ þakS w2 þ akuw2 ; where u = v w and S* = 1 I* R*. Hence, we obtain system (6), where now:
Uðw; v Þ ¼ ðk þ 2akI Þv w þ ðakS k 2akI Þw2 akv w2 akw3 ;
ð20Þ
and,
a ¼ kS ð1 þ 2aI Þ kI ð1 þ aI Þ;
b ¼ kI ð1 þ aI Þ:
We can now refer to system (7), by setting v (v1, v2) = (w, v)T; J = [aij], whose entries are given by
a11 ¼ a ðr þ lÞ;
a12 ¼ b;
a21 ¼ r;
a22 ¼ ðl þ cÞ;
where a and b are given in (21), and G = (g1, g2) (U, 0)T, where U is given by (20). Reasoning as in the proof of Theorem 1, we get equality (10), where now:
W ¼ a1 T v 1 Uðv 1 ; v 2 Þ;
ð21Þ
B. Buonomo, S. Rionero / Applied Mathematics and Computation 217 (2010) 4010–4016
4015
where T < 0, a1 is given by (11) and it is assumed that a11 – 0. Recall now that the solutions of (1) are bounded, so that we can consider uniformly bounded perturbations. That is, we assume that a positive constant m exists, m 6 2, such that
jv 1 j; jv 2 j < m: In this way, taking into account of (20), we get the estimate:
jWj 6 ma1 jTj
1 2 jk þ 2akI j v 21 þ v 22 þ jakS k 2akI jv 21 þ 2akmv 1 ; 2
that is,
jWj 6 d1
v 21 þ v 22
;
where,
d1 ¼ ma1 jTj½jk þ 2akI j þ jakS k 2akI j þ 2akm: Observe now that, from (9), it follows:
k1
v 21 þ v 22
< V < k2
v 21 þ v 22
;
ð22Þ
where
k1 ¼
AjTj ; 2
k2 ¼ jTj max A; 2 a211 þ a221 ; 2 a212 þ a222 :
ð23Þ
Hence, taking into account of (10), (22) and (23), it follows:
_ Vj ð7Þ 6 ðd3 þ d4 ÞV; where
d3 ¼
AT 2 ; k2
d4 ¼
d1 : k1
ð24Þ
Hence,
V 6 Vð0Þeðd3 d4 Þt : We can summarize the result of this section in the following theorem: Theorem 2. When the SIRS epidemic model (1) with the incidence force given by (18) admits a locally stable equilibrium, then such equilibrium is nonlinearly globally stable, provided that d3 > d4, where d3 and d4 are given by (24).
6. Conclusions In this paper we provide a direct application to epidemic models of a recent approach to nonlinear stability, introduced by Rionero in the context of L2-stability analysis for binary reaction–diffusion systems and here adapted to ODE systems. We find exact coincidence between the linear and nonlinear stability conditions for the endemic equilibrium and enlarge (with respect to what established in [10]) the class of incidence rates ensuring the nonlinear stability of endemic equilibria. Further investigations will include diffusion processes and higher dimensional epidemic systems, as SEIR-like models. Acknowledgements The present work has been performed under the auspices of the italian National Group for Mathematical Physics (GNFMINDAM). References [1] [2] [3] [4] [5]
R.M. Anderson, R.M. May, Infectious Diseases of Humans: Dynamics and Control, Oxford University Press, Oxford, 1991. S.N. Busenberg, K. Cooke, Vertically Transmitted Diseases, Biomathematics, vol. 23, Springer-Verlag, Berlin, 1993. V. Capasso, Global solution for a diffusive nonlinear deterministic epidemic model, SIAM J. Appl. Anal. 35 (1978) 274–284. V. Capasso, Mathematical Structures of Epidemic Systems, Lecture Notes in Biomath, vol. 97, Springer-Verlag, Berlin, 1993. V. Capasso, E. Grosso, G. Serio, I modelli matematici nella indagine epidemiologica Applicazione all’epidemia di colera verificatasi in Bari nel 1973, Annali sclavo 19 (1977) 193–208. [6] V. Capasso, G. Serio, A generalization of the Kermack-Mc Kendrick deterministic epidemic model, Math. Biosci. 42 (1978) 41–61. [7] A. d’Onofrio, Vaccination policies and nonlinear force of infection: generalization of an observation by Alexander and Moghadas (2004), Appl. Math. Comput. 168 (2005) 613–622. [8] H.W. Hethcote, The mathematics of infectious diseases, SIAM Rev. 42 (2000) 599–653.
4016
B. Buonomo, S. Rionero / Applied Mathematics and Computation 217 (2010) 4010–4016
[9] A. Korobeinikov, Lyapunov functions and global stability for SIR and SIRS epidemiological models with non-linear transmission, Bull. Math. Biol. 30 (2006) 615–626. [10] A. Korobeinikov, Global properties of infectious disease models with nonlinear incidence, Bull. Math. Biol. 69 (2007) 1871–1886. [11] Y. Jin, W. Wang, S. Xiao, An SIRS model with a nonlinear incidence rate, Chaos, Solitons Fractals 34 (2007) 1482–1497. [12] G.H. Li, Z. Jin, Global stability of an SEI epidemic model with general contact rate, Chaos, Solitons Fractals 23 (2005) 997–1004. [13] M.Y. Li, J.S. Muldowney, P. van den Driessche, Global stability of SEIRS models in epidemiology, Can. Appl. Math. Quart. 7 (1999) 409–425. [14] S. Rionero, L2 stability of solutions to a nonlinear binary reaction-diffusion system of P.D.E.s, Rend. Math. Acc. Lincei 16 (2005) 227–238. [15] S. Rionero, A nonlinear L2 stability analysis for two species dynamics with dispersal, Math. Biosci. Eng. 3 (2006) 189–204. [16] S. Rionero, A rigorous reduction f the L2-stability of the solutions to a nonlinear binary reaction-diffusion system of P.D.E.s to the stability of the solutions to a linear binary system of ODE’s, J. Math. Anal. Appl. 319 (2006) 377–397. [17] S. Rionero, Nonlinear L2-stability analysis for two-species population dynamics in spatial ecology under Neumann boundary data, Rend. Circ. Math. Palermo. Ser. II 78 (Suppl.) (2006) 273–283. [18] S. Rionero, On the nonlinear stability of the critical points of an epidemic SEIR model via a novel Lyapunov function, Rend. Acc. Sci. Fis. Nat. Napoli 75 (4) (2008) 115–129. [19] P. van den Driessche, J. Watmough, A simple SIS epidemic model with backward bifurcation, J. Math. Biol. 40 (2000) 525–540. [20] D. Xiao, S. Ruan, Global analysis of an epidemic model with nonmonotone incidence rate, Math. Biosci. 208 (2007) 129–419.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
On the Lyapunov stability for SIRS epidemic models with general nonlinear incidence rate Bruno Buonomo ⇑, Salvatore Rionero Department of Mathematics and Applications, University of Naples Federico II, via Cintia, I-80126 Naples, Italy
a r t i c l e
i n f o
Keywords: Epidemic model Lyapunov direct method Lyapunov function Nonlinear incidence rate
a b s t r a c t We study an epidemic model for infections with non permanent acquired immunity (SIRS). The incidence rate is assumed to be a general nonlinear function of the susceptibles and the infectious classes. By using a peculiar Lyapunov function, we obtain necessary and sufficient conditions for the local nonlinear stability of equilibria. Conditions ensuring the global stability are also obtained. Unlike the recent literature on this subject, here no restrictions are required about the monotonicity and concavity of the incidence rate with respect to the infectious class. Among the applications, the noteworthy case of a convex incidence rate is provided. Ó 2010 Elsevier Inc. All rights reserved.
1. Introduction In the seventies, Capasso and his coworkers stressed the importance to consider nonlinear incidence rates for some disease spreads, [3,5,6]. The case study was the investigation on a cholera epidemic spread. Since then, many scholars have proposed various nonlinear forms of the incidence rate (see e.g., the brief surveys contained in [7,20]). However, as it has been underlined in [10], the details of transmission of infectious diseases are generally unknown, and may depend on several factors. For this reason there is a wide and growing interest in studying epidemic models with general incidence rates. Examples are given in [10,12,13]. In this paper, we study an epidemic model for infections with non permanent acquired immunity (SIRS). We consider a general nonlinear incidence rate depending on susceptibles (S) and infectious (I) classes, say f(S, I). In particular, the stability of the equilibria states is performed. We introduce a Lyapunov function, which is the ODEs ‘adaptation’ of a Lyapunov function introduced by Rionero in the context of L2-stability analysis for binary reaction–diffusion systems, see e.g. [14–18]. This function allows us to obtain the following results: (i) Exact coincidence between the linear and nonlinear stability of equilibria (at least locally). (ii) The stability result requires the function f(S, I) to be increasing respect to the variable S and does not require any restriction about the monotonicity and the concavity of the function f, with respect to the variable I. Recently, the same model has been considered by Korobeinikov, [9,10]. By using the Lyapunov direct method, he found sufficient conditions on f(S, I) ensuring that there exists a unique endemic equilibrium and it is globally asymptotically stable. Such conditions are satisfied by several well known incidence rates as, e.g., the ones reflecting the effect of saturation with respect to the infectious I. However, the sufficient conditions found by Korobeinikov rule out some other incidence rates, as the ones considered, e.g. in [19,20], which are convex (i.e. concave up) respect to the variable I. Here, among the applications,
⇑ Corresponding author. E-mail addresses: [email protected] (B. Buonomo), [email protected] (S. Rionero). 0096-3003/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2010.10.007
B. Buonomo, S. Rionero / Applied Mathematics and Computation 217 (2010) 4010–4016
4011
the case of a convex incidence rate is provided. Hence, in view of the results (i) and (ii) above, we enlarge the class of incidence rates ensuring the nonlinear stability. The paper is organized as follows. In Section 2 we introduce the SIRS model. In Section 3 we present our main result, i. e. the nonlinear stability/ instability thresholds. Some applications are given in Section 4. In particular, we consider three cases: the classical bilinear incidence rate; the incidence rate with ‘psichological effect’, which is described by a non monotone (respect to I) functional; a quadratic incidence rate, which is an example of convex (respect to I) incidence rate. For this last case, in Section 5, we enhance the result of the previous sections, by obtaining sufficient conditions for the global stability. Concluding remarks are given in Section 6. 2. The SIRS model We consider the following SIRS epidemic model (for the general theory about epidemic models, see e.g., [1,2,4,8]):
8 _ > < S ¼ l f ðS; IÞ lS þ cR; I_ ¼ f ðS; IÞ ðr þ lÞI; > :_ R ¼ rI ðl þ cÞR;
ð1Þ
where the upper dot denotes the time derivative, d /dt. The state variables are the fractions in which a host population is divided: the susceptibles, S, the infectious, I, and the recovered, R. As in [9], we assume the equality between birth rate and death rate. This assumption allows us a better comparison with the result obtained in [9], although this restriction is unnecessary for our analysis. The parameters (all positive constants) have the following meanings: l is the birth/death rate, r is the recovery rate and c is the rate at which recovered individuals lose immunity and return to the susceptible class. The nonlinear incidence, f, is required to satisfy the following general assumptions:
8 > < f : ðS; IÞ 2 Q ! f ðS; IÞ 2 Rþ ; f 2 C 1 ðQ Þ; > : f ðS; 0Þ ¼ f ð0; IÞ ¼ 0; for all S P 0;
I P 0;
R2þ .
where Q ¼ ½0; 1 ½0; 1 Denoting by N(t) the size of the total population, i.e. setting N = S + I + R, from (1) we get: N_ ¼ l lN, so that for any initial N(t0) P 0 it follows: limt?1N(t) = 1. Hence, the plane S + I + R = 1, is an invariant manifold of system (1), which is attracting in the first octant, [20]. Here, we assume that the population is in equilibrium and investigate the dynamics of system (1) on such plane. Equilibria E(S*, I*, R*) of system (1) satisfy the system:
8 f ðS ; I Þ þ lS ¼ l þ cR ; > > > < f ðS ; I Þ ¼ ðr þ lÞI ; > rI ¼ ðl þ cÞR ; > > : S ¼ 1 I R ;
ð2Þ
Because of f(S, 0) = 0, for all S P 0, there exists the disease-free equilibrium:
E0 ð1; 0; 0Þ: The existence of one or more endemic equilibria (i.e. equilibria with positive components) will depend on the form of the incidence rate, as we will see in Section 4. The existence of the invariant manifold allow us to refer to the reduced binary ODE:
(
I_ ¼ f ðSðR; IÞ; IÞ ðr þ lÞI; R_ ¼ rI ðl þ cÞR;
ð3Þ
where, S(R, I) = 1 I R, and (I(t), R(t)) 2 R, where,
R ¼ fðI; RÞ 2 R2 : I; R P 0; I þ R 6 1g: 3. Nonlinear stability/ instability Let us consider the perturbations to the generic equilibrium, u = S S*, the perturbation system:
_ ¼ ðr þ lÞw þ u½uðv ; wÞ; w; w v_ ¼ rw ðl þ cÞv ;
v = R R* and w = I I*. From (2) and (3) we get ð4Þ
where,
u½uðv ; wÞ; w ¼ f ðuðv ; wÞ þ S ; w þ I Þ f ðS ; I Þ;
ð5Þ
4012
B. Buonomo, S. Rionero / Applied Mathematics and Computation 217 (2010) 4010–4016
and u(v, w) = v w. We expand f into Mac Laurin series in powers of u and ditions for expansion) and obtain:
v (we assume that f satisfies the required con-
_ ¼ ða ðr þ lÞÞw þ bv þ Uðw; v Þ; w v_ ¼ rw ðl þ cÞv ;
ð6Þ
where U includes all the terms involving u and
a¼
@ u
b¼
@ u
@w u¼0;w¼0 @ v u¼0;w¼0
v with powers higher than unity, U(0, 0) = 0, and,
@f @S þ @f ; @I S¼S ;I¼I S¼S ;I¼I @f @S : S¼S ;I¼I
In the following, we will prove a stability/ instability theorem. In order to underline the wide generality of our result, we will refer to the following general system:
v_ ¼ Jv þ Gðv Þ;
ð7Þ
where J = [aij], whose entries are given by:
a11 ¼ a ðr þ lÞ;
a21 ¼ r;
a12 ¼ b;
a22 ¼ ðl þ cÞ;
ð8Þ
T
and G = (g1, g2) (U, 0) . Observe that (7) reduces to system (6) by setting: v (v1, v2) = (w, v)T. Denote the principal invariants of the Jacobian matrix J as: T = a11 + a22, and A = a11a22 a12a21. Now we prove the stability theorem for the general system (7). As corollary, the same result will be true for system (6). qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Theorem 1. Let g1 and g2 be infinitesimal of higher order than v 21 þ v 22 . If T < 0, and A > 0, then the null solution of (7) is nonlinearly (locally) asymptotically exponentially stable. If T > 0 and A > 0, or A < 0, then it is unstable. Proof. The result will be obtained following the discussion given in [18]. We introduce the following Lyapunov function:
V¼
i 1 h 2 T A v 1 þ v 22 þ ða11 v 2 a21 v 1 Þ2 þ ða12 v 2 a22 v 1 Þ2 ; 2
ð9Þ
= 1, for T > 0, and = 1, for T < 0. Along the solution of (7), it turns out:
where
2 _ Vj ð7Þ ¼ AT
v 21 þ v 22
þ Wðv 1 ; v 2 Þ;
ð10Þ
where,
W ¼ T½ða1 v 1 a3 v 2 Þg 1 ðv 1 ; v 2 Þ þ ða2 v 2 a3 v 1 Þg 2 ðv 1 ; v 2 Þ; and,
a1 ¼ A þ a221 þ a222 ; a2 ¼ A þ a211 þ a212 ; a3 ¼ a11 a21 þ a12 a22 :
ð11Þ
In view of the hypothesis, two positive constants m0 and g0 exist, such that:
jg 1 j þ jg 2 j 6 m0
v 21 þ v 22
g0 þ1
;
so that
gþ1 jWj 6 m v 21 þ v 22 ; where m and g are positive constants. Consider now the circle:
Bð0; rÞ ¼
8 < :
2
ðv 1 ; v 2 Þ 2 R : v þ v < 2 1
2 2
jAjT 2 m
!21g 9 = : ;
It follows that: (i) For T < 0, A > 0, = 1, in any circle B the function V is positive definite and V_ is negative definite. Therefore the zero solution of (7) is nonlinearly locally asymptotically stable under the same conditions of linear stability; (ii) For T > 0, A > 0, = 1, both V and V_ are positive definite in any circle B, and hence the zero solution of (7) is (nonlinearly) unstable; (iii) For A < 0, = 1, the function V takes positive values in any circle B and V_ is positive definite there. Then, the zero solution of (7) is (nonlinearly) unstable. h We now specialize the result obtained in Theorem 1 to the SIRS epidemic model case:
B. Buonomo, S. Rionero / Applied Mathematics and Computation 217 (2010) 4010–4016
4013
Corollary 1. If
ðl þ cÞ @f @I
r lþc þ 1
S¼S ;I¼I
@f @S
@f @S
S¼S ;I¼I
S¼S ;I¼I
@f @I
S¼S ;I¼I
þ r þ l > 0; ð12Þ
ðr þ lÞ ðl þ cÞ < 0;
then the generic equilibrium E = (S*, I*, R*) is nonlinearly locally asymptotically stable. If one of the two inequalities (12) holds reversed, then E is nonlinearly unstable.
l + c)[a Proof. The proof follows from Theorem 1 by observing that, from (8), A=p( ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi (r + l)] rb, and T = a r c 2l, and that the boundedness of the solutions of (1) guarantees that U ¼ o w2 þ v 2 (see system (6)). h Set now: R0 ¼
1 @f : ðr þ lÞ @I S¼1;I¼0
ð13Þ
Because of the hypothesis f(S, 0) = 0, it follows that (@ f/@S)S=1,I=0 = 0. Hence, in terms of the equilibrium E0, the inequalities (12) can be written:
ðr þ lÞð1 R0 Þ > 0;
ðr þ lÞðR0 1Þ ðl þ cÞ < 0:
ð14Þ
From inequalities (14) and Theorem 1, we get the following result concerning the stability/instability of the trivial equilibrium E0: if R0 < 1, then E0 is nonlinearly stable. If R0 > 1, then E0 is unstable. 4. Incidence rates In this section we will apply our result by considering several different incidence rates. We will consider the class:
f ðS; IÞ ¼ SgðIÞ:
ð15Þ
(i) Principle of mass action. In the classical SIR model with vital dynamics (see [4]), it is postulated that the spread of the disease among the population occurs according to the principle of mass action. This means that the incidence rate is bilinear respect to S and I. In other words,
gðIÞ ¼ kI;
ð16Þ 1
where k is a positive constant. From (13), (15) and (16) we have R0 = k(r + l) . Then, if R0 > 1, system (1), where the incidence rate f is given by (15) and (16), admits the unique endemic equilibrium E (S*, I*, R*), where:
S ¼ R1 0 ;
R ¼ rl1 I ;
I ¼ lg1 ðR0 1Þ;
and g = k[1 rc(l + r)1(l + c)1] > 0. It is easy to check that the inequalities (12) are satisfied, provided that R0 > 1. As a consequence, Theorem 1 ensures that the equilibrium E, when it exists, is nonlinearly stable. (ii) Incidence rates with ‘psychological’ effect. When the function g(I) in (15) is non monotone, in the sense that it is increasing when I is small and decreasing when I is large, then it may represent the so called ‘psychological’ effect. Large number of infectious may push the population to reduce the number of contacts per unit time, so that the infection force first increases and then decreases as I increases (inhibitory effect). An example is given by the following functional, which has been proposed in [20]:
gðIÞ ¼
kI 1 þ aI2
;
ð17Þ
where k and a are positive constants. In (17), kI measures the infection force of the disease and 1/(1 + a I2) describes the inhibitory effect from the behavioral change of the susceptible individuals when the number of infectious individuals is very large. The functional (17) is non monotone. From (13), (15) and (17) it follows: R0 = k(r + l)1. It is shown in [20] that if R0 > 1, system (1), where the incidence rate f is given by (15) and (17), admits the unique endemic equilibrium E (S*, I*, R*). In [20] it is also shown that such equilibrium is a stable node. In other word, the inequalities (12) are satisfied. Hence, Theorem 1 ensures that the equilibrium E is nonlinearly stable. (iii) A quadratic incidence rate. Consider:
gðIÞ ¼ kIð1 þ aIÞ;
ð18Þ
4014
B. Buonomo, S. Rionero / Applied Mathematics and Computation 217 (2010) 4010–4016
where k and a are positive constants. As stressed in [19], the use of this functional means the inclusion of an increased rate of infection due to two exposures over a short time period. The single contacts lead to infection at the rate kIS, whereas the new infectious individuals arise from double exposures at a rate akI2S. From (13), (15) and (18) we have R0 = k(r + l)1. The dynamics of this model is very rich. Following the same discussion given in [11], it can be seen that for R0 < 1, system (1), where the incidence rate f is given by (15)–(18), may admit none, one, or two endemic equilibria. Backward bifurcation occurs. For R0 > 1, there exists a unique endemic equilibrium E. Stability switches are also possible, so that the stability/instability of the endemic equilibria depends on the fulfillment of additional conditions on the parameters. For example, in the case R0 > 1, it may be shown that the determinant of the Jacobian matrix is always positive, whereas the sign of the trace depends on the parameters. The following result has been stated. Set:
r k 4ð1 þ qÞ2 q¼ ; B ¼ ; B0 ¼ 2 ; k ðl þ cÞ l q ðq þ 2Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð4q þ 4 þ 2 4q2 þ 8q þ 4 þ B2 q2 Þðq þ 1Þ ; p0 ¼ 2B2 q pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi B þ 2q þ 4 þ 2 ðB þ q þ 2Þðq þ 2Þ : p2 ¼ B2
p¼
aðl þ cÞ
;
Proposition 1 [11]. If B 6 B0 and p < p0, or B > B0 and p < p2, then there exist an endemic equilibrium E, which is locally stable. As a consequence, Theorem 1 ensures that the hypotheses of Proposition 1 guarantee also the nonlinear stability of E. 5. Global nonlinear asymptotic stability As stressed in Section 1, in [10] it has been shown that the SIRS model (1) admits a unique endemic equilibrium, which is globally asymptotically stable in R2þ , provided that R0 > 1 and the incidence rate satisfies some sufficient conditions. Precisely, the following inequalities must be satisfied, [10]:
f ðS; I Þ < f ðS ; I Þ for all 0 < S < S ; f ðS; I Þ > f ðS ; I Þ for all S > S ; and,
I=I 6 f ðS; IÞ=f ðS; I Þ 6 1 for 0 < I 6 I ; 1 6 f ðS; IÞ=f ðS; I Þ 6 I=I
ð19Þ
for I P I :
It can be seen that inequalities (19) rule out several incidence rates, as the ones considered, e.g. in [19,20], which are convex (i.e. concave up) respect to the variable I. Here, the use of the Lyapunov direct method, through the peculiar Lyapunov function (9), allows to show that the conditions for the linear stability are necessary and sufficient for the local nonlinear stability of hyperbolic equilibria, according to Theorem 1. This result implies that no restrictions are required about the monotonicity and concavity of the incidence rate with respect to the infectious class. Furthermore, if one is able to ‘control’ the functional W in (10), then sufficient conditions for the global stability may be obtained. In this section, we find such conditions in the case of quadratic incidence rate, (15)– (18), which is convex (concave up) with respect to the variable I. We may start from (4). Observe that (5) is now:
u ¼ kI ð1 þ aI Þu þ kS ð1 þ 2aI Þw þ ðk þ 2akI Þuw þ þakS w2 þ akuw2 ; where u = v w and S* = 1 I* R*. Hence, we obtain system (6), where now:
Uðw; v Þ ¼ ðk þ 2akI Þv w þ ðakS k 2akI Þw2 akv w2 akw3 ;
ð20Þ
and,
a ¼ kS ð1 þ 2aI Þ kI ð1 þ aI Þ;
b ¼ kI ð1 þ aI Þ:
We can now refer to system (7), by setting v (v1, v2) = (w, v)T; J = [aij], whose entries are given by
a11 ¼ a ðr þ lÞ;
a12 ¼ b;
a21 ¼ r;
a22 ¼ ðl þ cÞ;
where a and b are given in (21), and G = (g1, g2) (U, 0)T, where U is given by (20). Reasoning as in the proof of Theorem 1, we get equality (10), where now:
W ¼ a1 T v 1 Uðv 1 ; v 2 Þ;
ð21Þ
B. Buonomo, S. Rionero / Applied Mathematics and Computation 217 (2010) 4010–4016
4015
where T < 0, a1 is given by (11) and it is assumed that a11 – 0. Recall now that the solutions of (1) are bounded, so that we can consider uniformly bounded perturbations. That is, we assume that a positive constant m exists, m 6 2, such that
jv 1 j; jv 2 j < m: In this way, taking into account of (20), we get the estimate:
jWj 6 ma1 jTj
1 2 jk þ 2akI j v 21 þ v 22 þ jakS k 2akI jv 21 þ 2akmv 1 ; 2
that is,
jWj 6 d1
v 21 þ v 22
;
where,
d1 ¼ ma1 jTj½jk þ 2akI j þ jakS k 2akI j þ 2akm: Observe now that, from (9), it follows:
k1
v 21 þ v 22
< V < k2
v 21 þ v 22
;
ð22Þ
where
k1 ¼
AjTj ; 2
k2 ¼ jTj max A; 2 a211 þ a221 ; 2 a212 þ a222 :
ð23Þ
Hence, taking into account of (10), (22) and (23), it follows:
_ Vj ð7Þ 6 ðd3 þ d4 ÞV; where
d3 ¼
AT 2 ; k2
d4 ¼
d1 : k1
ð24Þ
Hence,
V 6 Vð0Þeðd3 d4 Þt : We can summarize the result of this section in the following theorem: Theorem 2. When the SIRS epidemic model (1) with the incidence force given by (18) admits a locally stable equilibrium, then such equilibrium is nonlinearly globally stable, provided that d3 > d4, where d3 and d4 are given by (24).
6. Conclusions In this paper we provide a direct application to epidemic models of a recent approach to nonlinear stability, introduced by Rionero in the context of L2-stability analysis for binary reaction–diffusion systems and here adapted to ODE systems. We find exact coincidence between the linear and nonlinear stability conditions for the endemic equilibrium and enlarge (with respect to what established in [10]) the class of incidence rates ensuring the nonlinear stability of endemic equilibria. Further investigations will include diffusion processes and higher dimensional epidemic systems, as SEIR-like models. Acknowledgements The present work has been performed under the auspices of the italian National Group for Mathematical Physics (GNFMINDAM). References [1] [2] [3] [4] [5]
R.M. Anderson, R.M. May, Infectious Diseases of Humans: Dynamics and Control, Oxford University Press, Oxford, 1991. S.N. Busenberg, K. Cooke, Vertically Transmitted Diseases, Biomathematics, vol. 23, Springer-Verlag, Berlin, 1993. V. Capasso, Global solution for a diffusive nonlinear deterministic epidemic model, SIAM J. Appl. Anal. 35 (1978) 274–284. V. Capasso, Mathematical Structures of Epidemic Systems, Lecture Notes in Biomath, vol. 97, Springer-Verlag, Berlin, 1993. V. Capasso, E. Grosso, G. Serio, I modelli matematici nella indagine epidemiologica Applicazione all’epidemia di colera verificatasi in Bari nel 1973, Annali sclavo 19 (1977) 193–208. [6] V. Capasso, G. Serio, A generalization of the Kermack-Mc Kendrick deterministic epidemic model, Math. Biosci. 42 (1978) 41–61. [7] A. d’Onofrio, Vaccination policies and nonlinear force of infection: generalization of an observation by Alexander and Moghadas (2004), Appl. Math. Comput. 168 (2005) 613–622. [8] H.W. Hethcote, The mathematics of infectious diseases, SIAM Rev. 42 (2000) 599–653.
4016
B. Buonomo, S. Rionero / Applied Mathematics and Computation 217 (2010) 4010–4016
[9] A. Korobeinikov, Lyapunov functions and global stability for SIR and SIRS epidemiological models with non-linear transmission, Bull. Math. Biol. 30 (2006) 615–626. [10] A. Korobeinikov, Global properties of infectious disease models with nonlinear incidence, Bull. Math. Biol. 69 (2007) 1871–1886. [11] Y. Jin, W. Wang, S. Xiao, An SIRS model with a nonlinear incidence rate, Chaos, Solitons Fractals 34 (2007) 1482–1497. [12] G.H. Li, Z. Jin, Global stability of an SEI epidemic model with general contact rate, Chaos, Solitons Fractals 23 (2005) 997–1004. [13] M.Y. Li, J.S. Muldowney, P. van den Driessche, Global stability of SEIRS models in epidemiology, Can. Appl. Math. Quart. 7 (1999) 409–425. [14] S. Rionero, L2 stability of solutions to a nonlinear binary reaction-diffusion system of P.D.E.s, Rend. Math. Acc. Lincei 16 (2005) 227–238. [15] S. Rionero, A nonlinear L2 stability analysis for two species dynamics with dispersal, Math. Biosci. Eng. 3 (2006) 189–204. [16] S. Rionero, A rigorous reduction f the L2-stability of the solutions to a nonlinear binary reaction-diffusion system of P.D.E.s to the stability of the solutions to a linear binary system of ODE’s, J. Math. Anal. Appl. 319 (2006) 377–397. [17] S. Rionero, Nonlinear L2-stability analysis for two-species population dynamics in spatial ecology under Neumann boundary data, Rend. Circ. Math. Palermo. Ser. II 78 (Suppl.) (2006) 273–283. [18] S. Rionero, On the nonlinear stability of the critical points of an epidemic SEIR model via a novel Lyapunov function, Rend. Acc. Sci. Fis. Nat. Napoli 75 (4) (2008) 115–129. [19] P. van den Driessche, J. Watmough, A simple SIS epidemic model with backward bifurcation, J. Math. Biol. 40 (2000) 525–540. [20] D. Xiao, S. Ruan, Global analysis of an epidemic model with nonmonotone incidence rate, Math. Biosci. 208 (2007) 129–419.