Das SIR-Modell mit Covid

Was ist ein SIR Modell?

Bei mathematischer Modellierung geht es darum, eine schematische, häufig grob vereinfachte, Abstraktion eines Sachverhalts vorzunehmen. Im Beispiel von Epidemien, um die es heute gehen soll, muss eine Gruppe von Menschen in unterschiedliche Klassen aufgeteilt werden. Eine einfache Aufteilung geht beispielsweise so, dass man die Bevölkerung in Infizierte , Gesunde , und Genesene einteilt. Jedes Individuum gehört genau zu einer Gruppe. Im Englischen spricht man von Susceptibles, Infected und Recovered; kurz eben SIR. Im einfachsten Fall geht man davon aus, dass die Gesamtheit konstant bleibt und dass einmal Infected, die genesen, immun werden.

Vermutlich neu für die meisten von euch ist nun, dass es hier nicht um eine Momentaufnahme geht, sondern dass wir Funktionen als Lösungen eines dynamischen Systems suchen. Wir möchten zum Beispiel voraussagen können, wie viele zu einem Zeitpunkt $t$ infiziert sein werden, oder vielleicht wie lange es dauert, bis alle infizierten immun sind. Mathematisch heisst das, dass wir Funktionen $S$ , $I$ und $R$ in Abhängigkeit der Zeit $t$ suchen, die uns solche Fragen beantworten können.

$S(t)$ : Anzahl der Empfänglichen zum Zeitpunkt $t$ .
$I(t)$ : Anzahl der Infektiösen zum Zeitpunkt $t$ .
$R(t)$ : Anzahl der Genesenen zum Zeitpunkt $t$ .

Wie finden wir diese Funktionen?

Man sieht rasch, dass es schwierig ist, die Funktionen zu erraten. Jedoch: Wir können etwas über die zeitliche Änderung dieser Funktionen sagen.

Die Anzahl Personen, welche an einem Tag gesunden ist proportional zur Anzahl der infizierten:

\dot{R}(t) = \beta\cdot I(t)

Wenn es im Schnitt $k$ Tage dauert, bis man wieder fit ist, so hat $\beta$ die Grössenordnung $\frac{1}{k}$ .

Die Anzahl der Personen, die an einem Tag infiziert werden ist sowohl proportional zur Anzahl Infektiösen als auch zur Anzahl der Empfänglichen:

\dot{S}(t) = -\alpha\cdot S(t)\cdot I(t)

Den Proportionalitätsfaktor $\alpha$ kann man nicht so leicht abschätzen wie $\beta$ . $\alpha$ hängt unter anderem sicherlich von der Art der Krankheit und von der Anzahl der Kontakte zwischen den Individuen ab.

Nun ist klar, wie sich nach unserem Modell die Gruppe der Infizierten ändert:

\dot{I}(t) = \alpha\cdot S(t)\cdot I(t) - \beta\cdot I(t)

Vom Diskreten zum Kontinuierlichen

Methode

Wir haben es mit einem diskreten Modell zu tun, dass wir mit kontinuierlichen Methoden angehen werden. Damit wird es "geschmeidiger", kann aber auch nicht plausible Werte liefern wie Fliesskommazahlen. Oft ist es in der Praxis so, dass man die Lösungen numerisch bestimmen muss, dass man also sowieso kleine Zeitschritte einbaut und nicht analytisch die Lösung findet. Uns ist aber ohnehin klar, dass es beispielsweise $14.2857$ Gesunde nicht geben kann. Aber die Methode des Grenzwertübergangs liefert uns präzisere Daten.

Schliesslich wird jetzt noch einmal dieses Differentialgleichungssystem schlank im Überblick projeziert:

\begin{align} \dot{S} &= -\alpha S I\\ \dot{I} &= \alpha S I - \beta I\\ \dot{R} &= \beta I \end{align}

Plausibilität unseres SIR Modells

Einheitencheck

Stehe $t$ wie üblich für die Zeit. Betrachten wir das Differentialgleichungssystem

\begin{align} \dot{S} &= -\alpha S I\\ \dot{I} &= \alpha S I - \beta I\\ \dot{R} &= \beta I \end{align}

so sehen wir, dass $\beta$ die Einheit $[\frac{1}{t}]$ haben muss, da ja die linken Seiten die Einheit $[\frac{\mathrm{Individuen}}{t}]$ haben. Etwas Ähnliches haben wir uns vorher schon überlegt. Analog findet man, dass $\alpha$ die Einheit $[\frac{1}{\mathrm{Individuen}\cdot t}]$ hat.

Addieren der Gleichungen des Systems

Wir erhalten $(S+I+R)\dot{}=0$ . Das heisst, die Population ist konstant, wie wir angenommen haben.

Abschätzung der Startwerte und qualitative Verläufe

Ganz zu Beginn des Ausbruchs der Epidemie lässt sich unser System grob so formulieren:

\begin{align} \dot{S} &= -\alpha S I\\ \dot{I} &= \alpha S I \end{align}

Es kann ja noch gar niemand genesen sein, wenn die Krankheit beispielsweise "neu" ist. Das heisst wir haben $R = 0$ und damit $\dot{R} = 0$ . Ist zudem - und dies ist in der Praxis eigentlich immer gegeben - $S$ zu Beginn viel grösser als $I$ , so kann $S$ für kleine $t$ als konstant betrachtet werden. Ist nun $S$ annähernd konstant, dann ist $\dot{S} \approx 0$ und daher haben wir:

\begin{align} \dot{I} &= \gamma I \end{align}

wobei, $\alpha S =: \gamma$ .

Was für eine Funktion erfüllt diese Gleichung? Genau, das ist

I(t)=I_0\cdot e^{\gamma t}.

Deshalb haben wir zu Beginn exponentielles Wachstum.

Ignoriert man $R$ , wenn zum Beispiel der Genesungsprozess verhältnismässig lange dauert oder man nicht immun wird, dann haben wir

\begin{align} \dot{S} &= -\alpha S I\\ \dot{I} &= \alpha S I \end{align}

Da wir annehmen, dass die Gesamtpopulation konstant bleibt, gilt $S+I=1$ für alle Zeiten $t$ , also $S=(1-I)$ . Somit ergibt sich, da wir eigentlich in diesem Fall nur eine Gleichung brauchen,

\begin{align} \dot{I} &= \alpha (1-I) I \end{align}

was logistisches Wachstum bedeutet.

Übrigens ist diese Gleichung analytisch lösbar, indem man die Variablen trennt und anschliessend mit Hilfe einer Partialbruchzerlegung einfache Integrale bereit stellt. Wir werden das im AM erledigen. Die Neugierigen können ja bereits Stift und Papier zur Hand nehmen und selber rechnen.

Ein Geogebra-File zum Ausprobieren:

Open in GeoGebra

Das SIR Modell - jetzt quantitativ

Ich habe die Gleichungen im nächsten Kapitel in Geogebra bereit gestellt. Dort wird durch das Programm dynamisch der Verlauf berechnet und geplotet. Bei den Schiebereglern kann man die Parameter $\alpha$ und $\beta$ ändern. Es sei noch vermerkt, dass wir eigentlich nur die ersten beiden Gleichungen brauchen, die dritte Grösse (z.B. $R$ ) ergibt sich jeweils aus der konstanten Populationsgrösse. Nun folgen diverse statische Plots mit Python, die interessante Details zeigen.

Es ist also möglich, durch Ändern des Parameters $\alpha$ , den zeitlichen Verlauf der Krankheit innerhalb einer Population zu beeinflussen. Im ersten Beispiel ist der Peak der Infected markant und zwischenzeitlich sind mehr als $80%$ infiziert. Das muss man vermeiden! Im zweiten Beispiel wurde $\alpha$ reduziert auf einen Zentel, was zum Beispiel mit Massnahmen wie Social Distancing , regelmässig Händewaschen etc. erreicht werden kann. Bemerkenswert ist auch, dass im letzten Beispiel nicht alle Susceptibles krank werden und die Krankheit ausstirbt. Beachte auch, dass man für den zeitlichen Verlauf aber viel mehr Geduld braucht.

Daten aus der Schweiz (Wikipedia)

Auf dieser Site sind einige Zahlen/Daten aus der Schweiz zum COVID-19 unter anderem auch graphisch präsentiert. Insbesondere kann hier zum Beispiel der exponentielle und logistische Verlauf gefunden werden.

Parameterschätzung --- für den Kanton Bern anhand der März-Zahlen

Man kann nun zum Beispiel $I$ anpassen, indem man die ersten bekannten Messwerte von $I$ nimmt und ein Curve Fitting durchführt. In den ersten paar Tagen der Epidemie gilt wie oben gesehen näherungsweise $I(t)=I_0\cdot e^{\gamma\cdot t}$ . So erhält man eine Schätzung für $\gamma$ und kann daraus eine Schätzung für $\alpha$ ermitteln, da zu Beginn $\gamma=\alpha\cdot S_0$ .

Hier die Daten für den Kanton Bern im März (nachgewiesene Infektionen):

Datum im März	Zahl der nachgewiesenen Infektionen
14.	78
16.	123
18.	193
19.	282
20.	377
21.	418
23.	470
24.	532
25.	624

Leider hat der Kanton nicht für jeden Tag Daten, so dass ab und zu ein Tag ausgelassen wird. Ich habe ohne Rücksicht auf Verluste einfach diejenigen Werte kumuliert, die ich zur Verfügung hatte. Natürlich "fehlen" dann jeweils nach einem datenlosen Tag diese in der nächsten Aufsummierung. Für das erhobene fitting bedeutet das, dass der erhaltene Wert von $\gamma$ konservativ ausfällt; also eine milde Schätzung ist.

WolframAlpha beispielsweise kann so ein Fitting für uns tun:

FindFit[{{0,78},{2,201},{4,394},{5,676},{6,1013},{7,1431},{9,1901},{10,2433},{11,3067}}, 78*Exp[c t], {c}, t],

wobei ich wie erwähnt die nachgewiesenen Infektionen kumuliert habe.

Es findet für den Ansatz $f(t)=78\cdot e^{ct}$ den Wert $c=0.341643$ , mit unserer Notation

\gamma = 0.341643.

Die Daten habe ich aus der Zeit der Schulschliessung. Man kann heute, 31.Mai 2020, eine Reproduktionszahl von $R_0=1.22$ für Mitte März nachschlagen, als die Massnahmen bereits in Kraft waren. Vorher war diese Zahl bei etwa $3$ anzusetzen und aufgrund des Fittings noch höher. Beachte auch, dass das Fitting vermutlich anders ausgefallen wäre, wenn ich die beiden letzen Datenpunkte ausgelassen hätte. Was würde dies für unser Modell für $\alpha$ heissen?

Es ist $\gamma=\alpha\cdot S_0$ , also kriegen wir etwas pingelig $\alpha=\frac{\gamma}{S(0)-I(0)}=$ . Die Reproduktionszahl betrüge demnach $R_0=\frac{\alpha}{\beta}\cdot S_0=6.833\dots$ . Wollen wir wissen, welchem Wert $\alpha$ in etwa eine Rerpoduktionszahl von $R_0=3$ entspricht, so stellen wir nach $\alpha$ um und erhalten 0.15 gerundet. Der Wert $\alpha=0.15$ entspräche also ungefähr einem Wert der Reproduktionszahl von $R_0=3$ für unser Modell.

Man will natürlich die Reproduktionszahl unter $1$ bringen. Dies hiesse für $\alpha$ etwa 0.05. Der Verlauf wäre nun für ein $R_0$ knapp grösser $1$ .

Natürlich kann man jetzt mit schickeren Modellen arbeiten. Wenn ihr Lust habt, so hoffe ich gezeigt zu haben, wie man von diesem Punkt aus hier die Modelle auch verfeinern kann. Stichworte zum Vertiefen sind SEIR , SIRV , recursiv SIR , zellurläre Automaten ...

Die "echten" Daten kann man zum Beispiel auf der offiziellen Site COVID 19 Informationen Kanton Bern einsehen und mit den Modellen vergleichen. Aber auch bei offiziellen Seiten immer gut selber gucken, welche Daten wie verwendet wurden.

Abschliessend noch einige mathematische Bemerkungen zu unserem SIR Modell.

Mathematische Analyse des SIR Modells

Es folgen noch einige Rechnungen, welche auch die Begriffe Basisreproduktionszahl und Herdenimmunität mathematisch greifbar werden lassen.

Lambert'sche W-Funktion