SEIR-Modell

Epidemie-Verlauf im SEIR-Modell für eine Latenzzeit von 1 bis 4 Tagen. Im Vergleich dazu gestrichelt der Verlauf im SIR-Modell. Die Abweichung vom SIR-Modell vergrößert sich mit steigender Latenzzeit. Die Kurven zeigen die Neuinfektionen pro Tag in Abhängigkeit der Zeit in Tagen, relativ zur Populationsgröße. Parameter: Basisreproduktionszahl 2,0, infektiöse Zeit von 4 Tagen, 100 Infektiöse pro Mio. zum Startzeitpunkt.

Als SEIR-Modell bezeichnet man in der mathematischen Epidemiologie einen Ansatz zur Beschreibung der Ausbreitung von ansteckenden Krankheiten. Die Beschreibung ist näher am realen Verlauf als die des SIR-Modells, da hier berücksichtigt wird, dass ein Individuum nach seiner Infektion nicht sofort selbst infektiös ist. Im Gegensatz zu einem Individuum-basierten Modell ist die Beschreibung makroskopisch, d. h. die Population wird als Gesamtheit betrachtet.

Klassisches SEIR-Modell

Modelldefinition

Die Population von N Individuen sei zerlegt in die vier Kompartimente S, E, I, R, so dass

S + E + I + R = N . {\displaystyle S+E+I+R=N.}

Jedes Individuum kann die Prozedur

Susceptible (S) → Exposed (E) → Infectious (I) → Recovered (R)

durchlaufen. Die Ausbreitungsdynamik der Krankheit wird beschrieben durch das System

S = 1 N β S I , E = 1 N β S I α E , I = α E γ I , R = γ I , {\displaystyle {\begin{aligned}S'&=-{\tfrac {1}{N}}{\beta }SI,\\E'&={\tfrac {1}{N}}\beta SI-\alpha E,\\I'&=\alpha E-\gamma I,\\R'&=\gamma I,\end{aligned}}}

wobei

S ( t ) = d S ( t ) d t , E ( t ) = d E ( t ) d t , I ( t ) = d I ( t ) d t , R ( t ) = d R ( t ) d t {\displaystyle S'(t)={\frac {\mathrm {d} S(t)}{\mathrm {d} t}},\quad E'(t)={\frac {\mathrm {d} E(t)}{\mathrm {d} t}},\quad I'(t)={\frac {\mathrm {d} I(t)}{\mathrm {d} t}},\quad R'(t)={\frac {\mathrm {d} R(t)}{\mathrm {d} t}}}

die Ableitungen nach der Zeit sind. Betrachtet man die bezüglich der konstanten Population relativen Anteile

s := S N , e := E N , i := I N , r := R N , {\displaystyle s:={\frac {S}{N}},\quad e:={\frac {E}{N}},\quad i:={\frac {I}{N}},\quad r:={\frac {R}{N}},}

nimmt das System die Gestalt

s = β s i , e = β s i α e , i = α e γ i , r = γ i {\displaystyle {\begin{aligned}s'&=-\beta si,\\e'&=\beta si-\alpha e,\\i'&=\alpha e-\gamma i,\\r'&=\gamma i\end{aligned}}}

an. Hierbei handelt es sich um ein autonomes nichtlineares System von gewöhnlichen Differentialgleichungen, weshalb das Modell ein dynamisches System beschreibt. Jeder Zustand ohne Exponierte und Infektiöse stellt eine Ruhelage des dynamischen Systems dar, die krankheitsfreies Gleichgewicht genannt wird, zuweilen abgekürzt als DFE für disease free equilibrium. Sollte die effektive Reproduktionszahl den kritischen Wert 1 übersteigen, handelt es sich um eine instabile Ruhelage, aus deren Umfeld sich eine Epidemie entwickelt, woraufhin mit dem Durchlaufen der Epidemie eines der stabilen krankheitsfreien Gleichgewichte angestrebt wird.

Größe Einheit Erklärung
S(t) indv Anteil der Anfälligen, engl. susceptible. Noch nicht infiziert.
E(t) indv Anteil der Exponierten, engl. exposed. Infiziert, aber noch nicht infektiös.
I(t) indv Anteil der Infektiösen, engl. infectious.
R(t) indv Anteil der Erholten, engl. recovered oder resistant. Bzw. verstorben oder nach Symptomen in Quarantäne.
t d Zeit in Tagen, engl. time.
β 1/d Transmissionsrate. Der Kehrwert ist die mittlere Zeit zwischen Kontakten mit Übertragung.
γ 1/d Erholungsrate. Der Kehrwert ist die mittlere infektiöse Zeit.
α 1/d Übergangsrate. Der Kehrwert ist die mittlere Latenzzeit.

Die mittlere Latenzzeit ist die durchschnittliche Zeit, die ein Individuum in der Gruppe E der Exponierten verbringt; diese ist zu unterscheiden von der mittleren Inkubationszeit, denn der Beginn der Infektiosität muss nicht mit dem Beginn der Symptome übereinstimmen.

Für die Transmissionsrate (Übertragungsrate) ist auch die Bezeichnung Kontaktrate geläufig. Eine genauere Überlegung zerlegt diese in ein Produkt β = k b {\displaystyle \beta =kb} , wobei b {\displaystyle b} die Transmissionswahrscheinlichkeit und k {\displaystyle k} die eigentliche Kontaktrate ist.

Modellannahmen

Das Modell macht eine Reihe von Annahmen, welche die Wirklichkeit stark vereinfachen. Im Einzelnen:

  • Es findet eine gute Durchmischung der Bevölkerung statt, so dass jedes mit jedem anderen Individuum in Kontakt kommen kann. Andernfalls müsste man Kontaktnetzwerke oder räumliche Dynamiken beschreiben.
  • Es gibt hinreichend viele Infektiöse, so dass der Verlauf der Epidemie näherungsweise deterministisch stattfindet. Andernfalls müsste man einen stochastischen Prozess beschreiben.
  • Die Werte sind kontinuierlich. Das Modell beschreibt keine diskreten Anzahlen.
  • Sowohl die Dauer der latenten Phase als auch der infektiösen Phase ist exponentialverteilt (siehe den Abschnitt Infektionsalter-Modell unten). Die Transmissionsrate in der infektiösen Phase ist hierbei konstant.
  • Die Krankheit führt entweder zur vollständigen und lebenslangen Immunität oder zum Versterben.
  • Jedes anfällige Individuum ist gleich anfällig für die Krankheit.
  • Ein Nachwuchs neuer Anfälliger und natürliches Versterben tritt nicht auf. Das Modell enthält keine Beschreibung der demografischen Dynamik. Insbesondere ist die Größe der Bevölkerung konstant.

Beziehung zur Basisreproduktionszahl

Zwischen den Parametern des Modells und der Basisreproduktionszahl lässt sich eine Beziehung herstellen. Damit sich die Krankheit nicht weiter ausbreitet, muss e = 0 {\displaystyle e'=0} und i = 0 {\displaystyle i'=0} sein. Einsetzung dieser Bedingungen in die Differentialgleichungen führt zu α e = β s i {\displaystyle \alpha e=\beta si} und α e = γ i {\displaystyle \alpha e=\gamma i} , und somit β s = γ {\displaystyle \beta s=\gamma } . Diese Gleichgewichtsbedingungen entsprechen der effektiven Reproduktionszahl R q = 1 {\displaystyle R_{q}=1} . Weil aber per Definition R q := R 0 s {\displaystyle R_{q}:=R_{0}s} gilt, ergibt sich[1]

R 0 = β γ . {\displaystyle R_{0}={\frac {\beta }{\gamma }}.}

Gelegentlich wird β {\displaystyle \beta } alternativ so definiert, dass β N = R 0 γ {\displaystyle \beta N=R_{0}\gamma } gilt, mit der Bewandtnis, bei Benutzung der Großbuchstaben den Parameter N {\displaystyle N} nicht in den Gleichungen mitführen zu müssen.[2][3] Unterscheidungen dieser Art sind durch Normalisierung auf N = 1 {\displaystyle N=1} entfernbar, denn das Modell kodiert per se keine diskreten Zahlen.

Man definiert q := 1 s {\displaystyle q:=1-s} , das ist der immune Anteil der Population. Im Gleichgewicht

1 = R q = R 0 ( 1 q ) {\displaystyle 1=R_{q}=R_{0}\cdot (1-q)}

ist demnach

q = 1 1 R 0 . {\displaystyle q=1-{\frac {1}{R_{0}}}.}

Dieser Anteil, oberhalb dessen Herdenimmunität besteht, heißt kritische Immunisierungsschwelle. Diese Schwelle ist nur von der Basisreproduktionszahl abhängig, nicht aber vom gewählten Modell.

Der Parameter α {\displaystyle \alpha } , der den Übergang in das Zwischenstadium beschreibt, in der die Person infiziert, aber noch nicht ansteckend ist, spielt bei der Basisreproduktionszahl keine Rolle, solange man die Möglichkeit des Versterbens in der Gruppe der Exponierten vernachlässigt.[2] Falls man dies nicht tut, erhält man eine Formel für R 0 {\displaystyle R_{0}} wie unten im Abschnitt Einbeziehung demografischer Dynamik angegeben.

Die Beziehung der Parameter zur Basisreproduktionszahl ist also im Allgemeinen für jede Modifikation des Modells neu zu ermitteln. Eine systematische Methode der Bestimmung wurde im Jahr 1989 von Odo Diekmann, Hans Heesterbeek und Hans Metz gefunden,[4] und später von Pauline van den Driessche und James Watmough weiter ausgearbeitet.[5] Bei dieser Herangehensweise wird die Jacobi-Matrix des Systems am krankheitsfreien Gleichgewicht dergestalt zu einer Differenz zerlegt, dass der Minuend die Neuinfektionen charakterisiert und der Subtrahend die sonstigen Übergänge zwischen den Kompartimenten. Das Produkt aus dem Minuend und der Inverse des Subtrahenden bildet eine neue Matrix, von ihnen Next generation matrix genannt. Ihr Spektralradius ist die Basisreproduktionszahl.

Maximal möglicher Anteil an Infizierten

Die Basisreproduktionszahl bestimmt wesentlich den maximal möglichen Anteil von Infizierten an einer gegebenen Population. Bezeichnen wir den Gesamtanteil aller Infizierten mit j {\displaystyle j} , also j ( t ) := e ( t ) + i ( t ) {\displaystyle j(t):=e(t)+i(t)} , so folgt aus dem obigen Differentialgleichungssystem

d j d s = d e + d i d s = γ i + β s i β s i = γ β s 1 {\displaystyle {\frac {\mathrm {d} j}{\mathrm {d} s}}={\frac {\mathrm {d} e+\mathrm {d} i}{\mathrm {d} s}}={\frac {-\gamma i+\beta si}{-\beta si}}={\frac {\gamma }{\beta s}}-1} ,

oder mit der Basisreproduktionszahl R 0 = β γ {\displaystyle R_{0}={\frac {\beta }{\gamma }}} einfach

d j d s = 1 R 0 s 1. {\displaystyle {\frac {\mathrm {d} j}{\mathrm {d} s}}={\frac {1}{R_{0}s}}-1.}

Mit Umstellung dieser Gleichung nach d j = d s R 0 s d s {\displaystyle \mathrm {d} j={\frac {\mathrm {d} s}{R_{0}s}}-\mathrm {d} s} liefert eine Integration durch Trennung der Variablen

j ( t ) + s ( t ) 1 R 0 ln s ( t ) = const {\displaystyle j(t)+s(t)-{\frac {1}{R_{0}}}\ln s(t)={\text{const}}}

für alle t, wobei ln {\displaystyle \ln } der natürliche Logarithmus ist. Insbesondere gilt damit

j ( t ) + s ( t ) 1 R 0 ln s ( t ) = j ( 0 ) + s ( 0 ) 1 R 0 ln s ( 0 ) {\displaystyle j(t)+s(t)-{\frac {1}{R_{0}}}\ln s(t)=j(0)+s(0)-{\frac {1}{R_{0}}}\ln s(0)}

oder äquivalent

j ( t ) = j ( 0 ) + s ( 0 ) 1 R 0 ln s ( 0 ) + 1 R 0 ln s ( t ) s ( t ) =:   f ( s ( t ) ) . {\displaystyle j(t)=j(0)+s(0)-{\frac {1}{R_{0}}}\ln s(0)+\underbrace {{\frac {1}{R_{0}}}\ln s(t)-s(t)} _{=:\ f(s(t))}.}

Die Hilfsfunktion f ( x ) = ln x R 0 x {\displaystyle f(x)=\ln {\frac {x}{R_{0}}}-x} hat wegen f ( x ) = 1 R 0 x 1 {\displaystyle f'(x)={\frac {1}{R_{0}x}}-1} ein Maximum bei x = 1 R 0 {\displaystyle x={\frac {1}{R_{0}}}} . Der maximal mögliche Anteil j max {\displaystyle j_{\max }} der Infizierten wird also für s ( t ) = 1 R 0 {\displaystyle s(t)={\frac {1}{R_{0}}}} erreicht, d. h. j max {\displaystyle j_{\max }} hängt nur von der Basisreproduktionszahl und den Anfangswerten von j {\displaystyle j} und s {\displaystyle s} ab:

j max = j ( 0 ) + s ( 0 ) 1 R 0 ln s ( 0 ) 1 R 0 + 1 R 0 ln 1 R 0 = j ( 0 ) + s ( 0 ) + 1 R 0 ( ln 1 R 0 s ( 0 ) 1 ) . {\displaystyle j_{\text{max}}=j(0)+s(0)-{\frac {1}{R_{0}}}\ln s(0)-{\frac {1}{R_{0}}}+{\frac {1}{R_{0}}}\ln {\frac {1}{R_{0}}}=j(0)+s(0)+{\frac {1}{R_{0}}}\left(\ln {\frac {1}{R_{0}s(0)}}-1\right).}

Für eine neu auftretende Erkrankung wie eine Epidemie durch ein unbekanntes Virus gilt s ( 0 ) = 1 {\displaystyle s(0)=1} und j ( 0 ) = 0 {\displaystyle j(0)=0} , d. h. der maximal mögliche Anteil der Infizierten an der Population hängt in diesem Fall wie folgt von der Basisreproduktionszahl ab:

j max = 1 + 1 R 0 ( ln 1 R 0 1 ) . {\displaystyle j_{\text{max}}=1+{\frac {1}{R_{0}}}\left(\ln {\frac {1}{R_{0}}}-1\right).}

Diese Gleichung entspricht der unter SIR-Modell #Maximale Zahl der Infizierten angegebenen Gleichung für das SIR-Modell.

Ausmaß der Epidemie nach ihrem Abklingen

Im Zusammenhang mit der Basisreproduktionszahl steht auch, welcher Anteil der Population insgesamt infiziert wird, unter Annahme, die Epidemie würde ohne jegliche Quarantäne durchlaufen. Unter Heranziehung der Differentialgleichungen findet man

d s d r = s ( t ) r ( t ) = β s γ = R 0 s . {\displaystyle {\frac {\mathrm {d} s}{\mathrm {d} r}}={\frac {s'(t)}{r'(t)}}=-{\frac {\beta s}{\gamma }}=-R_{0}s.}

Für den Anfangswert r ( 0 ) = 0 {\displaystyle r(0)=0} ist demnach

s ( t ) = s ( 0 ) e R 0 r ( t ) . {\displaystyle s(t)=s(0){\mathrm {e} }^{-R_{0}r(t)}.}

Bei t {\displaystyle t\to \infty } ist nun e = i = 0 {\displaystyle e=i=0} , und daher s + r = 1 {\displaystyle s+r=1} . Daraus resultiert die Gleichung

1 q = s ( 0 ) e R 0 q , q := lim t r ( t ) . {\displaystyle 1-q=s(0){\mathrm {e} }^{-R_{0}q},\quad q:=\lim _{t\to \infty }r(t).}

Dieser Zusammenhang, Final size equation genannt, besitzt eine über das SEIR-Modell hinausgehende Bedeutung. Algebraische Umformung führt zur Gleichung

R q e R q = ( q 1 ) R 0 e ( q 1 ) R 0 = s ( 0 ) R 0 e R 0 . {\displaystyle -R_{q}{\mathrm {e} }^{-R_{q}}=(q-1)R_{0}{\mathrm {e} }^{(q-1)R_{0}}=-s(0)R_{0}{\mathrm {e} }^{-R_{0}}.}

Mit x := s ( 0 ) R 0 e R 0 {\displaystyle x:=-s(0)R_{0}\mathrm {e} ^{-R_{0}}} und y := ( q 1 ) R 0 {\displaystyle y:=(q-1)R_{0}} lautet die letzte Gleichung y e y = x {\displaystyle y\mathrm {e} ^{y}=x} , die mit der lambertschen W-Funktion W ( x ) {\displaystyle W(x)} nach y {\displaystyle y} umgestellt werden kann, d. h. y = W ( x ) {\displaystyle y=W(x)} , was wieder zurückersetzt und nach q {\displaystyle q} umgestellt

q = 1 + 1 R 0 W ( s ( 0 ) R 0 e R 0 ) {\displaystyle q=1+{\frac {1}{R_{0}}}W(-s(0)R_{0}{\mathrm {e} }^{-R_{0}})}

ergibt.

Für die praktische Berechnung des hier relevanten Teils der W-Funktion betrachtet man das quadratische Taylorpolynom der Funktion f ( w ) := w e w x {\displaystyle f(w):=w\mathrm {e} ^{w}-x} an der Stelle −1 und bestimmt von diesem die Nullstelle, auf welche noch einmalig die Fixpunktiteration w x e w {\displaystyle w\mapsto x\mathrm {e} ^{-w}} angewendet wird. Das Resultat ist

w = x exp ( 1 2 e x + 2 ) . {\displaystyle w=x\exp(1-{\sqrt {2\mathrm {e} x+2}}).}

Eine gute Näherung der W-Funktion erhält man nun als

W ( x ) φ n ( w ) , {\displaystyle W(x)\approx \varphi ^{n}(w),}

wobei φ n {\displaystyle \varphi ^{n}} die n-te Iteration des Newtonverfahrens

φ ( w ) = w f ( w ) f ( w ) = x e w + w 2 w + 1 {\displaystyle \varphi (w)=w-{\frac {f(w)}{f'(w)}}={\frac {x\mathrm {e} ^{-w}+w^{2}}{w+1}}}

ist. Für alle praktischen Bedürfnisse ist n = 4 {\displaystyle n=4} völlig ausreichend.

Exponentielle Anfangsphase

Am Anfang der Epidemie verläuft die Ausbreitung der Krankheit in guter Näherung exponentiell. Eine Untersuchung dieser anfänglichen Phase gelingt mit der folgenden systematischen Methode, die auch auf andere Modelle anwendbar ist. Der Vektor x = ( s , e , i , r ) {\displaystyle x=(s,e,i,r)} fasst die vier Anteile zusammen, so dass das autonome System in der abstrakten Form x = f ( x ) {\displaystyle x'=f(x)} dargestellt werden kann. An einem Zustand x 0 {\displaystyle x_{0}} nähert die Linearisierung

f ( x ) f ( x 0 ) + J ( x x 0 ) {\displaystyle f(x)\approx f(x_{0})+J(x-x_{0})}

das System an, wobei mit J := D f ( x 0 ) {\displaystyle J:=Df(x_{0})} die Jacobi-Matrix gemeint ist. Der Anfang der Epidemie liegt am krankheitsfreien Gleichgewicht x 0 = ( s 0 , 0 , 0 , 1 s 0 ) {\displaystyle x_{0}=(s_{0},0,0,1-s_{0})} , das als Ruhelage die Gleichung f ( x 0 ) = 0 {\displaystyle f(x_{0})=0} erfüllt, also eine Nullstelle des Vektorfeldes f {\displaystyle f} ist. Die Jacobi-Matrix besitzt in dieser Ruhelage den Wert

J = ( 0 0 β s 0 0 0 α β s 0 0 0 α γ 0 0 0 γ 0 ) . {\displaystyle J={\begin{pmatrix}0&0&-\beta s_{0}&0\\0&-\alpha &\beta s_{0}&0\\0&\alpha &-\gamma &0\\0&0&\gamma &0\end{pmatrix}}.}

Die vierte Zeile des Systems darf entfallen, da sie von den übrigen entkoppelt vorliegt, womit sich die Jacobi-Matrix auf drei Zeilen und Spalten verkleinert. Im linearisierten System ist zudem die erste Zeile entkoppelt, womit eigentlich die Betrachtung des kleinen Systems

( e i ) = ( α β s 0 α γ ) ( e i ) {\displaystyle {\begin{pmatrix}e'\\i'\end{pmatrix}}={\begin{pmatrix}-\alpha &\beta s_{0}\\\alpha &-\gamma \end{pmatrix}}{\begin{pmatrix}e\\i\end{pmatrix}}}

genügt. Die allgemeine Lösung des linearisierten Systems erhält man mit Standardmethoden. Man kann nachrechnen, dass die Jacobi-Matrix für R 0 s 0 1 {\displaystyle R_{0}s_{0}\neq 1} immer diagonalisierbar ist. Die Lösung ist demzufolge von der Form x ( t ) = x 0 + k c k e λ k t v k {\displaystyle \textstyle x(t)=x_{0}+\sum _{k}c_{k}\mathrm {e} ^{\lambda _{k}t}v_{k}} , wobei v k {\displaystyle v_{k}} der Eigenvektor zum Eigenwert λ k {\displaystyle \lambda _{k}} ist und die Konstanten c k {\displaystyle c_{k}} aus der Anfangsbedingung zu ermitteln sind.

Diese allgemeine Lösung gilt auch für sonderbare Anfangswerte. Einer der Summanden fällt jedoch rasch auf null ab, weshalb man ein wenig später für die Neuinfektionen eigentlich nur eine exponentielle Kurve c ( t ) = c 0 e λ t {\displaystyle c(t)=c_{0}\mathrm {e} ^{\lambda t}} beobachtet, womit sich die Lösungen ebenfalls zu Exponentialfunktionen vereinfachen. Das Verhältnis von Exponierten zu Infizierten strebt mit dem Abfallen des Summanden einen speziellen Wert an. Weil infolge i = λ i {\displaystyle i'=\lambda i} aufgrund c = s = β s 0 i {\displaystyle c=-s'=\beta s_{0}i} gilt, finden sich die natürlichen Anfangswerte

s ( 0 ) = s 0 e ( 0 ) i ( 0 ) , e ( 0 ) = 1 α ( λ + γ ) i ( 0 ) , i ( 0 ) = c 0 β s 0 . {\displaystyle s(0)=s_{0}-e(0)-i(0),\quad e(0)={\frac {1}{\alpha }}(\lambda +\gamma )i(0),\quad i(0)={\frac {c_{0}}{\beta s_{0}}}.}

Weiterhin ergibt sich nun, dass die Wachstumskonstante λ {\displaystyle \lambda } als Eigenwert der Jacobi-Matrix die quadratische Gleichung

α β s 0 = ( λ + α ) ( λ + γ ) {\displaystyle \alpha \beta s_{0}=(\lambda +\alpha )(\lambda +\gamma )}

erfüllen muss, womit eine Beziehung zwischen der Wachstumskonstante und den Parametern des Modells gefunden ist.[6] Wie bei jedem exponentiellen Wachstumsvorgang steht die Wachstumskonstante in eins-zu-eins-Beziehung zur anschaulicheren Vervielfachungszeit

T n = ln ( n ) λ , {\displaystyle T_{n}={\frac {\ln(n)}{\lambda }},}

speziell zur Verdopplungszeit T 2 {\displaystyle T_{2}} .

Die Linearisierung an der Ruhelage ist überdies für die Untersuchung der Stabilität von Bedeutung. Besitzt die Jacobi-Matrix keinen Eigenwert mit Realteil null, spricht man von einer hyperbolischen Ruhelage. Laut dem Linearisierungssatz, einem Grundresultat der Stabilitätstheorie dynamischer Systeme, besitzt das linearisierte System in der Nähe der hyperbolischen Ruhelage dasselbe qualitative Verhalten wie das ursprüngliche. Somit ist eine hyperbolische Ruhelage genau dann stabil, wenn sich sämtliche Eigenwerte links der imaginären Achse befinden. Weil die vorliegende Jacobi-Matrix allerdings den Eigenwert null enthält, müsste man weitergehende Mittel zur Untersuchung der Stabilität einsetzen. Eine Umgehung dieses Problems bietet das später beschriebene Modell mit demografischer Dynamik im speziellen Fall, wo die Geburtenrate mit der Sterberate übereinstimmt. Dann ist N konstant bezüglich der Zeit, womit man das äquivalente System definieren darf, in dem N eine explizite Konstante ist. Der Eigenwert null tritt nun nicht mehr auf. Damit einhergehend verbleibt aber aufgrund des Versterbens der Resistenten nur noch die vollständig anfällige Population als einziges krankheitsfreies Gleichgewicht. Ferner tritt im instabilen Fall das stabile endemische Gleichgewicht in Erscheinung, das im einfachen Modell nicht existiert.

Die Untersuchung zur Stabilität wurde von van den Driessche in allgemeiner Weise durchgeführt. Sie kam zum Resultat, dass das krankheitsfreie Gleichgewicht bei einem deterministischen Kompartiment-Modell unter gewissen natürlichen Voraussetzungen an die Modelldefinition erwartungsgemäß immer dann instabil ist, wenn die Basisreproduktionszahl über dem kritischen Wert 1 liegt, und zumindest lokal asymptotisch stabil, wenn sie darunter liegt. Ebenso findet man bei Erläuterungen dieser Untersuchung aufgrund der besagten Problematik zunächst die Beschränkung auf solche Modelle vor, wo ein einziges krankheitsfreies Gleichgewicht existiert.

Beispielrechnung

Es folgt eine Beispielrechnung zu einer Parameterbelegung, wie sie im März 2020 für die COVID-19-Pandemie in Deutschland abgeschätzt wurde.[7] Die Ausbreitung verläuft bezüglich einer Basisreproduktionszahl von 2,4, was der Unterlassung wesentlicher Quarantäne entspricht. Zur numerischen Lösung des Anfangswertproblems genügt das Euler-Verfahren.

from numpy import array as vector

# Explizites Euler-Verfahren
def euler_method(f, t0, x0, t1, h):
    t = t0; x = x0
    a = [[t, x]]
    for k in range(0, 1 + int((t1 - t0)/h)):
        t = t0 + k*h
        x = x + h*f(t, x)
        a.append([t, x])
    return a

def SEIR_model(alpha, beta, gamma):
    def f(t, x):
        s, e, i, r = x
        return vector([
            -beta*s*i,
            beta*s*i - alpha*e,
            alpha*e - gamma*i,
            gamma*i
        ])
    return f

def SEIR_simulation(alpha, beta, gamma, e0, i0, days, step=0.1):
    x0 = vector([1.0 - e0 - i0, e0, i0, 0.0])
    f = SEIR_model(alpha, beta, gamma)
    return euler_method(f, 0, x0, days, step)

def diagram(simulation):
    import matplotlib.pyplot as plot
    plot.style.use('fivethirtyeight')
    figure,axes = plot.subplots()
    figure.subplots_adjust(bottom = 0.15)
    axes.grid(linestyle = ':', linewidth = 2.0, color = "#808080")
    t,x = zip(*simulation())
    s, e, i, r = zip(*x)
    axes.plot(t, s, color = "#0000cc")
    axes.plot(t, e, color = "#ffb000", linestyle = '--')
    axes.plot(t, i, color = "#a00060")
    axes.plot(t, r, color = "#008000", linestyle = '--')
    plot.show()

def simulation1():
    N = 83200000 # Einwohnerzahl von Deutschland 2019/2020
    R0 = 2.4; gamma = 1/3.0
    return SEIR_simulation(
        alpha = 1/5.5, beta = R0*gamma, gamma = gamma,
        e0 = 40000.0/N, i0 = 10000.0/N, days = 140)

diagram(simulation1)
Dynamics of an infectious disease, according to the basic SEIR model

Die vier Anteile, jeweils abhängig von der Zeit in Tagen.
s in Blau, e in Gelb gestrichelt, i in Magenta, r in Grün gestrichelt.

Ersichtlich ist an diesem Beispiel, dass die Epidemie aufgrund der Aufheizung durch die Infektiösen auch noch nach dem Erreichen der kritischen Immunisierungsschwelle

q c = 1 1 R 0 = 58 % {\displaystyle q_{\text{c}}=1-{\frac {1}{R_{0}}}=58\,\%}

weiterläuft. Insgesamt würden sich 88 % der Bevölkerung mit der Krankheit anstecken. Die Epidemie ließe sich allerdings spätestens nach Erreichen der kritischen Immunisierungsschwelle durch eine ca. einmonatige strenge Quarantäne stoppen.

Modellerweiterungen

Einbeziehung demografischer Dynamik

Flussdiagramm

Unter der Annahme einer konstanten Geburtenrate ν {\displaystyle \nu } und konstanten Sterberate μ {\displaystyle \mu } wurde das erweiterte Modell

S = ν N 1 N β S I μ S , E = 1 N β S I α E μ E , I = α E γ I μ I , R = γ I μ R {\displaystyle {\begin{aligned}S'&=\nu N-{\tfrac {1}{N}}\beta SI-\mu S,\\E'&={\tfrac {1}{N}}\beta SI-\alpha E-\mu E,\\I'&=\alpha E-\gamma I-\mu I,\\R'&=\gamma I-\mu R\end{aligned}}}

aufgestellt. Gegenüber dem einfachen SEIR-Modell beschreibt dieses Modell auch einen langfristigen endemischen Verlauf, bei dem es zu einer Schwingung des Anteils der Anfälligen kommen kann, bis sich dieser ins endemische Gleichgewicht eingependelt hat. Für verschwindende Latenzzeit geht das Modell ins korrespondierende SIR-Modell über, wo dieses Phänomen bereits auftritt.

Sofern sich die Geburten- und Sterberate voneinander unterscheiden, ist die Populationsgröße

N ( t ) = S ( t ) + E ( t ) + I ( t ) + R ( t ) {\displaystyle N(t)=S(t)+E(t)+I(t)+R(t)}

keine Konstante mehr, sondern verläuft gemäß

N = ( ν μ ) N . {\displaystyle N'=(\nu -\mu )N.}

Die Beziehung zur Basisreproduktionszahl ist

R 0 = α β ( μ + α ) ( μ + γ ) . {\displaystyle R_{0}={\frac {\alpha \beta }{(\mu +\alpha )(\mu +\gamma )}}.}

Ist das Wachstum der Bevölkerung oder ihr Schwund über den betrachteten Zeitraum verschwindend gering, kann man die Geburten- und Sterberate zur Vereinfachung gleich setzen. In diesem Fall pendelt sich der Anteil der Anfälligen genau in den Gegenwert der kritischen Immunisierungsschwelle ein.

Zeitabhängige Transmissionsrate

Vorgänge wie Verhaltensänderungen, Quarantäne und Saisonalität bewirken eine Veränderung der Transmissionsrate. Diese Umstände finden ihre Berücksichtigung in der Modellierung der Transmissionsrate als zeitabhängige Funktion β ( t ) {\displaystyle \beta (t)} , wobei das übrige Modell identisch beibehalten wird.[8] Die einfachsten Ansätze für die Saisonalität nehmen die Transmissionsrate zum Beispiel als Sinusschwingung an, mit Berg in den kälteren und Tal in den wärmeren Monaten.[9]

Zu bemerken ist, dass das Differentialgleichungssystem mit der expliziten Zeitabhängigkeit kein autonomes mehr ist und damit nicht mehr direkt ein dynamisches System vorliegt. Man kann aus dem System allerdings durch Hinzunahme der Gleichung t = 1 {\displaystyle t'=1} künstlich ein autonomes gewinnen.

Interagierende Subpopulationen

Da soziale Kontakte, Schwere und Letalität der Erkrankung sowie die Wirksamkeit prophylaktischer Maßnahmen sich zwischen interagierenden Subpopulationen deutlich unterscheiden können, werden im „Interacting Subpopulation SEIR model“ individuelle SEIR Modelle für jede Subgruppe konstruiert und durch Kontakt-Interaktionen zwischen den Gruppen verknüpft[10]. Solche Modelle werden zum Beispiel zur Modellierung der Covid-19-Pandemie verwendet, um personalisierte, beschleunigte Impfstrategien zu entwickeln[11], welche eine Verkürzung der Pandemie und eine Reduktion von Fallzahlen und Todesfällen voraussagen.

Verallgemeinerungen

Infektionsalter-Modell

Die im SEIR-Modell vorausgesetzte Exponentialverteilung der Latenzzeit und der infektiösen Zeit mag recht künstlich erscheinen. Ein allgemeines, von dieser Voraussetzung befreites Modell beschreibt den zeitlichen Verlauf des Anteils der Anfälligen durch die Gleichung[12]

s ( t ) = s ( t ) 0 n ( a ) s ( t a ) d a . {\displaystyle s'(t)=s(t)\int _{0}^{\infty }n(a)s'(t-a)\,\mathrm {d} a.}

Hierbei ist n ( a ) {\displaystyle n(a)} die Infektionsrate im Infektionsalter a {\displaystyle a} . Mit dem Infektionsalter ist nicht das Alter eines Individuums gemeint, sondern die seit der Infektion vergangene Zeit. s ( t ) {\displaystyle s'(t)} ist die Zeitableitung von s ( t ) {\displaystyle s(t)} .

Das Modell ist wie folgt interpretierbar. Ursache für die Neuinzidenz s ( t ) {\displaystyle -s'(t)} sind in der Vergangenheit infizierte Infektiöse. Um a {\displaystyle a} in der Zeit zurück kam s ( t a ) {\displaystyle -s'(t-a)} hinzu, die sich nun im Infektionsalter a {\displaystyle a} befinden und daher zu n ( a ) s ( t a ) {\displaystyle -n(a)s'(t-a)} beitragen. Das Integral über die gesamte Vergangenheit ist schließlich noch mit s ( t ) {\displaystyle s(t)} zu multiplizieren, weil zum aktuellen Zeitpunkt nur noch der Anteil s ( t ) {\displaystyle s(t)} anfällig ist.

Per Definition von n ( a ) {\displaystyle n(a)} ergibt sich die Basisreproduktionszahl nun gemäß der Formel

R 0 = 0 n ( a ) d a . {\displaystyle R_{0}=\int _{0}^{\infty }n(a)\,\mathrm {d} a.}

Die normierte Funktion g ( a ) := n ( a ) / R 0 {\displaystyle g(a):=n(a)/R_{0}} lässt sich als Dichte einer Wahrscheinlichkeitsverteilung auffassen, die man Verteilung der Generationszeit, englisch generation interval distribution, nennt. Ihr Erwartungswert, die Generationszeit, ist der mittlere Abstand zwischen dem Zeitpunkt der Infektion des Spenders und dem Zeitpunkt der Infektion des Empfängers der Krankheit. Im Modell bleibt die Final size equation gültig, womit das Ausmaß der Epidemie allein von der Basisreproduktionszahl abhängt, nicht jedoch von der Dichte.

Zur Anfangsphase der Epidemie ergibt sich aus dem Modell für eine vollständig anfällige Population ( s ( t ) 1 {\displaystyle s(t)\approx 1} in guter Näherung) mit dem Ansatz s ( t ) = s ( 0 ) e λ t {\displaystyle s'(t)=s'(0)\mathrm {e} ^{\lambda t}} zudem die epidemiologische Euler-Lotka-Gleichung

1 = R 0 0 g ( a ) e λ a d a . {\displaystyle 1=R_{0}\int _{0}^{\infty }g(a)\mathrm {e} ^{-\lambda a}\,\mathrm {d} a.}

Diese besagt, bei Kenntnis der Dichte und der Wachstumskonstante λ {\displaystyle \lambda } ist die Basisreproduktionszahl ohne Kenntnis der Infektionsrate ermittelbar.

Um zu verstehen, in welchem Bezug die Modelle stehen, zerlegt man die Infektionsrate in das Produkt n ( a ) = β ( a ) p ( a ) {\displaystyle n(a)=\beta (a)p(a)} , wobei β ( a ) {\displaystyle \beta (a)} die Transmissionsrate ist, und p ( a ) {\displaystyle p(a)} die Wahrscheinlichkeit infektiös zu sein. Im Folgenden ist die Transmissionsrate konstant. Beginnt die infektiöse Zeit sofort und ist deren Dauer zufällig verteilt gemäß der Verteilungsfunktion F {\displaystyle F} , gilt p ( a ) = 1 F ( a ) {\displaystyle p(a)=1-F(a)} . Speziell für p ( a ) = e γ a {\displaystyle p(a)=\mathrm {e} ^{-\gamma a}} ergibt sich das SIR-Modell, wobei die mittlere infektiöse Dauer γ 1 {\displaystyle \gamma ^{-1}} nichts anderes als der Erwartungswert der Verteilung ist.

Dauert nun die Latenzzeit T 1 {\displaystyle T_{1}} zufällig lang gemäß Verteilungsfunktion F 1 {\displaystyle F_{1}} , und die danach beginnende infektiöse Zeit T 2 {\displaystyle T_{2}} zufällig lang gemäß Verteilungsfunktion F 2 {\displaystyle F_{2}} , erhält man

p ( a ) = P ( T 1 a < T 1 + T 2 ) = 0 a ( 1 F 2 ( a t ) ) d F 1 ( t ) . {\displaystyle p(a)=P(T_{1}\leq a<T_{1}+T_{2})=\int _{0}^{a}(1-F_{2}(a-t))\,\mathrm {d} F_{1}(t).}

Speziell für F 1 ( t ) = 1 e α t {\displaystyle F_{1}(t)=1-\mathrm {e} ^{-\alpha t}} und F 2 ( t ) = 1 e γ t {\displaystyle F_{2}(t)=1-\mathrm {e} ^{-\gamma t}} ergibt sich das SEIR-Modell. Die Erwartungswerte sind die mittlere Latenzzeit und die mittlere infektiöse Zeit.

Für einen nicht-zufälligen infektiösen Zeitraum mit Beginn im Infektionsalter a 1 {\displaystyle a_{1}} und Ende in a 2 {\displaystyle a_{2}} , also p ( a ) = χ [ a 1 , a 2 ] ( a ) {\displaystyle p(a)=\chi _{[a_{1},a_{2}]}(a)} , vereinfacht sich das Modell stattdessen zu

s ( t ) = β s ( t ) ( s ( t a 2 ) s ( t a 1 ) ) , {\displaystyle s'(t)=-\beta s(t)(s(t-a_{2})-s(t-a_{1})),}

eine retardierte Differentialgleichung. Die Beziehung zur Basisreproduktionszahl ist dann R 0 = β ( a 2 a 1 ) {\displaystyle R_{0}=\beta \cdot (a_{2}-a_{1})} .

Stochastische Varianten

CTMC-Modell

Waagerecht die Dichte der Verteilung des Zeitpunktes, an dem die Epidemie zu einer Wahl von Parametern im CTMC-Modell erlischt. Mit wachsender Basisreproduktionszahl trennt sich die Dichte in einen kleinen und einen großen Berg auf. Von einer Population von tausend Individuen sind zu Beginn vier infektiös.

Zur Beschreibung von Epidemie-Verläufen im Bereich weniger Infizierter wurde zu den klassischen Modellen jeweils eine direkte Entsprechung als zeitkontinuierlicher Markow-Prozess aufgestellt (CTMC: Continuous time Markov chain). In diesen stochastischen Modellen kann eine Epidemie vorläufig zum Erliegen kommen, obwohl die Reproduktionszahl größer als eins ist.

Für einen Markow-Prozess X {\displaystyle X} sei wie üblich

p a , b ( Δ t ) = P ( X ( t + Δ t ) = b X ( t ) = a ) {\displaystyle p_{a,b}(\Delta t)=P(X(t+\Delta t)=b\mid X(t)=a)}

die Übergangswahrscheinlichkeit vom Zustand a {\displaystyle a} in den Zustand b {\displaystyle b} . Für das SEIR-Modell ist X = ( S , E , I ) {\displaystyle X=(S,E,I)} und[13]

p ( s , e , i ) , ( x , y , z ) ( Δ t ) = o ( Δ t ) + { 1 N β s i Δ t , ( x , y , z ) = ( s 1 , e + 1 , i ) , α e Δ t , ( x , y , z ) = ( s , e 1 , i + 1 ) , γ i Δ t , ( x , y , z ) = ( s , e , i 1 ) , 1 ( 1 N β s i + α e + γ i ) Δ t , ( x , y , z ) = ( s , e , i ) , 0 , sonst. {\displaystyle p_{(s,e,i),(x,y,z)}(\Delta t)=o(\Delta t)+{\begin{cases}{\tfrac {1}{N}}\beta si\Delta t,&(x,y,z)=(s-1,e+1,i),\\\alpha e\Delta t,&(x,y,z)=(s,e-1,i+1),\\\gamma i\Delta t,&(x,y,z)=(s,e,i-1),\\1-({\tfrac {1}{N}}\beta si+\alpha e+\gamma i)\Delta t,&(x,y,z)=(s,e,i),\\0,&{\text{sonst.}}\end{cases}}}

Ein Verfahren zur Simulation dieser Markow-Prozesse ist der Gillespie-Algorithmus. Gemäß den Gesetzmäßigkeiten von Markow-Prozessen ist die Verweildauer im Zustand exponentialverteilt mit Verteilungsfunktion F ( t ) = 1 exp ( λ t ) {\displaystyle F(t)=1-\exp(-\lambda t)} , wobei die Ereignisrate λ {\displaystyle \lambda } die Summe der Einzelraten ist, das heißt

λ = 1 N β s i + α e + γ i . {\displaystyle \lambda ={\frac {1}{N}}\beta si+\alpha e+\gamma i.}

Mittels der Inversionsmethode ist für eine Standardzufallszahl u [ 0 , 1 ] {\displaystyle u\in [0,1]} die Verweildauer bzw. der Zeitschritt demnach

t = F 1 ( u ) = 1 λ ln ( 1 u ) . {\displaystyle t=F^{-1}(u)=-{\frac {1}{\lambda }}\ln(1-u).}

Eine zweite Standardzufallszahl wählt aus, welches der paarweise disjunkten Ereignisse stattfinden soll:

  • Eine Infektion geschieht mit Wahrscheinlichkeit β s i λ N {\displaystyle {\tfrac {\beta si}{\lambda N}}} ,
  • ein Übergang ins Infektiöse mit Wahrscheinlichkeit α e λ {\displaystyle {\tfrac {\alpha e}{\lambda }}} ,
  • eine Erholung mit der restlichen Wahrscheinlichkeit γ i λ {\displaystyle {\tfrac {\gamma i}{\lambda }}} .

DTMC-Modell

Epidemie-Verläufe im DTMC-Modell. Aufgetragen ist jeweils die Anzahl der Infektiösen. Zu Beginn sind es zehn von tausend. Die schwarze Kurve zeigt den Erwartungswert, die gestrichelte den Wert im klassischen Modell.

Für einen kleinen Zeitschritt Δ t {\displaystyle \Delta t} entsteht aus dem CTMC-Modell eine Näherung in Form eines zeitdiskreten Markow-Prozesses (DTMC: Discrete time Markov chain). Ist der Zeitschritt hinreichend klein, so dass mit ihm höchstens ein Ereignis stattfindet, darf man den Term o ( Δ t ) {\displaystyle o(\Delta t)} ignorieren und erhält die Übergangswahrscheinlichkeiten analog zum CTMC-Modell.[14]

Die Simulation gestaltet sich besonders einfach. Eine Standardzufallszahl wählt aus, welches der paarweise disjunkten Ereignisse im Zeitschritt stattfinden soll:

  • Eine Infektion geschieht mit Wahrscheinlichkeit 1 N β s i Δ t {\displaystyle {\tfrac {1}{N}}\beta si\Delta t} ,
  • ein Übergang ins Infektiöse mit Wahrscheinlichkeit α e Δ t {\displaystyle \alpha e\Delta t} ,
  • eine Erholung mit Wahrscheinlichkeit γ i Δ t {\displaystyle \gamma i\Delta t} ,
  • Verharren im Zustand mit der restlichen Wahrscheinlichkeit.

Weblinks

  • Modellrechner für den Verlauf einer Epidemie nach dem SEIR-Modell mit änderbaren Input-Variablen (englisch)

Siehe auch

Literatur

  • Matt J. Keeling, Pejman Rohani: Modeling Infectious Diseases in Humans and Animals. Princeton University Press, 2008.
  • Matthias an der Heiden, Udo Buchholz: Modellierung von Beispielszenarien der SARS-CoV-2-Epidemie 2020 in Deutschland. Robert Koch-Institut, Abteilung für Infektionsepidemiologie (20. März 2020). doi:10.25646/6571.2.

Einzelnachweise

  1. Institute for Disease Modeling: Originals vom 7. Mai 2020 im Internet Archive)  Info: Der Archivlink wurde automatisch eingesetzt und noch nicht geprüft. Bitte prüfe Original- und Archivlink gemäß Anleitung und entferne dann diesen Hinweis.@1@2Vorlage:Webachiv/IABot/www.idmod.org (abgerufen am 12. April 2020)
  2. a b Pauline van den Driessche: Reproduction numbers of infectious disease models. In: Infectious Disease Modelling, Band 2, KeAi Publishing, August 2017, S. 288–303. doi:10.1016/j.idm.2017.06.002. PMC 6002118 (freier Volltext).
  3. Odo Diekmann, Hans Heesterbeek, Tom Britton: Mathematical tools for understanding infectious disease dynamics. Princeton University Press, 2013, S. 35.
  4. Odo Diekmann, J. A. P. Heesterbeek, J. A. J. Metz: On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. In: Journal of Mathematical Biology, Band 28, 1990, S. 365–382. doi:10.1007/BF00178324.
  5. Pauline van den Driessche, James Watmough: Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. In: Mathematical Biosciences, Band 180, Nr. 1–2, November–December 2002, S. 29–48. doi:10.1016/s0025-5564(02)00108-6.
  6. Junling Ma: Estimating epidemic exponential growth rate and basic reproduction number. In: Infectious Disease Modelling, Band 5, 2020, S. 129–141. doi:10.1016/j.idm.2019.12.009.
  7. Stellungnahme der Deutschen Gesellschaft für Epidemiologie (DGEpi) zur Verbreitung des neuen Coronavirus (SARS-CoV-2). (PDF) Deutsche Gesellschaft für Epidemiologie, 18. März 2020, abgerufen am 26. März 2020.  Die der Beispielrechnung zugrunde liegende Version vom 19. März ist nicht mehr auffindbar. Hier ersatzweise die Version vom Vortag. Die aktuelle Version siehe hier (Memento des Originals vom 29. Mai 2020 im Internet Archive)  Info: Der Archivlink wurde automatisch eingesetzt und noch nicht geprüft. Bitte prüfe Original- und Archivlink gemäß Anleitung und entferne dann diesen Hinweis.@1@2Vorlage:Webachiv/IABot/www.dgepi.de.
  8. Gerardo Chowell, Cécile Viboud, Lone Simonsen, Seyed M. Moghadas: Characterizing the reproduction number of epidemics with early subexponential growth dynamics. In: Journal of the Royal Society Interface, Band 13, Nr. 123, 17. August 2016. doi:10.1098/rsif.2016.0659.
  9. M. Keeling, P. Rohani: Modeling Infectious Diseases in Humans and Animals. Abschnitt 5.2. (S. 159): Modeling forcing in childhood infectious diseases: Measles.
  10. Patrick Hunziker: Personalized-dose Covid-19 vaccination in a wave of virus Variants of Concern: Trading individual efficacy for societal benefit. In: Precision Nanomedicine. 24. Juli 2021, ISSN 2639-9431, doi:10.33218/001c.26101 (precisionnanomedicine.com [abgerufen am 18. August 2021]). 
  11. Patrick Hunziker: Vaccination strategies for minimizing loss of life in Covid-19 in a Europe lacking vaccines. Infectious Diseases (except HIV/AIDS), 1. Februar 2021, doi:10.1101/2021.01.29.21250747 (medrxiv.org [abgerufen am 18. August 2021]). 
  12. Fred Brauer: Mathematical epidemiology: Past, present, and future. In: Infectious Disease Modelling, Band 2, Nr. 2, 2017, S. 113–127. doi:10.1016/j.idm.2017.02.001. PMC 6001967 (freier Volltext).
  13. William Tritch, Linda J. S. Allen: Duration of a minor epidemic. In: Infectious Disease Modelling, Band 3, 2018, S. 60–73. doi:10.1016/j.idm.2018.03.002. PMC 6326226 (freier Volltext).
  14. Linda J. S. Allen: An Introduction to Stochastic Epidemic Models. In: Fred Brauer, Pauline van den Driessche, Jianhong Wu: Mathematical Epidemiology. In: Lecture Notes in Mathematics, Band 1945. Springer, Berlin/Heidelberg 2008, S. 81–130. doi:10.1007/978-3-540-78911-6 3.