新型コロナウイルスが東京など都市部で感染者数がそれほど多くない理由の仮説

# はじめに
　新型コロナウイルス感染症(COVID-19)の感染者数の推移を[都道府県別新型コロナウイルス感染者マップ](https://gis.jag-japan.com/covid19jp/)等で確認しますと、2020年3月7日現在、北海道で98名、相模原市で21名など局所的にかなりの数に上っている一方、東京都73名、大阪府41名、京都府8名と人の往来の多いはずの都市部で比較的少ない印象です。
　毎日何十万、何百万という人々が通勤・通学あるいは観光し、人口密度の高い場所で接触しながら活動している割に、都市部で感染者数が意外に少ないです。とても不思議ですよね。
　そこで、本記事ではその理由を数理モデルを使って説明できないか、ある仮説を立てて検証してみました。基本的な道具立ては、[先日投稿した記事](https://qiita.com/oki_mebarun/items/60592c703dd4e8d1ff72)をご参照ください。
# 用語
　[前回の記事](https://qiita.com/oki_mebarun/items/60592c703dd4e8d1ff72)でも書きましたが、基本的な用語をおさらいします。

+ S:感染症に対して免疫を持たない者(Susceptible)
+ E:感染症が潜伏期間中の者(Exposed)
+ I:発症者(Infectious)
+ R:感染症から回復し免疫を獲得した者(Recovered)

また、以下のパラメータを使います。

+ R0:一人の感染者が生み出す２次感染者数（基本再生産数）
+ lp:潜伏期間
+ ip:感染期間 

# 仮説
　さて、本記事で提案する仮説ですが、一言でいえば、**ベッドタウンによる浄化効果**です。まずは、下図をご覧ください。
![COVID-Tokyo.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597022/66e38916-aab7-2bb7-e890-96526d64fcc1.jpeg)
　まず、言葉を定義します。

+ S1: 郊外から都市部に通勤・通学する感染症に対して免疫を持たない者
+ S2: 都市部に住む感染症に対して免疫を持たない者
+ E1: 郊外に住む感染症が潜伏期間中の者
+ E2: 都市部に住む感染症が潜伏期間中の者
+ I1: 郊外に住む発症者
+ I2: 都市部に住む発症者
+ R1: 郊外に住む感染症から回復し免疫を獲得した者
+ R2: 都市部に住む感染症から回復し免疫を獲得した者

　図の流れですが、郊外にいるS1の人々が通勤・通学するとともに、都市部にいるS2の人々とともに同じ学校や職場で合流します。そこで、感染者Iから感染し、それぞれE1およびE2となります。E1およびE2からはそれぞれ、I1およびI2に遷移し、最終的にR1およびR2となります。
　ただし、郊外に住むI1のうち、一定比率r $(0 \leq r \leq 1)$ の人々のみが再び都市部に通勤・通学するものとします。残りの比率1-rの人々は自宅待機あるいは入院することで都市部には移動しないものとします。
　この流れから、以下の仮説を立てました。

+ **仮説:** 郊外（ベッドタウン）に住む人々が、感染しても郊外から都市部に再流入する比率rが十分小さければ、かえって、都市部の感染者数のピークが低くなるのではないか？ 

# SEIRモデルの修正
　上記の仮説を[SEIRモデル](https://ja.wikipedia.org/wiki/SEIR%E3%83%A2%E3%83%87%E3%83%AB)を修正して、数式で表してみましょう。

```math
\frac{dS_1}{dt}
 = - \frac{R_0}{ip} (I_2 + r I_1) \frac{S_1}{S_1 + S_2} , S_1 \geq 0
 \\
\frac{dS_2}{dt}
 = - \frac{R_0}{ip} (I_2 + r I_1) \frac{S_2}{S_1 + S_2}, S_2 \geq 0
 \\
\frac{dE_1}{dt}
 = -\frac{dS_1}{dt} - \frac{1}{lp} E_1 \\
\frac{dE_2}{dt}
 = -\frac{dS_2}{dt} - \frac{1}{lp} E_2 \\
\frac{dI_1}{dt}
 = \frac{1}{lp}E_1 - \frac{1}{ip} I_1 \\
\frac{dI_2}{dt}
 = \frac{1}{lp}E_2 - \frac{1}{ip} I_2 \\
\frac{dR_1}{dt}
 = \frac{1}{ip} I_1 \\
\frac{dR_2}{dt}
 = \frac{1}{ip} I_2 
```

最初の２つの式では、感染者（I2 + rI1）が生み出す２次感染者をS1とS2で按分しています。その他の式はSEIRモデルと変わりありません。

# Pythonで計算してみる
上記の仮説をPythonで計算してみましょう。
ライブラリをインポートします。

```Python
import numpy as np
import matplotlib.pyplot as plt
```

R0を計算する関数を定義します。

```Python
def COVID19R0(er):
    if np.random.rand() < er:
        # good environment
        if np.random.rand() < 0.8:
            R0 = 0
        else:
            R0 = np.random.randint(1,4)*1.0
    else:
        # bad environment
        R0 = np.random.randint(0,12)*1.0
    return R0
```

ODEを定義します。

```Python
def seir_eq6(v,t, keys):
    er = keys['er']
    lp = keys['lp']
    ip = keys['ip']
    r  = keys['r']
    #
    R0 = COVID19R0(er)
    #
    s1 = v[0]; s2 = v[1]
    e1 = v[2]; e2 = v[3]
    i1 = v[4]; i2 = v[5]
    r1 = v[6]; r2 = v[7]    
    #
    ds1 = - R0/ip * (i2 + r * i1) * (s1/(s1 + s2)) if s1 >= 0 else 0 # note: s1 >= 0
    ds2 = - R0/ip * (i2 + r * i1) * (s2/(s1 + s2)) if s2 >= 0 else 0 # note: s2 >= 0
    de1 = -ds1 - (1/lp) * e1
    de2 = -ds2 - (1/lp) * e2
    di1 = (1/lp)*e1 - (1/ip)*i1
    di2 = (1/lp)*e2 - (1/ip)*i2
    dr1 = (1/ip)*i1
    dr2 = (1/ip)*i2
    return [ds1, ds2, de1, de2, di1, di2, dr1, dr2]
```

ODEを解くための関数を定義します。

```Python
def my_odeint(deq, ini_state, tseq, keys):
    sim = None
    v = np.array(ini_state).astype(np.float64)
    dt = (tseq[1] - tseq[0])*1.0
    for t in tseq:
        dv = deq(v,t, keys)
        v = v + np.array(dv) * dt
        if sim is None:
            sim = v
        else:
            sim = np.vstack((sim, v))
    return sim
```

最後にシミュレーションと表示用のコードです。

```Python
ini_state=[430, 999, 0, 0, 0, 1, 0, 0]
t_max=180
dt=0.01
t=np.arange(0,t_max,dt)
#
keys = {'er':0.5, 'lp':5.5, 'ip':8, 'r':0.01 }
sim = my_odeint(seir_eq6, ini_state, t, keys)
#
plt.rcParams["font.size"] = 12
fig, ax = plt.subplots(figsize=(10,5))
ax.plot(t,sim[:,[4,5]]) # extract Infected1, Infected2
ax.set_xticks(np.linspace(0,t_max,19))
yw = 100; yn = int(500/yw)+1
ax.set_yticks(np.linspace(0,yw*(yn-1), yn))
ax.grid(which='both')
ax.legend(['Infected1(Outside)','Infected2(in City)'])
ax.set_xlabel('date')
plt.show()
```

init_stateのS1,S2部分を色々と変えて何パターンか計算してみます。比率rは仮に0.01としました。また、その他のパラメータは[前回の記事](https://qiita.com/oki_mebarun/items/60592c703dd4e8d1ff72)と同一です。

# シミュレーション結果

都市部の人口を1000人とし、0日目に1人感染者がいる状態(S2=999,I2=1)からスタートします。 1000人というのは計算しやすさのための仮の値です。

##### まずは、郊外からの流入S1が0の場合です。
![m6_2_s1_0_all.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597022/7eb6cba2-4be1-1d8c-5a7e-019a934721fb.png)

都市部の感染者I2はおよそ62日目でピークが約320名となりました。

##### 郊外からの流入S1が50%の場合です。
S1=500からスタートしました。
![m6_2_s1_50_all.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597022/69e99add-501b-418b-cd73-944d0f5d945c.png)
　都市部の感染者I2はおよそ102日目でピークが約260名となりました。つまり、ピークがおよそ40日遅れ、かつピーク感染者数が約19%減っています。
　ちなみに、[東京都の報道資料](https://www.metro.tokyo.lg.jp/tosei/hodohappyo/press/2018/03/20/12.html)によれば、平成27年の国勢調査の結果、

+ 東京都への流入人口は2,906,056人
+ 東京都の昼間就業者数8,006,342人、昼間通学者1,679,335人（合計9,685,677人）

とありますので、S1とS2の比率は430:1000が一つの目安です。

##### 郊外からの流入S1が80%の場合です。
S1=800からスタートしました。
![m6_2_s1_80_all.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/597022/470fe95e-1e11-5739-144a-26cb17d7653c.png)
　都市部の感染者I2はおよそ142日目でピークが約210名となりました。つまり、ピークがおよそ80日遅れ、かつピーク感染者数が約34%減っています。
　80%という数字は多いようですが、[こちらのページ](http://www.sangyo-rodo.metro.tokyo.jp/toukei/tourism/h30-jittai/)によれば、「平成30年に東京都を訪れた外国人旅行者数は約1,424万人、日本人旅行者数は約5億3,650万人となった。」らしいので、単純計算で、毎日150万人は流入することになります。

# 考察
以上から、郊外のベッドタウンから都市部に人口が流入しつつ、感染者が郊外に留まる（自宅待機・郊外で入院）ことによる感染のピーク変動に関して、次の傾向がシミュレーションから導けます。

+ 郊外からの流入が都市部人口の50%の場合、ピークがおよそ40日遅れ、かつピーク感染者数が約19%減った。
+ 郊外からの流入が都市部人口の80%の場合、ピークがおよそ80日遅れ、かつピーク感染者数が約34%減った。

したがって、**ベッドタウンによる浄化効果**はそれなりにあると思われます。

# さらに言えば・・・

+ 中国の武漢市やイタリアのロンバルディア州のように、都市封鎖をした場合、都市からの感染者流出は防げるでしょうが、**ベッドタウンによる浄化効果**がなくなることで却って、ピーク感染者数を短期に増大させ、医療崩壊を引き起こす恐れがあるかもしれません…。
+ 日本でも、北海道や相模原市など、特定の地域で急速に感染が増大していますが、他地域からの流入が少ないことが要因の一つかもしれません。その場合、地域内での人の交流が密になるのではと懸念されます。
+ 本記事で検証した効果の要因のもう一つの側面は、新規感染者がS1とS2で按分される点です。つまり、S1=0だと閉じた系でS2とI2が密に接さざる得ないため、より感染の機会が増えると思われます。その場合、外出を出来るだけ控える対策が必要でしょう。

# 参考リンク
下記のページを参考にさせて頂きました
[都道府県別新型コロナウイルス感染者マップ](https://gis.jag-japan.com/covid19jp/)
[東京都の昼間人口(従業地・通学地による人口)の概要](https://www.metro.tokyo.lg.jp/tosei/hodohappyo/press/2018/03/20/12.html)
[平成30年訪都旅行者数等実態調査](http://www.sangyo-rodo.metro.tokyo.jp/toukei/tourism/h30-jittai/)
[SEIRモデル](https://ja.wikipedia.org/wiki/SEIR%E3%83%A2%E3%83%87%E3%83%AB)
[感染症数理モデル事始め　PythonによるSEIRモデルの概要とパラメータ推定入門](https://qiita.com/Student-M/items/4e3e286bf08b7320b665)
[新型コロナウイルス対策としての休業の効果をSEIRモデルで検証](https://qiita.com/oki_mebarun/items/60592c703dd4e8d1ff72)