新型コロナウイルスのGUIシミュレーション　(SEIRモデル）

#はじめに

前回の記事に、SEIRモデルで新型コロナウイルスの挙動を予測するプログラムを掲載しました。
今回は、そのプログラムをGUI化したので、その内容を共有致します。

`前回の記事：SEIRモデルで新型コロナウイルスの挙動を予測してみた。`
リンク:https://qiita.com/kotai2003/items/ed28fb723a335a873061

![bandicam-2020-03-20-11-11-35-036.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/d935c6fa-b42b-ad83-e575-8d801b9609ca.gif)



#実行画面
##入力パラメータ一覧

現在、新型コロナウイルスの発症事例より、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  |

##実行例

例えば、[社会距離戦略](https://ja.wikipedia.org/wiki/%E7%A4%BE%E4%BC%9A%E8%B7%9D%E9%9B%A2%E6%88%A6%E7%95%A5)により、感染率が0.5から0.4に下がった場合、感染者のピークがどう変動するかシミュレーションで確認することが可能です。
###Case 1: 感染率　0.5
<img width="740" alt="IR0.5.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/374246d7-3cac-4ff4-64ef-0476d133d312.png">
###Case 2: 感染率 0.4 
感染者（Infections）のピークが下がり、そのタイミングが右側に移動しています。
このような計算により、政府の新型コロナウイルス感染症対策本部が2月23日に発表した「対策の目的（基本的な考え方）」の効果を確認することが可能です。

<img width="746" alt="IR0.4.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/cd829d62-a1c1-9859-faeb-1d0e5c99f094.png">
<img width="512" alt="counter.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/208980/730f36d1-ca4c-6384-f236-2ce34a790611.png">


#ソースコード
## main_routine.py

```Python
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

```Python


# 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]]
```

#参考資料
1. [SEIRモデルで新型コロナウイルスの挙動を予測してみた。](https://qiita.com/kotai2003/items/ed28fb723a335a873061)
2. [Epidemic analysis of COVID-19 in China by dynamical modeling](https://arxiv.org/abs/2002.06563) 
3. [【Python】Tkinterのテンプレート](https://qiita.com/kotai2003/items/fe9b5e59c7164a95ded8)
4. [Tkinterにmatplotlibグラフを埋め込む](https://qiita.com/kotai2003/items/45953b4d037a62b2042c)
5. [感染病の数学予測モデルの紹介　(SIRモデル）](https://qiita.com/kotai2003/items/3078f4095c3e94e5325c)
6. [新型コロナウイルス感染症への対応について](https://www.kantei.go.jp/jp/singi/novel_coronavirus/th_siryou/sidai_r020223.pdf)