[PYTHON] [Introduction to simulation] I tried playing by simulating corona infection ♬ Part 2

This time, it is a sequel to Simulation summarized on March 18. Actually, I felt like I was done, but when I thought about it, I had to dig deeper into the following two points. ① Infection transmission that cannot be explained by the SIR model ② Propagation when it occurs in a cluster In other words, in (1), looking at the transition of the number of infections in the previous particle density-dependent transition region, the peak number of infections is divided into several peaks, which cannot be explained by the SIR model. Also, real-life infections often occur in clusters and somehow do not appear to be monotonously transmitted. So, this time, we have implemented the following. ・ Transmission of infection in the critical density region ・ Infection transmission when the distribution of susceptibility carriers is distributed in a cluster The code is almost the same as the previous one, but it has been extended to give speed to the distribution and initial values of susceptibility holders. I'll put the code as a bonus.

• This time, I will omit the simulation to give speed.

・ Transmission of infection in the critical density region

When the susceptibility holders are randomly and uniformly distributed as in the previous time ① Infection range; Rc ② Infection probability; p ③ Density of susceptibility holders; d ④ Movement speed of susceptibility holder; mr Infection transmission and infection rate can be controlled by changing.

Typical infection transmission

The figure on the right is blue; number of susceptible carriers, red; infected, green; healer, and this graph is a familiar picture in the familiar SIR model. Parameters: ①Rc = 20x20, ②p = 0.05, ③d = 1000 / 200x200, ④mr = 1.0

Critical transmission

If the infection probability is lowered as follows, infection transmission will not occur soon. Infection transmission just before that can be realized with the following parameters, and it will be like the following Gif animation.

Parameters: ①Rc = 20x20, ②p = 0.0075, ③d = 1000 / 200x200, ④mr = 1.0

As can be seen from this figure, as shown in the figure on the right, this propagation is in a region that can no longer be represented by the SIR model, and it can be seen that propagation is occurring stochastically. In this figure, the rise is exponential, but after that, the infection spreads and heals by slowly pulling the tail.

Parameters: ①Rc = 13x13, ②p = 0.03, ③d = 1000 / 200x200, ④mr = 1.0

Critical transmission of infection could be achieved even if the infection range was reduced, and the result was as follows. In this animated gif, even at the start, the exponential increase disappears, it increases slowly, decreases shabby, and increases again. Thus, transmission of infection in this region depends on the distribution of susceptibility carriers and can be said to be an almost stochastic process.

・ Infection transmission when distributed in a cluster

In the actual distribution, it is assumed that the population density usually fluctuates in each village or town, and that the distribution fluctuates in each gathering or gathering of people. So, I tried to spread the infection when the distribution of susceptibility carriers is distributed in a cluster.

When clusters are densely present

This time, for the sake of simplicity, we first looked at infection transmission in clusters distributed on a 10x10 grid.

Parameters: ①Rc = 20x20, ②p = 0.5, ③d = 100x10 / 200x200, ④mr = 1.5

The meaning of d = 100x10 / 100x100 means that 100 clusters are placed in the area of 200x200 and 10 sensitivity holders are placed at each point on the grid as shown in the figure below. With this parameter, infection transmission is a region that can be expressed by the SIR model, which is the same as the uniform distribution.

Parameters: ①Rc = 20x20, ②p = 0.5, ③d = 100x10 / 200x200, ④mr = 1.25

When the speed of the random walk of the susceptibility carrier decreases, the probability of contact with the infected person decreases, and the transmission between clusters occurs depending on the infection by the susceptibility carrier as shown below, and the following is not possible. Regular propagation occurs.

Parameters: ①Rc = 20x20, ②p = 0.5, ③d = 100x10 / 200x200, ④ mr = 1.09

And the speed around 1.09 is the boundary of whether it is critical and propagates. And when the speed becomes lower than this, the propagation disappears. In other words, it will not be infected. Before infection, transmission stops at the first infected person or at least one of the adjacent clusters.

If the clusters are sparse

What if the clusters are larger and they are far apart? That is, the image is when large villages are scattered.

Parameters: ①Rc = 20x20, ②p = 0.5, ③d = 4x250 / 200x200, ④mr = 10

The meaning of d = 4x250 / 100x100 means that 200x200 is interspersed with clusters containing four 250 susceptibility carriers. In this case, it collapses a little due to the distribution of clusters, but it becomes a single peak for the time being.

Parameters: ①Rc = 20x20, ②p = 0.5, ③d = 4x250 / 200x200, ④mr = 7.5

In the case of this cluster distribution, this parameter area gives critical propagation in the sense of propagation between clusters.

When the cluster is sparse and the initial infected person is unevenly distributed on the lower side and the velocity distribution is uneven

It looks like a complicated setting when written in letters, but it seems to be possible in reality.

Parameters: ①Rc = 20x20, ②p = 0.5, ③d = 4x250 / 200x200, ④mr = 7.5, ⑤mx = 1, my = 0.5

With the last parameter ⑤, the speed is increased in the x direction and the speed in the y direction is halved.

Parameters: ①Rc = 20x20, ②p = 0.0075, ③d = 4x250 / 200x200, ④mr = 7.5, ⑤mx = 1, my = 0.75

If you place the early infected person in the lower left village from the beginning, it will be as follows. And the infection probability is also set to a small value of 0.0075 (the value of critical transmission of uniform distribution). Even if the velocity in the y direction is 0.75, it can be seen that the susceptible carriers belonging to the upper village are hardly infected. In this area, it can be seen that the transition of the number of infections is very structured and complicated.

Discussion

What is important is that there are areas where the number of beautiful infections does not increase or decrease as seen in the SIR model. And in the case of critical propagation (1) Infection does not spread if the movement speed of infected persons and susceptibility carriers is small ② If the infection rate is lowered, infection transmission will not occur. ③ If the density of susceptibility holders is low, infection will not occur. Therefore, it is dangerous to mix (travel) infected and susceptible individuals. The infection rate may be reduced to 0 by masking (infected persons) or washing hands (susceptible persons). Well, some people are asymptomatic, so all masks need to be worn. Furthermore, above all, the practice of measures to avoid three honeys (crowded, closely, closed), such as not gathering because infection can be prevented simply by lowering the density below a certain level. It seems that infection can be avoided by practicing the measures that are said in the streets. Finally, although not shown in the calculations, I understand that the concept of so-called herd immunity is almost meaningless because it changes by simply changing the density and infection probability in such calculations.

Summary

・ The SIR model does not seem to hold in the critical propagation region. ・ Cluster-like distribution of susceptibility carriers results in structural infections ・ When clusters are distributed in villages, transmission of infection between villages is unlikely to occur. ・ In the case of critical propagation, it does not propagate if the velocity is low. ・ The concept of herd immunity is meaningless (changes depending on infection probability, infection range, susceptibility carrier distribution, etc.)

・ If the infected person has a speed, the rate of increase in infection will naturally increase, but the simulation has not been completed.

bonus

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import random
import time

PARTICLE_NO = 1000 #Number of particles
ITERATION = 200 #Maximum number of loops: Stops when the number of infected people reaches 0
MIN_X, MIN_Y = -100.0, -100.0 #Minimum range at the start of search
MAX_X, MAX_Y = 100.0, 100.0 #Maximum range at the start of search
recovery=30 #Healed after a certain period of time
p=0.5 #0.03 #probability of infecion
rc=100 #121 #169 #225 Range of infection Radius of circle^2

start = time.time()

def plot_particle(sk,positions,elt,r,g,b):
    #fig, ax = plt.subplots()
    el_time = time.time()-start
    fig, (ax1, ax2) = plt.subplots(1, 2, sharey=False,figsize=(8*2, 8))
        
    for j in range(0,PARTICLE_NO):
        x=positions[j]["x"]
        y=positions[j]["y"]
        c=positions[j]["c"]
        s = 5**2
        ax1.scatter(x, y, s, c, marker="o")
    ax1.set_xlim([MIN_X, MAX_X])
    ax1.set_ylim([MIN_Y, MAX_Y])
    ax1.set_xlabel("x")
    ax1.set_ylabel("y")
    ax1.set_title("{:.2f}:InfectionRate;{:.2f} %".format(el_time,(PARTICLE_NO-b[-1])/PARTICLE_NO*100))
    
    ind = np.arange(len(elt))  # the x locations for the groups
    width = 0.3       # the width of the bars

    ax2.set_ylim([0, PARTICLE_NO])
    ax2.set_title("{:.2f}:red_{} green_{} blue_{}".format(el_time,r[-1],g[-1],b[-1]))
    rect1 = ax2.bar(ind, b,width, color="b")
    rect2 = ax2.bar(ind+width, g, width, color="g") #, bottom=b)
    rect3 = ax2.bar(ind+2*width, r,width, color="r") #, bottom=b)
    plt.pause(0.1)
    plt.savefig('./fig/fig{}_.png'.format(sk)) 
    plt.close()
    
#Particle position update function
def update_position(positions,velocity):
    x0 = []
    y0 = []
    for i in range(PARTICLE_NO):
        c=positions[i]["c"]
        t_time = positions[i]["t"]  #Initial value 0, infection time at infection
        k_time = time.time()-start  #elapsed time
        s = positions[i]["flag"]    #No infection 0, infection: 1
        if s == 1 and c == "red":   #If infected
            if k_time-t_time>recovery:  #Healed after a certain period of time
                #print("inside",i,s,c,k_time-t_time)
                c = "blue"
                positions[i]["c"] = "green"
                positions[i]["flag"] = 1   #However, the infection history remains
        if c == "red":  #Get location information if infected red
            x0.append(positions[i]["x"])
            y0.append(positions[i]["y"])
            #print("4",i,s,c,t_time)
    #print(x0,y0)   
    position = []
    for j in range(PARTICLE_NO):
        x=positions[j]["x"]
        y=positions[j]["y"]
        c=positions[j]["c"]
        s = positions[j]["flag"]
        t_time = positions[j]["t"]
        for k in range(len(x0)):
            if (x-x0[k])**2+(y-y0[k])**2 < rc and random.uniform(0,1)<p:
                if s ==0:
                    c = "red"
                    t_time = time.time()-start
                    s = 1
                    positions[j]["flag"]=s
                else:
                    continue
        vx = velocity[j]["x"]+1.085*random.uniform(-1, 1) #Coefficient is the magnitude of particle motility
        vy = velocity[j]["y"]+1.085*random.uniform(-1, 1)
        new_x = x + vx
        new_y = y + vy
        p_color = c
        s=s
        
        position.append({"x": new_x, "y": new_y, "c": p_color, "t": t_time,"flag":s})
        velocity.append({"x": vx, "y": vy})

    return position, velocity, x0

def count_brg(position):
    r=0
    g=0
    b=0
    for j in range(len(position)):
        if position[j]["c"] == "red":
            r += 1
        elif position[j]["c"] == "green":
            g += 1
        else:
            b += 1
    return r,g,b        

def main():
    #Start time measurement
    #start = time.time()
    xy_min, xy_max = -32, 32
    #Initial position of each particle,speed, personal best,global best and search space settings
    position = []
    velocity = []  #Expanded speed for use
    
    #Initial position,Initial speed
    #position.append({"x": random.uniform(MIN_X, MAX_X), "y": random.uniform(MIN_Y, MAX_Y), "c": "red", "t":0, "flag":1})
    position.append({"x": 0, "y": 0, "c": "red", "t":0, "flag":1}) #Place one initial infected person in the middle (0,0)
    velocity.append({"x": 0, "y": 0}) #The initial speed of the infected person is set to 0.
    for i in range(0,10): #Arrange 10 villages in the x direction
        for j in range(0,10): #Arrange the villages with mesh 10 in the xy direction
            for k in range(0,10): #Distribution of 10 susceptibility holders per village
                s=k+j*10+i*100;
                position.append({"x": 10+(-100+i*20)+random.uniform(MIN_X/100, MAX_X/100), "y":10+(-100+j*20)+ random.uniform(MIN_Y/100, MAX_Y/100), "c": "blue", "t": 0, "flag":0})
                velocity.append({"x": 0, "y": 0})
    print(len(position))
    sk = 0
    red=[]
    green=[]
    blue=[]
    elapsed_time = []
    while sk < ITERATION:
        position,velocity, x0 = update_position(position,velocity) ######
        r,g,b = count_brg(position)
        red.append(r)
        green.append(g)
        blue.append(b)
        el_time=time.time()-start
        elapsed_time.append(el_time)
        #print("{:.2f}:red_{} green_{} blue_{}".format(el_time,r,g,b))
        plot_particle(sk,position,elapsed_time,red,green,blue)
        if x0==[]:
            break
        sk += 1

    #Time measurement finished
    process_time = time.time() - start
    print("time:", process_time)
    
if __name__ == '__main__':
    main()