# Etivity 1.1.1: Modelling an infected cohort

In this etivity, you will code your first simple model and use it to answer a research question. To get used to the format for models in R, we'll start by giving you some parts of the code, and you need to fill in the missing numbers or code blocks by replacing the *#YOUR CODE#* placeholders. At the end of each etivity, we'll provide a solution file that contains the full code, annotated with comments. Make sure you understand each line, and soon you'll be able to write your own model from scratch!

Your task is to find out how long it takes for
a cohort of infected people to recover. As you saw in the video,
to answer this question you need to keep track of 2 populations: 
those that are infected (compartment I), and those that have recovered (compartment R). 
Infected people recover at a rate *gamma*. The differential equations describing this are:

\begin{align}
\frac{dI}{dt} & = -\gamma I \\
\frac{dR}{dt} & = \gamma I
\end{align}

Load the packages which you need for this etivity, by running the following cell:


```
library(deSolve)   # package to solve the model
library(reshape2)  # package to change the shape of the model output
library(ggplot2)   # package for plotting
```

To start, it is useful to code what we know about the situation we want to model.
We are looking at a cohort of 10$^6$ currently infected people, and no one has recovered so far. The average duration of infection is 10 days. The question we want to answer is how many people will recover from the infection over a 4-week period.

Given this data, fill in the correct values to the following variables, and run the cell:


```
initial_number_infected <- 1000000  # the initial infected population size
initial_number_recovered <- 0#YOUR CODE#  # the initial number of people in the recovered state
recovery_rate <- 0.1            # the rate of recovery gamma, in units of days^-1
follow_up_duration <- 4*7 #YOUR CODE#        # the duration to run the model for, in units of days
  
# Hint: the units of the recovery rate and the follow-up duration should be consistent.
```

Now, we combine this data into objects that are recognised by the deSolve package as **model input**. To do this, again run the code below.


```
# The initial state values are stored as a vector and each value is assigned a name.
initial_state_values <- c(I = initial_number_infected, 
                          R = initial_number_recovered)

# Parameters are also stored as a vector with assigned names and values. 
parameters <- c(gamma = recovery_rate)  
# In this case we only have one parameter, gamma.
```

Think about: what kind of information is stored in the *initial_state_values* and *parameters* vectors?

Additionally, we need to specify the time we want the model to run for. This depends on the question we want to answer. In the cell below, the duration you specified earlier is automatically filled in when you run it.


```
# The times vector creates a sequence of timepoints at which we want to calculate the number of people 
# in the I and R compartment.
times <- seq(from = 0, to = follow_up_duration, by = 1) 
```

Think about: what kind of information is stored in the *times* vector?

Check your answers by having a look at each of these vectors to familiarise yourself with the structure: in the following code cell we have typed the object names, so you just need to press "Run" to see what each of them contains.


```
initial_state_values
parameters
times
```


<dl class=dl-horizontal>
	<dt>I</dt>
		<dd>1e+06</dd>
	<dt>R</dt>
		<dd>0</dd>
</dl>




<strong>gamma:</strong> 0.1



<ol class=list-inline>
	<li>0</li>
	<li>1</li>
	<li>2</li>
	<li>3</li>
	<li>4</li>
	<li>5</li>
	<li>6</li>
	<li>7</li>
	<li>8</li>
	<li>9</li>
	<li>10</li>
	<li>11</li>
	<li>12</li>
	<li>13</li>
	<li>14</li>
	<li>15</li>
	<li>16</li>
	<li>17</li>
	<li>18</li>
	<li>19</li>
	<li>20</li>
	<li>21</li>
	<li>22</li>
	<li>23</li>
	<li>24</li>
	<li>25</li>
	<li>26</li>
	<li>27</li>
	<li>28</li>
</ol>



The next step is specifying the model. Using the example code in the introductory document, complete the following model function with the differential equations above.


```
cohort_model <- function(time, state, parameters) { 
  
    with(as.list(c(state, parameters)), {
      
      dI <- -gamma*I
      dR <- gamma*I
      
      return(list(c(dI, dR)))
    })
  
}
```

Now all there's left to do is solving this set of equations using the deSolve package. Fill in the following command, which calculates and stores the number of infected and recovered people at each timestep in the "output" dataframe. Don't forget to run it!


```
# Hint: if you can't remember what those arguments correspond to, just look up the ode help file:
??ode
```


```
output <- as.data.frame(ode(y = initial_state_values, 
                            times = times, 
                            func = cohort_model,
                            parms = parameters))
```

Printing the **model output** returns a dataframe with columns time (containing our times vector), I (containing the number of infected people at each timestep) and R (containing the number of recovered people at each timestep):


```
output
```


<table>
<thead><tr><th scope=col>time</th><th scope=col>I</th><th scope=col>R</th></tr></thead>
<tbody>
	<tr><td> 0          </td><td>1.000000e+06</td><td>      0.0   </td></tr>
	<tr><td> 1          </td><td>3.678794e+05</td><td> 632120.6   </td></tr>
	<tr><td> 2          </td><td>1.353351e+05</td><td> 864664.9   </td></tr>
	<tr><td> 3          </td><td>4.978697e+04</td><td> 950213.0   </td></tr>
	<tr><td> 4          </td><td>1.831559e+04</td><td> 981684.4   </td></tr>
	<tr><td> 5          </td><td>6.737923e+03</td><td> 993262.1   </td></tr>
	<tr><td> 6          </td><td>2.478741e+03</td><td> 997521.3   </td></tr>
	<tr><td> 7          </td><td>9.118772e+02</td><td> 999088.1   </td></tr>
	<tr><td> 8          </td><td>3.354606e+02</td><td> 999664.5   </td></tr>
	<tr><td> 9          </td><td>1.234090e+02</td><td> 999876.6   </td></tr>
	<tr><td>10          </td><td>4.539958e+01</td><td> 999954.6   </td></tr>
	<tr><td>11          </td><td>1.670156e+01</td><td> 999983.3   </td></tr>
	<tr><td>12          </td><td>6.144155e+00</td><td> 999993.9   </td></tr>
	<tr><td>13          </td><td>2.260306e+00</td><td> 999997.7   </td></tr>
	<tr><td>14          </td><td>8.315196e-01</td><td> 999999.2   </td></tr>
	<tr><td>15          </td><td>3.058986e-01</td><td> 999999.7   </td></tr>
	<tr><td>16          </td><td>1.125336e-01</td><td> 999999.9   </td></tr>
	<tr><td>17          </td><td>4.139872e-02</td><td>1000000.0   </td></tr>
	<tr><td>18          </td><td>1.522972e-02</td><td>1000000.0   </td></tr>
	<tr><td>19          </td><td>5.602676e-03</td><td>1000000.0   </td></tr>
	<tr><td>20          </td><td>2.061173e-03</td><td>1000000.0   </td></tr>
	<tr><td>21          </td><td>7.582818e-04</td><td>1000000.0   </td></tr>
	<tr><td>22          </td><td>2.789371e-04</td><td>1000000.0   </td></tr>
	<tr><td>23          </td><td>1.026131e-04</td><td>1000000.0   </td></tr>
	<tr><td>24          </td><td>3.773270e-05</td><td>1000000.0   </td></tr>
	<tr><td>25          </td><td>1.391230e-05</td><td>1000000.0   </td></tr>
	<tr><td>26          </td><td>5.151106e-06</td><td>1000000.0   </td></tr>
	<tr><td>27          </td><td>1.851387e-06</td><td>1000000.0   </td></tr>
	<tr><td>28          </td><td>6.704086e-07</td><td>1000000.0   </td></tr>
</tbody>
</table>



### Question: Based on the output, how many people have recovered after 4 weeks? What proportion of the total population does this correspond to?

Now, plot your model output in the following cell, with time on the x axis and the number of infected and recovered people on the y axis.  You can use the introductory document for help with this.


```
# First turn the output dataset into a long format, so that the number in each compartment at each timestep
# are all in the same column
output_long <- melt(as.data.frame(output), id = "time")                  

# Plot the number of people in each compartment over time
ggplot(data = output_long,                                               # specify object containing data to plot
        aes(x = time, y = value, colour = variable, group = variable))+  
       # in the long-format output dataset, the number in each compartment is in the "value" column and the 
       # compartment the number relates to is in the "variable" column. We are telling ggplot to assign a different 
       # colour to each compartment/group, which automatically generates a legend
  geom_line() +                                                          
       # we want to represent the data over time as lines. This command automatically looks to the specifications saved
       # in the ggplot command above to know which data to plot
  xlab("Time (days)")+                                                   # add label for x axis
  ylab("Number of people") +                                             # add label for y axis
  labs(title = paste("Number infected and recovered over time when gamma =",parameters["gamma"],"days^-1")) # add title
# Using the paste() command, we can combine sentences with the values stored in variables (here gamma)
```

### Question: Based on the plot, at what timepoint were infected and recovered individuals equal in number?

For the last part of the etivity, try varying *gamma* to see how it affects the output. For example, in the cell below change *gamma* to correspond to an average infectious period of 2 days and 20 days. What is the recovery rate in these 2 cases?

7 days


```
parameters <- c(gamma = 1.0)
output <- as.data.frame(ode(y = initial_state_values, 
                            times = times, 
                            func = cohort_model,
                            parms = parameters))

output_long <- melt(as.data.frame(output), id = "time")                  # turn output dataset into long format

ggplot(data = output_long,                                               # specify object containing data to plot
       aes(x = time, y = value, colour = variable, group = variable)) +  # assign columns to axes and groups
  geom_line() +                                                          # represent data as lines
  xlab("Time (days)")+                                                   # add label for x axis
  ylab("Number of people") +                                             # add label for y axis
  labs(title = paste("Number infected and recovered over time when gamma =",parameters["gamma"],"days^-1")) # add title              
# Since the initial number in each compartment and the timesteps of interest haven't changed, 
# these are the only parts of the code we need to rerun.

# Now, copy-paste your plot code from above here to visualise the output.
```

### Question: What changes do you observe in the transition to the recovered compartment if *gamma* is higher or lower? For example, how long does it take for everyone to recover in both cases?

higher gamma, faster recovery time

## Well done on writing your first model code. Now, check the solutions!
Once you are done with an etivity, you can open the "Solutions to etivities" folder in your Jupyter workspace. In this folder, you will find feedback and model code for all the coding etivities throughout the course. After completing an exercise, you should always carefully compare your answers and code with the solutions! It is especially important to get the coding right from the start, because we are building on the same modelling framework throughout the course. Additionally, the solution files sometimes contain further information that will help you deepen your understanding of the lesson.  
To make the most of the learning experience, we always recommend trying to complete the whole etivity before checking the solutions - but if you are stuck at any point, they can also help you to move on to the next part.


```

```

2f737c25fa23fdee9f90ca16369fe177a02388dc
ipynb
Jupyter Notebook
course-1/utf-8''IDM1.1.07 - Modelling an infected cohort.ipynb
RKiddle/IDM
f32e0038df97f67acc47097827ed1430d8e13478
MIT
course-1/utf-8''IDM1.1.07 - Modelling an infected cohort.ipynb
RKiddle/IDM
f32e0038df97f67acc47097827ed1430d8e13478
MIT
course-1/utf-8''IDM1.1.07 - Modelling an infected cohort.ipynb
RKiddle/IDM
f32e0038df97f67acc47097827ed1430d8e13478
MIT
Qwen/Qwen-72B
1. YES
2. YES