Sonntag, 8. März 2020

Coronavirus: Dynamik der Ausbreitung anhand des Standardmodells

Übersicht

Das Coronavirus beschäftigt seit Mitte Januar 2020 breite Bevölkerungsteile, besonders seit es immer näher kommt und in der Schweiz am 5. März die erste infizierte Tote gegeben hat. Im Dezember stellte man in China eine Häufung von Lungenentzündungen fest, die von einem neuartigen Virus verursacht schienen. Dieses nennt sich SARS-CoV-2, wobei SARS für "severe acute respiratory syndrome" steht und die Nummer 2 es vom Erreger der Epidemie von 2003 unterscheidet. Die Krankheit wiederum wird Covid-19 genannt.

Die vier Haupt-Klassen von Erregern, neben Pilzen, Milben und Prionen, sind:
  • Viren,
  • Bakterien,
  • Einzeller (Protozooen) und
  • Würmer (Helminthen).

Die Infektionsübertragung folgt verschiedenen Wegen, hier auch vier:

  1. Person --> Person (Grippe, HIV etc.),
  2. Person --> Umwelt --> Person (Typhus, Cholera)
  3. Reservoir --> Träger --> Person (Pest, Dengue, Malaria etc.)
  4. Reservoir --> Person (Tollwut, Brucellose etc.).

Als Laie würde man die Wege 1 und 4 als die bisher gefährlichsten ansehen, mit Malaria sicher als der hartnäckigsten und über die Zeit schädlichsten Infektion. Dies hängt mit der Häufigkeit des Trägers (Vektor) und seiner erfolgreichen Übertragung zusammen; Mücken gibt es fast unzählige.


Einer starke Vermutung nach soll das Virus von Fledermäusen auf den Menschen übertragen worden sein (Weg 3) und nun auch den gefährlichen Übertragungsweg von Mensch zu Mensch beschreiten (Weg 1). Der Übergang von Tierkrankheit zum Befall von Menschen nennt man Zoonose. Die meisten neunen Infektionen folgen diesem Muster, d.h. Ebola, West-Nil, Hanta, HIV usw.

Das Modell

Wir folgen für das Modell der Darstellung von Hethcote (1989) und nehme für uns keine besonderen Meriten in Anspruch ausser der Tatsache, die Lösungsfindung Schritt für Schritt darzustellen, ohne "wie man leicht nachvollziehen kann", "der Leser sofort erkennt", "allgemein bekannt ist" etc. Man kann diesen Blog auch als Auffrischung von Analysis II lesen.

Es gibt drei klassische Modelle aus der theoretischen Biologie, die sich als Zustandsdiagramme darstellen lassen. Die Zustände beschreiben den prozentualen Anteil der Bevölkerung (oder Kontrollgruppe), entweder ansteckbar und damit gesund sind (engl. susceptible S), ansteckend sind (infective I) und sich erholt und immunisiert haben (recoverd R). Die Infektiösen sind manifeste oder latente Befallene. Diese Modelle sind mit den Zustände S, I und R:
  1. SI: Ansteckung ohne Gesundung,
  2. SIS: Ansteckung  ohne Immunitätsbildung,
  3. SIR: Ansteckung  mit Immunitätsbildung.
Im Lauf der Zeit wurde noch ein weiterer Zustand in die Modelle aufgenommen, nämlich die "Exposed": Elemente dieser Gruppe sind bereits infiziert aber noch nicht infektiös. Diese Modelle sind dann SEIS, SEIR etc. 

Soweit man heute weiss, folgt das endemische Corona-Virus (SARS-CoV-2) leider dem zweiten Typus, denn es sind Wiederansteckungen beobachtet worden. Damit wird eine passive Immunisierbarkeit unwahrscheinlich. Wir gehen deshalb speziell auf das SIS-Modell ein.
Abb. 1: Zustandsdiagramme für SIR (oben) und SIS (unten). Die demographischen Variablen von Geburt und Sterblichkeit bei den Gesunden kann man auch ausblenden.



Die Gleichungen

Das SIR-Modell ist das komplexeste, weil es am meisten Varaiblen, nämlich deren drei, aufweist. SIS und SI haben zwei Variablen, die durch die Schliessbedingung, z.B. S+I+R=1, miteinander gekoppelt sind.

In Abb. 2 ist das SIS-Modell nochmals mit Parametern dargestellt. Betrachtet man das Kästchen S(t), dann kann man die Änderung des Inhalts S(t) als Bilanz von Zu- und Abfluss festhalten. Zufliessen ersten alle Neugeborenen, die als nicht-infektiös betrachtet werden und der Anteil γ der Infektiösen, die sich wieder erholt haben. Abfliessen die natürlichen Toten von S(t) mit der Sterberate μ, die wir mit der Geburtenrate gleichsetzen. Die Abflüsse durch Ansteckung bestimmen sich durch eine Ansteckungsrate λ, die auf das Produkt von Ansteckbaren S(t) und Infektiösen I(t) wirkt. Motiviert wird dieses Produkt durch die Feststellung, dass wenige Ansteckungen erfolgen, wenn wenige Infektiöse vorhanden sind, genau gleich, wie wenn fast alle schon angesteckt sind. Bei einer Aufteilung von 50/50 ist die Ansteckung am höchsten.

Die Änderung der Infektiösen in einem Zeitintervall dt ergibt sich ebenfalls aus Zu- und Abflüsse. Daraus kann man die entsprechenden Formeln ableiten.

Abb. 2: Das SIS-Modell mit seinen Parametern.
Formal lauten die Formel für die Änderungen dS/dt und dI/dt wie folgt:
Mit der Schliessungsbedingung (Gl. 3) kann man je eine Variable eliminieren, so dass alternativ folgende Darstellung gültig ist:




Im weiteren betrachten wir nur Gleichung 5; denn wenn I(t) bekannt ist, ist auch S(t)=1-I(t) gegeben. Es handelt sich hier um eine sogenannte gewöhnliche Differentialgleichung, bei der zu einer gesuchten Funktion I(t) nur Ableitungen nach genau einer Variablen t auftreten. Man sieht, dass zudem Summen von zwei Potenzen der gesuchten Funktion vorliegen. Glücklicherweise bildet unsere Gleichung eine Ausprägung der Bernoulli'schen Differentialgleichungen, für die wir eine geschlossene Lösunge herleiten. 

Sollte der Term (λ-γ-μ) null sein, dann folgt aus Gl. 5 nach ummöblieren:

mit der Anfangsbedingung I(0)=I0 ziemlich einfach.

Betrachten wir nun den allgemeineren Fall, d.h. λ<>γ+μ. Um der Übersichtlichkeit willen verzichten wir auf den Hinweis, dass die Variablen z.T. von der Zeit t abhängen. Wir werden mit zwei aufeinanderfolgenden Variablentransformationen die Differentialgleichung soweit vereinfachen, dass wir eine einfache, bekannte Lösung finden werden.

Die erste Transformation lautet I=1/Z. Das Differential muss man auch transformieren (Gl. 11). Damit, und der Vereinfachung der Parameter, ergibt sich die lineare Dgl. nach Gl. 12, formal:




Wir setzen nun eine weitere Transformation an und zwar wie folgt:
Die DGL nach Gleichung 17 besitzt die offensichtliche Lösung log(G), wobei der log der natürliche Logarithmus ist (aus Tabellen: Stammfunktion von 1/x ist log|x|). G bekommt man durch Exponentieren, sodann setzt man die Rücktransformationen ein und erhält die Funktion I(t). 
Die  Integrationskonstante c2 bestimmt man durch die Anfangswerte I(t=0)=I0. In Gl. 25 ist die Lösung mit allen Parametern. 
Nun haben wir die Formel für den Anteil der Infektiösen an der Gesamtbevölkerung. Allerdings sind vier Paramter empirisch zu bestimmen, nämlich λ, γ, μ und I0.

Die Parameter

Im Zustand (oder Reservoir) I(t) sind Neuansteckungen, Erholungsrate und Sterbereate pro Zeiteinheit sowie der Anfangswert von I zum Betrachtungszeitpunkt zu schätzen. Uns interessiert vor allem die Mortalität durch das Virus also M(t)=μ I(t).

Die von Hethcote (1989) verwendete Terminologie ist nicht die einzige. Heute wird viel von der Basic Reproduction Number R0 gesprochen, für die Schätzungen für verschiedene Infektionen vorliegen (Wikipedia contributors, 2020). Wie man aus Gl. 26 sieht, hat die Reproduktion mit dem Verhältnis von Kontakten zu Genesungen plus Todesfälle zu tun.
Für R0 wird ein Wert von 1.4 bis 3.8 geschätzt. Für die SARS-Infektion von 2003 gilt ein mittlerer Wert von 3.5 (Wikipedia contributors, 2020). Die saisonale Grippe liegt bei 1.3.

Das Gamma lässt sich als Inverse der Ansteckungsdauer verstehen. Bei einer Dauer von 14 Tagen nimmt man Gamma=1/14. Die Mortalität ist schätzungsweise (neu) 1.6%, wobei dieser Durchschnitt sehr inhomogene Sub-Gruppen zusamenfasst. Bei den Risikogruppen ist der Wert möglicherweise substantiell höher, bei Kindern etwa wesentlich tiefer. Diese Inhomogenität wird in diesem einfachen Modell nicht berücksichtigt. Bei der vorherigen SARS-Infektion schätz man die Mortalität auf 9.6%, bei MERS auf 34% und der porcinen Influenza auf 0.02%.

Der heutige Wert von I0 ist wegen der mögliche Dunkelziffer schwer zu schätzen. Wir nehmen einfach per 5.3.2020 an, es seien die publizierten Zahlen anwendbar, also ein Anteil von 0.0000265 positiv Getesteten. Er ändert nichts am prizipiellen Verlauf der Kurve, hat aber einen Einfluss, wann der starke Anstieg erfolgt. Je kleiner der Wert desto später der Anstieg.

Abb. 3: Verlaufe des Anteils Infektiösen mit verschiedenen Parameterschätzungen.

Die Unsicherheit

Die Paramter unterliegen grossere Unsicherheit. Wie man der Abb. 3 entnehmen kann, ist das Sigma oder R0 der wichtigste in diesem Modell. Je höher dieser Wert liegt, desto schneller und desto höher der Anteil der Infizierten. 

Die Parameter γ und μ wirken genau gleich, was man aus der Formel ersieht. Eine Verringerung dieser Werte führt zu einem steileren Anstieg und damit schnelleren Erreichen des asyptotischen Wertes. Sigma, oder R0, bestimmen den zu erreichenden Wert alleine. 

Sigma ist aber keine Naturkonstante des Virus sondern hängt wesentlich vom Verhalten der Menschen ab. Die Faktoren sind die Häufigkeit und die Qualität der Kontakte. Je weniger und je vorsichtiger desto besser. Weniger Kontakte bedeutet Ansammlungen zu vermeiden, Sicherheitsabstände einhalten und Berührungen meiden. Qualität bedeutet Emission von Tröpfchen kontrollieren und Hände gut waschen. 

Die theoretischen Verläufe ziehen die Möglichkeit von immer wirksameren Gegenmitteln nicht in Betracht. Die Wirksamkeit ist eine Funktion der Suche, die wiederum von der Zeit abhängt. Aus diesem Grund ist es ebenfalls dringlich, den Anstieg der Infektion so weit wie möglich in die Zukunft zu verschieben. Die Abb. 3. allerdings zeigt Zeiträume für Sigma=1.6, die für einen zu entwickelnden Impfschutz nicht ausreichen.

Modelle zeigen nur einen Auschnitt der Realität, die in Unkenntnis von vielen Faktoren potentiell wichtige ausser Acht lassen. Zudem zeigt Abb. 3, dass Paramtervariationen grosse Unterschiede verursachen. Wie Silver (2012, 214) schon bezüglich der Grippe von 2009 bemerkte, können Voraussagen ziemlich ungenau sein. Bemerkenswert ist die Tatsache, dass die geringe Anzahl Toter in der vermeintlichen Zone des Epidemieausbruchs nicht gut mit den obigen Parametern zusammenpassen (ausser die drakonischen Massnahmen haben so stark gewirkt). 

Im Nachhinein werden wir sicher klüger sein und Modell und Parameter auf seine relative Gültikeit beurteilen können.

Dennoch zeigt das Modell, dass die behördlicherseits getroffenen Vorkehrungen und Ermahnungen sicher richtig und wirksam sind. Möglicherweise sind sie aber zu zögerlich.

Erweiterung des Modells

Wir wollen die Situation untersuchen, bei der sich das Sigma von 3.8, sehr hohe Infektionsrate, auf 1.3 senkt. Dies sollte die getroffenen Massnahmen von geringerem Kontakt und grösserer Vorsicht, Händewaschen etc. darstellen. Dabei variieren wir die Geschwindigkeit der Senkung. Diese zwei Varianten stellen die Konsequenz der Umsetzung der Massnahmen dar.

Kucharski et al. (2020) haben R0-Werte anhand der Daten von Wuhan geschätzt. Sie gehen von einem R0 vor Reisebeschränkungen von 2.35 aus, das sich zwei Wochen später auf 1.05 gesenkt hat. Allerdings sind die 95%-Konfidenzgrenzen sehr hoch, was die Unsicherheit und die geringe Stichprobe spiegelt. Für 2.35 sind es 1.15 bis 4.77 und danach 0.41 bis 2.39. Unsere grobe Schätzungen liegen in den Intervallen. Auch die Zeitliche Dimension von 2 Wochen ist einigermassen abgebildet. Für den weiteren Verlauf ist wesentlich, ob sich das R0 unter 1 senkt, denn nur dann klingt die Infektion ab.

Die Funktion des Sigmas/R0 wählen wir aus der Familie der fallenden Funktionen. Da wir die logistische Funktion schon gebraucht haben (Gl. 24 und 25), nehmen wir etwas wie eine atanh-Funktion ("Arcus tangens hyperbolicus") ohne weitere Begründung ausser der Tatsache, dass wir sie so etwa von Hand gezeichnet hätten. Speziell also:
wobei p=5 (schnelles) oder 10 (langsames Abklingen) genommen wird. Graphisch sehen die Verläufe wie folgt aus:

Abb. 4: Verläufe von Sigma
Nun nehmen wir die Gleichung 12 zur Hand und ersetzen λ=σ(μ+γ), wobei μ=0.016 und γ=1/14 ist. Anstatt eine sicher sehr langwierige analytische Rechnung zu probieren, nehmen wir das klassische numerische Verfahren von Runge-Kutta mit 4 Stufen. Damit berechnen wir die schnelle und die langsame Variante Z(t), wobei wir dann Gl. 10 zur Bestimmung von I(t) anwenden. Sodann vergleichen wir das Resultat mit der Infektionsrate nach Gl. 25, wenn man mit einem konstaten Sigma von 1.3 gerechnet hätte. Graphisch:


Abb. 5: Die unterschiedlichen Verläufe in Funktion der Schnelligkeit der Senkung der Infektionsrate.


Was man aus der Abb. 5 gut erkennt sind die grossen zeitlichen Unterschiede, die durch die anfänglich hohen Sigma-Werte erzeugt werden. Ein Unterschied von rund 10 Tagen bei der Senkung produzierte eine Differenz von rund 80 Tagen. Dies würde im realen Leben für ein sehr dezidiertes Vorgehen plädieren. Der Unterschied stammt von den ersten Tagen, die ausschlaggebend sind( siehe Abb. 6).

Abb. 6: Die ersten Tage der unterschiedlichen Infektionsverläufe.

Und nun die Mortalität, das heikle Thema. Wir gehen von einer Gesamtbevölkerungszahl der Schweiz von 8.57 Mio. aus. Aufgrund der Infektionsrate σ und der Mortalitätsrate μ von 1.6% ergeben sich die folgenden Schaubilder. Ein wichtiger Punkt ist die zeitliche Verzögerung. Angenommen die Einführung eines wirksamen Impfschutzes daure 270 Tage, dann sieht man, dass bei der schnellen Senkung des Sigmas/R0 viele davon profitieren können, während es bei der langsamen Senkung für viele zu spät kommt. Das Wichtigste modellmässig gesprochen, ist dσ(t)/dt zu minimieren.

Abb. 7: Mortalität für die zwei Szenarien mit fallender Infektionsrate
Abb. 8: Tägliche Mortalität für die zwei Szenarien mit fallender Infektionsrate (Bevölkerung Schweiz 8,57 Mio Einwohner)
Der Bundesrat Berset, u.a. Schweizer Gesundheitsminister, hat mit seinem Ausspruch "Man würde sagen, spinnt er oder was?" zur Frage des früheren Treffens von Massnahmen an der Pressekonferenz vom 13. März 2020 gezeigt, dass seine Experten die Dringlichkeit und die prinzipiellen Aussagen dieser Analysen nicht beherzigten oder nicht als politsch umsetzbar betrachtet haben. Die Reaktion hätte sein müssen: die rigorosen Massnahmen der Länder zu übernehmen, die zeitlich vorauslaufen. 

Postscript

"Es irrt der Mensch solang' er strebt" schreibt Goethe im Faust. In diesem Sinne kann ich Fehler in der Darstellung nicht ausschliessen und eo ipso facto nicht haftbar gemacht werden.

Daten

https://github.com/CSSEGISandData/COVID-19/blob/master/csse_covid_19_data/csse_covid_19_daily_reports/03-05-2020.csv


Literatur

Hethcote H.W. (1989) Three Basic Epidemiological Models. In: Levin S.A., Hallam T.G., Gross L.J. (eds) Applied Mathematical Ecology. Biomathematics, vol 18. Springer, Berlin, Heidelberg

Kucharski, Adam J., et al. "Early Dynamics of Transmission and Control of COVID-19: a Mathematical Modelling Study." The Lancet. Infectious Diseases, 2020.


Silver, Nate (2012) The signal and the noise : why so many predictions fail-- but some don't.  Penguin Press, New York.

Wikipedia contributors. (2020, March 6). Basic reproduction number. In Wikipedia, The Free Encyclopedia. Retrieved 13:44, March 7, 2020, from https://en.wikipedia.org/w/index.php?title=Basic_reproduction_number&oldid=944253318