Implement the mathematical model "SIR model" of infectious diseases with OpenModelica

About the SIR model

When a state of emergency was declared due to the COVID-19 epidemic on April 7, 2020, Dr. Nishiura of Hokkaido University showed a graph showing the estimation when the contact rate was reduced by 20%, 60%, and 80%. Was there. Looking at what the "SIR model" used at this time is, there are examples of calculations in Python and R, and the same expression can be implemented in OpenModelica. When I tried it, I was able to simulate it, so I will write it as an article.

Please note that the results obtained from this calculation do not reflect the actual situation of COVID-19, and are only modeled by engineers who are not medical professionals with temporary parameters. Put.

In modeling, I referred to the following sites and articles.

** Introduction of mathematical prediction model for infectious diseases (SIR model) ** https://qiita.com/kotai2003/items/3078f4095c3e94e5325c ** Differential equations and mathematical epidemiology of infectious diseases ** https://www.ms.u-tokyo.ac.jp/~inaba/inaba_science_2008.pdf

Implementing a basic SIR model with OpenModelica

The SIR model is represented by the following three differential equations.

// * Quoted from "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 &:Infectable person\quad \text{(Susceptible)} \\
I &:Infected person\quad \text{(Infectious)} \\
R &:Those who died after infection or who acquired immunity\quad \text{(Removed)} \\
\beta &:Infection rate\quad \text{(The infectious rate)} \quad [1/day] \\
\gamma &:Exclusion rate\quad \text{(The Recovery rate)} \quad [1/day] \\
\end{align}

Prerequisites for this SIR model --A person who has acquired immunity will never be infected and will not lose immunity. --There is no inflow or outflow from the outside in the total population. No one died due to causes other than birth and infection.

// Quote up to here

Let's use this formula to implement the basic SIR model on OpenModelica.

Infection rates β and γ are defined as parameters, initial values of S, I and R are provided as inputs, and change values of S, I and R are provided as outputs. Here is the 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 using the cl_SIRsimple class

Let's create a model that gives the initial value as a constant and move it. The initial value is S = 999, I = 1, R = 0, assuming that the total population is 1000. The parameters β = 0.5 / 1000 and γ = 0.2. 01_SIR_Simple_Model.png

The following graph is obtained, which is equivalent to the SIR model calculated in other languages. Think of the unit time (s) as replacing day. 02_SIR_Simple_Graph.png

Create a class cl_SIRbetaInput that gives the infection rate externally

In order to simulate the Nishiura model, we will add an input interface for changing the β value, and also add an output of "number of newly infected people (increment per day)".

First, comment out the contact_rate that was a parameter so that it will be received from Input. Inew is the Output of the number of newly infected people.

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

Added Inew to the equation part. It means "the number of newly infected people = the reduction of the number of people who can be infected".

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

The whole code looks like this.

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;

Simulation results with varying contact rates

Using this, let's simulate by changing the contact rate on the 12th day. The downward triangle of the "SIR" icon is the output connector for the number of newly infected people. 03_SIR_betainput_model.png

I don't think the parameters have been adjusted sufficiently, but I think we were able to reproduce the graph of changes in the number of newly infected people in the "Nishiura model." 06_SIR_betainput_graph_60-80per.png

Summary

OpenModelica has been able to simulate an SIR model of an infectious disease. I also made a model that shows the effect of the vaccine, a model that gives a mortality rate to the number of recoverers, and a model that considers reinfection, so I would like to introduce it in another article.

Recommended Posts

Implement the mathematical model "SIR model" of infectious diseases with OpenModelica
Implement the mathematical model "SIR model" of infectious diseases with OpenModelica (example of repeating regulation and deregulation)
Implement the mathematical model "SIR model" of infectious diseases in OpenModelica (see the effect of the vaccine)
Implement the mathematical model "SIR model" of infectious diseases in OpenModelica (reflecting mortality rate and reinfection rate)
Introduction of mathematical prediction model for infectious diseases (SIR model)
Mathematical prediction model for infectious diseases (SIR model): Case study (1)
Mathematical model of infectious disease epidemics
Predict infectious disease epidemics with SIR model
[Introduction to SIR model] Predict the end time of each country with COVID-19 data fitting ♬
Calibrate the model with PyCaret
Analyze the topic model of becoming a novelist with GensimPy3
I tried to predict the number of domestically infected people of the new corona with a mathematical model
[Required subject DI] Implement and understand the mechanism of DI with Go
[Introduction to SIR model] Consider the fitting result of Diamond Princess ♬
Python-Simulation of the Epidemic Model (Kermack-McKendrick Model)
Validate the learning model with Pylearn2
Implement a model with state and behavior (3) --Example of implementation by decorator
I tried to predict the behavior of the new coronavirus with the SEIR model.
Align the size of the colorbar with matplotlib
Let's tune the model hyperparameters with scikit-learn!
Implement part of the process in C ++
Check the existence of the file with python
Run the interaction model with Attention Seq2 Seq
Make the SIR model calculation program GUI
Sum of variables in a mathematical model
The third night of the loop with for
The second night of the loop with for
Reading comprehension of "The Go Memory Model"
Implement a model with state and behavior
Count the number of characters with echo
I tried to create a model with the sample of Amazon SageMaker Autopilot