Solving mathematical models of infectious disease epidemics in Python

**Caution! This article does not contain more information than I have tried. ** **

Infectious disease mathematical model

In response to the recent COVID 19 situation, I think that various mathematical models of infectious diseases have been talked about, but let's solve it with Python odeint and visualize it with matplotlib as it is. It sounds awkward, but you shouldn't express any feelings about COVID 19 from the result of properly implementing something in this article.

So, let's say that the easiest SIR model to solve is

\begin{align}
\frac{\mathrm{d} S}{\mathrm{d} t} & = - \beta I S\\
\frac{\mathrm{d} I}{\mathrm{d} t} & = \beta I S - \gamma I\\
\frac{\mathrm{d} R}{\mathrm{d} t} & = \gamma I
\end{align}

I will leave it as. $ \ Beta $ is the infection rate (/ (person x hours)), the number of people who have not yet been infected per unit time per sick person, and $ \ gamma $ is the cure rate ( / Hour), the reciprocal of the time it takes to heal. $ S $ is the number of people who have not been infected and have no infectivity, $ I $ is the number of people who are infected and have the ability to infect others, $ R $ is quarantine, cure, death, etc. The number of people who have been infected with or have lost the ability to infect others.

Python code Since it is troublesome, I put a solid code without any explanation.

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider, Button

def solve(value, t, foi, rr):
	# fio = force of infection (beta)
	# rr = recovery rate
	S, I, R = value
	dSdt = - foi * I * S
	dIdt = foi * I * S - rr * I
	dRdt = rr * I
	return [dSdt, dIdt, dRdt]


span = 100
t = np.linspace(0, span, 256)
initial_value = [0.9, 0.1, 0] # initial values for S, I, R
solution = odeint(solve, initial_value, t, args = (0.5, 0.5))

fig, ax = plt.subplots(figsize=(12, 6))
plt.xlim(0, span)
plt.ylim(0, 1)
plt.subplots_adjust(bottom = 0.2)

data_S, = plt.plot(t, solution[:, 0], "blue", label = "Susceptible")
data_I, = plt.plot(t, solution[:, 1], "red", label = "Infectious")
data_R, = plt.plot(t, solution[:, 2], "orange", label = "Removed")

plt.legend(loc = "best")

rr_slider = Slider(plt.axes([0.2, 0.02, 0.65, 0.02]), "recovery rate", 0, 1, valinit = 0.5, valstep = 0.1)
foi_slider = Slider(plt.axes([0.2, 0.05, 0.65, 0.02]), "force of inf.", 0, 1, valinit = 0.5, valstep = 0.1)

def update(val):
	rr = rr_slider.val
	foi = foi_slider.val
	solution = odeint(solve, initial_value, t, args=(foi, rr))
	data_S.set_ydata(solution[:, 0])
	data_I.set_ydata(solution[:, 1])
	data_R.set_ydata(solution[:, 2])
	fig.canvas.draw_idle()

rr_slider.on_changed(update)
foi_slider.on_changed(update)

plt.show()

Execution result

so

Well play it. If you want to solve the differential equation while undulating the parameters, please do not hesitate to contact us.