M5858 Spojité deterministické modely I, cvičenie, 8.12.2021
Epidemiologické modely
Pri štúdiu šírenia epidémie môžeme vo všeobecnosti rozdeliť populáciu na štyri skupiny. Tieto
skupiny, respektíve počty jedincov v jednotlivých skupinách budeme značiť postupne symbolmi
S, E, I a R. Skupina S pozostáva zo zdravých jedincov náchylných k infekcii. Skupina E je
tvorená jedincami nachádzajúcimi sa v latentnom období teda tými, ktorí sú už síce nakazení,
ale ešte nie sú schopní nákazu šíriť. Skupinu označenú symbolom I tvoria nakazení jedinci, ktorí
sú zároveň šíriteľmi ochorenia. V skupine R sú zahrnutí jedinci, ktorí ochoreli, ale ochorenie
ďalej nešíria, teda napríklad vyliečení alebo izolovaní jedinci. Na základe konkrétnych predpokladov
o chovaní epidémie a populácie môžeme charakterizovať spôsob prechodu jedincov
medzi jednotlivými skupinami a tým konštruovať konkrétne modely. V každom z nasledujúcich
modelov budeme predpokladať, že veľkosť populácie sa v čase nemení. Zdôraznime, že všetky
uvažované systémy sú autonómne s polynomiálnou pravou stranou definovanou na celom Rn
,
čo znamená, že akýkoľvek počiatočný problém má jediné úplné riešenie.
Model SI
Začneme najjednoduchším modelom. Predpokladajme, že dĺžka latentného obdobia je zanedbateľná
a ignorujme vplyv izolácie či uzdravenia, teda v našom modeli nebudeme uvažovať
populačné skupiny E a R. Tieto predpoklady sú relevantné v počiatočných štádiách niektorých
ochorení. Vzhľadom k tejto skutočnosti je konštantnosť veľkosti populácie určená vzťahom
S + I = N. Tento model môžeme reprezentovať nasledujúcou schémou:
S
β
// I
kde β > 0 je tzv. koeficient šírenia nákazy, ktorý vyjadruje mieru prechodu jedincov zo skupiny
S do skupiny I. Môžeme si ho predstaviť ako súčin relatívneho počtu kontaktov v populácii
za jednotku času a pravdepodobnosti prenosu choroby pri kontakte z člena skupiny I na člena
skupiny S. Prírastok do skupiny I respektíve úbytok zo skupiny S za čas ∆t je
∆I = −∆S = ∆tβSI.
Po vydelení rovnice veličinou ∆t a limitným prechodom ∆t → 0 dostaneme nasledujúci systém
diferenciálnych rovníc:
S = −βSI,
I = βSI.
Všimnite si, že rovnosť S + I = 0 korešponduje s predpokladom S(t) + I(t) = N. Vzťah
S = N − I môžeme dosadiť do druhej rovnice, čo nás privedie k rovnici so separovanými
premennými:
I = β(N − I)I.
Nájdime riešenie tejto rovnice spĺňajúce počiatočnú podmienku I(0) = I0. Ak I0 = 0 alebo
I0 = N, potom má riešenie zrejme tvar I(t) ≡ 0 alebo I(t) ≡ N. V opačnom prípade postupujme
nasledujúcim spôsobom:
I = β(N − I)I,
I
(N − I)I
= β
t
0
ds ,
1
t
0
I (s)
(N − I(s))I(s)
ds =
t
0
β ds = βt.
Integrál na ľavej strane vypočítame nasledujúcim spôsobom:
t
0
I (s)
(N − I(s))I(s)
ds
I(s) = u t I(t)
I (s) ds = du 0 I0
=
I(t)
I0
1
(N − u)u
=
1
N
I(t)
I0
1
N − u
+
1
u
du
=
1
N
log
u
N − u
I(t)
I0
=
1
N
log
I(t)
N − I(t)
− log
I0
N − I0
=
1
N
log
I(t)
I0
·
N − I0
N − I(t)
.
Výsledok dosaďme a upravme:
1
N
log
I(t)
I0
·
N − I0
N − I(t)
= βt,
I(t)
I0
·
N − I0
N − I(t)
= eNβt
,
I(t) · (N − I0) = (N − I(t))I0eNβt
,
I(t) · N − I0 + I0eNβt
= NI0eNβt
,
I(t) =
NI0eNβt
N − I0 + I0eNβt
=
N
1 +
N
I0
− 1 e−βNt
.
Grafom tejto funkcie je logistická krivka. Všimnime si, že limt→∞ I(t) = N. Funkcia počtu
nových prípadov je v skutočnosti mierou zmeny funkcie I(t), teda jej deriváciou:
I (t) =
βN2 N
I0
− 1 e−βNt
1 +
N
I0
− 1 e−βNt
2 =
βN2 N
I0
− 1 eβNt
N
I0
− 1 + eβNt
2 .
Jej graf sa nazýva epidemická krivka. Typickým javom je počiatočné narastanie počtu nových
prípadov. Po dosiahnutí maxima počty nových prípadov neustále klesajú. Presvedčme sa o tom
výpočtom. Nájdeme stacionárne body prvej derivácie: I (t) =
β2
N3 N
I0
− 1 eβNt
eβNt
+
N
I0
− 1
2
− 2β2
N3 N
I0
− 1 eβNt
eβNt
+
N
I0
− 1 eβNt
N
I0
− 1 + eβNt
4
=
β2
N3 N
I0
− 1 eβNt
eβNt
+
N
I0
− 1 − 2β2
N3 N
I0
− 1 e2βNt
N
I0
− 1 + eβNt
3
=
β2
N3 N
I0
− 1 eβNt N
I0
− 1 − eβNt
N
I0
− 1 + eβNt
3 .
2
Stacionárny bod je riešením rovnice
N
I0
− 1 − eβNt
= 0 ⇔ eβNt
=
N
I0
− 1 ⇔ t =
1
βN
log
N
I0
− 1 .
Vzhľadom k tomu, že druhá derivácia je naľavo od tohto bodu kladná a napravo záporná, má
funkcia I (t) v tomto bode globálne maximum a naľavo od tohto bodu je rastúca a napravo je
klesajúca. Funkcia I(t) má v tomto bode inflexný bod a funkčnú hodnotu N
2
.
Model SIR
Teraz do našich úvah pridáme populačnú skupinu R, čo nás privedie k modelu SIR. Znázorniť
ho môžeme nasledujúcou schémou:
S
β
// I
ν // R
kde ν > 0 je koeficient popisujúci rýchlosť prechodu členov skupiny I do skupiny R. Rovnice
definujúce tento model majú nasledujúci tvar:
S = −βSI,
I = βSI − νI,
R = νI.
Fakt, že S +I +R = 0, opäť korešponduje s predpokladom S(t)+I(t)+R(t) = N o konštantnej
veľkosti populácie. Venujme sa nasledujúcemu počiatočnému problému:
S(0) = S0 > 0, I(0) = I0 > 0, R(0) = 0, S0 + I0 = N.
Všimnime si, že prvé dve rovnice systému nezávisia od premennej R. Môžeme ich preto študovať
samostatne ako autonómny systém v rovine:
S = −βSI,
I = βSI − νI.
Začneme nulklinami. S-nulklina je zjednotením priamok S = 0 a I = 0. Na druhej strane,
I-nulklina je zjednotením priamok S = ν/β a I = 0. To znamená, že každý jeden bod na osi S
je stacionárny bod a žiadne iné stacionárne body neexistujú. Ďalej si všimnime, že kladná časť
3
osi I je invariantná množina (v skutočnosti je trajektóriou riešenia S(t) = 0, I(t) = exp(−νt)).
V konečnom dôsledku môžeme povedať, že prvý kvadrant je pozitívne invariantná množina, čo
by sme v kontexte modelovaného javu aj očakávali. Dokonca pre akékoľvek riešenie spĺňajúce
S(0) > 0 a I(0) > 0 platí S(t) > 0 a I(t) > 0 pre každé t ∈ [0, ∞). Pre S ≤ ν/β je I ≤ 0,
čo znamená, že počiatočné podmienky v tejto oblasti vedú k riešeniam, ktorých zložka I(t)
pre t ∈ [0, ∞) klesá. Takže v takom prípade k prepuknutiu epidémie nedôjde. Ďalej budeme
predpokladať, že S0 > ν/β.
Pokúsme sa odvodiť tvar trajektórií. Za tým účelom budeme študovať vlastnosti riešení
rovnice
dI
dS
=
βSI − νI
−βSI
=
ρ − S
S
,
kde ρ = ν/β. Pre 0 < S < ρ riešenia tejto rovnice rastú a pre S > ρ klesajú. Navyše pre ich
druhú deriváciu platí
d2
I
dS2
= −
ρ
S2
< 0,
čo znamená, že tieto riešenia sú konvexné funkcie. Rovnicu s počiatočnou podmienkou I(S0) =
I0 môžeme dokonca vyriešiť a okrem iného nájsť hodnotu globálneho maxima:
I (S) =
ρ
S
− 1
S
S0
ds
I(s) = u S I(S)
I (s) ds = du S0 I0
S
S0
I (s) ds =
S
S0
ρ
s
− 1 ds
I(S)
I0
du = [ρ log s − s]S
S0
I(S) = ρ log S − S − ρ log S0 + S0 + I0
= ρ log
S
S0
− S + N.
V bode globálneho maxima platí
I(ρ) = ρ log
S
S0
− 1 + N.
Vráťme sa k pôvodnému trojrozmernému systému. Už sme sa presvedčili o tom, že riešenie
(S(t), I(t), R(t)) nami študovaného počiatočného problému spĺňa I(t) > 0 pre každé t ∈ [0, ∞).
4
Podľa poslednej rovnice je R (t) > 0 pre každé t ∈ [0, ∞), čo implikuje R(t) > 0 pre každé
t ∈ (0, ∞). Vďaka tomu pre každé nezáporné t platí 0 ≤ S(t), I(t), R(t) ≤ N, takže riešenie je
ohraničené, a preto je definované pre každé t ≥ 0. Teraz sa pokúsime vyšetriť správanie riešenia
pre t → ∞. Medzitým ešte raz formálne odvodíme niektoré zrejmé skutočnosti. Sčítanie νnásobku
prvej rovnice s βS-násobkom treťej rovnice nás privedie k rovnici
νS + βSR = 0,
ktorá má podobnú štruktúru ako exaktná diferenciálna rovnica. Vynásobme túto rovnicu funkciou
k(S, R), ktorú špecifikujeme neskôr:
k(S, R)νS + k(S, R)βSR = 0.
Ak existuje funkcia F(S, I, R) s vlastnosťou
∂F
∂S
= k(S, R)ν,
∂F
∂I
= 0,
∂F
∂R
= k(S, R)βS,
potom je táto funkcia zrejme prvý integrál nášho systému. Na to, aby táto funkcia existovala,
nám stačí
∂
∂R
(k(S, R)ν) =
∂
∂S
(k(S, R)βS).
Takže budeme postupovať rovnako ako v prípade exaktných rovníc. Označme M(S, R) = ν
a N(S, R) = βS. Keďže MR = 0 = β = NS, funkcia k(S, R) je netriviálna. Vzhľadom k
tomu, že (NS − MR)/M = β/ν = 1/ρ je funkcia premennej R (konštantná), môžeme hľadať
integračný faktor ako funkciu k(R), ktorá je riešením rovnice k = k/ρ. Jej partikulárne riešenie
je k = exp(R/ρ). Pozrime sa na funkciu F:
F(S, R) = νe
R
ρ dS = νSe
R
ρ + c(R),
FR = βSe
R
ρ + c (R) = βSe
R
ρ ⇒ c (R) = 0 ⇒ c(R) = c0 ∈ R.
Napríklad funkcia F(S, R) = νSe
R
ρ je prvým integrálom nášho systému, čo znamená, že jej
hodnoty sú pozdĺž riešení konštantné. Takže ak F(S0, 0) = νS0, potom pre každé t platí
F(S(t), R(t)) = νS0, teda
S(t) = S0e
−
1
ρ
R(t)
.
Ukážme si jednoduchší spôsob, ktorý nás privedie k rovnakému výsledku. Ak vydelíme prvú
rovnicu treťou, dostaneme lineárnu diferenciálnu rovnicu
dS
dR
= −
S
ρ
,
ktorej riešenie má tvar
S = S0e
−
1
ρ
R
,
čo po dosadení času t dá ten istý výsledok.
Keďže R(t) ≥ 0, platí 0 < S ≤ S0 < N a navyše na základe analýzy S-nulklín vieme, že
S(t) je klesajúca funkcia. Z poslednej rovnice tiež môžeme vyjadriť funkciu R(t):
R(t) = −ρ log
S(t)
S0
.
5
Ďalej platí
R (t) = ν(N − S(t) − R(t)) ≤ ν(N − R(t)),
R (t) + νR(t) ≤ νN, eνt
R (t)eνt
+ νR(t)eνt
≤ νNeνt
,
(R(t)eνt
) ≤ νNeνt
.
Poslednú rovnosť môžeme integrovať v medziach od 0 do t, čo nás privedie k nerovnosti
R(t)eνt
≤ N(eνt
− 1), e−νt
R(t) ≤ N(1 − e−νt
) < N.
V konečnom dôsledku tiež máme
I(t) = N − S(t) − R(t) < N,
keďže S(t) > 0 a R(t) ≥ 0.
Funkcia R(t) je rastúca a ohraničená, a preto existuje konečná limita R∞ = limt→∞ R(t).
Uvedomme si, že do skupiny R sa jedinec môže dostať len tak, že v nejakom momente ochorie.
Jedinci sa zo skupiny R už dostať nemôžu, takže hodnota R∞ vyjadruje, do akej miery sa
infekcia rozšíri. Pokúsme sa túto hodnotu vyčísliť. Vzhľadom k tomu, že
R (t) = νI(t) = ν(N − S(t) − R(t)) = ν(N − S0e
−
1
ρ
R(t)
− R(t)),
platí
lim
t→∞
R (t) = ν(N − S0e
−
1
ρ
R∞
− R∞).
Odbočme na chvíľu, aby sme sa mohli presvedčiť o tom, že existencia konečných limít limt→∞ R(t)
a limt→∞ R (t) implikuje nulovosť druhej limity. Tvrdenie dokážeme sporom, takže budeme
predpokladať, že spomínaná limita je nenulová. Vzhľadom k tomu, že R (t) = νI(t) > 0, má
zmysel uvažovať len prípad limt→∞ R (t) = L > 0. Zvoľme 0 < ε < L ľubovoľne. Existuje t0 ≥ 0
také, že pre každé t ≥ t0 platí R (t) ≥ ε. Integrovaním tejto nerovnosti v medziach od t0 do t
dostaneme R(t) ≥ ε(t − t0) + R(t0). To znamená, že pre t → ∞ funkcia R(t) rastie nad všetky
medze, čo je v spore s predpokladom o konečnosti limity limt→∞ R(t). Vzhľadom k tretej rovnici
nášho systému platí
lim
t→∞
R (t) = lim
t→∞
I(t) = 0.
Limitná hodnota R∞ je preto koreňom rovnice
F(x) = N − S0e
−
1
ρ
x
− x = 0.
Ukážeme, že uvedená rovnica má jediný kladný koreň a nájdeme jeho dolný odhad. Platí
F (x) =
S0
ρ
e
−
1
ρ
x
− 1 = 0 ⇔ x = R∗
= −ρ log
ρ
S0
,
Takže funkcia F (x) má jediný nulový bod R∗
> 0. Naľavo od tohto bodu je F (x) kladná, takže
F(x) rastie a napravo je F (x) záporná, takže F(x) klesá. Na základe vzťahov limx→±∞ F(x) =
−∞, F(0) = N − S0 > 0 a F(R∗
) = N − ρ > 0 môžeme usúdiť, že funkcia F(x) má práve dva
korene, jeden záporný a jeden kladný, pričom platí R∗
< R∞. Navyše máme aj horný odhad:
R∞ = N − S0e
−
1
ρ
R∞
< N − S0e
−
1
ρ
N
,
6
čo je dôsledkom toho, že F(N) < 0, teda R∞ < N.
Zhrňme všetky získané poznatky. Riešenie (S(t), I(t), R(t)) počiatočného problému S(0) =
S0, I(0) = I0 a R(0) = 0, kde S0 + I0 = N spĺňa nasledujúce limitné vzťahy:
lim
t→∞
R(t) = R∞, lim
t→∞
I(t) = 0, lim
t→∞
S(t) = N − R∞ = S∞,
pričom platí
ρ log
S0
ρ
< R∞ < N − S0e
−
1
ρ
N
, S0e
−
N
ρ < S∞ < ρ,
kde posledná nerovnosť je dôsledkom rovnosti S∞ = S0 exp(−R∞/ρ) a prvej nerovnosti. Rozsah
epidémie môžeme zmenšiť zväčšením parametra ρ, teda buď zväčšením parametra ν alebo
zmenšením parametra β. V praxi to znamená rýchlejšiu izoláciu a uzdravovanie nakazených
jedincov alebo obmedzenie kontaktov a zvýšenie odolnosti voči infekcii.
Model SIRS
Tento model sa od toho predošlého líši tým, že jedinci zo skupiny R sa opäť môžu stať náchylní
k infekcii. Schematicky ho môžeme znázorniť nasledujúcim obrázkom:
S
β
// I
ν // R
γ
hh
kde γ > 0 je koeficient popisujúci rýchlosť prechodu členov skupiny R do skupiny S. Rovnice
definujúce tento model majú nasledujúci tvar:
S = −βSI + γR,
I = βSI − νI,
R = νI − γR.
Fakt, že S +I +R = 0, opäť korešponduje s predpokladom S(t)+I(t)+R(t) = N o konštantnej
veľkosti populácie. Vzťah R = N − S − I môžeme dosadiť do prvej rovnice, čo nás privedie k
autonómnemu systému v rovine:
S = −βSI + γ(N − S − I),
I = βSI − νI.
Systém má jednu S-nulklinu definovanú rovnicou
I =
γ(N − S)
βS + γ
= −
γ
β
+
γ
β
γ
β
+ N
S + γ
β
a dve I-nulkliny definované rovnicami
I = 0, S =
ν
β
.
Tie sa pretínajú v dvoch stacionárnych bodoch:
[S1, I1] = [N, 0], [S2, I2] =
ν
β
,
γ(βN − ν)
β(ν + γ)
.
Z druhej rovnice vidno, že ak S ≤ ν/β, potom I ≤ 0 a I preto klesá. Nutnou podmienku
pre vypuknutie epidémie je σ = βN/ν > 1, čo budeme ďalej aj predpokladať. Číslo σ sa
7
nazýva infekčné kontaktné číslo alebo vnútorná reproduktívna rýchlosť infekcie. Dôsledkom
nášho predpokladu je fakt, že stacionárny bod [S2, I2] leží vo vnútri prvého kvadrantu.
Jacobiho matica systému má tvar
J(S, I) =
−βI − γ −βS − γ
βI βS − ν
.
Vzhľadom k tomu, že det J(N, 0) = −γ(βN − ν) < 0, je stacionárny bod [S1, I1] sedlom. Ďalej
platí
tr J(S2, I2) = −
γ(βN − ν)
ν + γ
− γ = −
γ(βN + γ)
ν + γ
< 0,
det J(S2, I2) = (ν + γ)
γ(βN − ν)
ν + γ
= γ(βN − ν) > 0,
(tr J(S2, I2))2
− 4 det J(S2, I2) = γ
γ(βN + γ)2
(ν + γ)2
− 4(βN − ν) .
To znamená, že bod [S2, I2] je v každom prípade stabilný a je to ohnisko pre malé γ a uzol pre
veľké γ. Všimnime si ešte, že
S2 + I2 =
ν
β
+
γ(βN − ν)
β(ν + γ)
=
γβN + ν2
β(ν + γ)
=
γN + ν2
β
ν + γ
<
γN + νN
ν + γ
= N,
čo znamená, že bod [S2, I2] leží vo vnútri trojuholníka K = {(x, y) ∈ R2
| x, y ≥ 0 x + y ≤ N}.
pokúsime sa ukázať, že tento trojuholník je pozitívne invariantná množina.
Vzhľadom k tomu, že os S je I-nulklinou, je to invariantná množina systému. Inými slovami,
S(t) = (S0 − N)e−γt
+ N, I(t) = 0
je riešenie. Na druhej strane, na ose I má prvá zložka vektorového poľa tvar γ(N − I), čo je
nezáporný výraz pre I ≤ N. Trajektórie začínajúce na ose I pre 0 ≤ I ≤ N teda vchádzajú do
trojuholníka K. Nakoniec sa pozrieme, ako sa mení veličina S + I pozdĺž riešení:
(S(t) + I(t)) = γ(N − S(t) − I(t)) − νI(t)
Na prepone trojuholníka K platí S(t) + I(t) = N, takže v týchto bodoch je posledný výraz
záporný, okrem stacionárneho bodu [N, 0]. Preto trajektórie začínajúce na prepone trojholníka
K do neho tiež vstupujú. Ukázali sme, že K je pozitívne invariantná množina. Navyše, je to
tiež kompaktná množina.
8
Pomocou Dulacovho kritéria ukážeme, že v trojuholníku K neexistuje uzavretá trajektória.
Počítajme:
∂
∂S
1
I
(−βSI + γ(N − S − I)) +
∂
∂I
1
I
(βSI − νI)
=
∂
∂S
−βS + γ
N − S − I
I
+
∂
∂I
(βS − ν)
= −β −
γ
I
.
Posledný výraz je vo vnútri trojuholníka K záporný, takže K neobsahuje uzavretú trajektóriu.
Zvoľme ľubovoľné riešenie (S(t), I(t)) spĺňajúce (S(0), I(0)) ∈ K a označme symbolom C+
jeho kladnú polotrajektóriu, teda C+
= {(S(t), I(t)) | t ≥ 0}. Vzhľadom k tomu, že K je
pozitívne invariantná množina, platí C+
⊆ K. V dôsledku toho má C+
kompaktný uzáver.
Predpokladajme, že ω-limitná nnožina Ω(C+
) neobsahuje žiadny stacionárny bod. Pripomeňme,
že Ω(C+
) ⊆ C+ ⊆ K. Potom podľa Poincarého-Bendixsonovej vety ([1], strana 11) je
Ω(C+
) cyklus, čo je v spore s tým, že v K žiadny cyklus nie je. Preto Ω(C+
) obsahuje aspoň
jeden stacionárny bod. Ak je to stabilný stacionárny bod [S2, I2], potom zrejme platí
Ω(C+
) = {[S2, I2]} a platí limt→∞(S(t), I(t)) = [S2, I2]. V opačnom prípade množina Ω(C+
)
obsahuje sedlo [S1, I1]. Podľa Vety 8 o štruktúre ω-limitnej množiny rovinného autonómneho
systému so stacionárnymi bodmi ([1], strana 11) je Ω(C+
) zjednotením stacionárnych bodov
a heteroklinických a homoklinických trajektórií medzi nimi. Heteroklinickú trajektóriu neobsahuje,
lebo v takom prípade by obsahovala aj bod [S2, I2]. Sedlo [S1, I1] nemá homoklinickú
trajektóriu, keďže riešenia, ktoré sa k nemu blížia pre t → ∞, ležia na osi S. V konečnom
dôsledku platí buď Ω(C+
) = {[S2, I2]} alebo Ω(C+
) = {[S1, I1]}, teda buď
lim
t→∞
(S(t), I(t)) = [S2, I2], alebo lim
t→∞
(S(t), I(t)) = [S1, I1],
pričom druhá možnosť nastane len v prípade I(0) = 0. Náš počiatočný problém S(0) = S0 > 0,
I(0) = I0 > 0, R(0) = R0 ≥ 0 a S0 + I0 + R0 = N sa teda ustáli na hodnote
ν
β
,
γ(βN − ν)
β(ν + γ)
,
ν(βN − ν)
β(ν + γ)
=
ν
β
,
γν(σ − 1)
β(ν + γ)
,
ν2
(σ − 1)
β(ν + γ)
=
ν
β
,
N − ν/β
1 + ν/γ
,
N − ν/β
1 + γ/ν
.
Podobne ako v prípade modelu SIR, rozsah epidémie môžeme zmenšiť buď zväčšením parametra
ν alebo zmenšením parametra β.
Literatúra
[1] KALAS, Josef, POSPÍŠIL, Zdeněk. Spojité Modely v Biologii. 1. vyd. Brno: Masarykova
univerzita, 2001.
9