Available online at www.sciencedirect.com
Physics Procedia 00 (2011) Physics Procedia 22 (2011) 20 –000–000 39 www.elsevier.com/locate/procedia
ICPST 2011
Analytic Comparison of Some Epidemic Models with Vaccination M. De la Sen a *, S. Alonso-Quesada a , A. Ibeas b a
b
Department of Electricity and Electronics, UPV/EHU, Leioa, 48080-Bilbao, Spain Department of Telecommunications and Systems Engineering, Universitat Autònoma de Barcelona, 08193-Bacelona, Spain
Abstract In this paper, we discuss the elementary properties of some simple SI, SR, SIR and SEIR epidemic models whose parameterizing functions (such as per-capita death rate, disease transmission, removal rate etc.) might be eventually time-varying but nonnecessarily time-integrable. Vaccination rules based of feedback, measuring the numbers of some of the partial populations defining the disease progress, are also discussed.
© Published Elsevier © 2011 2010 Published by by Elsevier B.V. B.V. Selection and/or peer-review under responsibility of Garry Lee. PACS:028987;028789 Keywords:Epidemic models; SI (susceptible/infectious); SR (susceptible/immune); SIR (susceptible/infectious/immune) ; SEIR (susceptible/ infected/infectious/immune); vaccination control.
1. Introduction Important control problems nowadays related to Life Sciences are the control of ecological models like, for instance, those of population evolution (Beverton-Holt model, Hassell model, Ricker model etc.) via the online adjustment of the species environment carrying capacity, that of the population growth or that of the regulated harvesting quota as well as the disease propagation via vaccination control. In a set of papers, several variants and generalizations of the Beverton-Holt model (standard time–invariant, timevarying parameterized, generalized model or modified generalized model) have been investigated at the levels of stability, cycle- oscillatory behavior, permanence and control through the manipulation of the carrying capacity (see, for instance, [1-5]). The design of related control actions has been proved to be important in those papers at the levels, for instance, of aquaculture exploitation or plague fighting. On the other hand, the literature about epidemic mathematical models is exhaustive in many books and papers . A non-exhaustive list of references is given in this manuscript, cf. [6-14] (see also the references listed therein). The sets of models include the most basic ones, [6-7]: SI- models where not removed- by – immunity population is assumed. In other words, only susceptible and infected populations are assumed. ____________ *Corresponding author. Tel.:+0-0034946012548; fax: +0-0034946012548 E-mail address:
[email protected]
1875-3892 © 2011 Published by Elsevier B.V. Selection and/or peer-review under responsibility of Garry Lee. doi:10.1016/j.phpro.2011.11.006
M. De la Sen/ Sen etPhysics al. / Physics Procedia 22 (2011) 20 – 39 Procedia 00 (2011) 000–000
-
SIR models, which include susceptible plus infected plus removed- by –immunity populations. SEIR- models where the infected population is split into two ones (namely, the “ infected” which incubate the disease but do not still have any disease symptoms and the “ infectious” or “ infective” which do have the external disease symptoms). Those models have also two major variants, namely, the so-called “pseudo-mass action models”, where the total population is not taken into account as a relevant disease contagious factor and the so-called “true-mass action models”, where the total population is more realistically considered as an inverse factor of the disease transmission rates. There are many variants of the above models, for instance, including vaccination of different kinds: constant [8], impulsive [12], discrete – time etc., incorporating point or distributed delays [12-13], oscillatory behaviours [14-18] etc. On the other hand, some ´ad- hoc´ variants of such models are known to become considerably simpler for the disease transmission among plants [67]. In [19], a control point of view of a vaccination strategy in continuous- time has been proposed for the true mass action (namely, the whole population numbers influence the rate of disease transmission) socalled SEIR (i.e. susceptible/infected/infectious and immune populations) epidemic model under constant whole population assumption. This model generalizes simpler SIR epidemic models where infected (i.e. those still without symptoms) and infectious (i.e. those already with disease symptoms) are not mutually distinguished. The vaccination strategy involves an auxiliary control being proportional to either the susceptible or to the whole population so that the unsuitable dynamics is removed and replaced for an asymptotically stabilizing term of the susceptible dynamics. The disease propagation is studied in a number of papers (see, for instance, [20-27] and references there in). In this paper, we first discuss three elementary epidemic models of respective types as follows: SI (susceptible/infectious), SR (susceptible/ immune) and SIR (susceptible /infected/ immune) models whose parametrizing functions (as for instance, per-capita death rate, disease transmission etc.) might be eventually time-varying but either timeintegrable or not. Furthermore, a more general SEIR (susceptible/infected/infectious/immune) model is also analyzed in this paper. Finally, the effect of some vaccination policies on the populations is discussed for this model. 2. A time-varying SI epidemic model Bernouilli proposed in 1760 a simple epidemic model where the infection is removed instantaneously so that all the population passes from susceptible to removed by immunity (the simplest SR model), [20]. The model was assumed in particular for instantaneous infective effect via inoculation of the smallpox. Since then a lot of investigation has been devoted to epidemic models including the incorporation of infected and infectious populations (SIR and SEIR epidemic models), the presence of delays in the disease transmission etc. The following simple one-parameter time-varying model is a generalization to the time- varying case of the simpler time- invariant SI (susceptible/infectious) epidemic model: N = S( t ) + I ( t ) (1) I ( t ) = β (t ) S (t ) I (t ) (2) = β (t ) S (t ) (N − S (t ) ) = β (t ) I (t ) (N − I (t ) )
(3) with initial conditions N ≥ S (0) = S 0 ≥ 0 , N ≥ I (0) = I 0 = N − S 0 ≥ 0 such that the total population N (t ) = S (t ) + I (t ) = S 0 + I 0 = N (0) = N 0 is constant for all time and the disease transmission function β : R 0+ → R 0+ with R 0+ : = R + ∪ { 0 } . Thus, one gets for N > 0: ⎛ 1 ⎞ dS (t ) ⎞ dI ( t ) ⎛ 1 d S (t ) 1 dI ( t ) 1 ⎟⎟ ⎟⎟ =− = − ⎜⎜ + = β (t ) dt (4) = ⎜⎜ + S (t ) ( N − S (t )) I ( t ) ( N − I ( t )) ⎝ I ( t ) N − I ( t ) ⎠ N ⎝ S (t ) N − S (t ) ⎠ N so that I (t ) I ( 0) S ( 0) S (t ) (5) ln − ln = ln − ln = N β ( t ) ; ∀ t ∈ R 0+ N − I (t ) N − I (0) N − S ( 0) N − S (t )
21
22
Sen et al. / Physics 22 (2011) 20 – 39 DeM. la De Sen/laPhysics Procedia 00Procedia (2011) 000–000
where β : R 0 + → cl R 0 + (i.e. with the image space being the closure of the nonnegative real numbers so that the + ∞ - point is added) is defined by β ( t ) :=
β ∈ L 1 ( R 0 + , R 0 + ) and β (t = ∞ ) = β
∞=
t
∫ β ( τ ) dτ ≤ β (t = ∞) = β 0
∞
with 0 ≤ β
∞<
∞ if
+ ∞ , otherwise. The solution of (1)-(2) is obtained from (5)
as follows: I (t ) =
N I ( 0)
I ( 0 ) + ( N − I ( 0) ) e − N β ( t )
S (t ) =
N (N − I (0) ) e − N β ( t )
=
N − I (0) I (t ) e − N β ( t ) ; ∀ t ∈ R 0 + I (0)
(6) I (0) + (N − I (0) ) e which are nonnegative for all time , so that (1)-(3) is a positive dynamic system ( see [15-17]), and have finite nonnegative limits as t → ∞ which is a global attractor of the trajectory- solution and it is also a globally asymptotically stable endemic ( in the sense that the disease propagates) equilibrium point: N I ( 0) (7.a) I (∞ ) = ∈ [0, N ] I ( 0) + ( N − I ( 0) ) e − N β ∞
S (∞ ) =
− N β (t )
N (N − I (0) ) e − N β I (0) + (N − I (0) )
∞
e− N β ∞
=
N − I (0) I (∞) e − N β I (0)
∞
= N − I (∞) ∈ [0, N
]
(7.b)
which becomes in particular if β ∞ = + ∞ the following strongly endemic ( in the sense that the whole population becomes infectious) equilibrium point: I (∞) = N ; S (∞) = 0 if I (0) ≠ 0 , and I (∞) = I (t ) = I (0) = 0 ; S (∞) = S (t ) = S (0) = N if I (0) = 0 ; ∀ t ∈ R 0 + (8) Thus, the solution of (1)-(2) is nonnegative for all time if the initial conditions are nonnegative and converges asymptotically to a stable equilibrium point for any given nonnegative initial conditions satisfying a constant population constraint N = S (0) + I (0) . Simple calculations yield: dI ( ∞ ) N 2 e −N β ∞ d S( ∞ ) d I( ∞ ) N 2 e −N β ∞ ; (9) = =− =− −N β ∞ ⎞ 2 dI ( 0 ) ⎛ I ( 0 ) + (N − I ( 0 ) ) e − N β ∞ ⎞ 2 dI ( 0 ) dI ( 0 ) ⎛ ⎜ ⎟ ⎜ I ( 0 ) + (N − I ( 0 ) ) e ⎟ ⎝ ⎠ ⎝ ⎠ dS ( ∞) what leads to = −1 and the following cases can occur: dI (∞) dI ( ∞) d S ( ∞) >0 , < 0 . Thus, if I (0) increases (decreases) then the susceptible 1) If 0 ≤ β ∞ < ∞ then dI (0) dI (0) limit decreases (increases) and the infected limit increases ( decreases).
2) If β ∞ = + ∞ or if N = 0 (leading to the trivial solution of (1)-(2)) then dI ( ∞) d S ( ∞) =0 , = 0 . Then, the equilibrium point coincides with the initial conditions. dI (0) dI (0) d S( ∞ ) dI ( ∞ ) = 1 >0, = −1 < 0 . Thus, if I (0) increases (decreases) then the limit 3) If β ∞ = 0 then dI ( 0 ) dI ( 0 ) infected increases (decreases) and the limit susceptible decreases (increases). The solution of (1)-(2) may be alternatively written as follows with given upper-bounding functions:
23
M. De la Sen/ Sen etPhysics al. / Physics Procedia 22 (2011) 20 – 39 Procedia 00 (2011) 000–000 ⎛
⎞
⎜ ⎝
⎟ τ ∈[ 0 , t ) ⎠
β ( t )⎜ N − max I (τ )⎟
e
⎛
⎞
⎜ ⎝
τ ∈[ 0 , t ) ⎠
β ( t ) ⎜ N − min I (τ )⎟
N −e
+e
t ∫0
β (τ
⎟
I ( 0 ) ≤ I (t ) = e
t ∫0
β (τ
)( N − I (τ ) ) dτ
⎛
⎞
⎜ ⎝
τ ∈[ 0 , t ) ⎠
β ( t )⎜ N − min I (τ )⎟
I ( 0) ≤ e
⎟
I (0) ; ∀ t ∈ R 0 +
(10)
t ⎛ ∫ β (τ ) ( N − I (τ ) ) dτ ⎞⎟ I (0) ≤ S (t ) = ⎜1 − e 0 N ⎜ ⎟ ⎝ ⎠
) ( N − I ( τ ) )dτ
⎛
⎞
⎜ ⎝
τ ∈[ 0 , t ) ⎠
β ( t )⎜ N − max I (τ )⎟
S( 0 ) ≤ N − e
⎟
I (0)
; ∀ t ∈ R 0 + (11) Also, one gets from (6):
dI ( t ) N 2 I ( 0 ) (N − I ( 0 )) e − N β ( t )β (t ) dI (t ) N 2e − N β ( t ) ; ; ∀ t ∈ R 0 + (12) = = 2 2 dI (0) ⎛ dt − N β (t ) ⎞ − N β (t ) ⎞ ⎛ ⎟ ⎜ I (0) + (N − I (0) ) e ⎜ I (0) + (N − I ( 0 ) ) e ⎟ ⎠ ⎝ ⎝ ⎠ Thus, the infected population at any time t increases (decreases) when I (0) is increased (decreased) .
The susceptible population behaves in the contrary sense. That is I ´ (t ) > I (t ) and S ´ (t ) < S (t ) at any time if I ´ (0 ) > I (0) and conversely. Also, the infected (susceptible) population is a monotone strictly increasing (decreasing) function for all time for any initial condition in (0, N ] . 3. A time-varying SR epidemic model
A simple SR (susceptible/ immune –also called removed by immunity-) time-varying epidemic model extending its time- invariant counterpart is (see [20]):
N ( t ) = −μ (t ) N ( t )
(13)
S ( t ) = − ( μ (t ) + β (t )) S (t ) R ( t ) = N (t ) − S (t ) = − μ (t ) R(t ) + β (t ) S (t )
(14)
(15) under initial conditions N (0) = N 0 = S (0) + R(0) ≥ 0 , N 0 ≥ S (0) = S 0 ; R(0) = R 0 = N 0 − S 0 ≥ 0 from which two differential equations are independent. The transmission function is β : R 0 + → R 0 + and μ : R 0+ → R 0+ is the per capita death ratio at time t . The unique solution of (13)-(15) is: − ∫ 0 ( μ (τ ) + β (τ )) d τ t
S (t ) = e
R(t ) = N (t ) − S (t ) = e
which
take
− ∫ 0 μ (τ ) d τ t
N (t ) = e
S (0)
t − ∫0
nonnegative
μ (τ ) dτ ⎛
⎜ N (0) − e ⎜ ⎝
real
values
t − ∫0
β (τ ) dτ
for
all
β ( t ):= ∫ 0 β ( τ )d τ ∈ cl R 0+ with 0 < μ (t ) ≤ μ t
μ
∞ :=
lim μ (t ) < ∞
t →+ ∞
if
N (0) ; ∀ t ∈ R 0 +
(16.a)
⎞ S (0) ⎟ ; ∀ t ∈ R 0 + ⎟ ⎠
(16.b)
time.
∞≤ ∞
Define
cl R
0+
∋ μ ( t ):= ∫ 0 μ ( τ )d τ t
;
and 0 < β (t ) ≤ β ∞ ≤ ∞ ; ∀ t ∈ R 0 + with
and only if μ ∈ L 1 ( R 0+ , R 0+
)
(otherwise, μ
∞= + ∞
) , and
24
β
Sen et al. / Physics 22 (2011) 20 – 39 DeM. la De Sen/laPhysics Procedia 00Procedia (2011) 000–000 ∞ :=
lim β (t ) < ∞ if and only if β ∈ L 1 ( R 0 + , R 0 + ) ( otherwise, β ∞ = + ∞ ). From (16), the
t →+ ∞
following properties hold: 1) If μ ∞ = + ∞ then N (∞ ) = S (∞) = R ( ∞ ) = 0 irrespective of the initial conditions so that the whole and the two partial populations asymptotically extinguish. 2) If β ∞ = + ∞ then S (∞) = 0 and R (∞ ) = N (∞) = e − μ ∞ N (0) so that the susceptible population asymptotically extinguishes and the whole one is asymptotically immune identical to a finite limit which depends on the initial total population and the value μ ∞ . If μ ∞ = + ∞ then such a limit is zero so that the whole population again asymptotically extinguishes as in the above case 1. 3) It follows from (13)-(15) that dN ( t ) / dN (t ) = − μ (t ) ≤ 0 , dS ( t ) / dS (t ) = − ( μ (t ) + β (t )) ≤ 0 dR ( t ) / dR(t ) = − μ (t ) ≤ 0 , dR ( t ) / dS (t ) = β (t ) ≥ 0 Thus, the whole population and the susceptible and immune populations are monotone decreasing functions so that they have finite limits. Also, the immune population increases (decreases) as the susceptible one increases (decreases). An equivalent result to (16.b) for the immune population is calculated directly from (15) and (16.a) as follows ∀ t ∈ R 0 + : R(t ) = e
− ∫ μ (τ ) d τ t 0
τ
R(0)
t − ∫ μ (τ ´) d τ ´ + ∫0 e 0 β
(t − τ )S (t − τ ) d τ
− ∫ 0 μ (τ ) d τ t
=e
t −τ ⎛ t − ∫τ μ (τ ´) d τ ´ ( μ (τ ´) + β (τ ´) ) d τ ´ ⎞⎟ −∫ + ⎜ ∫0 e 0 β ( t − τ )e 0 d τ S (0 ) ⎜ ⎟ ⎝ ⎠
R(0)
− ∫ 0 μ (τ ) d τ t
=e
R (0)
t −τ ⎛ t − ∫τ μ (τ ´) dτ ´ −∫ ( μ (τ ´) + β (τ ´) ) d τ ´ ⎞⎟ + ⎜ ∫0 e 0 β ( t − τ )e 0 d τ (N (0 ) − R (0) ) ⎜ ⎟ ⎝ ⎠
t = ⎛⎜ e − μ ( t ) − ∫ 0 e − (μ ( τ ) + μ ( t −τ )+β (t −τ ))β (t − τ ) dτ ⎞⎟ R ( 0 ) + ⎛⎜ ∫ t e − (μ ( τ ) + μ ( t −τ )+β (t −τ ) )β (t − τ ) dτ ⎞⎟ N (0 ) ⎝ ⎠ ⎝ 0 ⎠ = e − μ ( t ) R (0) + ⎛⎜ e − (μ ( t )+ β ( t )) ∫0t e μ (τ )+ β (τ )β (t − τ ) dτ ⎞⎟ S ( 0 ) (17) ⎝ ⎠
∞ β (τ ) β (t − τ R (∞ ) := e − μ ∞ R ( 0 ) + ⎛⎜ ∫ 0 e
) dτ
⎞⎟ S ( ∞ ) ⎠ ∞ β (τ ) −μ ∞ ⎛ β (t − τ ) dτ ⎞⎟ S ( ∞ ) ( N (0) − S (0)) + ⎜ ∫ 0 e =e ⎝ ⎠ ∞ = N ( ∞ ) + ⎛⎜ ∫ 0 e β ( τ ) β (t − τ ) dτ ⎞⎟ S ( ∞ ) − e − μ ∞ S ( 0 ) ⎝ ⎠ ⎝
(18)
4. A time-varying SIR epidemic model
A simple SIR (susceptible/infectious/immune) time-varying epidemic model extending its timeinvariant counterpart (Kermack-McKendrick model - 1927), [20] is: S ( t ) = − β (t ) S (t ) I (t ) (19) I ( t ) = β (t ) S (t ) I (t ) − γ (t ) I (t ) (20)
25
M. De la Sen/ Sen etPhysics al. / Physics Procedia 22 (2011) 20 – 39 Procedia 00 (2011) 000–000
R (t ) = γ (t ) I (t ) (21) γ : R 0+ → R 0+ is the removal rate with initial conditions where S (0) = S 0 ≥ 0 , I (0) = I 0 ≥ 0 , R(0) = R 0 ≥ 0 and where the new parameterizing function, right handsides of (19)-(21), yields a constant total population N = N ( 0 ) = N 0 = S 0 + I 0 + R 0 . The use of (21) in β (t ) R (t ) (19) yields S ( t ) = − S ( t ) which together with (19) yields the two equivalent expressions for γ (t ) the susceptible: t
S (t ) = e
−∫ 0 β (τ ) I (τ ) d τ
S (0) = e
(
)
t −∫ 0 β (τ ) R (τ ) /γ (τ ) d τ
S (0) ; ∀ t ∈ R 0 + (22) Eq. 20 might be rewritten equivalently as I ( t ) = (β (t ) S (t ) − γ (t )) I (t ) what leads to the following solution: ∫ (β (τ ) S (τ ) −γ (τ ) ) d τ I (t ) = e 0 I (0) ≥ 0 ; ∀ t ∈ R 0 + which is well-posed if t
∫ 0 (β (τ ) S (τ ) −γ (τ ) ) d τ
(23)
t
e or, equivalently, if ⎛ ∫ t (β ( τ ) S ( τ )−γ ( τ ) ) d τ ⎜ 0 ⎜ ⎝
e
I (0) ≤ N = S (0) + I (0) + R(0) ; ∀ t ∈ R 0 +
(24)
⎞ − 1⎟ I ( 0 ) ≤ S ( 0 ) + R( 0 ) ; ∀ t ∈ R 0 + ⎟ ⎠
(25)
∫ (β (τ ) S (τ ) −γ (τ ) ) d τ ≤ 1; which holds irrespective of any nonnegative values of the initial conditions if e 0 ∀ t ∈ R 0 + . That is guaranteed if β ( t ) ≤ γ (t ) / S (t ) ; ∀ t ∈ R 0 + . Focusing on the susceptible population given by (22), one gets that if β ( t ) ≤ γ (t ) / S (t ) ; ∀ t ∈ R 0 + , guaranteeing non- negativity for all time of the infectious population, then ⎛ −∫ t β (τ ) I (τ ) d τ ⎞ ⎜e 0 − 1⎟ S (0) ≤ 0 ≤ I (0) + R(0) ; ∀ t ∈ R 0 + (26) ⎜ ⎟ ⎝ ⎠ so that the susceptible population is also nonnegative for all time under the same sufficiency – type condition as the infectious one is nonnegative, that is, β ( t ) ≤ γ (t ) / S (t ) for all time. On the other hand, it also follows directly from integration of (22)-(23) through time that: t
R(t ) = R(0) + ∫ 0 γ (τ )I (τ ) dτ t
τ ⎛ t ∫ (β (τ ´) S (τ ´) −γ (τ ´) ) d τ ´ dτ = R (0) + ⎜ ∫ 0 γ (τ )e 0 ⎜ ⎝
⎞ ⎟ I ( 0) ; ∀ t ∈ R 0 + ⎟ ⎠
(27)
which is nonnegative for all time if the infectious population is also nonnegative for all time which is guaranteed if the disease transmission function is sufficiently small to satisfy the upper- bounding condition β ( t ) ≤ γ (t ) / S (t ) ; ∀ t ∈ R 0 + . Then, it follows as a global result that if β ( t ) ≤ γ (t ) / S (t ) ; ∀ t ∈ R 0 + then 0 ≤ S (t ), I (t ) , R(t ) ≤ N for all time. Since this condition is always guaranteed for t=0, it follows by complete induction that if the stronger condition β ( t ) ≤ γ (t ) / N ; ∀ t ∈ R 0 + holds then 0 ≤ S (t ), I (t ) , R(t ) ≤ N for all time for any set of well-posed initial conditions so that the given SIRmodel is positive. On the other hand, it follows by observing (27) that the condition β ( t ) ≥ γ (t ) / S (t ) ;
26
Sen et al. / Physics 22 (2011) 20 – 39 DeM. la De Sen/laPhysics Procedia 00Procedia (2011) 000–000
∀ t ∈ R 0 + implies that R (t ) → + ∞ if γ (t ) ∈ R + for all time what makes the epidemic mathematical model to be ill-posed. However, those conditions can fail on time intervals of finite measures remaining still the model well-posed. A necessary condition for a well-posed model is R (t ) ≤ N = S (0) + I (0) + R(0) ; ∀ t ∈ R 0 + for any given initial conditions what in view of (27) translates into the following constraint: t ∫0 γ
∫ 0 (β (τ ´) S (τ ´) −γ (τ ´) ) d τ ´ τ
(τ )e
dτ − 1 ≤ S (0) / I (0) ; ∀ t ∈ R 0 +
(28)
t ⎛ ⎞ − ∫ γ (τ ) d τ The worst case of (28) is when S (0) = 0 ⇒ ⎜ S (t ) = 0 ∧ I (t ) = e 0 I (0) ⎟ ; ∀ t ∈ R 0 + from (22) – ⎜ ⎟ ⎝ ⎠ (23) so that the necessary condition, irrespective of the disease transmission function, for well-posed model for the worst-case of (28), and then valid for any set of initial conditions, is
t ∫0 γ
(τ ) e
τ
− ∫ 0 γ (τ ´ ) d τ ´
dτ ≤ 1 ; ∀ t ∈ R 0 +
(29)
what is guaranteed in particular with an upper-bounding function of exponential order K e − α t of the left-hand-side of (29) if there exist constants α ∈ R 0+ , K (≤ α )∈ R + such that
α ≤ min (1 / t ) ln ( K / γ (t ) ) + ∫ 0 γ (τ ) dτ t
(30)
t ∈ R 0+
The combination of (19) and (20) leads to: I ( t ) = − γ (t ) I (t ) − S ( t ) so that
(31)
0 ≤ I (0) + S (0) − (I (t ) + S (t ) ) = ∫ 0 γ (τ ) I (τ ) d τ t
t = ∫ 0 R (τ ) d τ = R (t ) − R(0) ; ∀ t ∈ R 0 +
(32)
where (21) has been used and provided that the model is well-posed so that the infectious population is nonnegative for all time what implies: I (0) + S (0) ≥ I (t ) + S (t ) ; N (0) = I (0) + S (0) + R(0) = N (t ) = I (t ) + S (t ) + R(t ) (as expected); ∀ t ∈ R 0 + so that the joint susceptible plus infected population is a monotone decreasing real function independent of initial conditions and the whole population is constant [the constant property of the whole population was already known from simple inspection of (19)-(21)]. Equivalently, (20) may be integrated via integration by parts as follows by also using (19): t τ t τ − ∫ γ (τ ) d τ − ∫ γ (τ ´) d τ ´ −∫ γ (τ ) d τ − ∫ γ (τ ´) d τ ´ t τ =t I (t ) = e 0 I ( 0 ) − ∫0 e 0 S ( t − τ ) dτ = e 0 I ( 0) − e 0 S (t − τ ) τ = 0 t
− ∫0 e t
=e
τ
− ∫ 0 γ ( τ´) d τ ´
− ∫ 0 γ (τ ) dτ
γ ( τ ) S ( t − τ ) dτ ⎛
( I (0) − S (0) ) − ⎜⎜ ∫0t e
τ
− ∫ 0 γ ( τ´) d τ ´
t −τ
−∫0
β ( τ´) I ( τ´) d τ ´
⎞
γ (τ ´ ) dτ´ ⎟ S ( 0) + S (t )
(33) ⎟ ⎝ ⎠ ; ∀ t ∈ R 0 + . Then, under any of the conditions β ( t ) ≤ γ (t ) / S (t ) ; ∀ t ∈ R 0 + or β ( t ) ≤ γ (t ) / N ; ∀ t ∈ R 0 + guaranteeing that the mathematical model is positive, it follows that : e
27
M. De la Sen/ Sen etPhysics al. / Physics Procedia 22 (2011) 20 – 39 Procedia 00 (2011) 000–000 t
e
− ∫ 0 γ (τ ) d τ
τ
t
+ ∫0 e
− ∫ 0 γ ( τ´) d τ ´
e
t −τ
e
−∫0
β ( τ´) I ( τ´) d τ ´
γ (τ ´ ) dτ´
t − ∫ 0 γ (τ ) d τ
≤
I( 0 ) ; ∀ t ∈ R 0+ S( 0 )
(34)
irrespective of any well-posed set of initial conditions. Now, one can easily obtain from (19) to (21): S ( t ) + I (t ) = − R ( t ) = − γ ( t ) ( N − S (t ) − R ( t ) ) ; ∀ t ∈ R 0 + (35) Then R ( t ) = − γ ( t ) R ( t ) + γ ( t ) ( N − S (t ) ) ; ∀ t ∈ R 0 + (36) which has the following solution by using (22): t τ − ∫ γ (τ ) d τ − ∫ γ ( τ´) d τ ´ t R( t ) = e 0 R( 0 ) + ∫0 e 0 γ ( t −τ )( N − S ( t − τ ) ) dτ t
=e
− ∫ 0 γ (τ ) d τ
t
=e
− ∫ 0 γ (τ ) dτ
t
τ
R( 0 ) + ∫0 e t
R( 0 ) + ∫ 0 e
− ∫ 0 γ ( τ´) d τ ´ ⎛
τ
t−τ
⎜ N −e −∫ 0 ⎜ ⎝
β ( τ ´) I ( τ ´ ) d τ ´
t−τ
− ∫ 0 γ ( τ´) d τ ´ ⎛
⎜ S ( 0 ) + I ( 0 ) + R( 0 ) − e − ∫ 0 ⎜ ⎝
⎞ S ( 0 ) ⎟γ ( t −τ ) dτ ⎟ ⎠
β ( τ ´) I ( τ ´ ) d τ ´
⎞ S ( 0 ) ⎟γ ( t −τ ) dτ ⎟ ⎠
τ t−τ ⎡ ⎛ ⎛ t − ∫ τ γ ( τ´) d τ ´ ⎞ −∫ β ( τ ´) I ( τ ´ ) d τ ´ t − ∫ γ ( τ´) d τ ´ ⎜ 1−e 0 = ⎜ ∫0 e 0 γ (t − τ )d τ ⎟ I ( 0 ) + ⎢ ∫0 e 0 ⎜ ⎜ ⎟ ⎢ ⎝ ⎝ ⎠ ⎣
⎞ ⎟ γ (t − τ )dτ ⎟ ⎠
⎤ ⎥ S( 0 ) ⎥ ⎦
⎛ − ∫ t γ ( τ ) d τ t − ∫ τ γ ( τ´) d τ ´ ⎞ +⎜ e 0 + ∫0 e 0 γ (t − τ )d τ ⎟ R( 0 ) ≥ 0 , ∀ t ∈ R 0 + (37) ⎜ ⎟ ⎝ ⎠ If β ( t ) ≤ γ (t ) / S (t ) ; ∀ t ∈ R 0 + or β ( t ) ≤ γ (t ) / N ; ∀ t ∈ R 0 + for any set of well-posed initial conditions, this leads to the following integral constraint: τ t−τ ⎛ − ∫ 0 β ( τ ´) I ( τ ´ ) d τ ´ ⎞ t − ∫ 0 γ ( τ´) d τ ´ ⎜ ⎟ γ (t − τ )dτ ≥ 0 − e 1 e (38) ∫0 ⎜ ⎟ ⎝ ⎠ what follows directly by contradiction arguments by choosing initial conditions S (0) = N , I (0) = R (0) = 0 so that if (38) is false for some time t then R (t ) < 0 under the positivity constraint for the model (19)-(21) β ( t ) ≤ γ (t ) / S (t ) ; ∀ t ∈ R 0 + or β ( t ) ≤ γ (t ) / N ; ∀ t ∈ R 0 + what is impossible. Note that positivity, [28-30] is an important property, as it is that of positive realness of transfer matrices, of the solution in some kinds of problems concerning dynamic systems. Some simple concerns with the spreading or not of the disease are now discussed. 1) From (20) and (23), I (0) = 0 ⇒ I (0) = 0 ∧ I (t ) = I(t ) = 0; ∀ t ∈ R 0+ and the disease does not spread through time. 2) Assume that I (0) > 0 and S (t ) = γ (t ) /β (t ) ; ∀ t ∈ R 0 + (a necessary condition being β (t ) ∈ R ). Then , one gets I (t ) = 0 , so that I (t ) = I (0) and S (t ) = −γ (t ) I (t ) ; ∀ t ∈ R and
(
)
+
0+
the disease still does not spread out through time so that t t S (t ) = γ (t ) /β (t ) = S (0) − ∫ 0 γ (τ ) I (τ ) d τ = γ (0) /β (0) − I (0) ∫ 0 γ (τ ) d τ ; ∀ t ∈ R 0 +
from (19)-(20). The above equation is equivalent to the subsequent rule for the disease transmission function:
28
Sen et al. / Physics 22 (2011) 20 – 39 DeM. la De Sen/laPhysics Procedia 00Procedia (2011) 000–000
β ( t)=
γ (t γ (0 ) β
−1
)
(0 )− I (0) ∫ 0t γ (τ ) d τ
; ∀ t ∈ R 0+
(39)
Since the denominator of (39) cannot be negative at any time, it is necessary that ∫ 0t γ (τ ) dτ
≤
γ (0 ) ; β (0 )I (0)
∀ t ∈ R 0 + . if β (t ) ≤ γ (t ) / S (t ) , so that it is also required as a result γ ∈ L 1 ( R 0+ , R 0+ ) (needing in
addition γ ( t ) → 0 as t → +∞ ) with sufficiently small bound for the integral on R 0+ depending on the initial infectious population and the initial values of the average initial disease removal rate and initial disease transmission function value. Note that dI (t ) / dR( t ) = 0 if S ( t ) = γ (t ) / β (t ) In this case, it follows from (21) that the immune population is given by:
R( t ) = R( 0 ) + I ( 0 ) ∫ 0 γ (τ ) d τ ; ∀ t ∈ R 0 + t
(40)
which introduced an additional constraint on the average removal rate of the disease which together with the former one deriving from (39) yields the stronger necessary constraint for non-propagation of the disease under initial nonzero infected population which remains constant for all time: 1 γ (0 ) t ; ∀ t ∈ R 0+ (41) ∫ 0 γ (τ )d τ ≤ I ( 0 ) β (0) One has also from (20)-(21) that β ( t ) S (t ) − 1≤ 0 ; ∀ t ∈ R 0 + so that the infected population is a monotone decreasing d I (t ) / d R (t ) = γ ( t) function (eventually being constant) with respect to the immune one. 3) The disease still does not spread if I (0) > 0 and S (t ) ≤ γ (t ) /β (t ) ; ∀ t ∈ R 0 + what implies I( t ) ≤ 0 ; ∀ t ∈ R so that I (t ) ≤ I (0) ; ∀ t ∈ R . Then, 0+
γ ( t ) / β (t ) ≥ S ( t ) ≥
t S( 0 ) − ∫ 0 γ
0+
(τ ) I (τ ) d τ ; ∀ t ∈ R 0 +
so that , one has instead of (39),
β ( t )≤
γ (t
)
S (0 ) − I (0) ∫ 0 γ (τ ) d τ t
; ∀ t ∈ R 0+
(42)
which, by taking also into account (40), requires the necessary condition below: t ∫ 0 γ (τ )dτ <
S( 0 ) I( 0 )
; ∀ t ∈ R 0+
(43)
This condition is weaker than (41) but still guarantees that the disease does not spread with the infectious population being a monotone decreasing function including eventually to be a constant function defined by the initial conditions. 4) The disease spreads through time if I (0) > 0 and S (t ) ≥ γ (t ) /β (t ) : ∀ t ∈ R 0 + and , furthermore, S (t ) > γ (t ) /β (t ) ; ∀ t ∈ I T ⊂ R 0+ with the interval I T being some non necessarily →R is connected subset of R of infinite measure. Thus, I( t ) > 0 ; ∀ t ∈ IT so that I : R 0+
0+
monotone increasing. Then,
S ( t ) ≤ S ( 0 ) − ∫ 0 γ (τ ) I (τ ) d τ ≥ R t
p
:= γ ( t ) / β (t ) ; ∀ t ∈ R 0 +
S ( t ) ≤ S ( 0 ) − ∫ 0 γ (τ ) I (τ ) d τ > γ ( t ) / β (t ) ; ∀ t ∈ IT t
what is guaranteed under (43) if
0+
29
M. De la Sen/ Sen etPhysics al. / Physics Procedia 22 (2011) 20 – 39 Procedia 00 (2011) 000–000
β ( t )≥
γ (t
)
t S (0 ) − I ( 0) ∫ 0
γ (τ ) d τ
; ∀ t ∈ R 0+ ;
β ( t)>
γ (t
)
t S (0 ) − I (0) ∫ 0
γ (τ ) d τ
; ∀ t ∈ IT
(44)
and R p := γ (t ) /β (t ) is said to be the basic reproduction ratio of the disease at time t which allows its propagation. This condition implies also from (20)-(21) that d I (t ) / d R (t ) =
β ( t ) S (t ) − 1≥ 0 ; ∀ t ∈ R 0 + γ ( t)
so that the infected population is a monotone increasing function with respect to the immune one. Since the above derivative is strictly positive on the time interval IT then the infected population is a strictly monotone increasing function with respect to the immune on the time interval IT. 5) The susceptible population can be calculated through time independent of the other populations as follows. One gets combining (19) and (21) by using (23): β (t) β (t) d S (t ) / S (t ) = − dR(t ) = − R (t )dt γ ( t) γ (t) −∫ 0 β (τ ) I ( τ ) d τ t
⇒ S( t ) = S( 0 ) e
(
= S( 0 ) e
)
τ ∫ β ( τ´ ) S ( τ´)−γ ( τ´ d τ´ t − I ( 0 ) ∫ 0 β (τ ) e 0 dτ
; ∀ t ∈ R 0+
5. SEIR epidemic model
The following SEIR- model distinguishes as two separate populations the “infected” E(t) which do not still have external disease symptoms from the “infectious” I(t) (also so-called “ infective” ) which exhibit already such symptoms. Let S (t) be the “susceptible” population of infection at time t, E (t) the “infected” (i.e. those which incubate the illness but do not still have any symptoms) at time t, I (t ) is the “ infectious” (or “infective”) population at time t, and R (t) is the “removed by immunity ” (or “ immune”) population at time t. Consider the SEIR-type epidemic model: S (t ) I (t ) S ( t ) = − μ S (t ) + ω R (t ) − β +ν N (t ) (1 − V (t ) ) N (t )
S (t ) I (t ) − ( μ + σ ) E (t ) E ( t ) = β N (t ) I( t ) = − (μ + γ ) I (t ) + σ E (t )
R ( t )= − ( μ + ω ) R (t ) + γ (1 − ρ ) I ( t ) + ν N (t )V (t )
(45) (46) (47)
(48) subject to initial conditions S 0 = S ( 0 ) ≥ 0 , E 0 = E ( 0 ) ≥ 0 , I 0 = I ( 0 ) ≥ 0 and R 0 = R ( 0 ) ≥ 0 under the vaccination function V : R 0 + → R 0 + . The vaccination control is either the vaccination function itself or some appropriate four dimensional vector depending on it defined “ad –hoc” for some obtained equivalent representation of the SEIR- model as a dynamic system. In the above SEIR – model, N(t) is the total population, μ is the rate of deaths from causes unrelated to the infection, ω is the rate of losing immunity, β is the transmission constant (with the total number of infections per unity of time at time t S (t ) I (t ) ), σ −1 and γ −1 are, respectively, the average durations of the latent and infective being β N (t ) periods. All the above parameters are nonnegative. The parameter ω is that of rate of immunity lost since it makes the susceptible to increase and then the immune to decrease. The usual simplified SEIR- model is obtained with ν = μ and ρ = 0 . In that case N (t ) = S ( t ) + E (t ) + I (t ) + R ( t )
30
Sen et al. / Physics 22 (2011) 20 – 39 DeM. la De Sen/laPhysics Procedia 00Procedia (2011) 000–000
= μ (N ( t ) − S ( t ) − E (t ) − I (t ) − R ( t ) ) = 0 ; ∀ t ∈ R 0+ ⇒ N (t ):= S ( t ) + E (t ) + I (t ) + R ( t ) = N ( 0 ) = N 0 = N > 0
If ν > μ then the new-born lost of maternal immunity is considered in the model. If ν < μ then there is a considered mortality incidence by external causes to the illness. The parameter ρ ∈ ( 0 ,1] is the per- capita probability of dying from the infection. If either ν ≠ μ and ρ = 0 or ν = μ and ρ ≠ 0 , and otherwise, (ν − μ ) N (t ) occurs eventually on a set of zero measure only then the total population varies I (t ) = ργ through time as obtained by correspondingly summing- up both sides of (45)-(48). Furthermore, (45) and (48) and (46) and (47) might be separately summed up to obtain the evolution dynamics of the separate populations of joint susceptible and immune and joint infected and infectious. This leads to: N (t ) = (ν − μ ) N (t ) − ρ γ I (t ) (49) ⎛ S (t ) ⎞ S ( t ) + R ( t ) = − μ ( S (t ) + R( t ) ) + ⎜⎜ γ (1 − ρ ) − β ⎟⎟ I (t ) +ν N (t ) N (t ) ⎠
⎝
⎛ S (t ) ⎞ ⎟ I (t ) E ( t ) + I (t ) = − μ ( E (t ) + I ( t ) ) − ⎜⎜ γ − β N (t ) ⎟⎠ ⎝ Note that (49) is identically zero if ν − μ = ρ = 0 ,and
(50) (51)
t N ( t ) = e (ν −μ )t N ( 0 ) − ρ γ ∫ 0 e (ν −μ )(t − τ )I ( τ ) d τ
S (t ) + R( t ) = e − μ t ( S ( 0 ) + R( 0 ) ) ⎛ ⎛ S (τ ) ⎞ t ⎟ I (τ + ∫ 0 e − μ ( t −τ ) ⎜ ν N ( τ ) + ⎜⎜ γ (1 − ρ ) − β ⎜ N (τ ) ⎟⎠ ⎝ ⎝
⎞
)⎟⎟ dτ
(52)
⎠
E (t ) + I ( t ) = e − μ t ( E ( 0 ) + I ( 0 ) ) ⎛ S (τ ) ⎞ t ⎟ I (τ ) d τ − ∫ 0 e − μ ( t −τ ) ⎜⎜ γ − β N (τ ) ⎟⎠ ⎝ In order to further solve (52), an integration by parts is performed as follows:
(53)
t t ∫ 0 p (τ ) d q(t , τ ) = ∫ 0 p (τ ) q (t , τ ) dτ t t t = ∫ 0 N ( τ ) e − μ ( t −τ )d τ = N ( τ ) q ( t ,τ ) ] − ∫ 0 q ( t ,τ ) N ( τ ) d τ 0
∂ q (t , τ ) , and ∂t e − μ ( t −τ ) t q ( t ) : = ∫ 0 e − μ ( t −τ ) d τ = μ
(54)
where q (t , τ ) =
t
= q ( t ,τ )
0
t 0
=
1− e − μ t
μ
= q (t , t ) − q (t , 0 )
(55)
so that q ( t , t ) = 1 / μ and q ( t , 0 ) = e − μ t / μ and then using (49) in (54) yields
(
) − μ1 ∫ 0t e − μ (t −τ )((ν − μ ) N (τ ) − ρ γ I (τ ) ) d τ
t − μ ( t −τ ) dτ = N ( t ) − e − μ t N( 0 ) ∫0 N( τ )e μ
1
which, after grouping identical terms, leads to t − μ ( t −τ ) dτ ∫0 N( τ ) e
=
1 ⎛ t ⎜ N ( t ) − e − μ t N ( 0 ) + ρ γ ∫ 0 e − μ ( t −τ ) I (τ ) d τ ⎞⎟ ⎠
ν ⎝
(56)
31
M. De la Sen/ Sen etPhysics al. / Physics Procedia 22 (2011) 20 – 39 Procedia 00 (2011) 000–000
Thus, one gets:
S ( t ) + R( t ) − N ( t ) = − ( E( t ) + I ( t ) ) = e − μ t ⎛⎜ S ( 0 ) + R( 0 ) − N ( 0 ) + ∫ t e μ τ ⎛⎜ γ − β S (τ ) ⎞⎟ I (τ )d τ ⎞⎟ ⎜ ⎟ 0 ⎜ ⎟ ⎝
⎝
⎛ ⎛ S (τ ) t = − e − μ t ⎜⎜ E ( 0 ) + I ( 0 ) − ∫ 0 e μ τ ⎜⎜ γ − β N (τ ) ⎝ ⎝
⎞ ⎞ ⎟⎟ I (τ ) d τ ⎟ ⎟ ⎠ ⎠
N (τ ) ⎠
⎠
(57)
6. Vaccination control
If the control objective S ( t ) = γ N ( t ) / β for all time is achieved with a positive vaccination control in [ 0 , 1 ] , it is proven below that the whole population converges exponentially to the sum of the susceptible population plus the immune population while both the infectious and infective converge exponentially to zero. This is theoretically the ideal objective since the infection is collapsing as time increases while the susceptible plus the immune populations are approximately integrating the whole population for large time. Other alternative objective has been that the immune population be the whole one but this is a more restrictive practical objective since the whole susceptible population should asymptotically track the immune one even those of the susceptible who are not contacting the disease. Theorem 1. Assume that β > γ ≥ 0 and that the vaccination function is such that S ( t ) = γ N ( t ) / β ; ∀ t ∈ R 0 + with a vaccination control in [ 0 , 1 ] for all time. Then, the SEIR model (45)-(49) is positive for all time. Furthermore,
S ( t ) + R( t ) − N ( t ) = − ( E( t ) + I ( t ) ) = e − μ t ( S ( 0 ) + R( 0 ) − N ( 0 ) ) = − e − μ t ( E ( 0 ) + I ( 0 ) )
(58) for all time what implies the following constraint for the initial conditions: γ N( 0 ) γ ( E( 0 ) + I ( 0 ) + R( 0 ) ) S( 0 ) = = β β −γ As a result, β −γ R( t ) = N ( t ) − S ( t ) − e − μ t ( E ( 0 ) + I ( 0 ) ) = N ( t ) − e − μ t ( E (0 ) + I ( 0 ) ) β ⎛ ⎞ β −γ β −γ β−γ N ( t ) + e − μ t ⎜⎜ R ( 0 ) − S ( 0 ) ⎟⎟ ≤ N ( t ) ; ∀ t ∈ R 0+ β γ β ⎝ ⎠ β −γ and R( t ) → N ( t ) as t → ∞ .Furthermore, the following two limits exist: β lim ( S (t ) + R( t ) − N ( t ) ) = lim ( E( t ) + I ( t ) ) = 0 =
t →∞
t →∞
(59)
If , in addition, ν − μ = ρ = 0 then
N ( t ) = N ( 0 ) = N = lim ( S (t ) + R( t ) ) ; lim E( t ) = lim I ( t ) = 0 t →∞
t →∞
t →∞
(60)
Proof: The mathematical SEIR- model (45)-(49) is positive since the vaccination control is in [ 0 , 1 ] for all time so that no population takes negative values at any time. On the other hand, Eqs. 58-59 follow directly from (57) and S ( t ) = γ N ( t ) / β for all time. Eqs. 60 follow from (58)-(59) since ν − μ = ρ = 0
imply N(t) ≡ N(0) , N ( t ) ≡ 0 from (49).
An associate stability result follows: Theorem 2. Assume that ρ γ ≥ 0 .Then, the following properties hold:
32
Sen et al. / Physics 22 (2011) 20 – 39 DeM. la De Sen/laPhysics Procedia 00Procedia (2011) 000–000
(i) The SEIR model is globally stable if 0 ≤ ν ≤ μ and the vaccination law fulfils V : R 0 + → [ 0 , 1 ] . (ii) If S ( t ) = γ N ( t ) / β and ν > μ ≥ 0 then the following conditions are jointly necessary for global stability under Theorem 1 N ( 0 ) = ρ γ ∫ 0 e (ν − μ )τ I (τ ) d τ , lim N ( t ) = 0 ∞
μ < ν < μ + ρ γ ; ργ > 0
(ν − μ ) N ( t ) ;
t →∞
∀t ≥ t 0 ( finite )∈ R 0 + then global stability of the SEIRργ model (45)-(48) is guaranteed if V : R 0 + → [ 0 , 1 ] . If ν > μ ≥ 0 , V : R 0 + → [ 0 , 1 ] and
(iii) If ν > μ ≥ 0 and I (t ) =
I (t ) =
(ν − μ ) N ( t ) ργ
is replaced with the weaker condition I (t ) −
(ν − μ ) N ( t ) ργ
(
)
= o e − α t for some
α ∈ R + then the SEIR – model (45) –(48) is globally stable.
Proof: (i) If 0 ≤ ν ≤ μ and ρ γ ≥ 0 then N (t ) = (ν − μ ) N (t ) − ρ γ I (t ) ≤ (ν − μ ) N (t ) ≤ 0 ; ∀ t ∈ R 0 + so
that N ( t ) ≤ N ( 0 ) < ∞ ; ∀ t ∈ R 0 + . Since the SEIR – model is positive if V : R 0 + → [ 0 , 1 ] then all the populations are nonnegative and upper-bounded by N(0). (ii) On the other hand, the solution of (49) for any initial conditions is t N ( t ) = e (ν −μ )t ⎛⎜ N ( 0 ) − ρ γ ∫0 e −(ν −μ )τ I (τ ) d τ ⎞⎟ ⎠ ⎝
which is uniformly bounded for all time only if N ( 0 ) = ρ γ ∫0 e −(ν −μ )τ I (τ ) d τ since ν > μ ≥ 0 . Also, ∞
N ( t ) < ∞ ; ∀ t ∈ R 0 + only if N ( t ) ≤ 0 on a non-necessarily connected set of infinite Lebesgue measure. Thus, there is a finite sufficiently large finite time “ t” such that : (ν − μ ) N ( t ) = (ν − μ ) (S ( t ) + E( t ) + I ( t ) + R( t ) ) I (t ) ≥ ργ ργ ⎛ ν−μ⎞ ν−μ ⎟⎟ I ( t ) ≥ ⇔ ⎜⎜ 1 − ( S ( t ) + E( t ) + R( t ) ) ⇔ I ( t ) ≥ ν − μ ( S ( t ) + E( t ) + R( t ) ) ρ γ ρ γ μ + ρ γ −ν ⎝ ⎠ which requires the parametrical conditions ρ γ > 0 and μ < ν < μ + ρ γ . Since I(t) is of exponential
order of at most (− μ ) from Theorem 1 [Eq. (58)] then (S ( t ) + E( t ) + R( t )) is also of exponential of order of at most (− μ ) so that N(t) extinguishes exponentially as they do all the populations of susceptible, infected , infectious and immune. (ν − μ ) N ( t ) with ν > μ after some finite time t then N ( t ) = N ( t ) < ∞ ; ∀ t ≥ t (iii) If I (t ) = 0 0 0 ργ and the SEIR -model is positive since V : R 0 + → [ 0 , 1 ] . I (t ) −
(ν − μ ) N ( t ) ργ
(
)
Thus, global stability follows. If
= o e − α t replaces the above stronger condition I (t ) =
(ν − μ ) N ( t ) ργ
after a finite
time then N ( t ) is of exponential order (− α ) so that N ( t ) is uniformly bounded for all time and the global stability still holds.
Note that the case ν > μ is not feasible in practice for ρ γ = 0 since the population diverges . If ρ γ > 0 , it requires a collapsing effect of the illness on the population which is also unfeasible in practical situations. It is now discussed how the vaccination law is generated to keep simultaneously the SEIR- model positivity plus the tracking objective of Theorem 1 which requires positivity. The tracking
33
M. De la Sen/ Sen etPhysics al. / Physics Procedia 22 (2011) 20 – 39 Procedia 00 (2011) 000–000
objective S ( t ) = γ N ( t ) / β for all time is equivalent for all time to any of the subsequent equivalent identities below: N ( t ) = γ N ( t ) / β + E( t ) + I ( t ) + R( t ) β −γ β ⇔ ⎛⎜ β ⎞⎟ N ( t ) = E( t ) + I ( t ) + R( t ) ( E( t ) + I ( t ) + R( t ) ) ⇔ N( t ) = ⎝ ⎠ β −γ β −γ ⇔ R( t ) = N ( t ) − E( t ) − I ( t ) (61)
β
which requires as necessary condition β > γ ≥ 0 . Although unrelated to the physical problem at hand, the necessary condition will be also accomplished with β < 0 and γ ≤ 0 with S ( t ) = γ N ( t ) / β . From Theorem 1, Eqs. 50-51 imply that lim ( S (t ) + R( t ) − N ( t ) ) = lim ( E( t ) + I ( t ) ) = 0 The solution of (51) is: t →∞
t →∞
⎤ ⎛ S (τ ) ⎞ ⎜⎜ γ − β ⎟⎟ I (τ ) d τ ⎥ (62) N (τ ) ⎠ ⎢⎣ ⎥⎦ ⎝ Then, the solution of (48) matches (61) for all time if and only if: ⎡ ⎤ β −γ β −γ ⎛ S (τ ) ⎞ R( t ) = N ( t ) − E( t ) − I ( t ) = N ( t ) − e − μ t ⎢ E (0 ) + I ( 0 ) − ∫ t e μ τ ⎜⎜ γ − β ⎟⎟ I (τ ) d τ ⎥ 0 β β N (τ ) ⎠ ⎢⎣ ⎥⎦ ⎝ ⎡
E ( t ) + I (t ) = e − μ t ⎢ E (0 ) + I ( 0 ) − ∫0t e
μτ
= e − (μ + ω ) t × ⎛⎜ R( 0 ) + ∫ 0t e ( μ +ω )τ (γ (1 − ρ ) I ( τ ) + ν N (τ )V (τ ) ) d τ ⎞⎟ ⎠ ⎝
(63)
Define an everywhere time- differentiable auxiliary function h : R 0 + → R defined as h (t ) = h ( 0 ) + ∫ 0 (γ (1 − ρ ) I ( τ ) + ν N (τ )V (τ )) dτ t
such that h (t ) = γ (1 − ρ ) I ( t ) + ν N (t )V (t ) ⇔ V (t ) =
(
1
)
(64) h ( t ) − γ (1 − ρ ) I ( t ) ν N( t ) for all time so that the last right-hand – side additive term in (63) becomes after integration by parts: t e − ( μ +ω )t e ( μ +ω )τ h (τ ) d τ = e − ( μ +ω ) t ⎛⎜ e ( μ +ω )t h (t ) − h ( 0 )− ( μ +ω ) t e ( μ +ω )τ h (τ ) d τ ⎞⎟ (65)
∫0
∫0
⎝
The replacement of (65) into (63) yields: ⎡ β − γ ( μ +ω )t e N ( t ) − e ω t ⎢ E (0 )+ I ( 0 ) − ∫0t e μ τ β ⎣⎢
⎠
⎤ ⎛ S (τ ) ⎞ ⎜⎜ γ − β ⎟⎟ I (τ )d τ ⎥ N (τ ) ⎠ ⎥⎦ ⎝
= R( 0 ) + ∫ 0 e ( μ +ω )τ ( γ (1 − ρ ) I ( τ ) + ν N (τ )V (τ ) ) d τ t
= R( 0 ) + e ( μ +ω ) t h (t ) − h ( 0 ) − ( μ +ω ) ∫ 0 e ( μ +ω )τ h (τ ) d τ t
(66)
and equivalently, and since S ( t ) = γ N ( t ) / β for all time: h (t ) =
t − ( μ +ω )( t − τ ) β −γ h (τ ) dτ N ( t ) + e −( μ +ω ) t ( h ( 0 ) − R( 0 ) ) + ( μ +ω ) ∫ 0 e β
⎡ t μτ − e − μ t ⎢ E (0 ) + I ( 0 ) − ∫0 e ⎣⎢
=
⎤ ⎛ S (τ ) ⎞ ⎜⎜ γ − β ⎟⎟ I (τ ) d τ ⎥ N (τ ) ⎠ ⎥⎦ ⎝
(
β −γ t N ( t ) + ( μ +ω ) ∫ 0 e −( μ +ω )( t − τ ) h (τ ) dτ + e − μ t e −ω t ( h ( 0 ) − R( 0 ) ) − E (0 ) − I ( 0 ) β
)
(67)
34
Sen et al. / Physics 22 (2011) 20 – 39 DeM. la De Sen/laPhysics Procedia 00Procedia (2011) 000–000
generated from:
β −γ h (t ) = β
[(ν − μ ) N ( t ) − ρ γ I ( t )] − ( μ +ω )e −( μ +ω ) t ( h (0 ) − R( 0 ) ) + μ e −μ t (E (0 ) + I ( 0 ) )
t − ( μ +ω ) 2 ∫ 0 e −( μ +ω )(t − τ ) h (τ ) dτ + ( μ +ω ) h ( t )
(68)
so that
(β − γ )(ν − μ ) N ( t ) + γ ⎛⎜ γ ρ − 1⎞⎟ I ( t ) − ( μ +ω )e −( μ +ω ) t ( h ( 0 ) − R( 0 ) ) h ( t ) − γ (1 − ρ ) I ( t ) = ⎜ β ⎟ β ⎝ ⎠
t + μ e − μ t (E (0 ) + I ( 0 ) ) − ( μ +ω ) 2 ∫ 0 e −( μ +ω )(t −τ ) h (τ ) dτ + ( μ +ω ) h ( t )
(69)
The vaccination law which ensures the positivity of the mathematical SEIR- model (45)–(48) is generated as the subsequent saturation rule: ⎧ V ( t ) if V ( t ) ∈ [ 0 , 1 ] ⎪ V( t ) = ⎨ 1 if V ( t ) > 1 ⎪ 0 if V ( t ) < 0 ⎩
(70)
where V (t ) = h ( t ) − γ (1 − ρ ) I ( t ) Define the indicator function i(t) as follows:
(71)
i( t ) = 0 if V ( t ) ∈ [ 0 , 1 ] and i( t ) = 1 , otherwise Then, one has instead of (57)
(72)
S ( t ) + R( t ) − N ( t ) = − ( E( t ) + I ( t ) ) = e − μ t
⎞ ⎛ ⎞ ⎛ ⎜ S ( 0 ) + R( 0 ) − N ( 0 ) + ∫ t e μ τ ⎜ γ − β S (τ ) ⎟ I (τ ) i ( τ )d τ ⎟ ⎟ ⎜ 0 ⎟ ⎜ N (τ ) ⎠ ⎝ ⎠ ⎝ ⎛
= − e − μ t ⎜ E ( 0 ) + I ( 0 ) − ∫ 0t e ⎜ ⎝
⎞ S (τ ) ⎞ ⎟ ⎜ γ − β N (τ ) ⎟⎟ I (τ ) i (τ )d τ ⎟ ⎝ ⎠ ⎠
μ τ ⎛⎜
(73)
which coincides with (57) for all time if the indicator function is identically zero, that is, if h ( t ) is such that the auxiliary vaccination law (71) is in [ 0 , 1 ] for all time. Also, for any given real ε > 0 and T = T ( ε ) such that ⎛ N ( 0 ) − S ( 0 ) − R( 0 ) ⎞ ln ⎜ ⎟ ε ⎝ ⎠ and ∀ t ≥ T , one gets from (57) that T=
1
μ
⎛ S (τ ) ⎞ t − γ ⎟⎟ I (τ ) i ( τ ) d τ N ( t ) − S ( t ) − R( t ) ≤ ε + ∫ 0 e − μ (t −τ ) ⎜⎜ β (74) ⎝ N (τ ) ⎠ and the right-hand-side integral takes into account the tracking deterioration if there is a time interval of nonzero Lebesgue measure such that V ( t ) ≠ V ( t ) ; ∀ t ∈ R 0 + . The following result is important to
discuss stability when the vaccination law V ( t )∈ [ 0 , 1] but it is not identically equal to V ( t ) . In fact the positivity part of Theorem 1 still holds since the SEIR- model is positive since V ( t )∈ [ 0 , 1] ; ∀ t ∈ R 0 + and the whole population evolution is independent of the vaccination law according to (49).
35
M. De la Sen/ Sen etPhysics al. / Physics Procedia 22 (2011) 20 – 39 Procedia 00 (2011) 000–000
However, the whole susceptible plus immune does not asymptotically track the whole population. In summary, one has: Theorem 3. The vaccination law (68), (70)-(71) makes the SEIR – model (1-(4) positive and globally stable under Theorem 2. Furthermore ⎛ ⎞ N ( t ) − S ( t ) − R( t ) ≤ lim sup ∫ t e − μ (t −τ ) ⎜ β S (τ ) − γ ⎟ I (τ ) i ( τ ) d τ as t → ∞ ⎜ N (τ ) ⎟ 0 t →∞ ⎝ ⎠
An alternative vaccination law replaces the saturating function (70), subject to (68), by: 1 K N (t ) N ( t ) + K I (t ) I (t ) + K R ( t ) R * (t ) + K R d ( t ) R * (t ) V (t ) = ν N(t )
(
)
R * (t )= h (t ) N ( t ) ; R * ( 0 ) = R( 0 )
(75) (76)
7. Simulation results
This section illustrates through simulation examples the theoretical results stated in the previous Sections 5 and 6 for the SEIR controlled system. Notice that the SEIR model is the most general one being the SI, SR and SIR models, described in Sections 2, 3 and 4 respectively, reduced versions of this. Moreover, the SEIR model also experiences the richest dynamic behaviour and, for this reason, it is the model considered in the numerical examples. The first example in Section 7.1 is concerned with the saturated vaccination law described by (70) while the second one in Section 7.2 is related to the unsaturated modified vaccination law introduced in Eq.75 and close equations. The SEIR model is described by the following parameters: 1 1 1 = 25550 days, = 1.9 days, = 10 days, ρ = 0.15 , 1 = 12775 days β = 1.5 days −1 (77) μ σ ω ν while γ = σ . The initial conditions are given by, S(0) = 400, E(0) = 150, I(0) = 250 and R(0) = 200 individuals so that the total population at initial time is N(0) = 1000 individuals. The function h(t) for Eq. 1 (76) is defined by h(t ) = e −ct R (0) + N (t ) 1 − e −ct with 1 / c = 5 days. Notice that h(t ) → 1 as N (t ) t → ∞ and, therefore, R * ( t ) → N ( t ) . The evolution of the model without vaccination is represented in Figure 1 in order to compare this evolution with the ones associated with different vaccination policies.
(
(
))
500
450
400
R
Populations
350
300
S 250
200
I
150
E
100
50
0
5
10
15
20
25 time (days)
30
35
40
45
50
Figure 1. Evolution of populations without vaccination
Figure 1 displays the evolution of all the partial populations through time. In particular, the numerical simulation shows a number of infected and infectious individuals, 62 of each, which correspond in total to
36
Sen et al. / Physics 22 (2011) 20 – 39 DeM. la De Sen/laPhysics Procedia 00Procedia (2011) 000–000
a 19.03% of the total population at the end of the simulation. The vaccination policies introduced in this work are implemented in order to reduce this percentage of infected and infectious population. 7.1 Saturated vaccination law Firstly, we consider the saturated vaccination policy whose basic feature is that the vaccination effort is restricted to the interval [0, 1]. The nonzero controller parameters are selected as K r = K rd = 1 . The results are shown in Figures 2 and 3. 500
450
400
R
Populations
350
S
300
250
200
150
I
100
50
E 0
5
10
15
20
25 time (days)
30
35
40
45
50
Figure 2. Evolution of the populations with saturated vaccination Vaccination law
1
Populations
0.8
0.6
0.4
0.2
0
0
5
10
15
20
25 time (days)
30
35
40
45
50
Figure 3. Vaccination law for the saturated case
Figures 1 and 2 point out the similarity between the vaccination-free and the saturated vaccination policies: for both of them, the populations reach a steady-state where there exists a non-zero value of infected and infectious population. In particular, each population of infective and infectious is now of 61 individuals which make together a 19% of the total population. It is appreciated that the vaccination effort is not large enough to eradicate the illness as it is restricted to the interval [0, 1]. In fact, Figure 3 shows that the control law is always saturated to its largest value trying to cope with the infectious and infective populations. 7.2 Unsaturated modified vaccination law In this subsection, the saturation to unity of the previous vaccination law is removed under the restriction of all the populations being nonnegative for all time. The results are depicted in Figures 4, 5 and 6 below:
M. De la Sen/ Sen etPhysics al. / Physics Procedia 22 (2011) 20 – 39 Procedia 00 (2011) 000–000
900
800
R 700
Populations
600
500
400
300
200
I
100
S E
0
0
5
10
15
20
25 time (days)
30
35
40
45
50
Figure 4. Evolution of the populations under unsaturated vaccination 1000
N 900
800
Populations
700
R* R
600
500
400
300
200
100
0
5
10
15
20
25 time (days)
30
35
40
45
50
Figure 5. Reference immune trajectory tracking and convergence to the total population Unsaturated vaccination law
1400
1200
Vaccination, V
1000
800
600
400
200
0
0
5
10
15
20
25 time (days)
30
35
40
45
50
Figure 6. Unsaturated vaccination law
Figure 4 shows that the infectious and infective populations tend to zero when the vaccination is not saturated to unity. Thus, the vaccination law depicted in Figure 6 is powerful enough to eradicate the illness and make all the population immune. Furthermore, Figure 5 shows how the proposed control law is able to make the immune population perfectly track the desired one converging to the total population N. Moreover, Figure 6 also explains why the saturated control law is not enough to eradicate the illness since the unity value of the saturated vaccination function is very small in comparison with the vaccination depicted in Figure 6. Acknowledgements
37
38
Sen et al. / Physics 22 (2011) 20 – 39 DeM. la De Sen/laPhysics Procedia 00Procedia (2011) 000–000
The authors thank the Spanish Ministry of Education by its support of this work through Grant DPI2009-07197 and to the Basque Government by its support through Grants IT378-10 and SAIOTEK SPE09UN12.
References [1] M. De la Sen and S. Alonso-Quesada, “A control theory point of view on Beverton-Holt equation in population dynamics and some of its generalizations”, Applied Mathematics and Computation, Vol. 199, No. 2, pp. 464-481, 2008. [2] M. De la Sen and S. Alonso-Quesada, “Control issues for the Beverton-Holt equation in population in ecology by locally monitoring the environment carrying capacity: Non-adaptive and adaptive cases” , Applied Mathematics and Computation, Vol. 215, No. 7, pp. 2616-2633, 2009. [3] M. De la Sen and S. Alonso-Quesada, “Model-matching-based control of the Beverton-Holt equation in Ecology”, Discrete Dynamics in Nature and Society, Article number 793512, 2008. [4] M. De la Sen, “About the properties of a modified generalized Beverton-Holt equation in ecology models”, Discrete Dynamics in Nature and Society, Article number 592950, 2008. [5] M. De la Sen, “The generalized Beverton- Holt equation and the control of populations”, Applied Mathematical Modelling, Vol. 32, No. 11, pp. 2312-2328, 2008. [6] Epidemic Models: Their Structure and Relation to Data, Publications of the Newton Institute, Cambridge University Press, Denis Mollison Editor, 2003. [7] M. J. Keeling and P. Rohani, Modeling Infectious Diseases in Humans and Animals, Princeton University Press, Princeton and Oxford, 2008. [8] M. De la Sen, A. Ibeas and S. Alonso-Quesada, “Feedback linearization-based vaccination control strategies for true-mass action type SEIR epidemic models”, Nonlinear Analysis: Modelling and Control, Vol. 16, No. 3, pp. 283-314, 2011. [9] M. De la Sen, R.P. Agarwal, A. Ibeas and S. Alonso-Quesada, “On a generalized time-varying SEIR epidemic model with mixed point and distributed time-varying delays and combined regular and impulsive controls”, Advances in Difference Equations, Vol. 2010, Article ID 281612, 2010.doi:10.1155/2010/281612. [10] M. De la Sen, R.P. Agarwal, A. Ibeas and S. Alonso-Quesada, “On the existence of equilibrium points, boundedness, oscillating behaviour and positivity of SVEIRS epidemic model under constant ands impulsive vaccination”, Advances in Difference Equations, Vol. 2011, Article ID 748608, 2010.doi:10.1155/2011/748698. [11] H. Khan , R.N. Mohapatra, K. Vajravelu and S, J. Liao, “The explicit series solution of SIR and SIS epidemic models”, Applied Mathematics and Computation, Vol. 215, No. 2 , pp. 653-669, 2009. [12] X. Y. Song, Y. Jiang and H.M. Wei , “Analysis of a saturation incidence SVEIRS epidemic model with pulse and two time delays”, Applied Mathematics and Computation, Vol. 214, No. 2 , pp. 381-390, 2009. [13] T.L. Zhang, J.L. Liu and Z.D. Teng, “Dynamic behaviour for a nonautonomous SIRS epidemic model with distributed delays”, Applied Mathematics and Computation, Vol. 214, No. 2 , pp. 624-631, 2009. [14] B. Mukhopadhyay and R. Bhattacharyya, “Existence of epidemic waves in a disease transmission model with two- habitat population”, International Journal of Systems Science, Vol. 38, No. 9 , pp. 699-707, 2007. [15] M. De la Sen, “On positivity and stability of a class of time-delay systems” , Nonlinear Analysis- Real World Applications, Vol. 8, No. 3, pp. 749-768, 2007. [16] M. De la Sen, “On positivity of singular regular linear time-delay time-invariant systems subject to multiple internal and external incommensurate point delays”, Applied Mathematics and Computation, Vol. 190, No. 1, pp. 382-401, 2007. [17] M. De la Sen, “On impulsive time-varying systems with unbounded time-varying point delays: stability and compacteness of the relevant operators mapping the input space into the state and output spaces”, Rocky Mountain Journal of Mathematics, Vol. 37, No. 1, pp. 79-129, 2007. [18] S. Alonso-Quesada and M. de la Sen, “Robust adaptive control of discrete nominally stabilizable plants”, Applied Mathematics and Computation, Vol. 150, No. 2, pp. 555-583, 2004. [19] M. De la Sen and S. Alonso-Quesada, ´´A simple vaccination control strategy for the SEIR epidemic model´´, Proceedings of the 5th IEEE International Conference on Management of Innovation and Technology, pp. 1037-1044, 2010. [20] D.J. Dalley and J. Gani, Epidemic Modelling. An Introduction, Cambridge University Press, Cambridge, 1999.
M. De la Sen/ Sen et al. / Physics Procedia 22 (2011) 20 – 39 Physics Procedia 00 (2011) 000–000
[21] J. Zhang and J. Zhen, “The analysis of epidemic network model with infectious force in latent and infected period”, Discrete Dynamics in Nature and Society, Vol. 2010, Article ID 604329, 12 pages, 2010.doi: 1155/2010/604329. [22] X. Chen, F. Chen, Q. Su and N. Zhang, “On the stability property of the infection-free equilibrium of a viral infection model”, Discrete Dynamics in Nature and Society, Vol. 2010, Article ID 262414, 9 pages, 2010.doi: 1155/2010/262414. [23] K. Gao and H. Da-yin, “Effects of immunity on global oscillations in epidemic spreading in small-world networks”, Physics Procedia, Vol. 3, Issue 5, pp. 1801-1809, 2010, doi: 10.1016/j.phpro.2010.07.022. [24] P. Sun, X.B Cao, W.D. Du and C.L. Chen, “The effect of geographical distance of epidemic spreading”, Physics Procedia, Vol. 3, Issue 5, pp. 1811-1818, 2010, doi: 10.1016/j.phpro.2010.07.023. [25] C. Y. Xia, S.W. Wu, Z.X. Liu and Z.Q. Chen, “Influence of mobile agents on the spreading behaviour of SIS model”, Physics Procedia, Vol. 3, Issue 5, pp. 1825-1830, 2010, doi: 10.1016/j.phpro.2010.07.025. [26] H.F. Zhang and B.H. Wang, “Different methods for the threshold of epidemic on heterogeneous networks”, Physics Procedia,Vol. 3, Issue 5, pp. 1831-1837, 2010, doi: 10.1016/j.phpro.2010.07.026. [27] M. De la Sen, A. Ibeas, S. Alonso-Quesada and R. Nistal, “On the equilibrium points, boundedness and positivity of a Sveirs epidemic model under constant regular constrained vaccination”, Informatica, Vol. 22, No. 3, pp.339-370, 2011. [28] M. De la Sen, “A method for general design of positive real functions”, IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, Vol.45, No. 7, pp. 764-769, 1998. [29] M. De la Sen, “ Preserving positive realness through discretization”, Positivity, Vol. 1, Issue 1, pp. 31-45, 2002. [30] M. De la Sen and A. Bilbao- Guillerna, “Fractional hold circuits versus positive realness of discrete transfer functions”, Discrete Dynamics in Nature and Society, Vol. 2005, Issue 3, pp. 373-378, 2005. doi:10.1155/DDNS.2005.373.
39