Numeriˇ cno modeliranje ˇ sirjenja bolezni COVID - 19

UNIVERZA V LJUBLJANI

FAKULTETA ZA MATEMATIKO IN FIZIKO Finanˇ cna matematika – 1. stopnja

Jan ˇ Sifrer

Numeriˇ cno modeliranje ˇ sirjenja bolezni COVID - 19

Delo diplomskega seminarja

Mentorica: prof. dr. Marjetka Knez

Ljubljana, 2021

Kazalo

1. Uvod 4

2. Opis zaˇcetnega problema diferencialnih enaˇcb ter metod 5

2.1. Runga-Kutta metode 7

3. SIR model 8

3.1. Razˇsirjen SIR model 11

4. θ-SEIHRD model 12

4.1. Predpostavke modela 12

4.2. Matematiˇcen zapis modela 13

4.3. Rezultati modela 17

4.4. Ocenjevanje parametrov 18

4.5. Numeriˇcna simulacija in konfiguracija modela za primer Kitajske 20

5. Zakljuˇcek 27

Literatura 28

Numeriˇcno modeliranje ˇsirjenja bolezni COVID - 19

Povzetek

V tem delu diplomskega seminarja bomo govorili o temi, ki je nekako zaznamovala leti 2020 in 2021. Govorili bomo namreˇc o epidemiji koronavirusa COVID-19 in kako bi se ˇsirjenje le te dalo napovedati s pomoˇcjo numeriˇcnih modelov. Bolj natanˇcno, spoznali bomo tri razliˇcne modele. Najprej se bomo seznanili s SIR modelom, ki spada med najbolj osnovne modele za napovedovanje ˇsirjenja nalezljivih bolezni.

Pove nam, kako bi se epidemija ˇsirila, ˇce ne bi vanjo niˇc posegali. Opisuje prehode ljudi med tremi stanji - dovzetni, okuˇzeni in ozdravljeni. V nadaljevanju bomo spo- znali razˇsirjeni SIR model, kjer bomo upoˇstevali, da so se ob doloˇcenih trenutkih sprejeli doloˇceni ukrepi za zajezitev ˇsirjenja virusa COVID-19. Videli bomo, kako je to vplivalo na sam potek ˇsirjenja epidemije. Tretji model, kateremu bomo na- menili najveˇc besed, pa je θ-SEIHRD model, ki opisuje prehajanje ljudi med kar 9-timi razliˇcnimi stanji ob upoˇstevanju razliˇcnih ukrepov za zajezitev ˇsirjenja virusa COVID-19. Dodatna stanja prinesejo veˇc enaˇcb in parametrov, ki jih je potrebno doloˇciti, kar pa se izkaˇze za teˇzak problem, saj bomo videli, da so rezultati zelo odvisni od natanˇcnosti nastavljenih parametrov.

Numerical modeling of COVID-19 disease spread

Abstract

In this seminar we consider a topic that somehow marked the years 2020 and 2021.

Namely, we talk about the COVID-19 coronavirus epidemic and how its spread could be predicted with the help of various numerical models. More precisely, we get to know three different models. First, we describe the SIR model, which is one of the most basic models for predicting the spread of infectious diseases. It tells us how the epidemic would spread if nothing would be done about it. It describes the transitions of people between three conditions - susceptible, infected and cured. In the following, we present the extended SIR model, where we take into account that certain measures are taken at certain moments to curb the spread of the COVID- 19 virus. We show how this affects the very course of the spread of the epidemic.

The third model, to which we devote most words, is the θ -SEIHRD model, which describes the transition of people between 9 different conditions, taking into account various measures to curb the spread of the virus. Here we show that additional conditions bring more equations and parameters that need to be determined, which turns out to be a difficult problem as the results are highly dependent on the accuracy of the chosen parameters.

Math. Subj. Class. (2010): 34A12, 65L05, 65L06

Kljuˇcne besede: Runge-Kutta metoda, epidemija, COVID-19, numeriˇcna simula- cija, SIR model, θ-SEIHRD model.

Keywords: Runge-Kutta method, coronavirus epidemic, COVID-19, numerical si- mulation, SIR model, θ-SEIHRD model.

1. Uvod

Koronavirusi so skupina virusov, ki povzroˇcajo bolezni pri ljudeh in ˇzivalih. Pri ljudeh koronavirusi povzroˇcajo okuˇzbe dihal, ki pa so lahko blage - podobne na- vadnemu prehladu, lahko pa so tudi zelo hude, celo smrtne, kot so npr. SARS (povzroˇca ga virus SARS-CoV), MERS (povzroˇca ga virus MERS-CoV), COVID-19 (povzroˇca ga virus SARS-CoV-2). Najveˇcji problem pa predstavljajo ˇcloveˇstvu zato, ker za povzroˇcene bolezni ˇse ne obstajajo nobena odobrena zdravila.

COVID-19 je torej bolezen, ki je nekako zaznamovala leti 2020 in 2021. Prvi primer okuˇzbe so zaznali decembra 2019 v Kitajskem mestu Wuhan, od koder se je virus zaˇcel ˇsiriti po vsej Kitajski in ˇse po veliko drugih drˇzavah po svetu. Dne, 11.

marca 2020, je bilo okuˇzenih ˇze veˇc kot 118,000 ljudi v 114 razliˇcnih drˇzavah, od tega je bilo veˇc kot 90 % okuˇzb v samo dveh drˇzavah (na Kitajskem in v Juˇzni Koreji), zato je Svetovna zdravstvena organizacija (WHO) razglasila pandemijo. Virus se med ljudmi ˇsiri s pomoˇcjo aerosolov, ki jih v zrak izpuˇsˇcamo ˇze samo z govorjenjem, velik problem pa je tudi, da lahko ˇclovek virus ˇsiri ˇze nekaj dni, preden sploh dobi kakrˇsne koli simptome okuˇzbe. Pri najveˇc primerih okuˇzbe (80 % vseh primerov) so simptomi vroˇcina, suh kaˇselj in utrujenost ter izguba vonja in okusa. Pri nekaterih primerih (20 % vseh primerov) pa se razvijejo ˇse dodatni - hujˇsi simptomi, kot so zasoplost, izguba apetita, zmedenost, stalne boleˇcine ali pritiski v prsih, vroˇcina (nad 38 stopinj). Taki primeri so hospitalizirani in potrebujejo tretmaje kisika za premagovanje bolezni, 5 % hospitaliziranih pa potrebuje celo intenzivno nego - to so po navadi tisti primeri, ko okuˇzeni ne morejo veˇc sami dihati.

In ravno tu se kaˇze eden od problemov epidemij, s katerim se sooˇcajo strokovnjaki in politiki - ˇstevilo osebja in postelj v bolnicah je zgolj konˇcno mnogo. Kaj naj torej naredijo, da ne pride do kolapsa sistema? Natanˇcneje, da ne pride do tega, da bi nekdo potreboval medicinsko oskrbo, a je ne bi mogel dobiti, ker zanj ne bi bilo proste postelje ali osebja. Za pomoˇc pri razreˇsevanju tega problema je kar nekaj strokovnih skupin in raziskovalcev razvilo razliˇcne modele, s katerimi poskuˇsajo napovedati, kako se bo situacija razvijala v prihodnosti.

V delu diplomskega seminarja se bomo posvetili bolj matematiˇcnim modelom, ki temeljijo na numeriˇcnem reˇsevanju sistema navadnih diferencialnih enaˇcb. Bolj natanˇcno, najprej bodo opisane metode, ki se najveˇckrat uporabljajo v teh mate- matiˇcnih modelih. Nato pa bo opisan model, ki so ga uporabili na Kitajskem.

2. Opis zaˇcetnega problema diferencialnih enaˇcb ter metod Veliko znanstvenih in tudi tehnoloˇskih problemov se da opisati s pomoˇcjo mate- matiˇcnih modelov, ki temeljijo na reˇsevanju sistemov diferencialnih enaˇcb. Ker pa veˇcina takih problemov nima analitiˇcnih reˇsitev, jih reˇsujemo s pomoˇcjo numeriˇcnih metod. Zgodovinsko gledano sta prednika vseh danes uporabljenih metod za nu- meriˇcno reˇsevanje diferencialnih enaˇcb Eulerjeva metoda, ki jo je leta 1768 razvil Svicarski matematik Leonhard Euler, in Runge-Kutta metode, ki sta jih razvila Carlˇ Runge in Martin Kutta leta 1901 ([11]).

Numeriˇcne metode se da razdeliti na dva razreda:

(1) Diferenˇcne metode, kjer pribliˇzek k reˇsitvi iˇsˇcemo diskretno.

(2) Pribliˇzek k reˇsitvi iˇsˇcemo kot funkcijo v nekem izbranem podprostoru.

Najbolj osnovna naloga je poiskati reˇsitev t. i. zaˇcetnega problema diferencialnih enaˇcb oblike y^′ = f(x, y) pri danem zaˇcetnem pogoju y(a) = ya. ˇCe se posvetimo diferenˇcnim metodam, potem je reˇsitev sestavljena iz zaporedja toˇck a=x0 < x1 <

x₂ < · · · < x_M = b in pripadajoˇcega zaporedja vrednosti y₀, y₁, y₂, . . . , y_M, kjer je y_n izraˇcunan pribliˇzek za toˇcno reˇsitev vx_n, to je y_n≈y(x_n).

Diferenˇcne metode delimo na:

• Enoˇclenske metode: pribliˇzeky_n+1 je izraˇcunan iz pribliˇzka y_n.

• Veˇcˇclenske (k-ˇclenske) metode: pribliˇzek yn+1 je izraˇcunan iz pribliˇzkov yn, yn−1, yn−2, . . . , yn−k+1.

Glede na naˇcin izraˇcuna pribliˇzkov, pa jih delimo na:

• Eksplicitne metode: imamo direktno formulo za izraˇcun pribliˇzka y_n+1.

• Implicitne metode: pribliˇzek y_n+1 dobimo tako, da reˇsimo enaˇcbo, ki je obiˇcajno nelinearna.

Definicija 2.1. Sistem diferencialnih enaˇcb 1. reda je sestavljen iz diferencialnih enaˇcb in zaˇcetnih pogojev:

dy1

dx =f₁(x, y₁, y₂, . . . , y_m), y₁(a) = y_1,a, dy₂

dx =f₂(x, y₁, y₂, . . . , y_m), y₂(a) = y_2,a,

... ...

dy_m

dx =f_m(x, y₁, y₂, . . . , y_m), y_m(a) = y_m,a,

kjer so f_i : [a, b]×R^m →Rdane funkcije, y_i : [a, b]→R pa funkcije, ki jih moramo izraˇcunati. Tak problem tudi imenujemo zaˇcetni problem diferencialnih enaˇcb.

Zaradi same poenostavitve definirajmo ˇse vektorski funkciji

y=

⎡

⎢

⎢

⎣ y_m

... y2

y₁⎤

⎥

⎥

⎦

, y: [a, b]→R^m,

f=

⎡

⎢

⎢

⎣ f_m

... f₂ f₁⎤

⎥

⎥

⎦

, f: [a, b]×R^m →R^m.

Torej lahko zaˇcetni problem diferencialnih enaˇcb vektorsko zapiˇsemo kot y^′ =f(x, y),

y(a) = y_a =

⎡

⎢

⎢

⎣ y_m,a

... y_2,a y_1,a⎤

⎥

⎥

⎦ .

V tem diplomskem delu bomo take sisteme diferencialnih enaˇcb reˇsevali numeriˇcno s pomoˇcjo Runga-Kutta metode stopnje 4 (RK4). Preden pa nadaljujemo, pa se ˇse prepriˇcajmo, da pri primernih pogojih na neki okolici a res dobimo natanko eno reˇsitev prej zapisanega sistema. To izvemo iz eksistenˇcnega izreka za sisteme navadnih diferencialnih enaˇcb ([5]).

Izrek 2.2. Naj bodo fk, k = 1, . . . , m, funkcije m + 1 spremenljivk, ki so na paralelepipedu

|x−a| ≤b₀,

|y_k−y_k,a| ≤b_k, k = 1, . . . , m, zvezne in zadoˇsˇcajo Lipschitzevemu pogoju

|f_k(x, y₁, . . . , y_m)−f_k(x, η₁, . . . , η_m| ≤C

m

∑︂

j=1

|y_j −η_j|, k = 1, . . . , m.

Naj bo

M₀ = max{|f_k(x, y₁, . . . , y_m)|:|x−a| ≤b₀, |y_k−y_k,a| ≤b_k, k = 1, . . . , m}.

Tedaj ima problem

y^′_k=f_k(x, y₁, . . . , y_m) y_k(a) =y_k,a, k= 1, . . . , m, na intervalu

|x−a| ≤c= min {︃

b₀, b₁ M0

, . . . , b_n M0

, α mC

}︃

,0< α <1, natanko eno reˇsitev.

Dokaz izreka lahko najdemo v viru [5].

2.1. Runga-Kutta metode. Najbolj popularen numeriˇcni pristop za reˇsevanje zaˇcetnih problemov diferencialnih enaˇcb oziroma sistemov diferencialnih enaˇcb je pristop z uporabo Runga-Kutta metod. Te metode so namreˇc precej natanˇcne, nu- meriˇcno stabilne in se jih da preprosto izraˇcunati ter implementirati. Spadajo pa v razred diferenˇcnih metod.

Ideja teh metod je, da izraˇcunamo pribliˇzno vrednost funkcijefv stoˇckah na nekem intervalu [xn, xn+1] ⊂ [a, b]. To nam da vrednosti odvodov y v teh toˇckah, kar pa doloˇca premike v smeri tangente (pri eni diferencialni enaˇcbi). Konˇcen premik sestavimo iz uteˇzenega povpreˇcja dobljenih premikov. Za laˇzje razumevanje si kar poglejmo, kako izgleda Runga-Kutta metoda stopnje 4.

2.1.1. Runga-Kutta metoda stopnje 4. Za izraˇcun y_n+1, ki je izraˇcunan pribliˇzek za toˇcno reˇsitev v x_n+1, potrebujemo 4 parametre, ki jih oznaˇcimo s k_i, i= 1, . . . , 4.

Toˇckea =x₀ < x₁ < x₂ <· · ·< x_M =bso pri osnovni razliˇcici te metode razdeljene ekvidistantno, kar pomeni, da je razdalja med njimi konstantna, poimenujemo jo lahko tudi korak metode in jo oznaˇcimo s h = ^b−a_M . Parametre k_i izraˇcunamo s pomoˇcjo naslednjih ˇstirih enaˇcb:

k₁ =h·f(x_n, y_n) k₂ =h·f(x_n+ h

2, y_n+ 1 2·k₁) k₃ =h·f(x_n+ h

2, y_n+ 1 2·k₂) k₄ =h·f(x_n+h, y_n+k₃).

Pribliˇzek y_n+1 pa izraˇcunamo tako, da te parametre ustrezno uteˇzimo:

y_n+1 =y_n+1

6 ·(k₁+ 2·k₂+ 2·k₃+k₄).

To metodo bi v Matlabu zapisali na naslednji naˇcin.

1 f u n c t i o n [ x , y ] = RungeKutta ( fun , a , b , y0 , h )

2 % Runge=Kutta metoda r e d a 4 za r e s e v a n j e ( s i s t e m o v ) d i f e r e n c i a l n i h enacb p r v e g a r e d a

3 % Vhodni p o d a t k i : f u n k c i j a f ( x , y ) ; i n t e r v a l [ a , b ] , na katerem i s c e m o r e s i t e v ;

4 % z a c e t n i p o g o j y0 ( v e k t o r ) ; k o r a k h

5

6 x = a : h : b ;

7 y=z e r o s(l e n g t h( y0 ) ,l e n g t h( x ) ) ;

8 y ( : , 1 ) = y0 ;

9 f o r i = 2 :l e n g t h( x )

10 k1 = h*f u n ( x ( i=1) , y ( : , i=1) ) ;

11 k2 = h*f u n ( x ( i=1)+ 1/2*h , y ( : , i=1) + 1/2*k1 ) ;

12 k3 = h*f u n ( x ( i=1) + 1/2*h , y ( : , i=1) + 1/2*k2 ) ;

13 k4 = h*f u n ( x ( i=1) + h , y ( : , i=1) + k3 ) ;

14 y ( : , i ) = y ( : , i=1) + 1 / 6*( k1 + 2*k2 + 2*k3 + k4 ) ;

15 end

Primer 2.3. Zanima nas reˇsitev diferencialne enaˇcbe y^′ =−y−5 exp(−x) sin(5x) na intervalu [0, 3] pri zaˇcetnem pogoju y(0) = 1 in koraku h= 0.3. Problem bomo reˇsili s pomoˇcjo Runga-Kutta metode stopnje 4.

Spodnja koda prikazuje, kako v Matlabu definiramo vse potrebne vhodne podatke za Runge-Kutta metodo.

1 a =0;

2 b=3;

3 h = 0 . 3 ;

4 f u n = @( x , y ) (=y =5*exp(=x )*s i n( 5*x ) ) ;

5 y0 = 1 ;

6 [ x1 , y1 ] = RungeKutta ( fun , a , b , y0 , h ) ;

Metoda vrne le pribliˇzke v izbranih ekvidistantnih toˇckah. Aproksimacijo iskane funkcije y dobimo, ˇce te pribliˇzke interpoliramo z odsekoma linearno funkcijo. Na grafu to pomeni, da toˇcke poveˇzemo z daljicami. Na sliki 1 lahko vidimo primerjavo med toˇcno reˇsitvijo (rdeˇci graf), ki jo dobimo, ˇce problem reˇsimo analitiˇcno, in numeriˇcno izraˇcunanim pribliˇzkom (moder graf) s korakom h = 0.3 in korakom h= 0.1. Vidimo, da se krivulji kar dobro ujemata. Veˇcjo natanˇcnost pa dobimo, ˇce

zmanjˇsamo korak h. ♢

Slika 1. Primerjava analitiˇcne reˇsitve (rdeˇc graf) s reˇsitvijo izraˇcunano z Runga-Kutta metode stopnje 4 (moder graf) s korakom h= 0.3 (levo) in korakom h= 0.1 (desno).

3. SIR model

Modeli, ki nam dajejo moˇznost napovedovanja ˇsirjenja epidemije, so nekako osnova za uspeˇsno poznavanje ˇsirjenja le te in sprejemanje uˇcinkovitih ukrepov zoper nje.

Za laˇzje razumevanje, kako delujejo taki modeli, si najprej poglejmo najosnovnejˇsi model, ki se imenuje SIR model, ki opisuje, kako se ˇsiri epidemija, ˇce vanjo ne bi niˇc posegali. Natanˇcneje, poglejmo si, kako bi se ˇsiril COVID-19, ˇce ne bi bilo nobenih ukrepov. Sam model privzema, da so vsi ljudje v enem izmed naslednjih treh stanj:

• dovzetni za okuˇzbo (S),

• okuˇzeni (I),

• ozdravljeni (R).

Posebej pripomnimo ˇse, da za ta stanja velja N =S(t) +I(t) +R(t) za vsak t, kjer je N ˇstevilo ljudi, t pa nek poljuben ˇcas. ˇCe torej odvajamo ta izraz po ˇcasu t, dobimo 0 = ^dS_dt + ^dI_dt + ^dR_dt, kar torej pomeni, da sprememba enega stanja povzroˇci spremembo drugih dveh stanj.

Ostale predpostavke SIR modela pa so:

• hitrost ˇsirjenja virusa je odvisna od velikosti populacije,

• ni obdobja inkubacije - torej ob kontaktu oseba bodisi zboli, bodisi ostane zdrava,

• na zaˇcetku ni nihˇce imun na bolezen,

• ko postane imun, to za vedno ostane,

• bolezen ni usodna.

Ker bolezen ni usodna, lahko doloˇcimo verjetnost ozdravitve v posamezni enoti merjenja, kar oznaˇcimo zγ, njen inverz _γ¹ pa nam pove povpreˇcen ˇcas kuˇznosti. Mo- del dalje privzema, da je verjetnost, da okuˇzeni posameznik pride v stik z neokuˇzenim posameznikom nakljuˇcna, kar oznaˇcimo s parametrom c, verjetnost prenosa okuˇzbe ob kontaktu pa oznaˇcimo z a. ˇCe torej pomnoˇzimo obe konstanti, dobimo novo konstanto in jo oznaˇcimo zβ, tj. β =a·c. Ta konstantaβ nam pove, s kolikˇsno ver- jetnostjo okuˇzeni posameznik okuˇzi nekoga drugega v eni enoti merjenja. ˇCe torej pomnoˇzimo β s povpreˇcnim ˇcasom kuˇznosti _γ¹, dobimo ˇstevilo, kateremu pravimo tudi reprodukcijsko ˇsteviloR₀ =β·¹_γ, ki nam pove, koliko ljudi okuˇzena oseba okuˇzi v svojem kuˇznem obdobju ([12]).

Sedaj pa si poglejmo, kako s sistemom diferencialnih enaˇcb opiˇsemo povezave oz.

prehode med temi tremi stanji po osnovnem SIR modelu:

dS

dt =−βIS N , (1)

dI

dt = βIS N −γI, (2)

dR dt =γI.

(3)

Tak sistem diferencialnih enaˇcb bomo reˇsevali s pomoˇcjo numeriˇcnih metod oz. bolj natanˇcno z Runga-Kutta metodo stopnje 4.

Na sliki 2 lahko vidimo reˇsitev v primeru, ko sta koeficienta β in γ poznana (β = 0.4, γ = 0.1), uporabimo korak h= 0.1 in je v opazovani populaciji 5000 ljudi. Za zaˇcetno stanje izberemo S(0) = 4999, I(0) = 1, R(0) = 0.

Opazimo lahko eno pomembno stvar, in sicer da se vsaka taka epidemija prej ali slej konˇca. Pomembno je poudariti, da to ni zato, ker se vsi ljudje prekuˇzijo, temveˇc zato, ker se doseˇze skupinsko imunost. Na sliki 2 lahko vidimo, da je na koncu epidemije ˇse vedno nekaj ljudi v stanju S - skoraj 100 ljudi oz. skoraj 2 % populacije. Vendar se tukaj lahko vpraˇsamo, kaj bi to pomenilo oz. kakˇsno ceno bi morali plaˇcati, ˇce bi pustili epidemiji prosto pot. Najbolj problematiˇcna krivulja je seveda krivulja okuˇzenih (na sliki 2 oznaˇcena z rdeˇco barvo), saj je za nekatere okuˇzene potrebno ustrezno poskrbeti, torej je zelo pomembno, da v nekem trenutku to ˇstevilo okuˇzenih ne preseˇze zdravstvenih kapacitet. Torej nas zanima, kolikˇsno bo najveˇcje ˇstevilo okuˇzenih v nekem trenutku. To pa izraˇcunamo na naslednji naˇcin.

Slika 2. Potek epidemije brez ukrepanja.

Delimo enaˇcbo (2) z enaˇcbo (1) in to novo zvezo integriramo:

dI dS =

βIS N −γI

−^βIS_N =−1 + γN βS,

∫︂ dI dS dS =

∫︂

−1 + γN βS dS.

Dobimo

I =−S+γN

β ln(S) +C,

kjer je C neka poljubna konstanta. Od tod sledi, da velja zveza I(t) +S(t)− γN

β ln(S(t)) = C.

(4)

Konstanto C doloˇcimo iz zaˇcetnih pogojev C =I(0) +S(0)−γN

β ln(S(0)).

(5)

Iz (4) in (5) sledi, da se ˇstevilo okuˇzenih ob ˇcasu t izraˇza po formuli I(t) = I(0) +S(0)−S(t) + γN

β ln S(t) S(0).

Iz enaˇcbe (2) vidimo, da bo pri maksimalni vrednosti I_max imel S vrednost ^γN_β . To sledi iz tega, da bo globalni maksimum funkcije I doseˇzen, ko bo odvod ^dI_dt enak 0. To pa bo natanko tedaj, ko bo S = ^γN_β . Torej, ˇce zdruˇzimo vse poznane zveze, dobimo

I_max=I(0) +S(0) + γN β

(︂

ln(︂ γN βS(0)

)︂−1)︂

.

V naˇsem primeru nam sama simulacija (pri koraku h= 0.1) vrne, da je bilo najveˇc okuˇzenih 2017.36 ljudi, naˇs izraˇcun pa nam vrne vrednost 2017.38. Opazimo lahko, da je priˇslo do razlike samo na drugem decimalnem mestu, to pa je posledica nu- meriˇcnega raˇcunanja - ˇce bi pri simulaciji uporabili manjˇsi korak h, bi bila napaka manjˇsa.

3.1. Razˇsirjen SIR model. Do sedaj smo obravnavali zelo preprost model, kjer sta bila tako β kot tudi γ konstanti. Torej je bilo konstantno tudi reprodukcijsko ˇstevilo. Sedaj pa si poglejmo ˇse primer, kako bi epidemija potekala, ˇce bi se to reprodukcijsko ˇstevilo spreminjalo, ob predpostavki, da je ˇcas kuˇznosti ˇse vedno konstanten. Denimo, da bi epidemija potekala takole. Epidemija se zaˇcne ˇsiriti in R₀ ima na zaˇcetku vrednost 4. Nato se 20. dan epidemije sprejmejo razni ukrepi za zajezitev ˇsirjenja (obvezno noˇsenje mask, prepoved zbiranja, pouk na daljavo itd.) in s tem se reprodukcijsko ˇstevilo R0 zniˇza na 0.8. Epidemija se zaˇcne postopno umirjati. Nato pa se 50. dan epidemije ukrepi nekoliko omilijo in R0 se dvigne na vrednost 2. ˇStevilo okuˇzenih ponovno zaˇcne rasti, ampak ne tako intenzivno kot na zaˇcetku, kar lahko vidimo na sliki 3, ki prikazuje grafiˇcni prikaz ˇsirjenja te epidemije.

Slika 3. Potek epidemije z ukrepi za zajezitev ˇsirjenja.

Ce primerjamo sliki 2 in 3, lahko opazimo, da sta na sliki 3 narisana 2 valova inˇ tudi dolˇzina same epidemije se je podaljˇsala, vendar pa smo s tem pristopom dosegli eno zelo pomembno zadevo. Maksimalno ˇstevilo okuˇzenih smo zmanjˇsali iz okoli 2000 na zgolj 500, kar pomeni, da bi bilo tudi manj ljudi hospitaliziranih in bi zato laˇzje poskrbeli za njih.

Iz statistike o COVID-19 lahko vidimo, da je pri premagovanju te epidemije veljalo kar nekaj ukrepov, npr. imeli smo pouk na daljavo, zaprto je bilo vse, razen kritiˇcnih infrastruktur itd. Na zgornjem primeru pa smo videli, kaj povzorˇcajo spremembe ukrepov v smislu poteka same epidemije in bi bilo zato smiselno, da se pri razvijanju modelov, ki bi nam pomagali pri napovedovanju ˇsirjenja neke doloˇcene epidemije, v

modele vkljuˇci tudi znaˇcilnosti epidemije. Pri COVID-19 je ena izmed teh znaˇcilnosti definitivno hitro ˇsirjenje, kar zahteva veliko ˇstevilo ukrepov. Prav tako pa bi bilo mogoˇce zanimivo pogledati ˇse kakˇsno drugo stanje. Da torej ne bi ljudi delili samo na dovzetne, okuˇzene in ozdravljene, ampak bi dodali ˇse kakˇsno stanje. In ravno to so naredili raziskovalci, ki so razvili θ-SEIHRD model ([3]), ki bo opisan v naslednjem poglavju.

4. θ-SEIHRD model

Poglejmo si torej novi model, ki so ga razvili ˇspanski raziskovalci in z njegovo pomoˇcjo napovedovali potek ˇsirjenja bolezni COVID-19 na Kitajskem za obdobje med decembrom 2019 in marcem 2020 ([3]). Gre za novi θ-SEIHRD model, ki je izpeljan iz Be-CoDiS modela ([4]), s katerim se lahko opisuje ˇsirjenje neke bolezni po celem svetu. Be-CoDiS model je bil na primer uporabljen tudi za opisovanje ˇsirjenja Ebole med letoma 2014 in 2016 in za ˇsirjenje Ebole v Demokratiˇcni republiki Kongo med letoma 2018 in 2020. V obeh primerih je podajal zelo toˇcne napovedi ([3]).

Najprej si bomo ogledali matematiˇcni zapis modela in kakˇsen sistem diferencial- nih enaˇcb moramo reˇsiti. Nato bomo opisali, kakˇsne rezultate vraˇca model, nekaj besed bomo nato namenili izboru parametrov, na koncu pa bomo sam model ˇse konfigurirali za neke realne podatke in komentirali dobljene rezultate.

4.1. Predpostavke modela. Model predpostavlja, da je vsaka oseba v enem izmed naslednjih stanj:

• Obˇcutljiva (oznaˇcimo sS): Oseba ni okuˇzena z virusom.

• Izpostavljena (oznaˇcimo zE): Oseba je v obdobju inkubacije virusa in nima nobenih vidnih kliniˇcnih znakov okuˇzbe. Taka oseba lahko okuˇzi druge, vendar z manjˇso verjetnostjo kot oseba, ki je ˇze kuˇzna oziroma nalezljiva - to postane po inkubacijski dobi.

• Okuˇzena (oznaˇcimo z I): Po inkubacijski dobi oseba postane nalezljiva, to je prvo stanje bolezni. Naslednje stanje je, da oseba poˇcasi zaˇcne razvijati kliniˇcne znake okuˇzbe. Iz tega stanja gre lahko v tri druga stanja - I_u, H_R inH_D.

• Nalezljiva ampak nezaznana (oznaˇcimo z I_u): To lahko postane oseba po stanjuI. Oseba ˇse vedno lahko okuˇzi druge, ima blage kliniˇcne znake okuˇzbe, vendar je njen test negativen oziroma ga sploh ni opravila. Model predvideva, da so v tem stanju lahko samo osebe z blagimi simptomi, ne pa tudi tiste ki bodo zaradi simptomov umrle.

• Ozdravljena nezaznana oseba (oznaˇcimo z Ru): Po stanju Iu oseba ozdravi in gre zato v stanjeRu. V tem stanju ne more veˇc okuˇziti drugih, razvije pa tudi naravno imunost na virus.

• Hospitalizirana ali v domaˇci karanteni, ki bo ozdravela (oznaˇcimo zH_R): Sem uvrstimo osebe iz skupine I, ki imajo pozitivni test in bodo ozdravele. Tudi v tem stanju ˇse vedno lahko okuˇzijo druge. ˇCas hospitalizacije oznaˇcujemo s C₀, ki je merjen v dneh in predstavlja povpreˇcni ˇcas hospitalizacije.

• Ozdravljene zaznane osebe (oznaˇcimo z R_d): Po stanju H_R oseba pride v to stanje. V tem stanju ne more veˇc okuˇziti drugih, razvije pa tudi naravno imunost na virus.

• Hospitalizirana, ki bo umrla (oznaˇcimo z H_D): Sem pridejo osebe iz I, ki bodo izgubile boj z virusom.

• Umrla zaradi COVID-19 (oznaˇcimo zD).

Same predpostavke se najlaˇzje grafiˇcno predstavi na naˇcin, kot je prikazan na sliki 4. Model pa dopuˇsˇca tudi vkljuˇcitev nekaterih ukrepov za zajezitev ˇsirjenja bolezni, kot so na primer:

• Izolacija: Okuˇzeni ljudje nimajo kontakta z drugimi ljudmi, ampak samo z zdravstvenim osebjem, kateri uporabljajo ustrezno zaˇsˇcitno opremo.

• Karantena: Omejitev gibanja ljudi, da se omeji moˇzne stike z okuˇzenimi osebami.

• Sledenje: Uporaba naˇcinov za prepoznavanje potencialnih okuˇzenih oseb.

• Poveˇcanje ˇstevila zdravstvenih delavcev in postelj za bolnike.

Na sliki 4 lahko vidimo tudi, na katera stanja ukrepi vplivajo.

S E I HR Rd

I_u R_u

H_D D

Ukrepi

Slika 4. Grafiˇcni prikaz predpostavk modela.

θ-SEIHRD model je uporabljen za opisovanje ˇsirjenja bolezni znotraj doloˇcenih teritorijev oziroma obmoˇcij v nekem fiksnem ˇcasovnem intervalu. Na zaˇcetku mora uporabnik vnesti nekaj zaˇcetnih pogojev. Kakˇsni so ti pogoji, se bo videlo ˇze v naslednjem razdelku. Omenimo pa lahko ˇse to, da model v naˇsem primeru privzema za t = 0 datum 1. 12. 2019 ob 10 uri dopoldan (ob tej uri namreˇc WHO objavlja posodobljene podatke), ki je najverjetnejˇsi datum prve okuˇzbe. Simulacijo pa seveda lahko zaˇcnemo na kateremkoli t₀ ≥ 0. Model dalje privzema, da v neokuˇzenih obmoˇcjih ˇzivijo samo osebe, ki so v stanjuS, torej osebe, ki so dovzetne za okuˇzbo.

V okuˇzenih obmoˇcjih pa so osebe v stanjih S, E, I, I_u, H_R, H_D, R_d, R_u, D in nato prehajajo med temi stanji v ˇcasovnem intervalu [t₀, t₀+T_max], kjer jeT_max∈N in predstavlja najveˇcje ˇstevilo dni simulacije. Simulacija se ustavi, ˇce pridemo do najveˇcjega ˇstevila dni simulacije oziroma ˇce je ob koncu dneva t v vsakem stanju E, I, Iu, HR, HD manj kakor 1 oseba v vseh teritorijih. Prav tako pa ima uporabnik tudi moˇznost, da zaˇcene simulacijo z ukrepi za zajezitev ˇsirjenja virusa in na ta naˇcin preizkuˇsa uˇcinkovitost le teh. Sedaj pa ˇse definirajmo, kakˇsen sistem diferencialnih enaˇcb sestavlja model.

4.2. Matematiˇcen zapis modela. Privzemimo, da so osebe v doloˇcenem teritoriju samo v enem izmed stanj opisanih v razdelku 4.1: S, E, I, I_u, H_R, H_D, R_d, R_u, D.

Nato pa ˇse privzemimo, da so novo rojene osebe v stanju S. Z naslednjim sistemom

diferencialnih enaˇcb opiˇsemo spreminjanje vrednosti zgornjih stanj:

dS⁽ⁱ⁾

dt (t) = S⁽ⁱ⁾(t) N⁽ⁱ⁾

(︂

m⁽ⁱ⁾_E (t)β_E⁽ⁱ⁾E⁽ⁱ⁾(t) +m⁽ⁱ⁾_I (t)β_I⁽ⁱ⁾I⁽ⁱ⁾(t) +m⁽ⁱ⁾_I

u(t)β_I⁽ⁱ⁾

u(θ⁽ⁱ⁾(t))I_u⁽ⁱ⁾(t))︂

− S⁽ⁱ⁾(t) N⁽ⁱ⁾

(︂

m⁽ⁱ⁾_H

R(t)β_H⁽ⁱ⁾

RH_R⁽ⁱ⁾(t) +m⁽ⁱ⁾_H

D(t)β_H⁽ⁱ⁾

DH_D⁽ⁱ⁾(t))︂

−µ⁽ⁱ⁾_mS⁽ⁱ⁾(t) +µ⁽ⁱ⁾_n (︂

S⁽ⁱ⁾(t) +E⁽ⁱ⁾(t) +I⁽ⁱ⁾(t) +I_i⁽ⁱ⁾(t) +R⁽ⁱ⁾_d (t) +R⁽ⁱ⁾_u (t) )︂

, dE⁽ⁱ⁾

dt (t) = S⁽ⁱ⁾(t) N⁽ⁱ⁾

(︂

m⁽ⁱ⁾_E (t)β_E⁽ⁱ⁾E⁽ⁱ⁾(t) +m⁽ⁱ⁾_I (t)β_I⁽ⁱ⁾I⁽ⁱ⁾(t) +m⁽ⁱ⁾_Iu(t)β_I⁽ⁱ⁾_u(θ⁽ⁱ⁾(t))I_u⁽ⁱ⁾(t) )︂

+ S⁽ⁱ⁾(t) N⁽ⁱ⁾

(︂

m⁽ⁱ⁾_H

R(t)β_H⁽ⁱ⁾

RH_R⁽ⁱ⁾(t) +m⁽ⁱ⁾_H

D(t)β_H⁽ⁱ⁾

DH_D⁽ⁱ⁾(t))︂

−µ⁽ⁱ⁾_mE⁽ⁱ⁾(t)−γ_EE⁽ⁱ⁾(t) +τ₁⁽ⁱ⁾(t)−τ₂⁽ⁱ⁾(t), dI⁽ⁱ⁾

dt (t) = γ_EE⁽ⁱ⁾(t)−(µ⁽ⁱ⁾_m +γ_I⁽ⁱ⁾(t))I⁽ⁱ⁾(t), dIu⁽ⁱ⁾

dt (t) = (1−θ⁽ⁱ⁾(t))γ_I⁽ⁱ⁾(t)I⁽ⁱ⁾(t)−(µ⁽ⁱ⁾_m +γ_I⁽ⁱ⁾

u(t))I_u⁽ⁱ⁾(t), dH_R⁽ⁱ⁾

dt (t) = θ⁽ⁱ⁾(t)(︂

1− ω⁽ⁱ⁾(t) θ⁽ⁱ⁾(t)

)︂

γ_I⁽ⁱ⁾(t)I⁽ⁱ⁾(t)−γ_H⁽ⁱ⁾

R(t)H_R⁽ⁱ⁾(t), dH_D⁽ⁱ⁾

dt (t) = ω⁽ⁱ⁾(t)γ_I⁽ⁱ⁾(t)I⁽ⁱ⁾(t)−γ_H⁽ⁱ⁾

D(t)H_D⁽ⁱ⁾(t), dR⁽ⁱ⁾_d

dt (t) = γ_H⁽ⁱ⁾

R(t)H_R⁽ⁱ⁾(t)−µ⁽ⁱ⁾_mR⁽ⁱ⁾_d (t), dR⁽ⁱ⁾u

dt (t) = γ_I⁽ⁱ⁾_u(t)I_u⁽ⁱ⁾(t)−µ⁽ⁱ⁾_mR⁽ⁱ⁾_u (t), dD⁽ⁱ⁾

dt (t) = γ_H⁽ⁱ⁾

D(t)H_D⁽ⁱ⁾(t).

Razloˇzimo, kaj pomenijo oznake, ki so vkljuˇcene v model, a ne predstavljajo ne- znank.

• Indeks i∈ {1, . . . , N_C}, kjer je N_C ∈Nˇstevilo drˇzav/teritorijev/obmoˇcij.

• N⁽ⁱ⁾ je ˇstevilo ljudi v teritorijui pred zaˇcetkom pandemije.

• µ⁽ⁱ⁾n ∈[0, 1] predstavlja nataliteto, merjeno v enotah (dan⁻¹).

• µ⁽ⁱ⁾m ∈[0, 1] predstavlja mortaliteto, merjeno v enotah (dan⁻¹).

• ω⁽ⁱ⁾(t)∈[ω⁽ⁱ⁾, ω⁽ⁱ⁾]⊂[0, 1] je deleˇz umrlih oseb zaradi virusa, kjerω⁽ⁱ⁾inω⁽ⁱ⁾ predstavljata oceni za minimalni in maksimalni deleˇz umrlih v posameznem teritoriju.

• θ⁽ⁱ⁾(t)∈[ω⁽ⁱ⁾,1] je deleˇz okuˇzenih ljudi, ki so evidentirani v teritorijuiv ˇcasu t. Za manjˇso poenostavitev samega modela lahko privzamemo, da so vse osebe, ki so umrle zaradi COVID-19, bile evidentirane. Zato lahko reˇcemo, da jeθ⁽ⁱ⁾(t)≥ω⁽ⁱ⁾. Ocene spremenljivkeθ za pandemijo COVID-19 lahko za nekatere drˇzave (vkljuˇcno s Kitajsko) najdemo v literaturi [3].

• β_E⁽ⁱ⁾, β_I⁽ⁱ⁾, β_I⁽ⁱ⁾

u, β_H⁽ⁱ⁾

R, β_H⁽ⁱ⁾

D ∈ R⁺ so parametri, ki nam povedo, koliko ljudi oseba iz doloˇcenega stanja okuˇzi na dan. Merjeni so v enotah (dan⁻¹).

• β_I⁽ⁱ⁾_u(θ) ∈ R⁺ je parameter, ki pove, koliko ljudi okuˇzi oseba iz stanja I_u na dan, pri pogoju, da je deleˇz registrirano okuˇzenih enak θ.

• γ_E ∈(0,∞) je stopnja prehoda iz stanja E v stanje I, ki pa je enaka za vse teritorije. Merjen pa je v enotah (dan⁻¹).

• γ_I⁽ⁱ⁾(t)∈(0,∞) je stopnja prehoda iz stanjaI v stanja I_u, H_R, H_D v terito- riju iob ˇcasu t.

• γ_I⁽ⁱ⁾

u(t), γ_H⁽ⁱ⁾

R(t), γ_H⁽ⁱ⁾

D(t)∈(0,∞) oznaˇcujejo stopnjo prehoda iz stanj I_u, H_R, H_D v stanjaR_u, R_d, D v teritoriju i ob ˇcasu t.

• m⁽ⁱ⁾_E (t), m⁽ⁱ⁾_I (t), m⁽ⁱ⁾_I

u(t), m⁽ⁱ⁾_H

R(t), m⁽ⁱ⁾_H

D(t)∈[0, 1] so funkcije, ki predstavljajo uˇcinkovitost ukrepov za zajezitev ˇsirjenja bolezni, glede na doloˇceno stanje.

Merjene so v procentih in so doloˇcene posebej za vsak teritorij iin ˇcast.

• τ₁⁽ⁱ⁾ predstavlja ˇstevilo okuˇzenih oseb, ki pridejo v teritorij iiz drugih terito- rijev v enem dnevu.

• τ₂⁽ⁱ⁾ predstavlja ˇstevilo okuˇzenih oseb, ki zapustijo teritoriji v enem dnevu.

Zaˇcetni pogoji danega sistema diferencialnih enaˇcb pa so dani z

S⁽ⁱ⁾(t₀), E⁽ⁱ⁾(t₀), I⁽ⁱ⁾(t₀), I_u⁽ⁱ⁾(t₀), H_R⁽ⁱ⁾(t₀), H_D⁽ⁱ⁾(t₀), R⁽ⁱ⁾_d (t₀), R⁽ⁱ⁾_u (t₀), D⁽ⁱ⁾(t₀) Ce torej sliki 4 dodamo ˇse te parametre, ki jih doloˇˇ cimo pri matematiˇcnemu zapisu modela, lahko na sliki 5 vidimo grafiˇcen prikaz tega matematiˇcnega zapisa modela.

S E I HR Rd

I_u R_u

H_D D

Ukrepi

β⁽ⁱ⁾_E +β_I⁽ⁱ⁾+β⁽ⁱ⁾_I

u+ β⁽ⁱ⁾_H

R+β_H⁽ⁱ⁾

D

γE(i) θ⁽ⁱ⁾(1−^ω_θ(i)⁽ⁱ⁾)γ_I⁽ⁱ⁾ γ_H⁽ⁱ⁾

R

γ_I⁽ⁱ⁾

u

γ_H⁽ⁱ⁾

D

(1−θ⁽ⁱ⁾)γ_I⁽ⁱ⁾

θ^(i)ω_θ(i)⁽ⁱ⁾γ_I⁽ⁱ⁾

Slika 5. Grafiˇcni prikaz matematiˇcnega zapisa modela.

Opomba 4.1. Opazimo lahko, da 9. enaˇcba iz zgornjega sistema diferencialnih enaˇcb ni povezana z ostalimi enaˇcbami in jo lahko zato izraˇcunamo posebej. Tako dobimo

D⁽ⁱ⁾(t) =D⁽ⁱ⁾(t₀) +

∫︂ t t0

γ_H⁽ⁱ⁾

D(s)H_D⁽ⁱ⁾(s)ds.

(6)

Opomba 4.2. Zaradi manjˇse poenostavitve samega modela, bomo v tem diplom- skem delu obravnavali samo en terirorij, torej i= 1, in ga zato ne bomo pisali, prav tako pa ne bomo upoˇstevali natalitete in mortalitete, saj se izkaˇze, da ti dve koliˇcini (vsaj v krajˇsem ˇcasovnem obdobju) nimata velikega vpliva. Privzemimo torej, da je µn=µm = 0.

Z upoˇstevanjem prejˇsnje opombe dobimo naslednji sistem diferencialnih enaˇcb:

dS

dt(t) = S(t) N

(︂

m_E(t)β_EE(t) +m_I(t)β_II(t) +m_Iu(t)β_Iu(θ(t))I_u(t))︂

− S(t) N

(︂

mH_R(t)βH_RHR(t) +mH_D(t)βH_DHD(t) )︂

, dE

dt (t) = S(t) N

(︂

m_E(t)β_EE(t) +m_I(t)β_II(t) +m_Iu(t)β_Iu(θ(t))I_u(t))︂

+ S(t) N

(︂

m_HR(t)β_HRH_R(t) +m_HD(t)β_HDH_D(t))︂

−γ_EE(t) +τ₁(t)−τ₂(t), dI

dt(t) =γ_EE(t)−γ_I(t)I(t), dI_u

dt (t) = (1−θ(t))γ_I(t)I(t)−γ_Iu(t)I_u(t), dHR

dt (t) =θ(t) (︂

1− ω(t) θ(t) )︂

γI(t)I(t)−γH_R(t)HR(t), dHD

dt (t) =ω(t)γ_I(t)I(t)−γ_HD(t)H_D(t), dR_d

dt (t) =γ_HR(t)H_R(t), dR_u

dt (t) =γ_Iu(t)I_u(t), dD

dt (t) =γHD(t)HD(t).

Opomba 4.3. Tudi tukaj lahko opazimo, da enaˇcbe na 7., 8.in 9.mestu niso pove- zane z ostalimi enaˇcbami iz sistema in jih zato lahko izraˇcunamo posebej. Dobimo:

Rd(t) =Rd(t0) +

∫︂ t t0

γHR(s)HR(s)ds, R_u(t) =R_u(t₀) +

∫︂ t t0

γ_Iu(s)I_u(s)ds, D(t) =D(t₀) +

∫︂ t t0

γ_HD(s)H_D(s)ds.

Ostale neznane funkcije iz sistema diferencialnih enaˇcb pa bomo izraˇcunali s pomoˇcjo Runga-Kutta metode stopnje 4.

Za laˇzjo predstavo pa si tudi tukaj poglejmo, kako izgleda sam poenostavljen model in kakˇsne prehode med stanji dopuˇsˇca. Grafiˇcni prikaz vidimo na sliki 6.

S E I H_R R_d

I_u R_u

H_D D

β_E+β_I+β_Iu+ β_HR+β_HD

γ_E θ(1−^ω_θ)γI γHR

γ_Iu

γH_D

(1−θ)γ_I

θ^ω_θγ_I

Slika 6. Grafiˇcni prikaz modela, ki je posledica opombe 4.2.

4.3. Rezultati modela. Sedaj bomo opisali nekaj koliˇcin, ki jih lahko dobimo, ko enkrat izraˇcunamo reˇsitev sistema diferencialnih enaˇcb in zato te koliˇcine imenujemo kar rezultati modela.

• c_m(t): predstavlja komulativno ˇstevilo okuˇzenih z virusom COVID-19 v ˇcasu t. Podan je kot

cm(t) =HR(t) +HD(t) +Rd(t) +D(t), oziroma ga lahko izraˇcunamo tudi takole:

c_m(t) =c_m(t₀) +

∫︂ t t0

d(H_R+H_D +R_d+D)

dt (s)ds

=c_m(t₀) +

∫︂ t t0

θ(s)(︂

1−ω(s) θ(s)

)︂

γ_I(s)I(s)−γ_HR(s)H_R(s) +ω(s)γ_I(s)I(s)−γ_HD(s)H_D(s) +γ_HR(s)H_R(s) +γ_HD(s)H_D(s)ds

=c_m(t₀) +

∫︂ t t0

θ(s)γ_I(s)I(s)ds.

• d_m(t): predstavlja komulativno ˇstevilo umrlih zaradi COVID-19 v ˇcasu t, podan pa je z D(t).

• R₀ in R_e predstavljata reprodukcijski ˇstevili, torej koliko oseb v povpreˇcju okuˇzi okuˇzena oseba v svojem kuˇznem obdobju. Pri tem R₀ predstavlja t. i. bazno reprodukcijsko ˇstevilo, ko ne velja noben ukrep, odvisen pa je od sestave populacije in se ne spreminja skozi samo epidemijo. Oznaka R_e pa predstavlja efektivno reprodukcijsko ˇstevilo, ki se spreminja skozi samo simulacijo, saj je zelo odvisno od tega, koliko ljudi je ˇse v stanju S. Velja R_e(0) = R₀. Iz vira [3] na strani 7 lahko vidimo, da se koliˇcini izraˇzata po formulah

R0 = βIU(1−θ)

γ_Iu + βHR(θ−ω) γ_HR + βI

γ_I +βE

γ_E + βHDω γ_HD , (7)

pri ˇcemer so vse koliˇcine, ki so odvisne od ˇcasa, izvrednotene za ˇcas t = t₀, in

R_e= mIuβI_U(1−θ)

γ_Iu +mH_RβH_R(θ−ω)

γ_HR +mIβI

γ_I +mEβE

γ_E +mH_DβH_Dω γ_HD .

• Hos(t): predstavlja ˇstevilo hospitaliziranih ljudi. Podan pa je kot Hos(t) =H_D(t) +p(t)(H_R(t) + (R_d(t)−R_d(t−C₀))),

kjerp(t) predstavlja deleˇz hospitaliziranih oseb v ˇcasut, ˇsteviloC0pa je para- meter, ki nam pove povpreˇcen ˇcas hospitalizacije. Poznavanje te funkcije nam omogoˇca, da lahko napovedujemo ˇstevilo potrebnih postelj v bolniˇsnicah.

• MHos predstavlja maksimalno ˇstevilo hospitaliziranih oseb v obdobju simu- lacije. Izraˇcuna se ga kot MHos = max

t∈[t0,T]Hos(t).

• Γ_E(t), Γ_Iu(t) in Γ_H(t) predstavljajo ˇstevilo ljudi okuˇzenih v obdobju simu- lacije zaradi kontaktov z ljudmi iz stanj E, I_u in H = H_R+H_D. Izraˇcuna se jih po formulah ([3])

Γ_E(t) =

∫︂ t t0

m_E(s)β_EE(s)S(s) N ds, Γ_Iu(t) =

∫︂ t t0

m_Iu(s)β_IuI_u(s)S(s) N ds, ΓH(t) =

∫︂ t t0

(︂

mH_R(s)βH_RHR(s) +mH_D(s)βH_DHD(s) )︂S(s)

N ds.

4.4. Ocenjevanje parametrov. Da bodo rezultati modela dovolj natanˇcno napo- vedovali potek epidemije je potrebno, da pravilno in dovolj natanˇcno nastavimo parametre, ki nastopajo v modelu. Ravno ta del pa je eden teˇzjih problemov takega modeliranja in bo obravnavan v tem poglavju. Nekaj parametrov, ki se jih uporabi v simulacijah, je bilo vzetih iz literature. Ker pa je o virusu SARS-Cov-19 ˇse marsikaj neznanega, pa je bilo potrebno nekatere parametre oceniti na novo z empiriˇcnimi metodami ali pa s kakˇsnimi drugaˇcnimi pripomoˇcki za ocenjevanje parametrov.

4.4.1. Parametri ukrepov zoper ˇsirjenja COVID-19. Pri samem modelu lahko opa- zimo, da parametre, ki oznaˇcujejo, koliko oseb okuˇzi ena okuˇzena oseba v doloˇcenem stanju (torej β_E, β_I, β_Iu, β_HR, β_HD), mnoˇzimo s funkcijami, ki jim pravimo para- metri ukrepov zoper ˇsirjenja COVID-19. V naˇsem primeru jih izberemo na naslednji naˇcin ([3]):

m_E(t) = m_Iu(t) =m_I(t) = m_HR(t) = m_HD(t) =

=

⎧

⎪⎪

⎪⎪

⎪⎪

⎨

⎪⎪

⎪⎪

⎪⎪

⎩

(m₀−m₁) exp(︂

−κ₁(t−λ₀))︂

+m₁, t∈[t₀, λ₁], (m(λ₁)−m₂) exp(︂

−κ₂(t−λ₁))︂

+m₂, t∈(λ₁, λ₂], ...

(m(λq−1)−mq) exp (︂

−κq(t−λq−1) )︂

+mq, t∈(λq−1, ∞), kjer za vsakj ∈ {0, . . . , q}ˇstevilom_j ∈[0,1] predstavlja mero intenzivnosti kontrol- nih ukrepov (veˇcja kot je ˇstevilka, manjˇse je ˇstevilo novo okuˇzenih oseb),κ_j ∈[0,∞) predstavlja mero uˇcinkovitosti kontrolnih ukrepov (veˇcja kot je ˇstevilka, manjˇse je ˇstevilo novo okuˇzenih oseb), λj ∈ [t0, ∞) pa predstavlja prvi dan uvedbe novih

ukrepov. ˇStevilo λ₀ ∈ [0,∞) predstavlja uveljavitev ukrepa pred ˇcasom t₀, ˇce je seveda takrat ˇze veljal kakˇsen od ukrepov, ˇsteviloq ∈N pa predstavlja ˇstevilo vseh ukrepov oziroma njihovih sprememb v danem ˇcasovnem obdobju.

Vrednostiλ_j so v praksi znane, saj se lahko poiˇsˇce podatke, kdaj je kakˇsen teritorij oziroma drˇzava uvedla doloˇcen ukrep, prav tako lahko ocenimo vrednostimj, ostale parametre pa je potrebno oceniti na bolj kompliciran naˇcin ki bo predstavljen v razdelku 4.4.5.

4.4.2. Deleˇz umrlih oseb. Deleˇz umrlih oseb ω(t) je odvisen od ˇcasa t in podaja razmerje med umrlimi in skupnim ˇstevilom okuˇzenih (tako zaznanih okuˇzenih kot tudi nezaznanih). Na to razmerje lahko zelo vplivajo doloˇceni ukrepi (npr. zgodnje prepoznavanje okuˇzb, boljˇsi pogoji v zdravstvu itd.). V naˇsem primeru definiramo

ω(t) =mI(t)ω+ (1−mI(t))ω,

kjer je ω ∈ [0, 1] deleˇz umrlih oseb, ko ne velja noben ukrep (takrat je funkcija m_I(t) = 1), ω ∈ [0, 1] pa je deleˇz umrlih oseb, ko veljajo vsi ukrepi. Konstantni ω inω je potrebno ˇse oceniti, kar pa bo opisano v razdelku 4.4.5. Ker velja ω ≥ω, lahko zapiˇsemo ω =ω+δω, kjer jeδω ∈[0, 1−ω].

4.4.3. Stopnje prehodov iz enih stanj v druga. Oznaˇcimo z d_E, d_I, d_Iu, d_HR in d_HD povpreˇcen ˇcas, ki ga neka oseba preˇzivi v stanjuE, I, I_u, H_RinH_D brez upoˇstevanja ukrepov. Upoˇstevajmo pa ˇse:

• Prehod iz stanjaEv stanjeI je odvisen le od tipa bolezni, zato je konstanten.

• Na parameter γ_I(t) je moˇzno vplivati z doloˇcenimi ukrepi (npr. zgodnejˇse testiranje). Iz tega pa sledi, da je moˇzno vplivati tudi na parametre γ_Iu(t), γ_HR(t) in γ_HD(t), saj so potem osebe dlje ˇcasa v doloˇcenem stanju.

• Vrednostd_HR je niˇzja od vrednostid_HD, zato lahko zapiˇsemod_HD =d_HR+δ_R, kjer je δ_R>0.

• Zaradi same poenostavitve privzemimo ˇse, da je povpreˇcen ˇcas ozdravitve nezaznane okuˇzene osebe isti kot ˇcas ozdravitve zaznane osebe, zato jedIu = dHR. Privzemimo pa ˇse, da je ˇcas okrevanja konstanten.

Z upoˇstevanjem zgornjih predpostavk dobimo naslednje enaˇcbe, s katerimi ocenju- jemo vrednosti teh parametrov:

γ_I(t) = 1

d_I−g(t), γ_Iu(t) =γ_HR(t) = 1

d_I+g(t), γ_HD(t) = 1

d_Iu+g(t) +δ_R. Pri temg(t) =d_g(1−m_I(t)) predstavlja zmanjˇsanje ˇcasad_Izaradi doloˇcenih ukrepov v ˇcasu t, d_g pa predstavlja maksimalno ˇstevilo, za katero lahkod_I zmanjˇsamo.

4.4.4. Deleˇzi kontaktov okuˇzenih oseb. Od teh parametrov je najteˇzje doloˇciti pa- rameter, ki nam opisuje deleˇz kontaktov oseb v stanju I_u, zato bomo ta parameter poiskuˇsali oceniti s formulo

β_Iu(θ) =

⎧

⎪⎨

⎪⎩

β_I, ˇceθ =ω,

neka nenaraˇsˇcajoˇca funkcija, ˇceθ ∈(ω, 1),

βI, ˇceθ = 1,

kjer staβ_I inβ_I spodnja in zgornja meja. Zaradi same poenostavitve pa privzemimo, da velja enakost β_I = β_I. Prav tako pa privzemimo ˇse, da obstaja zveza med β_I, β_E, β_HR, β_HD in β

I. Dalje privzemimo ˇse, da so osebe v stanjih E, H_R, H_D

inI_u manj nalezljive kot osebe v stanjuI (zaradi manjˇse koncentracije virusa ali pa zaradi samoizolacije). Torej lahko parametre izraˇcunamo po formuli

β_E =C_Eβ_I, β_HR(t) = β_HD(t) =C_H(t)β_I inβ

I =C_uβ_I,

kjer so C_E, C_H(t) in C_u ∈ [0, 1] ([3]). V razdelku 4.5 bomo pokazali in natanˇcno razloˇzili, kako se doloˇci te parametre za primer Kitajske.

4.4.5. Ocenjevanje parametrov s pomoˇcjo multi-objektivne metode. Naj bo Ω druˇzina vseh parametrov, ki jih ˇzelimo oceniti s pomoˇcjo multi-objektivne metode. Doloˇciti moramo multi-objektivni optimizacijski problem, ki temelji na minimizaciji relativ- nih napak med registriranim komulativnim ˇstevilom okuˇzenih in mrtvih (oznaˇcimo jih z c_r in d_r) in komulativnim ˇstevilom okuˇzenih in mrtvih, ki jih vrne model ob uporabi parametrov iz Ω (oznaˇcimo jih z c^Ω_m in d^Ω_m). Ta optimizacijski problem zapiˇsemo na naslednji naˇcin

min{f_c(Ω), f_d(Ω)} pri pogoju Ω∈S ⊂Rⁿ,

kjer jeS mnoˇzica, ki vkljuˇcuje zgornje in spodnje meje ocenjevanih parametrov,nje ˇstevilo parametrov v druˇzini Ω, funkcijif_c inf_d pa sta definirani na naslednji naˇcin

f_c(Ω) = ∥c_r−c^Ω_r∥_L²_{([t0, t0+Tmax])}

∥c_r∥_L²_{([t0, t0+Tmax])} , f_d(Ω) = ∥d_r−d^Ω_r∥_L²_{([t0, t0+Tmax])}

∥d_r∥_L²_{([t0, t0+Tmax])} , kjer ∥g∥_L²_{([t0, t0+Tmax])} = (︂

∫︁t0+Tmax

t0 g(t)²dt)︂1/2

oznaˇcuje L² normo funkcije g na intervalu [t₀, t₀+T_max].

Tak optimizacijski problem se lahko reˇsuje z WASF-GA algoritmom (Weighting Achievement Scalarizing Function Genetic Algorithm), ki spada med t. i. genet- ske algoritme. Ideja tega algoritma je, da na zaˇcetku iz prostora Ω izbere nekaj nakljuˇcnih n-teric parametrov (v naˇsem primeru moramo s to metodo doloˇciti 7 parametrov, torej izbere nekaj sedmeric teh parametrov). Nato za vsako n-terico (sedmerico) izraˇcuna model in preveri, za koliko je reˇsitev oddaljena od globalnega minimuma. Nato izbere tiste parametre, ki so najbliˇzji globalnemu minimumu. V nadaljevanju med njimi naredi razne mutacije in kriˇzanja, tako da dobi nove n- terice (sedmerice) parametrov, na katerih postopek ponovi. Zaradi same teˇzavnosti takega programiranja se s tem problemom ne bomo ukvarjali in bomo samo pri- vzeli podatke, ki so podani v viru [3] in so predstavljeni v tabeli 1. Pripomnimo lahko samo ˇse to, da datum 29. 3. pomeni, da smo model prilagajali podatkom do vkljuˇcno 29. 3., podobno datum 25. 2. pomeni, da smo model prilagajali podatkom do vkljuˇcno 25. 2., torej smo do tega datuma nastavljali to sedmerico parametrov itd.

4.5. Numeriˇcna simulacija in konfiguracija modela za primer Kitajske. V tem razdelku bomo naredili nekaj numeriˇcnih simulacij samega modela in predstavili, kako uˇcinkovit je v napovedovanju ˇsirjenja bolezni COVID-19. Na sliki 7 lahko vidimo, kako se je spreminjalo ˇstevilo okuˇzenih ljudi na Kitajskem med 12. 1. 2020 in 28. 3. 2020. V obdobju med 11. 2. in 16. 2. lahko opazimo zelo velik porast ˇstevila okuˇzenih. Razlog za to je, ker v samem zaˇcetku epidemije ˇse niso imeli toˇcno doloˇcenih strategij, kako bodo beleˇzili pozitivne osebe. Zato se je oblast neke province na Kitajskem odloˇcila, da ne bo poroˇcala samo tistih okuˇzenih, ki so pozitivni na testih, temveˇc tudi tiste, za katere so zdravniki menili, da imajo zelo verjetno COVID-19. Kasneje je WHO izdala malo bolj natanˇcna navodila, kako naj

Zapis 29. 3. 25. 2. 8. 2. 29. 1.

β_I 0.2834 0.3653 0.3718 0.4012 C_E 0.3806 0.0800 0.0828 1.4306·10⁻⁵ Cu 0.4001 0.4053 0.4271 0.8559 δR 7.0000 9.7528 7.2773 12.6250 δ_ω 0.0195 0.0155 1.3265·10⁻⁵ 4.6360·10⁻⁶

ω 0.0161 0.0220 0.0220 0.0220 κ₁ 0.1090 0.1007 0.1415 0.2000

Tabela 1. Parametri doloˇceni s pomoˇcjo multi-objektivne metode ([3]).

se beleˇzi okuˇzene osebe (naj se upoˇsteva zgolj tiste, ki so pozitivni na testih) in s tem poenotila kriterij, kdaj je neka oseba pozitivna.

Ker pa iz 4.4.5 vemo, da se nekateri parametri ocenjujejo glede na preteklo pot epidemije, bi taki nenadni porasti, kot so tukaj, lahko bistveno vplivali na same parametre. Zato bomo sedaj najprej malo preoblikovali te podatke. Teh 15,141 novo okuˇzenih bomo s pomoˇcjo naslednje formule razdelili na obdobje pred 13. 2. 2020 po formuli

c_ar(t) =c_r(t) + 15141

∑︁t

τ=tcr(τ)

∑︁ta

τ=tc_r(τ),

kjer je c_ar(t) popravljeno ˇstevilo novo okuˇzenih v ˇcasu t, c_r(t) je poroˇcano ˇstevilo okuˇzenih, t = 12. 1. 2020 in t_a = 13. 2. 2020. Na sliki 7 lahko vidimo, kako sta se c_ar in c_r spreminjala v opazovanem obdobju.

Slika 7. ˇStevilo okuˇzenih oseb na Kitajskem

Sedaj, ko smo popravili podatke, pa si poglejmo, kako bomo nastavili parametre.

V vseh simulacijah bomo uporabili naslednje podatke, ki so dobljeni iz virov [7], [6], [8]:

• N = 1,400,812,636 ljudi,

• C₀ = 14 dni,

• p(t)≡1,

• d_E = 5.5 dni, iz ˇcesar sledi, da je γ_E = 1/5.5 dni⁻¹,

• d_I = 6.7 dni,

• d_g = 6 dni,

• d_Iu = 14−d_I = 7.3 dni.

Za βIu pa privzemimo tole linearno zvezo [3]

β_Iu(θ) =β

I+ β_I−β

I

1−ω(t)(1−θ), (8)

pri ˇcemer je β

I = CuβI konstanta. Iz podatkov WHO je razvidno, da je bilo 20.

2. 2020 na Kitajskem z novim koronavirusom okuˇzenih 2055 zdravstvenih delavcev od skupaj 74,651 okuˇzenih. Iz teh podatkov lahko sklepamo, da je deleˇz okuˇzenih kontaktov v stanjih H_R in H_D pribliˇzno 100₇₄₆₅₁²⁰⁵⁵ % = 2.75 % od vseh okuˇzenih.

Ce se sedaj spomnimo, kaj nam pove reprodukcijsko ˇstevilo (koliko oseb okuˇˇ zi ena oseba v svojem kuˇznem obdobju) in ˇce to ˇstevilo pomnoˇzimo s ˇstevilom okuˇzenih in delimo s povpreˇcnim ˇcasom kuˇznosti, dobimo ˇstevilo, ki nam pove, koliko oseb je bilo okuˇzenih v enem dnevu. Razmerje, ki je podano iz podatkov WHO, je torej isto, kakor ˇce bi gledali razmerje reprodukcijskih ˇstevil izvrednotenih v t = 20.2.2020, kjer bi imeli v ˇstevcu reprodukcijsko ˇstevilo stanjH_RinH_D, v imenovalcu pa reprodukcijsko ˇstevilo vseh stanj. Reprodukcijsko ˇstevilo za stanji H_R in H_D se izraˇcuna kot

R_HR,HD = β_HR(1− ^ω_θ)θ γHR

+ β_HD(^ω_θ)θ γHD

. Ce sedaj upoˇstevamo ˇseˇ βHR =βHD =CHβI, dobimo

R_HR,HD =C_Hβ_Iθ(︂

(1−ω θ) 1

γ_HR + ω θ

1 γ_HD

)︂

. Ob upoˇstevanju enaˇcbe (7) pa dobimo zvezo

C_Hβ_Iθ(︂

(1− ^ω_θ)_γ¹

HR +^ω_θγ¹

HD

)︂

C_Hβ_Iθ(︂

(1−^ω_θ)_γ¹

HR + ^ω_θγ¹

HD

)︂

+ ^β_γ^I

I +^β_γ^E

E + (1−θ)^β_γ^Iu

Iu

= 0.0275, iz ˇcesar sledi, da je

C_H =

0.0275 (︂βI

γI +^β_γ^E

E + (1−θ)^β_γ^Iu

Iu

)︂

(1−0.0275)β_Iθ(︂

(1−^ω_θ)_γ¹

HR + ^ω_θγ¹

HD

)︂. (9)

Iz literature [3] nato ˇse razberemo, da je bilo do 23. januarja 2020 kar 86 % vseh primerov nezaznanih. Po 8. februarju pa je ta deleˇz padel na 35 %. Zato glede na te deleˇze doloˇcimo

θ(t) =

⎧

⎪⎨

⎪⎩

1−0.86 = 0.14, ˇcet⩽24. 1. 2020,

linearna funkcija, ˇce 24.1. 2020⩽t⩽8. 2.2020, 1−0.35 = 0.65, ˇce 8. 2. 2020⩽t.

(10)