Implémenter le modèle mathématique «modèle SIR» des maladies infectieuses avec Open Modelica

À propos du modèle SIR

Lorsque la déclaration d'urgence due à l'épidémie de COVID-19 a été déclarée le 7 avril 2020, le Dr Nishiura de l'Université de Kita a montré un graphique montrant l'estimation lorsque le taux de contact diminuait de 20%, 60% et 80%. Était là. En regardant ce qu'est le "modèle SIR" utilisé actuellement, il existe des exemples de calculs en Python et R, et la même expression peut être implémentée dans OpenModelica. Quand je l'ai essayé, j'ai pu le simuler, donc je vais l'écrire sous forme d'article.

Veuillez noter que les résultats obtenus à partir de ce calcul ne reflètent pas la situation réelle du COVID-19, et ne sont modélisés que par des ingénieurs qui ne sont pas des professionnels de la santé avec des paramètres temporaires. Mettre.

En modélisation, j'ai fait référence aux sites et articles suivants.

** Introduction d'un modèle de prédiction mathématique pour les maladies infectieuses (modèle SIR) ** https://qiita.com/kotai2003/items/3078f4095c3e94e5325c ** Équations différentielles et épidémies mathématiques de maladies infectieuses ** https://www.ms.u-tokyo.ac.jp/~inaba/inaba_science_2008.pdf

Implémenter le modèle SIR de base avec OpenModelica

Le modèle SIR est représenté par les trois équations différentielles suivantes.

// * Extrait de "Introduction of Mathematical Prediction Model for Infectious Diseases (SIR Model)" *

\begin{align}

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

\end{align} 
\begin{align}
S &:Personne infectable\quad \text{(Susceptible)} \\
I &:Personne infectée\quad \text{(Infectious)} \\
R &:Ceux qui sont décédés après une infection ou qui ont acquis une immunité\quad \text{(Removed)} \\
\beta &:Taux d'infection\quad \text{(The infectious rate)} \quad [1/day] \\
\gamma &:Taux d'exclusion\quad \text{(The Recovery rate)} \quad [1/day] \\
\end{align}

Prérequis pour ce modèle SIR

// Citer jusqu'à ici

Utilisons cette formule pour implémenter le modèle SIR de base sur OpenModelica.

Les taux d'infection β et γ sont définis comme des paramètres, les valeurs initiales de S, I et R sont fournies en tant qu'entrées, et les valeurs de changement de S, I et R sont fournies en tant que sorties. Voici le 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 à l'aide de la classe cl_SIRsimple

Créons un modèle qui donne la valeur initiale sous forme de constante et déplacez-la. La valeur initiale est S = 999, I = 1, R = 0, en supposant que la population totale est de 1000. Les paramètres β = 0,5 / 1000 et γ = 0,2. 01_SIR_Simple_Model.png

Le graphique suivant est obtenu, ce qui équivaut au modèle SIR calculé dans d'autres langues. Considérez les heures de l'unité comme remplaçant le jour. 02_SIR_Simple_Graph.png

Créez une classe cl_SIRbetaInput qui donne le taux d'infection en externe

Afin de simuler le modèle Nishiura, nous ajouterons une interface d'entrée pour changer la valeur β, et ajouterons également une sortie "Nombre de personnes nouvellement infectées (incrémentation par jour)".

Commencez par mettre en commentaire le contact_rate qui était un paramètre afin qu'il soit reçu de Input. Inew est la sortie du nombre de personnes nouvellement infectées.

//  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)));

Ajout de Inew à la partie équation. Cela signifie "le nombre de personnes nouvellement infectées = réduction du nombre de personnes pouvant être infectées".

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

Le code entier ressemble à ceci.

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;

Résultats de simulation avec des taux de contact variables

En utilisant cela, simulons en modifiant le taux de contact le 12e jour. Le triangle descendant de l'icône "SIR" est le connecteur de sortie pour le nombre de personnes nouvellement infectées. 03_SIR_betainput_model.png

Je ne pense pas que les paramètres aient été suffisamment ajustés, mais je pense que nous avons pu reproduire le graphique de l'évolution du nombre de personnes nouvellement infectées dans le «modèle Nishiura». 06_SIR_betainput_graph_60-80per.png

Résumé

OpenModelica a pu simuler un modèle SIR d'une maladie infectieuse. J'ai aussi fait un modèle qui montre l'effet du vaccin, un modèle qui donne un taux de mortalité au nombre de guérisseurs, et un modèle qui considère la réinfection, je voudrais donc l'introduire dans un autre article.

Recommended Posts

Implémenter le modèle mathématique «modèle SIR» des maladies infectieuses avec Open Modelica
Implémenter le modèle mathématique «modèle SIR» des maladies infectieuses avec OpenModelica (exemple de régulation répétée et de relaxation)
Implémenter le modèle mathématique «modèle SIR» des maladies infectieuses dans OpenModelica (voir l'effet du vaccin)
Mettre en œuvre le modèle mathématique «modèle SIR» des maladies infectieuses dans OpenModelica (reflétant le taux de mortalité et de réinfection)
Introduction d'un modèle de prédiction mathématique pour les maladies infectieuses (modèle SIR)
Modèle de prédiction mathématique des maladies infectieuses (modèle SIR): étude de cas (1)
Modèle mathématique des épidémies de maladies infectieuses
Prédire les épidémies de maladies infectieuses avec le modèle SIR
[Introduction au modèle SIR] Prédire l'heure de fin de chaque pays avec l'ajustement des données COVID-19 ♬
Calibrer le modèle avec PyCaret
Analysez le modèle thématique pour devenir romancier avec GensimPy3
J'ai essayé de prédire le nombre de personnes infectées au niveau national de la nouvelle corona avec un modèle mathématique
[Objet obligatoire DI] Implémenter et comprendre le mécanisme de DI avec Go
[Introduction au modèle SIR] Considérez le résultat de l'ajustement de Diamond Princess ♬
Simulation Python du modèle épidémique (modèle Kermack-McKendrick)
Validez le modèle d'entraînement avec Pylearn2
Implémenter un modèle avec état et comportement (3) - Exemple d'implémentation par décorateur
J'ai essayé de prédire le comportement du nouveau virus corona avec le modèle SEIR.
Alignez la taille de la barre de couleurs avec matplotlib
Affinons les hyper paramètres du modèle avec scikit-learn!
Implémenter une partie du processus en C ++
Vérifier l'existence du fichier avec python
Attention Seq2 Exécutez le modèle de dialogue avec Seq
Transformez le programme de calcul du modèle SIR en une interface graphique
Somme des variables dans un modèle mathématique
La troisième nuit de la boucle avec pour
La deuxième nuit de la boucle avec pour
Implémenter un modèle avec état et comportement
Compter le nombre de caractères avec écho
J'ai essayé de créer un modèle avec l'exemple d'Amazon SageMaker Autopilot