APPM Lab 3: Zombies!

    • Introduction

As you already know, in the past month, zombies have overrun much of North America, including all major cities on both the East and the West coast. You have been selected as lead scientists for the newly formed Zombie Strategic Defense Program for two reasons. First, previous studies have shown that, since we are isolated from most major metropolitan areas (such as New York and Los Angeles), Boulder is one of the safest places to be in the event of a zombie apocalypse. Second, and perhaps more important, you are the only surviving humans that understand di erential equations! Your analysis and recommendations will inform the nation’s response to this zombie outbreak. Remember that this outbreak has the potential to spread rapidly. We are forecasting that you will probably survive, at least until nals week. Therefore, your written report (in .pdf form) and all supporting code must be uploaded to D2L by Friday, December 8, 2017 before 5:00 pm.

    • Mathematical and Biological Background

2.1    SIR Models

There is a long history of di erential equation models in epidemiology; we will consider a type referred to as a compartment model. One very simple model for the spread of an epidemic separates a population of N individuals into three classes: susceptible, infectious (also infected in this model), and recovered. We denote these classes by S, I, and R, respectively. Susceptible individuals become infected by contact with infectious individuals. This interaction can be modeled by the Law of Mass Action. Assuming that an infectious person transmits the disease to N other people per unit time and that this infectious person encounters a susceptible person on S=N of all encounters, then the rate of infection is ( N)(S=N)I = SI. These models often assume that the recovery rate is proportional to the number of infectious persons. If T is the duration of infection, then the constant of proportionality is the reciprocal of the duration of infection = 1=T . Assuming that people are born at constant rate and die of natural causes at constant rate , this simple model can be expressed as a system of di erential equations:


dS

=  N
SI
S;


dt












dI

=  SI
I
I;


dt











dR

=  I
R:


dt











This is commonly known as an SIR model. Notice that the total population is N(t) = S(t) + I(t) + R(t). This means that our solutions will grow without bound if > , that is, if the birth rate exceeds the death rate due to natural causes. To avoid this problem, we will assume that birth and death rates are negligible on the fast timescale of the epidemic (which is days). That will allow us to set = = 0 to obtain a constant population model. This is a \nice" assumption because, when the population is a constant (that is, N(t) = N0), we only need to solve for S(t) and I(t) because R(t) = N0 S(t) I(t).

1

2.2    The Jacobian Matrix

In Section 4, we will ask you to investigate the Jacobian matrix of this system so, we need to de ne the Jacobian matrix for you. You may recall the term Jacobian from Calculus 3. For a vector-valued function f : R2 ! R2 which takes input x = [x1; x2] 2 R2 and outputs f(x) = [f1(x); f2(x)] 2 R2, the Jacobian matrix, J, of f is de ned as

    • 3


@f1

@f1



J =
4
@x1

@x2
5
:
(1)


@x1

@x2





@f2

@f2




A two-dimensional system of autonomous,    rst-order di erential equations can be written as

x0 (t) = f (x(t)) ;

where x = [x1; x2] is a vector of the solution variables and f = [f1; f2] is a vector of right hand side functions. The Jacobian matrix for this system is also given by equation (1). In particular, when

=  = 0, the SIR model from the previous page can be written as


dS

= f1(S; I) =
SI;








dt





dI

= f2(S; I) =  SI   I:







dt


The Jacobian matrix for this model is
3
2
3
2




@f1
@f1
J =
@S

@I

@S2

@I2
4




@f
@f




I    S
5=4    5:
I    S

2.3    Linear Stability

After we nd an equilibrium point, x0, we can evaluate the Jacobian matrix at the equilibrium point. Denote this by A = Jjx0 . Then the linearized system at x0 is

x0 (t) = Ax(t):

Why is this useful? We can compute the eigenvalues of A = Jjx0 . These eigenvalues tell us about the stability of solutions.

    • If the eigenvalues of the Jacobian matrix (evaluated at the equilibrium point) all have negative real part, then the equilibrium point is asymptotically stable. We expect that, if a solution starts close enough to an asymptotically stable equilibrium point, then the solution will converge to that equilibrium point as t ! 1.

    • On the other hand, if at least one of the eigenvalues of the Jacobian matrix (evaluated at the equilibrium point) has positive real part, then the equilibrium point is unstable.

See page 437 in the course text for additional information about classifying equilibria.






2
    • SZR Model

Several assumptions in the SIR model above do not apply to the zombie outbreak. In particular,

    • Zombies are not giving birth.

    • Human/zombie encounters may result in zombie death (in addition to human infection). We assume that the average susceptible human kills N zombies per unit time.

    • Zombies do not recover at this time. Instead, deceased humans and zombies are placed in the ‘removed’ class, R. Anyone in class R can be revived as a zombie with rate R. (This assumption is unrealistic. It means that every removed person will eventually revive as a zombie. Don’t worry about that too much. After all, this is a zombie project.)

Susceptible humans are born, become converted into zombies, or die of natural causes. When a human interacts with a zombie, either the human becomes a zombie or the zombie is destroyed. Zombies can also be revived from the removed class. These peculiarities of the zombie outbreak can be incorporated into a susceptible, zombie, and removed (SZR) model:


dS
=S  ZS
S;

dt







dZ
=ZS+R
ZS;

dt






dR
=S+ZS
R;

dt







The SZR model is represented schematically in Figure 1. Notice the S, Z, and R compartments.


S






S

Susceptibles





ZS






Zombies




SZ


R






Removed


Figure 1: A diagram of the SZR model.

Since the zombie epidemic takes place in a few days, let’s assume that = 0 and = 0. This means that the population is constant. If N(0) = N0 is the total population at time zero, then we know that R(t) = N0 S(t) Z(t), so we only need to solve for S(t) and Z(t) and the model reduces to

dS
=   ZS;
(2a)
dt





dZ
= ZS+ (N0   S  Z)   ZS:
(2b)



dt



We will analyze equations (2a) and (2b) in Section 4.


3
    • Analysis

Consider the SZR model represented by equations (2a) and (2b). Assume that    ;    ;    > 0.

    (a) Find the nullclines (you may do that by hand) and  nd the equilibrium points of the system.

    (b) Find the Jacobian for the system and evaluate the Jacobian at each equilibrium point.

    (c) Find the eigenvalues of the Jacobian evaluated at each equilibrium point.

    (d) Use the eigenvalues to classify the equilibria as asymptotically stable or unstable. The stability of one equilibrium point is not obvious. Find a condition on parameters ; ; ; and N0 that makes one of the equilibrium points unstable. What are the implications for human survival?

    (e) Solve equations (2a) and (2b) numerically for 35 days to see how the system evolves. Use initial conditions S(0) = 59; 999 and Z(0) = 1. Use the parameter values in Table 1. Plot S(t), Z(t), and R(t) on the same graph. MatLab code is provided below for your convenience.

    (f) Explain (in sentences) what the code that we provided is doing.

    (g) What happens to the population in your simulation? You should nd that the population is not anywhere near an equilibrium point at 35 days. Will the simulation ever reach an equilibrium point? Explain. (Hint: Think carefully about your answer to (d).)


SZR Parameters

Symbol
Value
Units
Description





0.00002
[1 / humanoids / day]
average zombie kill rate per human per day

0.00003
[1 / humanoids / day]
average infection rate per zombie per day

0.000025
[1 / day]
revival rate
N0
60,000
[humanoids]
total population size

Table 1: Use these parameters for your simulation in Section 4.

To make units simple, we will refer to humans, zombies, and removed persons as \humanoids." The following MatLab function creates the right side of the model. We call this with ode45.


function dydt = szr(t, y, alpha, beta, gamma, N0)

    • Evaluates the right hand side of the SZR model for equations S(t) and Z(t).

    • Here, y(t) = [S(t); Z(t)], so y(1) = S(t) and y(2) = Z(t).
dydt = [-beta*y(1)*y(2);
beta*y(1)*y(2) + gamma*(N0 - y(1) - y(2)) - alpha*y(1)*y(2)];
end

To solve equations (2a) and (2b) in MatLab, you may use the following script.


% Assign parameter values

alpha = 0.00002; beta = 0.00003; gamma = 0.000025; N0 = 60000;

    • Set length of simulation tspan = [0 35];

    • Set initial conditions y0 = [59999; 1];

    • Solve system

[t,y] = ode45(@(t,y) szr(t,y,alpha,beta,gamma,N0), tspan, y0);


4
    • How Can We Fix This Zombie Thing?

5.1    An Antidote!

Before the East coast was overrun by zombies, scientists at the CDC were developing an antidote to return zombies to their human forms. They believe they have found one that works and thanks to one forward thinking epidemiologist, some of it has been sent to you in Boulder! Though side e ects have yet to be determined, the antidote should work if a human can manage to inject it into
    • zombie.

        (a) Modify equations (2a) and (2b) to include the antidote under the following assumptions:

            ▪ The antidote will returns former zombies to the susceptible class. (In other words, it is a temporary cure that does not grant permanent immunity).

            ▪ The antidote returns zombies to the susceptible class at rate Z where is a parameter with units [1 / day].

See Figure 2 for a schematic diagram of the model with an antidote.

S






S

Susceptibles




ZS



Z






Zombies




SZ


R






Removed


Figure 2: A diagram of the SZR model with an antidote.


    (b) Solve the system numerically with = 0:2 for 35 days. Use initial conditions S(0) = 59999 and Z(0) = 1 and the parameter values from Table 1. Continue to assume that = = 0. Plot S(t) and Z(t).

    (c) Is this an e ective strategy? Why or why not?

    (d) Find the equilibrium points of the model with antidote and classify their stability, when = 0:2 and the other parameters are given by Table 1. Use this information to explain your simulation with antidote.

5.2    A time dependent modi cation

Suppose that we modify our equation for Z to be of the form

dZ
= ZS+ (N0
S  Z)    ZS  K(S; t);
(3)
dt





where t is time measured in days. Suggest a reasonable function K(S; t) that models some time-dependent feature of the system. Discuss with your team how you can make the model more realistic or, come up with a strategy for more e ectively defeating Zombies. Make sure your function K explicitly depends on t! Be creative, there is no single correct answer here.

5

• specify your equation K(S; t) explicitly and describe what it means in real world terms.

    • Include plots of S, Z, and R versus time and show how the function K a ects the outcome.

    • Conclusion

Write a letter to the president that summarizes your discoveries in a few paragraphs. Was the antidote an e ective control strategy? Was your strategy e ective?

    • Report Guidelines

Your group will submit your project on D2L, in the appropriate dropbox (you can nd these under the \assessments" tab in D2L). Adhere to the following guidelines:

    • Do not put o nding a group, do this early. You should have a group set up within one week of the project assignment date.

    • Submit your project in pdf format. Contents of .zip les will not necessarily be graded. (Word documents not acceptable because equations are commonly jumbled around by D2L.)

    • Submit ALL code used for your project (.nb  les for Mathematica, .m  les for MatLab, etc).

    • Have only ONE group member submit the project. Having multiple people in your group submit the project to D2L will result in multiple grades, and we will take the LOWEST one.

    • Include the names of all group members working on the project.

Your report needs to accurately and consistently describe the steps you took in answering the questions asked. This report should have the look and feel of a technical paper. Presentation and clarity are very important. Here are some important items to remember:

    • Absolutely make sure your recitation number is on your submitted report.

    • Start with an introduction that describes what you will discuss in the body of your document. A brief summary of important concepts used in your discussion could be useful here as well.

    • Summarize what you have accomplished in a conclusion. No new information or new results should appear in your conclusion. You should only review the highlights of what you wrote about in the body of the report.

    • Always include units in your answers.

    • Always label plots and refer to them in the text. The main body of your paper should NOT include lengthy calculations. These should be included in an appendix, and referred to in the main body.

    • Labs must be typed. Including the equations in the main body (part of your learning ex-perience is to learn how to use an equation editor). An exception can be made for lengthy calculations in the appendix, which can be hand written (as long as they are neat and clear), and minor labels on plots, arrows in the text and a few subscripts.

    • Your report does not have to be long. You need quality, not quantity of work. Of course you cannot omit any important piece of information, but you need not add any extras.

    • DO NOT include printouts of computer software screens. This will be considered as garbage. You simply need to state which software you used in each step, and what it did for you.

6

    • You must include any plot that supports your conclusions or gives you insight in your inves-tigations.

    • Write your report in an organized and logical fashion. Section headers such as Introduction, Background, Problem Statement, Calculations, Results, Conclusion, Appendix, etc... are not mandatory, but are highly recommended. They not only help you write your report, but help the reader navigate through your paper, besides giving it a clearer look.


























































7