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.

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.

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.

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."

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.