Západočeská univerzita v Plzni Fakulta aplikovaných věd Katedra matematiky. Matematické modely v epidemiologii

Západoˇceská univerzita v Plzni Fakulta aplikovaných vˇed Katedra matematiky

Bakaláˇrská práce Matematick´ e modely v epidemiologii

Plzeˇn, 2015

Eva Beranov´ a

Prohl´ aˇ sen´ı Prohlaˇsuji, ˇze jsem bakaláˇrskou práci vypracovala samostatnˇe a v´ yhradnˇe s pouˇzit´ım citovan´ ych pramen˚ u. V Plzni dne 19. kvˇetna 2015 Eva Beranová

Podˇ ekov´ an´ı Ráda bych podˇekovala doc. Ing. Gabriele Holubové, Ph.D. za cenné rady, vˇecné pˇripom´ınky a vstˇr´ıcnost pˇri konzultac´ıch a vypracováván´ı bakaláˇrské práce.

Abstrakt Tato práce uvád´ı pˇrehled vybran´ ych epidemiologick´ ych model˚ u, dále ˇreˇs´ı modifikace základn´ıho SIR modelu a na pˇr´ıkladu porovnává teoretické v´ ysledky se znám´ ym pˇr´ıpadem ˇs´ıˇren´ı spálové ang´ıny. Zamˇeˇruje se na z´ıskán´ı informac´ı o vrcholu epidemie a porovnán´ı numerick´ ych v´ ysledk˚ u pro r˚ uzná vyjádˇren´ı SIR modelu. Pˇr´ınosem práce je pak pˇredevˇs´ım pˇreveden´ı SIR modelu na rovnici Abelova typu, z n´ıˇz bylo následnˇe vyjádˇreno ˇreˇsen´ı SIR soustavy v parametrickém tvaru.

Abstract This thesis contains an overview of selected epidemiological models, it also discusses the modifications of the basic SIR model and compares the theoretical results with known spread of scarlatinal tonsillitis. It focuses on obtaining information about the peak of the epidemic and comparing numerical results for different expressions of the SIR model. The main contribution of this work is the transfer of SIR model to Abel’s type differential equation from which the parametric solution of the SIR system was sebsequently expressed.

Kurzfassung ¨ Diese Arbeit enthält einen Uberblick von ausgewählten epidemiologischen Modellen, sie löst auch die Modifikationen von SIR-Modell und vergleicht theoretische Resultate mit dem bekannten Fall der Ausbreitung der Scharlach-

Angina. Sie orientiertet sich auf den Erwerb der Informationen u ¨ber die Spitze der Epidemie und den Vergleich von numerischen Ergebnissen f¨ ur ver¨ schiedene Außerungen von SIR-Modell. Der Beitrag der Arbeit ist vor allem ¨ die Ubertragung von SIR-Modell auf die abelsche Gleichung, aus der nachfolgend die Lösung des SIR-Systems in der Parameterform ausgedr¨ uckt war.

Obsah 1 Pˇ rehled vybran´ ych pojm˚ u a z´ akladn´ıch model˚ u 1.1 Základn´ı pˇredpoklady a pojmy . . . . . 1.1.1 Dynamické systémy . . . . . . . 1.1.2 Epidemiologie . . . . . . . . . . 1.2 Základn´ı SIR model . . . . . . . . . . . 1.3 SIR model s vitáln´ı dynamikou . . . . 1.4 Rozˇs´ıˇrené SIR modely . . . . . . . . . 1.4.1 Vakcinace . . . . . . . . . . . . 1.4.2 Karanténa . . . . . . . . . . . . 1.4.3 Obecná m´ıra kontaktu . . . . . 1.4.4 Dalˇs´ı modifikace . . . . . . . . . 1.5 SI model . . . . . . . . . . . . . . . . . 1.6 SIS model . . . . . . . . . . . . . . . . 1.7 SIRI a SIRS modely . . . . . . . . . . 1.8 SEIR a SEIS modely . . . . . . . . . . 1.9 SIHR model . . . . . . . . . . . . . . . 1.10 SITR model . . . . . . . . . . . . . . . 1.11 SIRD model . . . . . . . . . . . . . . . 1.12 Shrnut´ı pˇrehledu model˚ u . . . . . . . .

epidemiologick´ ych . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

2 Anal´ yza z´ akladn´ıho SIR modelu 2.1 Pˇreveden´ı modelu na rovnici vyˇsˇs´ıho ˇrádu . . . . . . . . . . . 2.2 Pˇreveden´ı modelu na soustavu rovnic s parametrem . . . . . . 2.2.1 Vrchol epidemie . . . . . . . . . . . . . . . . . . . . . . 2.3 Pˇreveden´ı modelu na rovnici Abelova typu se separovan´ ymi promˇenn´ ymi . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Porovnán´ı numerick´ ych v´ ysledk˚ u jednotliv´ ych systém˚ u pro r˚ uzn´ y software . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Základn´ı reprodukˇcn´ı ˇc´ıslo SIR modelu a celkov´ y rozsah epidemie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15 15 17 18 19 22 24 24 24 28 29 29 30 31 33 36 37 38 39 41 41 43 46 52 57 64

3 Porovn´ an´ı re´ aln´ ych dat s teoretick´ ym 3.1 Spálová ang´ına . . . . . . . . . . . . 3.1.1 Volba modelu . . . . . . . . . 3.2 Modelován´ı pr˚ ubˇehu nákazy . . . . . ´ 3.2.1 Uskal´ ı reáln´ ych dat . . . . . . A

modelem . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

66 66 66 68 70

72 A.1 Abelova rovnice . . . . . . . . . . . . . . . . . . . . . . . . . . 72 A.2 Bernoulliho rovnice . . . . . . . . . . . . . . . . . . . . . . . . 73

´ Uvod Epidemie naˇsi populaci velmi ovlivˇ nuj´ı, proto se lidstvo snaˇz´ı z´ıskat poznatky o ˇs´ıˇren´ı nemoc´ı a jejich pˇrenosu a modelovat pr˚ ubˇeh epidemi´ı budouc´ıch. V praxi z´ıskaná data jsou ˇcasto nekompletn´ı a natolik nepˇresná, ˇze je nejsme schopni interpretovat a chován´ı epidemie m˚ uˇze p˚ usobit chaoticky. Také odhad parametr˚ u a souˇcást´ı modelu je velmi tˇeˇzk´ y, proto se nab´ız´ı otázka, zda je pro nás matematické modelován´ı v epidemiologii v˚ ubec pˇr´ınosné. S jeho pomoc´ı jsme ovˇsem schopni porozumˇet skryt´ ym mechanism˚ um, které ovlivˇ nuj´ı ˇs´ıˇren´ı nemoci a navrhovat strategie, jak epidemie dostat pod kontrolu. Modely ˇcasto dokáˇz´ı rozeznat i chován´ı, které je z experimentáln´ıch dat nejasné. V´ yzkum infekˇcn´ıch chorob m˚ uˇze b´ yt rozdˇelen na deskriptivn´ı, analytick´ y, experimentáln´ı a teoretick´ y. Matematické modely jsou zaloˇzeny na populaˇcn´ı dynamice, chován´ı a pˇrenosu nemoci, vlastnostech infekˇcn´ıch agens a na dalˇs´ıch sociáln´ıch a psychologick´ ych faktorech. Pˇri matematickém modelován´ı pˇrenosu nemoc´ı je vˇzdy nutné zvolit kompromis mezi jednoduch´ ymi modely, které vynechávaj´ı vˇetˇsinu podrobnost´ı a jsou vytvoˇreny jen pro osvˇetlen´ı obecného kvalitativn´ıho chován´ı nákazy, a detailn´ımi modely, obvykle vytvoˇren´ ymi pro specifick´ y pˇr´ıpad a obsahuj´ıc´ı’ mi krátkodobou kvantitativn´ı pˇredpovˇed . Detailn´ı modely je obecnˇe sloˇzité nebo nemoˇzné ˇreˇsit analyticky a jejich vyuˇzit´ı pro teoretické u ´ˇcely je omezené, pˇrestoˇze jejich strategická hodnota m˚ uˇze b´ yt vysoká. Pˇri ˇreˇsen´ı reáln´ ych situac´ı je totiˇz pro zdravotnické odborn´ıky jednoduch´ y model nedostateˇcn´ y a je nutná numerická simulace modelu detailn´ıho. Spojen´ı epidemiologické teorie, biostatistiky a poˇc´ıtaˇcov´ ych simulac´ı pˇrispˇelo ke zlepˇsen´ı naˇsich vˇedomost´ı o mechanismech pˇrenosu epidemi´ı, k v´ yvoji epidemiologie a ke vzniku efektivnˇejˇs´ıch metod kontrolován´ı infekˇcn´ıch cho13

rob. Modelován´ı tˇechto nemoc´ı odhalilo mnoˇzstv´ı zaj´ımav´ ych jev˚ u a odliˇsnost´ı v jejich chován´ı a zapˇr´ıˇcinilo vznik celé ˇrady podnˇet˚ u pro dalˇs´ı v´ yzkum. Tato práce je ˇclenˇena následovnˇe: kapitola 1 slouˇz´ı k seznámen´ı se s vybran´ ymi pojmy z teorie dynamick´ ych systém˚ u a epidemiologie a uvád´ı pˇrehled základn´ıch epidemiologick´ ych model˚ u. V kapitole 2 je pak provedena anal´ yza základn´ıho SIR modelu, jsou uvaˇzovány moˇznosti pˇreveden´ı modelu do jiné formy a dále jsou srovnávány numerické v´ ysledky tˇechto nov´ ych systém˚ u v softwarech Matlab a Mathematica. Kapitola 3 poukazuje na komplikovanost modelován´ı v epidemiologii srovnán´ım reáln´ ych dat z lokáln´ı epidemie spálové ang´ıny s teoretick´ ymi hodnotami z´ıskan´ ymi z jednoduchého modelu.

14

Kapitola 1 Pˇ rehled vybran´ ych pojm˚ u a z´ akladn´ıch epidemiologick´ ych model˚ u V této kapitole pˇredstav´ıme struˇcn´ yu ´vod do modelován´ı v epidemiologii, pˇredpoklady, z kter´ ych modelován´ı vycház´ı, vysvˇetl´ıme nˇekteré pojmy z epidemiologie a teorie dynamick´ ych systém˚ u, pˇredstav´ıme vybrané modely a shrneme moˇznosti, které matematické modelován´ı v epidemiologii nab´ız´ı. Pokud nen´ı uvedeno jinak, informace v této kapitole byly ˇcerpány z [2]. Ilustraˇcn´ı obrázky byly vytvoˇreny v programu MATLAB pomoc´ı funkce ode45.

1.1

Z´ akladn´ı pˇ redpoklady a pojmy

Modely b´ yvaj´ı pojmenovány pomoc´ı zkratek anglick´ ych názv˚ u oznaˇcuj´ıc´ıch jednotlivé modelové tˇr´ıdy - skupiny populace: • S susceptible - náchyln´ı jedinci Tito jedinci nemaj´ı proti nemoci ˇzádnou imunitu a mohou se snadno nakazit. • I infectious - infekˇcn´ı jedinci Nakaˇzen´ı, kteˇr´ı mohou po celou dobu onemocnˇen´ı ˇs´ıˇrit nemoc dál mezi náchylné. • R recovered - imunn´ı jedinci Prodˇelán´ım nemoci, pˇr´ıpadnˇe oˇckován´ım, z´ıskali imunitu, která m˚ uˇze b´ yt trvalá, nebo doˇcasná. 15

• D dead - zemˇrel´ı jedinci Vyskytuj´ı se pˇri modelován´ı smrteln´ ych nemoc´ı. • E exposed - infikovan´ı jedinci Vyskytuj´ı se pˇri modelován´ı nemoc´ı s latentn´ım obdob´ım, v nˇemˇz je osoba infikovaná, ovˇsem bez pˇr´ıznak˚ u nemoci, a schopnost pˇrenáˇset nemoc na náchylné jedince maj´ı tito lidé jen ˇcásteˇcnˇe, nebo v˚ ubec ne. • H hibernator Infekˇcn´ı jedinci, kteˇr´ı z˚ ustávaj´ı doma a nemoc tedy dále neˇs´ıˇr´ı. • T treatment - léˇcen´ı jedinci ˇ ast infekˇcn´ıch jedinc˚ C´ u, které se dostalo léˇcby, a zkrátila se tak doba onemocnˇen´ı. Matematick´ ym modelem epidemiologického jevu pak budeme rozumˇet nelineárn´ı systém obyˇcejn´ ych diferenciáln´ıch rovnic. Poˇcet rovnic odpov´ıdá poˇctu pouˇzit´ ych typ˚ u modelov´ ych tˇr´ıd. V naˇsem pˇr´ıpadˇe z tohoto systému z´ıskáváme pouze numerickou aproximaci ˇreˇsen´ı. Nezávislou promˇennou tˇechto kompartmentov´ ych model˚ u je ˇcas t. Souˇcet pouˇzit´ ych modelov´ ych tˇr´ıd mus´ı pokr´ yt celou populaci N. Tˇr´ıdy je dále moˇzné dˇelit podle poˇzadovan´ ych kritéri´ı (pohlav´ı, vˇek apod.). Poˇzadujeme, aby velikost jednotliv´ ych tˇr´ıd byla dostateˇcnˇe velká a m´ıˇsen´ı jej´ıch ˇclen˚ u by tak mohlo b´ yt povaˇzováno za homogenn´ı. Tento pˇredpoklad nen´ı platn´ y na zaˇca´tku epidemie pˇri prvotn´ım vypuknut´ı nákazy, kdy se vyskytuje pouze nˇekolik málo infikovan´ ych jedinc˚ u a pˇrenos infekce je stochastickou událost´ı závislou na kontaktech mezi jednotlivci. Pro modelován´ı poˇcáteˇcn´ı fáze nákazy se proto vyuˇz´ıvá jin´ y pˇr´ıstup neˇz ten, kter´ y je uveden´ y v této práci. Pˇredpokládáme, ˇze pr˚ ubˇeh epidemie je deterministick´ y, tud´ıˇz ˇze chován´ı populace je urˇceno jej´ı minulost´ı a pravidly popisuj´ıc´ımi model. Slabinou tohoto modelován´ı je skuteˇcnost, ˇze vn´ımáme populaci jako celek a pˇredpokládáme jej´ı homogennost, tedy ˇze uvaˇzujeme jedin´ y typ pr˚ umˇerného jedince, kter´ y se v populaci vyskytuje. Abychom zahrnuli kaˇzdého jednotlivce a odliˇsnosti jeho chován´ı od pr˚ umˇeru, potˇrebovali bychom daleko sloˇzitˇejˇs´ı model. D˚ uleˇzitou veliˇcinou je základn´ı reprodukˇcn´ı ˇc´ıslo R0 , jeˇz m˚ uˇzeme definovat jako poˇcet druhotn´ ych infekc´ı zp˚ usoben´ ych jedin´ ym infekˇcn´ım ˇcinitelem, 16

kter´ y vstoupil mezi ˇcleny kompletnˇe náchylné populace. Epidemie propukne právˇe tehdy, kdyˇz R0 > 1. Vyˇsˇs´ı R0 vede k rychlejˇs´ımu nástupu epidemie a k vˇetˇs´ımu poˇctu nemocn´ ych (tj. k závaˇznˇejˇs´ı epidemii celkovˇe), urˇcuje proporci zdrav´ ych jedinc˚ u na vrcholu epidemie (1/R0 ) a prevalenci nemoci na vrcholu epidemie. Hodnotu R0 je moˇzné odhadnout po skonˇcen´ı epidemie nebo bˇehem jej´ıho pr˚ ubˇehu. Odborn´ıci se snaˇz´ı navrhovat v praxi takové postupy, které povedou ke sn´ıˇzen´ı reprodukˇcn´ıho ˇc´ısla pod hraniˇcn´ı hodnotu a t´ım dojde k zamezen´ı vzniku epidemie.

1.1.1

Dynamick´ e syst´ emy

Následuj´ıc´ı vybrané pojmy z teorie dynamick´ ych systém˚ u byly ˇcerpány z [9]. Autonomn´ı soustavu obyˇcejn´ ych diferenciáln´ıch rovnic prvn´ıho ˇra´du x01 = f1 (x1 , x2 , ..., xn ), x02 = f2 (x1 , x2 , ..., xn ), .. . 0 xn = fn (x1 , x2 , ..., xn ), kde f1 , f2 , ..., fn jsou reálné funkce n reáln´ ych promˇenn´ ych, m˚ uˇzeme vektorovˇe zapsat jako x0 = f (x), kde x = (x1 , x2 , ..., xn ) ∈ Rn a f : Rn → Rn . Soustavu lze interpretovat jako matematický model nˇejakého reálného dynamického systému se závislou promˇennou x a nezávislou promˇennou t, jehoˇz stav je v ˇcase t urˇcen uspoˇra´danou n-tic´ı ˇc´ısel x(t) = (x1 (t), x2 (t), ..., xn (t)). Potom x0 (t) =

dx = (x01 (t), x02 (t), ..., x0n (t)) dt

udává rychlost zmˇeny stavu systému, která je zadána prostˇrednictv´ım funkce f , jej´ıˇz hodnota závis´ı pouze na okamˇzitém stavu systému. Jestliˇze funkce f závis´ı také na ˇcase t, tj. f = f (t, x), pak rovnici 17

x0 = f (t, x) naz´ yváme neautonomn´ı. Jestliˇze je funkce f lineárn´ı v x, tj. f (x, t) = A(t)x + b(t), kde A(t) je maticová funkce a b(t) vektorová funkce, naz´ yváme soustavu lineárn´ım modelem. Plat´ı f (x+y) = f (x)+f (y) a f (λx) = λf (x), ∀x, y ∈ Rn , ∀λ ∈ R. Cauchyovou (poˇcáteˇcn´ı) u ´lohou rozum´ıme u ´lohu, kde hledáme ˇreˇsen´ı x 0 soustavy x = f (t, x) splˇ nuj´ıc´ı pˇredem zadanou podm´ınku x(t0 ) = x0 . Dan´ y bod [t0 , x0 ] je tzv. poˇcáteˇcn´ı podm´ınka. yváme equilibriem (rovnováˇzn´ ym stavem) soustavy, jestliˇze Bod x∗ naz´ f (x ) = 0. Stav systému se v tomto bodˇe v ˇcase nemˇen´ı, tj. (x∗ )0 = 0. Equilibrium je stabiln´ı, pokud jemu bl´ızká“ ˇreˇsen´ı z˚ ustanou bl´ızká“ i ve veˇske” ” rém následuj´ıc´ım ˇcase, tj. pro vˇsechna okol´ı O equilibria x∗ existuje okol´ı O1 takové, ˇze kaˇzdé ˇreˇsen´ı x(t) s x(0) = x0 v O1 je definované a setrvává v O pro vˇsechna t > 0. Pokud nav´ıc m˚ uˇzeme zvolit O1 tak, ˇze limt→∞ x(t) = x∗ , jedná se o asymptotickou stabilitu. ∗

1.1.2

Epidemiologie

Modelován´ı pˇrenosu a pr˚ ubˇehu nemoc´ı se neobejde bez znalosti základn´ıch pojm˚ u z epidemiologie. Následuj´ıc´ı pˇrehled vycház´ı pˇredevˇs´ım z [8].

• epidemiologie - vˇeda zab´ yvaj´ıc´ı se studiem v´ yskytu chorob a faktor˚ u, které jej ovlivˇ nuj´ı • infekˇcn´ı agens - p˚ uvodce nákazy • patogen - chorobn´ y ˇcinitel, obvykle choroboplodn´ y zárodek • hostitel - organismus, v jehoˇz tˇele se nacházej´ı patogeny • infekce - proniknut´ı choroboplodn´ ych zárodk˚ u do organismu • epidemie - náhlé vypuknut´ı a hromadné ˇs´ıˇren´ı infekˇcn´ı nemoci v populaci, jej´ıˇz znaˇcná ˇca´st je postiˇzena 18

• pandemie - rozsáhlá epidemie • endemie - v´ yskyt infekˇcn´ıho onemocnˇen´ı ohraniˇcen´ y na urˇcitou oblast • prevalence - poˇcet vˇsech pˇr´ıpad˚ u urˇcitého onemocnˇen´ı vztaˇzen´ y obvykle na 100 000 obyvatel a kalendáˇrn´ı rok Obdob´ı, kdy se patogen snaˇz´ı mnoˇzit a hostitel se jej snaˇz´ı zbavit, se naz´ yvá dynamick´ y stav. Interakce hostitel-patogen prob´ıhá r˚ uznˇe na základˇe konkrétn´ıho druhu patogenu a hostitele a jej´ım v´ ysledkem nemus´ı vˇzdy b´ yt vyvolán´ı nemoci u hostitele.

1.2

Z´ akladn´ı SIR model

Tento model vytvoˇrili A. G. McKendrick a W. O. Kermack roku 1927. Jde o kompartmentov´ y model zaloˇzen´ y na relativnˇe jednoduch´ ych pˇredpokladech o rychlosti toku mezi r˚ uzn´ ymi tˇr´ıdami ˇclen˚ u populace. Neuvaˇzuj´ı se zde demografické zmˇeny (natalita, mortalita), z´ıskaná imunita je povaˇzována za trvalou a nemoc nen´ı smrtelná. Populace je homogenn´ı a pravdˇepodobnost nakaˇzen´ı je konstantn´ı. Vˇsechny dalˇs´ı modely m˚ uˇzeme vn´ımat jen jako modifikaci tohoto základn´ıho modelu, ilustrovaného obrázkem 1.1 a popsaného rovnicemi dS = −βSI, dt dI = βSI − γI, dt dR = γI, dt

(1.2.1) (1.2.2) (1.2.3)

s poˇca´teˇcn´ımi podm´ınkami: S(0) = S0 ∈ R+ , S0 + I0 + R0 = N,

I(0) = I0 ∈ R+ , R+ = h0; +∞).

R(0) = R0 ∈ R+ ,

Konstanty β a γ udávaj´ı m´ıry pˇrechodu mezi tˇr´ıdami a m˚ uˇzeme je uvaˇzovat jako pravdˇepodobnosti, tedy 0 ≤ β ≤ 1 a 0 ≤ γ ≤ 1 . M´ıra infekce β popisuje pravdˇepodobnost nakaˇzen´ı pˇri kontaktu mezi náchyln´ ymi a infekˇcn´ımi jedinci S a I, respektive ˇr´ıká, ˇze pr˚ umˇern´ y ˇclen populace se za jednotku ˇcasu dostane do kontaktu dostateˇcného k infikován´ı s βN dalˇs´ımi jedinci. Jej´ı vliv na model je moˇzné pozorovat na obrázku 1.2. Je tˇeˇzké ji urˇcit, protoˇze 19

Obrázek 1.1: Graf pro základn´ı SIR model (β = 0.003, γ = 0.3). nezávis´ı jen na studované nemoci, ale také na sociáln´ıch faktorech a chován´ı jednotlivc˚ u. Konstanta γ udává pr˚ umˇernou rychlost uzdraven´ı a pˇrechod z tˇr´ıdy infekˇcn´ıch jedinc˚ u I mezi imunn´ı R, infekˇcn´ı jedinci tedy opouˇstˇej´ı 1 tˇr´ıdu I rychlost´ı γI za jednotku ˇcasu a doba trván´ı infekce d je . Vliv této γ konstanty na model je ilustrovan´ y obrázkem 1.3. Model dává smysl pouze dokud jsou S a I nezáporné. Pokud jedna z tˇechto tˇr´ıd dosáhne nuly, je rozumné modelován´ı ukonˇcit. Vˇsimnˇeme si, dS dI ˇze <0a > 0 právˇe tehdy, kdyˇz S je vˇetˇs´ı neˇz γ/β. Hodnota I tud´ıˇz dt dt roste tak dlouho, dokud je S > γ/β, ale jakmile S klesne pod hodnotu γ/β, I klesá téˇz, aˇz dosáhne nuly. Jestliˇze je na zaˇcátku modelován´ı poˇcet náchyln´ ych jedinc˚ u v ˇcase nula S0 menˇs´ı neˇz γ/β, poˇcet infekˇcn´ıch jedinc˚ u klesne k nule a k epidemii nedojde. Aby epidemie probˇehla, mus´ı platit, ˇze S0 > γ/β. Poˇcet infikovan´ ych jedinc˚ u I pak roste, aˇz dosáhne maxima pro S = γ/β a zaˇcne klesat k nule. β je hraniˇcn´ı hodnota, odpov´ıdaj´ıc´ı základn´ımu reprodukˇcn´ımu γ ˇc´ıslu R0 . Pokud je R0 > 1, vypukne epidemie. Odvozen´ı závislosti celkového rozsahu epidemie na R0 je uvedeno v kapitole 2.5. Veliˇcina S0

20

Soustavu rovnic (1.2.1), (1.2.2) a (1.2.3) m˚ uˇzeme pˇrevést i na soustavu diferenˇcn´ıch rovnic Sn+1 = Sn − βSn In , In+1 = In + βSn In − γIn , Rn+1 = Rn + γIn ,

(1.2.4) (1.2.5) (1.2.6)

kde Sn = S(n), tj. S v ˇcase t = n, In = I(n), Rn = R(n) a plat´ı S0 > 0, I0 > 0, R0 > 0 . Tento diskrétn´ı model je vhodn´ y po porovnán´ı simulovan´ ych dat se statistick´ ymi daty nebo v poˇca´tc´ıch epidemie, obdobnˇe lze pˇrepsat vˇetˇsinu spojit´ ych model˚ u. Studium diskrétn´ıch model˚ u vˇsak nen´ı obsahem této práce.

Obrázek 1.2: Porovnán´ı graf˚ u SIR modelu pro r˚ uzné hodnoty β (vlevo β = 0.001, vpravo β = 0.004, γ = 0.3 u obou graf˚ u).

Obrázek 1.3: Porovnán´ı graf˚ u SIR modelu pro r˚ uzné hodnoty γ (vlevo γ = 0.1, vpravo γ = 0.7, β = 0.003 u obou graf˚ u).

21

1.3

SIR model s vit´ aln´ı dynamikou

Na rozd´ıl od základn´ıho SIR modelu uvaˇzuje model s vitáln´ı dynamikou i natalitu a mortalitu, ty jsou vˇsak v rovnováze, a proto nedocház´ı ke zmˇenˇe celkové populace N. Vyuˇz´ıvá se pˇredevˇs´ım pro endemické nemoci, které p˚ usob´ı po delˇs´ı ˇcas a demografické zmˇeny jsou d´ıky tomu nezanedbatelné. Model vypadá následovnˇe: dS = −βSI + µ(N − S), dt dI = βSI − (γ + µ)I, dt dR = γI − µR, dt

(1.3.1) (1.3.2) (1.3.3)

kde µ je stejnomˇerná m´ıra narozen´ı a m´ıra u ´mrt´ı, jej´ıˇz vliv na model m˚ uˇzeme pozorovat na obrázku 1.4. Jsou zde pouˇzity stejné hodnoty pro konstanty γ a β jako v pˇr´ıpadˇe základn´ıho SIR modelu na obrázku 1.1, nav´ıc pˇribyla konstanta µ s hodnotou 0.1. Ta zp˚ usobuje prodlouˇzen´ı epidemie, nebot’ v populaci neustále pˇrib´ yvaj´ı náchyln´ı jedinci. Pˇri pouˇzit´ı niˇzˇs´ı hodnoty µ (viz obrázek 1.5) nebude celkov´ y rozd´ıl model˚ u tolik patrn´ y, pˇresto si m˚ uˇzeme vˇsimnout poklesu poˇctu imunn´ıch jedinc˚ u a pˇr´ır˚ ustku jedinc˚ u náchyln´ ych.

22

Obrázek 1.4: Graf pro SIR model s vitáln´ı dynamikou (β = 0.003, γ = 0.3, µ = 0.1).

Obrázek 1.5: Graf pro SIR model s vitáln´ı dynamikou (β = 0.003, γ = 0.3, µ = 0.01).

23

1.4

Rozˇ s´ıˇ ren´ e SIR modely

Do modelu je moˇzné zahrnout dalˇs´ı jevy, které tvoˇr´ı v´ ysledky reálnˇejˇs´ımi a které lze libovolnˇe kombinovat.

1.4.1

Vakcinace

Oˇckován´ım pˇred vypuknut´ım epidemie je ˇca´st náchylné populace pˇresunuta pˇr´ımo mezi imunn´ı jedince. V tomto pˇr´ıpadˇe staˇc´ı pouˇz´ıt základn´ı SIR model a pouze upravit poˇca´teˇcn´ı velikosti jednotliv´ ych tˇr´ıd. Pokud oˇckován´ı zasahuje do pr˚ ubˇehu epidemie, pouˇz´ıvá se následuj´ıc´ı u ´prava modelu: dS = −(β − ψS)SI, dt dI = (β − ψS)SI − γI, dt dR = γI + ψS, dt

(1.4.1) (1.4.2) (1.4.3)

m´ıra vakcinace je závislá na poˇctu náchyln´ ych osob ψS. Vakcinace zp˚ usob´ı pomalejˇs´ı nástup epidemie a celkovˇe sn´ıˇz´ı poˇcet postiˇzen´ ych jedinc˚ u, doba trván´ı epidemie se ale prodlouˇz´ı (viz obrázek 1.6). V reálu ovˇsem nem˚ uˇzeme oˇcekávat, ˇze by oˇckován´ı mˇelo 100% u ´spˇeˇsnost.

1.4.2

Karant´ ena

Karanténa se vyuˇz´ıvá ke sn´ıˇzen´ı mnoˇzstv´ı nákaz osob prostˇrednictv´ım izolace infekˇcn´ıch nebo náchyln´ ych jedinc˚ u. V´ ysledkem je pokles hodnoty prevalence na vrcholu epidemie a sn´ıˇzen´ı poˇctu zasaˇzen´ ych jedinc˚ u. Karanténu m˚ uˇzeme modelovat nˇekolika r˚ uzn´ ymi zp˚ usoby. Jednou z moˇznost´ı je uvaˇzován´ı ne´ uplné karantény, kdy v ˇcase tk dojde k v´ yznamnému sn´ıˇzen´ı m´ıry kontaktu mezi náchyln´ ymi a infekˇcn´ımi jedinci. V praxi se jedná pˇredevˇs´ım o uzavˇren´ı ˇskol v dobˇe epidemie, tuto situaci ilustruje obrázek 1.7. Jak je z obrázku vidˇet, dopad epidemie se sn´ıˇz´ı, ale prodlouˇz´ı se jej´ı trván´ı. Dalˇs´ı moˇznost´ı je pˇresunout v ˇcase tk ˇca´st z infekˇcn´ıch jedinc˚ u (obrázek 1.8) nebo z náchyln´ ych jedinc˚ u (obrázek 1.9) do u ´plné izolace (viz. [6]). Naˇs´ım c´ılem je samozˇrejmˇe poslat do karantény co nejvˇetˇs´ı ˇcást jedinc˚ u 24

v co nejkratˇs´ım ˇcase. V tomto pˇr´ıpadˇe se do modelu obvykle pˇridává tˇr´ıda Q, Q0 = −δ, kde δ udává poˇcet jedinc˚ u, kter´ y pˇrejde za jednotku ˇcasu mezi imunn´ı R. V pˇr´ıpadˇe karantény náchyln´ ych jedinc˚ u m˚ uˇzeme tˇr´ıdu R chápat jako jedince vylouˇcené z procesu nákazy. Podoba rovnice pro R0 do ukonˇcen´ı karantény je následuj´ıc´ı: R0 = γI + δ.

Obrázek 1.6: Graf pro SIR model s vakcinac´ı (β = 0.003, γ = 0.3, ψ = 0.000003).

25

Obrázek 1.7: Graf pro SIR model s poklesem hodnoty β v ˇcase 1 na ˇctvrtinu (β1 = 0.003, β2 = 0.00075, γ = 0.3).

Obrázek 1.8: Graf pro SIR model s karanténou (β = 0.003, γ = 0.3). V ˇcase 2 jsme pˇresunuli do karantény 150 infekˇcn´ıch jedinc˚ u.

26

Obrázek 1.9: Graf pro SIR model s karanténou (β = 0.003, γ = 0.3). V ˇcase 0.5 jsme pˇresunuli do karantény 300 náchyln´ ych jedinc˚ u.

27

1.4.3

Obecn´ a m´ıra kontaktu

Realistiˇctˇejˇs´ım modelem je pˇredpokládat m´ıru kontaktu mezi skupinami S a I jako neklesaj´ıc´ı funkci celkové populace N. Tento pˇr´ıstup je vhodn´ y pˇredevˇs´ım pro pohlavnˇe pˇrenosné choroby. Vytvoˇr´ıme pˇredpoklad, ˇze pr˚ umˇern´ y ˇclen populace se dostane do C(N) 0 kontakt˚ u za jednotku ˇcasu (C (N ) ≥ 0) a definujeme β(N ) =

C(N ) . N

(1.4.4)

Standardn´ımu v´ yskytu kontakt˚ u odpov´ıdá C(N ) = λ (poˇcet kontakt˚ u nezávis´ı na velikosti populace N ), pˇri masovém v´ yskytu kontakt˚ u C(N ) = βN (poˇcet kontakt˚ u roste lineárnˇe s velikost´ı populace N ). V praxi se m˚ uˇzeme setkat také s Michaelis-Mentenov´ ym typem interakce v podobˇe C(N ) =

aN . 1 + bN

(1.4.5)

Pro kontakt ve mˇestech stˇredn´ı velikosti je vhodná forma C(N ) = λN a ,

a = 0.05.

(1.4.6)

Protoˇze je nyn´ı velikost populace vyuˇz´ıvána v modelu, mus´ıme do nˇej zahrnout i rovnici pro N 0 . To nám umoˇzn ˇuje jednoduˇse zahrnout zemˇrelé jedince do modelu, kde ˇca´st f z γI ˇclen˚ u populace pˇrejde mezi imunn´ı, zb´ yvaj´ıc´ı ˇca´st (1 - f ) na nemoc zemˇre. Parametr f tedy udává základn´ı m´ıru pˇreˇzit´ı nemoci. Rovnici pro R0 nemus´ıme uvádˇet, nebot’ R je stanoveno pomoc´ı S, I a N, které známe. Model bude vypadat následovnˇe: dS = −β(N )SI, dt dI = β(N )SI − γI, dt dN = −(1 − f )γI. dt

(1.4.7) (1.4.8) (1.4.9)

Pokud by se f rovnalo 1, celková populace by z˚ ustala konstantn´ı (oznaˇcme N = K, plat´ı K 0 = 0) a systém rovnic (1.4.7), (1.4.8) a (1.4.9) by pˇreˇsel zpˇet do základn´ıho modelu (1.2.1) a (1.2.2) s konstantou β = β(K). 28

1.4.4

Dalˇ s´ı modifikace

• Zlepˇsuj´ıc´ı se léˇcba - V základn´ım modelu pˇredpokládáme uzdraven´ı konstantn´ı rychlost´ı γ. V pr˚ ubˇehu epidemie se ale léˇcba vyv´ıj´ı: nemoc je dˇr´ıve správnˇe urˇcena a vhodné léky jsou nasazeny rychleji. M´ıra uzdraven´ı tak lineárnˇe stoupá: γ = tν.

(1.4.10)

• Denn´ı reˇzim - Pˇres den je pravdˇepodobnost nakaˇzen´ı vˇetˇs´ı, zat´ımco v noci klesá. Infekˇcn´ı konstantu tedy m˚ uˇzeme zamˇenit za sinusoidn´ı kˇrivku závislou na ˇcase, jej´ımˇz posunut´ım zajist´ıme, ˇze koeficient bude vˇzdy kladn´ y: β = ϑ sin(λt) + κ.

(1.4.11)

• Sezónnost - V základn´ım modelu pˇredpokládáme, ˇze kontakt mezi jedinci je stále stejn´ y, ve skuteˇcnosti se vˇsak v pr˚ ubˇehu roku mˇen´ı. Typicky napˇr. u dˇetsk´ ych nemoc´ı, kdy pravdˇepodobnost nákazy vzr˚ ustá se ˇskoln´ım rokem a klesá o prázdninách. Tyto sezónn´ı zmˇeny m˚ uˇzeme do modelu zahrnout u ´pravou β(t) = β0 (1 + α cos(2πt)),

(1.4.12)

kde β0 je pr˚ umˇerná m´ıra pˇrenosu infekce, α je amplituda sezónn´ı odchylky, 0 ≤ α ≤ 1, a ˇcas t je mˇeˇren v roc´ıch. Pˇredpokládáme tedy, ˇze β(t) je periodická funkce (s periodou 1 rok) a model naz´ yváme sezónnˇe buzen´ y.

1.5

SI model

Nejjednoduˇsˇs´ım modelem je SI model, kter´ y vyuˇz´ıvá pouze tˇr´ıd náchyln´ ych a infekˇcn´ıch jedinc˚ u S a I. Modelován´ı konˇc´ı onemocnˇen´ım celé populace (viz obrázek 1.10) a obvykle se pouˇz´ıvá v pˇr´ıpadˇe, ˇze nás zaj´ımá pouze nástup nemoci. V praxi se jedná o modelován´ı poˇcáteˇcn´ıch stadi´ı onemocnˇen´ı horn´ıch cest d´ ychac´ıch. dS = −βSI, dt dI = βSI. dt 29

(1.5.1) (1.5.2)

Obrázek 1.10: Graf pro SI model (β = 0.001).

1.6

SIS model

SIS model popisuje nemoci, které nezanechávaj´ı imunitu proti opakovanému infikován´ı. Jednotlivci tedy pˇrecházej´ı z tˇr´ıdy náchyln´ ych do tˇr´ıdy infekˇcn´ıch a odtud zpˇet mezi náchylné: dS = −βSI + γI, dt dI = βSI − γI. dt

(1.6.1) (1.6.2)

Zmˇenu pr˚ ubˇehu nemoci pro r˚ uzné parametry m˚ uˇzeme pozorovat na obrázku 1.11. Nemoci tohoto typu jsou obvykle bakteriáln´ıho p˚ uvodu. Dále se model pouˇz´ıvá pro velkou ˇca´st pohlavnˇe pˇrenosn´ ych chorob (kapavka), encefalitidy a nemoci zp˚ usobené hl´ısty (paraziti jako roup dˇetsk´ y a ˇskrkavka dˇetská p˚ usob´ıc´ı v tenkém stˇrevˇe, svalovec stoˇcen´ y zp˚ usobuj´ıc´ı trichinelózu, vlasovec m´ızn´ı zp˚ usobuj´ıc´ı tzv. slon´ı nemoc apod.). β Základn´ı reprodukˇcn´ı ˇc´ıslo R0 se rovná (N − γ/β) a rozhoduje o tom, γ zda se nemoc prosad´ı nebo vymiz´ı. Pokud S = γ/β, mluv´ıme o tzv. endemickém equilibriu a nemoc pˇretrvává, pˇri S = N − γ/β se jedná o tzv. beznákazové equilibrium.

30

Obrázek 1.11: Porovnán´ı graf˚ u SIS modelu pro r˚ uzné hodnoty parametr˚ u (vlevo nahoˇre β = 0.001, γ = 0.3, vpravo nahoˇre β = 0.001, γ = 0.5, vlevo dole β = 0.003, γ = 0.3, vpravo dole β = 0.003, γ = 0.5).

1.7

SIRI a SIRS modely

Pˇredpoklad trvalé imunity, vyuˇz´ıvan´ y v SIR modelu, nesplˇ nuje celá ˇrada nemoc´ı, napˇr. chˇripka, spála, záˇskrt nebo tuberkulóza. Z´ıskaná imunita po prodˇelán´ı nemoci nemus´ı b´ yt u ´plná a nˇekteˇr´ı jedinci se mohou opˇetovnˇe nakazit. To vyjadˇruje model SIRI s konstantou α, která udává m´ıru reinfekce jako d˚ usledek pouze ˇca´steˇcné imunity tˇr´ıdy R: dS = −βSI, dt dI = βSI − γI + αRI, dt dR = γI − αRI. dt

(1.7.1) (1.7.2) (1.7.3)

Jak si m˚ uˇzeme vˇsimnout na obrázku 1.12, ne´ uplná imunita epidemii prodluˇzuje a zhorˇsuje jej´ı dopad na populaci. 31

Obrázek 1.12: Graf pro SIRI model (β = 0.003, γ = 0.3, α = 0.0005). Imunita m˚ uˇze b´ yt také pouze doˇcasná a postupnˇe se m˚ uˇze sniˇzovat. To zachycuje SIRS model: dS = −βSI + θR, (1.7.4) dt dI = βSI − γI, (1.7.5) dt dR = γI − θR, (1.7.6) dt kde θ je proporcionáln´ı m´ıra ztráty imunity. Tento pˇr´ıpad ilustruje obrázek 1.13. V porovnán´ı s obrázkem 1.12 pro SIRI model nen´ı epidemie tak závaˇzná a na vrcholu vykazuje niˇzˇs´ı prevalenci. Model m˚ uˇzeme dále upravit zaveden´ım pevnˇe dané konstantn´ı délky doˇcasné imunity ω, po které se imunn´ı jedinec vrát´ı mezi náchylné, do podoby dS = −βS(t)I(t) + γI(t − ω), dt dI = βS(t)I(t) − γI(t), dt dR = γI(t) − γI(t − ω). dt 32

(1.7.7) (1.7.8) (1.7.9)

Obrázek 1.13: Graf pro SIRS model (β = 0.003, γ = 0.3, θ = 0.1).

1.8

SEIR a SEIS modely

Nemoci jako napˇr. tuberkulóza, syfilis, borelioza a nemoci vyvolávané herpetick´ ymi viry (opary, plané neˇstovice, Kaposiho sarkom, mononukleózy, ˇsestá nemoc) se vyznaˇcuj´ı tzv. latentn´ım stadiem, kdy je hostitel sice nakaˇzen, ale jeho tˇelo neobsahuje dostatek patogenn´ıch agens, aby mohl nemoc ˇs´ıˇrit dál, a neobjevuj´ı se ˇza´dné pˇr´ıznaky nemoci. Proto vznikly tyto modely obsahuj´ıc´ı dalˇs´ı tˇr´ıdu E jako mezistupeˇ n mezi náchyln´ ymi a infekˇcn´ımi jedinci. Rychlost, s jakou nemoc naplno propukne a infikovan´ y jedinec se stane infekˇcn´ım, vyjadˇruje konstanta σ, pr˚ umˇerná latentn´ı doba je 1/σ. dS dt dE dt dI dt dR dt

= −βSI,

(1.8.1)

= βSI − σE,

(1.8.2)

= σE − γI,

(1.8.3)

= γI.

(1.8.4)

U nˇekter´ ych nemoc´ı se vyskytuje bˇehem latentn´ıho stadia urˇcitá n´ızká m´ıra infekˇcnosti. To m˚ uˇze b´ yt modelováno pˇredpokládán´ım infekˇcnosti sn´ıˇzené 33

vynásoben´ım koeficientem ε bˇehem stadia latence: dS dt dE dt dI dt dR dt

= −βS(I + εE),

(1.8.5)

= βS(I + εE) − σE,

(1.8.6)

= σE − γI,

(1.8.7)

= γI.

(1.8.8)

Tento model je nˇekdy uvádˇen také jako SLIR model, s tˇr´ıdou L s redukovanou infekˇcnost´ı, kterou nahrad´ıme tˇr´ıdu E. Základn´ı reprodukˇcn´ı ˇc´ıslo R0 je u tohoto modelu rovno

βN βN +ε . γ σ

Obrázek 1.14: Graf pro SEIR model (β = 0.003, γ = 0.3, σ = 0.1).

34

Proti nˇekter´ ym nemocem se u nakaˇzeného jedince nevytváˇr´ı imunita, a proto se po prodˇelán´ı nemoci vrát´ı zpˇet do tˇr´ıdy náchyln´ ych:

dS = −βSI + γI, dt dE = βSI − σE, dt dI = σE − γI. dt

(1.8.9) (1.8.10) (1.8.11)

Dopadem delˇs´ı latentn´ı doby na pr˚ ubˇeh epidemie je pomalejˇs´ı poˇca´teˇcn´ı rychlost epidemie a delˇs´ı trván´ı epidemie, jak si m˚ uˇzeme vˇsimnout, porovnáme-li obrázek 1.14, znázorˇ nuj´ıc´ı graf pro SEIR model, s grafem základn´ıho SIR modelu na obrázku 1.1, kde uvaˇzujeme nulovou latenci. Obdobnˇe m˚ uˇzeme porovnávat i obrázek 1.15, ilustruj´ıc´ı SEIS model, s grafem SIS modelu (viz obrázek 1.11 vlevo dole), kde byly pouˇzity stejné hodnoty konstant β a γ. Celkov´ y poˇcet nakaˇzen´ ych na délce latence nezávis´ı.

Obrázek 1.15: Graf pro SEIS model (β = 0.003, γ = 0.3, σ = 0.1).

35

1.9

SIHR model

Tento model popisuje situaci, kdy se nemocn´ı sami od sebe pˇresouvaj´ı do izolace v podobˇe tˇr´ıdy H a sniˇzuje se tak kontakt mezi náchyln´ ymi a infekˇcn´ımi jedinci a t´ım i ˇs´ıˇren´ı nákazy. dS dt dI dt dH dt dR dt

= −βSI,

(1.9.1)

= βSI − γI − θI,

(1.9.2)

= θI − δH,

(1.9.3)

= γI + δH.

(1.9.4)

Konstanta θ urˇcuje rychlost, s jakou se infekˇcn´ı jedinci po nakaˇzen´ı pˇresunou do domác´ı izolace, konstanta δ udává rychlost pˇrechodu jedinc˚ u v domác´ı izolaci mezi imunn´ı. Izolace zp˚ usob´ı poˇca´teˇcn´ı pomalejˇs´ı nár˚ ust epidemie, niˇzˇs´ı poˇcet onemocnˇen´ı, ale také delˇs´ı dobu trván´ı epidemie. To m˚ uˇzeme pozorovat, srovnáme-li obrázek 1.16 s grafem SIHR modelu s obrázkem 1.1, na kterém je model základn´ı.

Obrázek 1.16: Graf pro SIHR model (β = 0.003, γ = 0.3, θ = 0.5, δ = 0.2).

36

1.10

SITR model

Pokud existuje léˇcba pro infekˇcn´ı osoby, m˚ uˇzeme pˇredpokládat, ˇze ˇcást infekˇcn´ıch jedinc˚ u α za jednotku ˇcasu je vybrána pro léˇcen´ı a léˇcba sn´ıˇz´ı jejich infekˇcnost o pod´ıl δ. M´ıra pˇresunu z tˇr´ıdy léˇcen´ ych je κ. To vede k SITR modelu (viz obrázek 1.17), kde jedinec m˚ uˇze b´ yt bud’ náchyln´ y → infekˇcn´ı → imunn´ı, nebo náchyln´ y → infekˇcn´ı → léˇcen´ y → imunn´ı: dS dt dI dt dT dt dR dt

= −βS(I + δT ),

(1.10.1)

= βS(I + δT ) − (γ + α)I,

(1.10.2)

= αI − κT,

(1.10.3)

= γI + κT.

(1.10.4)

Obrázek 1.17: Graf pro SITR model (β = 0.003, γ = 0.3, α = 0.5, δ = 0.75, κ = 0.5).

37

1.11

SIRD model

SIRD model se vyuˇz´ıvá pˇri modelován´ı nemoc´ı, které mohou b´ yt smrtelné. Zavedeme tˇr´ıdu zemˇrel´ ych jedinc˚ u D a konstantu δ, která udává m´ıru u ´mrt´ı nakaˇzen´ ych jedinc˚ u. T´ım z´ıskáme soustavu dS dt dI dt dR dt dD dt

= −βSI,

(1.11.1)

= βSI − γI,

(1.11.2)

= γI − δI,

(1.11.3)

= δI.

(1.11.4)

Vlastn´ı pr˚ ubˇeh epidemie se oproti SIR modelu nemˇen´ı, pouze se kv˚ uli u ´mrt´ım sniˇzuje poˇcet imunn´ıch jedinc˚ u v populaci a stejnˇe tak klesá celková velikost populace N (viz obrázek 1.18).

Obrázek 1.18: Graf pro SIRD model (β = 0.003, γ = 0.3, δ = 0.1). SIRD model lze modifikovat podobnˇe jako SIR model, napˇr´ıklad zahrnut´ım vakcinace nebo karantény, zaj´ımavé jsou ovˇsem pˇredevˇs´ım u ´pravy, které zasahuj´ı do tˇr´ıdy D. Kupˇr´ıkladu léˇcba infikovan´ ych sniˇzuje u ´mrtnost o konstantu ζ. Rovnice (1.11.1) a (1.11.2) z˚ ustávaj´ı stejné, rovnice (1.11.3) 38

a (1.11.4) se zmˇen´ı na dR = γI − (δ − ζ)I, dt dD = (δ − ζ)I. dt

(1.11.5) (1.11.6)

´ Umrt´ ı m˚ uˇzeme do modelu zahrnout i bez pˇridán´ı zvláˇstn´ı tˇr´ıdy D a vyjádˇrit je pomoc´ı soustavy rovnic (1.4.7), (1.4.8) a (1.4.9). V této soustavˇe vyuˇz´ıváme parametr f s t´ım, ˇze pˇredpokládáme m´ıru pˇreˇzit´ı nemoci pˇrinejmenˇs´ım rovnu f. V [5, str. 688] lze naj´ıt následuj´ıc´ı tvrzen´ı: Vˇ eta 1. Plat´ı, ˇze limt→∞ S = S∞ pro epidemii s reprodukˇcn´ım ˇc´ıslem R0 a s m´ırou pˇreˇzit´ı pˇrinejmenˇs´ım f nen´ı menˇs´ı neˇz hodnota S∞ pro epidemii bez u ´mrt´ı s reprodukˇcn´ım ˇc´ıslem R0 /f . D´ıky této vˇetˇe jsme schopni z´ıskat velmi dobrou aproximaci S∞ pro model s u ´mrt´ımi. Velikost epidemie pro smrtelnou nemoc tedy nen´ı stanovena pˇresnˇe, ale horn´ı hranici velikosti epidemie udává koneˇcná velikost epidemie bez u ´mrt´ı s vyˇsˇs´ım reprodukˇcn´ım ˇc´ıslem, coˇz je pro modelován´ı dostaˇcuj´ıc´ı.

1.12

Shrnut´ı pˇ rehledu model˚ u

Z hlediska pˇrenosov´ ych mechanism˚ u zahrnuj´ı modely v epidemiologii ˇradu rozliˇcn´ ych faktor˚ u. Jsme schopni zohlednit kontaktn´ı pˇrenáˇsen´ı onemocnˇen´ı (z ˇclovˇeka na ˇclovˇeka), vertikáln´ı pˇrenáˇsen´ı onemocnˇen´ı (z rodiˇc˚ u na potomky) a vektorové pˇrenáˇsen´ı onemocnˇen´ı (z mezihostitele, napˇr. z krys, komár˚ u nebo kl´ıˇst’at, na ˇclovˇeka), zaˇclenit do modelu inkubaci, latentn´ı periodu onemocnˇen´ı, u ´mrtnost, sezónnost, izolaci, karanténu, vakcinaci, vakcinaci se ztrátou imunity, léˇcbu, nákazu v rámci skupiny nebo mezi v´ıcero skupinami i u skupin s odliˇsnou populaˇcn´ı dynamikou. Sofistikovanˇejˇs´ı modely zahrnuj´ı vˇek, pohlav´ı, prostorovou strukturu rozloˇzen´ı populace, cestován´ı, psychologické struktury, odliˇsné pravdˇepodobnosti nakaˇzen´ı pro r˚ uzné vˇekové kategorie a zjednoduˇsené demografické jevy (natalita, mortalita, migrace). Koeficienty v modelech mohou b´ yt konstantn´ı nebo závislé na ˇcase. Pˇredchoz´ı v´ yˇcet, kter´ y byl obsahem této kapitoly, tedy pˇredstavuje pouze omezenou ˇcást moˇznost´ı, zvláˇstˇe vzhledem k tomu, ˇze jsme schopni v´ yˇse uvedené alternativy rozˇs´ıˇren´ı model˚ u navzájem kombinovat a vytváˇret tak detailnˇejˇs´ı a sloˇzitˇejˇs´ı modely.

39

Vˇetˇsina deterministick´ ych model˚ u je zaloˇzena na soustavˇe obyˇcejn´ ych diferenciáln´ıch rovnic, ale pro komplikovanˇejˇs´ı modely jsou vyuˇz´ıvány napˇr. parciáln´ı diferenciáln´ı rovnice prvn´ıho nebo druhého ˇra´du, diferenciáln´ı rovnice se zpoˇzdˇen´ım a diferenciáln´ı rovnice s impulsy. Nˇekdy jsou vyuˇz´ıvány i modely diskrétn´ı, pˇredevˇs´ım pro stavy pˇred vypuknut´ım epidemie, pro modelován´ı podle sesb´ıran´ ych statistick´ ych dat nebo u model˚ u s vˇekovou strukturou. Jsou tvoˇreny soustavou diferenˇcn´ıch rovnic, vzhledem k povaze statistick´ ych u ´daj˚ u je ˇcasov´ y krok obvykle volen jako jeden den. Pˇri studiu deterministick´ ych model˚ u nás obvykle zaj´ımá pˇredevˇs´ım korektnost model˚ u a jejich ˇreˇsen´ı, pˇretrváván´ı nemoc´ı, existence a stabilita stacionárn´ıch ˇreˇsen´ı, které charakterizuj´ı ˇs´ıˇren´ı nemoc´ı a jejich endemickost, existence a stabilita periodick´ ych ˇreˇsen´ı, které popisuj´ı oscilaci pˇrenosu nemoci, a v´ yskyt bifurkac´ı a chaotického chován´ı.

40

Kapitola 2 Anal´ yza z´ akladn´ıho SIR modelu 2.1

Pˇ reveden´ı modelu na rovnici vyˇ sˇ s´ıho ˇ r´ adu

Vˇeta v této kapitole byla vyslovena v [3], z této publikace vycház´ı i idea d˚ ukazu, kter´ y je zde podrobnˇeji rozpracován. Vˇ eta 2. Systém rovnic (1.2.1), (1.2.2) a (1.2.3) je ekvivalentn´ı s nelineárn´ı rovnic´ı druhého ˇrádu

β

R00 = c βR0 e− γ R −γR0 , β

(2.1.1) β

kde c = S0 e γ R0 . Pro neznámé S a I plat´ı S = c e− γ R , I = N − S − R. D˚ ukaz. Nejprve zderivujeme (1.2.1) podle ˇcasu a do z´ıskané rovnice S 00 = −βS 0 I − βSI 0

(2.1.2)

dosad´ıme za I pod´ıl z (1.2.1), tj. I=−

S0 , βS

(2.1.3)

odkud dostaneme (S 0 )2 S = − βSI 0 . S 00

41

(2.1.4)

Rovnici (2.1.4) vydˇel´ıme souˇcinem −βS a po u ´pravˇe z´ıskáme vztah " 0 2 # 00 1 S S I0 = − − . (2.1.5) β S S Nyn´ı uprav´ıme rovnici (1.2.2). Za levou stranu dosad´ıme vztah (2.1.5) a za I dosad´ıme (2.1.3), vynásob´ıme −β a z´ıskáme 0 2 S0 S 00 S − = S 0β − γ . (2.1.6) S S S Do rovnice (1.2.3) dosad´ıme za I ze vztahu (2.1.3), tj. γ S0 R =− , βS 0

(2.1.7)

ˇc´ımˇz máme z (1.2.1) i z (1.2.3) eliminovanou promˇennou I. Vztah (2.1.7) zintegrujeme podle ˇcasu na γ R + c1 = − (ln |S|) β

(2.1.8)

a uprav´ıme do podoby β

S = c e− γ R ,

(2.1.9)

kde c je kladná integraˇcn´ı konstanta, kterou urˇc´ıme z podm´ınek R(0) = R0 , S(0) = S0 , tj. β

c = S0 e γ R0 .

(2.1.10)

Rovnici (2.1.9) pak zderivujeme podle ˇcasu na S 0 = −c

β 0 − βγ R . R e γ

(2.1.11)

Podruhé zderivujeme rovnici (2.1.7) podle ˇcasu a do z´ıskaného vztahu " 0 2 # 00 γ S S R00 = − (2.1.12) − β S S dosad´ıme z rovnice (2.1.6) a roznásob´ıme, tj. R00 = −γS 0 + 42

γ2 S0 . β S

(2.1.13)

γ2 S0 nahrad´ıme na základˇe β S vztahu (2.1.7) souˇcinem −γR0 a z´ıskáme nelineárn´ı rovnici druhého ˇrádu Ted’ uˇz jen dosad´ıme za S 0 z (2.1.11) a souˇcin β

R00 = c βR0 e− γ R −γR0 .

(2.1.14)

Modelován´ı epidemie pomoc´ı jedné nelineárn´ı ODR 2. ˇra´du zkracuje ˇcas numerické simulace oproti bˇeˇznˇe uˇz´ıvané soustavˇe ODR (viz str. 62). Zároveˇ n toto vyjádˇren´ı pomáhá dokázat vˇetu 3 v následuj´ıc´ı sekci.

2.2

Pˇ reveden´ı modelu na soustavu rovnic s parametrem

Vˇeta v této kapitole byla opˇet vyslovena v [3], z této publikace vycház´ı i idea d˚ ukazu, kter´ y je zde podrobnˇeji rozpracován. ˇ sen´ı SIR soustavy lze vyjádˇrit parametricky pomoc´ı parametru u Vˇ eta 3. Reˇ ve tvaru: Z u dξ , (2.2.1) t = u0 ξ(−βN − γ ln(ξ) + cβξ) S = cu, γ I = ln(u) − cu + N, β γ R = − ln(u), β

(2.2.2) (2.2.3) (2.2.4)

β R0 kde c = S0 e γ . Poznámka. Plat´ı t = 0 pro u = u0 a t > 0 pro u < u0 . D˚ ukaz. Zavedeme funkci u(t) pˇredpisem β

u = e− γ R .

(2.2.5)

Dále si vyjádˇr´ıme prvn´ı a druhou derivaci u podle ˇcasu: β u0 = − R0 u, γ β β u00 = − R00 u − u0 R0 . γ γ 43

(2.2.6) (2.2.7)

Plat´ı R

0

R00

γ u0 , = − βu β 0 0 00 γ u + γR u = − . β u

(2.2.8) (2.2.9)

Nyn´ı vyuˇzijeme vˇetu 2 a dosazen´ım do (2.1.1) z´ıskáme 0

β 0γ u 00 γ u0 γ u − γu β u − (γ − cβu) = 0, − β u βu

(2.2.10)

po u ´pravách a vynásoben´ım u2 dostaneme rovnici uu00 − (u0 )2 + uu0 (γ − cβu) = 0.

(2.2.11)

Nyn´ı si zavedeme novou funkci φ pˇredpisem φ=

dt du

(2.2.12)

a odtud u0 =

1 , φ

u00 = −

(2.2.13)

1 dφ 1 dφ =− 3 . 2 φ dt φ du

(2.2.14)

Dosazen´ım do (2.2.11) z´ıskáme rovnici 1 dφ −u 3 − φ du

2 1 1 + u (γ − cβu) = 0 φ φ

(2.2.15)

a po roznásoben´ı φ3 , vydˇelen´ım −u a následné u ´pravˇe dostaneme diferenciáln´ı rovnici Bernoulliho typu dφ 1 + φ = (γ − cβu)φ2 . du u

(2.2.16)

Obecn´ ym ˇreˇsen´ım této rovnice (viz pˇr´ıloha A.2) je φ(u) =

1 , u(C − γ ln u + cβu) 44

(2.2.17)

kde C je libovolná integraˇcn´ı konstanta. Pokud toto ˇreˇsen´ı zintegrujeme podle u, z´ıskáme Z Z 1 dt du = du, (2.2.18) du u(C − γ ln u + cβu) resp. pokud pouˇzijeme urˇcité integrály a pˇreznaˇc´ıme integraˇcn´ı promˇenné, dostáváme Z u Z t 1 1 dτ = dξ, (2.2.19) u0 ξ(C − γ ln(ξ) + cβξ) t0 Z u dξ t − t0 = , (2.2.20) u0 ξ(C − γ ln(ξ) + cβξ) za t0 m˚ uˇzeme zvolit nulu bez u ´jmy na obecnosti. Z (2.1.9) vyjádˇr´ıme S = cu

(2.2.21)

a zlogaritmován´ım (2.2.5) z´ıskáme γ R = − ln(u). β

(2.2.22)

Pˇredpis pro I z´ıskáme následuj´ıc´ım zp˚ usobem: nejprve vynásob´ıme rovnici (2.2.17) souˇcinem βu a uprav´ıme, tedy β βu = , u0 C − γ ln u + cβu

(2.2.23)

dále dosazen´ım ze vztahu (2.2.8) dostaneme rovnici β β = , C − γ ln u + cβu R0 γ

(2.2.24)

γ β = , R0 C − γ ln u + cβu

(2.2.25)

−β po u ´pravˇe −

vyuˇzijeme vztahu (2.1.7) a uprav´ıme levou stranu rovnice, tj. −

βS β = − , 0 S C − γ ln u + cβu 45

(2.2.26)

a následnˇe s vyuˇzit´ım vztahu (2.1.3) a po u ´pravˇe z´ıskáme C γ ln(u) − cu − . β β

I =

(2.2.27)

Souˇctem (2.2.21), (2.2.22) a (2.2.27) z´ıskáváme C S+I +R = − , β

(2.2.28)

odtud C = −βN.

(2.2.29)

Z´ıskali jsme parametricky vyjádˇrené neznámé funkce S, I a R a ˇcas t. ˇ C´ıselnou hodnotu u v ˇcase nula urˇcuj´ı poˇca´teˇcn´ı podm´ınky a parametry modelu, respektive β

u0 = e− γ R0 .

(2.2.30)

V´ yhodou parametrického vyjádˇren´ı ˇreˇsen´ı je skuteˇcnost, ˇze udává analytické ˇreˇsen´ı popisuj´ıc´ı dynamick´ y v´ yvoj SIR systému pro dané poˇca´teˇcn´ı podm´ınky S0 , I0 a R0 a libovolné hodnoty konstant β a γ. Pˇresnost ˇreˇsen´ı ovlivˇ nuje pouze integrál, kter´ y nemus´ıme b´ yt schopn´ı spoˇc´ıtat analyticky.

2.2.1

Vrchol epidemie

Funkce I(t) je diferencovatelná a v ˇcase, ve kterém dosahuje svého maxima (dále budeme znaˇcit tmax ), mus´ı m´ıt nulovou prvn´ı derivaci. Plat´ı 0 = I 0 (tmax ) = βS(tmax )I(tmax ) − γI(tmax ),

(2.2.31)

po u ´pravách z´ıskáme vztah S(tmax ) = γ/β.

(2.2.32)

Z parametrického vyjádˇren´ı SIR modelu m˚ uˇzeme jednoduˇse vyjádˇrit maximáln´ı poˇcet nemocn´ ych i ˇcas, ve kterém tento stav nastane. V´ıme, ˇze nejv´ıce pacient˚ u je pro S(tmax ) = γ/β. Z (2.2.2) pak z´ıskáme hodnotu parametru u, se kterou urˇc´ıme maximum nemocn´ ych jako umax =

γ1 . βc

46

(2.2.33)

Dosazen´ım (2.2.33) do vztahu (2.2.3) z´ıskáme nejvyˇsˇs´ı poˇcet infekˇcn´ıch jedinc˚ u jako γ γ γ Imax = I(tmax ) = ln − + N. (2.2.34) β βc β Toho je dosaˇzeno v ˇcase Z

umax

tmax = u0

dξ . ξ(−βN − γ ln(ξ) + cβξ)

(2.2.35)

Pˇripomeˇ nme, ˇze β

β

u0 = e− γ R0 ,

c = S 0 e γ R0 ,

umax =

γ β

βS0 e γ R0

,

(2.2.36)

a tedy Imax

γ ln = β

γ γ γ − ln (S0 ) − R0 − + N, β β β

 γ β R βS0 e γ 0 β − R e γ 0

Z tmax =



 

dξ

β

ξ −βN − γ ln(ξ) + S0 e γ R0 βξ

.

Hodnotu tohoto integrálu jsme schopni pˇribliˇznˇe stanovit za pouˇzit´ı numerick´ ych ˇreˇsiˇc˚ u. Poˇcet infekˇcn´ıch jedinc˚ u na vrcholu epidemie klesá s rostouc´ım poˇctem imunn´ıch jedinc˚ u na poˇcátku epidemie (obrázek 2.1), analyticky: dImax d γ γ γ γ = ln − ln (N − R0 − I0 ) − R0 − + N dR0 dR0 β β β β γ 1 γ = − − −1= − 1. β N − R0 − I0 βS0 Protoˇze mus´ı platit, ˇze S0 > v závislosti na R0 klesaj´ıc´ı.

γ dImax (viz str. 20), < 0 a funkce Imax je β dR0

Hodnota Imax klesá i pˇri zvˇetˇsuj´ıc´ım se pod´ılu parametr˚ u γ a β (obrázek γ 2.2), jak m˚ uˇzeme opˇet analyticky dokázat. Oznaˇcme =: δ, β 47

Obrázek 2.1: Poˇcet infekˇcn´ıch jedinc˚ u na vrcholu epidemie Imax v závislosti na R0 , N = 1020, I0 = 20, γ = 0.3, β = 0.003.

dImax d = [δ ln (δ) − δ ln (S0 ) − R0 − δ + N ] dδ dδ δ = ln(δ) + − ln(S0 ) − 1 δ γ = ln (δ) − ln(S0 ) = ln − ln(S0 ). β γ γ Jak jiˇz bylo zm´ınˇeno, S0 > , a tedy ln(S0 ) > ln . Proto je vˇzdy β β γ dImax < 0 a funkce Imax je v závislosti na klesaj´ıc´ı. dδ β γ Závislost Imax na hodnotách R0 i souhrnnˇe zachycuje obrázek 2.3. β ˇ Cas, ve kterém epidemie dosáhne svého vrcholu a infekˇcn´ıch jedinc˚ u je nejv´ıce, rovnˇeˇz závis´ı na zastoupen´ı imunn´ıch jedinc˚ u v populaci. Obecnˇe v´ıce imunn´ıch jedinc˚ u na zaˇcátku epidemie znamená prodlouˇzen´ı epidemie a pozdˇejˇs´ı dosaˇzen´ı jej´ıho vrcholu, záleˇz´ı ovˇsem i na poˇctu jedinc˚ u v ostatn´ıch tˇr´ıdách v ˇcase 0 a také na parametrech γ a β. Podrobnˇejˇs´ı anal´ yza funkce tmax (R0 ) by tedy byla velice sloˇzitá, pˇribliˇznou pˇredstavu o jej´ım chován´ı si m˚ uˇzeme udˇelat na základˇe numerick´ ych experiment˚ u na obrázc´ıch 2.4 a 2.5.

48

Obrázek 2.2: Poˇcet infekˇcn´ıch jedinc˚ u na vrcholu epidemie Imax v závislosti γ na pod´ılu , S0 = 800, I0 = 20, R0 = 200. β

Obrázek 2.3: Zmˇena poˇctu infekˇcn´ıch jedinc˚ u na vrcholu epidemie Imax γ v závislosti na pod´ılu a na R0 , N = 1020, I0 = 20. β

49

ˇ tmax dosaˇzen´ı maximáln´ıho poˇctu infekˇcn´ıch jedinc˚ Obrázek 2.4: Cas u Imax v závislosti na R0 , N = 1000, I0 = 1, zleva doprava γ = 0.3, 0.3, 0.1, 0.9, β = 0.001, 0.01, 0.003, 0.003.

50

ˇ tmax dosaˇzen´ı maximáln´ıho poˇctu infekˇcn´ıch jedinc˚ Obrázek 2.5: Cas u Imax v závislosti na R0 , N = 1000, γ = 0.3, β = 0.003, zleva doprava I0 = 1, 10, 50, 100, 400, 600.

51

2.3

Pˇ reveden´ı modelu na rovnici Abelova typu se separovan´ ymi promˇ enn´ ymi

Vˇ eta 4. Systém rovnic (1.2.1), (1.2.2) a (1.2.3) je ekvivalentn´ı s rovnic´ı Abelova typu β dv = v 2 + βv 3 , (2.3.1) dψ γ 0 kde v = RR00 . Pro neznámé funkce v závislosti na v plat´ı: S = β1 v1 + γ , h i R = − βγ ln γ + v1 − ln (βS0 ) − βγ R0 , I = N − S − R. D˚ ukaz. Z rovnice (1.2.3) vyjádˇr´ıme I a vztah následnˇe zderivujeme podle ˇcasu, tedy R0 , γ R00 . = γ

I =

(2.3.2)

I0

(2.3.3)

Do rovnice (1.2.1) dosad´ıme za I vztah (2.3.2), do rovnice (1.2.2) dosad´ıme za I a I 0 vztahy (2.3.2) a (2.3.3): β S 0 = − SR0 , γ R00 β = SR0 − R0 . γ γ

(2.3.4) (2.3.5)

Z (2.3.5) vyjádˇr´ıme souˇcin βS =

R00 +γ R0

(2.3.6)

a zderivujeme podle ˇcasu na R000 βS = 0 − R 0

R00 R0

2 .

(2.3.7)

Rovnici (2.3.4) vynásob´ıme konstantou β, do pravé strany dosad´ıme ze vztahu (2.3.6) a levou stranu nahrad´ıme vztahem (2.3.7), z´ıskáme tedy rovnici R000 − R0

R00 R0

2

=

R00 +γ R0 52

β 0 − R . γ

(2.3.8)

Pro zjednoduˇsen´ı zavedeme substituci R0 = ψ,

R00 = ψ 0 ,

R000 = ψ 00

(2.3.9)

a dosad´ıme do (2.3.8), vynásob´ıme ψ 2 a z´ıskáme následuj´ıc´ı diferenciáln´ı rovnici druhého ˇra´du (ψ 0 )2 − ψψ 00 =

β 2 0 ψ ψ + βψ 3 . γ

(2.3.10)

S pomoc´ı transformace w= ψ0 =

1 , w

dt 1 = 0, dψ ψ ψ 00 = −

(2.3.11)

1 dw , w3 dψ

(2.3.12)

se (2.3.10) zmˇen´ı na 2 1 dw 1 β 1 +ψ 3 = ψ 2 + βψ 3 w w dψ γ w

(2.3.13)

a po vynásoben´ı s w3 a u ´pravˇe z´ıskáme ψ

dw β + w = ψ 2 w2 + βψ 3 w3 . dψ γ

(2.3.14)

Zavedeme si novou funkci v pˇredpisem v = wψ =

ψ , ψ0

(2.3.15)

kterou dosad´ıme do rovnice (2.3.14) a uprav´ıme na ψ

β dw + w = v 2 + βv 3 . dψ γ

(2.3.16)

Vyuˇzijeme-li rovnost dw d ψ +w =ψ dψ dψ

1 ψ0

1 d + 0 = ψ dψ

1 ψ 0 ψ

,

(2.3.17)

z´ıskáme diferenciáln´ı rovnici Abelova typu dv β = v 2 + βv 3 . dψ γ 53

(2.3.18)

´ Upravou rovnice (2.3.6) z´ıskáme pˇredpis pro S : 1 R00 1 1 S= +γ . +γ = β R0 β v

(2.3.19)

Pˇredpis pro R z´ıskáme z (2.3.15) s vyuˇzit´ım vztahu (2.1.1) z vˇety 2: R0 R0 v = 00 = , β β R S0 e γ R0 βR0 e− γ R −γR0 1 , v = β β S0 β e γ R0 e− γ R −γ β β 1 + γ = S0 β e γ R0 e− γ R . v Po u ´pravách rovnici zlogaritmujeme, tedy ! 1 + γ β v ln = − R, β γ S0 β e γ R0

(2.3.20) (2.3.21) (2.3.22)

(2.3.23)

a vyjádˇr´ıme R: γ 1 β R=− ln γ + − ln (βS0 ) − R0 . β v γ

(2.3.24)

Také soustavu rovnic (1.3.1), (1.3.2) a (1.3.3), tedy soustavu pro SIR model s vitáln´ı dynamikou, m˚ uˇzeme pˇrevést na Abelovu rovnici. Tato transformace je provedena v [3, str.190]: b µ dv = (a + )v 3 + (c + )v 2 , dψ ψ ψ µ kde a = β( +1), γ

b = µ(µ+γ−βN ),

54

c=

β , γ

ψ = R0 + µR,

(2.3.25) v =

ψ . ψ0

ˇ sen´ı rovnice Abelova typu Reˇ Rovnici (2.3.1) m˚ uˇzeme ˇreˇsit metodou separace promˇenn´ ych. Je-li pravá strana rovnice dv β 2 = v + βv 3 dψ γ rovna nule, z´ıskáváme konstantn´ı ˇreˇsen´ı 1 v=− , γ

v = 0.

Pokud je ale r˚ uzná od nuly, m˚ uˇzeme provést následuj´ıc´ı u ´pravu:

β 2 v γ

Vztah zintegrujeme, tedy Z

dv = dψ. + βv 3

dv = β 2 v + βv 3 γ

(2.3.26)

Z 1dψ,

(2.3.27)

a pouˇzijeme rozklad na parciáln´ı zlomky: Z Z γ2 γ γ3 − + dv = 1dψ. βγv + β βv βv 2

(2.3.28)

Nyn´ı uˇz snadno z´ıskáme ˇreˇsen´ı −

γ2 γ2 γ ln(v) + ln(γv + 1) − = ψ + c. β β βv

Ze vztahu (2.3.21) v´ıme, ˇze 1

v0 := v(t0 ) = S0 β e

β R γ 0

−β R γ 0

e

= −γ

a ze vztah˚ u (1.2.3) a (2.3.9)

ψ0 := ψ(t0 ) = R0 (0) = γI0 . 55

1 S0 β − γ

(2.3.29)

Odtud z´ıskáme hodnotu integraˇcn´ı konstanty c jako γ2 1 γ2 γ γS0 β − γ 2 c = −γI0 − ln + ln +1 − , β S0 β − γ β S0 β − γ β po u ´pravˇe γ2 c= β

βI0 βS0 ln(S0 β) + 1 − − . γ γ

(2.3.30)

Nyn´ı proved’me následuj´ıc´ı u ´pravy: dR = ψ, dt dR(v) dv = ψ(v), dv dt

(2.3.31) (2.3.32)

po dosazen´ı za R a ψ β d γ 1 R 0 − ln γ + − ln βS0 e γ dv = dv β v 2 γ2 γ γ ln(γv + 1) − − c dt. − ln(v) + β β βv Derivován´ım a dalˇs´ımi u ´pravami z´ıskáváme 2 γ γ γ2 γ 1 dv = − ln(v) + ln(γv + 1) − − c dt, βv 2 γ + v1 β β βv γ βv 2 (γ+ v1 ) γ2 β

ln(γv + 1) −

γ2 β

ln(v) −

γ βv

−c

dv = 1dt.

(2.3.33) (2.3.34)

Pokud vztah (2.3.34) zintegrujeme a pˇreznaˇc´ıme integraˇcn´ı promˇenné, dostáváme γ βξ 2 (γ+ 1ξ )

v

Z

v0

Z

v

v0

Z

t

1dτ,

(2.3.35)

γ dξ = t − t0 , (γξ + 1)(γ 2 ξ ln(γ + 1ξ ) − γ − cβξ)

(2.3.36)

γ2 β

ln(γξ + 1) −

γ2 β

ln(ξ) −

γ βξ

−c

dξ = t0

kde za t0 m˚ uˇzeme volit nulu bez u ´jmy na obecnosti a v0 a c jsou známy.

56

Z´ıskali jsme tedy vyjádˇren´ı neznám´ ych funkc´ı S, I a R a ˇcasu t v závislosti na parametru v : Z v γ dξ, t = 1 2 v0 (γξ + 1)(γ ξ ln(γ + ξ ) − γ − cβξ) 1 1 S = S = +γ , β v I = N − S − R, γ 1 β R = − ln γ + − ln (βS0 ) − R0 . β v γ Toto parametrické vyjádˇren´ı ˇreˇsen´ı SIR soustavy je obdobné jako ve vˇetˇe 3. M˚ uˇzeme ovˇsem oˇcekávat v´ yraznˇe vyˇsˇs´ı v´ ypoˇcetn´ı sloˇzitost, proto nebude vyuˇzito v následuj´ıc´ı kapitole srovnávaj´ıc´ı numerické v´ ysledky pro r˚ uzné systémy.

2.4

Porovn´ an´ı numerick´ ych v´ ysledk˚ u jednotliv´ ych syst´ em˚ u pro r˚ uzn´ y software

Tato kapitola se vˇenuje porovnán´ı r˚ uzn´ ych zp˚ usob˚ u vyjádˇren´ı SIR modelu a shrnut´ı jejich v´ yhod a negativ. Dále srovnává vhodnost vyuˇzit´ı program˚ u Matlab a Mathematica z hlediska kvality a ˇcasové nároˇcnosti v´ ypoˇct˚ u. Oba programy byly spuˇstˇeny na stejném poˇc´ıtaˇci a za stejn´ ych podm´ınek pro dan´ y systém pˇetkrát, ˇcasová nároˇcnost uvedená v tabulce 2.1 vznikla jako pr˚ umˇer tˇechto pˇeti hodnot. Ve vˇsech pˇr´ıpadech byly zvoleny následuj´ıc´ı hodnoty parametr˚ u a poˇct˚ u obyvatel: β = 0.003, γ = 0.3, S0 = 930, I0 = 20, R0 = 50. Pro porovnán´ı v´ ysledk˚ u byly stanoveny absolutn´ı a relativn´ı chyby v´ ypoˇct˚ u. Absolutn´ı chyba pˇredstavuje absolutn´ı hodnotu rozd´ılu pˇresného a pˇribliˇzného ˇreˇsen´ı, relativn´ı chyba je pak pomˇer absolutn´ı chyby k absolutn´ı hodnotˇe pˇresného ˇreˇsen´ı. Za pˇresné ˇreˇsen´ı byly povaˇzovány hodnoty neznám´ ych funkc´ı vypoˇcten´ ych ze základn´ıho SIR modelu pˇredstavovaného systémem ODR v daném softwaru. Pouˇzité verze program˚ u: • Wolfram Mathematica 8 • MathWorks MATLAB R2011b 57

K numerickému ˇreˇsen´ı diferenciáln´ıch rovnic v Matlabu byla vyuˇzita funkce ode45. Ta je zaloˇzena na explicitn´ı Runge-Kuttovˇe metodˇe 4. a 5. ˇrádu, na Dormand-Princeovˇe metodˇe. Jedná se o jednokrokov´ y ˇreˇsiˇc. V Mathematice byla vyuˇzita funkce NDSolve, která sama vol´ı vhodnou metodu v závislosti na typu rovnice, pokud ji sami neurˇc´ıme.

Model jako soustava diferenci´ aln´ıch rovnic Pˇrestoˇze z SIR modelu popsaného systémem diferenciáln´ıch rovnic nelze snadno vyˇc´ıst nˇekteré hodnoty, napˇr. ˇcas dosaˇzen´ı vrcholu epidemie, je jeho numerické ˇreˇsen´ı vypoˇcteno velmi rychle. Oba programy dosáhly obdobn´ ych v´ ysledk˚ u (viz obrázek 2.6, pro stanoven´ı velikosti relativn´ı chyby byly za pˇresné povaˇzovány v´ ysledky z Matlabu), v pˇr´ıpadˇe Matlabu probˇehl v´ ypoˇcet o nˇeco rychleji neˇz v Mathematice (viz tabulka 2.1).

ˇ sen´ı modelu v parametrick´ Reˇ e formˇ e Pˇrestoˇze pˇreveden´ı systému diferenciáln´ıch rovnic na ˇreˇsen´ı modelu v parametrické formˇe pˇrináˇs´ı mnoho v´ yhod, pˇredevˇs´ım jednoduˇsˇs´ı vyjádˇren´ı nˇekter´ ych mezn´ıch hodnot, je numerická simulace v´ ysledk˚ u v tomto pˇr´ıpadˇe v´ yraznˇe ˇcasovˇe nároˇcnˇejˇs´ı. Parametrizace je neuniformn´ı. Parametr pro numerické simulace je rozumné volit od hodnoty, pˇri které se poˇcet infekˇcn´ıch jedinc˚ u bl´ıˇz´ı nule (v ukázkovém pˇr´ıkladˇe byla zvolena hodnota 0.000046), do hodnoty u0 . V pˇr´ıpadˇe Matlabu, ˇc´ım v´ıc se vzdalujeme od ˇcasu t0 , t´ım potˇrebujeme jemnˇejˇs´ı krok, abychom z´ıskali dostateˇcn´ y poˇcet bod˚ u pro dobrou pˇredstavu o pr˚ ubˇehu epidemie, viz obrázek 2.7 ilustruj´ıc´ı rozloˇzen´ı vykreslovan´ ych bod˚ u pˇri konstantn´ım kroku. Pokud bychom zvolili velmi mal´ y pevn´ y krok pro cel´ y v´ ypoˇcet, mnoˇzstv´ı z´ıskan´ ych bod˚ u by bylo pro naˇsi pˇredstavu o pr˚ ubˇehu epidemie dostateˇcnˇe velké, ale ˇcasová nároˇcnost v´ ypoˇct˚ u by se na bˇeˇzném poˇc´ıtaˇci v´ yraznˇe zvˇetˇsila aˇz do ˇra´du dn˚ u. Vyˇsˇs´ı hustoty spoˇcten´ ych bod˚ u v grafu na obrázku 2.8 a relativnˇe hezkého ˇcasu v tabulce 2.1 bylo dosaˇzeno rozdˇelen´ım v´ ypoˇctu do 4 ˇca´st´ı s r˚ uzn´ ymi kroky parametru. Preciznˇejˇs´ım dˇelen´ım by se ˇcasová nároˇcnost mohla dále sniˇzovat. Mathematica sice zvládla vykreslen´ı grafu za zlomek doby oproti Matlabu, u vyˇsˇs´ıch ˇcas˚ u byl ovˇsem proveden nedostateˇcn´ y poˇcet v´ ypoˇct˚ u a skuteˇcn´ y pr˚ ubˇeh funkc´ı byl lineárnˇe aproximován (viz obrázek 2.9). 58

SIR model jako soustava ODR ˇreˇsen´ı SIR modelu v parametrické formˇe SIR model pˇres jednu ODR 2. ˇra´du

Matlab Mathematica 1.25 s 1.60 s 2.5 h 6.20 s 1.05 s 0.36 s

Tabulka 2.1: Porovnán´ı ˇcasové nároˇcnosti v´ ypoˇct˚ u odliˇsn´ ych zp˚ usob˚ u vyjádˇren´ı SIR modelu.

Obrázek 2.6: Porovnán´ı v´ ysledk˚ u z´ıskan´ ych pro SIR model jako soustava ODR v Matlabu a v Mathematice. 59

Obrázek 2.7: SIR model jako soustava rovnic s parametrem s pevn´ ym krokem ui+1 = ui + 0.0001.

Obrázek 2.8: SIR model v Matlabu: soustava ODR 1. ˇra´du, ˇreˇsen´ı v parametrické formˇe a rozd´ılnost dosaˇzen´ ych v´ ysledk˚ u. 60

Obrázek 2.9: SIR model v Mathematice: soustava ODR 1. ˇra´du, ˇreˇsen´ı v parametrické formˇe a rozd´ılnost dosaˇzen´ ych v´ ysledk˚ u.

61

Model pˇ res ODR 2. ˇ r´ adu Pokud vypoˇcteme mnoˇzstv´ı imunn´ıch jedinc˚ u R v ˇcase t z nelineárn´ı obyˇcejné diferenciáln´ı rovnice 2. ˇrádu a na základˇe tˇechto hodnot dopoˇcteme náchylné jedince S a infekˇcn´ı jedince I, z´ıskáme v´ ysledky nejrychleji z uveden´ ych zp˚ usob˚ u (viz tabulka 2.1). Mathematica bude v tomto pˇr´ıpadˇe rychlejˇs´ı neˇz Matlab. Grafy s pr˚ ubˇehem epidemie a rozd´ılnost v´ ysledk˚ u v porovnán´ı se systémem ODR 1. ˇra´du zachycuj´ı obrázky 2.10 a 2.11.

Obrázek 2.10: SIR model pˇres jednu ODR 2. ˇra´du v Matlabu a rozd´ılnost v´ ysledk˚ u v porovnán´ı se soustavou ODR.

62

Obrázek 2.11: SIR model pˇres jednu ODR 2. ˇra´du v Mathematice a rozd´ılnost v´ ysledk˚ u v porovnán´ı se soustavou ODR.

Shrnut´ı numerick´ ych v´ ysledk˚ u

Vypoˇcten´ı hodnot funkc´ı S, I a R dosáhneme nejrychleji pˇri numerické simulaci SIR modelu vyjádˇreného pˇres jednu obyˇcejnou diferenciáln´ı rovnici druhého ˇrádu v programu Mathematica. Pokud bychom chtˇeli vyuˇz´ıt program Matlab, je také rozumné pouˇzit´ı vyjádˇren´ı pˇres jednu ODR 2. ˇra´du, pˇr´ıpadnˇe soustavy tˇr´ı ODR 1. ˇrádu (v tomto pˇr´ıpadˇe z´ıskáme hodnoty rychleji v Matlabu neˇz v Mathematice).

Pokud bychom povaˇzovali v´ ysledky soustavy ODR, která je pro popis pr˚ ubˇehu epidemi´ı vyuˇz´ıvána nejˇcastˇeji, za pˇresné, bude se relativn´ı chyba pˇri vyuˇzit´ı jiného zp˚ usobu popisu epidemie uvedeného v práci pohybovat u Matlabu v ˇrádu tis´ıcin a v Mathematice v ˇrádu miliontin.

63

2.5

Z´ akladn´ı reprodukˇ cn´ı ˇ c´ıslo SIR modelu a celkov´ y rozsah epidemie

Obsah této kapitoly vycház´ı z [1, str. 353]. Obecné informace o základn´ım reprodukˇcn´ım ˇc´ısle R0 byly uvedeny na stranˇe 17, v souvislosti s SIR modelem bylo R0 zm´ınˇeno na stranˇe 20. D˚ uleˇzit´ ym pˇredpokladem je skuteˇcnost, ˇze v ˇcase 0 tvoˇr´ı témˇeˇr celou β β populaci náchyln´ı jedinci, tedy ˇze S0 ≈ N a R0 = S0 ≈ N . γ γ Ze soustavy rovnic (1.2.1), (1.2.2) a (1.2.3) m˚ uˇzeme vynechat posledn´ı rovnici, nebot’ poˇcet imunn´ıch jedinc˚ u je jasnˇe dan´ y známou velikost´ı populace a poˇcty náchyln´ ych a infekˇcn´ıch jedinc˚ u. Z´ıskáme tak soustavu dS = −βSI, dt dI = βSI − γI, dt

(2.5.1) (2.5.2)

s poˇca´teˇcn´ı podm´ınkou S0 + I0 = N . Seˇcten´ım (2.5.1) a (2.5.2) z´ıskáme dS dI + dt dt

= −γI.

(2.5.3)

Souˇcet (S + I ) je nezáporná klesaj´ıc´ı funkce. Oznaˇcme si lim (S + I) = I∞ + S∞ .

t→∞

(2.5.4)

Pˇredpokládejme, ˇze poˇcet infekˇcn´ıch jedinc˚ u klesá k nule a epidemie skonˇc´ı, tj. lim I(t) = I∞ = 0,

t→∞

a limita ve (2.5.4) je tedy rovna S∞ . Integrován´ım (2.5.3) z´ıskáme Z ∞ − (S + I)0 dt = S0 + I0 − S∞ − I∞ = N − S∞ .

(2.5.5)

(2.5.6)

0

Rovnici (2.5.1) vydˇel´ıme funkc´ı S a zintegrujeme: Z ∞ 0 Z ∞ S dt = −β I dt, S 0 0 Z ∞ ln S0 − ln S∞ = β I dt. 0

64

(2.5.7) (2.5.8)

Do pravé strany rovnice dosad´ıme ze vztahu (2.5.3) a následnˇe z (2.5.6), tj. Z ∞ Z β ∞ β I dt = − (S + I)0 dt = [N − S∞ ] . (2.5.9) β γ 0 γ 0 Z´ıskáváme v´ yslednou rovnici S0 β S∞ S∞ ln = N 1− = R0 1 − , S∞ γ N N

(2.5.10)

S0 a S∞ jsou urˇceny sérologickou studi´ı mˇeˇren´ım imunitn´ıch reakc´ı ze vzork˚ u krve ˇclen˚ u populace pˇred a po epidemii. Rovnice (2.5.10) je transcendentn´ı rovnic´ı udávaj´ıc´ı závislost celkového rozsahu epidemie na R0 , v angliˇctinˇe se naz´ yvá final size relation. Celkov´ y rozsah epidemie, nebo-li poˇcet jedinc˚ u infikovan´ ych v pr˚ ubˇehu epidemie, je roven N − S∞ . Hodnota (1 − S∞ /N ) se oznaˇcuje jako m´ıra napaden´ı.

65

66

Kapitola 3 Porovn´ an´ı re´ aln´ ych dat s teoretick´ ym modelem V této kapitole budou porovnána data z´ıskaná ze 4 tˇr´ıdn´ıch knih z obdob´ı prob´ıhaj´ıc´ı epidemie spálové ang´ıny zaˇcátkem ˇskoln´ıho roku 2014/2015 na druhém stupni jedné ze základn´ıch ˇskol v Sokolovˇe.

3.1

Sp´ alov´ a ang´ına

Spálová ang´ına je infekˇcn´ı onemocnˇen´ı zp˚ usobené streptokoky. Projevuje se horeˇckou, vyráˇzkou a akutn´ım zánˇetem krˇcn´ıch mandl´ı a okoln´ı lymfatické tkánˇe. M˚ uˇze m´ıt m´ırn´ y i velmi závaˇzn´ y pr˚ ubˇeh, typicky se vyskytuje v menˇs´ıch epidemi´ıch v dˇetsk´ ych kolektivech. Inkubaˇcn´ı doba je vˇetˇsinou 2-5 dn˚ u, postiˇzen´ y b´ yvá izolován a léˇcen antibiotiky (10-14 dn˚ u). Pˇri zanedbán´ı léˇcby existuje riziko pozdˇejˇs´ıch komplikac´ı. Tyto a dalˇs´ı informace si m˚ uˇzete pˇreˇc´ıst v [8] a v [10].

3.1.1

Volba modelu

Data, která máme k dispozici z tˇr´ıdn´ıch knih, pˇredstavuj´ı ˇzáky v domác´ı izolaci, tedy tˇr´ıdu H, a protoˇze spálová ang´ına je nemoc s latentn´ım obdob´ım, potˇrebujeme také tˇr´ıdu E. Model obsahuj´ıc´ı obˇe tyto tˇr´ıdy nebyl v pˇrehledu základn´ıch epidemiologick´ ych model˚ u (kapitola 1) uveden, proto bylo zapotˇreb´ı vytvoˇrit model nov´ y. Spojen´ım SEIR (kapitola 1.8.5, str. 34) a SIHR (kapitola 1.9, str. 36) modelu vnikl následuj´ıc´ı model: 67

dS dt dE dt dI dt dH dt dR dt

= −βS(I + εE),

(3.1.1)

= βS(I + εE) − σE,

(3.1.2)

= σE − γI − θI,

(3.1.3)

= θI − δH,

(3.1.4)

= γI + δH.

(3.1.5)

Parametr β pˇredstavuje pravdˇepodobnost nakaˇzen´ı pˇri kontaktu náchylné osoby s infekˇcn´ım jedincem, konstanta γ udává pr˚ umˇernou rychlost uzdraven´ı, koeficient ε sniˇzuje infekˇcnost jedinc˚ u v latentn´ım obdob´ı oproti infekˇcn´ım jedinc˚ um, konstanta σ vyjadˇruje rychlost pˇrechodu z infikovan´ ych jedinc˚ u mezi infekˇcn´ı, konstanta θ urˇcuje rychlost, s jakou se infekˇcn´ı jedinci po propuknut´ı nemoci pˇresunou do domác´ı izolace, a konstanta δ udává rychlost pˇrechodu jedinc˚ u v domác´ı izolaci mezi imunn´ı. Rovnice (3.1.1) ˇr´ıká, ˇze ˇca´st náchyln´ ych jedinc˚ u se nakaz´ı kontaktem s infikovan´ ymi a infekˇcn´ımi jedinci. Tato ˇcást pak pˇrejde mezi infikované a ti se po uplynut´ı obdob´ı latence pˇresunou mezi infekˇcn´ı (rovnice (3.1.2)). Rovnice (3.1.3) zachycuje tento pˇresun infikovan´ ych osob mezi infekˇcn´ı a také pˇresun ˇca´sti infekˇcn´ıch jedinc˚ u do domác´ı izolace, zbytek osob ve tˇr´ıdˇe z˚ ustane aˇz do odeznˇen´ı nemoci a potom se pˇresune mezi imunn´ı. Jedinci, kteˇr´ı se léˇcili doma, se po skonˇcen´ı izolace pˇresunou také mezi imunn´ı (rovnice (3.1.4)). Pˇredpokládáme tedy, ˇze izolace trvá aˇz do doby, kdy jedinec nen´ı schopn´ y nemoc dál pˇrenáˇset. Rovnice (3.1.5) ukazuje, jak pˇrib´ yvá imunn´ıch jedinc˚ u. Spálová ang´ına nezanechává trvalou imunitu, ale vzhledem ke krátkému ˇcasovému u ´seku nepˇredpokládáme, ˇze by se nˇekdo nakazil opakovanˇe. Pˇrestoˇze se jedná o zjednoduˇsenou pˇredstavu o pr˚ ubˇehu onemocnˇen´ı, mˇeli bychom uˇz s t´ımto modelem pˇri dobˇre zvolen´ ych hodnotách parametr˚ u dosáhnout relativnˇe dobr´ ych v´ ysledk˚ u.

68

3.2

Modelov´ an´ı pr˚ ubˇ ehu n´ akazy

Obrázek 3.1: Poˇcty chybˇej´ıc´ıch ˇzák˚ u v pr˚ ubˇehu epidemie spálové ang´ıny – data z tˇr´ıdn´ıch knih. Ve tˇr´ıdách, ve kter´ ych se vyskytla spálová ang´ına a pro které bylo provedeno srovnán´ı teoretick´ ych a reáln´ ych hodnot, bylo celkem 99 ˇza´k˚ u, pro modelován´ı tedy byly stanoveny poˇcty ˇza´k˚ u v modelov´ ych tˇr´ıdách následovnˇe: S0 = 97,

E0 = 2,

I0 = 0,

H0 = 0,

R0 = 0.

Hodnoty parametr˚ u jsme pak na základˇe informac´ı, které o spálové ang´ınˇe a chován´ı ˇzák˚ u máme, urˇcili takto: • γ=

1 , délka nemoci je totiˇz obvykle aˇz 2 t´ ydny, 14

1 • ε = , tento parametr byl odhadnut tak, aby se vypoˇctené hodnoty co 2 nejlépe shodovaly s reáln´ ymi daty, nebot’ podklady pro jeho stanoven´ı nejsou nikde uvádˇeny, 1 • σ = , protoˇze latentn´ı stadium nemoci trvá pr˚ umˇernˇe 4 dny, 4 • θ=

9 , vˇetˇsina dˇet´ı totiˇz z˚ ustává doma jiˇz pˇri prvn´ıch pˇr´ıznac´ıch, 10 69

• δ=

1 , nebot’ ˇzaći obvykle chybˇeli 10 dn´ı, 10

4/5 ve vˇsedn´ı dny, v souladu s daty totiˇz pˇredpokládáme pomˇernˇe • β= 99 vysokou pravdˇepodobnost nakaˇzen´ı pˇri kontaktu s pˇrenaˇseˇcem nemoci, pro jednoduchost ovˇsem uvaˇzujeme stále pln´ y poˇcet jedinc˚ u ve tˇr´ıdách, • β = 0 o v´ıkendech, protoˇze pˇredpokládáme, ˇze se dˇeti mimo ˇskolu nest´ ykaj´ı a tedy se od sebe nemohou nakazit. Graf pro model s tˇemito parametry je na obrázku 3.2.

Obrázek 3.2: Modelován´ı pr˚ ubˇehu epidemie spálové ang´ıny. Tˇr´ıdn´ı knihy obsahuj´ı pouze u ´daje pro vˇsedn´ı dny (viz obrázek 3.1). Aby byl zˇretelnˇejˇs´ı trend chybˇen´ı ˇzák˚ u na obrázku 3.3, byly doplnˇeny oˇcekávané poˇcty dˇet´ı v domác´ı izolaci o v´ıkendech. Ty byly stanoveny jako aritmetick´ y pr˚ umˇer dané páteˇcn´ı a pondˇeln´ı hodnoty, v pˇr´ıpadˇe nutnosti zaokrouhlen´ y na celé ˇc´ıslo dle pravidel zaokrouhlován´ı. Ve ˇctyˇrech tˇr´ıdách obvykle pr˚ umˇernˇe chyb´ı 3 aˇz 5 ˇzák˚ u, proto bylo k vypoˇcten´ ym hodnotám funkce H pro lepˇs´ı porovnán´ı s reáln´ ymi daty pˇriˇc´ıtáno ˇc´ıslo 4, které pˇredstavuje ˇza´ky chybˇej´ıc´ı z jiného d˚ uvodu neˇz je onemocnˇen´ı spálovou ang´ınou. Obrázek 3.3 tedy ukazuje doplnˇená reálná data, jejich aproximaci polynomem 7. stupnˇe a teoretické poˇcty ˇzák˚ u ve tˇr´ıdˇe H (zv´ yˇsené 70

o 4). Stupeˇ n polynomu byl zvolen tak, aby byl dobˇre patrn´ y trend chybˇen´ı ˇza´k˚ u a zároveˇ n nedoˇslo k jeho pˇr´ıliˇsnému rozvlnˇen´ı kv˚ uli velk´ ym v´ ykyv˚ um reáln´ ych dat. Pˇrestoˇze jsou hodnoty ˇcasto znaˇcnˇe rozd´ılné, m˚ uˇzeme ˇr´ıci, ˇze jsme d´ıky modelován´ı z´ıskali relativnˇe dobrou pˇredstavu o pr˚ ubˇehu epidemie. D˚ uvody odliˇsnosti v´ ysledk˚ u jsou podrobnˇeji rozebrány v následuj´ıc´ı sekci.

Obrázek 3.3: Poˇcty chybˇej´ıc´ıch ˇza´k˚ u doplnˇené o oˇcekávané hodnoty ˇza´k˚ u v domác´ı izolaci o v´ıkendech, aproximace tˇechto dat polynomem 7. stupnˇe pro lepˇs´ı znázornˇen´ı trendu chybˇen´ı ˇza´k˚ u a teoretické hodnoty chybˇej´ıc´ıch ˇza´k˚ u dle SEIHR modelu.

3.2.1

´ Uskal´ ı re´ aln´ ych dat

Rozd´ılnost teoretick´ ych a experimentáln´ıch dat je zp˚ usobena t´ım, ˇze u takto n´ızkého poˇctu jedinc˚ u je d˚ uleˇzité chován´ı kaˇzdého z nich. Do modelu nejsme schopni zahrnout vˇsechny situace, které reálnˇe nastávaj´ı, napˇr. ˇze nemocné dˇeti se na jeden den do ˇskoly vrát´ı kv˚ uli d˚ uleˇzité p´ısemce anebo na v´ıce dn´ı proto, ˇze rodiˇce nesehnali hl´ıdán´ı. Stejnˇe tak nejsme schopni do modelu zahrnout r˚ uznorodost kontaktn´ıch struktur na základˇe sociáln´ıch vazeb (vyˇsˇs´ı m´ıra kontaktu ve skupinkách kamarád˚ u, jedinci vylouˇcen´ı z kolektivu apod.), která je u takto malého poˇctu dˇet´ı v´ yznamná. Data z tˇr´ıdn´ıch knih také neodhal´ı, jestli d˚ uvodem absence ˇza´ka byla právˇe spálová ang´ına, pˇr´ıpadnˇe jiná nemoc nebo rodinné d˚ uvody. Model také neuvaˇzuje moˇznost, ˇze by se dˇeti nakazily mimo ˇskoln´ı prostˇred´ı a ˇze se pravdˇepodobnost nakaˇzen´ı v pr˚ ubˇehu 71

epidemie mˇen´ı. Takto bychom mohli jmenovat i dalˇs´ı faktory ovlivˇ nuj´ıc´ı kvalitu modelu. Velk´ y vliv sehrála skuteˇcnost, ˇze od pondˇel´ı 27. 10. do stˇredy 29. 10. mˇely dˇeti prázdniny. Mnohé z nich jely s rodiˇci na dovolenou, coˇz je d˚ uvodem vysokého poˇctu absenc´ı ve stˇredu, ˇctvrtek a pátek v t´ ydnu pˇred prázdninami, a ve ˇctvrtek a v pátek v prázdninovém t´ ydnu. Proto na obrázku 3.1 nedocház´ı ke konci epidemie k takovému poklesu chybˇej´ıc´ıch ˇzák˚ u, jak bychom oˇcekávali.

72

Pˇ r´ıloha A

A.1

Abelova rovnice

V této sekci bylo ˇcerpáno z [7]. Abelova rovnice je obyˇcejná diferenciáln´ı rovnice, kubická v neznámé funkci, pojmenovaná podle N. H. Abela. Obvykle uvaˇzujeme Abelovu rovnici prvn´ıho nebo druhého ˇrádu. Nelineárn´ı Abelova diferenciáln´ı rovnice prvn´ıho ˇra´du má podobu dy = p(x)y 3 + q(x)y 2 + r(x)y + s(x), p(x) 6= 0, (A.1.1) dx a hraje d˚ uleˇzitou roli v mnoha fyzikáln´ıch a technick´ ych aplikac´ıch. Nen´ı tˇeˇzké nalézt jej´ı partikulárn´ı ˇreˇsen´ı, uˇz sloˇzitˇejˇs´ı je nalezen´ı ˇreˇsen´ı obecného. Pˇresto existuje celá ˇrada postup˚ u, jak obecné ˇreˇsen´ı z´ıskat, napˇr. pomoc´ı substituce E(x) , u = y − y1 (x) R 2 E(x) = e [3p(x)y1 +2q(x)y1 +r(x)]dx , kde y1 (x) je nám známé partikulárn´ı ˇreˇsen´ı, a následnou transformac´ı rovnice (A.1.1) do podoby du p(x)E 2 (x) + = −E(x) [3p(x)y1 (x) + q(x)] . dx u

(A.1.2)

Pokud y1 = −

q(x) , 3p(x)

pravá strana rovnice (A.1.2) je nulová a obecné ˇreˇsen´ı rovnice (A.1.1) lze z´ıskat integrac´ı diferenciáln´ı rovnice se separovan´ ymi promˇenn´ ymi.

73

A.2

Bernoulliho rovnice

Bernoulliho rovnic´ı se rozum´ı nelineárn´ı diferenciáln´ı rovnice dy + p(x)y = q(x)y n , dx kterou se zab´ yval matematik Jacob Bernoulli. Jej´ı obecné ˇreˇsen´ı bylo objeveno jiˇz koncem 17. stolet´ı a má tvar  " #1/(1−n) R (1−n) R p(x)dx   (1 − n) e q(x)dx + C 1   R pro n 6= 1,   (1−n) p(x)dx  e y=      R   C2 e [q(x)−p(x)]dx pro n = 1, kde C1 a C2 jsou integraˇcn´ı konstanty.

74

Literatura [1] Brauer, F., Castillo-Chavez, C.: Mathematical models in population biology and epidemiology. Springer, New York, druhé vydán´ı 2012. [2] Brauer, F., Van den Driessche, P., Wu, J.: Mathematical epidemiology. Springer, Berl´ın, 2008. [3] Harko, T., Lobo, F.S.N., Mak, M.K.: Exact analytical solutions of the Susceptible-Infected-Recovered (SIR) epidemic model and of the SIR model with equal death and birth rates. Applied Mathematics and Computation 236, 2014, str. 184–194. [4] Ma, Z., Li, J.: Dynamical modeling and analysis of epidemics. World Scientific Publishing, Singapur, 2009. [5] Brauer, F.: Age-of-infection and the final size relation. Mathematical biosciences and engineering, Volume 5, Number 4, October 2008, str. 681690. [6] Neely, D.: Quarantine for an infectious disease. [online]. [cit. 2014-1211]. Dostupné na: http://www.westminster.edu/acad/math/dept/pdf/neely capstone paper.pdf [7] Mak, M.K., Harko, T.: New method for generating general solution of Abel differential eqution. Computers and mathematics with applications 43, 2002, str. 91-94. [8] Vokurka, M., Hugo, J.: Praktický slovn´ık medic´ıny. Maxdorf, Praha, sedmé vydán´ı 2004. [9] Drábek, P.: Základy matematické teorie nelineárn´ıch dynamických ˇ systém˚ u. Editaˇcn´ı stˇredisko VSSE, Plzeˇ n, 1990. [10] Spálová ang´ına. WikiSkripta. [online]. [cit. 2015-04-17]. Dostupné z: www.wikiskripta.eu

75

Západočeská univerzita v Plzni Fakulta aplikovaných věd Katedra matematiky. Matematické modely v epidemiologii

Recommend Documents