OpenModelicaで感染症の数理モデル「SIRモデル」を実装する

SIRモデルについて

2020年4月7日にCOVID-19の流行に伴う緊急事態宣言がされた際、北大の西浦先生が接触率を20%減、60%減、80%減した時の試算をグラフで示されていました。
この時使われた「SIRモデル」とはどういうものか調べてみると、Python や R で計算されている例があり、同じ式を OpenModelica で実装できるのではと試してみたところ、シミュレーションできましたので記事にしておきます。

なお、この計算で得られた結果はCOVID-19の実際の状況を反映しているものではなく、医療従事者でもないエンジニアが仮のパラメータでモデル化しただけのものであることをお断りしておきます。

モデル化にあたり、以下のサイト、記事を参考にさせていただきました。

感染病の数学予測モデルの紹介　(SIRモデル）
https://qiita.com/kotai2003/items/3078f4095c3e94e5325c
微分方程式と感染症数理疫学
https://www.ms.u-tokyo.ac.jp/~inaba/inaba_science_2008.pdf

OpenModelica で基本の SIRモデルを実装

SIRモデルは次の3つの微分方程式で表されます。

// 「感染病の数学予測モデルの紹介　(SIRモデル）」から引用

\begin{align}

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

\end{align} 

\begin{align}
S &: 感染可能者 \quad \text{(Susceptible)}　\\
I &: 感染者 \quad \text{(Infectious)}　\\
R &: 感染後死亡者、もしくは免疫を獲得した者 \quad \text{(Removed)}　\\
\beta &: 感染率\quad \text{(The infectious rate)} \quad [1/day] \\
\gamma &:除去率\quad \text{(The Recovery rate)} \quad [1/day]　\\
\end{align}

このSIRモデルの前提条件
- 一度免疫を獲得した者は2度と感染することなく、免疫を失うこともない。
- 全体人口で外部からの流入及び流出はない。なお出生者及び感染以外の原因で死亡する者もいない。

// 引用ここまで

この式を使って OpenModelica 上で基本になる SIRモデルを実装してみます。

パラメータとして感染率 β、γ を定義し、入力として S、I、R の初期値、出力として S、I、R の変化値を設けました。
コードを示します。

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;

cl_SIRsimple クラスを使ってシミュレーション

初期値を定数で与えるモデルを作成して動かしてみます。
初期値は全人口を 1000人として、S = 999、I = 1、R = 0 とします。
パラメータ β = 0.5/1000、γ = 0.2 としました。

以下のグラフが得られ、他の言語で計算されている SIRモデルと同等の結果となりました。
単位時間 (s) は day と置き換えて考えてください。

感染率を外部から与えるクラス cl_SIRbetaInput を作成

西浦モデルをシミュレーションするために、β値を変化させるための入力インターフェイスを追加し、さらに「新規感染者数（1日あたりの増分）」の出力も追加します。

まず、パラメータになっていた contact_rate をコメントアウトし、Inputから受け取るようにします。Inew が新規感染者数の Output。

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

等式部分に Inew を追加。
「新規感染者数 = 感染可能者の減分」ということです。

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

コード全体はこうなります。

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;

接触率を変化させたシミュレーション結果

これを使って、12日目に接触率を変化させてシミュレーションしてみます。
「SIR」アイコンの下向き三角が新規感染者数の出力コネクタです。

パラメータの合わせ込みが十分ではないと思いますが、「西浦モデル」の新規感染者数の推移のグラフが再現できたのではないかと思います。

まとめ

OpenModelica で感染症の SIRモデルのシミュレーションができました。
ワクチンの効果が見られるモデルや回復者数に死亡率を与えるモデル、再感染を考慮したモデルも作ってみましたので、別記事でご紹介したいと思います。