BLog |
 |
|
|||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
Modelle zur Ausbreitung des n-Coronavirus-2019
In der mathematischen Epidemiologie wurden verschiedene Modelle entwickelt, um den dynamischen Verlauf der Ausbreitung von Krankheiten zu beschreiben. Als Basis gilt das SI-Modell. Dabei handelt es sich um einen Satz von zwei gewöhnlichen Differentialgleichungen, die die zeitlichen Verläufe der Anzahl der gesunden (Suszipierfähigen) und der Anzahl der erkrankten (Infektiösen) Individuen in einer Population beschreibt. Das SI-Modell ist Was ist der Zweck der Übung?Als die Ausbreitung in China in der dritten Januar-Woche richtig an Fahrt aufnahm, gab es schnell Stimmen, die das anfängliche exponentielle Wachstum ad infinitum bis hin zur Auslöschung der Menschheit in Jahresfrist extrapolierten. Jetzt wissen wir ja von Albert Einstein, daß nur zwei Dinge unendlich sind, nämlich das Universum und die menschliche Dummheit, und weiter, daß er sich beim ersten nicht wirklich sicher ist. Ich war gerade Anfang Januar wieder glücklicher Vater geworden und wollte natürlich wissen, ob das Schicksal meiner Kinder schon besiegelt ist, oder ob ihnen trotz allem noch eine glückliche Zukunft offen steht. Als mathematisch Gebildeter war mir natürlich sofort klar, daß die Fortschreibung der Exponentialfunktion einfach Quark war und ist, denn in der realen Welt gibt es nunmal die logistische Begrenzung und die macht sich immer schon relativ früh in den jeweiligen Kurvenverläufen bemerkbar, zunächst unscheinbar - man muß schon genauer hinsehen, dann aber soweit, daß nur unendliche Dummheit jemanden davon abhalten könnte, diese zu sehen. Die für uns Bürger spannenden Fragen waren, sind und bleiben, welches Ausmaß die Corona-Ausbreitung annehmen wird, und wie lange es dauern könnte. Was das voraussichtliche Ausmaß angeht, wollen wir uns mit der Größenordnung begnügen, und zur Dauer der Corona-Welle sollte uns reichen zu wissen, ob es Monate oder Jahre dauert bis man sein Leben wieder in einen halbwegs normalen Gang bringen kann. Eine Information wie, die Infektionsrate in D verdoppelt sich nur noch alle 9 Tage anstelle von 5 Tagen in der Vorwoche, und jeder Infizierte infiziert im Schnitt nur noch eine weitere Person, ist zwar beruhigend, aber noch nicht wirklich das was man wissen will. Dank Internet stehen uns heutzutage Informationen und Daten in Hülle und Fülle zur Verfügung, an denen wir unser mathematisches Geschick erproben können. Da wären die WHO-Situationsberichte, und die Daten für die China-Zeitreihe hatte ich noch mühsam händisch daraus zusammengetragen. Mittlerweile unterhält das Center for Systems Science and Engineering (CSSE) an der Johns Hopkins University in ihrem GitHub-Repository tagesaktualisierte Daten zur Ausbreitung von Covid-19 in vielen Ländern und Regionen, und man kann sich die Daten als CSV-Tabelle herunterladen. Dafür habe ich ein Tool programmiert mit dem man die Zeitreihen für spezifizierte Länder bzw. Regionen automatisch extrahieren und in eine spaltenorientierte TSV-Datei exportieren kann, und ich habe dann an den Daten für Deutschland, Italien und Brasilien das logistische Wachstumsgesetz erprobt. Es zeigt sich nun, daß der berichtete Verlauf der Corona-Ausbreitung in den genannten Ländern zwar grob mit der Logistischen Funktion in Einklang gebracht werden kann, aber wenn man genauer hinsieht, gibt es systematische Abweichungen, und diese werden bei der Fortschreibung mit jedem Berichtstag größer, so daß nur unendliche Dummheit uns davon abhalten könnte, diese nicht zu sehen. Die einfache Logistische Funktion stößt an die Grenzen ihrer AnwendbarkeitIn der linearen Auftragung zeigen sich von Tag zu Tag zunehmende systematische Abweichungen im unteren Bogen, und die Lage des Wendepunktes der Anpassung ist unsicher. Ferner liegt die angepaßte Kurve bei den letzten Punkten ständig zu niedrig, was eine Extrapolation fragwürdig macht. 
In der halblogarithmischen Auftragung erkennt man einerseits die beim Fortschreiben der Kurvenanpassung von Tag zu Tag zunehmenden Abweichungen im unteren Bereich und ferner einen anderen merkwürdigen Effekt speziell in der Deutschland-Zeitreihe. 
Der Verlauf der berichteten bestätigten Fallzahlen läßt sich zwanglos in Segmente gruppieren. Sehen wir hier möglicherweise den Einfluß unterschiedlicher Test-Regime mit schrittweise erhöhter Schlagzahl? Jedenfalls muß man davon ausgehen, daß das Coronavirus keine 8-Tage-Woche kennt. So oder so fängt die Logistische Funktion sowas natürlich nicht ein. Die Dunkelziffer
Hier begegnen wir mal wieder der berühmten Dunkelziffer, die uns immer wieder an den Kopf geworfen wird. Angefangen von
Selbstverständlich werden je nach Test-Regime nicht alle Infizierten erfaßt, und in einer solcherart dynamischen Situation hinkt man selbstverständlich mit der logistischen Herausforderung hinterher, den Bedarfsstellen eine genügende Anzahl Test-Kits bereitzustellen. Die Klugscheißerfrage,
Bei der Zeitreihenanalyse des Corona-Ausbruchs in China, hatte ich die Logistische Funktion mit einem zeitlich lineareren Glied + Versatz modifiziert. Die Annahme dazu ist, daß sich die Dunkelziffer ganz grob in demselben Ausmaß entwickelt wie die gemeldeten Fallzahlen, und damit wäre sie also inhärent im Verlauf über einen zwar unbekannten wenngleich Zurück an die ModellwerkbankWas ist so toll an der Logistischen Funktion? Nun, sie ist die Basis von unzähligen Wachstumsverläufen in der Natur, und so auch in der Epidemiologie. Nur heißt sie dort anders, nämlich SI-Modell. Das SI-ModellWie bereits einleitend gesagt, beschreiben die beiden gewöhnlichen Differentialgleichungen die zeitlichen Verläufe der Anzahl der gesunden (Suszipierfähigen) und der Anzahl der erkrankten (Infektiösen) Individuen in einer Population.
Hierin sind S und I strickt gegenläufig und deren Summe ist konstant: S + I = y0 + y1 = a1 + a2 = c. Unter Verwendung der Konstanten c und mit den Setzungen L = I, y = y1 und a1 = c vereinfacht sich das SI-Modell zur einfachen Logistischen Differentialgleichung (LDE):
Diese wiederum hat eine analytische Lösung, nämlich die Logistische Funktion (LF):
Die Modelle SI, LDE und LF sind also vom Prinzip her identisch, und an und für sich müssen wir uns damit nicht mehr länger aufhalten. Nur noch soviel, der Grenzwert in den L bei der Kurvenanpassung der gemeldeten Fallzahlen einmündet ist gleich S0, also die Anzahl der suszipierfähigen UND Jetzt ist es aber tatsächlich so, daß ein Großteil der erkannten Fälle mit einer gewissen Latenz aus der Infektionskette ausscheidet (engl. Remove), und zwar durch Isolierung, Gesundung und Immunisierung oder unglücklicherweise auch durch den Tod. Das SI-Modell wird um die dritte Differentialgleichung R erweitert und man kommt zum SIR-Modell. Das SIR-Modell
Bei diesem Modell wird der Verlauf von R (= y2) über die Parameter a0...a5 an die gemeldeten Fallzahlen angepaßt. Unter der Voraussetzung, daß die allermeisten infektiösen Fälle (d.h. die schwer erkrankten Infizierten) nach einer Latenz zum Melde- und dann zum Isolierfall werden, hätte man sich auf elegante Weise der Dunkelziffer entledigt, denn die bleibt in S und I stecken und sollte jedenfalls bei der Kurvenanpassung nicht weiter stören. Der Parameter a3 ist der Kehrwert der Latenzzeit, und diese läßt sich laut RKI sogar gut abschätzen, jedenfalls geht das RKI bei seinen Modellen von 2,6 Tagen aus, die vergehen bis ein infektiöser Fall aus der Infektionskette herausgenommen werden kann.
Um Kurvenanpassungen von Systemen mit gewöhnlichen Differentialgleichungen an die CSSE@JHU-Covid-19-Zeitreihen zu ermöglichen, habe ich nun in mein Kommandozeilen-Programm
# xcssecovid -i -m SIR -f012 -a0 0.6 -a3 0.38 -a4 8.5 -z 75 Germany time_series_covid19_confirmed_global.csv - # Germany # Model: SIR Differential Equations # S dy0/dt = -a0/a1·y0·y1 || y0(a5) = a1-a2-a4 # I dy1/dt = a0/a1·y0·y1 - a3·y1 || y1(a5) = a2 <- a4/a3 # R dy2/dt = a3·y1 || y2(a5) = a4 # # a0 = 0.59163 ± 0.848 % # a1 = 217146 ± 3.2512 % # a2 = 72.9125 ± 10.772 % # a3 = 0.38 # a4 = 8.5 # a5 = 35 # # ChiSqr = 938988.7
Die Option Die Gesamtzahl der betroffenen Individuen in Deutschland betrüge demnach a1 + a2 + a4 = 217145. Die Anzahl, die sich davon bei Erreichen des Grenzwertes infiziert haben würden und gemeldet worden wären, nämlich 134228 muß man aus den Kurven ablesen, und diese SIR-Kurven sehen nun schon erheblich brauchbarer aus als die vergleichbaren LF(SI)-Kurven weiter oben.
Wir sehen allerdings noch eine systematische Abweichung im unteren Bogen, bei den Fallzahlen bis ca. 15000. Dennoch ist die Abweichung im oberen Bereich verschwunden, und damit trauen wir uns nun wieder eine Vorhersage zu. Bei Fortsetzung des derzeitigen Eindämmungsregimes sollte eigentlich das Schlimmste in 3-4 Wochen überstanden worden sein. Die Kurvenanpassung muß jedoch in den nächsten Tagen fortgeschrieben werden, damit aktuelle Entwicklungen in die Vorhersage einfließen können. Das ERF-ModellDabei handelt es sich um ein Modell, das vom Institute for Health Metrics and Evaluation (IHME) an der University of Washington State in Seattle zur Vorhersage der Auswirkungen der Coronavirus-Ausbrüche auf die Gesundheitssysteme in den einzelnen Staaten der USA und weiteren Ländern verwendet wird. In einem Fachartikel wurde das Modell offengelegt, nämlich eine verschobene Fehlerfunktion (Error Function - ERF). Die Fehlerfunktion ist so wie die Logistische Funktion eine S-Funktion, sie ist allerdings steiler. Die Zielmenge ist -1 bis +1 bei Funktionsargumenten von -∞ bis +∞. Zur Verwendung als Modellfunktion zur Beschreibung der zeitlichen Verläufe der Corona-Todesfälle in den USA hat das IMHE sie um +1 nach oben und um einen variablen Parameter auf der Zeitachse nach rechts verschoben.
Die verschobene Fehlerfunktion habe ich auch als Modell in mein Tool
# xcssecovid -i -m ERF -z 75 Germany time_series_covid19_confirmed_global.csv - # Germany # Model: Shifted Error Function # y = 0.5·a1·(1 + erf(a0·(t - a2))) # # a0 = 0.0797187 ± 2.0207 % # a1 = 132419 ± 2.2802 % # a2 = 68.8025 ± 0.41395 % # # ChiSqr = 621844.3 Selbst wenn die Fehlerfunktion ein Ad-hoc-Modell zu sein scheint, muß man ihr doch folgendes lassen, sie paßt, wackelt und hat Luft. Die Anpassung gelingt mit einem erheblich geringerem 𝜒2 von nur 2/3 der vergleichbaren Anpassung des SIR-Modells. Die vorhergesagten Grenzwerte sind 132419 zu 134228 quasi identisch. Der Wendepunkt befindet sich beim Tag 68, und weil die Fehlerfunktion um den Wendepunkt herum symmetrisch ist, wäre der Auslauf vorhersagbarer.
Es bleibt ein gewisses Unbehagen, denn eine Funktion, die einfach so mal aus einer Vielzahl möglicher S-Funktionen ausgewählt wurde weil sie gut paßt, ist eigentlich kein Modell. Die Fehlerfunktion kommt unter gewissen Randbedingungen in der Wärmeleitungsgleichung vor. Die Ausbreitung des Coronavirus in Analogie zur Wärmeleitung? Da käme uns ja die Fermi-Dirac-Statistik als passenderes Analogon in den Sinn, und deren Verteilungsfunktion ist auch eine sigmoide, nämlich die Grundform der Logistischen Funktion, da ist sie ja schon wieder diese LF. Das GLF-Modell
Das SI-Modell (= LF-Modell, Logistische Funktion) und das verfeinerte SIR-Modell haben den großen Vorteil, daß sie die zugrundeliegenden Vorgänge der Virus-Ausbreitung erklären. Die Ergebnisse mit dem ERF-Modell zeigen uns allerdings, daß die LF nicht steil genug ist. Moment mal, da gibt es die Verallgemeinerte Logistische Funktion, bei der die Steilheit mit dem Exponenten 1/𝜈 eingestellt werden kann, und wir kämen damit wieder auf einen theoretisch fundierteren Boden. Ich habe auch die GLF als Modell in
Durch Ausprobieren findet man, daß man mit der Setzung von a3 = 0,3 für die Deutschland-Zeitreihe die beste Kurvenanpassung über die anderen Parameter a0...a2 erhält:
# xcssecovid -i -m GLF -a3 0.3 -z 75 Germany time_series_covid19_confirmed_global.csv - # Germany # Model: Generalised Logistic Function # y = a1/(1 + exp(-a0·a3·(t - a2)))^(1/a3) # # a0 = 0.381885 ± 2.2646 % # a1 = 159938 ± 2.6699 % # a2 = 58.3897 ± 0.21511 % # a3 = 0.3 # # ChiSqr = 609164.1 Die Anpassung gelingt mit einem nochmals geringerem 𝜒2 im Vergleich zur Anpassung des ERF-Modells. Der Grenzwert liegt mit kanpp 160000 Fällen jedoch deutlich höher, und es wird wohl auch noch etwas länger dauern bis er sich einstellt.
Der Exponent 1/a3 = 1/0,3 = 3,33... läßt sich so verstehen, daß sich im zugrundeliegenden Verlauf mehrere logistische Wirkungen potenzieren. Darin steckt die mit der Zeit zunehmende Testhäufigkeit inklusive Meldeverzögerung, darin steckt die zunehmend raschere Erkennung und Isolierung von Infizierten. Darin steckt das Social Distancing, der Public Shutdown bis hin zum Lockdown, und schließlich auch die unterschiedlichen Stadien und Anstiege der Infektionswellen in den Regionen. Der dem Verlauf inhärente Grenzwert ist kleiner und wird schneller erreicht und dazu ist die Kurve am Wendepunkt steiler, vorher stärker und im Auslauf schwächer gekrümmt. Das SEIR-ModellIn diesem Model wird berücksichtigt, daß nach der Infektion eines Individuums ein gewisser Zeitraum vergeht, bis es selber andere anstecken kann. Individuen in diesem Stadium waren also bereits dem Virus ausgesetzt (exponiert), sind aber selber noch nicht infektiös. In das SIR-Modell wird also eine weitere Differentialgleichung für die Exponierten eingefügt, und man erhält das SEIR-Modell.
# xcssecovid -i -m SEIR -z 75 Germany time_series_covid19_confirmed_global.csv - # Germany # Model: SEIR Differential Equations # S dy0/dt = -a0/a1·y0·y2 + a8/y0 || y0(a7) = a1-a2-a5-a6 # E dy1/dt = a0/a1·y0·y2 - a3·y1 || y1(a7) = a2 <- a6/a4/a3 # I dy2/dt = a3·y1 - a4·y2 || y2(a7) = a5 <- (1 - a4)·a3·a2 # R dy3/dt = a4·y2 || y3(a7) = a6 # # a0 = 0.510309 ± 2.6952 % # a1 = 157426 ± 1.9613 % # a2 = 334.427 ± 11.125 % # a3 = 0.4 # a4 = 0.1 # a5 = 120.394 # a6 = 17 # a7 = 35 # a8 = 0 # # ChiSqr = 773805.1 
Damit ergibt sich eine sehr gute Kurvenanpassungen für den Fall-Verlauf in Deutschland, fast gleichwertig zu der mit der Generalisierten Logistischen Funktion. Zu den Parametern (s. auch SEIR and SEIRS models):
Die Parameter mit den %-Angaben sind Variablen in der Kurvenanpassung, für die anderen müssen
Mit dem Parameter a1 = Der Punkt ist, daß im großen statistischen Schnitt nur ein Bruchteil der Gesamtbevölkerung, nämlich die virtuelle anfällige Population, überhaupt vom Virus getroffen werden kann, und der Rest befindet sich in einer mehr oder weniger guten Deckung. Nach dem vorliegenden Ergebnis befinden sich im Schnitt ca. 160000 Individuen ausserhalb der Deckung, also sozusagen auf dem Präsentierteller. Im Grunde kann jeder vom Virus getroffen werden, nur die Unvorsichtigen und Dummen immer zuerst. Auch in einem guten Western, trifft es immer diejenigen zuerst, die sich für unverwundbar halten und ihren Kopf partout in der Schußlinie exponieren wollen. Wie gesagt, das ist nur meine Meinung. Ein a0 von 0,51 ist die Wahrscheinlichkeit (51 %) mit der ein Infektiöser auf eine noch gesunde Person aus der virtuellen Population (die ausserhalb der Deckung) von ca. 160000 trifft und ansteckt. a2 = 334 war die Anzahl der vom Corona-Virus getroffenen (die Exponierten) aber noch nicht infektiösen Individuen am Tag a7 = 35. a3 = 0,4 ist der Kehrwert der Latenzzeit (manchmal auch ungenau als Inkubationsrate bezeichnet), und beim vorliegenden Modell bedeutet diese, daß es im Schnitt 1/0,4 = 2,5 Tage dauert bis ein Corona-Exponierter infektiös wird. a4 = 0,1 ist die Rate mit der infektiöse Individuen aus der Infektionskette ausscheiden, sei es durch Isolierung, Gesundung oder unglücklicherweise den Tod - jedenfalls sie sind Raus. a4 und a3 kann man wegen Parameterkorrelation bei der Kurvenanpassung nicht zugleich variabel halten. Durch Ausprobieren gelingt eine Minimierung der Abweichungen bei der 7-Tage-Simulationsretrospektive bei einer Ausscheiderate von 0,1, d.h. es dauert im Schnitt ca. 10 Tage. Hier habe ich gewisse Zweifel ob das stimmen kann. Tatsache ist, daß die Kurvenanpassung mit 0,1 die besten Resultate ergibt. a5 und a6 sind einfach nur Randwerte, die sich mittelbar aus a2, a3 und a4 ergeben. Wenn man sie auf 0 setzt, ändert das so gut wie nichts am Ergebnis, die Kurvenanpassung dauert aber viel länger (mehr Iterationen). a7 ist der Tag ab dem es mindestens 17 gemeldete Infektionsfälle gab. In Deutschland war das der Tag 35, der Fastnachtsdienstag, 25.2.2020, und am Aschermittwoch ging es dann richtig los. a8/y0 ist eine Art Diffusionsrate. a8 ist der Faktor aus einer Diffusionskonstanten d und der Gesamtbevölkeung TP, nämlich a8 = d·TP. Die Rate mit der Individuen entlang des Konzentrationsgefälles aus der Gesamtbevölkeung (TP) in den virtuellen Anteil der anfälligen Population (VSP) diffundieren ergibt sich damit zu a8/VSP = d·TP/VSP. FazitDie besten Kurvenanpassungen gelingen mit der Generalisierten Logistischen Funktion dicht gefolgt vom SEIR-Modell. Am SEIR-Modell lassen sich jedoch die Vorgänge in der Infektionskette besser erklären. Selbst wenn wir in den verschiedenen Ländern die jeweils ersten Infektionswellen überstanden haben, können weitere folgen. Wir müssen solange es keine Impfung gibt vorsichtig und in Deckung bleiben. Copyright © Dr. Rolf Jansen - 2020-04-05 13:11:25 Diskussion auf Twitter: 1247543331418071040 |
||