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. (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.
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.
Then define the variables as a whole.
In addition, the following are introduced as parameters.
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.
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
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.
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.
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 $ **.
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)
Let's take a look at the calculation result.
The distribution of basic reproduction numbers is as follows. The average number of basic reproductions was $ E [R_0] = 2.49 $.
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.
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.
Eventually, approximately 90% will be infected, resulting in herd immunity termination.
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. 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.
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% **.
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.
(Quoted from Paper by G.M. Szabo)
I referred to the following page.
Recommended Posts