Griepepidemie Project Modelleren C

Den Dolech 2, 5612 AZ Eindhoven Postbus 513, 5600 MB Eindhoven www.tue.nl

Auteur

Marco ten Eikelder & Hugo Bink ID (resp.): 0748594 & 0739250 Begeleider: dr. J.A.C. Resing Opdrachtgever: dr. J.S.H. van Leeuwaarden Faculteit: W&I Vakcode: 2WH03 Datum

Februari 2012 - Maart 2012

Where innovation starts

Griepepidemie Project Modelleren C

Inhoudsopgave Titel

Griepepidemie Project Modelleren C

1

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

2

Inleiding . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

3

Deterministische discrete modellen . . . . . . . . . . . . .

4

3.1

SIR model . . . . . . . . . . . . . . . . . . . . . . .

4

3.2

SIR model 2 . . . . . . . . . . . . . . . . . . . . .

6

Deterministische continue modellen . . . . . . . . . . . .

8

4.1

SIR model . . . . . . . . . . . . . . . . . . . . . . .

8

4.2

SEIRD model . . . . . . . . . . . . . . . . . . . . . 10

4.3

SIR model 2 . . . . . . . . . . . . . . . . . . . . . 13

4.4

SEIR-model (met geboorte en sterfte) . . . . . . . 17

4

5

6

7

Reproductiegetal . . . . . . . . . . . . . . . . . . . . . . . 20 5.1

SIR-model . . . . . . . . . . . . . . . . . . . . . . 20

5.2

SEIR-model . . . . . . . . . . . . . . . . . . . . . . 20

Probabilistische modellen . . . . . . . . . . . . . . . . . . 22 6.1

Eenvoudig model SI . . . . . . . . . . . . . . . . . 22

6.2

SIR model . . . . . . . . . . . . . . . . . . . . . . . 25

6.3

Het Reed frost model . . . . . . . . . . . . . . . . 33

Conclusie . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

Bibliografie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

Where innovation starts

Technische Universiteit Eindhoven University of Technology

1

Summary

Recent outbreaks of several (animal) diseases like the avian influenza (2003) in the Netherlands and SARS worldwide makes it very important to understand the spread of an infection. Mathematical models can help us to understand the dynamics of such a disease. They play an extreme important role in understanding and predicting the behavior of an infection. Results from mathematical epidemiology can and have been used in other science as well. One can think of the spread of information or the spread of computer viruses. There are many different ways to model an infection. Deterministic and stochastic models can be used. Stochastic models of epidemics are more realistic than the deterministic ones, but deterministic models can also help us to understand the behavior of an infection very well. There are two ways of analyzing the model deterministic: continuous or discrete; the time can be a continuous or a discrete variable. We have looked at many different models. A common used model for the spread of an infection is the SIR (Susceptible, Infected, Recovered) model. All particles in the model are in one of these three states. Infected particles can infect susceptible ones and infected ones can become recovered. This model can be analyzed deterministic and stochastic. When the model is continuous deterministic, we obtain a system of differential equations. Sometimes these systems can be solved analytic; it depends on the assumptions made. To solve these (linear) systems we use eigenvalues and eigenvectors. If the system is not analytically solvable, we use computer programs to simulate. Even if the system does not have an analytic solution, we can still learn something interesting from the model. For example we can plot the solutions as curves in the phase plane. A deterministic discrete model is almost the same as the continuous version. The difference is that we do not obtain a system of differential equations but we get a system of recurrence equations. Stochastic models are the most realistic and maybe the most interesting models to understand the behavior of a disease. We have used a stochastic model to predict how long it takes before a disease dies out. The basic reproduction number is (one of) the most important parameters in epidemiology. It can be used to predict if an infectious disease can invade in a population. The basic reproduction number is the number of cases one case generates on average over the course of their infectious period. If this number is smaller than one, the disease will die out. If the number is bigger than one the disease can spread out in a population. When the basic reproduction number almost equals one, it is a critical case. This is the most interesting part, which we will investigate in the second part of our survey.

2 Griepepidemie Project Modelleren C


2

Inleiding

Griep is een influenzavirus dat zich gemakkelijk verspreidt van persoon tot persoon. Enkele symptomen van griep zijn hoofdpijn, spierpijn over je gehele lichaam en koorts. Bij een griepepidemie worden veel mensen tegelijkertijd ziek. Dit kan vervelende gevolgen hebben en het is dus belangrijk zo’n epidemie zoveel mogelijk tegen te gaan. Het mathematisch modelleren van de verspreiding van een virus is erg belangrijk om de verspreiding van een virus te begrijpen. Vooral de duur en omvang van de epidemie zijn relevante gegevens. Door hier inzicht in te krijgen kunnen methoden worden ontwikkeld om de verspreiding af te remmen. Denk hierbij aan vaccinatie of quarantaine. In dit het eerste deel van het project zullen we verschillende modellen die een griepepidemie beschrijven onderzoeken. We zullen deterministische (zowel discrete als continue) modellen en stochastische modellen bekijken. Hoewel stochastische modellen realistischer zijn, zijn deterministische modellen ook interessant en nuttig. Het reproductiegetal is een belangrijk kengetal binnen de epidemiologie. Dit zullen we kort aan bod laten komen.



3

Deterministische discrete modellen

Een van de manieren om de verspreiding van een griepvirus te analyseren is door te kijken naar recurrente betrekkingen. Een van de aannames bij dit model is dat de tijd als een discrete variabele wordt gezien. We bekijken het aantal deeltjes dat in een bepaalde toestand zit op verschillende tijdstippen. De toestanden waarin de deeltjes kunnen zitten zijn: vatbaar/susceptible (S), blootgesteld/exposed (E), geïnfecteerd/infected (I), hersteld/recovered (R) en dood/dead (D). Het aantal deeltjes in een bepaalde toestand hangt onder andere af van het aantal deeltjes in die toestand op het vorige tijdstip of de vorige tijdstippen.

3.1

SIR model

We bekijken een SIR-model, dit is een model dat het aantal deeltjes in een bepaalde toestand weergeeft. We bekijken het aantal deeltjes in de toestanden vatbaar/susceptible (S), geïnfecteerd/infected (I) en hersteld/recovered (R). Een van de aannamen van het model is dat de herstelde deeltjes niet meer vatbaar zijn. Ook gaan we er vanuit dat geen nieuwe deeltjes in het model kunnen komen en geen deeltjes weg kunnen gaan (geen geboorte, sterfte en migratie). We definiëren: N = totaal aantal deeltjes Sn = het aantal vatbare deeltjes op tijdstip n In = het aantal zieke deeltjes op tijdstip n Rn = het aantal resistente deeltjes op tijdstip n. We nemen nu aan dat de verandering van het aantal deeltjes dat in een andere/de volgende toestand zal komen, proportioneel is aan het aantal deeltjes in de huidige toestand. Laat α de fractie zijn van het aantal vatbare deeltjes dat in de volgende tijdseenheid geïnfecteerd raakt, en laat β de fractie zijn van het aantal geïnfecteerde deeltjes dat in de volgende tijdseenheid hersteld is. Het stelsel recurrente betrekkingen dat bij dit model hoort is dan:  Sn+1 = Sn − αSn ∀n ∈ N      In+1 = In + αSn − In β ∀n ∈ N (1)      Rn+1 = Rn + In β ∀n ∈ N. We merken op dat Sn+1 + In+1 + Rn+1 = Sn + In + Rn ∀n ∈ N zodat het totaal aantal deeltjes binnen het model constant blijft. Het (lineaire) stelsel recurrente betrekkingen



zullen we nu oplossen. Allereerst schrijven we het stelsel in matrixvorm:      Sn+1 1−α 0 0 Sn  In+1  =  α 1 − β 0  I n  . Rn+1 0 β 1 Rn

(2)

Net als in het vorige model (SI) kunnen we het stelsel oplossen met behulp van eigenwaarden en eigenvectoren. Merk op dat de vierkante matrix een benedendriehoeksmatrix is, zodat we de eigenwaarden direct van de diagonaal kunnen aflezen: λ1 = 1 − α, α T λ2 = 1−β en λ3 = 1. Als bijbehorende eigenvectoren vinden we: v1 = (− −α+β β , − β , 1) , v2 = (0, −1, 1)T en v3 = (0, 0, 1)T . De algemene oplossing is nu   3 Sn X  In  = ci · λin · vi ∀n ∈ N, c1 , c2 , c3 ∈ R, (3) Rn 1 Met de begincondities S0   Sn         In          Rn

= S0 , I0 = 0, R0 = 0 vinden we als oplossing: = S0 (1 − α)n = −

S0 α((1 − α)n − (1 − β)n ) α−β

= −

S0 (−α + α(1 − β)n + β − (1 − α)n β) . α−β

(4)

In figuur 1 zijn voor S0 = 100, I0 = R0 = 0 en α = 0.15, β = 0.1 de grafieken van Sn , In en Rn weergegeven.

Figuur 1: Grafieken van Sn , In en Rn weergegeven



We zien in de grafiek dat op het tijdstip n = 8 de meeste deeltjes geïnfecteerd zijn. We hadden dit ook kunnen berekenen door de discrete formules als continu te beschouwen. Als we dan differentiëren krijgen we In0

S0 α((1 − α)n log(1 − α) − (1 − β)n log(1 − β)) . =− α−β

(5)

In0 gelijkstellen aan nul levert log n=−

log (1−α) log (1−β)

log (1 − α) − log (1 − β)

=

1−α 1−β

log(1 − β) log . log(1 − α)

(6)

Als we hier de waarden invullen zoals in de grafiek in figuur 1 krijgen we n ≈ 7, 58.

3.2

SIR model 2

Net als in het vorige model bekijken het aantal deeltjes in de toestanden vatbaar (S), geïnfecteerd (I) en hersteld (R). In het vorige model was de verandering van het aantal vatbare deeltjes alleen afhankelijk van het aantal vatbare deeltjes; dus niet van het aantal geïnfecteerde deeltjes. Omdat dit in de werkelijkheid niet zo zal zijn, veranderen we deze aanname. Net als in het vorige model definiëren we: N = totaal aantal deeltjes Sn = het aantal vatbare deeltjes op tijdstip n In = het aantal zieke deeltjes op tijdstip n Rn = het aantal resistente deeltjes op tijdstip n.

Stel een deeltje komt per tijdsperiode gemiddeld in contact met γ andere deeltjes. In het geval van een griepepidemie bedoelen we met dit contact dat een geïnfecteerd deeltje dicht genoeg bij een ander deeltje in de buurt komt om deze te besmetten. Natuurlijk zullen alleen de vatbare deeltjes geïnfecteerd kunnen raken. Echter zullen vatbare deeltjes niet altijd besmet raken als ze in contact komen met een geïnfecteerd deeltje. Zij τ de kans dat een vatbaar deeltje besmet raakt wanneer deze in contact komt met een geïnfecteerd deeltje. Dan is γ τ =: α˜ het gemiddeld aantal deeltjes dat een geïnfecteerd deeltje per tijdseenheid zal besmetten. Definieer nu α = Nα˜ . Omdat besmettingen echter alleen effect hebben voor de vatbare deeltjes zal ieder geïnfecteerd deeltje α˜ SNn = αSn vatbare deeltjes besmetten raken op tijdstip n. Als we verder aannemen dat ieder tijdstip een fractie β van het aantal geïnfecteerde deeltjes herstelt kunnen we een stelsel met recurrente betrekkingen opstellen:



         

Sn+1 = Sn − αSn In

        

Rn+1 = Rn + In · β

∀n ∈ N

In+1 = In + αSn In − In · β ∀n ∈ N (7) ∀n ∈ N

Sn + In + Rn = N

∀n ∈ N

Om de grafieken Sn , In en Rn te plotten voor specifieke waarden van α, β, S0 , I0 en R0 schrijven we een Java-programma. Het programma berekent uit de beginwaarden eerst S1 , I1 en R1 , vervolgens S2 , I2 en R2 etc.. In figuur 2 zijn de grafieken van Sn , In en Rn weergegeven voor α = 10−3 , β = 2 · 10−2 , S0 = 200, I0 = 1 en R0 = 0. We zien

200ææææææææææ

ææææ ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ææ ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ææ ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì ææ ì ì ì ì ì ì ì ì ì ì ì ì ì ì ì æ ì ì ì ì ì ì ì æ ì ì ì ì ì ì ì æ ì ì ì ì ì ì ì æ ì ì ì ì ì ì æ ì ì ì ì 150 ì ì ì æ ì ì ì ì ì ì æ ì ì ì ì àààààààà ì ì æ ì àà àààà ì ì ì àà ì æ àà ì àà ì ì ì àà ì æ àà ì àà ì ì àà ì ì æà àà ì ì ààì ì 100 ìààà æàà ì ì ààà ì ì àæ ààà ì ì ì ààà àæ ì ì ààà ì à æ ì ààà ì ààà ì à æ ì ààà ì àààà ì à æ ì àààà ì à ì àààà æ ì 50 àààà à ì ì ààààà æì à ààààà ì æì à àààààà ì æ à àààààà ì æ ì ààààààà à ìæ àààààààà à ì æ ì àààààààààà à ì ææ à ì àààààààààààà à ì æ ì àààààààààààààààà ì æææ ì ààààààààààààààààààà ààà ì ì à ì æ à æ ì ææææææææ ì ì àààì ì àì ì àì ì àì ì àì æææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææ ì àì ì àì ì àì ì àì àì ìì ì ì ì

50

100

150

200

æ

Sn

à

In

ì

Rn

n

Figuur 2: Grafieken van Sn , In en Rn weergegeven voor vaste α en β dat het aantal vatbare deeltjes alsmaar afneemt en uiteindelijk 0 wordt. Het aantal geïnfecteerde deeltjes neemt eerst toe en vervolgens weer af. Als we schrijven In+1 − In = αSn − β zien we dat In niet verandert als Sn = βα . Dus als Sn = βα dan bereikt In voor die n haar maximum.



4

Deterministische continue modellen

Een andere manier om een griepvirus te analyseren, is met behulp van differentiaalvergelijkingen. Het model bevat nu geen willekeurigheid; dezelfde beginsituatie zal exact hetzelfde resultaat opleveren. De tijd is nu een continue variabele in tegenstelling tot de recurrente betrekkingen. Net als bij de recurrente betrekkingen verdelen we de (vaste) groep deeltjes weer in een aantal groepen.

4.1

SIR model

In dit model bekijken we het aantal vatbare (S), geïnfecteerde (I) en herstelde (R) deeltjes. We nemen aan dat alle deeltjes in het begin vatbaar zijn. Het model berust nu op de aanname dat de verandering van het aantal deeltjes dat in een andere/de volgende toestand zal komen, proportioneel is aan het aantal deeltjes in de huidige toestand. Hiervoor definiëren we de volgende parameters: α = de intensiteit om van S naar I te gaan β = de intensiteit om van I naar R te gaan. In figuur 3 is het model weergegeven. We vinden dan het volgende stelsel differentiaal-

Figuur 3: Schematische weergave SIR model

vergelijkingen:

Merk op dat geldt:


 d S(t)   = −αS(t)   dt        d I (t) = αS(t) − β I (t) dt      d R(t)    = β I (t)    dt

(8)

d S(t) d I (t) d R(t) + + = 0, dt dt dt

(9)


dus het totaal aantal deeltjes in alle toestanden is constant! We gaan het gevonden stelsel (lineaire) differentiaalvergelijkingen exact oplossen. Uit S 0 (t) = αS(t) volgt direct S(t) = c1 · e−αt met c1 ∈ R. Als we deze oplossing invullen in de tweede differentiaalvergelijking vinden we: d I (t) = c1 · α · e−αt − β I (t). (10) dt De homogene oplossing Ih (t) van deze differentiaalvergelijking is: Ih (t) = c2 · e−βt

met c2 ∈ R.

(11)

Als particuliere oplossing proberen we een oplossing van de vorm: I p (t) = d · e−αt

met d ∈ R.

Invullen geeft: −dαe−αt = c1 αe−αt − dβe−αt ⇒ d =

(12) c1 α . β −α

(13)

De oplossing is dan: I (t) = c2 · e−βt +

c1 α · e−αt β −α

met c1 , c2 ∈ R.

Nu I (t) bekend is, kan ook R(t) gevonden worden: Z t Z t c1 α −αs −βs + R(t) = β I (s)ds = β c2 · e ·e ds β −α 0 0 c1 1 − e−tα + c2 1 − e−βt . = β −α Met de begincondities: S(0) = zodat de oplossing is:  S(t)         I (t)         R(t)

S0 en I (0) = R(0) = 0 vinden we c1 = S0 en c2 =

(14)

(15) (16) S0 α α−β ,

= S0 · e−αt =

S0 α e−αt − e−βt β −α

= S0 +

(17)

S0 βe−αt − αe−βt α−β

In figuur 4 zijn de grafieken van S(t), I (t) en R(t) weergegeven voor: S(0) = 100, I (0) = R(0) = 0, α = 0.4 en β = 0.2. We zien dat uiteindelijk alle deeltjes hersteld zijn.



100

SHtL

80

IHtL 60

RHtL 40

20

0 0

5

10

15

20

25

30

t

Figuur 4: Het aantal deeltjes in de verschillende toestanden over de tijd, bij vaste parameters

4.2

SEIRD model

We verdelen de bevolking nu in 5 categoriëen: vatbaar (S), blootgesteld (E), geïnfecteerd (I), hersteld (R) en dood (D). In figuur 5 zijn de verschillende categorieën weergegeven. We nemen aan dat de gehele bevolking in het begin vatbaar is. De toestand daarna is blootgesteld en dan geinfecteerd. Vervolgens kunnen geïnfecteerde deeltjes weer herstellen of dood gaan. Een persoon die hersteld is van het virus, kan niet opnieuw in aanraking komen met het virus. Het model berust weer op de aanname dat

Figuur 5: De verschillende categoriëen de verandering van het aantal deeltjes dat in een andere/de volgende toestand zal komen, proportioneel is aan het aantal deeltjes in de huidige toestand. We definiëren de volgende parameters: α β γ δ

= de intensiteit om van S naar E te gaan = de intensiteit om van E naar I te gaan = de intensiteit om van I naar R te gaan = de intensiteit om van I naar D te gaan.



We vinden dan het volgende stelsel differentiaalvergelijkingen:  d S(t)     dt       d E(t)      dt      d I (t)  dt       d R(t)      dt         d D(t) dt

= −αS(t) = αS(t) − β E(t) = β E(t) − (γ + δ)I (t)

(18)

= γ I (t) = δ I (t).

Merk op dat weer (net als in het SIR model) geldt: d S(t) d E(t) d I (t) d R(t) d D(t) + + + + =0 dt dt dt dt dt

(19)

zodat we vinden dat de het totaal aantal deeltjes binnen het model constant is: S(t) + E(t) + I (t) + R(t) + D(t) = Constant.

(20)

Het gevonden stelsel (18) kunnen we ook in matrixvorm schrijven, waarin we de vierkante matrix A noemen:  0      S(t) −α 0 0 0 0 S(t) S(t)  E(t)   α −β    E(t)  0 0 0      E(t)     I (t)  =  0     (21) β −(γ + δ) 0 0      I (t)  =: A  I (t)  .  R(t)   0      0 γ 0 0 R(t) R(t) D(t) 0 0 δ 0 0 D(t) D(t) Dit stelsel differentiaalvergelijkingen kunnen we oplossen m.b.v. eigenwaarden en eigenvectoren. Merk op dat de matrix A een benedendriehoeksmatrix is, zodat we de eigenwaarden direct van de diagonaal af kunnen lezen: λ1 = −α, λ2 = −β, λ3 = −(γ + δ) en λ4 = λ5 = 0. De bijbehorende eigenvectoren zijn:  (α−β)(α−γ −δ)          − 0 0 0 0 βδ  α(α−γ −δ)  β−γ −δ            δα   γ0+δ  0 0 βδ          v1 =   , v2 =  − αδ  −γ δ  , v3 = − γ δ  , v4 = 0 , v5 = 0 . (22)         1 γ 0   δ δ δ 1 0 1 1 1



De algemene oplossing van het stelsel is nu van de vorm:   S(t)  E(t)  5   X  I (t)  = ci · eλi ·t · vi , waarbij ci nader te bepalen constanten zijn.    R(t)  i=0 D(t)

(23)

Het stelsel differentiaalvergelijkingen kan ook eenvoudiger worden opgelost. Net als in het SIR model volgt uit S 0 (t) = αS(t) direct dat S(t) = c1 · e−αt met c1 ∈ R. Als we deze oplossing invullen in de tweede differentiaalvergelijking vinden we: d E(t) = c1 · α · e−αt − β E(t). dt

(24)

Op eenzelfde manier als bij het SIR model kan de oplossing van E(t) gevonden worden. Merk op dat E(t) van dit model gelijk is aan I (t) van het vorige SIR model: E(t) = c2 · e−βt +

c1 α · e−αt β −α

met c1 , c2 ∈ R.

(25)

Ook kan zo I (t) gevonden worden en vervolgens de functies R(t) en D(t). Met de begincondities: S(0) = S0 , E(0) = I (0) = R(0) = D(0) = 0 vinden we de volgende oplossing:  S(t) = S0 · e−αt         S0 α   e−αt − e−βt E(t) =    β −α        e−tα e−tβ e−t (γ +δ)    + − I (t) = S0 αβ   (−α + β)(−α + γ + δ) (−α + β)(β − γ − δ) (β − γ − δ)(−α + γ + δ)        e−tα β e−tβ α R(t) = S0 γ − − (−α + β)(−α + γ + δ) (−α + β)(β − γ − δ)        −t (γ +δ) αβ  e γ   + S0 γ + S0 ·    (γ + δ)(β − γ − δ)(−α + γ + δ) γ +δ        e−tα β e−tβ α    D(t) = S0 δ − −   (−α + β)(−α + γ + δ) (−α + β)(β − γ − δ)    −t (γ +δ) αβ  e δ   + S0 δ + S0 ·  (γ + δ)(β − γ − δ)(−α + γ + δ) γ +δ (26) 12 Griepepidemie Project Modelleren C


Uiteindelijk zijn alle deeltjes in de toestanden hersteld (R) en dood (D). Immers, ga na dat geldt: lim S(t) = lim E(t) = lim I (t) = 0. (27) t→∞

t→∞

t→∞

Het aantal deeltjes in de limiet in de toestanden hersteld (R) en dood (D) is een soort gewogen gemiddelde: lim R(t) = S0 ·

t→∞

δ γ , lim D(t) = S0 · . γ + δ t→∞ γ +δ

(28)

In figuur 6 zijn de grafieken van S(t), E(t), I (t), R(t) en D(t) weergegeven voor: S(0) = 100, E(0) = I (0) = R(0) = D(0) = 0, α = 0.4, β = 0.2, γ = 0.6 en δ = 0.2. Merk op dat 0.2 0.6 = 75 en lim D(t) = 100 · = 25. inderdaad geldt: lim R(t) = 100 · t→∞ t→∞ 0.6 + 0.2 0.6 + 0.2 SHtL

100

EHtL 80

IHtL RHtL

60

DHtL 40

20

0 0

5

10

15

20

25

30

t

Figuur 6: Het aantal deeltjes in de verschillende toestanden over de tijd, bij vaste parameters

4.3

SIR model 2

In het vorige model (SEIRD) werd aangenomen dat de verandering van het aantal deeltjes dat in een andere/de volgende toestand zal komen, proportioneel is aan het aantal deeltjes in de huidige toestand. In het model is dus de verandering van het aantal vatbare deeltjes alleen afhankelijk van het aantal vatbare deeltjes. In de werkelijkheid zal deze verandering ook afhangen van het aantal blootgestelde of geïnfecteerde deeltjes. We nemen dit effect mee in het volgende model. We bekijken nu het aantal deeltjes in de volgende 3 toestanden: vatbaar (S), geïnfecteerd (I) en hersteld (R). We definiëren dus de verandering van het aantal nieuwe infecties als α · S · I , waarbij α een infectieparameter is. Wel nemen we weer aan dat de geïnfecteerde deeltjes herstellen met een 13 Griepepidemie Project Modelleren C


gelijke kans op elk tijdstip. Het aantal nieuwe herstelde deeltjes is dan dus β I , waarin β een herstel parameter is. Het model is schematisch weergegeven in figuur 7. We

Figuur 7: Schematische weergave van het model

vinden dan het volgende stelsel differentiaalvergelijkingen:  d S(t)   = −αS(t) · I (t)   dt      d I (t) = αS(t) · I (t) − β I (t)  dt         d R(t) = β I (t). dt

(29)

Het stelsel differentiaalvergelijkingen is een niet-lineair stelsel, en heeft geen exacte algebraïsche oplossing. Toch kunnen we wel iets zeggen over het model. Ten eerste geldt net als in het vorige model: d S(t) d I (t) d R(t) + + = 0, dt dt dt

(30)

S(t) + I (t) + R(t) = Constant.

(31)

zodat we weer vinden: Met Wolfram Mathematica 7.0 kunnen we toch een benaderende grafiek maken van S(t), I (t) en R(t). In figuur 8 is het aantal vatbare, geïnfecteerde en herstelde deeltjes weergegeven over de tijd bij vaste parameters S(0) = 499, I (0) = 1, R(0) = 0, α = 10−3 en β = 2 · 10−1 .



SHtL

500

400

IHtL

300

RHtL 200

100

20

40

60

80

100

t

Figuur 8: Het aantal vatbare, geïnfecteerde en herstelde deeltjes weergegeven over de tijd bij vaste parameters We zien inderdaad dat, zoals verwacht, er op den duur geen geïnfecteerde deeltjes meer zijn; alle geïnfecteerde deeltjes zijn aan het eind van de epidemie hersteld. Merk op dat het bij dit model noodzakelijk is dat I (0) > 0. Immers als I (0) = 0, dan is d S(t) d I (t) d R(t) = = = 0. Alle deeltjes blijven voor altijd vatbaar. Dit zou een dt dt dt onrealistisch en niet interessant model zijn. Het aantal geïnfecteerde deeltjes verandert niet ( d Idt(t) = 0) als S(t) = βα . In figuur 8 is te zien dat het aantal vatbare deeltjes inderdaad gelijk aan

β α

=

2·10−1 10−3

= 200 als het aantal geïnfecteerde deeltjes maximaal is.

Oplossingen als banen in het fasevlak We bekijken de oplossingen van dit SIR-model als banen in het fasevlak. Merk op dat de differentiaalvergelijkingen van S(t) en I (t) onafhankelijk zijn van die van R(t). Omdat het aantal deeltjes in het model constant blijft, kunnen we R(t) vinden uit: S(t)+ I (t)+ R(t) = N . We zijn geïnteresseerd hoe de oplossingsbaan zich gedraagt in S I -fasevlak (S als horizontale as en I als verticale as). We nemen aan dat er in het begin geen deeltjes in de toestand R zijn: S(0) = S0 , I (0) = I0 en R(0) = 0. Dus op t = 0 geldt S + I = N , ofwel elke oplossing begint op de lijn: S + I = N . We weten dat S + I ≤ N voor alle t, de oplossing zal dus altijd binnen de driehoek ingesloten door de lijnen: S + I = N , S = 0 en I = 0 liggen. Ook weten we dat S 0 (t) < 0 (immers α > 0), zodat de banen naar links zullen gaan. Het teken van I (t) hangt af van S(t):  >0 als S(t) > βα      0 (32) I (t) = αS(t)I (t) − β I (t) = I (t)(αS(t) − β) = 0 als S(t) = βα      <0 als S(t) < βα 15 Griepepidemie Project Modelleren C


We moeten nu twee gevallen bekijken: βα > N en βα < N ; we maken dus onderscheid of βα groter of kleiner is dan het aantal deeltjes waaruit de bevolking bestaat. Bekijk het geval: βα < N . De isoclinen (curves waarlangs het vectorveld geheel horizontaal of verticaal is) kunnen gevonden worden door te stellen: 0 S (t) = 0 (33) I 0 (t) = 0 De isoclinen zijn dus de lijnen: S = 0, I = 0 en S = βα . In figuur 9 zijn 2 verschillende banen in het SI fase vlak weergegeven (voor het geval: βα < N ). Een baan start in

Figuur 9: Het SI fase vlak met daarin 2 banen van het geval: βα < N . De ene start in het punt P en de andere in het punt Q, beide op de lijn: S + I = N het punt P. Omdat de S−coördinaat van het punt P kleiner is dan βα , is I 0 (t) < 0 zal de baan naar beneden gaan. Eerder hadden we al gevonden dat elke baan naar links gaat, zodat we vinden dat deze baan naar linksonder gaat. De andere baan die begint in het punt Q (S Q > βα ) zal vanwege analoge argumenten de baan eerst naar linksboven gaan, totdat de lijn S = βα bereikt wordt. Hierna zal de baan naar linksonder gaan. Merk op dat de S−coördinaat van het punt waar een baan de lijn I = 0 snijdt (S ∗ ), het aantal personen is dat het virus niet krijgt. We willen (S ∗ ) bepalen. Hiervoor leggen we een relatie tussen S en I . We krijgen dI dt dS dt


=

dI αS I − β I β = = −1 + , dS −αS I αS

(34)


ofwel (scheiden van variabelen) d I = −1 +

β d S. αS

Integreren aan beide kanten levert: Z Z β d I = −1 + d S, αS

(35)

(36)

waaruit volgt: I = −S +

β log S + C α

Op t = 0 geldt dus: C = N − I (S) = −S +

β α

met C een constante, C ∈ R.

(37)

log S0 , invullen in vergelijking (37) geeft:

β β β S log S + N − log S0 = N − S − log . α α α S0

(38)

Het aantal personen dat niet geïnfecteerd raakt, S ∗ , volgt dus uit de vergelijking I (S) = 0.

4.4

SEIR-model (met geboorte en sterfte)

Het is natuurlijk niet logisch dat deeltjes alleen dood kunnen gaan als ze ziek zijn (zoals we in het SEIRD-model aannamen). Daarom gaan we er in dit model van uit dat deeltjes in iedere toestand (S, E, I en R) kunnen sterven. We gaan ervan uit dat deeltjes in iedere tijdsperiode evenveel kans hebben om te overlijden, onafhankelijk van de toestand waarin ze zich bevinden. De intensiteit van de sterfte is µ. Ook dit is natuurlijk niet helemaal realistisch, maar om het model niet meteen heel ingewikkeld te maken zullen we dit toch even aannemen. Daarnaast voegen we ook geboorte toe (intensiteit λ). We nemen aan dat alle deeltjes die geboren worden in toestand S terechtkomen. Verder gelden de overgangsintensiteiten zoals aangegeven in de figuur op de volgende pagina.



Figuur 10: Schematische weergave van het SEIR-model

We kunnen nu de volgende differentiaalvergelijkingen opstellen bij dit model:  d S(t)   = −αS(t)I (t) + λ − µS(t)   dt        d E(t)   = αS(t)I (t) − (µ + β)E(t)   dt   d I (t)   = β E(t) − (γ + µ)I (t)    dt         d R(t) = γ I (t) − µR(t). dt

(39)

Het optellen van de differentiaalvergelijkingen van het stelsel levert d (S(t) + E(t) + I (t) + R(t)) = λ − µ (S(t) + E(t) + I (t) + R(t)) . dt

(40)

Als deze uitdrukking groter/kleiner dan 0 is, zal het totaal aantal deeltjes binnen het model toenemen/afnemen. In figuur 11 is voor de vaste parameters α = 10−3 , β = 2·10−2 , γ = 10−2 , λ = 50 en µ = 7·10−2 het aantal vatbare, geïnfecteerde, blootgestelde en herstelde deeltjes weergegeven. Als beginwaarde hebben we gekozen: S(0) = 5000, I (0) = 5, E(0) = R(0) = 0.



SHtL

5000

EHtL 4000

IHtL 3000

RHtL 2000

1000

50

100

150

200

250

t

Figuur 11: Het aantal vatbare, geïnfecteerde, blootgestelde en herstelde deeltjes weergegeven over de tijd bij vaste parameters



5

Reproductiegetal

Het reproductiegetal R0 geeft het verwachte aantal besmettingen van één geïnfecteerd deeltje weer. Voor het getal geldt de volgende formule: R0 = τ · c¯ · d.

(41)

In deze formule is τ de kans dat een ziekte wordt overgedragen wanneer een geïnfecteerd deeltje in contact komt met een vatbaar deeltje. In het geval van een griepepidemie bedoelen we met dit contact dat een geïnfecteerd deeltje dicht genoeg bij een ander deeltje in de buurt komt om deze te besmetten. Verder is c¯ is het gemiddeld aantal keren contact tussen geïnfecteerde en vatbare deeltjes en d de gemiddelde tijdsperiode dat een deeltje geïnfecteerd is. Wanneer breekt er nu een epidemie uit?

5.1

SIR-model

We bekijken eerst het SIR-model met differentiaalvergelijkingen zoals in vergelijking (29). Er kan een epidemie komen als het aantal geïnfecteerde deeltjes toeneemt: αS(t) α d I (t) > 0 ⇔ αS(t) · I (t) − β I (t) > 0 ⇔ αS(t) − β > 0 ⇔ > 1 ⇒ N · > 1. dt β β (42) Merk op dat we I (t) kunnen wegdelen. Immers I (t) 6= 0 omdat er anders geen epidemie meer zou zijn. Omdat β de intensiteit is om van I naar R te gaan en d de gemiddelde tijdsperiode dat een deeltje geïnfecteerd is geldt β = d −1 . In de beginsituatie is het aantal vatbaren ongeveer gelijk aan N . Bij benadering geldt dus α = R0 > 1. τ cd ¯ = N (43) β Er kan dus een epidemie uitbreken als het reproductiegetal groter is dan 1. Ook intuïtief is eenvoudig te begrijpen dat een epidemie uitbreekt als ieder geïnfecteerd deeltje gemiddeld meer dan één vatbaar deeltje besmet. Als het reproductiegetal kleiner is dan 1 zal een geïnfecteerd deeltje gemiddeld minder dan één vatbaar deeltje besmetten en zal dus geen epidemie uitbreken. Verder geldt hoe groter het reproductiegetal, des te besmettelijker is de ziekte.

5.2

SEIR-model

We gaan uit van het SEIR-model met differentiaalvergelijkingen zoals in vergelijking (39) en zullen nu uitrekenen wat de waarde van R0 is zodat we kunnen bepalen wanneer 20 Griepepidemie Project Modelleren C


een epidemie zal uitbreken. We bepalen deze waarde op de manier zoals beschreven in [4]. We zien eenvoudig in dat we een stationair punt hebben voor S(t) = µλ en E(t) = I (t) = R(t) = 0. De Jacobi-matrix van de differentiaalvergelijkingen ziet er als volgt uit:   −α I (t) − µ 0 −αS(t) 0  α I (t) −(µ + β) αS(t) 0   . (44)  0 β −(γ + µ) 0  0 0 γ −µ Als we het gevonden stationaire punt invullen in de Jacobi-matrix en vervolgens de eigenwaarde uitrekenen en kleiner dan nul stellen krijgen we −µ(−(µ + β)(γ + µ)µ + αβλ) < 0.

(45)

Als we dit uitwerken vinden we (R0 =)


αβλ > 1. µ(β + µ)(γ + µ)

(46)


6

Probabilistische modellen

Een probabilistisch model van de verspreiding van een griepvirus neemt mee dat er wat variatie en willekeurigheid zit in (tenminste één van) de parameters en variabelen. De resultaten van het model zullen dus kansverdelingen zijn.

6.1

Eenvoudig model SI

We beginnen met een eenvoudig model met slechts twee toestanden: vatbaar/susceptible (S) en geïnfecteerd/infected (I). Hierbij heeft elk deeltje in toestand S kans om geïnfecteerd te raken en dus in toestand I te komen. Allereerst bekijken we één afzonderlijk deeltje. Er wordt aangenomen dat een vatbaar deeltje met kans p op het volgende tijdstip geïnfecteerd wordt, en dus met kans 1− p nog steeds vatbaar is. Evenzo is een geïnfecteerd deeltje op het volgende tijdstip weer vatbaar met kans q en dus nog steeds geïnfecteerd met kans 1 − q. In figuur 12 is een schematische weergave van de toestanden weergegeven met de bijbehorende overgangskansen. Merk op dat de volgende toestand van een willekeu-

Figuur 12: Schematische weergave van de toestanden met overgangskansen

rig deeltje alleen afhangt van de huidige toestand en niet van de vorige toestanden; een deeltje springt dus volgens een Markovketen. We zijn geïnteresseerd in de kans dat een vatbaar/geïnfecteerd deeltje na n stappen vatbaar/geïnfecteerd is. Noem X i de toestand van een deeltje op tijdstip i, dus X i = S of X i = I . We berekenen de kans dat een vatbaar deeltje op tijdstip t = 0 op tijdstip t = n weer/nog steeds vatbaar is: P(X n = S|X 0 = S).

(47)

De 1-staps overgangsmatrix is:

1− p p q 1−q

(48)

waarin Pi j = P(X k+1 = j|X k = i) met i het rijnummer en j het kolomnummer. Definieer de n-staps overgangskans met Pi(n) j = P(X k+n = j|X k = i) en de n-staps overgangsma22 Griepepidemie Project Modelleren C


trix met

P (n) ,

dan geldt

P (n)

=

Pn

P(X n = S|X 0 = S) = P(n) 11

n 1− p p = . q 1−q (n−1) P11 · (1 − p) + P(n) 12 · q

= (n)

(n)

P11 +P12 =1

=

=

(n−1) P11 · (1 − p − q) + q (n−2) P11 · (1 − p − q) + q · (1 − p − q) + q (n−2) q + q(1 − p − q) + (1 − p − q)2 · P11

= . . .

q·

=

n−1 X

n (1 − p − q)k + P(0) 11 · (1 − p − q)

k=0 voor p+q6 =0

=

=

p+q 1 − (1 − p − q)n + · (1 − p − q)n q· 1 − (1 − p − q) p+q q p + · (1 − p − q)n p+q p+q

(1) Merk op dat inderdaad P(0) 11 = 1 en P11 = 1 − p. De kans dat een deeltje na n stappen geïnfecteerd is als het nu vatbaar is, P(n) 12 , kan eenvoudig gevonden worden uit: (n) (n) P11 + P12 = 1. De andere twee overgangskansen kunnen gevonden worden door in (n) de formules van P(n) 11 en P12 de 1-staps overgangskansen p en q te verwisselen. Het getal 1 − p −q bepaalt de convergentiesnelheid van de rij overgangskansen: hoe kleiner 1 − p − q des te groter de convergentiesnelheid. We krijgen de volgende n-staps overgangsmatrix:   q p p p + (1 − p − q)n − (1 − p − q)n  p+q p+q p+q p+q (49)   q q q p n n − (1 − p − q) + (1 − p − q) p+q p+q p+q p+q

De evenwichtssituatie/stabiele toestand van P(n) 11 is: p q q π11 := lim P(n) = lim + · (1 − p − q)n = . (50) 11 n→∞ n→∞ p + q p+q p+q q Ga na dat ook geldt: π21 = lim P(n) = . Het maakt dus niet uit voor de limiet wat 21 n→∞ p+q de beginsituatie is. Definieer π1 := lim P(X n = S) en π2 := lim P(X n = Z ). Deze n→∞

n→∞

limieten kunnen gevonden worden door het volgende stelsel op te lossen:  π1 = π1 (1 − p) + π2 · q  π2 = π2 · p + π2 (1 − q)  π1 + π2 = 1 23 Griepepidemie Project Modelleren C

(51)


We vinden als oplossing: π1 =

q p+q

en π2 =

p p+q .

We bekijken de verspreiding van een griepvirus binnen een (constante) groep deeltjes. Hiervoor definiëren we: Hierin is: N = totaal aantal deeltjes Sn = het aantal vatbare deeltjes op tijdstip n In = het aantal geïnfecteerde deeltjes op tijdstip n. We gebruiken nu dat bij benadering Sn = P(X n = S) · N en In = P(X n = I ) · N . Hierbij is X i de toestand van een willekeurig deeltje op tijdstip i. Voor grote waarden van N zal deze benadering goed kloppen, maar in werkelijkheid zullen de gevonden waarden natuurlijk afgerond moeten worden op gehele getallen. Nu is het bijbehorende stelsel recurrente betrekkingen  Sn+1 = Sn (1 − p) + In · q ∀n ∈ N      In+1 = Sn · p + In (1 − q) ∀n ∈ N (52)      Sn + In = N ∀n ∈ N Het stelsel (52) kunnen we ook als volgt schrijven:  1− p q S  Sn+1 = · n ∀n ∈ N In+1 p 1−q In  Sn + In = N ∀n ∈ N

(53)

Het stelsel recurrente betrekkingen kan opgelost worden m.b.v. eigenwaarden en ei1− p p . De karakteristieke vergelijking is: genvectoren. We definiëren A := q 1−q 1− p−λ p = λ2 + λ( p + q − 2) + (1 − p − q) = 0. det(A − λI ) = q 1−q −λ Er zijn dus twee eigenwaarden, λ1 = 1 en λ2 = 1 − p − q. De eigenruimte bij λ1 is E 1 = h(q, p)i en bij λ2 E 1− p−q = h(1, −1)i. We vinden dan de volgende oplossing: Sn q 1 n n = c1 · λ1 · + c2 · λ2 · ∀n ∈ N, c1 , c2 ∈ R, (54) In p −1 waarbij c1 en c2 bepaald kunnen worden uit de beginconditie. S0 = S0 geeft S0 = c1 q +c2 S0 p − I0 q S0 + I0 en c2 = . Dit en I0 = I0 geeft I0 = c1 p − c2 . Hieruit vinden we c1 = p+q p+q levert dan de volgende oplossing:  S0 q + I0 q S0 p − I0 p  Sn = + · (1 − p − q)n ∀n ∈ N    p+q p+q (55)   S0 p + I0 p I0 q − S0 q   In = + · (1 − p − q)n ∀n ∈ N p+q p+q 24 Griepepidemie Project Modelleren C


n n n n Merk op dat Sn = S0 P11 + I0 P21 en In = S0 P12 + I0 P22 . Inderdaad klopt het dat Sn en In op n = 0 gelijk zijn aan respectievelijk S0 en I0 . Merk op dat Sn en In alleen convergeren als 1 − p − q = det A ≤ 1. Ze convergeren dan naar een ’gewogen’ gemiddelde van S0 en I0 . Als we in het stelsel (52) links en rechts de limiet nemen voor n → ∞, waarbij we definiëren: S∞ := lim Sn en I∞ := lim In , vinden we n→∞

n→∞

S∞ · p = I∞ · q. We constateren dat deze limieten overeenkomen met de limieten van stelsel 55. In figuur 13 zijn Sn en In weergegeven voor verschillende (continue) waarden van n. Hierbij is p = 0.15, q = 0.1, S0 = 100 en I0 = 0. De limieten zijn nu dus: S∞ = 40 en I∞ = 60.

Figuur 13: Sn en In weergegeven voor verschillende waarden van n

6.2

SIR model

We bekijken weer het SIR model en hebben dus 3 verschillende toestanden waarin de deeltjes zich kunnen bevinden: vatbaar (S), geïnfecteerd (I) en hersteld (R). Een geïnfecteerd deeltje en een vatbaar deeltje zullen resulteren in twee geïnfecteerde deeltjes met intensiteit α (dus een geïnfecteerd deeltje zal een vatbaar deeltje besmetten met intensiteit α). En een geïnfecteerd deeltje zal herstellen met intensiteit β. Er geldt dus P(infectie in (t, t + dt)) = αdt. De parameters α en β hangen dus niet van de tijd af. We nemen weer aan dat het totaal aantal deeltjes constant is, zeg N . We hebben dus de volgende twee overgangen: α

S + I −→ 2I 25 Griepepidemie Project Modelleren C

(56)

Technische Universiteit Eindhoven University of Technology β

I −→ R.

(57)

In figuur 14 is een schematisch model weergegeven.

Figuur 14: Schematische weergave van het SIR-model Zij X een stochast, X = (X 1 , X 2 ) ∈ N2 , laat X 1 , het aantal vatbare deeltjes en X 2 , het aantal geïnfecteerde deeltjes zijn. De verandering van de kans (tegen de tijd t) dat er S d vatbare deeltjes en I geïnfecteerde deeltjes zijn, P(X (t + dt) = (S, I )) kan uitgedrukt dt worden in 3 kansen: 1. De kans dat er S + 1 vatbare deeltjes en I − 1 geïnfecteerde deeltjes zijn op tijdstip t, P(X (t) = (S + 1, I − 1)) 2. De kans dat er S vatbare deeltjes en I + 1 geïnfecteerde deeltjes zijn op tijdstip t, P(X (t) = (S, I + 1)) 3. De kans dat er S vatbare deeltjes en I geïnfecteerde deeltjes zijn op tijdstip t, P(X (t) = (S, I )). We vinden de volgende relatie: P(X (t + dt) = (S, I )) = P(X (t) = (S + 1, I − 1))(α(S + 1)(I − 1)dt + o(dt)) + P(X (t) = (S, I + 1))(β(I + 1)dt + o(dt))

(58)

+ P(X (t) = (S, I ))(1 − ((αS I + β I )dt + o(dt))) ofwel P(X (t + dt) = (S, I )) − P(X (t) = (S, I )) = P(X (t) = (S + 1, I − 1))(α(S + 1)(I − 1)dt + o(dt)) + P(X (t) = (S, I + 1))(β(I + 1)dt + o(dt)) − P(X (t) = (S, I ))((αS I + β I )dt + o(dt)). (59) 26 Griepepidemie Project Modelleren C


Links en rechts delen door dt levert P(X (t + dt) = (S, I )) − P(X (t) = (S, I )) = P(X (t) = (S + 1, I − 1))(α(S + 1)(I − 1) + dt + P(X (t) = (S, I + 1))(β(I + 1) + − P(X (t) = (S, I ))((αS I + β I ) +

o(dt) dt )

o(dt) dt ).

(60)

Als we nu zowel links als rechts de limiet nemen voor dt → 0 krijgen we d P(X (t) = (S, I )) = P(X (t) = (S + 1, I − 1))α(S + 1)(I − 1) dt + P(X (t) = (S, I + 1))β(I + 1)

(61)

− P(X (t) = (S, I ))(αS I + β I ). Er ontstaat dus een niet-lineaire differentiaalvergelijking. Een bestaat geen exacte algebraïsche oplossing. Toch kunnen we wel iets zeggen over het model. We zullen twee dingen onderzoeken: de kans dat een bepaald aantal deeltjes niet geïnfecteerd raakt en de verwachte tijd van de duur van de epidemie. In figuur 15 (zie volgende pagina) zijn de overgangsintensiteiten schematisch weergegeven voor N = 4. Het aantal vatbare deeltjes wordt weergegeven op de horizontale as, het aantal geïnfecteerde deeltjes op de verticale as. Het aantal resistente deeltjes wordt niet weergegeven. Als S en I beide groter zijn dan 0 zijn er meerdere overgangen mogelijk. De eerste mogelijkheid is dat een geïnfecteerd deeltje een vatbaar deeltje besmet. In de figuur gaan we dan dus naar linksboven. Merk op dat dit niet mogelijk is als we al bovenin de figuur zitten, alle deeltjes zijn dan geïnfecteerd. De andere mogelijkheid is dat een geïnfecteerd deeltje herstelt, zodat dit deeltje uit het model verdwijnt. In de figuur gaan we dan één stap naar onder. Definieer Y S ∗ (S, I ) als de gebeurtenis dat bij start in het punt (S, I ) de eindtoestand zal zijn (S ∗ , 0). Uit Figuur 15 volgt Y0 (0, I ) = 1. Ook geldt natuurlijk Y S (S, 0) = 1. Als we in het punt (S, I ) zijn, kunnen we nooit eindigen in het punt (S ∗ , 0) voor S ∗ > S. We kunnen immers nooit meer vatbare deeltjes krijgen. Als we in het punt (S, I ) zijn, βI I gaan we met kans αSαS I +β I naar het punt (S − 1, I + 1) en met kans αS I +β I naar het punt


o(dt) dt )


Figuur 15: Schematische weergave van overgangsintensiteiten met N = 4. De overgangsintensiteit van (S, I ) naar (S −1, I +1) is αS I en de overgangsintensiteit van (S, I ) naar (S, I − 1) is β I .

(S, I − 1) (volgt ook uit figuur 15). We vinden nu  β αS P(Y S ∗ (S, I )) = αS+β P(Y S ∗ (S − 1, I + 1)) + αS+β P(Y S ∗ (S, I − 1)) voor S, I ≥ 1          P(Y0 (0, I )) =1 voor I ≥ 1         P(Y S ∗ (S ∗ , 0)) = 1    P(Y S ∗ (S, I ))         P(Y S ∗ (S, 0))       S, I ≥ 0

=0

voor S ∗ > S

=0

voor S ∗ 6= S

,S+ I ≤ N (62)

Merk ten eerste op dat geldt: P(Y S ∗ (S , I )) = ∗

β αS + β

I

,

(63)

om vanuit het punt (S ∗ , I ) het punt (S ∗ , 0) te bereiken moeten we immers I keer naar beneden gaan. Des te groter β des te groter deze kans; dit wisten we van tevoren al. Nu zijn we geïnteresseerd in P(Y S ∗ (S ∗ + 1, I )). Om vanuit het punt (S ∗ + 1, I ) in het punt (S ∗ , 0) te komen, zijn meerdere routes mogelijk. We bekijken deze kans eerst voor 28 Griepepidemie Project Modelleren C


I = 0, 1, 2. Voor I = 0 zijn en blijven we in het verkeerde punt, dus P(Y S ∗ (S ∗ + 1, 0)) = 0.

(64)

Voor I = 1 gebruiken we het gevonden stelsel vergelijkingen en vinden α(S ∗ + 1) P(Y S ∗ (S + 1, 1)) = α(S ∗ + 1) + β ∗

β αS ∗ + β

2

β αS ∗ + β

3

+

β P(Y S ∗ (S ∗ + 1, 0), (65) α(S ∗ + 1) + β

+

β P(Y S ∗ (S ∗ + 1, 1). (66) ∗ α(S + 1) + β

en voor I = 2 α(S ∗ + 1) ∗ P(Y S (S + 1, 2)) = α(S ∗ + 1) + β ∗

Hieruit vinden we de volgende recurrente betrekking  I +1 ∗ + 1) β β α(S    P(Y S ∗ (S ∗ + 1, I )) = + P(Y S ∗ (S ∗ + 1, I − 1)  ∗ ∗ ∗  α(S + 1) + β αS + β α(S + 1) + β         

voor I ≥ 0 P(Y S ∗ (S ∗ + 1, 0)) = 0 . (67)

We hebben dus een recurente betrekking van de vorm   a I = c1 c2I +1 + c3 a I −1 voor I ≥ 0 

(68)

a0 = 0.

De oplossing hiervan is c1 c22 (c2I − c3I ) , (69) c2 − c3 dit is na te gaan door de oplossing gewoon in te vullen in de recurrente betrekking. Na vereenvoudigen vinden we voor P(Y S ∗ (S ∗ + 1, I )) de volgende uitdrukking I I ! β β − (1 + S ∗ )β αS ∗ + β α(S ∗ + 1) + β P(Y S ∗ (S ∗ + 1, I )) = . (70) αS ∗ + β aI =

Op analoge wijze kan eventueel P(Y S ∗ (S ∗ + 2, I )) gevonden worden. In figuur 16 is voor α = 0.5, β = 5, S = 14, I = 7 en variabele S ∗ de kans P(Y S ∗ (S, I )) weergegeven. We zien in deze figuur alleen de kans voor S ∗ < 14, immers voor S ∗ = 14 is deze kans heel erg klein en voor S ∗ > S = 14 is de kans gelijk aan 0.


Technische Universiteit Eindhoven University of Technology P@YS* H14,7LD 0.20

0.15

0.10

0.05

0

1

2

3

4

5

6

7

8

9

10

11

12

13

S*

Figuur 16: P(Y S ∗ (S, I )) voor verschillende S ∗

Wat is de verwachte tijd totdat geen enkel deeltje meer geïnfecteerd is? Hiervoor definiëren we Z (S, I ) als de tijd die het duurt totdat geen enkel deeltje meer geïnfecteerd is als er S vatbare en I geïnfecteerde deeltjes zijn (Z (S, I ) is een stochast). Dan is E[Z (S, I )] de verwachte tijd tot geen enkel deeltje meer geïnfecteerd is, als er S vatbare deeltjes en I geïnfecteerde deeltjes zijn. Als we beginnen met enkel vatbare en resistente deeltjes geldt natuurlijk E[Z (S, I )] = E[Z (S, 0)] = 0. De kans dat vanuit het punt αS en de kans dat een stap vanuit (S, I ) een stap naar linksboven wordt gedaan is αS+β

(S, I ) naar beneden wordt gedaan is

β αS+β .

Merk op dat de verwachte tijd om in het

is. Om E[Z (S, I )] te berekenen, kunnen we de kansen van punt (S, I ) te blijven de overgangen combineren met de verwachte tijden tot elke mogelijke overgang. We vinden nu de volgende vergelijkingen 1 αS I +β I

 1 αS I  E[Z (S, I )] = + E[Z (S − 1, I + 1)]    αS I + β I αS I + β I    βI   + E[Z (S, I − 1)] voor S, I ≥ 1   αS I + β I      βI 1 E[Z (0, I )] = + E[Z (0, I − 1)] voor I ≥ 1   βI βI         E[Z (S, 0)] = 0       S, I ≥ 0, S+I ≤N


(71)


ofwel

 1 αS   + E[Z (S − 1, I + 1)] E[Z (S, I )] =   αS I + β I αS + β    β   E[Z (S, I − 1)] voor S, I ≥ 1 +    αS + β     1 E[Z (0, I )] = + E[Z (0, I − 1)] voor I ≥ 1   βI         E[Z (S, 0)] = 0       S, I ≥ 0, S+I ≤N

(72)

Als we bijvoorbeeld beginnen met één geïnfecteerd en één vatbaar deeltje is de kans β α α+β dat het geïnfecteerde deeltje het vatbare deeltje besmet. Met kans α+β zal het geïnfecteerde deeltje genezen voordat het vatbare deeltje is besmet. We krijgen dus 1 α β + E[Z (0, 2)] + E[Z (1, 0)] α+β α+β α+β α β 3α + 2β(1 + β) 1 1 1 1 + + + = . = α + β α + β β 2β α+β β 2β(α + β)

E[Z (1, 1)] =

(73)

Als er aan het begin slechts geïnfecteerde en resistente deeltjes aanwezig zijn is E[Z (S, I )] = E[Z (0, I )] eenvoudig te bepalen. Er is iedere keer namelijk maar één overgang mogelijk. We krijgen dan dus I I X 1 1X1 E[Z (0, I )] = = . (74) iβ β i i=1

i=1

Voor S = 1 geldt E[Z (1, I )] =

α β 1 + E[Z (0, I + 1)] + E[Z (1, I − 1)] α+β α+β α+β

I +1 1 α 1X1 1 α = + + + E[Z (0, I )] + E[Z (1, I − 2)] α+β α+ββ i α(I − 1) + β(I − 1) α + β i=1

= ... =

I +1 k 1 α 1 XX 1 1 + ... + + + E[Z (1, 0)] αI + βI α·1+β ·1 α+ββ i k=1 i=1

I +1 I +1 k 1 X1 α 1 XX 1 = + α+β i α+ββ i i=1

k=1 i=1

(75)



Een eenvoudige gesloten formule voor deze recurrente betrekking is niet makkelijk te vinden. Wel kunnen we iets zeggen over het gedrag van E[Z (1, I )] voor I groot. Er I X 1 ≈ log I . Merk ook op dat geldt namelijk i 1

I +1 X k X 1

I +1 X I +1 X 1

I +1 I +1 X X I +2 1 − (I + 1). = (I + 2 − i) = − 1 = (I + 2) = i i i i i k=i k=1 i=1 i=1 i=1 i=1 i=1 (76) We vinden dus dat voor I groot het volgende geldt E[Z (1, I )] ≈

I +1 X 1

1 α 1 log(I + 1) + · (I + 1) log(I + 1) − (I + 1). α+β α+ββ

(77)

Echter kunnen we voor specifieke α, β, S en I wel E[Z (S, I )] berekenen met behulp van een computersimulatie. We hebben een recursief computerprogramma geschreven dat E[Z (S, I )] berekent. In de figuren 17 en 18 zijn E[Z (1, 1)] en E[Z (6, 2)] grafisch weergegeven voor verschillende waarden van α en β. Α 10 5 0

1.5

1.0

E@ZH1,1LD 0.5

0.0 0

5

Β 10

Figuur 17: Schematische weergave van E[1, 1] voor verschillende waarden van α en β



0

Β 5

10

3

2

E@ZH6,2LD 1 0

5 0 10

Α

Figuur 18: Schematische weergave van E[6, 2] voor verschillende waarden van α en β

6.3

Het Reed frost model

In dit model nemen we aan dat de tijd een discrete variabele is. Weer bekijken we het aantal vatbare(S) en geïnfecteerde (I) deeltjes. We gaan ervan uit dat geïnfecteerde deeltjes maar voor één tijdseenheid geïnfecteerd zijn en daarna direct immuun worden. Voor n = 0, 1, 2, ..., laat Sn , In respectievelijk het aantal vatbare en het aantal geïnfecteerde deeltjes op tijdstip n zijn (Sn , In ∈ N). We nemen aan dat elk vatbaar deeltje een even grote kans, p, heeft om op een bepaald tijdstip geïnfecteerd te raken door een individueel geïnfecteerd deeltje. Als we aannemen dat infecties onafhankelijk van elkaar plaatsvinden, dan heeft elk vatbaar deeltje kans (1 − p) In om niet geïnfecteerd te raken op tijdstip n. Als een deeltje geïnfecteerd is geraakt op tijdstip n, kan het andere deeltjes gaan infecteren maar pas vanaf tijdstip n + 1. Merk op dat de kans dat een vatbaar deeltje geïnfecteerd raakt op tijdstip n gelijk is aan 1 − (1 − p) In . De kans dat er In+1 deeltjes geïnfecteerd raken is dan (1 − (1 − p) In ) In+1 . Als er In+1 deeltjes geïnfecteerd raken, zijn er Sn+1 deeltjes die niet geïnfecteerd raken. De kans dat er Sn+1 deeltjes niet geïnfecteerd raken is ((1− p) In ) Sn+1 = (1− p) In Sn+1 . Merk op dat het aantal geïnfecteerde deeltjes op tijdstip n + 1 gelijk is aan het aantal vatbare deeltjes op het vorige tijdstip, dus tijdstip n, min het aantal vatbare deeltjes op tijdstip n + 1: In+1 = Sn − Sn+1 . Na



itereren vinden we: S0 + I0 = S1 + I1 + I0 = S2 + I2 + I1 + I0 = ... = Sn +

n X

Ij.

(78)

j=0

De kans op In+1 geïnfecteerde deeltjes op tijdstip n wordt nu gegeven door Sn P(In+1 |Sn , In ) = (1 − (1 − p) In ) In+1 (1 − p) In Sn+1 Sn+1

(79)

Sn waarbij het aantal mogelijkheden is om Sn+1 vatbare deeltjes uit Sn te trekken. Sn+1 Dit is een binomiale verdeling met parameters Sn en 1 − (1 − p) In . Met deze gegevens schrijven we een simulatie in Java. De simulatie houdt het aantal geïnfecteerde deeltjes en het aantal vatbare deeltjes bij. In figuur 19 zien we voor een aantal simulaties het aantal vatbare deeltjes uitgezet tegen de tijd. Sn 5000

4000

3000

2000

1000

0 0

5

10

15

20

n

Figuur 19: Het aantal vatbare deeltjes tegen de tijd uitgezet. Hierbij is steeds de kans p = 5 · 10−4 , als begin aantallen zijn gekozen: S0 = 5000 en I0 = 1.



In figuur 20 zien we bij de getoonde aantallen vatbare deeltjes van figuur 19 het aantal geïnfecteerde deeltjes. We merken op dat het aantal vatbare deeltjes inderdaad In 1200

1000

800

600

400

200

0 0

5

10

15

20

n

Figuur 20: Het aantal geïnfecteerde deeltjes tegen de tijd uitgezet. Hierbij is steeds de kans p = 5 · 10−4 , als begin aantallen zijn gekozen: S0 = 5000 en I0 = 1. afneemt naarmate de verstreken tijd toeneemt. De simulatie stopt als het aantal geïnfecteerde deeltjes 0 is geworden (of als de maximale simulatietijd verstreken is, maar dat is hier niet het geval). Immers als In = 0 veranderen Sn en In niet meer. Het feit dat sommige grafieken eerder stoppen (geen waarden meer aannemen voor n groter dan een bepaalde waarde) komt dus doordat het aantal geïnfecteerde deeltje voor die n gelijk aan 0 is. Bij één van de simulaties werd in de eerste stap I1 = 0, echter zijn de grafieken van deze niet goed zichtbaar in de figuren.



7

Conclusie

Van een griepepidemie kunnen vele verschillende modellen gemaakt worden, zowel probabilistisch als deterministisch (discreet en continu). In dit verslag hebben we een aantal van deze modellen geanalyseerd. Deze modellen zijn erop gebaseerd dat deeltjes op elk tijdstip in een bepaalde toestand zijn. Van een aantal modellen kan uitgerekend worden hoeveel deeltjes in elke toestand op elk tijdstip zijn, dit kan met behulp van eigenwaarden en eigenvectoren. Bij het modelleren van een griepepidemie speelt het reproductiegetal een belangrijke rol. Dit getal geeft het verwachte aantal besmettingen van één geïnfecteerd deeltje weer. Het reproductiegetal is een zeer belangrijk kengetal binnen de epidemiologie. Hoe groter het reproductiegetal is des te besmettelijker is de ziekte. Dit getal kan berekend worden uit de Jacobi-matrix van de differentiaalvergelijkingen. In het eindverslag zullen we deze en nieuwe modellen verder onderzoeken. Ook gaan we dan het reproductiegetal nader analyseren en kijken naar verschillende strategieën om te vaccineren.



Bibliografie [1] SIR Models of Epidemics, F. Débarre, 2010, pp. 1-2. [2] Notes On R0 , Holland James Jones, Tech. report, Department of Anthropological Sciences, Stanford University, 2007, pp. 1-5. [3] SIR Model of Epidemics, T. Tassier, Epidemics and Development Policy, Fordham University NY, 2005, pp. 2-4. [4] Basic Reproduction Number for Infectious Dynamics Models and the Global Stability of Stationary Points, J. Guardiola, A. Vecchio, Wseas Transactions on Biology and Biomedicine, Napoli, Italy, 2003, p. 2. [5] Using Probabilistic Models to Infer Infection Rates in Viral Outbreaks, Auteur onbekend, 2007, pp. 5-7.


Griepepidemie Project Modelleren C

Recommend Documents