はじめに

最近はコロナコロナでなんだか世紀末みたいになってますね。ふとPythonでこのコロナの感染拡大について分析できないかなあ、と思いたち本稿を執筆しました。あくまでツールとして分析手法を提示するだけの記事のため、この記事が何かを結論付けることはありませんのでご了承ください。この記事を呼んで感染症分析やグラフ理論について興味を持っていただけたら幸いです。是非本稿でPythonのnetworkxを動かしてみて下さい☆

目次

• SIRモデル
• 接触ネットワーク
• Watts-Strogatzモデル

SIRモデル

SIRとはSusceptible Infectious Removedの略で、S（健康な人）・I（病気を移す可能性のある感染者）・R（病気を移す可能性のない感染者）のように人を３つのグループに分けて感染動態を表現する方法です。まず実際に式を見てみましょう。

S + I → 2I

I → R

これがSIRモデルとなります。化学反応式と同じイメージをして下さい。健康な人Sは感染者Iと会うことである一定の確率で感染者Iになります。これが一番目の式です。また、感染者Iは完治したり死亡したり隔離されるなどの様々な要因によって病気を移す可能性のない感染者Rとなります。これが二番目の式です。（賢い方はお気づきでしょうが再感染を考慮していません）

次にここに時間発展の概念を与えます。時間tにおけるSをS(t)、IをI(t)、RをR(t)とします。上の式ですが、化学反応と同様にSはSとIの積に比例して定率でIに移行し、Iは定率でRに移行するという仮定となっております。つまり、これを連続力学系で表すと

\frac{dS}{dt}(t)=-βS(t)I(t)

\frac{dI}{dt}(t)=βS(t)I(t)-γI(t)

\frac{dR}{dt}(t)=γI(t)

となります。ここでβは感染率、γは回復や隔離をまとめて表現した定数となります。以上がSIRモデルの説明となります。ここまでの内容をグーグルコラボで実装しました。このコードで実行したのが下の結果です。定数を変えていじってみて下さい。

#モジュールのインポート
import numpy as np
from scipy.integrate import odeint
from scipy.optimize import minimize
import matplotlib.pyplot as plt
%matplotlib inline

#式の定義V[0]はSでV[1]はIが入る
def SIR_EQ(v, t, beta, gamma):
    return [-beta*v[0]*v[1], beta * v[0] * v[1] - gamma * v[1], gamma * v[1]]

#tの最大値とtの微小変化値
t_max = 120
dt = 0.01
#βγの値を決める
beta_const = 0.001
gamma_const = 0.45
#SIRの初期値
S_0=499
I_0=1
R_0=0
ini_state = [S_0,I_0,R_0] #[S[0], I[0], R[0]]
#格納
times =np.arange(0,t_max, dt)
args  =(beta_const, gamma_const)

#SIRmodelを解く
result = odeint(SIR_EQ, ini_state, times, args)
#plot
plt.plot(times,result)
plt.legend(['Susceptible','Infectious', 'Removed'])

例　500人の中で１人の感染者からβ=0.001、γ=0.1で感染拡大

例　500人の中で１人の感染者からβ=0.001、γ=0.45で感染収束
Rの最終的な値が総感染者数と考えると、保因者の隔離の重要性がわかりますね。

接触ネットワーク

もう一度SIRモデルの仮定を考えて見ましょう

• S；感染する可能性がある人。基礎疾患の有無などは無視
• I；人に感染させるポテンシャルのある人。潜伏期などは無視
• R；Iから移行する感染することも感染させることもない人。回復・隔離・死亡などが含まれる
• SとIは化学反応のようにランダム接触する

以上がSIRモデルの仮定となります。しかし、実社会では地理的にも人間関係的にも接触はランダムではありません。この仮定には無理があります。ここで接触ネットワークの考え方を導入します。

まず人をノード(まる○)、接触をエッジ(線ー）と表すことにします。下の図を見て下さい。

人と人の接触をこのようなグラフで表し、感染アルゴリズムを考えていきます。
まず、この中からランダムに人(ノード)を一つ選択します。今回真ん中の人(ノード)を選択したとします。
$ $選ばれた人(ノード)がSの場合；このノードが確率 $βnIdt$ でIに変えます。このnは隣接するIの数で今回の例ではn=2となります。 $dt$ については時間ステップをかけていると考えて下さい。
$ $選ばれた人(ノード)がIの場合；このノードが確立 $γdt$ でRに変えます。
選ばれた人(ノード)がRの場合；何もしない。
以上のアルゴリズムを繰り返します。このようにネットワークにおける感染症をシミュレーションします。（追記；ノード一つあたりから伸びるエッジの本数を平均次数と言います。本稿ではkと表します）

Watts-Strogatzモデル

では、どのようなネットワークの上でこのアルゴリズムを走らせれば良いのでしょうか？これはWatts-Strogatzモデルを使うのが良いとされています。あるノードから最近傍のK(平均次数)個のノードと結び付けます。ノード数n=20、平均次数k=8の場合下の図の様になります。

ここから確率pでネットワークをランダムに組み替えます。p=0.01,0.1,1で組み替えてみましょう。

↓p=0.01

↓p=0.1

↓p=1

p=0が規則正しいレギュラーなネットワークで、p=1が完全にランダムなネットワークとなります。では、これらのL(平均ノード距離)、C(クラスター係数)を求めてみましょう。L(平均ノード距離)とは何人たどれば任意の人に会えるかみたいなイメージです（日本だと友達7人たどれば有名人に会えるなんていいますよね）。C(クラスター係数)は友達の友達が友達である確率のようなものです。自分とエッジで結ばれているノードとエッジで結ばれているノードが自分と直接エッジで繋がっている確率というものがC(クラスター係数)です。上の例だと次のようになります。

ノード数n=20なのであまり違いがわかりませんね。では、n=1000でやってみましょう。結果は下のようになります。

確率pが大きくなるに連れてLとCが減少することが分かりますね。面白いことにp=0.1,0.01ではCの値をある程度の大きさに保ったままLが大きく減少していますね。これらをSmall Worldと呼んだりします。次にこれらのネットワークからSIRモデルを考えてみましょう。以下にこの章で使ったコードを記載します。

#google colabにてインストール
!pip install networkx matplotlib
!pip install -q ndlib
!pip install -q bokeh

#モジュールのインポート
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
import ndlib.models.epidemics as ep
import ndlib.models.ModelConfig as mc
from bokeh.io import output_notebook, show
from ndlib.viz.bokeh.DiffusionTrend import DiffusionTrend

#watts_strogatzモデルの定義(ノード数,平均次数,組み換え確率)
G=nx.watts_strogatz_graph(1000,8,1)
nx.draw(G) #nx.draw_random(G)
plt.show()
print(nx.info(G))
#平均ノード間距離
print(nx.average_shortest_path_length(G))
#クラスター係数
print(nx.average_clustering(G))

SIR + Watts-Strogatzモデル

それでは、最後にSIRモデルとWatts-Strogatzモデルを実際に合わせてみましょう。下の様に実装します。

#google colabにてインストール
!pip install networkx matplotlib
!pip install -q ndlib
!pip install -q bokeh

#モジュールのインポート
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
import ndlib.models.epidemics as ep
import ndlib.models.ModelConfig as mc
from bokeh.io import output_notebook, show
from ndlib.viz.bokeh.DiffusionTrend import DiffusionTrend

#watts_strogatzモデルの定義(ノード数,平均次数,組み換え確率)
G=nx.watts_strogatz_graph(1000,8,0.01)
nx.draw(G) #nx.draw_random(G)
plt.show()
print(nx.info(G))
#平均ノード間距離
print(nx.average_shortest_path_length(G))
#クラスター係数
print(nx.average_clustering(G))

model = ep.SIRModel(G)
config = mc.Configuration()
config.add_model_parameter('beta',0.1)
config.add_model_parameter('gamma',0.1)
config.add_model_parameter('percentage_infected',0.001)
model.set_initial_status(config)

iterations=model.iteration_bunch(75)
trends=model.build_trends(iterations)
output_notebook()

viz = DiffusionTrend(model,trends)
p = viz.plot(width=400,height=400)
show(p)

ここでのpercentage_infectedというのは集団での最初の感染者の割合です。ここでは1000人の集団で１人の感染者からの増加を考えたいので、percentage_infected=0.001として計算します。$β=0.1$,$γ=0.1$,$n=1000,k=8$で計算してみましょう。上からp=0,0.01,0.1,1での出力結果となります。
↓p=0

↓p=0.01

↓p=0.1

↓p=1

ランダム性が増すにつれて感染爆発のタイミングが早く急激に起こることが改めて分かりますね。若者のライブやクラブなどはランダム性の増加につながるので辞めましょうね（こじつけ）

皆さんもこのコードのパラメータをいじって分析してみてください。この分析手法で医学的な観点で分析した記事についてはまた今度出すかもしれません、、

参考資料

2010年、日本人工知能学会 湯浅友幸・白山晋、感染症流行予測におけるネットワーク構造の影響分析
山下雄史のyoutube研究室SIRモデルとネットワーク理論