Jens Hahn - 01/12/2015
Content taken from:
Scipy Docs at http://scipy-cookbook.readthedocs.org/items/Zombie_Apocalypse_ODEINT.html
Munz et al. (2009): http://mysite.science.uottawa.ca/rsmith43/Zombies.pdf
Differential equations can be used to describe the time-dependent behaviour of a variable.
$$\frac{\text{d}\vec{x}}{\text{d}t} = \vec{f}(\vec{x}, t)$$
In our case the variable stands for the number of humans in a infected (zombies) or not infected population.
Of course they can also be used to describe the change of concentrations in a cell or other continuous or quasi-continuous quantity.
In general, a first order ODE has two parts, the increasing (birth, formation,...) and the decreasing (death, degradation, ...) part:
$$\frac{\text{d}\vec{x}}{\text{d}t} = \sum_{}\text{Rates}_{\text{production}} - \sum_{}\text{Rates}_{\text{loss}}$$You probably already know ways to solve a differential equation algebraically by 'separation of variables' (Trennung der Variablen) in the homogeneous case or 'variation of parameters' (Variation der Konstanten) in the inhomogeneous case. Here, we want to discuss the use of numerical methods to solve your ODE system.
Let's have a look at our equations:
Number of susceptible victims $S$: $$\frac{\text{d}S}{\text{d}t} = \text{P} - \text{B}\times S \times Z - \text{d}\times S$$
Number of zombies $Z$: $$\frac{\text{d}Z}{\text{d}t} = \text{B}\times S \times Z + \text{G}\times R - \text{A}\times S \times Z$$
Number of people "killed" $R$: $$\frac{\text{d}R}{\text{d}t} = \text{d}\times S + \text{A}\times S \times Z - \text{G}\times R$$
Parameters:
P: the population birth rate
d: the chance of a natural death
B: the chance the “zombie disease” is transmitted (an alive person becomes a zombie)
G: the chance a dead person is resurrected into a zombie
A: the chance a zombie is totally destroyed by a human
Before we start the simulation of our model, we have to define our system.
We start with our static information:
In [1]:
import numpy as np
# 1. initial conditions
S0 = 500. # initial population
Z0 = 0 # initial zombie population
R0 = 0 # initial death population
y0 = [S0, Z0, R0] # initial condition vector
# 2. parameter values
P = 0 # birth rate
d = 0.0001 # 'natural' death percent (per day)
B = 0.0095 # transmission percent (per day)
G = 0.0001 # resurect percent (per day)
A = 0.001 # destroy percent (per day)
# 3. simulation time
start = 0.0 # days
end = 15. # days
number_time_points = 1000.
t = np.linspace(start, end, number_time_points) # time grid, 1000 steps or data points
In the second step, we write a small function f, that receives a list of the current values of our variables x and the current time t. The function has to evaluate the equations of our system or $\frac{\text{d}\vec{x}}{\text{d}t}$, respectively. Afterwards, it returns the values of the equations as another list.
Important
Since this function f is used by the solver, we are not allowed to change the input (arguments) or output (return value) of this function.
In [2]:
# function 'f' to solve the system dy/dt = f(y, t)
def f(y, t):
Si = y[0]
Zi = y[1]
Ri = y[2]
# the model equations (see Munz et al. 2009)
f0 = P - B*Si*Zi - d*Si
f1 = B*Si*Zi + G*Ri - A*Si*Zi
f2 = d*Si + A*Si*Zi - G*Ri
return [f0, f1, f2]
Last but not least, we need to import and call our solver. The result will be a matrix with our time courses as columns and the values at the specified time points. Since we have a values for every time point and every species, we can directly plot the results using matplotlib.
In [5]:
# zombie apocalypse modeling
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.integrate import odeint
# solve the DEs
result = odeint(f, y0, t)
S = result[:, 0]
Z = result[:, 1]
R = result[:, 2]
# plot results
plt.figure()
plt.plot(t, S, label='Humans')
plt.plot(t, Z, label='Zombies')
plt.plot(t, R, label='Dead Humans')
plt.xlabel('Days from outbreak')
plt.ylabel('Population')
plt.title('Zombie Apocalypse - No Init. Dead Pop.; No New Births.')
plt.ylim([0,500])
plt.legend(loc=0)
Out[5]:
When you want to simulate biochemical reactions it is in principle not different to the things you've learned right now, but there is one important thing to consider. People have thought a lot about the best ways to describe certain kinds of chemical reactions and their kinetics. Let's have a look at one of the 'rates' from above: We wrote the birth rate of humans as $$v_{\text{birth}} = \text{P}$$ that means, there is a constant addition of human beeings in our system, or a constant 'influx'.
Then, we have the 'natural' death of humans $$v_{\text{death}} = \text{d}\times S$$
The extent of the rate is dependent on the current number of humans S, makes sense, right? The more human beeings there are, the more can die of 'natural' cause.
Finally, there is the killing of zombies, this rate is not only dependent on the number of zombies Z (the more there are, the more can be killed), but also on the number of humans S, the more there are, the more can kill zombies:
$$v_{\text{kill}} = \text{A}\times S \times Z$$
This kind of 'rate law' we call mass action, there is always a parameter multiplied by one or more variables.
The other very important kinetic rate law is the Michaelis-Menten kinetic, you've probably heard about it. This rate law was made for describing an enzymatic reactions, when the rate itself cannot reach unlimited extent. The velocity or flux is limited by the maximal rate of the enzyme: $$v_{\text{MM}} = \text{V}_\text{max}\times \frac{S}{\text{K}_{\text{m}} + S}$$
Equations: $$\frac{\text{d}R}{\text{d}t} = \alpha\times R - \beta\times R \times F$$ $$\frac{\text{d}F}{\text{d}t} = -\gamma\times F + \delta\times\beta\times R \times F$$
Parameters:
$\alpha$: natural growing rate of rabbits, when there's no fox around
$\beta$: natural dying rate of rabbits, due to predation
$\gamma$: natural dying rate of fox, when there's no rabbit
$\delta$: factor describing how many caught rabbits let create a new fox
a) Find a reasonable parameter set
b) Plot the time courses using matplotlib
c) Plot the phase plot (What is a phase plot?)
a) Use Michaelis-Menten kinetics to simulate this short reaction chain: $$\text{A} \longrightarrow \text{B} \longrightarrow \text{C}$$