[PYTHON] Ich habe versucht, die Anzahl der im Inland infizierten Menschen der neuen Korona mit einem mathematischen Modell vorherzusagen

Überblick

Ich fand The Lancets Artikel [^ 1], der die Anzahl infizierter Menschen mit einem mathematischen Modell vorhersagt. Deshalb implementierte ich ein mathematisches Modell, um die Anzahl infizierter Menschen in Japan vorherzusagen, und verglich es mit anderen Infektionskrankheiten.

Zusammenfassung

Die Anzahl der Infizierten wurde vom SEIR-Modell vorhergesagt. ――Die Grundreproduktionszahl der neuen Korona liegt bei etwa 2,9, was etwas höher ist als bei Influenza.
Wenn die Infektionskontrolle nicht wirksam ist, liegt die Spitzenzahl der im Inland infizierten Personen bei etwa GW. ――Die Zahl der Todesfälle im Inland wird im nächsten Jahr über 10.000 liegen (wenn die Letalitätsrate 0,01% beträgt).
Die Einwanderungsbeschränkungen für Ausländer sind bedeutungslos.

Inhaltsverzeichnis

[Zweck dieser Zeit](# 1-Zweck dieser Zeit)
[Mathematisches Modell der Infektionskrankheit](# 2-Mathematisches Modell der Infektionskrankheit)
[Parameterschätzung](# 3-Parameterschätzung)
[In Python implementiert und illustriert](# 4-Implementiert und in Python illustriert)
[Ergänzung](# 5-Ergänzung)

1. Zweck dieser Zeit

Prognostizieren Sie die Anzahl zukünftiger häuslicher Infektionen des neuen Koronavirus (COVID-19, im Folgenden als neue Korona bezeichnet) mithilfe eines mathematischen Modells.
Vergleichen Sie die Grundreproduktionszahlen der neuen Korona mit anderen Infektionskrankheiten.
Implementierung wissenschaftlicher Berechnungen mit Scipy.

2. Mathematisches Modell einer Infektionskrankheit

"Formulieren Sie den Prozess, wie sich eine Infektionskrankheit ausbreitet, wie lange sich eine infizierte Person entwickelt und schwer wird." [^ 2]

Im versicherungsmathematischen Modell von Infektionskrankheiten wird die Population nach dem Infektionsstadium unterteilt, und die Zunahme oder Abnahme jeder Gruppe wird durch eine Differentialgleichung ausgedrückt. Ich werde das einfache SEIR-Modell und die diesmal verwendete erweiterte Version erläutern.

2.1 SEIR-Modell

$ S $: Anzahl der Personen ohne Immunität (anfällig)
$ E $: Anzahl der Personen in der Zeit nach der Exposition (exponiert)
$ I $: Anzahl der Personen während der Infektionsperiode (Infektiös)
$ R $: Anzahl der Personen, die sich von einer Infektion erholt und Immunität erlangt haben (wiederhergestellt)
$ \ beta $: Infektionsrate, $ \ gamma $: Wiederherstellungsrate (Isolationsrate), $ \ sigma $: Wahrscheinlichkeit der Infektiosität, $ N $: Gesamtbevölkerung

Dann werden die folgenden simultanen Differentialgleichungen erhalten.

\begin{align}
\frac{dS(t)}{dt} &= -\beta \frac{I(t)}{N}S(t) \\

\frac{dE(t)}{dt} &= \beta \frac{I(t)}{N}S(t) - \sigma E(t) \\

\frac{dI(t)}{dt} &= \sigma E(t) - \gamma I(t) \\

\frac{dR(t)}{dt} &= \gamma I(t) \\

N &= S(t) + E(t) + I(t) + R(t)
\end{align}

Es ignoriert jedoch die Geburtenrate und die Sterblichkeitsrate und sagt, dass es nie wieder infiziert wird, da es sich nach dem Einsetzen erholt und Immunität erlangt.

Wenn Sie einen geeigneten Anfangswert angeben und diesen lösen, können Sie die Anzahl der infizierten Personen in der Zukunft vorhersagen. (Ich werde die Abbildung [^ 3] als Referenz setzen)

Ausführliche Erklärungen und Implementierungen finden Sie in Qiitas Artikel [^ 4].

Hier wird $ R_0 = \ beta / \ gamma $ als Grundreproduktionsnummer bezeichnet und "(alle Personen sind anfangs anfällig. Es bedeutet "die Anzahl der pro infizierter Person produzierten sekundär infizierten Personen (im Zustand des Habens)". [^ 5] ** Der Punkt ist, dass die Infektiosität umso größer ist, je größer $ R_0 $ ist. ** ** **

Sie können auch $ R_0 $ verwenden, um zu berechnen, wie eine weit verbreitete Impfung eine Epidemie verhindern kann.

Zusätzlich zu SEIR können verschiedene Modelle erstellt werden, indem die gesamte Bevölkerung nach Wohnort und Alter unterteilt wird und die Merkmale von Infektionskrankheiten in Differentialgleichungen wiedergegeben werden. Zum Beispiel gibt es SIR, das E nicht berücksichtigt, SEIS, das den Erwerb der Immunität nicht berücksichtigt, und MSEIR, das die passive Immunität berücksichtigt [^ 6].

2.2 Erweitertes SEIRD-Modell

Da das obige SEIR zu einfach ist, werden wir es diesmal auf ein Modell ausweiten, das den Zu- und Abfluss der Bevölkerung und die Sterblichkeitsrate im Ausland unter Bezugnahme auf das Papier [^ 1] berücksichtigt.

$ R_0 : Grundlegende Reproduktionsnummer ( R_0 = \ beta / \ gamma $)
$ T_E : Durchschnittlicher Zeitraum nach der Exposition (Tage) ( \ gamma ^ {-1} = T_I $)
$ T_I : Durchschnittliche Infektiositätsdauer (Tage) ( \ sigma ^ {-1} = T_E $)
$ N_ {J, I} $: Anzahl der Reisenden von Japan nach Übersee (Personen / Tag)
$ N_ {I, J} $: Anzahl der Reisenden aus Übersee nach Japan (Personen / Tag)
$ f $: Tödlichkeit
$ D (t) $: Anzahl der Todesfälle (Personen)

Dann werden auch die folgenden simultanen Differentialgleichungen erhalten.

\begin{align}
\frac{dS(t)}{dt} &= - \frac{R_0}{D_I} \frac{I(t)}{N(t)} S(t) - N_{J,I} + 
 \left(1 - \frac{E(t) + I(t)}{N(t)} \right)N_{I,J} \\

\frac{dE(t)}{dt} &= \frac{R_0}{D_I} \frac{I(t)}{N(t)} S(t) - \frac{E(t)}{T_E} + \frac{E(t) + I(t)}{N(t)} N_{I,J}\\

\frac{dI(t)}{dt} &= \frac{E(t)}{T_E} - \frac{I(t)}{T_I}\\

\frac{dR(t)}{dt} &= (1-f) \frac{I(t)}{T_I}\\

\frac{dD(t)}{dt} &= f \frac{I(t)}{T_I} \\

N(t) &= S(t) + E(t) + I(t) + R(t)
\end{align}

In Bezug auf $ N_ {J, I} $ und $ N_ {I, J} $ wird jedoch angenommen, dass alle, die aus Japan ins Ausland gehen, zu $ S $ gehören, und diejenigen, die aus Übersee nach Japan kommen, zu $ S, E $ gehören. Ich werde.

Wie in der Arbeit [^ 1] werden andere Parameter als $ R_0 $ durch entsprechende Werte ersetzt. $ T_I und T_E $ werden unter Bezugnahme auf SARS Corona auf $ T_E = 6.0 und T_I = 2.4 $ gesetzt. $ N_ {J, I}, N_ {I, J} $ wird auf $ N_ {J, I} = 40000, N_ {I, J} = 70000 $ bezogen auf die Anzahl der Reisenden im Februar 2019 gesetzt.

3. Parameterschätzung

Der unbekannte Parameter des Modells ist $ R_0 $. Es gibt verschiedene Schätzmethoden wie die Methode des minimalen Quadrats, die Schätzung der maximalen Wahrscheinlichkeit, die MAP-Schätzung, die Bayes'sche Schätzung [^ 7], aber dieses Mal wird die Wahrscheinlichkeit als instationärer Poisson-Prozess [^ 8] berechnet und $ R_0 $ wird durch die Schätzung der maximalen Wahrscheinlichkeit berechnet. ..

3.1. Wahrscheinlichste Schätzung

Jeder liebt die wahrscheinlichste Schätzung. Drücken Sie zunächst die Wahrscheinlichkeit (Wahrscheinlichkeit) der beobachteten Daten als Funktion von $ R_0 $ aus und finden Sie $ R_0 $, das sie maximiert. Wenn die Anzahl der infizierten Personen als instationärer Poisson-Prozess unter Bezugnahme auf das Papier [^ 1] modelliert und die Wahrscheinlichkeit berechnet wird,

L(R_0) = \prod_{d=D_s}^{D_e} \frac{e^{-\lambda_d}\lambda_d^{x_d}}{x_d!}

$ D_s $: erster Beobachtungstag, $ D_s $: letzter Beobachtungstag, $ x_d $: Anzahl infizierter Personen bei $ d $

\lambda (t, R_0) = E(t) + I(t) \\
\lambda _d (R_0) = \int_{d-1}^{d} \lambda (u, R_0) du

Wird besorgt. $ \ lambda $ repräsentiert die durchschnittliche Anzahl infizierter Personen, wenn sie der Poisson-Verteilung folgen.

Ohne Begriffe, die nicht mit $ R_0 $ in Zusammenhang stehen, von der logarithmischen Wahrscheinlichkeit kann $ R_0 $ unten geschätzt werden.

\newcommand{\argmin}{\mathop{\rm arg~min}\limits}

\argmin_{R_0 > 0} - \log L(R_0) = \argmin_{R_0 > 0} \sum_{d=D_s}^{D_e} (\lambda_d(R_0) - x_d \log \lambda_d(R_0))

4. In Python implementiert und illustriert

Implementierung der SEIR-Klasse

import scipy
from scipy.integrate import odeint
from scipy.optimize import minimize
from sklearn.metrics import mean_absolute_error, mean_squared_error


class SIR(object):

    def __init__(self, r_0, t_i, dt, init_S, init_I, init_R):
        self.r_0 = r_0
        self.t_i = t_i
        self.dt = dt
        self.init_state_list = [init_S, init_I, init_R]

    def _get_param(self, params, key, default, t=None):
        """Entspricht zeitlich variierenden Parametern
        """
        if isinstance(params, dict):
            if key in list(params.keys()):
                param = params[key]
                if isinstance(param, list):
                    return param[np.clip(int(t / self.dt), 0, len(param)-1)]
                elif isinstance(param, np.ndarray):
                    return param
                else:
                    return default
            else:
                return default
        else:
            return default

    def _ode(self, state_list, t=None, params=None):
        """Definieren Sie simultane Differentialgleichungen
        """
        r_0 = self._get_param(params, 'r_0', self.r_0, t=t)
        t_i = self._get_param(params, 't_i', self.t_i, t=t)
        S, I, R = state_list
        N = S + I + R
        dstate_dt = list()
        dstate_dt[0] = - (r_0 / t_i) * (I / N) * S
        dstate_dt[1] = (r_0 / t_i) * (I / N) * S - I / t_i
        dstate_dt[2] = I / t_i
        return dstate_dt

    def solve_ode(self, len_days=365, params=None):
        """Lösen Sie Differentialgleichungen
        """
        t = np.linspace(0, len_days, int(len_days / self.dt), endpoint=False)
        args = (params,) if params else ()
        return odeint(self._ode, self.init_state_list, t, args=args)


class CustomizedSEIRD(SIR):

    def __init__(self, r_0=None, t_e=6.0, t_i=2.4, n_i_j=70000, n_j_i=40000, f=0.0001, dt=1,
            init_S=126800000, init_E=0, init_I=1, init_R=0, init_D=0):
        self.r_0 = r_0
        self.t_e = t_e
        self.t_i = t_i
        self.n_i_j = n_i_j
        self.n_j_i = n_j_i
        self.f = f
        self.dt = dt
        self.init_state_list = [init_S, init_E, init_I, init_R, init_D]

    def _ode(self, state_list, t=None, params=None):
        """Definieren Sie simultane Differentialgleichungen
        """
        r_0 = self._get_param(params, 'r_0', self.r_0, t=t)
        t_e = self._get_param(params, 't_e', self.t_e, t=t)
        t_i = self._get_param(params, 't_i', self.t_i, t=t)
        n_i_j = self._get_param(params, 'n_i_j', self.n_i_j, t=t)
        n_j_i = self._get_param(params, 'n_j_i', self.n_j_i, t=t)
        f = self._get_param(params, 'f', self.f)
        S, E, I, R, D = state_list
        N = S + E + I + R
        dstate_dt = list()
        dstate_dt.append(- (r_0 / t_i) * (I / N) * S - n_j_i + (1 - ((E + I) / N)) * n_i_j)
        dstate_dt.append((r_0 / t_i) * (I / N) * S - E / t_e + ((E + I) / N) * n_i_j)
        dstate_dt.append(E / t_e - I / t_i)
        dstate_dt.append((1 - f) * I / t_i)
        dstate_dt.append(f * I / t_i)
        return dstate_dt

    def _calc_neg_log_likelihood_r0(self, r_0, X):
        """Log-Wahrscheinlichkeit(R_Nur der Teil, der sich auf 0 bezieht)Berechnen
        """
        solution = self.solve_ode(len_days=len(X), params=dict(r_0=r_0))
        lambda_arr = solution[int(1/self.dt)-1::int(1/self.dt), 2] # I
        return - np.sum(- lambda_arr + X * np.log(lambda_arr))
    
    def _calc_error(self, r_0, X, metric=mean_absolute_error):
        """Berechnen Sie den durchschnittlichen absoluten Fehler und den durchschnittlichen quadratischen Fehler
        """
        solution = self.solve_ode(len_days=len(X), params=dict(r_0=r_0))
        e_arr = solution[int(1/self.dt)-1::int(1/self.dt), 2]  # I
        return metric(X, e_arr)

    def exec_point_estimation(self, init_r_0, X, project='mle'):
        """Punktschätzungsparameter
        """
        if project == 'mle':
            result = minimize(self._calc_neg_log_likelihood_r0, init_r_0, args=(X,), method='Nelder-Mead')
        elif project == 'lad':
            result = minimize(self._calc_error, init_r_0, args=(X, mean_absolute_error), method='Nelder-Mead')
        elif project == 'ls':
            result = minimize(self._calc_error, init_r_0, args=(X, mean_squared_error), method='Nelder-Mead')
        else:
            print(f'Invalid project: {project}')
            return None
        if self.r_0 is None:
            self.r_0 = result.x[0]
        return result

    def exec_map(self):
        """MAP-Schätzung
        """
        pass

    def exec_mcmc(self):
        """MCMC-Methode
        """
        pass

if __name__ == '__main__':
    # 1/16 bis 3/Anzahl der Infizierten bis zu 8
    X = np.array([
        1, 1, 1, 1, 1, 1, 1, 1, 2, 3,  # 1/25
        4, 4, 7, 8, 14, 17, 20, 20, 20, 23,  # 2/4
        35, 45, 86, 90, 96, 161, 163, 203, 251, 259,  # 2/14
        338, 414, 520, 616, 705, 728, 743, 743, 838, 851, # 2/24
        862, 891, 919, 938, 947, 961, 980, 999, 1035, 1056, # 3/5
        1113, 1157, 1190 # 3/8
    ])
    model = CustomizedSEIRD()
    result = model.exec_point_estimation(init_r_0=1, X=X, project='mle')
    print(result.x[0])

2.8806152343750036

Die "Anzahl infizierter Personen" ist schwer zu handhaben, und nur diejenigen, die das Krankenhaus besuchen, werden tatsächlich gezählt, und niemand kennt das wahre $ x_d $ unter Berücksichtigung der diagnostischen Kriterien und der Sensitivitätsspezifität des Tests. Hmm. Hier wird es basierend auf der Anzahl der Infizierten [^ 9] einschließlich Kreuzfahrtschiffen geschätzt.

Vorhersage zukünftiger häuslicher Infektionen durch erweitertes SEIR-Modell (I ist die Anzahl der Infektionen)

Es ist ein Diagramm, wenn sich das durch die blaue Linie geschätzte $ R_0 $ in Zukunft nicht ändert. Wenn wir so weitermachen, wird die Anzahl der infizierten Menschen ihren Höhepunkt erreichen, sobald das GW vorbei ist. Es ist ein Diagramm, wenn sich die rote und die blaue Linie nach 3/9 $ R_0 $ ändern. Da $ R_0 $ "die Anzahl der pro infizierter Person produzierten sekundär infizierten Personen" ist, ist sie aufgrund von Gegenmaßnahmen gegen Infektionen und Bewegungseinschränkungen geringer. Wenn $ R_0 $ klein genug gemacht werden kann, kann es wie eine grüne Linie enthalten sein.

Die folgende Abbildung ist eine vom Ministerium für Gesundheit und Soziales veröffentlichte Zahl [^ 9], aber die Form des Diagramms der Anzahl der infizierten Personen stimmt genau überein.

厚生省コロナ感染者数.png

Wenn Sie beispielsweise $ N_ {i, j} $ ändern, können Sie auch berechnen, ob sich die Anzahl der infizierten Personen aufgrund von Einwanderungsbeschränkungen von Ausländern ändert. Bei Berechnung mit demselben Modell ändert sich die Anzahl der infizierten Personen nicht, selbst wenn $ N_ {i, j} = 0 $ nach 3/9 ist und $ R_0 $ gleich ist. Dies bedeutet, dass Einwanderungsbeschränkungen bedeutungslos sind (wenn dieses Modell korrekt ist).

Vergleichen wir als nächstes die Grundreproduktionszahl und die Letalität mit anderen Infektionskrankheiten [^ 10]. (Die Letalitätsrate der neuen Korona beträgt etwa 0,1% bis 5% der Letalitätsrate jedes Landes.)

Es scheint, dass die Grundreproduktionszahl und die Letalitätsrate der neuen Korona nicht signifikant höher sind als bei anderen Infektionskrankheiten. (Für das Risiko von Infektionskrankheiten ist es notwendig, die Schwere der Symptome und das Vorhandensein oder Fehlen von Impfstoffen zu berücksichtigen, und ich denke, dass niedrige Grundreproduktionszahlen und niedrige Letalitätsraten nicht bedeuten, dass sie sicher sind, aber die Details sind Fragen Sie jemanden, der gegen Infektionskrankheiten resistent ist.

5. Ergänzung

Ministerium für Gesundheit, Arbeit und Soziales

Japanische Gesellschaft für Infektionskrankheiten

WHO

CDC

Coole Website, auf der Sie die Anzahl der infizierten Personen überprüfen können (Johns Hopkins CSSE)

** Ich bin kein Spezialist für Infektionskrankheiten, kein Gesundheitsamt oder Statistiker. Die tatsächliche Infektionskrankheit ist nicht so einfach wie dieses Modell, daher ist diese Vorhersage wahrscheinlich falsch. Wir sind nicht verantwortlich für Schäden, die durch den Inhalt dieses Artikels verursacht werden. ** ** **

[^ 2]: Infektionskrankheiten mit einem mathematischen Modell stoppen

[^ 4]: Beginn des mathematischen Modells für Infektionskrankheiten: Überblick über das SEIR-Modell von Python und Einführung in die Parameterschätzung

[^ 5]: Vorhersage von Infektionskrankheiten: quantitative Probleme in mathematischen Modellen von Infektionskrankheiten

[^ 7]: Höchstwahrscheinlich Schätzung, MAP-Schätzung, Bayes-Schätzung in Fujii 4. Dan gelernt

[^ 8]: SIR-Modell und instationärer Poisson-Prozess

[^ 9]: Über eine neue Coronavirus-Infektion