[PYTHON] Predict infectious disease epidemics with SIR model

Introduction

Nowadays, it is said on TV and other media that "80% of behavioral regulations are necessary to resolve the new corona pneumonia." It seems that the basis of this 80% number is a classic model equation called the SIR model. I'm curious, but I implemented the SIR model and checked the effect of the self-restraint rate.

• The author has only the knowledge of Wikipedia about the SIR model, so there may be some mistakes. The source code for this article is just for ** an intuitive understanding of how well behavioral regulation can control infections **. Since the model here is very simple, ** I think it is difficult to predict and take measures against the actual number of infected people **.

SIR model

** The SIR model ** is a classic that deterministically describes the short-term epidemic process of [infectious diseases](https://ja.wikipedia.org/wiki/infectious diseases). [Model equation](https://ja.wikipedia.org/w/index.php?title=model equation & action = edit & redlink = 1). The name is for [Model](https://ja.wikipedia.org/wiki/Mathematical model)

--Sensitible holders (** S ** usceptible) --Infected (** I ** negative) --Immune holders (** R ** recovered or quarantine ** R ** emoved org / wiki / removed)))

Named after the acronym>. The prototype model is [W.O.Kelmac](https://ja.wikipedia.org/w/index.php?title=William Ogilvy Kermak & action = edit & redlink = 1) (English version //en.wikipedia.org/wiki/William_Ogilvy_Kermack)) and [AG McKendrick](https://ja.wikipedia.org/w/index.php?title=Anderson Gray McKendrick & action = edit & redlink = 1) (English version) proposed in a 1927 paper [1 ] / wiki / SIR model # cite_note-1).

However, * β *> 0 represents the infection rate and * γ *> 0 represents the recovery rate (quarantine rate) (the reciprocal 1 / * γ * represents the average infection period).

Exhibitor: [wikipedia]([https://ja.wikipedia.org/wiki/SIR%E3%83%A2%E3%83%87%E3%83%AB](https://ja.wikipedia.org/ wiki / SIR model)))

The SIR model is a three-dimensional simultaneous differential equation of S, R, I, $ S(t)+I(t)+R(t)=const. $ Therefore, it can be replaced with a two-variable equation. Furthermore, assuming * S >> I * and setting γ = 1 and R = βS, it can be replaced with the equation of I only. Here, the following equation is used as the governing equation. $ \frac{dI}{dt}=I(t)(R-1) $ Here, R is a value called [basic reproduction number](https://ja.wikipedia.org/wiki/basic reproduction number), and it is said that COVID-19 is about 2 to 3.

Literature → https://www.biorxiv.org/content/10.1101/2020.01.25.919787v1

If you write down the above formula as a forward difference, $ I(t+Δt)=I(t)+I(t)(R-1)Δt $ It will be. In the following source code, we will simulate based on this formula.

Source code

import

Use matplotlib for graph visualization.

from matplotlib import pyplot

Parameters

R = 2.5
activity = 0.2
x0 = 200
infection_days = 14
dt = 1/infection_days
stop_day = 30
start_day = 150
total_day = 180

Variable initialization

x = [x0]
new = [0]

Variable update

for day in range(1,total_day):
	p = 1
	if day >= stop_day:
		p = activity
	if day >= start_day:
		p = 1
	x.append(x[day-1]+x[day-1]*(R*p-1.0)*dt)
	new.append(x[day-1]*R*p*dt)

p represents the activity rate, which is 1 until day becomes stop_day, ʻactivity from stop_daytostart_day`, and 1 after that. Multiplying this value by R represents the basic reproduction number at the time of self-restraint. For example, if p = 0.5, then R × p = 2.5 × 0.5 = 1.25.

Graph display

pyplot.plot(x,label='Infected persons')
pyplot.plot(new,label='New infected persons')
pyplot.xlabel('days')
pyplot.ylabel('Persons')
pyplot.legend()
pyplot.show()

Transition of infected persons when behavior is reduced by 80% (ʻactivity = 0.2`)

The number of infected people decreased rapidly from the self-restraint start date (stop_day), but the activity start date (start_day) was reached before it was completely eradicated, and the number of infected people gradually increased after that date. I have. You can see that even if the number of infected people decreases, it is usually useless.

Transition of infected persons when behavior is reduced by 60% (ʻactivity = 0.4`)

Since the graph becomes hard to see, I set start_day = 180 (no reactivation). It's a graph that you often see on TV. If ʻactivity = 0.4`, then R × p = 1, which means that ** an average of 1 person is infected per infected person **. Under the premise of R = 2.5, it is necessary to maintain the status quo with a 60% reduction and to reduce the number of infected people by 80% in consideration of uncertainty.

at the end

Again, the source code for this article is just for ** an intuitive understanding of how behavioral regulation can control infections **. I think the model is too simple and tight for quantitative evaluation.

APPENDIX code full

from matplotlib import pyplot

#Basic reproduction number (number of newly infected persons per infected person during normal activities)
R = 2.5
#Behavior rate during the self-restraint period (1 is normal)
activity = 0.2
#Number of initial infections
x0 = 200
#Days to quarantine or recovery
infection_days = 14
#Computational time step
dt = 1/infection_days
#Self-restraint start date (activity rate changes from 1 to activity after this date)
stop_day = 30
#Day to end self-restraint and resume normal activities
start_day = 150
#Total simulation days
total_day = 180

#Current number of infected people
x = [x0]
#Number of newly infected people
new = [0]

#Value update (difference method)
for day in range(1,total_day):
	p = 1
	if day >= stop_day:
		p = activity
	if day >= start_day:
		p = 1
	x.append(x[day-1]+x[day-1]*(R*p-1.0)*dt)
	new.append(x[day-1]*R*p*dt)
	
#graph display
pyplot.plot(x,label='Infected persons')
pyplot.plot(new,label='New infected persons')
pyplot.xlabel('days')
pyplot.ylabel('Persons')
pyplot.legend()
pyplot.show()