Einfache Epidemie-Modelle

Bisher habe ich den Parameter R₀ in unüblicher Weise definiert. Beim RKI heißt es: "Die Reproduktionszahl beschreibt, wie viele Menschen eine infizierte Person im Mittel ansteckt." ...

... "Durch Infektionsschutzmaßnahmen lässt sich die Reproduktionszahl verringern. Man spricht von einer zeitabhängigen Reproduktionszahl R(t). Es gilt:

• Wenn R > 1, dann steigende Anzahl täglicher Neuinfektionen,
• Wenn R = 1, dann konstante Anzahl täglicher Neuinfektionen
• Wenn R < 1, dann sinkende Anzahl täglicher Neuinfektionen.

"

Ich möchte statt von R(t) von Rₜ sprechen. Ich würde es praktisch finden, diese Definition nur für den Anfang einer Epidemie gelten zu lassen. So ist sie ein Maß für den Erfolg die Gegenmaßnahmen. Die Anzahl der Neuerkrankungen geht am Ende auch ohne Maßnahmen zurück, wenn ein größerer Teil der Bevölkerung immun geworden ist.

Ich sage also jetzt Rₜ, wo ich früher R₀ gesagt habe. Außerdem gibt es noch eine Änderung. Um die Bedingung "Wenn R=1, dann ..." zu erfüllen musste ich einen Faktor b einbauen. Damit haben wir ein geändertes Modell, das "m/b" heißen soll im Unterschied zu "m/a", dem bisherigen Modell.

Modelldefition m/b

Das Modell ist nach wie vor durch das Flußdiagramm

characterisiert. Zur Bedeutung der Buchstaben S, E, I1, I2, Q, R, K siehe m/a

Das Modell ist definiert durch das folgende System von Differentialgleichungen

Hier ist nur der Faktor b neu. Er ist der Ausgleich dafür, dass ein durchschnittlicher Infizierter nicht über die ganze infektiöse Zeit andere anstecken kann. Er wird nämlich zu einem Anteil q nach 6 Tagen isoliert.

Parameter

• Für die Faktoren γ, die die durchschnittliche Verweildauer in den einzelnen Kompartimenten bedeutet, siehe m/a.
• Rₜ ändert sich mit der Zeit. Ich zerlegung die Zeit in 10 Tage breite Intervalle und vor und bestimme einen passenden Wert. Rₜ(t) ist also eine Treppenfunktion. Die Werte werden in einem automatisierten Verfahren durch Anpassung an die offiziellen Fallzahlen bestimmt.
• q hat am Anfang einer Epidemie (S ≈ N) gar keinen Einfluss. auf die Anpassung des Modells an die Daten. Bei hoher Dunkelziffer steigt die Zahl der Immunen schneller. Ich wähle q = 0,8.
• Anfang braucht das Modell eine Angabe über die Zahl der Infizierten. Ich habe Diese Zahl durch halbautomatische Anpassung an die Daten im März auf 1500 gesetzt.

Rₜ in Scenarios

Interesse an einem Modell wie m/b und an der Funktion Rₜ wird man deshalb haben, weil es eine Einschätzung der zukünftigen Entwicklung erlaubt. Wenn auch nicht sehr genau. Man kann davon ausgehen, daß es bei einer weitgehenden Lockerung zu einem schnellen Anstieg der Erkrankungen kommen würde.

Das Argument ist, dass die Krankheit sich schnell in der deutschen Gesellschaft ausbreiten kann, wie im März geschehen und dass der starke Rückgang der Neuerkrankungen im April auf die Quarantänemaßnahmen zurückzuführen ist. Die Werte von Rₜ, mit denen das Modell das vergangene Geschehen simulieren kann, können die Maßnahmen der entsprechende Zeit characterisieren. So kann man eine Vorstellung entwickeln, ob weitere Lockerungen möglich sind. Das Modell m/b ist ein minimalistisch einfaches Modell, doch es hat denselben Kern wie viele solche Modelle.

S′ = -β⋅S⋅I/N

wo β eine Konstante, S die Zahl der Ansteckbaren, I die Zahl der Ansteckenden und N die (konstante) Populationsgröße ist. Man kann argumentieren, dass es für einen Ansteckungsvorgang eine Person aus I braucht. Die hat nun im Mittel β Kontakte am Tag, mit irgendwelchen aus der Bevölkerung. Zu einem Anteil von S/N führen diese Kontakte zu einer Ansteckung. S′ ist aber das Negative der Zahl der Neuansteckungen pro Tag. In β verbirgt sich Rₜ. Es muß noch so normiert werden, dass über die gesamte Zeit der Infektiösität Rₜ die Zahl der Ansteckungen ist, die von einem Infizierten ausgeht.

Man stellt also die Frage: Wenn wir durch bestimmte Maßnahmen ein bestimmtes Rₜ für die Zukunft einstellen würden, was hätte das für Folgen?

Ich habe solche Szenarien mit dem m/b gerechnet.

Die Interpretation ist: Eine dauerhafte Überschreitung der Schwelle 1,0 für Rₜ wird später zu einer "Welle" an Erkrankungen führen, gegen die die "Welle" vom März nur ein Kräuseln war. Nach einer solchen Welle bleibt bei Rₜ=1,1 oder Rₜ=1,2 der Anteil der Immunen mäßig und die Maßnahmen müssen noch für viele Jahre aufrecht erhalten werden.

Die Abbildungen zeigen auch ungefähr stationäres Verhalten bei Rₜ=1.

Einmal kam eine Gruppe von 1000 Reisenden in Ihre Heimatprovinz zurück. Sie hatten sich alle auf dem Abschiedsball Mit Covid-19 angesteckt. In Ihrer Heimatprovinz war noch niemand erkrankt und noch niemand immun. Sie mischten sich unter das Volk und begannen die Grundbevölkerung zu infizieren. Die infektiöse Zeit der Einzelnen betrug 10 Tage.

Es gab zunächst 147 Neuansteckungen pro Tag. Aber nach einigen Tagen erkannte der Gesundheitsdienst sie und sie wurden isoliert. Das geschah im Mittel nach 6 Tagen und bei 80% der Gruppe. Also nach 6 Tagen blieben von den 1000 nur noch 200, die weiter für 4 Tage ansteckend waren. Im Effekt gab es jetzt pro Tag 147⋅0,2 ≈ 29 Ansteckungen. Wenn man über die infektiöse Zeit summiert

6⋅147 + 4⋅147⋅0,2 = 999,6 ≈ 1000

Jedes Gruppenmitglied hat also im Durchschnitt eine Neuinfektion erzeugt. Deshalb ist hier

R₀ = 1,0

Dies gilt unter der Voraussetzung, dass nach 6 Tagen 80% der Infizierten unter Quarantäne kommen. Wenn das wegfiele, hätte die Reisegruppe 1470 Neuinfektionen hervorgerufen und R₀ wäre 1,47.

Die Zahl 147 erhält man als 10/(6+4⋅0,2) ⋅ 0,1.

Die Situation der Reisegruppe ist eine besondere. Im allgemeinen hängt die Zahl der Neuerkrankungen noch von den Faktoren "Anzahl der Infizierenden" und "Anteil der Immunen in der Bevölkerung" ab. R₀ ist gerade so definiert, dass man von diesen Faktoren absieht. R₀ ist Ausdruck der Biologie des Erregers und der Kontaktfreude der Bevölkerung. Es kann sich mit der Zeit ändern.

Warum hatte Reisegruppe 1000 Mitglieder? Es schien grausam, von einem Patienten zu erzählen, der zu 80% in Quarantäne kommt.

Mein Modell hat bisher um den Faktor 1,47 zu hohe R₀ benutzt. Diese Zahlen gaben die Ansteckungsrate ohne Berücksichtigung der Quarantänemaßnahmen. Ich werde das durch den Einbau eines Normierungsfaktors korrigieren. Dann sollte im Anfang einer Epidemie bei R₀=1 eine konstante Zahl an Neuerkrankungen pro Tag vorhergesagt werden. Mir scheint das der Sprachgebrauch zu sein. Gibt es eine offizielle Definition?

Heute will ich auf die Parameter des Modells eingehen und wie sie bestimmt wurden.

Verfahren

Die meisten Parameter werden durch Anpassung an die Fallzahlen des RKI bestimmt. Aus den täglich erneuerten Tabellen--Files Lassen sich die gemeldeten Fallzahlen ablesen. Deren kumulative Aufsummierung ergibt eine mit der Zeit monoton steigende Reihe von Zahlen. Um diese mit dem Modell vergleichen zu können, berechne ich eine zusätzliche Funktion K, die kumulative Zahl der erkannten Kranken. Ihre Ableitung ist der positive Term in der Ableitung von Q

dK/dt = I₁/γ₁.

das Modell wird mit verschiedenen Parametersätzen berechnet und jeweils SSR, sum of squared residuals , gebildet, durch Berechnung der tägllchen Differenz zwischen K und den Daten des RKI. SSR ist eine Bewertung des Parametersatzes. Man versucht SSR zu minimieren und durch Wahl eines geeigneten Parametersatzes "das Modell an die Daten zu fitten".

Ich habe mir durch Variation der Parameter nach subjektiver Überlegung einen Eindruck verschafft von der Wirkung der Parameter auf die Kurve K(t) und ihre Abweichung von den Daten. Wenn ich eine halbwegs passable Anpassung erreicht hatte, habe ich mathematische Methoden (Least Squares ) zum lokalen Feintuning benutzt.

q

q bedeutet den Anteil (0 ≤ q ≤ 1) an Erkrankten, die vom Gesundheitssystem isoliert werden. Mit der Parameteranpassung hatte ich hier keinen Erfolg. Es schien, dass der Parameter geringen Einfluss hatte und deshalb aus diesen Fallzahlen nicht abgelesen werden kann.

Ich habe angenommen dass fast alle Fälle symptomatisch werden und dass Deutschland effektiv im Testen solcher Patienten ist. Außerdem werden Fälle erkannt bevor sie Symptome zeigen, weil Kontakte zurückverfolgt werden. Das entspricht, meine ich, einem

q = 0.8

R₀

soll die Zahl an Ansteckungen sein, die von einem Infizierten ausgeht, der in eine Population von ansteckbaren kommt. Der Parameter im Modell hat diese Eigenschaft: in der Gleichung für dS/dt kann man für diese Überlegung die beiden Klassen der Infektösen zusammenfassen, I = I₁ + I₂. Diese Zahl repräsentiert in der Gleichung die Zahl derer von denen Infektionen ausgehen. Hier ist I = 1 zu setzen, denn es geht nur um eine Person. Ferner ist in der beschriebenen Situation S = N. Sei γ = γ₁ + γ₂ die Dauer der Infektiösität, dann bleibt

dS/dt = -R₀/γ,

eine Konstante. Wenn wir die von 0 bis γ integrieren kommt -R₀ heraus . Das Negative davon ist die Zahl der hier neu infizierten.

R₀ hat also eine reale Bedeutung. Man kann es z.B in isolierten Heimen beobachten, wo man Ansteckungsketten verfolgen kann. Es ist klar, dass R₀ nicht nur eine feststehende Eigenschaft des Erregers ist, sondern von Art und Häufigkeit von Kontakten abhängt.

Es liegt deshalb nahe, R₀ als zetlich veränderlich anzusehen. Die Abbildung zeigt anhand der Daten, dass ab Tagg 20 eine Veränderung im Sinne einer Abflachung eintritt. Ich habe für das Modell die Zeit in Intervalle zerlegt. Das erste umfasst die ersten 20 Tage und es folgen weitere, jeweils 10 Tage breit. In jedem Intervall kann anderer R₀--Wert benutzt werden. Geeignete Werte können durch automatische Optimierung gefunden werden, wenn die Länge der Datenreihen ausreicht.

I₀

Dieser Parameter hat schon einen Einfluss auf die Kurve K(t), das aber abhängig vom Parameter R₀ . Anscheinend konnte man über den einen wie den anderen eine Anpassung erreichen. Ich habe deshalb das Intervall 2.5 ≤ I₀ ≤ 4.5 in 50 Teilintervalle zerlegt und für jedes die drei R₀--Werte durch die Optimierungssoftware bestimmen lassen. Ich benutzte die Daten vom 17.4.20 und teilte die Zeit in drei Intervalle, t<20, 20≥t<30 und t≥30.

Man sieht hier für vorgegebene Werte von I₀ das Ergebnis der Optimierung der drei R₀--Werte für die drei Zeitintervalle. SSR zeigt ein deutliches Minimum bei

I₀ =3000.

Diesen Wert halte ich fest.

Mein Modell orientiert sich an dem von
an der Heiden und Buchholz, das aber nicht in allen technischen Einzelheiten verfügbar gemacht wurde. Mein Modell ist jedenfalls einfacher, weil ich die Belastung des Gesundheitssystems nicht betrachte.

Die Grundstruktur des Modells ist eine Einteilung der deutschen
Bevölkerung in Klassen folgender Art



Die Klassen sind definiert als

S

Die ansteckbaren (susceptable)

E

Die frisch infizierten, noch nicht ansteckenden (exposed)

I₁

Die ansteckenden (infectious) der ersten Phase

I₂

Die ansteckenden (infectious) der zweiten Phase

Q

Die erkannten Erkrankten, die in Quarantäne sind

R

die die nach einer Erkrankung immun oder tot sind (recovered)

Für die einzelnen PatientInnen verläuft die Krankheit im Mittel so: nach der Ansteckung sind sie zunächst selbst noch nicht ansteckend (E, 3Tage). Dann beginnt bereits ihre Infektiösität, aber es dauert noch 2 Tage, bis sich Symptome zeigen. Danach dauert es im Mittel 4 Tage bis sie getestet und isoliert wird (I₁, 6Tage). Diese Erfassung und Isolierung gelingt nur in einem Teil der Fälle (z.B. 80%). Die restlichen bleiben infektiös (I₂, 4Tage). In Quarantäne bleiben die Patenten im Schnitt 8Tage.Diese Zeitangaben habe ich aus an der Heiden und Buchholz entnommen.

Im Modell werden die Anzahlen an Personen in den einzelnen Klassen als (reelle) Zahlen betrachtet, die als Zustandsvariablen bezeichnet werden. Sie werden ebenso bezeichnet wie die Klassen. Ihre Summe ist die Zahl der Einwohner

S + E + I₁ + I₂ + Q + R = N = 83 Millionen

und wird hier als konstant angenommen. Die Umverteilung zwischen den Klassen wird durch das folgende Gleichungssystem beschrieben:

Hier sind verschiedene Parameter unbekannt. R₀⁢, q und die verschiedenen γ. Die letzteren sind aber die mittleren Verweildauern in Tagen, die ich gerade erwähnt habe.

Die unabhängige Varable t ist die Zeit in Tagen. Deshalb sind die Verweildauern γ in Tagen angegeben.

Das Modell zu berechnen, bedeutet, ab einer Zeit t₀⁢ die Werte der Zustandsvariablen für alle Zukunft so anzugeben, dass das Differentialgleichungssystem erfüllt ist. Voraussetzung ist. dass die Startwerte, das sind die Werte der Zustandsvariablen bei t₀⁢, bekannt sind. Mein Startpunkt soll der 1.3.2020 sein. Ich nehme an dass es diesem Tag eine gewisse Zahl I₀⁢ von Infizierten existiert, die ich, mangels Kenntnis, zu gleichen Teilen auf E, I₁ und I₂ verteile. Die Startwerte sind also

S = N - I₀⁢
E = I₀/3
I₁ = I₀/3
I₂ = I₀/3
⁢Q = 0
R = 0


Nun ist, bis auf drei freie Paramter R₀, q und I₀ das Modell definiert. Zur Anpassung der Parameter an die Daten des RKI werde ich in Kürze berichten.