SIS ja SIRVI-malli

SIS-malli

Lähde löytyy täältä.

SIS-malli kuvaa tilannetta, jossa taudin sairastaminen ei tuotakaan pysyvää immuniteettia, vaan henkilö palaakin takaisin lähtötilanteeseen eli taudille alttiiksi.
\begin{equation}I(t+\triangle t) = I(t)+\beta \triangle t S(t) I(t) - \gamma \triangle t I(t), \end{equation}
kun kolmio-t lähestyy nollaa, saadaan differentiaaliyhtälö
\begin{equation}\frac{dI}{dt} =\frac{\beta SI}{N} - \frac{\gamma I}{N}. \end{equation}
sekä vastineeksi
\begin{equation} \frac{dS}{dt}=-\frac{\beta SI}{N}+\frac{\gamma I}{N}.\end{equation}

Oletus I+S=N saadaan käännettyä muotoon S=N-I, sekä sijoittamalla saatu muoto edelliseen differentiaalityhtälöön saadaan
\begin{equation}\frac{dI}{dt} \beta (N-I)I - \gamma I. \end{equation}
sekä edelleen muotoilemalla
\begin{equation}\frac{dI}{dt} = (\beta N - \gamma)I(1-\frac{\beta I}{N-\gamma}), \end{equation}
joka on siis logistinen kasvumalli.

Edellinen differentiaaliyhtälö voidaan edelleen ratkaista. Ratkaisun pituuden ja monimutkaisuuden vuoksi en käy ratkaisumenetelmää yksityiskohtaisesti tässä läpi, täydellinen ratkaisu löytyy tästä Khanin, Sabbirin ja Sadiqin artikkelista. Differentiaaliyhtälön ratkaisu lähtee ylempänä muotoillun yhtälön integroinnista. Kun muuttujan I differentiaaliyhtälö on saatu ratkaistua, saadaan muuttujan S ratkaisu yksinkertaisesti käyttämällä identiteettiä S = N-I, jolloin ratkaisuiksi saadaan
\begin{equation}S=N - \frac{\beta N - \gamma}{\beta - (\beta N - \gamma)\frac{\beta N - \gamma - i_0 \beta}{(\beta N - \gamma)i_0}e^{-(\beta N - \gamma )t}}. \end{equation}
sekä muuttujan I yhtälöksi
\begin{equation}I=\frac{\beta N - \gamma}{\beta - (\beta N - \gamma)\frac{\beta N - \gamma - i_0 \beta}{(\beta N - \gamma)i_0}e^{-(\beta N - \gamma )t}}. \end{equation}

Skipataan tässä yksinkertaisemmassa mallissa tasapainopisteiden haku sekä niiden stabiiliusominaisuuksien analyysi.  Aiheesta löytyy erinomainen esitelmä täältä.
 
Tartuntoja kuvaavien differentiaaliyhtälöiden kuvaaja on

SIS-mallissa perustartuntaluku on muotoa
\begin{equation}R_0 = \frac{\beta N}{\gamma}, \end{equation}
eli jos R_0 > 1, epidemia käynnistyy. Tällöin sairastuneella on vähintään yksi kontakti, ja tauti voi siten lähteä leviämään.

SIRVI-malli 

SIRVI-malli lähtee pitkälti samoista oletuksista kuin kaksi edellistäkin, sillä erotuksella, ettei rokotus tällä kertaa annakaan pysyvää immuniteettiä. SIRVIn asetelma muistuttaa siis vahvasti koronarokotteen vaikutuksia, ja on siksi hyvin ajankohtainen. Koronarokote ei siis antanut kestävää immuniteettia tautia vastaan, estänyt ketään sairastamasta tai levittämästä tautia, mutta vahvisti hetkellisesti immuniteettia tautia vastaan. SIRVIin löytyy erinomaisen perusteellinen lähde täältä, ja käytän sitä myös tämän blogipostauksen pohjana. Kaikki olennainen tieto on peräisin Turkyilmazoglun artikkelista, enkä ota tästä kunniaa itselleni.

SIRVI-malli muodostuu yhtälöistä:

\begin{equation} \frac{dS}{dt}=-\beta SI -vS\end{equation}

\begin{equation} \frac{dI}{dt}=\beta SI + \alpha VI -\gamma I\end{equation}

\begin{equation} \frac{dR}{dt}=\gamma I \end{equation}

\begin{equation} \frac{dV}{dt}=v S - \alpha VI \end{equation}

missä alfa, beta, gamma ja v ovat rokotettujen tartuntamäärä, epäiltyjen tartuntamäärä, parantumismäärä ja epäiltyjen rokotusmäärä.

Hetkellä t=0 annettujen funktioiden muuttujille pätee:

\begin{equation} S_0 = 1-\tau \end{equation}

\begin{equation} I_0 = \tau \end{equation}

\begin{equation} R_0 = 0 \end{equation}

\begin{equation} V_0 = 0 \end{equation}.

Tällöin on voimassa myös yhtälö

\begin{equation}S+I+R+V = 1 \end{equation}.
Mallin perustartuntaluku on
\begin{equation}R_0=\frac{S_0 \beta -V_0 \alpha}{\gamma}\end{equation}
sekä tasapainopisteet
\begin{equation}(S,I)=(0,0)\end{equation}
\begin{equation}(V,R)=(0 \leq V, R \leq 1)\end{equation}

SIRVI-mallin Jacobin matiiris on muotoa
\begin{equation}J=\begin{bmatrix}-v & -SB & 0 & 0 \\ 0 & V \alpha + S \beta - \gamma & 0 & 0 \\ v & -V \alpha & 0 & 0 \\ 0 & \gamma & 0 & 0  \end{bmatrix}\end{equation}

SIRVI-mallin stabiilisuusominaisuuksista voidaan ainakin sanoa, että nollasta eroavat ominaisarvot
\begin{equation}\lambda_1 = -v \end{equation}
sekä
\begin{equation}\lambda_2 = V \alpha + S \beta - \gamma, \end{equation}
mistä huomataan, että stabiilisuus pätee rokotetulla ja taudista vapaalla populaatiolla kun
\begin{equation}V \leq \frac{\gamma}{\alpha}\end{equation}

Kuvaaja on siis muotoa

Vaikka tässä käsiteltiinkin vain kolmea erityistä epidemiologista mallia, riittää vastaavia sovellutuksia maailmassa vaikka kuinka paljon. Epidemiologia on yksi sovelletun matematiikan merkittävimmistä haaroista.

Käytetyt koodit, perusten aiempaan SIR-mallissa käytettyyn Python-koodiin:

# SIS-malli
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# Total population, N.
N = 1000
# Initial number of infected and recovered individuals, I0 and R0.
I0 = 1
# Everyone else, S0, is susceptible to infection initially.
S0 = N - I0
# Contact rate, beta, and mean recovery rate, gamma, (in 1/days).
beta, gamma = 0.09, 2.0
# A grid of time points (in days)
t = np.linspace(0, 160, 160)

# The SIR model differential equations.
def deriv(y, t, N, beta, gamma):
    S, I = y
    dSdt = -beta * S * I / N + gamma * I/N
    dIdt = beta * S * I / N - gamma * I/N
    return dSdt, dIdt

# Initial conditions vector
y0 = S0, I0,
# Integrate the SIR equations over the time grid, t.
ret = odeint(deriv, y0, t, args=(N, beta, gamma))
S, I = ret.T

# Plot the data on three separate curves for S(t), I(t) and R(t)
fig = plt.figure(facecolor='w')
ax = fig.add_subplot(111, facecolor='#dddddd', axisbelow=True)
ax.plot(t, S/1000, 'b', alpha=0.5, lw=2, label='Susceptible')
ax.plot(t, I/1000, 'r', alpha=0.5, lw=2, label='Infected')
ax.set_xlabel('Time /days')
ax.set_ylabel('Number (1000s)')
ax.set_ylim(0,1.2)
ax.yaxis.set_tick_params(length=0)
ax.xaxis.set_tick_params(length=0)
ax.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = ax.legend()
legend.get_frame().set_alpha(0.5)
for spine in ('top', 'right', 'bottom', 'left'):
    ax.spines[spine].set_visible(False)
plt.show()

#################################################################

# SIRVI-malli
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# Total population, N.
N = 1000
# Initial number of infected and recovered individuals, I0 and R0.
I0, R0, V0 = 1, 0, 0
# Everyone else, S0, is susceptible to infection initially.
S0= N - I0 - R0 -V0
# Contact rate, beta, and mean recovery rate, gamma, (in 1/days).
beta, gamma, v ,alpha= 0.2, 2.10, 0.03, 0.03
# A grid of time points (in days)
t = np.linspace(0, 160, 160)

# The SIR model differential equations.
def deriv(y, t, N, beta, gamma, v, alpha):
    S, I, R, V = y
    dSdt = -beta * S * I/N - v * S
    dIdt = beta * S * I/N + alpha*V*I/N - gamma * I/N
    dRdt = gamma * I/N
    dVdt = v * S-alpha*V*I/N
    return dSdt, dIdt, dRdt, dVdt

# Initial conditions vector
y0 = S0, I0, R0, V0
# Integrate the SIR equations over the time grid, t.
ret = odeint(deriv, y0, t, args=(N, beta, gamma,v,alpha))
S, I, R, V = ret.T

# Plot the data on three separate curves for S(t), I(t) and R(t)
fig = plt.figure(facecolor='w')
ax = fig.add_subplot(111, facecolor='#dddddd', axisbelow=True)
ax.plot(t, S/1000, 'b', alpha=0.5, lw=2, label='Susceptible')
ax.plot(t, I/1000, 'r', alpha=0.5, lw=2, label='Infected')
ax.plot(t, R/1000, 'g', alpha=0.5, lw=2, label='Recovered with immunity')
ax.plot(t,V/1000, 'y', alpha=0.5, lw=2, label = 'vaccinated' )
ax.set_xlabel('Time /days')
ax.set_ylabel('Number (1000s)')
ax.set_ylim(0,1.2)
ax.yaxis.set_tick_params(length=0)
ax.xaxis.set_tick_params(length=0)
ax.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = ax.legend()
legend.get_frame().set_alpha(0.5)
for spine in ('top', 'right', 'bottom', 'left'):
    ax.spines[spine].set_visible(False)
plt.show()

SIR- ja SIRV-mallit sekä rokotuksen vaikutuksen mallinnus

SI-mallin johtaminen

Aiheen johdannoksi esitellään Charnovin monisteiden esittelemä tapa tartuntatautimallin määrittelemiseksi. Aloitetaan seuraavasti tunnetuksi oletetusta kaavasta
\begin{equation}I(t+\triangle t) = I(t)+ \beta \triangle t S(t) I(t)\end{equation}
missä t kuvaa aikaa, I(t) tartuntojen määrää ja S(t) alttiiden määrää hetkellä t sekä beta kontaktien määrää aikayksikössä, eli leviämisintensiteettiä. Viemällä I(t) funktion vasemmalle puolelle, jakamalla molemmat puolet kolmio-t:llä, kun kolmio-t lähestyy nollaa, saadaan tartuntojen määrän muutosta kuvaava differentiaaliyhtälö
\begin{equation} \frac{dI}{dt}=\beta SI, \end{equation}
joka on tuttu myöhemmistä aiheen variaatioista. Myöhempiä malleja ei johdeta yhtä perusteellisesti, vaan käytännössä aloitetaan valmiista ja annetuiksi oletetuista differentiaaliyhtälöistä.

SIR-malli

SIR-malli kuvaa epidemian dynamiikkaa olettamalla olevaksi kolme tilaa: alttiiden S, todetun tautitapauksen I sekä parantuneen R. Tässä mallissa ei oteta huomioon syntymiä eikä kuolemia. Kerran saatu immuniteetti ei häviä, eikä samaan tautiin voi sairastua toistamiseen. Näistä oletuksista on erilaisia vapaampia muotoiluja, kuten SIRI-malli, jossa henkilölle ei muodostu kestävää immuniteettiä, SIRD-malli, jossa henkilön kuolema otetaan huomioon, sekä joukko erilaisia rokotusmalleja, joista tässä käsiteltävässä immuniteetti on pysyvä. Vaihtoehtoisia muotoiluja oletuksineen käsitellään erikseen myöhemmin.

Yksinkertainen SIR-malli muodostuu differentiaaliyhtälöistä
\begin{equation} \frac{dS}{dt}=\frac{-\beta SI}{N} \end{equation}

\begin{equation} \frac{dI}{dt}=\frac{\beta SI}{N} -\gamma I\end{equation}

\begin{equation} \frac{dR}{dt}=\gamma I. \end{equation}
SIR-mallin parametreista beta kuvaa tartuntaintensiteettiä ja gamma poistumis- tai parantumisintensiteettiä riippuen miten mallia tulkitaan. Kahden parametrin avulla voidaan edelleen määritellä perustarttuvuusluku
\begin{equation} R_0 = \frac{\beta}{\gamma}, \end{equation}
joka karkeasti tarkoittaa odotettujen tartuntojen määrän. Seuraavaksi huomataan, että sijoittamalla perustarttuvuusluku tartuntojen muutosta kuvaavaan differentiaaliyhtälöön, saadaan
\begin{equation}\frac{dI}{dt}=(R_0 \frac{S}{N}-1)\gamma I, \end{equation}
mistä edelleen huomataan, että kun pätee
\begin{equation}R_0 S > N, \end{equation}
on infektioiden muutos samalla kasvava, eli tartuntojen määrä lisääntyy ja epidemia on kiihtymisvaiheessa. Käänteisessä tapauksessa tartuntojen määrä vähenee ja epidemia on asettumassa. Sama voitaisiin todeta toisin, eli olettamalla likimain jokaisen populaatioon kuuluvan yksilön olevan altis eli S = N, seuraisi siitä, että epidemian käynnistymiseen riittää ehto
\begin{equation}R_0 > 1, \end{equation}
sekä talttumiseen käänteinen tulos.

SIR-mallin kuvaajassa näkyvät eri ryhmien muutosten suunnat sekä kasvunopeudet. Seuraavan SIR-mallin kuvaajaan löytyy koodi täältä.

Kuvaajasta huomataan, että alttiiden määrä vähenee taudin leviämisen sekä siitä parantuneiden kasvamisen mukana. Kaikki muutosnopeudet tasaantuvat ajan mittaan. Samaista koodia sovelletaan mahdollisuuksien mukaan myös saman mallin eri sovelluksissa, kuten seuraavassa SIRV-mallissa.

SIRV-malli

Määritellään johdannoksi rokotusmalliin seuraava tulos. Olkoon S_0 alttiiden prosentuaalinen osuus hetkellä 0, olkoon lisäksi p rokotettujen osuus väestöstä, ja olkoon lisäksi voimassa S_0=(1-p). Jos oletetaan rokotuksen tuoma immuniteetti pysyväksi, voidaan kirjoittaa edellisen kohdan perustarttuvuuslukua hyödyntäen seuraava ehto epidemian käynnistymiselle:
\begin{equation} R_0 (1-p) > 1. \end{equation}
Annetun ehdon avulla voidaan määritellä rokotettavien minimäärä p* epidemian estämiseksi, jolloin voidaan kirjoittaa
\begin{equation} p^* = 1-\frac{1}{R_0}.\end{equation}
SIRV-malli on perusteiltaan kuin aiemmin kuvailtu yksinkertainen SIR-malli, mutta siinä on otettu huomioon rokotuksen vaikutus. Rokotus antaa tässä mallissa pysyvän immuniteetin, eikä tautiin kuolleita oteta tässäkään huomioon.
\begin{equation} \frac{dS}{dt}=\frac{-\beta SI}{N} -vS\end{equation}

\begin{equation} \frac{dI}{dt}=\frac{\beta SI}{N} -\gamma I\end{equation}

\begin{equation} \frac{dR}{dt}=\gamma I \end{equation}

\begin{equation} \frac{dV}{dt}=v S \end{equation}
Tässä mallissa parametrit tarkoittavat samaa kuin edellisessä mallissa, mutta uutena SIRV-malliin on lisätty rokotusintensiteettiä kuvaava parametri v. Kun muut parametrit ovat samat kuin ensimmäisessä ja rokotusintensiteetti v = 0.005, saadaan kuvaaja

Kuvasta siis näemme, että rokotuksilla saadaan taudin leviäminen hallintaan ja immuniteetti puolelle samasta määrästä ihmisiä, jotka saivat saman sairastamalla taudin. Infektoiden määrän kasvuvauhtia saatiin näin loivennettua ja huippu madallettua.

SIRV-mallissa käytetty Python-koodi, muokattu aiemmin linkatusta esimerkistä:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

N = 1000
I0, R0, V0 = 1, 0, 0
S0 = N - I0 - R0 -V0
beta, gamma, v = 0.2, 1./10, 0.005
t = np.linspace(0, 160, 160)

def deriv(y, t, N, beta, gamma, v):
    S, I, R, V = y
    dSdt = -beta * S * I / N - v * S
    dIdt = beta * S * I / N - gamma * I
    dRdt = gamma * I
    dVdt = v * S
    return dSdt, dIdt, dRdt, dVdt

y0 = S0, I0, R0, V0

ret = odeint(deriv, y0, t, args=(N, beta, gamma,v))
S, I, R, V = ret.T

fig = plt.figure(facecolor='w')
ax = fig.add_subplot(111, facecolor='#dddddd', axisbelow=True)
ax.plot(t, S/1000, 'b', alpha=0.5, lw=2, label='Susceptible')
ax.plot(t, I/1000, 'r', alpha=0.5, lw=2, label='Infected')
ax.plot(t, R/1000, 'g', alpha=0.5, lw=2, label='Recovered with immunity')
ax.plot(t,V/1000, 'y', alpha=0.5, lw=2, label = 'vaccinated' )
ax.set_xlabel('Time /days')
ax.set_ylabel('Number (1000s)')
ax.set_ylim(0,1.2)
ax.yaxis.set_tick_params(length=0)
ax.xaxis.set_tick_params(length=0)
ax.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = ax.legend()
legend.get_frame().set_alpha(0.5)
for spine in ('top', 'right', 'bottom', 'left'):
    ax.spines[spine].set_visible(False)
plt.show()