Inspired by Jason'sAlmost Looks Like WorkInspired by this, I'm going to show some models of viral transmission. It is important to note that this model does not reflect reality, so do not mistake it for a horrific West African epidemic. Instead, it should be viewed more as some sort of fictionalized zombie outbreak. So, let's get to the point.
That'sSIRmodel, where the letters S, I, and R reflect the different states that individuals may be in during a zombie epidemic.
- S represents the susceptible group, i.e., the number of potentially possible shifts in healthy individuals.
- I represents the infected population, i.e., the number of zombies.
- R stands for removal, i.e. the number of zombies that exit the game due to death, or the number that turn back to humans after infection. But there is no healer for vs. zombies, so let's not fool ourselves (there is still a healer if we want to apply the SIR model to flu infections).
- As for beta (beta) and gamma (gamma).
- Beta (beta) indicates the degree of contagiousness of the disease, which can be contracted simply by being bitten.
- Gamma(gamma) denotes the rate of moving from zombie to dead, depending on the average work rate of zombie hunters, of course, this is not a perfect model, please be patient with me.
- S′=?βIS tells us the rate at which healthy people become zombies, and S′ is the derivative with respect to time.
- I′ = βIS?γI tells us how the infected increase, and the rate at which walkers enter the removal state (pun intended).
- R′ = γI is just added to (gamma I), which is negative in the previous equation.
The model above does not consider the spatial distribution of S/I/R, which is corrected below!
One approach is to partition Sweden and the Nordic countries into grids where each cell can infect neighboring cells, as described below:
Where for the unit, and are the four units around it. (Don't get brain fatigue from diagonal units, we need our brains to not get eaten).
Initialize some stuff.
import numpy as np import math import as plt %matplotlib inline from matplotlib import rcParams import as mpimg rcParams[''] = 'serif' rcParams[''] = 16 rcParams[''] = 12, 8 from PIL import Image
The right beta and gamma values are capable of destroying a large portion of the country.
beta = 0.010 gamma = 1
Remember the definition of the derivative? When the derivative is known, and assuming Δt is small, after rearranging it, it can be used to approximate the next value of the predicted function, which we have already stated as u′(t).
Initialize some stuff.
import numpy as np import math import as plt %matplotlib inline from matplotlib import rcParams import as mpimg rcParams[''] = 'serif' rcParams[''] = 16 rcParams[''] = 12, 8 from PIL import Image
The right beta and gamma values are capable of destroying a large portion of the country.
beta = 0.010 gamma = 1
Remember the definition of the derivative? When the derivative is known, and assuming Δt is small, after rearranging it, it can be used to approximate the next value of the predicted function, which we have already stated as u′(t).
This method is calledEuler method (math.)The code is as follows:
def euler_step(u, f, dt): return u + dt * f(u)
We need the function f(u). The friendly numpy provides concise array operations. I'll probably review it in another post because they're so powerful they need more explanation, but for now this will do the trick:
def f(u):
S = u[0]
I = u[1]
R = u[2]
new = ([-beta*(S[1:-1, 1:-1]*I[1:-1, 1:-1] +
S[0:-2, 1:-1]*I[0:-2, 1:-1] +
S[2:, 1:-1]*I[2:, 1:-1] +
S[1:-1, 0:-2]*I[1:-1, 0:-2] +
S[1:-1, 2:]*I[1:-1, 2:]),
beta*(S[1:-1, 1:-1]*I[1:-1, 1:-1] +
S[0:-2, 1:-1]*I[0:-2, 1:-1] +
S[2:, 1:-1]*I[2:, 1:-1] +
S[1:-1, 0:-2]*I[1:-1, 0:-2] +
S[1:-1, 2:]*I[1:-1, 2:]) - gamma*I[1:-1, 1:-1],
gamma*I[1:-1, 1:-1]
])
padding = np.zeros_like(u)
padding[:,1:-1,1:-1] = new
padding[0][padding[0] < 0] = 0
padding[0][padding[0] > 255] = 255
padding[1][padding[1] < 0] = 0
padding[1][padding[1] > 255] = 255
padding[2][padding[2] < 0] = 0
padding[2][padding[2] > 255] = 255
return padding
Importing and down-sampling population density maps for the Nordic countries to get faster results
from PIL import Image
img = ('')
img = (([0]/2,[1]/2))
img = 255 - (img)
imgplot = (img)
imgplot.set_interpolation('nearest')
Population density map for the Nordic countries (without Denmark)
The S-matrix, or susceptible individuals, should approximate the population density. The initial value of infected individuals is 0, and we use Stockholm as the first source of infection.
S_0 = img[:,:,1] I_0 = np.zeros_like(S_0) I_0[309,170] = 1 # patient zero
Since no one has died yet, set the matrix to 0 as well.
R_0 = np.zeros_like(S_0)
Then initialize the simulation duration, etc.
T = 900 # final time dt = 1 # time increment N = int(T/dt) + 1 # number of time-steps t = (0.0, T, N) # time discretization # initialize the array containing the solution for each time-step u = ((N, 3, S_0.shape[0], S_0.shape[1])) u[0][0] = S_0 u[0][1] = I_0 u[0][2] = R_0
We need to customize a color table so that the infection matrix can be displayed on the map.
import as cm
theCM = cm.get_cmap("Reds")
theCM._init()
alphas = ((0, 1, ))
theCM._lut[:-3,-1] = alphas
Sit back and enjoy below...
for n in range(N-1): u[n+1] = euler_step(u[n], f, dt)
Let's do a little more image rendering and make it into a gif, everyone loves gifs!
from images2gif import writeGif
keyFrames = []
frames = 60.0
for i in range(0, N-1, int(N/frames)):
imgplot = (img, vmin=0, vmax=255)
imgplot.set_interpolation("nearest")
imgplot = (u[i][1], vmin=0, cmap=theCM)
imgplot.set_interpolation("nearest")
filename = "outbreak" + str(i) + ".png"
(filename)
(filename)
images = [(fn) for fn in keyFrames]
gifFilename = ""
writeGif(gifFilename, images, duration=0.3)
()