素人が常微分方程式を解いて、コロナの新規感染者数を考えてみた

「常微分方程式って美味しいの?」からスタートするSIRモデルの解法

最近 (2021/11/14現在)、急激な感染者数の減少で若干「勝利」を感じつつある COVID19 の動きを、ド素人の私が今更ながら確認したいと思います。

MATLAB の常微分方程式を解く機能を使ってみて検証

感染症の数理モデルの基本として SIR モデルというものがある。
このあたりの専門家なら基本の"き"のモデルだそうだ。もちろん、COVID19 向けにもう少し手の込んだモデルが提案されているのは言うまでもない。
モデルとしては以下の微分方程式として与えられるらしい。まぁ式を見てみれば、何となくモデルの"気持ち"は分からんでもないが、素人の私にはこれがどんな感じになるのか？見当もつかない。よし、MATLAB で解いてみよう。

各パラメータの意味とか、定義については興味が有る人は、[感染症の数理モデルと対策] (https://www.naika.or.jp/jsim_wp/wp-content/uploads/2020/11/nichinaishi-109-11-article_4.pdf)など読むと良い。
基本再生数

R_0 =\beta \frac{N}{\gamma }

感染率

\begin{array}{l}
\beta =m\times p/N\\
\textrm{where}\;p:\textrm{probability}\;\textrm{per}\;\textrm{contact}\;\textrm{on}\;\textrm{each}\;\textrm{day},\;m:\textrm{number}\;\textrm{of}\;\textrm{contacs}\;\textrm{per}\;\textrm{day}
\end{array}

回復率

\gamma =\frac{1}{\textrm{Days}\;\textrm{for}\;\textrm{recovery}},\;\textrm{(e.g.,}\;\textrm{1/14)}

日本の状況を再現すべく、人口を1億人として、感染症の数理モデルと対策を参考に以下を適当に設定:

• $\displaystyle N=10^8$
• $\displaystyle R_0 =2.5$
• $\gamma^{-1} =10$
• $\beta =R_0 \frac{\gamma }{N}=2.5\times 10^{-9}$

ちなみに

微分方程式を解くと言った時、私の頭に最初に浮かんだのは**学生時代にならったような美しい解 (解析解)**を得る事だったが、そんなものが常に求まるわけでは無いので、プロたちは **"解く= 数値解"**らしい。つまり計算機のパワーを使って、指定した時間期間の中での挙動が分かればそれで良しとするのだ。

知識ゼロからの実装

ツール

• MATLAB 本体のみ
• MATLABによる常微分方程式の解法 のコースを受講 (2-3 hr, 超おススメ)

初期条件 & 設定

常微分方程式を解くので、各パラメータの初期値と、方程式を解く期間 tRange を設定。
タイムステップは "日" にしています。初期の感染者数は I=1000 で設定。
Y0が列ベクトルで与えられているのに注意。

Code

N = 1e8; % Population
Y0 = [N-1000; 1000; 0]; % [S I R] 初期値
tRange = [0, 365]; % 1年間

解く!

ソルバーのことは良く知りませんが、とりあえず ode45 なるソルバーが、汎用的で最初に試すものらしいです。あとはデフォルトでよしなにしてくれます。

Code

[tSol, YSol] = ode45(@MySIR,tRange,Y0);

解説:

• tSol: 解いた後の時間軸
• YSol: 数値解 (tSol に対応)
• @MySIR: 微分方程式の情報を含んだ関数 MySIR の関数ハンドル
  "関数"の中身が気になるところだが、これは後述！とりあえず流れが先。

結果の処理 & 可視化

Code

S = YSol(:,1);
I = YSol(:,2);
R = YSol(:,3);
plot(tSol,S,tSol,I,tSol,R);
legend("Susceptible", "Infected", "Recovered","Location","best");
xlabel("Days");
ylabel("Number of people");

S, I, R の取り出し方がポイント。プロットは以下になります:

新規感染者の感じが本物に似ている！

メトリクス確認

ピーク人数、達する日までの期間などを調べる...

Code

[maxI, idx] = max(I)
maxDay = tSol(idx)

解釈:

• 最大で 2,000万人が1日で新規感染
• 1,000人の感染が確認された時点から、3か月弱がピーク
  おっそろしいウィルスだこと...

アニメーション機能を使ってシミュレーションっぽくしてみる

animatedlineなる関数を使うと、簡単にアニメが作成できる。そして、最近の MATLAB には一度実行したものを再生する機能も付いている!

Code

h = animatedline("Color",'m');
axis([0,365,0,3e7]);
xlabel("Days"); ylabel("# of Infected");
title("COVID 19 using SIR Model");
grid on;

for k=1:length(S)
    addpoints(h,tSol(k),I(k));
    drawnow
end

凄い！凄すぎるぞ Live Editor!
更に凄いのは、右下のアイコンから動画ファイルをそのまま作成できてしまうこと。これは結構最新機能。
このカーブの形状、見覚え有り! 最近 (2021/11/14) もこんな感じで低下中。

ODE(常微分方程式)の関数の中身

さて、本丸...と言っても、数行なのでQiita読者には辛くもないはず。

Code

function dYdt = MySIR(t,Y)
    % Dependant variables
    S = Y(1);
    I = Y(2);
    R = Y(3);

    % Parameters    
    Beta = 2.5e-9;
    Gamma = 0.1;

    % Diferentials
    dSdt = -Beta*S*I;
    dIdt = Beta*S*I - Gamma*I;
    dRdt = Gamma*I;
    
    % Vector of the Diferentials
    dYdt = [dSdt;dIdt;dRdt];
end

解説

• dSdt = -Beta*S*I: $\frac{dS}{dt} = \beta SI$
• dIdt = Beta*S*I - Gamma*I: $\frac{dI}{dt} = \beta SI$
• dRdt = Gamma*I: $\frac{dR}{dt} = \gamma I$
  を作成して、dYdt として列ベクトルで微分計算を返す関数を作成。ここは初期値Y0の与え方に対応!
  引数がt,Y だが、第一引数が独立変数、第二引数が従属変数を渡す。三つ目以降に何かを入れても良いが、使われないらしい。

まとめ

~~個人的には、COVID19 でのシミュレーションを予想ととらえて、"xx教授の予想は全く当たらない！" のような巷の声があり、非常に心を痛めていました。~~シミュレーションは私の感覚では絶対値で語るものではなく、オーダーやトレンドを追いかけるものです。そこで今回は、ウィルスの栄枯盛衰のカーブがどんな形状になるのか、一番基本的なSIRモデルで確認してみました。
素人の私が直ぐに微分方程式を解くことできたのは、MathWorks のオンライントレーニングのお陰です。包括契約がある企業や大学などだと無料で受講できる場合があるので、是非お試しください。