ID Patterns - Practice

In DSAIDE: Dynamical Systems Approach to Infectious Disease Epidemiology (Ecology/Evolution)

#load various variable definitions that are the same for each app
source('startup_script.R')
currentrmdfile = knitr::current_input() 
appsettings = get_settings(currentrmdfile,appdocdir,packagename)

Overview {#shinytab1}

This app allows you to explore a model that tracks the same compartments as the Characteristics of ID model. If you haven't yet explored the Characteristics of ID model, I suggest you try out that one first. The model for this app adds a few more processes. It includes natural births and deaths of hosts, seasonal variation in transmission, and waning immunity.

Read about the model in the "Model" tab. Then do the tasks described in the "What to do" tab.

The Model {#shinytab2}

Model Overview

This model has the same compartments as the Characteristics of ID model:

• S - uninfected and susceptible individuals
• P - individuals who are infected and do not yet show symptoms. Those individuals can potentially be infectious
• A - individuals who are infected and do not show symptoms. Those individuals can potentially be infectious
• I - individuals who are infected and show symptoms. Those individuals are likely infectious, but the model allows to adjust this, including no infectiousness.
• R - recovered/removed individuals. Those individuals have recovered and are immune. They can lose their immunity in this model.
• D - individuals who have died due to the disease.

We include the following processes in this model:

• Susceptible individuals (S) can become infected by pre-symptomatic (P), asymptomatic (A) or symptomatic (I) hosts. The rates at which infections from the different types of infected individuals (P, A and I) occur are governed by 3 parameters, b~P~, b~A~, and b~I~.
• All infected individuals first enter the presymptomatic stage. They remain there for some time (determined by rate g~P~, the inverse of which is the average time spent in the presymptomatic stage). A fraction f of presymptomatic hosts move into the asymptomatic category, and the rest become symptomatic infected hosts.
• Asymptomatic infected hosts recover after some time (specified by the rate g~A~). Similarly, the rate g~I~ determines the duration the symptomatic hosts stay in the symptomatic state. For symptomatic hosts, two outcomes are possible. Either recovery or death. The parameter d determines the fraction of hosts that die.
• Recovered individuals are initially immune to reinfection. They can loose their immunity at rate w and return to the susceptible compartment.
• New susceptibles enter the system/model at a fixed rate n. From each compartment (apart from the dead compartment), hosts "leave" after some time (i.e. they die) at rate m. The inverse of this rate is the average lifespan of a host.
• The rates of transmission, b~P~, b~A~ and b~I~ can vary seasonally/annually, as described more below.

Note that we only track people that die due to the disease in our $D$ compartment. All hosts dying due to other causes just "exit the system" and we don't further keep track of them (though we could add another compartment to "collect" and track all individuals who died from non-disease-related causes.)

Model Implementation

The flow diagram and equations describe the model implemented in this app:

knitr::include_graphics(here::here('inst/media',appsettings$modelfigname))

$$b_P^s = b_P(1+s \sin(2\pi t / T))$$ $$b_A^s = b_A(1+s \sin(2\pi t /T))$$ $$b_I^s = b_I(1+s \sin(2\pi t /T))$$ $$\dot S = n - S (b_P^s P + b_A^s A + b_I^s I) + wR - m S $$ $$\dot P = S (b_P^s P + b_A^s A + b_I^s I) - g_P P - m P$$ $$\dot A = f g_P P - g_A A - m A$$ $$\dot I = (1-f) g_P P - g_I I - m I $$ $$\dot R = g_A A + (1-d) g_I I - wR - m R$$ $$\dot D = d g_I I $$

Since we do not track people dying due to non-disease causes, all the "m - arrows" are not pointing to another compartment, instead of those individuals just "leave the system". Similarly new susceptibles enter the system (are born) from "outside the system".

Modeling seasonality

Generally, when running the kind of models we investigate here, the variables change (and each gets its own equation) and the parameters remain fixed during the simulation. However, sometimes one wants parameters that change over time. In this example, we want the transmission rate to vary seasonally. The model allows you to let the transmission rate parameters vary as the simulation progresses. We model seasonality with a sin-function which has a period of 1 year. At time $t=0$, the function is zero, it reaches its max of 1 at 3 months, goes back to 0 at 6 months, reaches its minimum of -1 at 9 months and is back to zero at 12 months.

The strength of this seasonal variation, i.e. the amplitude of the sinusoidal variation, is controlled by the parameter s. Reasonable values for s are between 0 and 1. For values larger than 1, the model produces negative transmission rates, which are of course nonsensical.

Note that for even more flexibility, we could have included another parameter that can shift the sin-function along the time axis such that the minimum doesn't always have to be at t=9 but instead any other time. For simplicity, that's not done here, but would likely be good to have for modeling any real disease (e.g. influenza which has its minimum of transmission during the summer, which is around month 6-8 in the northern hemisphere).

The parameter T is set depending on the time units chosen for the model. For example if you want to run the model in units of days, the underlying simulation code will set T=365, similarly, for weeks it will be T=52. This ensures that the seasonal variation always has a period of a year.

Instead of defining time-varying parameters, one could treat each such parameter as a variable and give it their own compartment/equation. But often, it is easier and makes conceptually more sense to still consider them as parameters, and implement them in the way shown here, with explicit time-dependence.

What to do {#shinytab3}

The tasks below are described in a way that assumes everything is in units of __MONTHS__ (rate parameters, therefore, have units of inverse months). If any quantity is not given in those units, you need to convert it first (e.g. if it says a year, you need to convert it to 12 months).

Some of the simulations might take a few seconds to run. Be patient.

#this is the running counter for the records which starts at 1 
rc=1

#empty object, will hold all outcomes
alloutcomes = NULL

# Task 1
tid = 1
tasktext = "Start with 1000 susceptible, 1 initially infected and presymptomatic host, simulation duration 10 years (convert to months). Assume that only symptomatic individuals transmit, at rate _b~I~_ = 0.002 and that there is no seasonal variation (s=0). Assume that the duration of the symptomatic period is 1 month long. (Hint: The parameter _g~I~_ is the inverse of this period.) Assume that the duration of the presymptomatic period is approximately 6 days long. (Make sure you convert units correctly. First convert from days to months. You can assume that a month is roughly 30 days. Then take the inverse to get the rate.) Assume that there are no asymptomatic infections. You can, therefore, set the rate of recovery of asymptomatics, _g~A~_ to anything, it doesn't matter because nobody will be asymptomatic. Assume nobody dies due to disease, immunity does not wane, no births, and no natural deaths. If you did it correctly, you should get a single outbreak with around 203 susceptibles left at the end."

# Record for task 1
nrec = 1 # number of items to record
out_records = c("Nothing")
out_types = rep("None",nrec)
out_notes = c("")
outcomes = data.frame( TaskID = rep(tid,nrec),
                       TaskText = rep(tasktext,nrec),
                      RecordID = paste0('T',tid,'R',(1:nrec)),
                      Record = out_records, 
                      Type = out_types, 
                      Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task 


# Task 2
tid = 2
tasktext = "Next, turn on births and deaths. It's always good to check parts of a model. Here, we want to first look at births and deaths in the absence of disease. Therefore, set initial number of presymptomatic to 0, keep the number of susceptible at 1000. Set the natural mortality rate, _m_, to a value that corresponds to an average lifespan of 41 years. Recall that the unit of _m_ needs to be in 1/months. You need to convert lifespan to months, then take the inverse to get the death rate. Round to the first significant digit (i.e. 0.00X). Set the birth rate to _n=4_ per month. Also increase the simulation duration to 500 years (make sure to convert to months). (Depending on the speed of your computer, it will take several seconds for the simulation to finish). Set the initial number of pre-symptomatic individuals to 0. This means there are no infected, and thus no transmission can happen. We are focusing on only the susceptibles. Run the simulation, record the number of susceptibles at the end. Repeat this for different starting values of S, you should always end up at the same final value. It is possible to compute this final value from the equations. In the absence of disease, all infection related states (P, A, I) are zero and the equation for the susceptibles becomes _dS/dt=n-mS_. S settles down to a state where there is no more change, the so-called steady state. At that state, we have _dS/dt=0=n-mS_. You can now solve this equation for _S_ as a function of the parameters _m_ and _n_. Do that, then stick in the values for the parameters you just used to run the simulation, and convince yourself that the _S_ you get doing the math and the value you get from running the simulation are the same (up to some rounding error)."

# Record for task 2
nrec = 1 # number of items to record
out_records = c("Number susceptibles at end of simulation")
out_types = rep("Rounded_Integer",nrec)
out_notes = rep("Report the rounded integer",nrec)
outcomes = data.frame( TaskID = rep(tid,nrec),
                       TaskText = rep(tasktext,nrec),
                      RecordID = paste0('T',tid,'R',(1:nrec)),
                      Record = out_records, 
                      Type = out_types, 
                      Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task 

# Task 3
tid = 3
tasktext = "Try out different values for the birth and death rates, specifically, set _n=2_ and _m=0.002_ and then _n=6_ and _m=0.003_. Use the equation from the previous task to compute the value you expect to see _S_ go to. Then check with the simulation. For each setting, think about what you expect to see, based on your intuition and based on the theory/math, i.e. the equation from the previous task, and compare with the simulation results. Depending on the value of _S_ you start with, you might have to increase the simulation time for the model to reach the value that matches the equation."

# Record for task 3
nrec = 2 # number of items to record
out_records = c("Number susceptible at end of simulation for _n=2_ and _m=0.002_",
                    "Number susceptible at end of simulation for _n=6_ and _m=0.003_")
out_types = rep("Rounded_Integer",nrec)
out_notes = rep("Report the rounded integer",nrec)
outcomes = data.frame( TaskID = rep(tid,nrec),
                       TaskText = rep(tasktext,nrec),
                      RecordID = paste0('T',tid,'R',(1:nrec)),
                      Record = out_records, 
                      Type = out_types, 
                      Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task 



# Task 4
tid = 4
tasktext = "Set the birth and death rate to _n=4_ and _m=0.002_ and the initial number of susceptibles to the steady state value you computed above. We do that because we don't want a mixing of underlying population growth/decline dynamics on top if the ID dynamics. That would make it more complicated to understand what's going on. (You get to do that later). Ensure that if you run the simulation now, nothing changes, i.e. the susceptibles remain at the same number. Now introduce an infected individual (_P0=1_), with other values as set in task 1. Set the simulation duration to 200 years. What do you expect to see? Run the simulation, compare expectations with results. If things work right, you should see around 501 susceptibles at the end of the simulation."

nrec = 1 # number of items to record
out_records = c("Nothing")
out_types = rep("None",nrec)
out_notes = c("")
outcomes = data.frame( TaskID = rep(tid,nrec),
                       TaskText = rep(tasktext,nrec),
                      RecordID = paste0('T',tid,'R',(1:nrec)),
                      Record = out_records, 
                      Type = out_types, 
                      Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task 


# Task 5
tid = 5
tasktext = "It is possible to compute the steady state values for the model variables at the endemic scenario, i.e. when the number of infected is non-zero. This computation is a bit messy for the model above because it has the extra _P_ and _A_ states that make the math more complicated (you can try for yourself as a challenge after you've done the simple model.) To keep the math simpler, we consider a simplified approximation of the above model. The way we have set our parameters so far, individuals spend a very short time in _P_ compared to _I_. We can thus approximately ignore the _P_ compartment and assume that individuals go directly from _S_ to _I_. Further, for our parameter choices, nobody enters _A_ and _D_, so we can ignore those compartments too. This produces the following simpler SIR model:

$$ dS/dt = n -   b_I IS - m S $$ 
$$ dI/dt = b_I SI - g_I I - m I $$ 
$$ dR/dt =  g_I I  - m R $$

We can follow the procedure above for computing the steady state for _S_ to compute the endemic steady state by setting the left side of the equations to 0 and solving for the variables _S_, _I_ and _R_ as functions of model parameters only. Try it. You should get _S=(g~I~ + m)/b~I~ _ and similar equations for _I_ and _R_. For more on steady state calculations, see e.g. [@vynnycky10] or [@keeling08] (but note that each of those references uses their own notation which is not quite the same as used here.) Compare the model simulations with the equations you found. Substitute the parameter values from the task above into the equations you found for _S_, _I_ and _R_ and confirm that you get approximately the same results as you get from running the model to steady state.

To get an exact match between simulation and math would require solving the full model (i.e. the one that includes _P_ and _A_ as shown in the _Model_ tab. But those equations are getting big and ugly. If you want to give it a try, I recommend using software that can solve such equations analytically. The free [Maxima software](http://maxima.sourceforge.net/) might be a good option. If you have access to it, Mathematica or Maple are other good choices."

# Record for task 5
nrec = 1 # number of items to record
out_records = c("Nothing")
out_types = rep("None",nrec)
out_notes = c("")
outcomes = data.frame( TaskID = rep(tid,nrec),
                       TaskText = rep(tasktext,nrec),
                      RecordID = paste0('T',tid,'R',(1:nrec)),
                      Record = out_records, 
                      Type = out_types, 
                      Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task 


# Task 6
tid = 6
tasktext = "The patterns we've seen so far are oscillations/cycles (i.e. repeated outbreaks) that eventually settle down to a steady state. These oscillations are often referred to as _intrinsic_, purely produced by the interplay between depletion of susceptibles due to infection and replenishment due to birth. Waning immunity is another mechanism by which susceptibles can be replenished and one can thus get repeat outbkreaks. To investigate that, set everything back to the values from task #1. Now, assume that immunity wanes within on average 10 months and set simulation time for 20 years. Record the number of susceptible and infected at the end of the simulation. Double the _duration_ of immunity. What do you expect? What do you see? Run again and record the number susceptible and infected at steady state. Try a few more waning immunity values and observe how the osciallations and the steady state is affected (you might have to change the simulation time to ensure you reach steady state)."

# Record for task 6
nrec = 4 # number of items to record
out_records = c("Number of susceptible at end of simulation, 10 months immunity",
                    "Number of infected at end of simulation, 10 months immunity",
                    "Number of susceptible at end of simulation, 20 months immunity",
                    "Number of infected at end of simulation, 20 months immunity")
out_types = rep("Rounded_Integer",nrec)
out_notes = rep("Report the rounded integer",nrec)
outcomes = data.frame( TaskID = rep(tid,nrec),
                       TaskText = rep(tasktext,nrec),
                      RecordID = paste0('T',tid,'R',(1:nrec)),
                      Record = out_records, 
                      Type = out_types, 
                      Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task 

# Task 7
tid = 7
tasktext = "We can again compare the steady state values we get from the simulation with those we get from a model. To do so, we follow the approach above. The model we now solve is one that doesn't include births and deaths but does including waning immunity. The equations are:

$$ dS/dt = w R  -   b_I IS $$ 
$$ dI/dt = b_I SI - g_I I $$ 
$$ dR/dt =  g_I I  - w R $$

You can set the left sides to zero and solve for the steady states of _S_, _I_ and _R_. By solving the second equation for _S_, you will find that its steady state is independent of _w_, which you should have also seen in your simulations. _I_ and _R_ are not uniquely determined, you'll realize you can't solve them as before. Going deeper into what this means is beyond what we want to discuss here. For more on steady states, see e.g. [@keeling08] or some mathematical biology textbook."

# Record for task 7
nrec = 1 # number of items to record
out_records = c("Nothing")
out_types = rep("None",nrec)
out_notes = c("")
outcomes = data.frame( TaskID = rep(tid,nrec),
                       TaskText = rep(tasktext,nrec),
                      RecordID = paste0('T',tid,'R',(1:nrec)),
                      Record = out_records, 
                      Type = out_types, 
                      Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task 


# Task 8
tid = 8 
tasktext = "Now combine births and deaths with waning immunity. Explore how the two different mechanisms interact and influence the observed ID dynamics by choosing different values for waning immunity, births and deaths and exploring how those different values impact both the cycles and the steady state you see."

# Record for task 8
nrec = 1 # number of items to record
out_records = c("Nothing")
out_types = rep("None",nrec)
out_notes = c("")
outcomes = data.frame( TaskID = rep(tid,nrec),
                       TaskText = rep(tasktext,nrec),
                      RecordID = paste0('T',tid,'R',(1:nrec)),
                      Record = out_records, 
                      Type = out_types, 
                      Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task 


# Task 9
tid = 9
tasktext = "So far, the cycles were due to the intrinsic dynamics of the system. Now we'll explore what happens if some outside phenomenon - e.g. the weather - influences some model parameters. In our model, the transmission parameters can vary sinusoidally with a period of 1 year. Set everything as in task #1, do a quick run to make sure you only get one outbreak. Now run the model with increasing seasonality by setting s = 0.001, 0.01, 0.1, and 1. Even though you are letting the transmission vary seasonally, you will only see a single outbreak, though the dynamics changes slightly. Why is that so? Why don't you see patterns like before? Think about the concept of _resource replenishment_ and what that means here."

# Record for task 9
nrec = 1 # number of items to record
out_records = c("Nothing")
out_types = rep("None",nrec)
out_notes = c("")
outcomes = data.frame( TaskID = rep(tid,nrec),
                       TaskText = rep(tasktext,nrec),
                      RecordID = paste0('T',tid,'R',(1:nrec)),
                      Record = out_records, 
                      Type = out_types, 
                      Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task 



# Task 10
tid = 10
tasktext = "I hope you figured out above that without replenishment of resources (i.e. susceptibles), you can't get multiple outbreaks, even with seasonally changing transmission. let's introduce replenishment by assuming immunity wanes within 10 months. Start with no seasonal variation, i.e. _s=0_. You should see one decent outbreak with a peak at around 10 months, another small one with a peak at around 30 months, and not much else (you might want to use plotly so you can easier zoom in/out). Now set seasonal variation strength to _s=0.5_. You should see that the outbreaks keep going now. Approximately how many months are there between outbreak peaks? (Remember in what way we let the transmission strength vary and make sure what we specify in the model and the time between peaks makes sense)."  

# Record for task 10
nrec = 1 # number of items to record
out_records = c("Months between peaks (approximate/rounded)")
out_types = rep("Rounded_Integer",nrec)
out_notes = rep("Report the rounded integer",nrec)
outcomes = data.frame( TaskID = rep(tid,nrec),
                       TaskText = rep(tasktext,nrec),
                      RecordID = paste0('T',tid,'R',(1:nrec)),
                      Record = out_records, 
                      Type = out_types, 
                      Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task 


# Task 11
tid = 11 
tasktext = "Explore how changing the rate of waning immunity interacts with seasonality to affect patterns. To that end, run simulations with _s=0.5_ and _w=0.2_ and _w=0.05_. You can further explore by switching off waning immunity and turning on births and deaths instead. You can also have both mechanisms present and explore how they interact with seasonality. And of course, you can set parameters to allow for asymptomatic infections and deaths, ore have an underlying growing or declining population in the absence of disease and let the disease run on top of it. There is lots to explore. While it can be fun running simulations with random parameters, I suggest you learn more if you are deliberate in your explorations. After you set the model inputs to some values and before you run a simulation, think about what you expect. Then run and compare results with expectations. If they don't agree, try to figure out what's going on. This way, you will gain more and more insight into how the different components and processes influence the observed outcome - the basics of doing science."

# Record for task 11
nrec = 1 # number of items to record
out_records = c("Nothing")
out_types = rep("None",nrec)
out_notes = c("")
outcomes = data.frame( TaskID = rep(tid,nrec),
                       TaskText = rep(tasktext,nrec),
                      RecordID = paste0('T',tid,'R',(1:nrec)),
                      Record = out_records, 
                      Type = out_types, 
                      Note = out_notes)
alloutcomes = rbind(alloutcomes,outcomes)
rc = rc + nrec #increment record counter by number of outcomes to record for this task 

#save the fully filled task table to a tsv file
alloutcomes$QuizID = paste0(packagename,"_",appsettings$appid)
alloutcomes$AppTitle = appsettings$apptitle
alloutcomes$AppID = appsettings$appid
#remove a few variables from the data frame
savedoutcomes <- dplyr::select(alloutcomes,QuizID,AppID,AppTitle,TaskID,TaskText,RecordID,Record,Type,Note)     
write.table(savedoutcomes, paste0(appsettings$appid,"_tasktable.tsv"), append = FALSE, sep = "\t", row.names = F, col.names = TRUE)

# Take all the text stored in the table and print the tasks and items to record
write_tasktext(alloutcomes)

Further Information {#shinytab4}

• This app (and all others) are structured such that the Shiny part (the graphical interface you see and the server-side function that goes with it) calls an underlying R script (or several) which runs the simulation for the model of interest and returns the results.
• For this app, the underlying function(s) running the simulation are called r appsettings$simfunction. You can call them directly, without going through the shiny app. Use the help() command for more information on how to use the functions directly. If you go that route, you need to use the results returned from this function and produce useful output (such as a plot) yourself.
• You can also download all simulator functions and modify them for your own purposes. Of course to modify these functions, you'll need to do some coding.
• For examples on using the simulators directly and how to modify them, read the package vignette by typing vignette('DSAIDE') into the R console.
• For some papers discussing ID patterns and seasonality, see e.g. [@dowell2001; @altizer2006].

References

DSAIDE documentation built on Aug. 24, 2023, 1:07 a.m.