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.

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).

{\displaystyle {\begin{aligned}{\frac {dS}{dt}}(t)&=-\beta S(t)I(t)\{\frac {dI}{dt}}(t)&=\beta S(t)I(t)-\gamma I(t)\{\frac {dR}{dt}}(t)&=\gamma I(t)\end{aligned}}}

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
Parameters Description
R Basic reproduction number (number of newly infected persons per infected person during normal activities)
activity Behavior rate during the self-restraint period (1 is normal)
For example, 1 for 80% self-restraint-0.8=0.It will be 2.
x0 Number of initial infections
infection_days Days to quarantine or recovery
dt Computational time step
stop_day Self-restraint start date (activity rate changes from 1 to activity after this date)
start_day Day to end self-restraint and resume normal activities
total_day Total simulation days

Variable initialization

x = [x0]
new = [0]
Variable name Description
x[day] t=Number of infected people on day
new[day] t=Number of new infections on day

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

figure1.png

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

figure2.png

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()

Recommended Posts

Predict infectious disease epidemics with SIR model
Mathematical model of infectious disease epidemics
Implement the mathematical model "SIR model" of infectious diseases with OpenModelica
Predict hot summers with a linear regression model
Solving mathematical models of infectious disease epidemics in Python
Introduction of mathematical prediction model for infectious diseases (SIR model)
COVID-19 simulation with python (SIR model) ~~ with prefectural heat map ~~
Mathematical prediction model for infectious diseases (SIR model): Case study (1)
Model fitting with lmfit
Day 71 I tried to predict how long this self-restraint will continue with the SIR model
Regression with linear model
Implement the mathematical model "SIR model" of infectious diseases with OpenModelica (example of repeating regulation and deregulation)
[Introduction to infectious disease model] I tried fitting and playing ♬