(c) I. Dahn, 2020

Das SI-Modell und das SIR Modell in der COVID-19-Pandemie in Deutschland im April 2020

Ingo Dahn (dahn@dahn-research.eu)

Kontext

Wie gut beschreiben mathematische Modelle reale Prozesse?

Die Ausbreitung der COVID-19-Pandemie in Deutschland konnte bis Ende März 2020 gut durch ein exponentielles Modell beschrieben werden. Ende März verschlechterte sich die Passung des exponentiellen Modells zu den gemeldeten Infektionsdaten zusehends. Zu diesem Zeitpunkt traten außerdem in Deutschland umfangreiche Kontaktbeschränkungen in Kraft mit dem Ziel, die Ausbreitung der Pandemie zu verlangsamen.

In diesem Jupyter-Notebook untersuchen wir die Passung der einfachsten Pandemie-Modelle - des SI-Modells und des SIR-Modells - zu den Daten der Pandemie in Deutschland im April 2020, wie sie vom Robert-Koch-Institut zur Verfügung gestellt wurden. Zunächst wird das System von Differentialgleichungen zusammengestellt, das diese Modelle beschreibt. Die zu untersuchende Frage ist, wie die Parameter dieses Differentialgleichungssystems gewählt werden müssen, damit die realen Daten möglichst gut approximiert werden und welche Qualität der Approximation so erreicht werden kann.

Für das SI-Modell bestimmen wir die Parameter mittels logistischer Regression gegen die Lösungsformel des SI-Differentialgleichungssystems. Für das SIR-System steht eine solche Lösung nicht zur Verfügung - hier verwenden wir eine numerische Näherungslösung, die mit Hilfe des Runge-Kutta-Verfahrens berechnet wird. Wir diskutieren dabei auch Möglichkeiten zur Verbesserung der Genauigkeit bzw. zur Reduzierung der erforderlichen Rechenzeit.

Ein abschließendes Fazit fasst die Ergebnisse zusammen.

In der ausführbaren Ansicht dieses Notebooks in CoCalc bzw. im CoCalc-Player können die Eingabezellen editiert und neu berechnet werden um die Verfahren zu prüfen oder sie auf andere Daten anzuwenden.

Dabei müssen die Zellen immer in der gegebenen Reihenfolge ausgeführt werden!

Falls Sie dabei bessere Approximationen oder Fehler finden, so lassen Sie es mich wissen.

Diese Seite wird unter der Creative Commons Lizenz CC BY-NC-SA 4.0 veröffentlicht.

Koblenz im Mai 2020

Dr. Ingo Dahn

Modelle

SI- und SIR-Modell teilen eine Bevölkerung in zwei bzw.drei Gruppen ein:

1. $S$: die noch nicht infizierten, aber anfälligen,
2. $I$: die infizierten und daher ansteckenden und,
3. $R$: die genesenen (oder gestorbenen), und daher nicht mehr ansteckenden Personen.

Differentialgleichungen

Bezeichnen wir die Anzahl der für das Modell wesentlichen Personen mit $N$, so nehmen die Modelle an, dass $$N = S + I + R.$$

Die Beziehung zwischen den Größen $S,I$ und $R$ wird im SIR-Modell durch die folgenden Differentialgleichungen beschrieben. $$ \begin{align} \frac{dS}{dt} &= -c \frac{S}{N} I &= -\frac{c}{N} S I\ \frac{dI}{dt} &= c \frac{S}{N} I - w I &= \frac{c}{N} S I - w I\ \frac{dR}{dt} &= w I \end{align} $$

$\frac{dS}{dt} = -c \frac{S}{N} I = -\frac{c}{N} S I$

$\frac{dI}{dt} = c \frac{S}{N} I - w I = \frac{c}{N} S I - w I$

$\frac{dR}{dt} = w I$

Dabei ist $c$ die Infektionsrate und $w$ die Genesungsrate. Das SI-Modell ist der Spezialfall des SIR-Modells der Genesungen und Todesfälle nicht berücksichtigt, d. h. $w=0$ bzw. $R$ ist konstant und damit $S+I=N$. Für das SI-Modell haben die Differentialgleichungen die Lösung $$I(t)=a\frac{N}{a+(N-a)\cdot e^{-ct}}.$$

Dabei ist $a=I(0)$.

Für das vollständige Gleichungsssystem des SIR-Modells ist keine Lösung in geschlossener Form bekannt.

Daten

Im Folgenden wird statt mit Kalender-Daten, mit Tagen seit dem 24.2.2020, dem Beginn der täglichen Datenemeldungen des Robert-Koch-Instituts (Tag 0) gerechnet. Die folgende Tabelle ermöglicht eine Umrechnung.

Es stehen taggenaue Meldungen des RKI zu

• gemeldeten Infektionen
• geschätzte Genesungen (auf volle 100 gerundet)
• gemeldete Todesfälle

zur Verfügung. Der April 2020 umfasst die Tage 37-66.

Die Zahl der Genesenen entwickelt sich - seit einem Sprung am 8.4.2020 - kontinuierlich. Diesen Sprung, der durch eine Änderung der Zählweise entstand, korrigieren wir, indem wir die Zahlen der Genesenen vor dem 8.4. (Tag 44 der Datenerfassung) um 30% nach oben korrigieren

Das folgende Diagramm zeigt die vom RKI gemeldeten Daten (einschließlich unserer Korrektur für die Genesenen), wobei die Werte normiert wurden. In absoluten Zahlen waren die maximalen Werte in dieser Zeit

• für die gemeldeten Infektionen (kumuliert): 159 119
• für die gemeldeten Todesfälle (kumuliert): 6 288
• für die geschätzte Zahl der Genesenen (kumuliert): 123 500

Hochrechnung

Im April 2020 wurden überwiegend Personen mit Symptomen getestet. Die Zahl der Test ist nur für die einzelnen Wochen bekannt und schwankt sehr stark:

• 30.3-5.4.: 408 348
• 6.4.-13.4.: 379 233
• 14.4.-21.4.: 330 027
• 21.4.-28.4.: 467 137

Informationen über Ergebnisse von randomisierten Tests in dieser Zeit liegen nicht vor, was eine Hochrechnung erheblich erschwert. Schätzungen über den Dunkelfaktor $d$ - das Verhältnis der Zahl der Infizierten zur Zahle der gemeldeten Infektionen - schwanken stark (A. Kekule (17.3.20): 4, L. H. Wieler: 2, Universität Göttingen (2.4.2020): 7).

Bewertung der Modelle

Wir berechnen Funktionen, die die beobachteten Daten möglichst gut erklären sollen. Wie gut eine Datenreihe dataP durch eine Funktion $f$ approximiert wird messen wir mit dem relativen Residuum rR2(dataP,f) wie folgt. Ist $\vec{vd}$ der Vektor der beobachteten Daten und $\vec{vf}$ der Vector der entsprechenden Funktionswerte, so definieren wir $rR2(\vec{vd},f)=\frac{|\vec{vd}-\vec{vf}|}{|\vec{vd}|}$. Je geringer dieser Wert desto genauer die Approximation.

SI-Modell

Zahl der Infektionen

In der logistischen Phase, ab Anfang April, kann das Wachstum der Zahl der gemeldeten Infektionen immer schlechter durch Exponentialfunktionen beschrieben. das Virus stößt - aus welchen Gründen auch immer - auf Faktoren, die seine Ausbreitung behindern (z.B. Senkung der Reproduktionsrate durch Einschränkung sozialer Kontakte, hoher Grad der Durchseuchung der Bevölkerung).

Ein solches Wachstum kann durch eine logistische Funktion $a\cdot \frac{N}{a+(N-a)\cdot e^{-c x}}$ beschrieben werden. Dabei ist $a$ der Anfangswert bei $x=0$, $N$ ist eine angenommene Obergrenze, der sich die Zahl der Infizierten asymptotisch annähert (Kapazitätsgrenze) und $c$ ist ein Parameter der die Geschwindigkeit dieser Annäherung beschreibt. $c$ und $N$ sind Parameter in der logiszischen Differentialgleichung $$\dot{I}=\frac{c}{N}\cdot S \cdot I = \frac{c}{N} \cdot (N-I) \cdot I.$$

Die logistische Funktion ergibt sich als Lösung der logistischen Differentialgleichung die ihrerseits aus dem SI-Modell der Entwicklung von Epidemien abgeleitet ist.

Die Zahl der Infizierten ergibt sich aus dem $d$-fachen der vom RKI gemeldeten Infizierten einschließlich Genesenen und Toten, da das SI-Modell diese nicht berücksichtigt: $$I_{SI} = d\cdot I_{RKI}$$

Wir testen auch, wie gut sich das SI-Modell zur Beschreibung der Zahl der Toten anwenden lässt, die wir unverändert aus den RKI-Daten übernehmen:

$$D_{SI} = D_{RKI}$$

Dies ist von Interesse, da wir annehmen, dass die Entwicklung der Zahl der Toten die Tendenz der Ausbreitung der Pandemie - mit einer Verzögerung von etwa 10 Tagen zwische Meldung der Krankheit und Tod - besser widerspiegelt als die Zahl der gemeldeten Infektionen, die zusätzlich vom Umfang des Testens abhängig ist. Dabei ist zu erwarten, dass für die Differentialgleichungen zur Beschreibung der Todesfälle andere Parameter $a,N,c$ gewählt werden müssen, als zur Beschreibung der gemeldeten Infektionen.

Die Entwicklung der Zahl der gemeldeten Infektionen im März ließ sich gut durch eine Exponentialfunktion mit einem Anfangswert $a=400$ und einem Exponenten $0.14x$ beschreiben. Im logistischen Modell ist der Exponent $cx$, wobei $N$, als $d$-faches des Maximalwerts der gemeldeten Infektionen, in der Größenordnung von mehreren $1 000 000 = 10^6$ zu erwarten ist.

Die Approximation der Daten $I_{SI}$ mit Hilfe der logistischen Regression durch die Funktion find_fit erweist sich als schwierig. Ob und wie gut dies funktioniert hängt stark von den gewählten Anfangswerten ab. Mit den Anfangswerten $a=400, N=1.6\cdot 10^{6}, c=0.14$ erhalten wir kein brauchbares Ergebnis: Versuchen Sie es elbst:

Da find_fit uns kein brauchbares Ergebnis liefert, wählen wir eine andere Regressionsmethode. Unsere globale Regression geht für alle Parameter der logistischen Differentialgleichung von einem Bereich aus, in dem wir geeignete Parameter vermuten. Dieser Bereich wird glichmäßig in Teilbereiche aufgeteilt, deren Anzahl durch eine natürliche Zahl - die Granularität gran gegeben ist. Für jedes Tupel von Parametern aus dem so definierten Gitter berechnen wir das relative Residuum bezüglich unserer Datenreihe. Nachdem dies für alle Tupel berechnet wurde, geben ein Tupel zurück, für das das relative Residuum minimal ist sowie den Wert des Residuums.

Da wir 3 Parameter, $a,N,c$, suchen haben wir die Berechnungen für $g^3$ Tupel durchzuführen wenn wir mit gran=$g$ rechnen. Bei $g=10$, das wir im Folgenden verwenden werden, sind dies 1 000 Tupel, bei $g=100$ sind es 1 000 000 Tupel - die Rechenzeit ver-1000-facht sich!

Wir suchen $a \in [200,10^4], N\in [10^6,2\cdot 10^6], c\in [0.1,0.2]$ mit einer Granularität von 10. Experimentieren Sie mit anderen Werten!

Wir erhalten einen Fehler von etwa 2%. Sie können diesen Fehler - bei längerer Rechenzeit - verringern indem Sie die Granularität erhöhen.

Als Gegenmittel gegen das Anwachsen der Rechenzeit können Sie aber auch die Größe der zu durchsuchenden Intervalle verkleinern. Sie können etwa zunächst mit der Granularität 10 rechnen und dann noch einmal mit der Granularität 10 für die Intervalle mit $\frac{1}{10}$ der bisherigen Länge um die im ersten Schritt erhaltenen Parameterwerte. Dann rechnen Sie, aufgrund der kleineren Intervalle, mit dem 100fachen der bisherigen Genauigkeit, aber die Rechenzeit dafür hat sich nur verdoppelt, nicht ver-1000-facht! (Das Teilgebiet der Mathematik, dass den Ressourcenbedarf von Algorithmen untersucht ist die Komplexitätstheorie)

Mit der Funktion g_regression, die in der folgenden Zelle definiert wird, können Sie dieses Verfahren depthmal durchführen. Die Funktion gibt die gefundenen Parameter und das erreichte relative Residuum zurück. Versuchen Sie es!

Wir berechnen analog ein Modell für die Zahl der Toten. Bei Vorgaben mit Toleranzen von $\pm 50$ für $a$, $\pm 5000$ für $N$ und $\pm 0.01$ für $c$ und 10 Schritten für jedes Argument erhalten wir eine Approximation mit einem Residuum von 3.5 %. Bei 2-maliger Anwendung der Regression sinkt der Relative Fehler auf 3.1 % mit den Werten $a=19, N=8 200, c=0.1108$.

SIR-Modell

Für das SIR-Modell multiplizieren wir die Zahl der vom RKI gemeldeten genesenen mit dem Dunkelfaktor um die Zahl der Genesenen $H_{SIR}$ für das SIR-Modell abzuschätzen.

$$H_{SIR} = d\cdot H_{RKI}$$

Die Zahl der Verstorbenen übernehmen wir unverändert.

$$D_{SIR} = D_{RKI}$$

Die Zahl der Infizierten ergibt sich aus dem $d$-fachen der vom RKI gemeldeten Infizierten ohne Genesene und Tote:

$$I_{SIR} = d\cdot (I_{RKI}-H_{RKI}-D_{RKI})$$

Die Zahl $R$ der nicht infektiösen für das SIR-Modell ist die Summe aus der Zahl der Toten und der Zahl der Genesenen des Modells.

$$R_{SIR} = D_{SIR}+H_{SIR}$$

Schließlich ergibt sich die Zahl $I$ der Infizierbaren des Modells durch Abzug der Infizierten und der nicht infektiösen von der Gesamtzahl:

$$S_{SIR} = N - I_{SIR} - R_{SIR}$$

Beachten Sie, dass die Datenreihen $I_{SIR}$ und $R_{SIR}$ sich aus den Angaben des RKI berechnen lassen, während $S_{SIR}$ von dem zu schätztenden Wert des Parameters $N$ des Modells abhängt. Statt der Parameter $a,N,c$, für die wir im SI-Modell mögliche Werte $a=19, N=1 600 000, c=0.1108$, erhalten haben, benötigen wir im SIR-Modell die Parameter $N,c,w$ für die Berücksichtigung von Genesungen und Sterbefällen. Statt des geschätzten Anfangswertes $a$ des SI-Modells verwenden wir die Werte der Datenreihen zu Beginn der Periode.

Wir suchen Werte für diese Parameter, so dass die Lösungen des SIR-Differentialgleichungssystems die Werte der Datenreihen $S_{SIR},I_{SIR}$ und $R_{SIR}$ von den Funktionen $S,I$ und $R$ möglichst gut approximiert werden. Für $N$ können wir dann die Kapazitätsgrenze des berechneten Modells nehmen.

Wir definieren die Differentialgleichungen und ihre numerischen Lösungen mit dem Runge-Kutta-Verfahren. Dabei werden die Anfangswerte zum Beginn der Periode vorgegeben und das Runge-Kutta-Verfahren berechnet schrittweise den Verlauf der Lösungen des Differentialgleichungssystems.

Das Runge-Kutta-Verfahren zur numerischen Lösung von Differentialgleichungen liefert statt eines Funktionsausdrucks drei Datenreihen. Dies erfordert Anpassungen unserer Werkzeuge zur Fehlerabschätzung und Regression.

Die Funktion regression_l_g erhält als Parameter die realen Datenreihen für die Infizierten $I$ sowie für die Genesenen und Verstorbenen $R$, die zu durchsuchenden Bereiche für die Parameter $N,c,w$ und die Granularität $gran$ der Suche.

Die dritte Datenreihe $S$ wird aus den beiden gegebenen und dem aktuellen wert des Parameters $N$ mit Hilfe der Funktion makeS berechnet und bei Änderungen dieses Parameters angepasst. Wir suchen nach Parameterwerten für die das Maximum des relativen Fehlers der numerischen Lösung des Differentialgleichungssystems verglichen mit den aktuellen Datenreihen möglichst klein wird.

Wir berechnen die Wertereihen des Modells für S,I und R mit unserer Funktion regression_l_g. Werden dabei für jede der Parameter $N,c,w$ 10 Schritte vorgesehen, so muss dabei das Runge-Kutta-Verfahren etwa 120 000 Werte berechnen - es kann also einige Minuten dauern. Wir rechnen deshalb hier nur mit 5 Schritten. Sie können das Ergebnis verbessern indem Sie den letzten Parameter von regression_l in Zeile 3 der folgenden Zelle erhöhen.

Mit der Funktion g_regression_l können wir versuchen, die in einem Regressionsschritt gefundenen Parameter in weiteren Regressionsschritten feiner zu variieren. Wie oft dies erfolgt wird durch den Parameter depth festgelegt.

Wir können unterschiedliche Strategien zur Reduzierung der Rechenzeit anwenden:

Mit $gran=8, depth=2$ rechnen wir in jedem Regressionsdurchlauf genauer, brauchen aber weniger Durchläufe als mit $gran=2,depth=8$. Im ersteren Fall sind etwa 131 000 Datenpunkte zu berechnen, im letzteren Fall nur ca. 19 000. Bildlich gesprochen fischen wir im ersteren Fall mit einem feineren Netz und im letzteren Fall dafür öfter.

Fazit

Durch logistische Regression haben wir die Zahl der COVID-19-Infizierten (hochgerechnet) und Toten in Deutschland im April 2020 durch SI-Modelle mit einer relativen Genauigkeit von 3% darstellen können. Für das SIR-Modell, das die Zahlen der infizierbaren, Infizierten und Genesenen bzw. Verstorbenen berechnet, erreichten wir mit numerischen Lösungen der SIR-Differentialgleichungen nach dem Runge-Kutta-Verfahren eine relative Genauigkeit von 10%. Für die Regression wurden Algorithmn verwendet, die einen Anfangsbereich systematisch nach Werten mit einem minimalen Residuum absuchen. Es wurden Wege diskutiert um die Genauigkeit der Näherung zu steigern oder die benötigte Rechenzeit zu senken. Der Leser wird ermutigt, im Notebook selbst mit diesen Algorithemn zu experimentieren um bessere Näherungen zu finden.