Über das SIR-Modell

Als die Notfallerklärung aufgrund der COVID-19-Epidemie am 7. April 2020 erklärt wurde, zeigte Dr. Nishiura von der Kita University eine Grafik, die die Schätzung zeigt, wenn die Kontaktrate um 20%, 60% und 80% reduziert wurde. War dort. In Bezug auf das derzeit verwendete "SIR-Modell" gibt es Beispiele für Berechnungen in Python und R, und der gleiche Ausdruck kann in OpenModelica implementiert werden. Als ich es ausprobiert habe, konnte ich es simulieren, also werde ich es als Artikel schreiben.

Bitte beachten Sie, dass die Ergebnisse dieser Berechnung nicht die tatsächliche Situation von COVID-19 widerspiegeln und nur von Ingenieuren modelliert werden, die keine Mediziner mit temporären Parametern sind. Stellen.

Bei der Modellierung habe ich auf die folgenden Websites und Artikel verwiesen.

** Einführung eines mathematischen Vorhersagemodells für Infektionskrankheiten (SIR-Modell) ** https://qiita.com/kotai2003/items/3078f4095c3e94e5325c ** Differentialgleichungen und mathematische Epidemien von Infektionskrankheiten ** https://www.ms.u-tokyo.ac.jp/~inaba/inaba_science_2008.pdf

Implementieren Sie das grundlegende SIR-Modell mit OpenModelica

Das SIR-Modell wird durch die folgenden drei Differentialgleichungen dargestellt.

// * Zitiert aus "Einführung des mathematischen Vorhersagemodells für Infektionskrankheiten (SIR-Modell)" *

\begin{align}

\frac{dS}{dt} &= -\beta SI \\
\frac{dI}{dt} &=  \beta SI -\gamma I \\
\frac{dR}{dt} &=  \gamma I \\

\end{align}

\begin{align}
S &:Infizierbare Person\quad \text{(Susceptible)}　\\
I &:Infizierte Person\quad \text{(Infectious)}　\\
R &:Diejenigen, die nach einer Infektion starben oder Immunität erlangten\quad \text{(Removed)}　\\
\beta &:Infektionsrate\quad \text{(The infectious rate)} \quad [1/day] \\
\gamma &:Ausschlussrate\quad \text{(The Recovery rate)} \quad [1/day]　\\
\end{align}

Voraussetzungen für dieses SIR-Modell ――Eine Person, die Immunität erworben hat, wird nie wieder infiziert und verliert keine Immunität.

In der Gesamtbevölkerung gibt es keinen Zu- oder Abfluss von außen. Niemand starb aufgrund anderer Ursachen als Geburt und Infektion.

// Bis hierher zitieren

Verwenden wir diese Formel, um das grundlegende SIR-Modell auf OpenModelica zu implementieren.

Die Infektionsraten β und γ werden als Parameter definiert, Anfangswerte von S, I und R werden als Eingaben bereitgestellt und Änderungswerte von S, I und R werden als Ausgaben bereitgestellt. Hier ist der Code.

class cl_SIRsimple
  extends Modelica.Blocks.Icons.Block;
  parameter Real contact_rate = 0.5 / 1000 "Infectious rate (1/day)";
  parameter Real gamma_const = 0.2 "Recovery rate (1/day)";
  Modelica.Blocks.Interfaces.RealInput Si "Connector of initial S(Susceptible)" annotation(
    Placement(visible = true, transformation(origin = {-100, 60}, extent = {{-20, -20}, {20, 20}}, rotation = 0), iconTransformation(origin = {-120, 60}, extent = {{-20, -20}, {20, 20}}, rotation = 0)));
  Modelica.Blocks.Interfaces.RealInput Ii "Connector of initial I (Infectious)" annotation(
    Placement(visible = true, transformation(origin = {-100, 0}, extent = {{-20, -20}, {20, 20}}, rotation = 0), iconTransformation(origin = {-120, 0}, extent = {{-20, -20}, {20, 20}}, rotation = 0)));
  Modelica.Blocks.Interfaces.RealInput Ri "Connector of initial R (Removed)" annotation(
    Placement(visible = true, transformation(origin = {-100, -62}, extent = {{-20, -20}, {20, 20}}, rotation = 0), iconTransformation(origin = {-120, -60}, extent = {{-20, -20}, {20, 20}}, rotation = 0)));
  Modelica.Blocks.Interfaces.RealOutput Sy "Output connector of S(Susceptible)" annotation(
    Placement(visible = true, transformation(origin = {100, 60}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {110, 60}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
  Modelica.Blocks.Interfaces.RealOutput Iy "Output connector of I (Infectious)" annotation(
    Placement(visible = true, transformation(origin = {100, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {110, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
  Modelica.Blocks.Interfaces.RealOutput Ry "Output connector of R (Removed)" annotation(
    Placement(visible = true, transformation(origin = {100, -60}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {110, -60}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
initial equation
  Sy = Si;
  Iy = Ii;
  Ry = Ri;
equation
  der(Sy) = -contact_rate * Sy * Iy;
  der(Iy) = contact_rate * Sy * Iy - Iy * gamma_const;
  der(Ry) = Iy * gamma_const;
  annotation(
    Icon(graphics = {Text(origin = {-14, 11}, extent = {{-56, 33}, {78, -51}}, textString = "SIR")}));
end cl_SIRsimple;

Simulation mit der Klasse cl_SIRsimple

Lassen Sie uns ein Modell erstellen, das den Anfangswert als Konstante angibt, und ihn verschieben. Der Anfangswert ist S = 999, I = 1, R = 0, vorausgesetzt, die Gesamtbevölkerung beträgt 1000. Die Parameter β = 0,5 / 1000 und γ = 0,2.

Es wird das folgende Diagramm erhalten, das dem in anderen Sprachen berechneten SIR-Modell entspricht. Stellen Sie sich die Zeiteinheit (en) als Ersatztag vor.

Erstellen Sie eine Klasse cl_SIRbetaInput, die die Infektionsrate extern angibt

Um das Nishiura-Modell zu simulieren, fügen wir eine Eingabeschnittstelle zum Ändern des β-Werts sowie eine Ausgabe von "Anzahl neu infizierter Personen (Inkrement pro Tag)" hinzu.

Kommentieren Sie zunächst die contact_rate aus, die ein Parameter war, damit sie von Input empfangen wird. Inew ist die Ausgabe der Anzahl neu infizierter Personen.

//  parameter Real contact_rate = 0.2/1000 "The infectious rate (1/day)";
parameter Real gamma_const = 0.2 "The Recovery rate (1/day)";

Modelica.Blocks.Interfaces.RealInput contact_rate "Connector of infectious rate" annotation(
  Placement(visible = true, transformation(origin = {0, 100}, extent = {{-20, -20}, {20, 20}}, rotation = -90), iconTransformation(origin = {0, 120}, extent = {{-20, -20}, {20, 20}}, rotation = -90)));
Modelica.Blocks.Interfaces.RealOutput Inew "Output connector of New I(Infectious)" annotation(
  Placement(visible = true, transformation(origin = {0, -100}, extent = {{-10, -10}, {10, 10}}, rotation = -90), iconTransformation(origin = {0, -110}, extent = {{-10, -10}, {10, 10}}, rotation = -90)));

Neu zum Gleichungsteil hinzugefügt. Es bedeutet "die Anzahl der neu infizierten Personen = Verringerung der Anzahl der Personen, die infiziert werden können".

equation
  der(Sy) = -contact_rate * Sy * Iy;
  der(Iy) = contact_rate * Sy * Iy - Iy * gamma_const;
  der(Ry) = Iy * gamma_const;
  Inew = -der(Sy);

Der ganze Code sieht so aus.

class cl_SIRbetaInput
  extends Modelica.Blocks.Icons.Block;
    //  parameter Real contact_rate = 0.2/1000 "The infectious rate (1/day)";
    parameter Real gamma_const = 0.2 "The Recovery rate (1/day)";
    Modelica.Blocks.Interfaces.RealInput Si "Connector of initial S(Susceptible)" annotation(
      Placement(visible = true, transformation(origin = {-100, 60}, extent = {{-20, -20}, {20, 20}}, rotation = 0), iconTransformation(origin = {-120, 60}, extent = {{-20, -20}, {20, 20}}, rotation = 0)));
    Modelica.Blocks.Interfaces.RealInput Ii "Connector of initial I(Infectious)" annotation(
      Placement(visible = true, transformation(origin = {-100, 0}, extent = {{-20, -20}, {20, 20}}, rotation = 0), iconTransformation(origin = {-120, 0}, extent = {{-20, -20}, {20, 20}}, rotation = 0)));
    Modelica.Blocks.Interfaces.RealInput Ri "Connector of initial R(Removed)" annotation(
      Placement(visible = true, transformation(origin = {-100, -62}, extent = {{-20, -20}, {20, 20}}, rotation = 0), iconTransformation(origin = {-120, -60}, extent = {{-20, -20}, {20, 20}}, rotation = 0)));
    Modelica.Blocks.Interfaces.RealOutput Sy "Output connector of S(Susceptible)" annotation(
      Placement(visible = true, transformation(origin = {100, 60}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {110, 60}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
    Modelica.Blocks.Interfaces.RealOutput Iy "Output connector of I(Infectious)" annotation(
      Placement(visible = true, transformation(origin = {100, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {110, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
    Modelica.Blocks.Interfaces.RealOutput Ry "Output connector of R(Removed)" annotation(
      Placement(visible = true, transformation(origin = {100, -60}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {110, -60}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
    Modelica.Blocks.Interfaces.RealInput contact_rate "Connector of infectious rate" annotation(
      Placement(visible = true, transformation(origin = {0, 100}, extent = {{-20, -20}, {20, 20}}, rotation = -90), iconTransformation(origin = {0, 120}, extent = {{-20, -20}, {20, 20}}, rotation = -90)));
    Modelica.Blocks.Interfaces.RealOutput Inew "Output connector of New I(Infectious)" annotation(
      Placement(visible = true, transformation(origin = {0, -100}, extent = {{-10, -10}, {10, 10}}, rotation = -90), iconTransformation(origin = {0, -110}, extent = {{-10, -10}, {10, 10}}, rotation = -90)));
  initial equation
    Sy = Si;
    Iy = Ii;
    Ry = Ri;
  equation
    der(Sy) = -contact_rate * Sy * Iy;
    der(Iy) = contact_rate * Sy * Iy - Iy * gamma_const;
    der(Ry) = Iy * gamma_const;
    Inew = -der(Sy);
    annotation(
      Icon(graphics = {Text(origin = {-14, 11}, extent = {{-56, 33}, {78, -51}}, textString = "SIR")}, coordinateSystem(initialScale = 0.1)));
end cl_SIRbetaInput;

Simulationsergebnisse mit unterschiedlichen Kontaktraten

Simulieren wir damit, indem wir die Kontaktrate am 12. Tag ändern. Das nach unten gerichtete Dreieck des Symbols "SIR" ist der Ausgangsanschluss für die Anzahl der neu infizierten Personen.

Ich denke, die Parameter wurden nicht ausreichend angepasst, aber ich denke, wir konnten das Diagramm der Änderungen der Anzahl neu infizierter Personen im "Nishiura-Modell" reproduzieren.

Zusammenfassung

OpenModelica konnte ein SIR-Modell einer Infektionskrankheit simulieren. Ich habe auch ein Modell erstellt, das die Wirkung des Impfstoffs zeigt, ein Modell, das die Sterblichkeitsrate der Anzahl der Wiederhersteller angibt, und ein Modell, das eine Reinfektion berücksichtigt. Daher möchte ich es in einem anderen Artikel vorstellen.