# Simulating Disease Systems: Examples

This R script will simulate several different disease models. We will begin with a very simple model that tracks only the numbers of suscpeptible (uninfected) and infected people without considering recovery.

Since you will use the function **ode** in the **deSolve** package to do the numerical simulation, make sure that package is installed by running

`library(deSolve)`

If you get the message “Error in library(deSolve) : there is no package called ‘deSolve’” then you should run the following code, which installs the package from a repository on U of T’s servers (you will have to have an internet connection for this).

`install.packages('deSolve', repos='http://probability.ca/cran')`

To learn more about the **deSolve** package, and the **ode** function in particular, you can bring up the R documentation for the function by running either

```
help('ode') # or, alternatively
?ode
```

A pdf document that provides very detailed information about the package and examples of how to use it can be viewed by running

`vignette("deSolve")`

If you look through the help documentation for **deSolve**, you will see that you need to define your own function to calculate the rate of change for the state variables. In order to work with the built-in numerical integrators of **ode**, that function has to take specific arguments, in a specific order, and return the derivatives in a specific way. If any of these requirements are not met, **ode** will fail.

**ode** requires that the function that calculates the rate of change of state variables take three arguments. Arguments are the elements that are passed to the function (for example, for \(f(x) = x^2\), \(f\) is a function that takes an argument \(x\), \(f(x,y) = x^2 + y^2\) is a function that takes arguments \(x\) and \(y\), and so on). For **ode**, the function must be of the form **function(t, state, parameters)**, where **t** is the variable that is being integrated over (in this case, time), **state** is a vector of state variables, and **parameters** is a vector of parameters required by the model.

#### S-I model without recovery

Let’s consider the following simple system: \[ \begin{align} \frac{dS}{dt} &= b (S + I) - \beta S I - \mu S, \\ \frac{dI}{dt} &= \beta S I - (\mu + \nu) I. \end{align} \] In this model, \(b\) is the birth rate, and we have assumed that both suscpetible and infected individuals have the same birth rate. \(\beta\) is the transmission rate per infected person. \(\mu\) is the baseline mortality rate, the mortality rate due to causes other than infection, and \(\nu\) is the mortality rate caused by infection. Notice: in the absence of any infection (\(I = 0\)), the susceptible population will grow exponentially.

To encode these dynamics in a model, let’s begin by defining the parameter vector. I will create a “named vector” in R using the **c** function. This will allow me to reference the parameters by name so I don’t have to keep track of the order of the parameters in the vector.

```
parameters <- c(b=1, B=0.05, m=0.25, v=1)
parameters
```

```
## b B m v
## 1.00 0.05 0.25 1.00
```

We can also define a vector holding the initial conditions. Again, we can create a named vector to make keeping track of things easier. We will assume that we have a population where there are many susceptible individuals, but only a single infected individual.

```
initial <- c(S=100, I=1)
initial
```

```
## S I
## 100 1
```

Now let’s define the dynamical system. We will give the function an informative name so we can easily remember what it does. Inside of the function, we will reference the parameter vector and assign parameters to new variables that will exist only inside the function. For simplicity of expression, I will also assign the elements of the **state** vector to variables.

```
si_model <- function(t, state, parameters) {
## assign parameter values to parameter variables
## the function 'unname()' just removes the name of the parameter - this is unnecessary, it just cleans things up a bit
b <- unname(parameters['b'])
B <- unname(parameters['B'])
m <- unname(parameters['m'])
v <- unname(parameters['v'])
## assign state variables to names
S <- unname(state['S'])
I <- unname(state['I'])
## compute the rates of change
dSdt <- b * (S + I) - B * S * I - m * S
dIdt <- B * S * I - (m + v) * I
## return as a list object, with the first element of the list being the derivatives. The order of derivatives must be the same as the order in the initial condition vector!
return(list(c(dSdt, dIdt)))
}
```

To simulate the model, we need to call **ode**, supplying the proper information to each of the arguments of that function. **ode** is a function taking five arguments: **y**, the vector of initial conditions; **times** a vector indicating when output should be recorded; **func**, the name of function that will compute the derivatives; **parms**, a vector of parameters; and **method**, which determines which numerical integration routine to use. For this class, you will almost certainly only need the default method, so we do not need to supply anything for **method**.

```
## When should output be recorded?
times <- seq(0, 100, 0.1)
## Simulate the system!
output <- ode(y=initial, times=times, func=si_model, parms=parameters)
## Examine the last few rows of data
tail(output)
```

```
## time S I
## [996,] 99.5 25 75
## [997,] 99.6 25 75
## [998,] 99.7 25 75
## [999,] 99.8 25 75
## [1000,] 99.9 25 75
## [1001,] 100.0 25 75
```

We can use plotting to visualize the dynamics of the system. In particular, let’s plot both variables on the same axis, so we can easily compare the two. You see the initial outbreak of disease when the majority of the population is susceptible, followed by a decline towards an endemic equilibrium where both susceptible and infected individuals coexist.

```
## the expression 'output[,"time"]' tells R to show you all of the data contained in the column called 'time'
## 'lwd' specifies the line width
## ylim specifies the ylimits to be the maximum and minimum values of S and I
plot(output[,"time"], output[,"S"], type='l', xlab='Time', ylab='Number of individuals', lwd=2, ylim=range(output[,c("S","I")]))
## add infected individuals to the plot as a red line
lines(output[,"time"], output[,"I"], col="red", lwd=2)
## add an informative legend to the top right corner of the plot
legend(x='topright', c("Susceptible", "Infected"), fill=c("black","red"), bty="n")
```

Of course, problems can arise. For example, if you try to run the code below, you will get an error message “Error in checkInput(y, times, …) : argument ‘func’ is missing, with no default,” because you forgot to tell **ode** what function to use to compute the derivatives.

`output <- ode(y=initial, times=times, parms=parameters)`

The following code will also generate the error “Error in checkInput(y,…) : ‘func’ must be a function or character vector.” This error is caused because you didn’t pass the arguments to **ode** in the correct order: **ode** is expecting that the order of arguments is **(y, times, func, parameters)**. If you do not explicity specify what you are assigning to each argument using an equals sign (that is, instead of typing **ode(y= my_y, times=my_t, …)** you just type

**ode(**), you

*my_y*,*my_t*, …)*must*put the arguments in the order the function is expecting.

`output <- ode(initial, parameters, times, si_model)`

For purposes of comparison, the following code will run just fine, even though the arguments are not in the order **ode** is expecting, because we have explicitly assigned everything to the correct function argument:

```
output <- ode(func=si_model, times=times, y=initial, parms=parameters)
tail(output)
```

```
## time S I
## [996,] 99.5 25 75
## [997,] 99.6 25 75
## [998,] 99.7 25 75
## [999,] 99.8 25 75
## [1000,] 99.9 25 75
## [1001,] 100.0 25 75
```

The following code has the same problem: the arguments are not given in the correct order. If you run this, you will get the error “The number of derivatives returned by func() (2) must equal the length of the initial conditions vector (1001).” This is because you passed ‘times’ to the argument ‘y’, and ‘initial’ to the argument ‘times’.

`output <- ode(times, initial, si_model, parameters)`

Another common problem is that the order of the variables in the initial conditions must match the order of the derivatives in the **return** statement inside your derivative function. For example, if you had put the initial conditions in the wrong order, the results would look very different:

```
wrong_initial <- c(I=1, S=100)
output <- ode(y=wrong_initial, times=times, func=si_model, parms=parameters)
tail(output)
```

```
## time I S
## [996,] 99.5 15 2.094e+34
## [997,] 99.6 15 2.257e+34
## [998,] 99.7 15 2.432e+34
## [999,] 99.8 15 2.622e+34
## [1000,] 99.9 15 2.826e+34
## [1001,] 100.0 15 3.046e+34
```

In general, these errors can be avoided by always being *very* careful. A more challenging problem to diagnose arises when the numerical integration works and generates a numerical solution, but one that seems incorrect. I will give a simple example here. Imagine that you want to explore what happens as you change the mortality rate caused by the disease, so you modify the parameter vector, lowering the value of \(\nu\), the mortality rate caused by the disease.

```
parameters2 <- c(b=1, B=0.05, m=0.25, v=0.5)
output2 <- ode(y=initial, times=times, func=si_model, parms=parameters2)
```

This code runs just fine, but if you look at the final values the system reaches, you notice that the value of the susceptible population size is unchanging, and reasonable (\(S = 20\)), but the infected population is huge, and still increasing (\(I = 1.34 \times 10^{13}\))! This seems pretty strange - what is going on?

`tail(output2)`

```
## time S I
## [996,] 99.5 20 1.180e+13
## [997,] 99.6 20 1.210e+13
## [998,] 99.7 20 1.240e+13
## [999,] 99.8 20 1.272e+13
## [1000,] 99.9 20 1.304e+13
## [1001,] 100.0 20 1.337e+13
```

Remember that I pointed out that, in the absence of infection, the susceptible host population would grow exponentially. In an exponentially growing population, if the growth rate is positive, the susceptible population will grow towards infinity and there is no equilibrium except 0. So, there are two possibilities when we are considering a population of both susceptible and infected individuals: either infection *will* regulate total population growth, leading to an equilibrium (as we saw initially), or it will not, and the total population size will continue to grow exponentially. That is what we are seeing here: the total population (\(S + I\)) is growing exponentially, and the huge number of infected individuals means that everyone gets infected, but not right away: the influx of new births keeps some number of susceptible hosts around, who quickly become infected.

A bit of analysis would have indicated to us that setting the mortality rate caused by infection too low could be problematic. In particular, we could have found the equilibria of this model by setting both \(\frac{dS}{dt} = 0\) and \(\frac{dI}{dt} = 0\) and solving the resulting set of algebraic equations for \(\hat{S}\) and \(\hat{I}\). In particular, setting \(\frac{dI}{dt} = 0\), you can immediately solve for \(\hat{S}\), the equilibrium number of susceptible hosts: \[ \hat{S} = \frac{\mu + \nu}{\beta}.\] Setting \(\frac{dS}{dt} = 0\) and plugging in this value for \(\hat{S}\), you can arrive at an expression for \(\hat{I}\): \[ \hat{I} = \frac{(b-\mu) \hat{S}}{\beta \hat{S}-b} = \frac{(b-\mu)(\mu+\nu)}{\beta (\mu + \nu - b)}. \] Note that \(\hat{S}\) will always be positive, but \(\hat{I}\) will only be positive if \(b - \mu\) and \(\mu + \nu - b\) have the same sign. Plugging in the parameter values we were just considering, \(b - \mu = 1 - 0.25 = 0.75\) and \(\mu + \nu - b = 0.25 + 0.5 - 1 = -0.25\). Because these have opposite signs, \(\hat{I} < 0\)! Obviously this equilibrium is biologically unfeasible and dynamically unattainable (because we are working in continuous time, to get \(I < 0\), \(I\) must pass through 0, but when \(I = 0\), \(\frac{dI}{dt} = 0\), meaning that \(I\) can never become negative if it starts out positive.) This suggests that there will be no equilibrium of the model, which is exactly what we just observed in the numerical simulation results.

Even with very complicated models, a small bit of analysis can sometimes go a long way towards helping understand your numerical results.

#### S-I model with recovery, but no immunity

How would we need to change the model in order to allow individuals to recover from infection, but not to be immune to getting reinfected? A simple model for this case would be \[ \begin{align} \frac{dS}{dt} &= b (S + I) - \beta S I - \mu S + \gamma I, \\ \frac{dI}{dt} &= \beta S I - (\mu + \nu + \gamma) I. \end{align} \] This model is identical to the one above, except we have added a new parameter \(\gamma\) that quantifies the rate of recovery from infection. I need to modify the parameter vector, but the initial conditions stay the same.

```
parameters_2 <- c(b=1, B=0.05, m=0.25, v=1, g=0.25)
initial <- c(S=100, I=1)
```

I also need to modify the function encoding the dynamical system:

```
si_model_2 <- function(t, state, parameters) {
## assign parameter values to parameter variables
## the function 'unname()' just removes the name of the parameter - this is unnecessary, it just cleans things up a bit
b <- unname(parameters['b'])
B <- unname(parameters['B'])
m <- unname(parameters['m'])
v <- unname(parameters['v'])
g <- unname(parameters['g']) ## add recovery rate as a parameter
## assign state variables to names
S <- unname(state['S'])
I <- unname(state['I'])
## compute the rates of change
dSdt <- b * (S + I) - B * S * I - m * S + g * I
dIdt <- B * S * I - (m + v + g) * I
## return as a list object, with the first element of the list being the derivatives. The order of derivatives must be the same as the order in the initial condition vector!
return(list(c(dSdt, dIdt)))
}
```

Now I can numerically simulate this new model using the same code as above. However, what if you made a mistake and passed the wrong parameter vector to this new model? Remember that the parameter vector used for the first model was called ‘parameters’ and this new parameter vector is called ‘parameters_2’. Notice that you do not get any error messages, but when you inspect ‘output_2’ using the function **tail** to look at the last few rows of data, you see that both the \(S\) and \(I\) population sizes are recorded as NA. This is typically what happens when a parameter is missing! The numerical seems to run, but when you actually look at the results, everything is NA. If you get results like these, the first thing to check is to make sure that all of the parameters are defined properly. R is case-sensitive, so upper- versus lower-case names matter! This is probably the most common mistake.

```
output_2 <- ode(y=initial, times=times, func=si_model_2, parms=parameters)
tail(output_2)
```

```
## time S I
## [996,] 99.5 NA NA
## [997,] 99.6 NA NA
## [998,] 99.7 NA NA
## [999,] 99.8 NA NA
## [1000,] 99.9 NA NA
## [1001,] 100.0 NA NA
```

Let’s run the model with the correct parameter vector, and then plot the results along with the results of the first model so we can directly compare them. You can see that the epidemic peak (the maximum number of infected individuals) is the same for both models, but the model with recovery settles into an equilibrium with more susceptible *and* more infected hosts. On the one hand, this seems counterintuitive: shouldn’t increasing recovery *reduce* the amount of infection? The reason for this prediction is that epidemics tend to fade out, either to extinction or to some equilibrium, because the disease exhausts the suscpetible population. Recovery “refills” the susceptible pool, allowing both the susceptible and infected populations to be larger at equilibrium.

```
## Run the new model
output_2 <- ode(y=initial, times=times, func=si_model_2, parms=parameters_2)
## Run the old model
output <- ode(y=initial, times=times, func=si_model, parms=parameters)
## Plot the number of susceptibles in the old model as a black line
plot(output[,"time"], output[,"S"], type='l', lwd=2, xlab='Time', ylab='Number of individuals', ylim=range( c( output[,2:3], output_2[,2:3] ) )) ## set ylim so that the none of the lines are cut off
## Plot the number of infecteds in the old model as a red line
lines(output[,'time'], output[,"I"], col="red", lwd=2)
## Plot the number of susceptibles in the new model as a dashed black line
lines(output_2[,"time"], output_2[,"S"], lwd=2, lty=2)
## Plot the number of infecteds in the new model as a dashed red line
lines(output_2[,"time"], output_2[,"I"], lwd=2, lty=2, col="red")
## Add a legend
legend(x='topright', c('S: old', 'S: new', 'I: old', 'I: new'), lty=c(1,2,1,2), lwd=2, col=c("black","black","red","red"))
```

#### S-I-R model with recovery to an immune class

For this model, we assume that individuals that recover from infection cannot get sick again - they are immune to the disease. The model for this scenario would look like \[ \begin{align} \frac{dS}{dt} &= b (S + I + R) - \beta S I - \mu S, \\ \frac{dI}{dt} &= \beta S I - (\mu + \nu + \gamma) I,\\ \frac{dR}{dt} &= \gamma I - \mu R \end{align} \] We have added a new state variable, \(R\), which is the total number of immune hosts. These hosts give birth to new susceptible hosts at a rate \(b\), just like susceptible and infected hosts, and they suffer a constant background mortality rate of \(\mu\).

Before we numerically simulate this system, let’s take a moment to think about what we might expect the dynamics to look like. If the recovery rate is very low (much lower than the background or disease-induced mortality rate, for example), then the dynamics will look like those of the first model. However, if the recovery rate is high, infection may no longer be able to regulate the host population size, and the population will grow exponentially.

Let’s see if our intuition is correct. I wills start with the a parameter vector where the recovery rate is very low. I also need to modify the initial conditions to include immune hosts. We will assume that the initial number of immune hosts is 0.

```
## low recovery rate
parameters_sir <- c(b=1, B=0.05, m=0.25, v=1, g=0.01)
## high recovery rate
parameters_sir_2 <- c(b=1, B=0.05, m=0.25, v=1, g=0.25)
initial_sir <- c(S=100, I=1, R=0)
```

I also need to modify the function encoding the dynamical system:

```
sir_model <- function(t, state, parameters) {
## assign parameter values to parameter variables
## the function 'unname()' just removes the name of the parameter - this is unnecessary, it just cleans things up a bit
b <- unname(parameters['b'])
B <- unname(parameters['B'])
m <- unname(parameters['m'])
v <- unname(parameters['v'])
g <- unname(parameters['g']) ## add recovery rate as a parameter
## assign state variables to names
S <- unname(state['S'])
I <- unname(state['I'])
R <- unname(state['R']) ## add in the new state variable
## compute the rates of change
dSdt <- b * (S + I + R) - B * S * I - m * S
dIdt <- B * S * I - (m + v + g) * I
dRdt <- g * I - m * R ## and keep track of its dynamics
## return as a list object, with the first element of the list being the derivatives. The order of derivatives must be the same as the order in the initial condition vector!
return(list(c(dSdt, dIdt, dRdt)))
}
```

Simulating the model with the recovery rate set to a very low value, we can see that the final dynamical state of the system is similar to that of the first model. There is a slightly higher number of infected individuals because recovered individuals contribute to population growth, increasing the size of the susceptible class upon which the infection depends.

```
output_sir <- ode(y=initial_sir, times=times, func=sir_model, parms=parameters_sir)
## compare the numerical simulation results for this model
tail(output_sir)
```

```
## time S I R
## [996,] 99.5 25.2 85.91 3.436
## [997,] 99.6 25.2 85.91 3.436
## [998,] 99.7 25.2 85.91 3.436
## [999,] 99.8 25.2 85.91 3.436
## [1000,] 99.9 25.2 85.91 3.436
## [1001,] 100.0 25.2 85.91 3.436
```

```
## with those of the first model
tail(output)
```

```
## time S I
## [996,] 99.5 25 75
## [997,] 99.6 25 75
## [998,] 99.7 25 75
## [999,] 99.8 25 75
## [1000,] 99.9 25 75
## [1001,] 100.0 25 75
```

Simulating the model with the recovery rate set much higher, you see that the dynamics are very different: the population is growing exponentially, leading to exponentially increasing infected and recovered classes.

```
output_sir_2 <- ode(y=initial_sir, times=times, func=sir_model, parms=parameters_sir_2)
## compare the numerical simulation results for this model
tail(output_sir)
```

```
## time S I R
## [996,] 99.5 25.2 85.91 3.436
## [997,] 99.6 25.2 85.91 3.436
## [998,] 99.7 25.2 85.91 3.436
## [999,] 99.8 25.2 85.91 3.436
## [1000,] 99.9 25.2 85.91 3.436
## [1001,] 100.0 25.2 85.91 3.436
```

#### Model of SARS

As a final example, I will show how to encode the SARS model you saw earlier in the semester. This model comes from Gumel, A. B. et al. (2004) Modelling strategies for controlling SARS outbreaks. *Proceedings of the Royal Society B* **271** 2223-2232. The paper and appendix can be found by following the links below. You will have be on campus to access the appendix. Paper Appendix

The original model consists of six classes of individuals: susceptible (\(S(t)\)), asymptomatic (\(E(t)\)), quarantined (\(Q(t)\)), symptomatic (\(I(t)\)), isolated (\(J(t)\)), and recovered individuals in a population of \(N = S(t) + E(t) + Q(t) + I(t) + J(t) + R(t)\) individuals. I am going to augment the model to also track the number of people dying of SARS (\(D(t)\)).

The population of susceptibles increases through both birth and migration, at the rate \(\Pi\). The population is reduced through natural mortality at the rate \(\mu\). Susceptible hosts become infected with SARS by interacting with any of the four classes of infected individuals. The parameter \(\beta\) is the basic transmission rate, which captures both the rate of contact between susceptible and infected hosts and the probability of transmission upon contact. Asymptomatic, quarantined, and isolated individuals may also transmit SARS, but at a modified rate, as captured by the parameters \(\epsilon_E\), \(\epsilon_Q\), and \(\epsilon_J\). If these parameters are all less than one, then transmission from any of these classes is less than transmission from the infected class, as you might expect. The authors assume frequency-dependent (rather than density-dependent) transmission, so the entire transmission rate is divided by the total number of individuals in the population. The dynamics of the susceptible population are given by: \[ \frac{dS}{dt} = \Pi - \frac{S \left(\beta I + \epsilon_E \beta E + \epsilon_Q \beta Q + \epsilon_J \beta J \right)}{N} - \mu S.\]

Asymptomatic individuals are those who are carrying (and possibly transmitting) SARS, but who do not yet show any symptoms of infection. New asymptomatic individuals may arrive into the population through migration at a rate \(p\), but are mainly increasing due to transmission. The asymptomatic population decreases due to quarantine at a rate \(\gamma_1\), development of symptoms at a rate \(\kappa_1\), and natural death. The dynamics are thus given by: \[ \frac{dE}{dt} = p + \frac{S \left(\beta I + \epsilon_E \beta E + \epsilon_Q \beta Q + \epsilon_J \beta J \right)}{N} - \left(\gamma_1 + \kappa_1 + \mu\right) E.\]

Asymptomatic individuals may also be quarantined. The model assumes that individuals who are actually infected will be quaratined (there is no possibility for accidental quarantine of susceptible individuals). There is no migration of quarantined individuals, but quarantined individuals may leave the class due to development of symptoms at a rate \(\kappa_2\) and natural mortality \(\mu\). The dynamics are therefore \[ \frac{dQ}{dt} = \gamma_1 E - \left(\kappa_2 + \mu \right) Q.\]

The number of symptomatic individuals increases due to disease progression from asymptomatic individuals at a rate \(\kappa_1\). Symptomatic individuals may be isolated at a rate \(\gamma_2\), they may recover from infection at a rate \(\sigma_1\), or they may die at the elevated rate \(d_1 + \mu\) due to infection-induced mortality and natural mortality. The dynamics of symptomatic individuals are then \[ \frac{dI}{dt} = \kappa_1 E - \left(\gamma_2 + \sigma_1 + d_1 + \mu \right) I. \]

Isolated individuals have developed clinical symptoms but have been isolated from the general population. The number of isolated individuals increases through the isolation of symptomatic individuals and the progression of asymptomatic quarantined individuals. The number decreases due to death at the rate \(d_2 + \mu\), where \(d_2 < d_1\) because isolated individuals have a better level of care and may receive partially-effective treatment. Some individuals also recover to an immune class at a rate \(\sigma_2\). The dynamics are then \[ \frac{dJ}{dt} = \gamma_2 I + \kappa_2 Q - \left(\sigma_2 + d_2 + \mu \right) J. \]

The number of recovered individuals increases through recovery of infected and isolated individuals and decreases due to background mortality: \[ \frac{dR}{dt} = \sigma_1 I + \sigma_2 J - \mu R. \]

Finally, the number of people dying of SARS is just: \[ \frac{dD}{dt} = d_1 I + d_2 J.\]

This model can be coded into R, with no small amount of difficulty, as

```
sars_model <- function(t, state, parameters) {
P <- unname(parameters['P']) # growth of susceptible population
B <- unname(parameters['B']) # transmission rate
e_E <- unname(parameters['e_E']) # relative transmission from asymptomatic individuals
e_Q <- unname(parameters['e_Q']) # relative transmission from quarantined individuals
e_J <- unname(parameters['e_J']) # relative transmission from isolated individuals
m <- unname(parameters['m']) # background mortality rate
p <- unname(parameters['p']) # migration rate of asymptomatic individuals
g_1 <- unname(parameters['g_1']) # rate of quarantine
K_1 <- unname(parameters['K_1']) # rate of symptom development for asymptomatic individuals
K_2 <- unname(parameters['K_2']) # rate of symptom development for quarantined individuals
g_2 <- unname(parameters['g_2']) # rate of isolation
s_1 <- unname(parameters['s_1']) # rate of recovery for symptomatic individuals
d_1 <- unname(parameters['s_1']) # rate of SARS-death for symptomatic individuals
s_2 <- unname(parameters['s_2']) # rate of recovery for isolated individuals
d_2 <- unname(parameters['d_2']) # rate of SARS-death for isolated individuals
S <- unname(state['S'])
E <- unname(state['E'])
Q <- unname(state['Q'])
I <- unname(state['I'])
J <- unname(state['J'])
R <- unname(state['R'])
D <- unname(state['D'])
# total population size
N <- S + E + Q + I + J + R
dSdt <- P - S * B * (I + e_E * E + e_Q * Q + e_J * J) / N - m * S
dEdt <- p + S * B * (I + e_E * E + e_Q * Q + e_J * J) / N - (g_1 + K_1 + m) * E
dQdt <- g_1 * E - (K_2 + m) * Q
dIdt <- K_1 * E - (g_2 + s_1 + d_1 + m) * I
dJdt <- g_2 * I + K_2 * Q - (s_2 + d_2 + m) * J
dRdt <- s_1 * I + s_2 * J - m * R
dDdt <- d_1 * I + d_2 * J
list(c(dSdt,dEdt,dQdt,dIdt,dJdt,dRdt,dDdt))
}
```

To evaluate the effect of different types of quaratine, we will simulate the system for different values of the parameters for the SARS outbreak in the Greater Toronto Area. These parameter values are taken from the paper and online appendix to the original paper, which are posted on the course website. **Especially if you will need to estimate model parameters on the basis of empirical data for your project, read the Appendix!** I will do three simulations: a simulation without any quaratine or isolation, a simulation with only quaratine, and a simulation with only isolation, to see how each affects the dynamics of SARS mortality over time. The solid line shows the cumulative number of SARS deaths when no quarantine or isolation is used; the dashed line shows the same when quarantine but not isolation, and the dotted line shows the case of isolation but no quaratine.

```
# take special care to make sure that the initial conditions are in the same order as the derivative in the function - this is an easy mistake to make when coding large models
# these are the parameters for the GTA
initial_sars <- c(S=4000000, E=6, Q=0, I=1, J=0, R=0, D=0)
times <- seq(0, 100, 0.1)
## parameters with no quarantine or isolation
parameters_sars_1 <- c(P=136, p=0.06, m=0.000034, K_1=0.1, K_2=0.125, s_1=0.0337, s_2=0.0386, d_1=0.0079, d_2=0.0068, g_1=0, g_2=0, B=0.2, e_J=0.36, e_E=0, e_Q=0)
## simulate the model
out_1 <- ode(y=initial_sars, times=times, func=sars_model, parms=parameters_sars_1)
## parameters with quarantine but no isolation
parameters_sars_2 <- c(P=136, p=0.06, m=0.000034, K_1=0.1, K_2=0.125, s_1=0.0337, s_2=0.0386, d_1=0.0079, d_2=0.0068, g_1=0.2, g_2=0, B=0.2, e_J=0.36, e_E=0, e_Q=0)
## simulate the model
out_2 <- ode(y=initial_sars, times=times, func=sars_model, parms=parameters_sars_2)
## parameters with isolation but no quarantine
parameters_sars_3 <- c(P=136, p=0.06, m=0.000034, K_1=0.1, K_2=0.125, s_1=0.0337, s_2=0.0386, d_1=0.0079, d_2=0.0068, g_1=0, g_2=0.2, B=0.2, e_J=0.36, e_E=0, e_Q=0)
## simulate the model
out_3 <- ode(y=initial_sars, times=times, func=sars_model, parms=parameters_sars_3)
## Plot the cumulative number of SARS deaths for each case.
## Use lty=1 to make the first line solid, lty=2 to make the second line dashed, and lty=3 to make the line dotted
plot(out_1[,c('time','D')], type='l', lty=1, xlab='Days since initial infection', ylab='Cumulative number of SARS deaths')
lines(out_2[,c('time','D')], lty=2)
lines(out_3[,c('time','D')], lty=3)
```