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

#「常微分方程式って美味しいの?」からスタートするSIRモデルの解法
最近 (2021/11/14現在)、急激な感染者数の減少で若干「勝利」を感じつつある COVID19 の動きを、**ド素人の私**が今更ながら確認したいと思います。

## MATLAB の常微分方程式を解く機能を使ってみて検証
感染症の数理モデルの基本として SIR モデルというものがある。
このあたりの専門家なら基本の"き"のモデルだそうだ。もちろん、COVID19 向けにもう少し手の込んだモデルが提案されているのは言うまでもない。  
モデルとしては以下の微分方程式として与えられるらしい。まぁ式を見てみれば、何となくモデルの"気持ち"は分からんでもないが、素人の私にはこれがどんな感じになるのか？見当もつかない。よし、MATLAB で解いてみよう。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/473841/e991e179-1a2d-b158-de3f-18238e4290cd.png",width="500">

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

```math
R_0 =\beta \frac{N}{\gamma }
```
**感染率**

```math
\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}
```

**回復率**

```math
\gamma =\frac{1}{\textrm{Days}\;\textrm{for}\;\textrm{recovery}},\;\textrm{(e.g.,}\;\textrm{1/14)}
```
  
日本の状況を再現すべく、人口を1億人として、[感染症の数理モデルと対策](https://www.naika.or.jp/jsim_wp/wp-content/uploads/2020/11/nichinaishi-109-11-article_4.pdf)を参考に以下を**適当に**設定:

   -  $\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による常微分方程式の解法](https://matlabacademy.mathworks.com/jp/details/solving-ordinary-differential-equations-with-matlab/odes) のコースを受講 (2-3 hr, 超おススメ)

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

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

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

```matlab:Code
[tSol, YSol] = ode45(@MySIR,tRange,Y0);
```
**解説**:
- `tSol`: 解いた後の時間軸
- `YSol`: 数値解 (`tSol` に対応)
- `@MySIR`: 微分方程式の情報を含んだ関数 `MySIR` の関数ハンドル
"関数"の中身が気になるところだが、これは後述！とりあえず流れが先。

## 結果の処理 & 可視化

```matlab: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` の取り出し方がポイント。プロットは以下になります:
 ![result.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/473841/aa031e67-5155-d3c6-8056-a9caff2d7964.jpeg)

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

###メトリクス確認
ピーク人数、達する日までの期間などを調べる...

```matlab:Code
[maxI, idx] = max(I)
maxDay = tSol(idx)
```

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

## アニメーション機能を使ってシミュレーションっぽくしてみる
`animatedline`なる関数を使うと、簡単にアニメが作成できる。そして、最近の MATLAB には一度実行したものを再生する機能も付いている!

```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
```

![2021-11-14_20-51-48.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/473841/df2c3a60-f271-924b-c722-830f3a5b446e.gif)

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

# ODE(常微分方程式)の関数の中身
さて、本丸...と言っても、数行なのでQiita読者には辛くもないはず。

```matlab: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 のオンライントレーニングのお陰です。包括契約がある企業や大学などだと無料で受講できる場合があるので、是非お試しください。