【シミュレーション入門】コロナ感染をシミュレートして遊んでみた♬

数値計算の分野では、昔から自由運動する粒子のシミュレーションやブラウン運動する粒子群の集団運動などのシミュレーションがある。
今回は、これをコロナウィルス感染のシミュレーションに応用してみようと思う。

似たようなシミュレーションとしてPSOというのがあるが今回はこれらのpythonコードを参考にした。
【参考】
・粒子群最適化(PSO)のパラメータについ
・粒子群最適化法(PSO)を救いた
###やったこと
・コード解説
・感染率依存性
・粒子密度依存性
・粒子運動依存性
###・コード解説
コード全体は以下においた。
・collective_particles/snow.py
コードの主要な部分を説明する。
利用ライブラリは以下のとおり

import numpy as np
import matplotlib.pyplot as plt
import random
import time

初期値の定義は以下のとおり

PARTICLE_NO = 1000 # 粒子数
ITERATION = 200 # 最大ループ回数　感染者が０になると止まる
MIN_X, MIN_Y = -100.0, -100.0 # 探索開始時の範囲最小値
MAX_X, MAX_Y = 100.0, 100.0 # 探索開始時の範囲最大値
recovery=30 #一定時間経過したら治癒
p=0.03 #probability of infecion

グラフ表示は、上記の結果を見るとわかるように、
ax1:粒子の位置情報
ax2:青；正常数、赤；感染数、緑；治癒数の頻度グラフ
を以下のように表示している。
引数は、グラフ番号、位置情報、経過時間、r,g,bの値

def plot_particle(sk,positions,elt,r,g,b):
    el_time = time.time()-start
    fig, (ax1, ax2) = plt.subplots(2, 1, sharey=False,figsize=(8, 16))
        
    for j in range(0,PARTICLE_NO):
        x=positions[j]["x"]
        y=positions[j]["y"]
        c=positions[j]["c"]
        s = 5**2  #粒子サイズ
        ax1.scatter(x, y, s, c, marker="o")
    ax1.set_xlim([MIN_X, MAX_X])
    ax1.set_ylim([MIN_Y, MAX_Y])
    ax1.set_xlabel("x")
    ax1.set_ylabel("y")
    ax1.set_title("{:.2f}:InfectionRate;{:.2f} %".format(el_time,(PARTICLE_NO-b[-1])/PARTICLE_NO*100)) #累計感染率
    ind = np.arange(len(elt))  # the x locations for the groups
    width = 0.3       # the width of the bars
    ax2.set_ylim([0, PARTICLE_NO])
    ax2.set_title("{:.2f}:red_{} green_{} blue_{}".format(el_time,r[-1],g[-1],b[-1]))
    rect1 = ax2.bar(ind, b,width, color="b")
    rect2 = ax2.bar(ind+width, g, width, color="g")
    rect3 = ax2.bar(ind+2*width, r,width, color="r")
    plt.pause(0.1)
    plt.savefig('./fig/fig{}_.png'.format(sk)) 
    plt.close()

今回の主なロジックは、粒子の位置情報更新関数に押し込めた。
すなわち、粒子オブジェクトの属性を以下のとおり定義している。

属性	値
位置情報；	(x,y)
感染属性；	(blue,red,green)
感染した時間（初期値からの経過）；	t_time
感染履歴flag；	s（感染１、未感染０）

position.append({"x": new_x, "y": new_y, "c": p_color, "t": t_time,"flag":s})

上記の変数の時間変化を以下の関数で求めます。
感染判定は以下の式で行っています。
半径２０の円の中に入ると接触したものとして、一定確率pで感染するかどうかを評価します。つまり、この半径と感染確率pが感染の強さを決めています。

if (x-x0[k])**2+(y-y0[k])**2 < 400 and random.uniform(0,1)<p:

上記のように粒子オブジェクトを定義したので、それぞれの変数を取り出したり代入しています。
もう一つ、簡単のために治癒は一定時間経過すると自動的に治癒するようにしました。ここも確率で治癒させてもいいのですが、ここでは事象を複雑化するだけなので、詳細化は不要です。
また、粒子運動（人の移動）は、平均的にはランダムウォークとしました。人はそんなに動かないという近似です。一方、分子運動のシミュレーションみたいに自由運動させることもできますが、やっていません。
そして、感染は治癒の人は感染させないという仮定です。

# 粒子の位置更新関数
def update_position(positions):
    x0 = []
    y0 = []
    for i in range(PARTICLE_NO):
        c=positions[i]["c"]
        t_time = positions[i]["t"]  #初期値０，感染は感染時時間
        k_time = time.time()-start  #経過時間
        s = positions[i]["flag"]    #感染なし０，感染：１
        if s == 1 and c == "red":   #感染済な場合
            if k_time-t_time>recovery:  #一定時間経過したら治癒
                c = "blue"
                positions[i]["c"] = "green"
                positions[i]["flag"] = 1   #ただし、感染履歴ありのまま
        if c == "red":  #感染redなら位置情報取得
            x0.append(positions[i]["x"])
            y0.append(positions[i]["y"])
    position = []
    for j in range(PARTICLE_NO):
        x=positions[j]["x"]
        y=positions[j]["y"]
        c=positions[j]["c"]
        s = positions[j]["flag"]
        t_time = positions[j]["t"]
        for k in range(len(x0)):
            if (x-x0[k])**2+(y-y0[k])**2 < 400 and random.uniform(0,1)<p:
                if s ==0:
                    c = "red"
                    t_time = time.time()-start
                    s = 1
                    positions[j]["flag"]=s
                else:
                    continue
        vx = 1*random.uniform(-1, 1)
        vy = 1*random.uniform(-1, 1)
        new_x = x + vx
        new_y = y + vy
        p_color = c
        s=s
        
        position.append({"x": new_x, "y": new_y, "c": p_color, "t": t_time,"flag":s})
    return position, x0

main関数は以下のとおり、
ポイントは、初期値で１粒子を感染；red、flag;1などとしている。
その他の粒子は、感染無しで乱数で配置している。最初の粒子の位置を真ん中に配置するというアイデアもあるが、いろいろな配置からの感染伝播が見たいので任意の位置とした。count_brg()関数は単純にカウントしている。

def main():
    # 各粒子の初期位置, 速度, 
    position = []
    velocity = []  #速度は今回使わない
    # 初期位置, 初期速度
    position.append({"x": random.uniform(MIN_X, MAX_X), "y": random.uniform(MIN_Y, MAX_Y), "c": "red", "t":0, "flag":1})
    for s in range(1,PARTICLE_NO):
        position.append({"x": random.uniform(MIN_X, MAX_X), "y": random.uniform(MIN_Y, MAX_Y), "c": "blue", "t": 0, "flag":0})
    sk = 0    
    red=[]
    green=[]
    blue=[]
    elapsed_time = []
    while sk < ITERATION:
        position, x0 = update_position(position)
        r,g,b = count_brg(position)
        red.append(r)
        green.append(g)
        blue.append(b)
        el_time=time.time()-start
        elapsed_time.append(el_time)
        plot_particle(sk,position,elapsed_time,red,green,blue)
        if x0==[]:
            break
        sk += 1

###・感染率依存性
まず、感染率pを変化させると感染伝播の様子や全体としての累計感染率はどう変化するだろうか。
上記の例は感染率p=30%のものであり、累計感染率=100%であり、確実に感染伝播が発生しており、津波のように伝播しているのがわかる。また、感染ピークは404、67.5secであった。この場合、感染伝播速度のような群速度が定義できそうであるが、なんとなくやめました。
一方、感染率p=5%では、累計感染率=92.20%まで下がり、感染伝播の様子も以下のように青が目立つようになり、感染ピーク235と抑えられ、ピーク位置も１７８secまで伸びた。

さらに、感染伝播が発生するかどうかの臨界的な感染率=3%の場合、以下のとおりとなった。すなわち、津波のような感染伝播は消え、恐る恐る伝播するイメージである。感染ピークも117と低くなり、260secと長くなった。累計感染率=52%まで減少した。なお、この感染率の場合、この計算結果も途中で終了しているが、ほとんど感染せずに初期の段階で消滅することもあった。すなわち、手洗いやマスクなどをうまく使って感染率を一定以下（今回のシミュレーションだと３％以下）に下げられればそれだけで感染伝播をなくすことも可能だという意味です。

これらの結果は、ある程度はマスクやハグやキスや咳エチケットなどのコミュニケーションの方法をコントロールすることにより、長引くかもしれないが累計感染率を下げることが出来ることを示していると考察できる。そして、このシミュレーションをよく見ると、一度治癒した人がいる領域にある青い粒子（未感染者）は感染から守られており、かつ感染伝播が内向きに伝播することはない。すなわち、この領域の集団は、感染に関して堅牢であり滅多なことでは感染するリスクを回避できる状態である。すなわち、集団免疫を獲得したと言える。この集団免疫は簡単に言うと治癒者が集団の未感染者の密度を実質的に下げるため感染伝播を防ぐと言い換えることもできる。
###・粒子密度依存性
もう一つの興味は、通常イベント中止とか集会中止などを起こしている、集まらないという常識が正解なのかどうかということである。
ここでは感染確率３０％、リカバリ時間３０secの条件下で、粒子密度を変化して振る舞いの違いを見てみた。
結果は、感染率の粒子密度依存性を表に表すと以下のようになった。

粒子密度	30	40	60	80	100	120	140	160	180	200
累計感染率％	3.33	27.5	8.33	40	33	49.17	67.86	70	91.67	96.5
すなわち６０以下の密度では感染伝播は発生せず、８０／２００Ｘ２００辺りから感染伝播が発生し始めるため累計感染率が大きくなり始め、２００辺りでほぼ１００％程度に到達する。
例として、以下に遷移領域である１４０のときのシミュレーションを示す。
特徴は、上記の密度が十分大きいときは単一ピークであったが、この程度の密度だと感染者のピークはよりなだらかであり、粒子分布の密度ゆらぎを反映して、いくつかのピークを示すことである。そして、感染伝播の様子は密度が小さい領域を乗り越えるために、苦労しており、時には間が空きすぎている場合には感染伝播ができないことが見える。
すなわち、集まらないという方針は全体の密度を下げることにより累計感染率を下げる効果があるばかりでなく、局所的にも集まらないという方針で行動することにより、感染の確率は下げられるということである。
###・粒子運動依存性
次に、出かけないというのは正しいかどうかを検証しよう。
この効果は一応、粒子運動が大きいと感染伝播や累計感染率にどういう影響があるか見た。
粒子密度120のとき、結果は以下の表の通りとなった。

粒子運動	0	1	2	3	4	8	16
累計感染率％	10	49.17	54.17	83.33	79.17	65.83	74.17
この評価から、静止している場合のこの密度ではほぼ感染はしない。一方、粒子運動があると累計感染率は増加し、激しいほど累計感染率は上がる傾向にあることがわかる。これは単純に上記で記載したように感染伝播を阻む隙間を、運動により乗り越えていると想像できる。
実際のシミュレーションの代表的な例は以下のとおりである。
このシミュレーションを見ると想像通り、少し大きな空間も乗り越えて感染伝播しているのがわかる。
つまり、出かけないというのも感染リスクを下げる行為であることがわかる。
この運動する粒子モデルについてもさらに現実の世界を模したいろいろなシチュエーションも考えてシミュレートできるがきりがないので、今回はここまでとする。
###まとめ
・コロナ感染をシミュレーションして遊んでみた
・人が集まらず、手洗いやマスクなどをうまく使って感染率を一定以下（今回のシミュレーションだと３％以下）に下げられればそれだけで感染伝播を発生させないことも可能という結論を得た
・集まらない、マスクや手洗いなどにより感染率を下げると、感染の終息には時間がかかるが、累計感染率を下げられることがわかった。
・出かけないの効果検証として、運動する場合は累計感染率が増加することがわかった
・高密度な場合、感染伝播は水滴を落とした波紋の広がりに似ている
・低密度だが感染伝播が存在する場合、治癒者と未感染者が共存する領域が広がって来るが、この領域は感染に対して堅牢であり、集団免疫な状態と言える

・このシミュレーションはもっと大規模で行うと、実際の街を模したものへの拡張も可能であるので、もう少し火力を上げて計算してみようと思う
ちなみに、今回はJetson-nanoを利用したが、1000までは問題なく計算できた。