はじめに
前回の記事に、SEIRモデルで新型コロナウイルスの挙動を予測するプログラムを掲載しました。
今回は、そのプログラムをGUI化したので、その内容を共有致します。
前回の記事：SEIRモデルで新型コロナウイルスの挙動を予測してみた。
リンク:https://qiita.com/kotai2003/items/ed28fb723a335a873061
実行画面
入力パラメータ一覧
現在、新型コロナウイルスの発症事例より、SEIRのパラメータを推定する研究論文が多数発表されています。今回は、2月16日に発表された論文に掲載されたパラメータ推定値で、SEIRモデルを計算してみます。（参考文献 2)
| Parameter | 中国本土（湖北省除く） | 湖北省（武漢除く） | 武漢 |
|---|---|---|---|
| 人口数 N (million) | 1340 | 45 | 14 |
| 感染率 [beta] | 1.0 | 1.0 | 1.0 |
| Latency period (days) | 2 | 2 | 2 |
| infectious_period (days) | 6.6 | 7.2 | 7.4 |
| E_0 | 696 | 592 | 318 |
| I_0 | 652 | 515 | 389 |
実行例
例えば、社会距離戦略により、感染率が0.5から0.4に下がった場合、感染者のピークがどう変動するかシミュレーションで確認することが可能です。
Case 1: 感染率 0.5
Case 2: 感染率 0.4
感染者（Infections）のピークが下がり、そのタイミングが右側に移動しています。
このような計算により、政府の新型コロナウイルス感染症対策本部が2月23日に発表した「対策の目的（基本的な考え方）」の効果を確認することが可能です。
ソースコード
main_routine.py
import tkinter as tk
from tkinter import ttk
from tkinter import Menu
from tkinter import messagebox
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
from calcSEIR import SEIR_EQ
class Application(tk.Frame):
def __init__(self,master):
super().__init__(master)
self.pack()
self.master.geometry("1000x600")
self.master.title("SEIR Epidemic Model Simulator")
self.create_widgets()
def create_widgets(self):
#Canvas Frame
self.canvas_frame = tk.Frame(self)
self.canvas_frame.configure(width=600, height=480)
self.canvas_frame.grid(row=0, column=0)
self.canvas_frame.grid(padx = 20, pady=20)
#Label Frame for Input Parameters
self.frame_param = tk.LabelFrame( self )
self.frame_param.configure( text=' Input Paramaters ' )
self.frame_param.grid( row=0, column=1 )
self.frame_param.grid( padx=20, pady=20 )
#1. Population
#Label_population
self.label_popu = tk.Label( self.frame_param)
self.label_popu.configure(text ='Population (Million)')
self.label_popu.grid(row =0, column = 0)
#Scale population
self.var_popu = tk.DoubleVar() #scale variable
self.scale_popu = tk.Scale( self.frame_param)
self.scale_popu.configure(orient="horizontal")
self.scale_popu.configure(from_=1, to= 1350)
self.scale_popu.configure(variable=self.var_popu)
self.scale_popu.grid(row=0, column=1)
#2. Infection Rate
# Label_Infection_Rate
self.label_IR = tk.Label( self.frame_param )
self.label_IR.configure( text='Infection Rate' )
self.label_IR.grid( row=1, column=0 )
# Scale Infection_Rate
self.var_IR = tk.DoubleVar() # scale variable
self.scale_IR = tk.Scale( self.frame_param )
self.scale_IR.configure( orient="horizontal" )
self.scale_IR.configure( from_=0.1, to=2 , resolution=0.1)
self.scale_IR.configure( variable=self.var_IR )
self.scale_IR.grid( row=1, column=1 )
#3. Latency Period
# Label_
self.label_LP = tk.Label( self.frame_param )
self.label_LP.configure( text='Latency Period (days)' )
self.label_LP.grid( row=2, column=0 )
# Scale
self.var_LP = tk.IntVar() # scale variable
self.scale_LP = tk.Scale( self.frame_param )
self.scale_LP.configure( orient="horizontal" )
self.scale_LP.configure( from_=1, to=14 , resolution=0.1)
self.scale_LP.configure( variable=self.var_LP )
self.scale_LP.grid( row=2, column=1 )
# 3.5 Infection Period
# Label_
self.label_IP = tk.Label( self.frame_param )
self.label_IP.configure( text='Infections Period (days)' )
self.label_IP.grid( row=3, column=0 )
# Scale
self.var_IP = tk.IntVar() # scale variable
self.scale_IP = tk.Scale( self.frame_param )
self.scale_IP.configure( orient="horizontal" )
self.scale_IP.configure( from_=1, to=14, resolution=0.1 )
self.scale_IP.configure( variable=self.var_IP )
self.scale_IP.grid( row=3, column=1 )
#4 E_0
self.label_E0 = tk.Label( self.frame_param )
self.label_E0.configure( text='E(t=0)' )
self.label_E0.grid( row=4, column=0 )
#Entry
self.Entry_E0 = tk.Entry(self.frame_param)
self.Entry_E0.grid(row=4, column=1)
self.Entry_E0.insert(tk.END,"696")
#5 I_0
self.label_I0 = tk.Label( self.frame_param )
self.label_I0.configure( text='I(t=0)' )
self.label_I0.grid( row=5, column=0 )
# Entry
self.Entry_I0 = tk.Entry( self.frame_param )
self.Entry_I0.grid( row=5, column=1 )
self.Entry_I0.insert( tk.END, "652" )
#6 R_0
self.label_R0 = tk.Label( self.frame_param )
self.label_R0.configure( text='E(t=0)' )
self.label_R0.grid( row=6, column=0 )
# Entry
self.Entry_R0 = tk.Entry( self.frame_param )
self.Entry_R0.grid( row=6, column=1 )
self.Entry_R0.insert( tk.END, "0" )
#7 Time
self.label_time = tk.Label(self.frame_param)
self.label_time.configure( text = 'Time [days]')
self.label_time.grid(row=7, column=0)
self.var_time = tk.IntVar() # scale variable
self.scale_time = tk.Scale( self.frame_param )
self.scale_time.configure( orient="horizontal" )
self.scale_time.configure( from_=10, to=500, resolution=1 )
self.scale_time.configure( variable=self.var_time )
self.scale_time.grid( row=7, column=1 )
#Basic Reproduction Number
# Label Frame result
self.frame_basicR0 = tk.LabelFrame( self )
self.frame_basicR0.configure( text=' Basic Reproduction Number ' )
self.frame_basicR0.grid( row=2, column=1 )
self.frame_basicR0.grid( padx=20, pady=20 )
self.label_basicR0 = tk.Label(self.frame_basicR0)
self.label_basicR0.grid(row = 0, column=0)
self.label_basicR0.configure(text = ' R0 is ')
self.msg_basicR0 = tk.Message(self.frame_basicR0)
self.msg_basicR0.grid(row=0, column=1)
self.msg_basicR0.configure(text ='')
# Button
##Label Frame for Buttons
# Label Frame for Input Parameters
self.frame_button = tk.LabelFrame( self )
self.frame_button.configure( text=' Operation ' )
self.frame_button.grid( row=2, column=0 )
self.frame_param.grid( padx=20, pady=20 )
#button
# Plot (Rungekutta. Plot..Canvas..)
self.button_plot = tk.Button( self.frame_button )
self.button_plot.configure( text='Calculate & Plot' )
self.button_plot.grid( column=0, row=1 )
self.button_plot.configure( command=self.plotCalc )
self.button_plot.configure(width = 20, height=2)
# Quit Button
self.button_quit = tk.Button( self.frame_button )
self.button_quit.config( text='Quit' )
self.button_quit.grid( column=2, row=1 )
self.button_quit.configure( command=self.quit_app )
self.button_quit.configure( width = 15, height=2 )
## Event Call Back
def plotCalc(self):
# parameters
self.t_max = self.var_time.get() # days
self.dt = 0.01
# initial_state
self.N_pop = 1e6*self.var_popu.get()
self.E_0 = int(self.Entry_E0.get())
self.I_0 = int(self.Entry_I0.get())
self.R_0 = int(self.Entry_R0.get())
self.S_0 = self.N_pop - (self.E_0 + self.I_0 + self.R_0)
self.ini_state = [self.S_0, self.E_0, self.I_0, self.R_0] # [S[0],E,[0], I[0], R[0]]
# 感染率
self.beta_const = self.var_IR.get() # 感染率
# 暴露後に感染症を得る率
self.epsilon_const = 1 / self.var_LP.get()
# 回復率や隔離率
self.gamma_const = 1 / self.var_IP.get()
#Basic Reproduction number in SEIR model
self.basicR0 = self.beta_const/self.gamma_const +self.beta_const/self.epsilon_const
self.msg_basicR0.configure( text=str(self.basicR0) )
#https://www.fields.utoronto.ca/programs/scientific/10-11/drugresistance/emergence/fred1.pdf
# numerical integration
self.times = np.arange( 0, self.t_max, self.dt )
self.args = (self.beta_const, self.epsilon_const, self.gamma_const, self.N_pop)
# Numerical Solution using scipy.integrate
# Solver SEIR model
self.result = odeint(SEIR_EQ, self.ini_state, self.times, self.args )
## Plotting
# Generate Figure instance
self.fig = plt.Figure()
#Generate Axe instance
#ax1
self.ax1 = self.fig.add_subplot(111)
self.ax1.plot(self.times, self.result)
self.ax1.set_title('SEIR Epidemic model')
self.ax1.set_xlabel('time [days]')
self.ax1.set_ylabel('population [persons]')
self.ax1.legend(['Susceptible', 'Exposed', 'Infectious', 'Removed'])
self.ax1.grid()
#Link to Axe instance to Canvas
#Then show Canvas onto canvas_Frame
self.canvas = FigureCanvasTkAgg( self.fig, self.canvas_frame )
self.canvas.draw()
self.canvas.get_tk_widget().grid(column=0, row=0)
def quit_app(self):
self.Msgbox = tk.messagebox.askquestion( "Exit Applictaion", "Are you sure?", icon="warning" )
if self.Msgbox == "yes":
self.master.destroy()
else:
tk.messagebox.showinfo( "Return", "you will now return to application screen" )
def main():
root = tk.Tk()
app = Application(master=root)#Inherit
app.mainloop()
if __name__ == "__main__":
main()
calcSEIR.py
# Define differential equation of SEIR model
'''
dS/dt = -beta * S * I / N
dE/dt = beta* S * I / N - epsilon * E
dI/dt = epsilon * E - gamma * I
dR/dt = gamma * I
[v[0], v[1], v[2], v[3]]=[S, E, I, R]
dv[0]/dt = -beta * v[0] * v[2] / N
dv[1]/dt = beta * v[0] * v[2] / N - epsilon * v[1]
dv[2]/dt = epsilon * v[1] - gamma * v[2]
dv[3]/dt = gamma * v[2]
'''
def SEIR_EQ(v, t, beta, epsilon, gamma, N ):
return [-beta * v[0] * v[2] / N ,beta * v[0] * v[2] / N - epsilon * v[1],
epsilon * v[1] - gamma * v[2],gamma * v[2]]