[PYTHON] Let's explain the difference in how the new coronavirus spreads between Japan and other countries with a cluster countermeasure model

Introduction

The new coronavirus infection (COVID-19) is spreading from China to Europe and the Americas and is rampant. In Japan as well, as of March 27, there are 1,254 infected people in Japan, and the governor of Tokyo is requesting restrictions on weekend outings in Tokyo as well. However, compared to the G7 countries where medical care seems to be relatively well-developed, the number of people in Japan is small, with 85594 in the United States, 80589 in Italy, 43938 in Germany, 9155 in France, 11658 in the United Kingdom, and 4043 in Canada. It stands out. In this article, as one hypothesis to explain this reason, we examined whether it could be explained by ** a mathematical model focusing on cluster countermeasures **. The flow of examination is the problems of the conventional model, the examination of the new model, the simulation, and the consideration.

Conclusion

First of all, from the conclusion,

  1. If the outbreak of clusters is left unattended, the total number of infected people will increase explosively. ** Clusters with 10 initial reproductions have 113, clusters with 20 initial reproductions have 3493, and clusters with 30 initial reproductions. Produces 58,437 infected people **.
  2. The effect of cluster countermeasures can be expected as soon as possible and the larger the cluster size, ** there is a possibility that the total number of infected people can be halved or more if isolated within 20 days from the time of cluster occurrence **. is there.
  3. However, ** If 50 days have passed since the cluster occurred, the effect of the cluster measures will almost disappear **.
  4. Therefore, if the number of infected people who cannot follow the route is increasing, it is rational to ban uniform going out at an early stage.
  5. The factors that make the difference between the G7 countries in Europe and the United States and Japan are: ** Moved to cluster measures early, restrained clusters with high initial production by refraining from large-scale events, and in addition to religious facilities. It is possible that there are few meetings **.

Limitations of conventional model (SEIR model)

In the article below, I have calculated using a model called the SEIR model, which calculates the transition of infectious diseases.

The SEIR model is simply S: Susceptible for infectious diseases, E: Exposed for infectious diseases, I: Infectious, R: Infectious diseases It is a model that expresses how the distribution of a certain group changes among the four people who recovered from the disease and acquired immunity (Recovered) by a differential equation. It seems that it is often used in the trend analysis of infectious diseases, and it often appears in the explanation materials of the Ministry of Health and Welfare in Japan, the Ministry of Health in the United Kingdom, and the US CDC. However, if you try to calculate the population group S based on a large group such as the whole country or the whole prefecture, it will not be a digit at all **. For example, in the calculation of the SEIR model examined in the above article, assuming a population of 1000 people, it peaks in 60 days, and at the peak, more than 30% of the population becomes infected, but in terms of population ratio Even if you look at Italy, where the number of infected people is extremely high, the current number of infected people is only about 0.1% of the population of 60.48 million. There are various possible reasons, but the following are assumed, for example.

So, let's think about ** Is there a model that can explain the facts more from the actually observed data **?

Consider a cluster-based SEIR model

Therefore, let's check the phenomenon that was actually observed. The figure below is the result calculated by Let's examine the convergence time from the global trend of the effective reproduction number of the new coronavirus. R0_爆発的感染が観測された地域+推定.png

The figures in the graph are calculated based on a certain standard for the effective basic reproduction number (** below, referred to as Rp. **). The dotted line is an approximate expression (log-linear) that approximates the graph. The trends that can be read from this graph are as follows.

So ** let's obediently model this observed trend **. In particular,

If you write it in a mathematical formula, it will look like this.

Rp(t) = Rp(t_0) \cdot 2^{-\frac{t-t_0}{T_h}},  T_h=7.5[days]

Incorporate this effective reproduction number Rp (t) into the SEIR model to make it a ** cluster-based SEIR model **.

\begin{eqnarray}
\frac{dRp}{dt} &=& - \frac{ln2}{T_h}Rp \\
\frac{dS}{dt} &=& -\frac{1}{ip} Rp \cdot I \\
\frac{dE}{dt} &=& -\frac{dS}{dt} - \frac{1}{lp} E \\
\frac{dI}{dt} &=& \frac{1}{lp} E - \frac{1}{ip} I \\
\frac{dR}{dt} &=& \frac{1}{ip} I
\end{eqnarray}

The point is that, unlike the original SEIR model, it does not depend on the initial value of ** S (t) **. In the following, I would like to examine the effect of cluster countermeasures on this model, that is, the effect of discovering clusters and discovering and isolating infected people originating from clusters **.

Try to calculate with Python

Let's calculate the above cluster-based SEIR model in Python. Import the library.

import numpy as np
import matplotlib.pyplot as plt

Define the ODE. Here, in order to calculate the effect of cluster countermeasures, after a certain time $ T_c $, it is calculated with $ Rp = 0 $.

#define differencial equation of seir model
def seir_eq7(v, t, keys):
    Th = keys['Th']
    lp = keys['lp']
    ip = keys['ip']
    Tc = keys['Tc']
    #
    Rp = v[0];
    s  = v[1];
    e  = v[2];
    i  = v[3];
    r  = v[4];
    #
    if t >= Tc:
        Rp = 0
    #
    dRp = - np.log(2)/Th * Rp
    ds  = - Rp/ip * i
    de  = - ds - (1/lp) * e
    di  = (1/lp)*e - (1/ip)*i
    dr  = (1/ip)*i
    return [dRp, ds, de, di, dr]

Define a function to solve the ODE.

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
    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))
    return sim

This is the code to simulate by changing the initial cluster value $ Rp (t_0) $ in various ways. The parameters are $ lp = 5, ip = 8, Th = 7.5 $.

# case 7.1
def calcsim(Rp0, keys):
    # solve seir model
    #          Rp  S  E  I  R
    ini_state=[Rp0, 0, 0, 1, 0]
    t_max=180
    dt=0.01
    t=np.arange(0,t_max,dt)
    #
    sim = my_odeint(seir_eq7, ini_state, t, keys)
    #
    plt.rcParams["font.size"] = 12
    fig, ax = plt.subplots(figsize=(10,5))
    ax.plot(t,sim[:,[2,3,4]]) # extract Rp, E I R
    ax.set_xticks(np.linspace(0,t_max,19))
    ax.grid(which='both')
    ax.legend([ 'Exposed', 'Infected','Recoverd'])
    ax.set_xlabel('date')
    ax.set_ylim(0,)
    plt.show()
    print("R(tmax):{}".format(sim[-1,4]))

# try different Rp0
keys = {'lp':5, 'ip':8, 'Th':7.5, 'Tc':200 }
calcsim(100, keys)
calcsim(30, keys)
calcsim(20, keys)
calcsim(10, keys)
# try different Tc
keys = {'lp':5, 'ip':8, 'Th':7.5, 'Tc':20 }
calcsim(10, keys)
keys = {'lp':5, 'ip':8, 'Th':7.5, 'Tc':10 }
calcsim(10, keys)

This is the code to simulate by changing the cluster countermeasure time $ T_c $ in various ways.

# case 7.2
def calcTctoRp(Rp0):
    #          Rp  S  E  I  R
    ini_state=[Rp0, 0, 0, 1, 0]
    t_max=180
    dt=0.01
    t=np.arange(0,t_max,dt)
    #
    lp = 5
    ip = 8
    keys = {'lp':lp, 'ip':ip, 'Th':7.5, 'Tc':(5+0) }
    #
    rslt = []
    for i in range(0,60):
        keys['Tc'] = lp + i
        sim = my_odeint(seir_eq7, ini_state, t, keys)
        rslt.append([keys['Tc'], sim[-1:,4]]) # (i, R(tmax))
    #
    rslt = np.array(rslt)
    ymax = max(rslt[:,1])
    plt.rcParams["font.size"] = 12
    fig, ax = plt.subplots(figsize=(10,5))
    #
    ax.plot([0 , 0],[0,ymax],'m:')
    ax.plot([lp,lp],[0,ymax],'b:')
    ax.plot([lp+ip,lp+ip],[0,ymax],'g:')
    #
    ax.plot( rslt[:,0], rslt[:,1], 'r')
    ax.legend([ '0 day', 'lp days','lp + ip days','Total infected'], loc='lower right')
    ax.grid(which='both')
    ax.set_xlabel('cluster shutdown date since cluster occured')
    ax.set_ylabel('R(tmax)')
    #
    for tat in [10,13,20,30,40]:
        idx = [i for i in range(len(rslt[:,0])) if rslt[i,0] >= tat][0]
        print("R_fin with Tc:{} is {}".format(rslt[idx,0], rslt[idx,1]))
    #
    ax.set_xlim(-1,)
    ax.set_ylim(0,)
    plt.show()
    
calcTctoRp(100)
calcTctoRp(30)
calcTctoRp(20)
calcTctoRp(10)

simulation result

Let's take a look at the calculation result. The result of changing some patterns of the initial reproduction number $ Rp (t_0) $ and changing the cluster countermeasure time $ T_c $ by fixing $ Rp (t_0) $ is shown.

1. Effect of initial reproduction number Rp (t_0)

When the initial reproduction number Rp (t_0) = 10

Eventually 113 people will be infected. m7_1_Rp0_10_Tc_200.jpg

When the initial reproduction number Rp (t_0) = 20

Eventually 3493 people will be infected. m7_1_Rp0_20_Tc_200.jpg

When the initial reproduction number Rp (t_0) = 30

Eventually 58,437 people will be infected. m7_1_Rp0_30_Tc_200.jpg

2. Impact of cluster countermeasure time Tc

When the initial reproduction number Rp (t_0) = 10, the cluster countermeasure time Tc is changed and verified. Of course, the sooner the cluster measures are taken, the more effective it is expected. In the above simulation, 113 people were eventually infected without clustering measures.

When the cluster countermeasure time Tc = 10 [days]

If the infected person can be quarantined within 10 days of the cluster outbreak, the final infected person can be reduced from 113 to 21. m7_1_Rp0_10_Tc_10.jpg

When the cluster countermeasure time Tc = 20 [days]

If the infected can be quarantined within 20 days of the cluster, the final number of infected can be reduced from 113 to 64. m7_1_Rp0_10_Tc_20.jpg

3. Relationship between cluster countermeasure time Tc and final number of infected persons

Let's graph the relationship between the cluster countermeasure time Tc and the final number of infected people. Dotted lines are drawn at three points: when a cluster occurs, when the incubation period lp has elapsed, and when the incubation period + infection period (lp + ip) has elapsed.

When the initial reproduction number Rp (t_0) = 10

m7_2_Rp0_10.jpg

When the initial reproduction number Rp (t_0) = 20

Compared to when Rp (t_0) = 10, the effect of cluster measures is relatively increased. m7_2_Rp0_20.jpg

When the initial reproduction number Rp (t_0) = 30

Compared to when Rp (t_0) = 20, the effect of cluster measures is relatively increased. m7_2_Rp0_30.jpg

Consideration

From the above, the following trends can be derived from the simulation regarding the transition of the number of infected persons and the effect of the cluster countermeasure time based on the cluster-based SEIR model.

Even in the cluster-based SEIR model that does not assume the initial value of + S, the curve has almost the same shape as the normal SEIR model.

Furthermore ...

Reference link

I referred to the following page.

  1. SEIR model
  2. Outline of cluster support strategy (provisional version on March 10, 2020)

Recommended Posts

Let's explain the difference in how the new coronavirus spreads between Japan and other countries with a cluster countermeasure model
Verify the effect of leave as a countermeasure against the new coronavirus with the SEIR model
[Python] Explain the difference between strftime and strptime in the datetime module with an example
How to get the date and time difference in seconds with python
About the difference between "==" and "is" in python
How to use argparse and the difference between optparse
Let's simulate the effect of introducing a contact tracking app as a countermeasure against the new coronavirus
Let's take a look at the infection tendency of the new coronavirus COVID-19 in each country and the medical response status (additional information).
What is the difference between a symbolic link and a hard link?
Save the pystan model and results in a pickle file
I tried to find out the difference between A + = B and A = A + B in Python, so make a note