[PYTHON] The theory that the key to controlling infection with the new coronavirus is hyperdispersion of susceptibility

Introduction

The SEIR model or SIR model is often used in the prediction and analysis of the spread of the new coronavirus infection (COVID-19), but recently based on the scale-free network [Paper by GM Szabo](https://arxiv. org / pdf / 2004.00067.pdf) and Paper by Y. Ohsawa have been published and are attracting attention. To put it simply, a scale-free network means that "some vertices are connected to many other vertices by edges and have a large degree, while most of the others have a few vertices. It is only connected to and has a small degree. ”([Complex network](https://ja.wikipedia.org/wiki/%E8%A4%87%E9%9B%91%E3%83%8D%] It is a network structure with E3% 83% 83% E3% 83% 88% E3% 83% AF% E3% 83% BC% E3% 82% AF)), but this property was found in a cluster group survey. ** The basic reproduction number of a small number of people is large, while the basic reproduction number of the majority of people is small ** ([New Corona Cluster Countermeasure Expert (2)]](https: // togetter. It is very similar to com / li / 1492002)) and seems to have a high affinity. scale-free-network.jpg (Quoted from Paper by G.M. Szabo) So, in this article, inspired by this study, ** how does the overall spread of infection affect when individual basic reproduction numbers are determined according to a wide-tailed distribution? I would like to extend ** to the SEIR model and verify it.

Definition of terms

Perform a simulation based on the SEIR model. First, define groups and individuals.

Next, we define a random variable for individual j.

Since it is a probability, it naturally satisfies the following.

  1. \forall j, 0 \leq s_j, e_j, i_j, r_j \leq 1
  2. \forall j, s_j + e_j + i_j + r_j = 1

Then define the variables as a whole.

In addition, the following are introduced as parameters.

Basic reproduction number and overdispersion

The basic reproduction number $ R_0 $ is defined as "the number of secondary infections produced by one infected person (in the early stages of infection) during the infection period". The reason for the early stage of infection is that the value fluctuates due to changes in the SEIR ratio over time. If it depends on time, it is written as $ R_t $. According to the Survey by the Cluster Countermeasures Group of the Ministry of Health, Labor and Welfare, the basic reproduction number is sampled as shown in the graph below. クラスター対策班のR0_2.jpg

The parameter $ e_r $ above shows the probability of being in a poorly ventilated environment on this graph. As can be seen from this graph, ** the basic reproduction number of a small number of people is large, while the basic reproduction number of the majority of people is small **. This distribution of basic reproduction numbers is realized by the following approximation code.

def COVID19R0(er):
    if np.random.rand() < er:
        # good environment
        if np.random.rand() < 0.8:
            R0 = 0
        else:
            R0 = np.random.randint(1,4)*1.0
    else:
        # bad environment
        R0 = np.random.randint(0,12)*1.0
    return R0

It should also be noted that ** basic reproduction numbers, by definition, have implications for both the time dimension and the number of infected persons **. For example

  1. One person continued to infect one person a day for eight days.
  2. One person infected eight people a day and did not subsequently develop a secondary infection.

In both cases, $ R_0 = 8 $, but 1 implies ** the strength of the infectivity of the individual **, and 2 is rather ** people are 3 dense, etc. It can be said that it implies the sensitivity ** of being in the scene. Now, let's set up a mathematical model that calculates the SEIR model based on the distribution of basic reproduction numbers. Consider the following two types of models.

  1. A model that assumes overdispersion of infectivity. ** Consider overdispersion as a factor on the part of the affected person **.
  2. A model that assumes hyperdispersion of susceptibility. ** Consider overdispersion as a factor on the part of the recipient **.

Mathematical model 1 (overdispersion of infectivity)

A model that assumes overdispersion of infectivity can be written as follows.

\begin{eqnarray}
\beta_j &=& \frac{R_{0j}}{ip N} \\
\frac{ds_j}{dt} &=& - s_j \cdot \sum_{k \neq j} \beta_k i_k \\
\frac{de_j}{dt} &=& -\frac{ds_j}{dt} - \frac{1}{lp} e_j \\
\frac{di_j}{dt} &=& \frac{1}{lp} e_j - \frac{1}{ip} i_j \\
\frac{dr_j}{dt} &=& \frac{1}{ip} i_j 
\end{eqnarray}

Especially the second line is important, but if you sum this formula with respect to $ j $,

\begin{eqnarray}
\frac{dS}{dt} &=& - S \cdot BI  + \sum_j s_j i_j \beta_j \\
\frac{dS}{dt} &=& \sum_j \frac{ds_j}{dt} \\
BI &=& \sum_k \beta_k i_k 
\end{eqnarray}

It will be. If $ \ forall k, \ beta_k = \ beta $,

\begin{eqnarray}
\frac{dS}{dt} &=& - \beta S \cdot I  + \sum_j s_j i_j \beta_j \\
I &=& \sum_j i_j \\
\end{eqnarray}

Therefore, only the second term is different from the normal SEIR model, but this term indicates that ** I will not make myself a secondary infected person **, considering the transition of the SEIR model. It can be ignored because it can be said that it is almost $ s_j i_j = 0 $. In other words, it is a model ** that results in a normal SEIR model when ** $ \ beta_j $ is homogenized.

Mathematical model 2 (hyperdispersion of sensitivity)

A model that assumes hyperdispersion of susceptibility can be written as follows.

\begin{eqnarray}
\beta_j &=& \frac{R_{0j}}{ip N} \\
\frac{ds_j}{dt} &=& - s_j \beta_j \cdot \sum_{k \neq j}  i_k \\
\frac{de_j}{dt} &=& -\frac{ds_j}{dt} - \frac{1}{lp} e_j \\
\frac{di_j}{dt} &=& \frac{1}{lp} e_j - \frac{1}{ip} i_j \\
\frac{dr_j}{dt} &=& \frac{1}{ip} i_j 
\end{eqnarray}

Focusing on the second line, you can see that $ \ beta_k $ has moved to $ \ beta_j $ for the mathematical model 1 above. As before, this model is a model that can be reduced to a normal SEIR model by equalizing ** $ \ beta_j $ **.

Try to calculate with Python

Now, let's calculate mathematical model 1 and mathematical model 2 using Python. First, the premise of calculation is as follows.

Consider a group of + $ N = 200 $.

Import the library.

import numpy as np
import matplotlib.pyplot as plt

A function that calculates the basic reproduction number and generates a distribution. This time, I chose $ er = 0.6 $ so that the average value is around $ E [R_0] = 2.5 $. Also, $ R_ {0j} $ is sorted in descending order for $ j $.

def COVID19R0(er):
    if np.random.rand() < er:
        # good environment
        if np.random.rand() < 0.8:
            R0 = 0
        else:
            R0 = np.random.randint(1,4)*1.0
    else:
        # bad environment
        R0 = np.random.randint(0,12)*1.0
    return R0

N  = 200
er = 0.6
R0 = np.array([COVID19R0(er) for i in range(N)])
R0 = np.sort(R0)[::-1]
plt.hist(R0,bins=10)
plt.title("Ro distribution with E[Ro]: {}".format(np.average(R0)))
plt.xlabel('R0')
plt.ylabel('num. of people')
plt.show()

A function that solves ordinary differential equations.

def my_odeint(deq, ini_state, tseq, keys):
    sim = None
    v = np.array(ini_state).astype(np.float64)
    dt = (tseq[1] - tseq[0])*1.0
    c = 0
    for t in tseq:
        dv = deq(v,t, keys)
        v = v + np.array(dv) * dt
        if sim is None:
            sim = v
        else:
            sim = np.vstack((sim, v))
        pg = 100.* t / tseq[-1]
        if pg >= c:
            c = c + 10
            print("progress: {}%".format(pg))
    return sim

This is a calculation function of mathematical model 1.

#define differencial equation of seir model
def seir_eq10_1(v, t, keys):
    N  = keys['N']
    b  = keys['b']
    lp = keys['lp']
    ip = keys['ip']
    #
    ds = np.zeros(N)
    de = np.zeros(N)
    di = np.zeros(N)
    dr = np.zeros(N)
    #
    BI = np.sum([ b[k]*v[2][k] for k in range(N)])
    #
    for j in range(N):
        s  = v[0][j];
        e  = v[1][j];
        i  = v[2][j];
        r  = v[3][j];
        #
        ds[j] = - s * (BI - b[j]*i) #Do not self-infect
        de[j] = - ds[j] -  (1/lp) * e
        di[j] = (1/lp)*e - (1/ip) * i
        dr[j] = (1/ip)*i
    return [ds, de, di, dr]

This is a calculation function of Mathematical Model 2.

#define differencial equation of seir model
def seir_eq10_2(v, t, keys):
    N  = keys['N']
    b  = keys['b']
    lp = keys['lp']
    ip = keys['ip']
    #
    ds = np.zeros(N)
    de = np.zeros(N)
    di = np.zeros(N)
    dr = np.zeros(N)
    #
    CI = np.sum([ v[2][k] for k in range(N)])
    #
    for j in range(N):
        s  = v[0][j];
        e  = v[1][j];
        i  = v[2][j];
        r  = v[3][j];
        #
        ds[j] = - s * b[j] * (CI - i) #Do not self-infect
        de[j] = - ds[j] -  (1/lp) * e
        di[j] = (1/lp)*e - (1/ip) * i
        dr[j] = (1/ip)*i
    return [ds, de, di, dr]

This is a display function.

def showSim_01(n, sim, keys):
    N = keys['N']
    plt.rcParams["font.size"] = 12
    fig, ax = plt.subplots(figsize=(10,5))
    
    x = range(N)
    
    ax.plot(x,sim[n+0]) # extract S
    ax.plot(x,sim[n+1]) # extract E
    ax.plot(x,sim[n+2]) # extract I
    ax.plot(x,sim[n+3]) # extract R
    ax.grid(which='both')
    ax.legend([ 'Susceptible','Exposed', 'Infected', 'Recoverd'])
    ax.set_xlabel('person no.')
    ax.set_ylim(0,1)
    plt.show()

def showSim_02(t, sim, keys):
    N = keys['N']
    plt.rcParams["font.size"] = 12
    fig, ax = plt.subplots(figsize=(10,5))
    
    n = len(t)
    sdat = [np.average(sim[4*i+0]) for i in range(n)]
    edat = [np.average(sim[4*i+1]) for i in range(n)]
    idat = [np.average(sim[4*i+2]) for i in range(n)]
    rdat = [np.average(sim[4*i+3]) for i in range(n)]

    ax.plot(t,sdat) # extract S
    ax.plot(t,edat) # extract E
    ax.plot(t,idat) # extract I
    ax.plot(t,rdat) # extract R
    ax.grid(which='both')
    ax.legend([ 'Susceptible','Exposed', 'Infected', 'Recoverd'])
    ax.set_xlabel('date')
    ax.set_ylim(0,)
    plt.show()

This is a simulation function. The initial state is uniform for each individual.

def calcsim( R0, keys):
    # solve seir model
    N = keys['N']
    #           S  　　　　E 　　　　　I 　　　　　　R
    ini_state=[np.ones(N), np.zeros(N), np.zeros(N), np.zeros(N)]
    ini_state=[np.ones(N)*0.99, np.zeros(N), np.ones(N)*0.01, np.zeros(N)]
    t_max=180
    dt=0.01
    t=np.arange(0,t_max,dt)
    keys['b'] = R0/(N*keys['ip'])
    #
    sim = my_odeint(seir_eq10, ini_state, t, keys)
    #    
    return sim,t

This is the code to execute and display the mathematical model 1. It's very heavy (-_-;).

seir_eq10 = lambda v, t, keys: seir_eq10_1(v, t, keys)
keys = {'N':N, 'lp':5, 'ip':8, 'R0':2.5 }
sim,t = calcsim(R0, keys)
showSim_01(np.int(180/0.01*0.95), sim, keys)
showSim_02(t, sim, keys)

This is the code to execute and display the mathematical model 2.

seir_eq10 = lambda v, t, keys: seir_eq10_2(v, t, keys)
keys = {'N':N, 'lp':5, 'ip':8, 'R0':2.5 }
sim,t = calcsim(R0, keys)
showSim_01(np.int(180/0.01*0.95), sim, keys)
showSim_02(t, sim, keys)

simulation result

Let's take a look at the calculation result.

Distribution of basic reproduction numbers

The distribution of basic reproduction numbers is as follows. The average number of basic reproductions was $ E [R_0] = 2.49 $. m10_R0dist.jpg

Mathematical model 1 (overdispersion of infectivity)

First, in Mathematical Model 1, the horizontal axis is each individual, and the vertical axis is the probability of (s, e, i, r) after 95% time has passed. case1_m10_1_N_200_at17100.jpg

How, ** There is almost no individual difference due to overdispersion of infectivity! ** The result is. Next, let's look at the transition of SEIR as a group like a normal SEIR model. case1_m10_2_N_200.jpg

Eventually, approximately 90% will be infected, resulting in herd immunity termination.

Mathematical model 2 (hyperdispersion of sensitivity)

Next, in Mathematical Model 2, the horizontal axis is each individual, and the vertical axis is the probability of (s, e, i, r) after 95% time has passed. case2_m10_1_N_200_at17100.jpg Very different from the previous result, ** the influence of individual R0 distribution is very strong! ** The result is. Individuals with larger $ R_0 $ tend to be more susceptible to infection, and individuals with smaller $ R_0 $ tend to be less susceptible to infection. Next, let's look at the transition of SEIR as a group. case2_m10_2_N_200.jpg

The end result is ** about 37% infected and terminated by herd immunity !! **. Since the result calculated by Mathematical Model 1 is 90%, ** the number of infected people was suppressed by 59% **.

Consideration

From the above, the following trends can be derived from the simulation regarding the calculation by the SEIR model considering the overvariance of the basic reproduction number.

Furthermore ...

initial-distribution.jpg (Quoted from Paper by G.M. Szabo)

Reference link

I referred to the following page.

  1. [Complex Network](https://ja.wikipedia.org/wiki/%E8%A4%87%E9%9B%91%E3%83%8D%E3%83%83%E3%83%88% E3% 83% AF% E3% 83% BC% E3% 82% AF)
  2. Propagation and mitigation of epidemics in ascale-free network
  3. COVID-19 Should be Suppressed by Mixed Constraints -- from Simulations on Constrained Scale-Free Networks
  4. Estimate the peak infectivity of the new coronavirus
  5. The corona turmoil may not end forever due to false positive PCR tests
  6. New Corona Cluster Countermeasure Expert (2)
  7. “Definition of terms for new corona cluster countermeasures” reported by specialists for new corona cluster countermeasures

Recommended Posts

The theory that the key to controlling infection with the new coronavirus is hyperdispersion of susceptibility
The story that the private key is set to 600 with chmod
I tried to predict the behavior of the new coronavirus with the SEIR model.
I tried to automatically send the literature of the new coronavirus to LINE with Python
I tried to visualize the characteristics of new coronavirus infected person information with wordcloud
Quantify the degree of self-restraint required to contain the new coronavirus
Plot the spread of the new coronavirus
To make sure that the specified key is in the specified bucket in Boto 3
How to judge that the cross key is input in Python3
Create a bot that posts the number of people positive for the new coronavirus in Tokyo to Slack
Did the number of store closures increase due to the impact of the new coronavirus?
I tried to streamline the standard role of new employees with Python
[Introduction to Python] What is the method of repeating with the continue statement?
Folding @ Home on Linux Mint to contribute to the analysis of the new coronavirus
Estimate the peak infectivity of the new coronavirus
Python> set> Convert with set ()> dictionary is only key> I was taught how to convert the values of dictionary to set / dir ({}) / help ({}) / help ({} .values)
Let's simulate the transition of infection rate with respect to population density with python
Hackathon's experience that it is most important to understand the feelings of the organizer
I tried to display the infection condition of coronavirus on the heat map of seaborn
A new form of app that works with GitHub: How to make GitHub Apps
The basis of graph theory with matplotlib animation
Factfulness of the new coronavirus seen in Splunk
GUI simulation of the new coronavirus (SEIR model)
Verify the effect of leave as a countermeasure against the new coronavirus with the SEIR model
Here is one of the apps with "artificial intelligence" that I was interested in.
Since the stock market crashed due to the influence of the new coronavirus, I tried to visualize the performance of my investment trust with Python.
Two solutions to the problem that it is hard to see the vector field when writing a vector field with quiver () of matplotlib pyplot
Add information to the bottom of the figure with Matplotlib
Try to get the contents of Word with Golang
Let's test the medical collapse hypothesis of the new coronavirus
Script that changes the length of the sound with REAPER
Use hash to lighten collision detection of about 1000 balls in Python (related to the new coronavirus)
The sound of tic disorder at work is ... I managed to do it with the code
When a character string of a certain series is in the Key of the dictionary, the character string is converted to the Value of the dictionary.
I tried to predict the infection of new pneumonia using the SIR model: ☓ Wuhan edition ○ Hubei edition
A simple version of government statistics (immigration control) that is easy to handle with jupyter
Try to visualize the nutrients of corn flakes that M-1 champion Milkboy said with Python
[VLC] How to deal with the problem that it is not in the foreground during playback
[Python] The status of each prefecture of the new coronavirus is only published in PDF, but I tried to scrape it without downloading it.
If the people of Tokyo become seriously ill with the new coronavirus, they may be taken to a hospital in Kagoshima prefecture.
I want to output while converting the value of the type (e.g. datetime) that is not supported when outputting json with python