##はじめに
現在、COVID-19の感染が拡大している。感染症が人間社会の中においてどのように広がっていく過程やその仕組みを理解することはコロナ禍を生きる私たちにとって重要なものである。感染症の数理モデルとして従来から用いられているSIRモデルやその拡張のSEIRモデルがよく知られているが、現在報告されているCOVID-19の特徴はこのモデルで再現できていない部分がある。本稿では、SEIRモデルを基礎として、COVID-19の現在報告されている特徴を挙げ、それらを考慮したより精緻な感染モデルを紹介する。

##SEIRモデル
まず、本稿で提案するモデルの基礎となるSEIRモデルについて簡単に説明する。モデルの名称は、以下に示すモデルの４つの変数の名前の頭文字をとったものである。

　　　図１　SEIRモデルの状態遷移図

SEIRモデルは図1の遷移図に示すように、四つの状態を順に人が遷移していく過程を常微分方程式でモデル化したものである[1]。

##COVID-19の特徴
次に本稿で新たに考慮する三つの特徴を述べる。ただし、COVID-19の感染拡大からまだ日も浅く、現時点で報告されていることをもとにしているため、必ずしもCOVID-19の特徴として正しいとは言えないことを留意してほしい。
　１つ目が、人によって発症した際の症状の重さが異なる点である。正確な割合は不明だが無症状または軽症である人が多く存在する。一方で重症化する人もおり、死に至る可能性も小さくない。この特徴により、無症状の者が無自覚で感染を拡大させる可能性があり、この点をSEIRモデルでは考慮されていない。
　２つ目は、まだ確実な情報ではないが免疫の減退である。抗体が感染し回復した後数ヶ月で失われている事例が報告されている[2]。これもまた免疫を獲得後、半永久的にそれが保持されていると仮定しているSEIRモデルでは表すことができない。
　３つ目は、接触率の時間変化である。感染拡大が始まると、人々が社会的距離をとるようになり、SEIRモデルが仮定している通常の人間社会とは違った状況が発生する。しかし、社会的距離をとるという施策も人々の意識次第であり、その強度は変化する。ある調査を見ても感染者数と自粛率の間には相関がある[3]。この変化の影響は小さくないはずであり、SEIRモデルに新たに組み込むべき特徴といえる。

##拡張SEIRモデル
　前章で挙げた３つの特徴をSEIRモデルに組み込むことを考える。まず、モデル内の人間の状態遷移図を描いてみると図のようになる。図2に示した通り、既存のSEIRモデルとの違いは、発症者を症状の重さによって分類している点、感染症による死者が出る点、及び免疫獲得者が一定期間ののち免疫を失う点である。

　　　図2　拡張SEIRモデルの状態遷移図

　これをもとに常微分方程式を立式すると、以下のようになる。ここでは、発症者のうち軽症・無症状者も感染拡大に影響することを考慮している。

a : 基本接触率（基本感染率）
b : 免疫消失率
c : 発症率
d : 回復率
e : 回復率
f : 死亡率
C : 接触率の時間変化係数

dS/dt=-aCS(E+ I_m )+bR　　　
dE/dt=aCS (E+ I_m )-cE
(dI_m)/dt=mcE- dI_m
(dI_s)/dt=scE-eI_s-fI_s
dR/dt=dI_m+eI_s-bR
dD/dt=fI_s

さらに、前章でも述べた接触率の変化をモデルに組み込むためにパラメータである接触率も時間変化するように上記の常微分方程式に以下の式も加える。この式では、接触率が発症者の数が増えれば増えるほど減少し、回復者が増加するほど増加する様子を表している。しかし、両者の人数はモデル内で規模が異なり、パラメータのgとhの設定が重要である。

dC/dt= -g(I_m+ I_s+D)+ hR

##結果と考察

　　　図3　実行結果

　　　図4　自粛率の推移

下の図4は前述した他者との接触率の変化を表している。

抗体が一定期間で消えること、また死者が存在することをモデルに組み込むことにより感染者数が振動する様子を確認できた。これによりそれぞれの人々が感染者数を参考に自らの行動を抑制、促進した場合の感染者数の推移について確認できた。
図4を見ても自粛率が時間とともに変化している様子を見ることが出来る。

ここでCについて簡単に説明する。
Cは上述したとおり「接触率の時間変化係数」を表すものである。
このCは軽症者数、重傷者数、死者数に-gを掛けたものと、抗体獲得者数にhを掛けたものである。単位は％。

そしてこのg,hは以下のようになっている。

dC/dt= -g(I_m+ I_s+D)+ hR（再掲）

g = 5/S0*100  
h = 1/S0*100

5と1はそれぞれの人数比を調整するためのものである。S0は人数の初期値を表していて、これで「軽症者数、重傷者数、死者数」や「抗体獲得者数」を割ることにより、全体の中でそれらの人がどの程度いるのか割合を示すようにしている。それを100で割ることで100分率にしている。

しかし、この5,1は人数比を調整するために恣意的に決めた値である。
この恣意的にバランスを取る点を改善するために数式を変更し、パラメータhをなくしたものが以下になる。

dC/dt= -g(cE - dI_m -eI_s - fI_s)

g = 5/(S0^0.5)/100
＊S0は人口の初期値

この上の式に置いて、
・ cE は発症者数の増加を
・ -dI_m は軽症者数の減少を
・ -eI_s - fI_s は重傷者数の減少を
表している。
これによりCは感染者が増加していればマイナスの値をとり、減少していればプラスの値を取る。

実行結果は以下の通り。

　　　図5　パラメータ変更後の実行結果

　　　図6　パラメータ変更後の自粛率

図5を見ると時間経過とともに人との接触量8割程度で収束していることが確認できる。

これにより軽症者、重傷者、死者の絶対数ではなく単位時間前からの増加量、つまり何人いるのかではなく、どれだけ増えたかによって人の行動が変わる状況を示した。

この結果を見ると自粛率振動の様子が見えなくなっている。
図3と同様に振動すると予想したが、結果としてはより早い段階で収束した。これは自粛率の変化を発症者数のみで表現したことによるものと考えられ、図らずも最適な自粛が行われた結果（つまり、急速な感染拡大を抑えウイルスと共存する状態）を見ることが出来た。

感染者数の変化だけを理由にし、意思決定をした場合には早い段階で収束することが確認できた。

このモデルの問題点としてgの式を(S0^0.5)としたが、この0.5は恣意的に決めた値であるということがる。この0.5は新規に感染する人と症状が亡くなる人の和は、全体の人数の0.5乗程度であるという感覚にで決められた値である。この値をいかにして変数にするかという観点が課題になる。

##おわりに
今回は感染者数の変化を元に自粛率を変化させる方法を模索し、その方法として微分方程式を用いたSEIRモデルの拡張として自粛率を変数としたモデルを提案した。しかし、今回のモデルにはパラメータの設定等に改善の余地が残されている。今後もモデルの更新を行っていく必要があるだろう。

最後に感染拡大の数理モデルは、私たちの社会において感染症がどのように流行するのかという仕組みやその概観を理解することを手助けしてくれる。しかしそれはあくまで仮定の上に成り立っているモデルであり、実社会をそっくりそのまま再現できているわけではない。それを踏まえた上でどのように実社会をより良くしていくかが大事である。

##使用したコード
参考までに今回使用したコードについても記載する。

import pandas as pd
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
%matplotlib inline
import seaborn as sns
sns.set
import math

def equation(v,t, a,b,c,d,e,f,g,h,m,s):
    dSdt  = b*v[4] - a * v[6] * v[0] * (v[1] + v[2])
    dEdt  = a * v[6] * v[0] * (v[1] + v[2]) - c*v[1]
    dImdt = m * c * v[1] - d * v[2]
    dIsdt = s * c * v[1] - e * v[3] - f * v[3]
    dRdt  = d * v[2] + e * v[3] - b * v[4]
    dddt  = f * v[3] 
    dadt  = -g * (c*v[1] - d * v[2] - e * v[3] - f * v[3])   #h * (v[4]) - g * ( v[2] + v[3])    
    return [dSdt,dEdt,dImdt,dIsdt,dRdt,dddt,dadt]

def graphize(v,disp):
    plt.figure(figsize=(20,20))
    # Main graph
    plt.subplot(2,1,1)
    for i ,name in disp:
        plt.plot(v[:,i],label = name)
    plt.plot(v[:,2]+v[:,3],label = "Im + Is")      # if you merge some variable, use this code
    plt.ylabel("Number of people")
    plt.title("Epidemic Model")
    plt.legend()
    # Parameter graph
    plt.subplot(2,1,2)
    plt.plot(v[:,6],label = "C")
    plt.xlabel("Day")  
    plt.title("Self restraint parameter")
    plt.legend()
    
def Table(V0,v):
    # Variable Table
    Variable = pd.DataFrame([V0,[math.floor(np.max(v[:,i])) for i in range(7)],[math.floor(np.min(v[:,i])) for i in range(7)]],
                             columns = ["S","E","Im","Is","R","D","C"],index = ["init","max","min"])
    return Variable
    
def solve(V0,dspan,params):
    t = np.arange(1,dspan,1)
    v  = odeint(equation,V0,t,args=tuple(params))
    return v
    
def main(V0,dspan,params,disp):
    v = solve(V0,dspan,params)
    graphize(v,disp)
    return Table(V0,v)

# information

Ro = 3        # basic repuroduction number
l = 7         # incubation period
i = 10        # recovery period 
r = 90        # antibody duration

# initial value

S0 = 10000    # population
E0 = 10       # expose潜伏期間
Im0 = 0       # minor symptom
Is0 = 0       # sivire symptom
R0 = 0        # recovery (antibody acquisition)
D0 = 0        # dead
C0 = 1       # contact opportunity rate

# parameters

a = Ro/i/S0   # basic infection rate 
b = 1/r       # antibody disappearance rate
c = 1/l       # incidence rate
d = 1/(i-2)   # recover rate (minor) 
e = 1/i       # recover rate (severe)
f = 0.04      # dead rate
g = 5/(S0**0.5)/100  #1/S0/100/2 # self rastraint rate
h = 5/S0/100  # relief rate
m = 0.8       # minor symptom rate
s = 0.2       # severe symptom rate

# simulate configuration

dspan = 365
disp = zip([0,1,4,5],["S","E","R","D"])
V0 = [S0,E0,Im0,Is0,R0,D0,C0]
params = [a,b,c,d,e,f,g,h,m,s]

pd.DataFrame(params,index=["a","b","c","d","e","f","g","h","m","s"],columns=["Parameter Value"]).T

T =  main(V0,dspan,params,disp)
T