On February I decided to attend the Model Thinking Coursea’s MOOC. After learning the basics of ‘contagion models’, I decided to dig a little more into the topic of SIR models  (Susceptible-Infected-Removed).
Being a total noob of this field I surfed the internet: I was nicely surprised to find the paper of Munz et al.,2009 where the SIR approach has been used to describe the zombie apocalypse.
I am a big fan of zombie movies  (The Walking Dead is one of my preferred TV serie) so I decided to  translate the basic model by Munz et al. in R, with some further simplifications introduced by the authors.
OK, let’s start with a little bit of theory about an outbreak of a zombie infection:

• We know from the zombie films that we are all infected and so we are all in the Susceptible class (S).  
• When one in the S class deceases, for natural causes or by a zombie attack, moves to the Removed class (R). 
• People in the R class can resurrect and become zombie (parameter ζ): from this moment a zombie (now in the Z class) start to wander craving for for human flesh and trying to bite the susceptibles. 
• As explained the paper, people in class S can avoid zombification through an altercation with a zombie by defeating it, cutting its head or destroying its brain and so resisting the infection (parameter α). 
• In the unfortunate case that a Susceptible succumbs to an attack by an undead, the infection transmits from the zombie to the Susceptible, the latter will die and  briefly resurrect as member of the zombie class (transmission parameter β).

The flowchart of the SZR model s described in the image below.

Flowchart of the SZR model

Thefollowing system of nonlinear differential equation is the formalization of the SZR model:


To solve the model in R I’ll use the deSolve package by Soetaert & Petzoldt, thanks to this great piece of software, the formalization is pretty just straightforward, like writing the equations with a pencil on a sheet of paper:

require(deSolve)

SZR = function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    dS = -beta*S*Z
    dZ =  beta*Z*S + zeta*R - alpha*S*Z
    dR =  alpha*S*Z - zeta*R
    
    return(list(c(dS, dZ, dR)))
  })
}

times = seq(from= 0, to= 10, by= 0.01)

state = c(S= 500, Z= 0.0, R= 1.0)

parameters = c(alpha= 0.03, beta= 0.0445, zeta=0.2)

output = lsoda(y=state, times=times, func=SZR, parms= parameters)

As written in the code above, the initial state is represented by 500 Susceptibles and just one Removed (case 0). The parameters used in the model are α=0.03, β=0.0455, ζ=0.2, as proposed by Witkowski and Blais, 2013 analyzing the evolution of the infection from Shaun of the Dead and Night of the Living Dead (hey, data collecting can be fun!)

The results of the modelled outbreak is shown below:

The outbreak of the zombie infection, the solution of the SZR basic (and simplified)	model

No hope here: this model predicts the doomsday with zombies quickly infecting everyone in a very short time  resulting in the collapse of the whole society.

Witowsky & Blais are very critic about the basic model proposed by Munz et al. and indeed three others models are proposed in the paper. I just agree with the critics: the ζR assumption is not reasonable (as far we can use the word ‘reasonable’ talking about a zombie outbreak): with this parametrization, when zombies are killed they can be recycled from the R group to the Z group again. Instead, when a zombie is defeated by destroying its brain it cannot resurrect anymore.

But for this first try, i think it is just enough. Lessons learned so far:

1. the SIR models are fun.
2. the SIR models with zombies are even funnier (at least for me, a zombie geek).
3. the deSolve package for R is really powerful but surprisingly easy to use.

In the next post I will try to implement a different and more realistic model for the outbreak of the zombie infection and see the differences with the one just proposed here. Stay tuned, and just in case you meet a zombie… just run and ‘don’t look back’.