新型コロナウイルスの感染者数予測～簡単なモデルから～

今中国、そして世界各国で猛威をふるいはじめている新型コロナウイルスの感染者数を、シンプルな微分方程式系でモデリングしてみたいと思います。伝染病関係の数理モデリングをするのは今回が初めてで、色々おかしい部分があると思いますので予めご容赦ください。

注)素人が実施したこの感染者数予測の信頼性は皆無です。

なお、文献としては以下を参考にしています。

微分方程式で数学モデルを作ろう

数理モデル

以下の状況を仮定します。

• 新型コロナウイルスは中国の湖北省（人口5902万人）で発生
• 人から人へのみ伝染する
• 湖北省から人の出入りはない(渡航禁止)
• 当局により感染が発覚した人間は即刻除外される（除外された人からは伝染しない）
• 当局が発見できる感染者は、全体の感染者のうちの一部である（割合γ）

また、次のような視点からモデルを組み立てます。

『感染者I』と『未感染者S』が接触することにより感染するので、『未感染者の減少速度』は『感染者I』と『未感染者S』の数の積に比例する。またこの比例定数を感染速度rとする。

$$ \frac{dS}{dt}=-rSI $$

しかし、『感染者I』のうちある一定割合（γ）は除外されるので、その分『感染者I』は減少する。

$$ \frac{dI}{dt}=rSI - \gamma I $$

上述の通り、『感染者I』のうちある一定割合（γ）は『除外者R』となる

$$ \frac{dR}{dt}=\gamma I $$

上記3つの式を連立させることで、微分方程式系を解くことができます。

パラメーターフィッティング

上記微分方程式系には2つのパラメーターが存在します。

• 全体の感染者数に対する、当局が発見できる感染者の割合γ
• 感染速度r

今回は簡略化のため、γ=0.05を仮定します。申し訳ないのですがこの値に特に根拠はありません。ただ、感染を申告すると隔離されてしまう状況だとすると隠す人が多いのではないでしょうか。

あとは、実データに合うように感染速度rを決定します。Wikipediaから湖北省の感染者推移データが入手できます。

「2019年-2020年中国武漢における肺炎の流行」の版間の差分 - Wikipedia

感染者数の推移

なお、今回のモデル上では上記記事にあるような『当局が発見した感染者数』＝『除外者数』であることに注意してください。

除外者の時系列データをまとめると以下のようになります。

day	R（除外者数）
0	1
23	27
26	44
28	59
34	41
35	41
36	41
37	41
38	41
39	41
40	45
41	62
42	198
43	198
44	270
45	375
46	444
47	549
48	729
49	1052
50	1413
51	2714

これについて、上記微分方程式系をMatlabのode45ソルバーで解き、実データのプロットに対してフィッティングしたグラフは以下のようになります。

○が実データ、実線がode45による計算値です。おおよそ実データと合っています。 このとき、r=3.65E-9でした。

感染者数の予測

パラメーターが確定できたので、感染者数の予測を行います。r=3.65E-9, γ=0.05で300日間のシミュレーションを行ったときの感染者数のグラフは以下のようになります。

上記グラフによると、今後感染者数が急激に増大し、120日～130日後（凡そ4月～5月頃）にピークを迎え、2500万人が感染していることになります（湖北省の人口の半分です）。

対策の検討

政府による諸対策はこのモデル上においては以下のような効果を生み出します。

• 病床数を増やす→γを増加させる
• 感染者を検知できる仕組みづくり→γを増加させる
• マスクをつける→rを減少させる

例えば、γを増加させることの数理モデル上の影響を確認してみます。

上記グラフの青がγ=0.05、赤がγ=0.15です。これだけでも感染者数の抑制効果が大きいことがわかります（r=3.65E-9のままです)。

今回はここまでにします。 最後にコードを載せて終わります。

微分方程式系(odefcn.m)

function dydt = odefcn(t,y,r,gamma)
dydt = zeros(3,1);
dydt(1) = -r*y(1)*y(2);
dydt(2) = r*y(1)*y(2) - gamma*y(2);
dydt(3) = gamma*y(2);

%% 実データ
tdata = [0 23 26 28 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51];
Rdata = [1 27 44 59 41 41 41 41 41 41 45 62 198 198 270 375 444 549 729 1052 1413 2714];

r = 3.65E-9;
gamma = 0.05;

%% ode45による数値計算
tspan = [0 300];
y0 = [59020000 1 0];
[t1,y1] = ode45(@(t,y) odefcn(t,y,r,gamma), tspan, y0);

%% 実データと計算値の比較
plot(t1,y1(:,3))
hold on
xlabel('day')
ylabel('除外者(R)')
plot(tdata,Rdata,'o')
xlim([0.0 60.0])
ylim([0 10000])
hold off

%% 感染者数のグラフ
plot(t1,y1(:,2))
xlabel('day')
ylabel('感染者数(I)')
hold off


%%　gammaの値を増加
r = 3.65E-9;
gamma = 0.15;

%% ode45による数値計算
tspan = [0 300];
y0 = [59020000 1 0];
[t2,y2] = ode45(@(t,y) odefcn(t,y,r,gamma), tspan, y0);

%% gammaの影響をプロットで比較
plot(t1,y1(:,2),t2,y2(:,2))
xlabel('day')
ylabel('感染者数(I)')