　　　　On a Thread of the Web: 緊急事態宣言と接触制限モデル（１）

北海道大学の西浦博さんが4月3日にマスコミで示していた図がある。対策がない状態の新規感染者数が指数関数的な増大をしているときに，人の接触を20%減らした場合と80%減らした場合の新規感染者数の変化を示した図だ。これにより彼は欧米に近い外出制限の必要性を訴えていた。残念ながら，これは専門家会議や政府の共通のコンセンサスとはならず，やがて4月7日の7都府県に対する5月6日までの緊急事態宣言につながっていく。

ここではその図がどのようにして得られたものかを，素人が持っている簡単な道具と知識で理解することを目指す。報道されたニュースや伝わってくる情報は鵜呑みにはせず，自分で考えることが必要だから。集団としては人口1400万人の東京都を想定する。日本経済新聞におけるこの内容に関する報道の文脈でも東京を対象としていることがうかがえる。東京の4月3日における新規感染者数は14人（新規感染数累計は753人）であり，人口の1ppm条件がちょうど達成された頃なので，これ以降をSIRモデルなどの単純な微分方程式系で扱うことは可能だろう。

今，探してもみあたらないのだが，西浦さんはこの報道の後に解説のための動画を公開していた。それによればしばらく前のドイツなみの再生産数を仮定し，現時点から15日目以降に制限措置を行った場合をシミュレートしているとした。なお，彼のモデルでは，計算開始から15日目が現在であり，30日目が制限措置の開始時点となる。

今から3週間前のドイツでは新規感染数累計の7日増倍率が7倍になっていた。仮に6.3倍（1日に1.3倍）を仮定すると，基本再生産数は$R_0 = 1 + 5/7 * \log 6.3 \approx 2.3$となる*。一方現在の東京では新規感染数累計の7日増倍率が2.5倍（1日に1.14倍）程度なので，基本再生産数は$R_0 = 1 + 5/7 * \log 2.5 \approx 1.65$（倍加時間と基本再生産率を参照）である。だからこのドイツ並の仮定はちょっとどうかと思うが，いちおうこれを採用する。
（*西浦さんはドイツなどの欧州の平均基本再生産数が$R_0=2.5$であると話していた）

つまり，現在100人の新規感染者数はこのまま放置すると15日後のt=30（1.3^15≒50）に5000人以上に達するという状況を西浦さんは想定していると思われる。今の東京の水準をそのままあてはめた場合は，15日後のt=30（1.14^15≒7）に東京の新規感染者数は700人程度であることに注意しよう。

SIRモデルに接触制限措置の効果を含めた微分方程式系をMathematicaで書くと次のようになる。nは集団人口（1400，以下人数の単位は1万人当りとなる），βは感染率（0.44），αは感染期間（5），νは感染数の初期値（0.0002），pは制限措置の程度（p=10→制限なし，p=9→80%に制限，p=6→20%に制限）である。なお，接触制限措置は，βにかけるシグモイド関数（幅が1日のオーダー）によってモデル化して階段関数は用いない。t=15が4月3日現在（新規感染者数=100）に対応する。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
f[n_,β_,α_,ν_,p_]:=
NDSolve[{S'[t]==-(p+(10-p)*Tanh[2*(30-t)])/10*β/n*J[t]*S[t],
J'[t]==(p+(10-p)*Tanh[2*(30-t)])/10*β/n*J[t]*S[t]-J[t]/α, R'[t]==J[t]/α, S[0]==n,J[0]==ν,R[0]==0}, {S,J,R}, {t,0,100}]

β=0.44 sol = f[1400, β, 7, 0.0002, 10];
s[t_] := S[t] /. sol[[1, 1]];
i[t_] := J[t] /. sol[[1, 2]];
r[t_] := R[t] /. sol[[1, 3]];
cf1[t_] = (10+0*Tanh[2*(30-t)])/10*β*i[t]*s[t]/1400;

sol = f[1400, β, 7, 0.0002, 9];
s[t_] := S[t] /. sol[[1, 1]];
i[t_] := J[t] /. sol[[1, 2]];
r[t_] := R[t] /. sol[[1, 3]];
cf2[t_] = (9+1*Tanh[2*(30-t)])/10*β*i[t]*s[t]/1400;

sol = f[1400, β, 7, 0.0002, 6];
s[t_] := S[t] /. sol[[1, 1]];
i[t_] := J[t] /. sol[[1, 2]];
r[t_] := R[t] /. sol[[1, 3]];
cf3[t_] = (6+4*Tanh[2*(30-t)])/10*β*i[t]*s[t]/1400;

Plot[{cf1[t], cf2[t], cf3[t]}, {t, 15, 45},
PlotRange -> {0, 2}, GridLines ->
{{29, 30, 31}, {0.1, 0.4, 0.5, 0.6, 0.7, 0.8}}]
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

計算の結果，次のグラフが得られた。西浦さんのものと概ね一致する。