[PYTHON] SIR model that even junior high school students can understand

Recently, I often hear the term SIR model related to coronavirus. It is known that the SIR model fits well in the graph of infected people such as influenza, and it is also useful for the countermeasures against this coronavirus.

I happen to be taking related lessons, so it may be about 100th brew, but I'd like to summarize it easily so that even junior high school students can understand it. Finally, let's simulate with Python code.

What is a SIR model?

The SIR model is a mathematical model that expresses how infectious diseases spread. S stands for uninfected person, I stands for infected person (currently infected person, not the total number of infected people), and R stands for recovered person. In the SIR model, the recoverer is immune and will no longer be infected. It is expressed by the following three formulas without using a small product.

\Delta S = - \beta S I

\Delta I = \beta S I - \nu I

\Delta R = \nu I

If you know chemistry, think of it as similar to the rate equation (law of mass action). $ \ Delta S $, $ \ Delta I $, and $ \ Delta R $ represent changes in $ S $, $ I $, and $ R $, respectively. Specifically, how much $ S $, $ I $, and $ R $ will increase or decrease tomorrow. For example, tomorrow's $ S $ will be $ S + \ Delta S $. $ \ Beta $ is the infection rate and $ \ nu $ is the recovery rate, which is determined by what kind of illness you have. $ 1 / \ nu $ is approximately the number of days it will take to recover.

Now, let's add both sides of these three formulas.

\Delta S + \Delta I + \Delta R = 0

You can see that it is. $ \ Delta S + \ Delta I + \ Delta R $ represents the sum of $ S $, $ I $, and $ R $ that will increase or decrease tomorrow. That's $ 0 $, which means that this model isn't a weird model that will increase or decrease the total population. This time, let's say $ S + I + R = 1 $. For example, if $ I = 0.1 $, then 10% of the inhabitants of that country are $ I $ (infected).

Meaning of SIR model

Next, consider the meaning of the formula. The first formula was as follows.

\Delta S = - \beta S I

This formula shows that tomorrow S (uninfected) will be reduced by $ \ beta S I $. Decreasing means turning into an infected person. The reason why the amount to be reduced is determined by the product of $ S $ and $ I $ can be imagined by substituting a concrete value.

In a country where everyone is infected and there are no uninfected people, $ S = 0 $, so $ SI = 0 . Similarly, in a country where there are no infected people ( I = 0 $), $ SI = 0. It's $. In countries where uninfected and half infected are $ SI = 0.5 \ times 0.5 = 0.25 $, which is the largest. In other words, $ SI $ represents something like the encounter rate between infected and uninfected.

Summarizing the above, we can see that this formula expresses that the product of the encounter rate and the infection rate $ \ beta $ is newly infected. Next, let's look at the third formula.

\Delta R = \nu I

This formula means that as many people as the number of infected people multiplied by the recovery rate will recover. It seems a little strange that the number of recovery depends on the proportion of infected people, but the SIR model assumes this. Finally, let's look at the second formula.

\Delta I = \beta S I - \nu I

Considering the meanings of the first and third formulas, you can see that the number of infected people increases by the amount of (newly infected people)-(newly recovered people). When (newly infected person) <(newly recovered person), this is negative, so we can see that the number of infected people decreases.

Basic reproduction number

I will write it in my free time

Simulation in Python

I will write the explanation in my spare time

sir.py

import matplotlib.pyplot as plt
import numpy as np

beta = 0.5
nu = 0.3
period = 100

def next_state(s, i, r):
    delta_s = - beta*s*i
    delta_i = beta*s*i - nu*i
    s = s + delta_s
    i = i + delta_i
    if s < 0:
        s = 0
    if i > 1:
        i = 1
    r = 1 - s - i
    return s, i, r

def main():
    results = []
    i = 0.01
    r = 0
    s = 1 - i - r

    l_s = [s]
    l_i = [i]
    l_r = [r]
    results.append([s,i,r])
    r0 = beta*s/nu
    
    print(s, i, r)
    for t in range(period):
        s, i, r = next_state(s, i, r)
        results.append([s,i,r])
        l_s.append(s)
        l_i.append(i)
        l_r.append(r)
    print("R0:{}".format(r0))
    plt.figure()
    plt.title("R0:{}".format(format(r0,".3f")))
    plt.xlabel('Time')
    plt.ylabel('Rate')

    x = np.linspace(0, 100, period+1)
    plt.plot(x, l_s, label="S")
    plt.plot(x, l_i, label="I")
    plt.plot(x, l_r, label="R")
    plt.legend()
    plt.show()
    
main()

Execution result