基本のSIRモデルに死亡率、再感染率を追加してみる
既出記事、「OpenModelicaで感染症の数理モデル「SIRモデル」を実装する」で作成したクラス cl_SIRbetaInput に、死亡率、再感染率を加えたフルバージョン版を作成してみます。
※最初の投稿で死亡率の数式で間違いを見つけたため修正しています。
※※（2020/5/13修正）再度見直しました。
このモデルの前提条件
- 感染者から回復者に移行したうち、一定の割合で死亡に至る
- 一度免疫を獲得した回復者でも時間の経過とともに一定の割合で免疫を失って感染可能者に戻る
- 回復者から感染可能者に戻るのは生存者のみ
なお、以下のシミュレーションは全て仮のパラメータでモデルの検証をしただけであり、実際の数値を反映・予想したものではないことをお断りしておきます。
クラス cl_SIRbetaInput に死亡率を加えたクラス cl_SIRreinfect を作成する。
パラメータとして「dying_rate = 死亡率」を設定可能とします。
デフォルト値は0.01（1%）とします。
parameter Real dying_rate = 0.01 "Dying rate";
追加になる等式部分は以下のとおり
I（感染者）から R（回復者）に移行する時、回復者のうち死亡率に比例して死亡者となり、残りが「生存回復者」となる。
equation
der(Rdead) = Iy * gamma_const * dying_rate;
// der(Rdead) = der(Ry) * dying_rate;
Rlive = Ry - Rdead;
さらに再感染率を加える
再感染を考慮した数式については、「Rとパンデミックの数理モデル 新型コロナウィルス(COVID-19)研究を例に」を参考にさせていただきました。
ρを再感染率とすると下記の数式になります。
\begin{align}
\frac{dS}{dt} &= -\beta SI + ρ R\\
\frac{dI}{dt} &= \beta SI -\gamma I \\
\frac{dR}{dt} &= \gamma I - ρ R\\
\end{align}
パラメータとして「reinfect_rate = 再感染率」を設定可能とします。
デフォルト値は0.001（0.1%）とします。
parameter Real reinfect_rate = 0.001 "Reinfection rate";
追加になる等式部分は以下のとおり
I（感染者）から R（回復者）に移行する時、回復者のうち死亡率に比例して死亡者となり、残りが「生存回復者」となる。
「生存回復者」のうち、再感染率に比例してS（感染可能者）に移行し、その分 R（回復者）が減る。
equation
der(Sy) = (-contact_rate * Sy * Iy) + Rlive * reinfect_rate;
der(Iy) = contact_rate * Sy * Iy - Iy * gamma_const;
der(Ry) = Iy * gamma_const - Rlive * reinfect_rate;
Inew = contact_rate * Sy * Iy;
コード全体はこうなります。
class cl_SIRreinfect
extends Modelica.Blocks.Icons.Block;
parameter Real gamma_const = 0.2 "Recovery rate (1/day)";
parameter Real reinfect_rate = 0.001 "Reinfection rate";
parameter Real dying_rate = 0.01 "Dying rate";
Modelica.Blocks.Interfaces.RealInput Si "Connector of initial S(Susceptible)" annotation(
Placement(visible = true, transformation(origin = {-120, 60}, extent = {{-20, -20}, {20, 20}}, rotation = 0), iconTransformation(origin = {-120, 62}, extent = {{-20, -20}, {20, 20}}, rotation = 0)));
Modelica.Blocks.Interfaces.RealInput Ii "Connector of initial I (Infectious)" annotation(
Placement(visible = true, transformation(origin = {-120, 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 = {-120, -60}, extent = {{-20, -20}, {20, 20}}, rotation = 0), iconTransformation(origin = {-120, -62}, extent = {{-20, -20}, {20, 20}}, rotation = 0)));
Modelica.Blocks.Interfaces.RealOutput Sy "Output connector of S(Susceptible)" annotation(
Placement(visible = true, transformation(origin = {110, 60}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {110, 62}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
Modelica.Blocks.Interfaces.RealOutput Iy "Output connector of I (Infectious)" annotation(
Placement(visible = true, transformation(origin = {110, 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 = {110, -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, 120}, 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 = {-60, -110}, extent = {{-10, -10}, {10, 10}}, rotation = -90), iconTransformation(origin = {-60, -110}, extent = {{-10, -10}, {10, 10}}, rotation = -90)));
Modelica.Blocks.Interfaces.RealOutput Rlive "Output connector of R (Live)" annotation(
Placement(visible = true, transformation(origin = {60, -110}, extent = {{-10, -10}, {10, 10}}, rotation = -90), iconTransformation(origin = {60, -110}, extent = {{-10, -10}, {10, 10}}, rotation = -90)));
Modelica.Blocks.Interfaces.RealOutput Rdead "Output connector of R (Dead)" annotation(
Placement(visible = true, transformation(origin = {0, -110}, 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;
Rdead = 0;
equation
der(Sy) = (-contact_rate * Sy * Iy) + Rlive * reinfect_rate;
der(Iy) = contact_rate * Sy * Iy - Iy * gamma_const;
der(Ry) = Iy * gamma_const - Rlive * reinfect_rate;
Inew = contact_rate * Sy * Iy;
der(Rdead) = Iy * gamma_const * dying_rate;
// der(Rdead) = der(Ry) * dying_rate;
Rlive = Ry - Rdead;
annotation(
Diagram,
Icon(coordinateSystem(initialScale = 0.1), graphics = {Text(origin = {-62, 53}, extent = {{-52, 29}, {174, -133}}, textString = "SIR\nReI")}));
end cl_SIRreinfect;
死亡率、再感染率を加えたシミュレーション結果
死亡率のみを設定してシミュレーション
作成した cl_SIRreinfect を使って、モデル mdl_SIRmodel_reinfect を作成し、シミュレーションします。
ここではワクチンは使用しないものとします。
わかりやすくするために、死亡率をかなり高めの20%に設定してシミュレーションします。
結果、死亡者数は全体の18%くらいに収斂しています。
再感染率のみを設定してシミュレーション
次に、再感染率のみ設定したモデルでシミュレーションしてみます。
再感染率は高めの 1%で試します。
100日までのシミュレーションでは、感染可能者数が増え続けていてこの先がどうなるかが見えません。500日までシミュレーションしてみると、S、I、Rともに振動しながら徐々に収束していきます。
I（感染者数）を見ると、第2波、第3波があるのがわかります。
死亡率、再感染率の両方を設定してみる
上2つのパラメータを同時に設定してみます。
※※（2020/5/13再修正して元に戻しました）
シミュレーションしてみると、再感染しながら死亡者数が徐々に増えつつ 59%くらいが死に至り、I（感染者）、RLive（生存回復者）がゼロに収束していきながら感染可能者数が40%弱に収束しています。
数式としては合っているはずですが、直感的なグラフとは違うような感じがします。
もし数式の立て方に問題があるようでしたら、ご指摘ください。
まとめ
OpenModelica でモデル化することで、基本のSIRモデルに少しずつ要素を加えてシミュレーションすることが容易にできたと思います。
実際はこのシミュレーションより、死亡率、再感染率ともに一桁以上低いと思われますが、より正確な実データが得られれば、パラメータの同定ができるかもしれません。
関連記事
OpenModelicaで感染症の数理モデル「SIRモデル」を実装する
OpenModelicaで感染症の数理モデル「SIRモデル」を実装する（ワクチンの効果を見る）
感染病の数学予測モデルの紹介 (SIRモデル）