Epidemie

Stochastyczne i deterministyczne modele epidemii

Zjawiska opisywane na poziomie "mikro" przez procesy Markowa zachowują się na poziomie "makro" jak funkcje deterministyczne, będące rozwiązaniami równań różniczkowych. Można to zaobserwować na przykładzie prostego (ale, niestety, bardzo aktualnego w chwili pisania tego artykułu) modelu początkowej, wykładniczej fazy epidemii.

Rozważmy bardzo dużą populację. Niech $| I(t)$ oznacza liczbę zarażonych w chwili $t.$ Zakładamy, że w krótkim okresie czasu $h$ każdy z nich zaraża średnio $|αh$ innych ludzi, ale przy tym z prawdopodobieństwem $|βh$ zdrowieje lub umiera i przestaje zarażać. Deterministyczny model jest następujący:

I(t+ h) = I(t) +(α − β)I(t)h + o(h) dla h 0.

(1)

(Traktujemy $I(t)$ jako wielkość ciągłą. Symbol $o(h)$ oznacza nieskończenie małą rzędu wyższego niż $h,$ czyli funkcję spełniającą równość $|limh 0o(h)/h = 0.$ ) Oczywiście (1) sprowadza się do zwyczajnego równania różniczkowego:

I(t) d-----= (α − β)I(t). d t

Popatrzmy na to samo zjawisko z bliska (w skali mikro). Potraktujmy $I(t)$ jako zmienną losową o wartościach całkowitoliczbowych. Ewolucja procesu $|I(t)$ jest opisana równaniami

P(I(t +h) = i + 1 I(t) = i) = αih+ o(h) dla h 0, P(I(t +h) = i− 1 I(t) = i) = βih+ o(h) dla h 0, P(I(t +h) = i I(t) = i) = 1− (α+ β)ih + o(h) dla h 0.

(2)

Zakładamy, że ${I(t) t⩾ 0}$ jest procesem Markowa, to znaczy ewolucja procesu po chwili $t$ zależy tylko od $|I(t),$ czyli stanu w chwili $t,$ a nie od wcześniejszej "historii" procesu. Jest to skokowy proces Markowa: w losowych momentach $T1 < T2 < ... < Tn < ...$ proces rośnie lub maleje o 1. Żeby lepiej zrozumieć strukturę tego procesu i sposób jego symulowania, przyjrzyjmy się czasom pomiędzy skokami, $Wn+1 = Tn+1− Tn,$ gdzie

Tn+1 = inf{t > Tn I(t) /= I(Tn)} (rys 1).

Okazuje się (co zostanie uzasadnione na końcu artykułu), że $Wn+1$ jest zmienną losową o rozkładzie wykładniczym z parametrem $(α +β )in$ . Ponadto

$pict$

Łatwo w to uwierzyć, patrząc na równania (2).

Przebieg procesu $| I(t)$ symulujemy zatem, losując kolejno czasy $| Tn$ i stany $|I(Tn)$ zgodnie z dwiema powyższymi regułami. Chemicy, biolodzy i epidemiolodzy znają tę metodę pod nazwą algorytmu Gillespie'go. W rzeczywistości taki opis skokowego procesu Markowa podał znakomity amerykański matematyk czeskiego pochodzenia, Joseph Doob. Przykładowy przebieg krótkich symulacji przedstawiamy w poniższej tabelce i na rysunku 1.

Popatrzmy teraz na dłuższe przebiegi symulacji i porównajmy trajektorie procesu losowego z rozwiązaniem równania różniczkowego, czyli funkcją wykładniczą $|I(t) = i0 exp{(α − β)t}.$ Przyjęliśmy następujące wartości parametrów: $α = 2,β= 1,i0 = I(0) = 1$ (śledzimy rozwój epidemii od pierwszego zarażonego). Na rysunku 2 przedstawionych jest 20 niezależnych trajektorii procesu Markowa. Spośród tych 20 trajektorii 12 wpadło dość wcześnie w stan 0. Jest to stan "pochłaniający": jeśli $I(Tn) = 0,$ to dla każdego $|t⩾ Tn$ mamy $I(t) = 0.$ Pogrubiona, kolorowa linia to wykres funkcji wykładniczej.

Uderzającą cechą procesu losowego opisanego równaniami (2) jest zgodność trajektorii w makroskali z rozwiązaniem równania różniczkowego, czyli z funkcjami wykładniczymi $I(t) = exp{(α − β)(t −c)}.$ Poszczególne realizacje procesu różnią się jednak wyraźnie przesunięciem w fazie $|c.$ To widać na rysunku 3, na którym oś pionowa została przedstawiona w skali logarytmicznej. Rysunek 4 pokazuje te same 20 trajektorii w początkowym okresie rozwoju. Wspomniane wyżej przesunięcia w fazie wynikają z losowych fluktuacji na początku. Rezultatem tych losowych fluktuacji jest to, że 12 spośród przedstawionych trajektorii jest "pochłoniętych" przez zero. Dalszy przebieg procesu jest niemal deterministyczny.

Model SIR

Rozważmy populację złożoną z $|ℓ$ osobników. Niech $I(t)$ oznacza liczbę zarażonych w chwili $t,$ zaś $|R(t)$ łączną liczbę uodpornionych i zmarłych. Liczba osobników narażonych na zakażenie jest równa $|S(t) = ℓ− I(t)− R(t).$ Osobnik typu $|S$ może przejść do kategorii $|I,$ a stąd do kategorii $R$ (stąd nazwa "model SIR"). Schematycznie:

|--|--α--|----β-- |-| -S-| -I-| R--

Zakładamy, że w krótkim okresie czasu $|h$ każdy osobnik $I$ zaraża średnio $α (S(t)/ℓ)h$ innych ludzi, a z prawdopodobieństwem $|βh$ zdrowieje lub umiera i przestaje zarażać. Zauważmy, że jeśli $|S(t)/ℓ≃ 1,$ to dostajemy model fazy wykładniczej, przedstawiony w zadaniu poprzednim.

Klasyczny model deterministyczny jest następujący:

S(t)- I(t +h) = I(t)+ (α ℓ − β) I(t)h + o(h) dlah 0, S(t +h) = S(t) −α S(t)-I(t)h+ o(h) dla h 0, ℓ

(3)

oraz $R(t) = ℓ− I(t)− S(t).$ Oczywiście (3) jest układem równań różniczkowych zwyczajnych:

d-I(t)- S(t)- dS(t)- S(t)- dt = (α ℓ − β)I(t), d t = −α ℓ I(t).

(4)

W skali mikro traktujemy $I(t)$ i $S(t)$ jako zmienne losowe o wartościach całkowitoliczbowych. Stan układu jest parą $(I,S)t = (I(t),S(t)) = (i,s).$ Ewolucja procesu jest opisana równaniami

s P((I,S)t+h = (i + 1,s− 1) (I,S)t = (i,s)) = α-ih+ o(h) dla h 0, ℓ P((I,S)t+h = (i− 1,s) (I,S)t = (i,s)) = β ih + o(h) dla h 0, s P((I,S)t+h = (i,s) (I,S)t = (i,s)) = 1− (α-+ β) ih + o(h) dla h 0. ℓ

(5)

Symulacja procesu opisanego powyższymi równaniami również może zostać przeprowadzona przy użyciu algorytmu Gillespie'go. Tym razem jednak musimy monitorować dwie wielkości - liczbę zainfekowanych $|(I)$ oraz liczbę osobników podatnych na zakażenie $(S)$ . Jeśli liczby te wynoszą $i$ i $s,$ odpowiednio, to czas oczekiwania na kolejne "zdarzenie" (czyli zainfekowanie nowej osoby lub wyzdrowienie/śmierć osoby zarażonej) ma rozkład wykładniczy z parametrem $α |ℓsi +β i.$ Kiedy już dojdzie do zdarzenia, to z prawdopodobieństwem $- | -`sisi+βi `$ jest to nowe zarażenie (oraz wyzdrowienie/śmierć w przeciwnym przypadku). Nowe stany odpowiadające tym dwóm rodzajom zdarzeń to, odpowiednio, $|(i + 1,s −1)$ oraz $(i− 1,s).$ Oczywiście symulacja kończy się, kiedy $|i = 0.$

Rys. 5 Trajektorie procesu $I t$ w modelu SIR

Na rysunku 5 widać kilka trajektorii procesu opisanego równaniami (5), wraz z zaznaczonym kolorem rozwiązaniem równania różniczkowego (4). Przyjęliśmy następujące wartości parametrów: $α = 2,$ $|β = 1,i0 = I(0) = 1,ℓ = 5000.$

Rys. 6 Liczba wszystkich dotychczas zainfekowanych. Przerywaną linią oznaczono rozmiar populacji. We wszystkich symulowanych przebiegach po wygaśnięciu epidemii około 1000 osób pozostaje w stanie $S$ (są to osoby, które nie zostały zarażone)

Początkowo trajektorie zachowują się tak, jak w modelu poprzednim: najpierw sporo losowych fluktuacji, później następuje faza wykładniczego wzrostu. Jeszcze później widać zupełnie inne zjawisko: hamowania wzrostu wskutek spadku liczby narażonych $|S(t),$ aż do zupełnego wygaśnięcia epidemii.

Dla różnych realizacji tego samego procesu SIR widać duże przesunięcia w fazie, wynikające z losowego charakteru początkowego fragmentu. Losowe fluktuacje (głównie początkowe) mają też pewien (niewielki) wpływ na maksymalną liczbę zarażonych (maksima poszczególnych krzywych).

W opracowaniach dotyczących rozwoju epidemii często przedstawia się również liczbę wszystkich osób, które zostały zainfekowane do danego momentu (tzn. liczbę zdarzeń polegających na zarażeniu nowej osoby). Liczba ta odpowiada procesowi $ℓ −S(t),$ widocznemu na rysunku 6. Warto zwrócić uwagę na niewielkie fluktuacje liczby wszystkich osób, które zostały zakażone do czasu wygaśnięcia epidemii.

Przedstawione w tym artykule modele rozwoju epidemii są skrajnie uproszczone i nie nadają się do ilościowego opisu rzeczywistego zjawiska. Niemniej nawet takie modele pozwalają trochę zrozumieć mechanizm epidemii na poziomie jakościowym.

Dlaczego czasy oczekiwania na skok mają rozkład wykładniczy?

Aby udowodnić, że zmienna $W$ ma rozkład wykładniczy z parametrem $|λ,$ wystarczy pokazać, że zmienna $e−λW$ ma rozkład jednostajny na odcinku $|[0,1],$ tzn. że dla dowolnego $t > 0$ : $P(W < t) = P(e −λW > e −λt) = 1 −e−λt.$

Przyjrzyjmy się czasowi pierwszego skoku, $|W1 = T1 = inf{t I(t) /= i0}.$ Nietrudno uwierzyć, że

P(W1 ⩽ t+h W1 > t) = P(I(t+h) /= i0 I(t) = i0)+o(h) = (α+β )i0h+o(h),

gdzie druga równość wynika z (2). Jeśli $F(t) = P(W1 ⩽ t),$ to dostajemy stąd równanie $|F t+h-−F-t--= (α +β )i0h + o(h), 1− F t$ zatem $1− F t |ddt log(1− F(t)) =-1− F-t --= −(α + β)i0t,$ skąd wnioskujemy, że $|1− F(t) = exp{ −(α +β )i0t}.$ Identyczne rozumowanie stosuje się do $|W n+1$ i prowadzi ono do analogicznego rezultatu, z $|i 0$ zastąpionym przez $|in.$