Multi Pathogen Dynamics - Documentation

#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 with 2 types of pathogens which can serially or simultaneously infect the hosts. Only one type of host is included. Read about the model in the "Model" tab. Then do the tasks described in the "What to do" tab.

This app assumes that you have worked through several of the previous ones, at least those in the Basics and Reproductive Number sections.

Learning Objectives

The Model {#shinytab2}

Model Overview

This model tracks susceptible hosts, hosts infected with either pathogen 1, pathogen 2 or both, and individuals recovered from infection with pathogen 1, pathogen 2 or both. As usual, infected hosts are assumed to be infectious.

The following compartments are included:

The included processes/mechanisms are the following:

The biological meaning of the assumptions made here is that there is no immunological interaction between pathogens. One could alter the model to allow for this, e.g. to include short-term immunological cross-protection or similar features.

Model Implementation

The flow diagram and the set of differential equations for the mathematical model implemented in this app are as follows:

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

$$\dot S = - (b_{1} (I_1+I_{1X}) + b_{2} (I_2+I_{2X}) + b_{12}I_{12}) S $$ $$\dot I_1 = (b_{1} (I_1+I_{1X}) + ab_{12} I_{12})S - (g_1 + b_{2} (I_2+I_{2X}) + b_{12} I_{12}) I_1$$ $$\dot I_2 = (b_{2} (I_2+I_{2X}) + (1-a) b_{12} I_{12})S - (g_2 + b_{1}(I_1 + I_{1X}) + b_{12} I_{12}) I_2$$ $$\dot I_{12} = (b_{2} (I_2+I_{2X}) + b_{12} I_{12}) I_1 + (b_{1}(I_1 + I_{1X}) + b_{12} I_{12}) I_2 - g_{12} I_{12}$$ $$\dot R_1 = g_1 I_1 - (b_2 (I_2 + I_{2X}) + b_{12} I_{12}) R_1$$ $$\dot R_2 = g_2 I_2 - (b_1 (I_1 + I_{1X}) + b_{12} I_{12}) R_2$$ $$\dot I_{1X} = (b_1 (I_1 + I_{1X}) + b_{12} I_{12}) R_2 - g_{1} I_{1X}$$ $$\dot I_{2X} = (b_2 (I_2 + I_{2X}) + b_{12} I_{12}) R_1 - g_{2} I_{2X}$$ $$\dot R_{12} = g_{1} I_{1X} + g_{2} I_{2X} + g_{12} I_{12} $$

What to do {#shinytab3}

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

#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 = "Set the initial number of infected with pathogen 1 to _I~1~_ = 1, all other infected at zero. Susceptibles at 1000.
Set all transmission rates to 0.001, all recovery rates at 0.5. Set fraction of pathogen 1 infections to _a=0.5_.
Run the simulation for 100 days. Check that you only get a single outbreak with pathogen 1, with a peak of around 154 infected. Now set initial number of infected with pathogen 1 to 0, with pathogen 2 to 1. Run the simulation and check that you get exactly the same size outbreak with pathogen 2. The number susceptible left at the end of each outbreak is the same. This finding should not surprise you: In the presence of only a single pathogen, you can only get an outbreak caused by that pathogen. For both pathogens to infect hosts, both need to be present."

nrec = 1 # number of items to record
out_records = c("Number susceptible left 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 2
#########################
tid = tid + 1
tasktext = "Leave everything as before but now start with 1 infected for both pathogen 1 and 2. You should see outbreaks for both pathogens with the same dynamics (in fact, you won't see the curves for _I~1~_, _R~1~_, etc. since they are covered by the curves for _I~2~_, _R~2~_.) If you compare the peak number of infected for either _I~1~_ or _I~2~_ to the value you found in the previous task, you will notice it is lower (around 80). The two pathogens are competing for susceptible hosts, thus neither can cause quite as big an outbreak when both pathogens are present than when only one of them is present. However, the 2 outbreaks together lead to more total infections, which you can confirm by looking at the final number of susceptible and comparing to the previous task."

nrec = 1 # number of items to record
out_records = c("Number susceptible left 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 = tid + 1
tasktext = "It basically never the case that two pathogens are the same in their transmission rates, recovery rates, starting values, etc. Let's mimic a more realistic scenario by adjusting the settings for pathogen 1. Double _b~1~_ to 0.002, keep the remaining parameters as before. Now pathogen 1 is twice as tranmissible as pathogen 2, thus we should expect it to dominate. Run the simulation and look at the peak infected for pathogen 1 and 2 to confirm that this is the case. As you have encountered many times previously, things are not linear, i.e. the fact that pathogen 1 is twice as transmissible than pathogen 2 does not lead to a peak for pathogen 1 that is twice as large as that for pathogen 1, in fact the difference is much larger, almost 100-fold and pathogen 1 completely dominates."
nrec = 2 # number of items to record
out_records = c("Maximum number infected with pathogen 1", "Maximum number infected with pathogen 2")
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 = tid + 1
tasktext = "Now also double the rate of recovery for pathogen 1 to _g~1~_ = 1. Recall that in this simple SIR model, for just a single pathogen _i_ (if the other one isn't present) the reproductive number is _S~0~b~i~_/_g~i~_. This means that pathogen 1 and 2 have now the same reproductive number (pathogen 1 has both a doubled transmission and recovery rate.) Given that they have the same reproductive number, do you expect to see outbreaks for the two pathogens of the same size? Once you contemplated what to expect, check by running the simulation."
nrec = 2 # number of items to record
out_records = c("Maximum number infected with pathogen 1", "Maximum number infected with pathogen 2")
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 5
#########################
tid = tid + 1
tasktext = "You should have found in the previous task that despite the same reproductive number values, pathogen 1 causes a larger outbreak. If this is confusing, recall that _R~0~_ says nothing about the speed of spread of a pathogen, only its infectiousness (how many new infections are produced by one infected host). If there is competition for hosts, a pathogen that spreads faster (pathogen 1 here) is advantaged. If none of that sounds familiar, you might want to revisit the reproductive number apps. We can give pathogen 2 a bit of a boost by letting it start with more initial infected. Let's see what that does. Set _I~2~_ = 10, keep everything else the same. Run the simulation. You'll find that the peaks for the infected for the 2 pathogens are getting closer. It's like in a race where runner 1 is faster, but runner 2 gets to start a distance ahead of runner 1. The initial difference in starting distance, the speed of the 2 runners, and the length of the race all determine the overall result. Feel free to explore a bit more how different initial values and transmission/recovery rates affect the competition between the 2 pathogens."
nrec = 2 # number of items to record
out_records = c("Maximum number infected with pathogen 1", "Maximum number infected with pathogen 2")
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 6
#########################
tid = tid + 1
tasktext = "Next, we'll allow for a difference between the 2 pathogens by assuming that if a host is double infected, they are more likely to transmit one of the pathogens to an uninfected host. Set all values as in task 1 (I~1~ = I~2~ = 1, all transmission rates 0.001, all recovery rates 0.5). Then set _a_ = 0.9. This means that a double-infected host will mainly transmit pathogen 1 to other hosts. Thus we should expect pathogen 1 to have an advantage. Run the simulation to confirm this."
nrec = 2 # number of items to record
out_records = c("Maximum number infected with pathogen 1", "Maximum number infected with pathogen 2")
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 = tid + 1
tasktext = "Change the rate at which double infected hosts recover from 0.5 to _g~12~_ = 5. Keep everything else unchanged. Run the simulation and look at the maximum for each pathogen. You will find that the values are much closer compared to the previous task. Can you figure out why? Think about the advantage that pathogen 1 has (all infections caused by double-infected hosts are with pathogen 1) and how that advantage is influenced as you increase the recovery rate _g~12~_."
nrec = 2 # number of items to record
out_records = c("Maximum number infected with pathogen 1", "Maximum number infected with pathogen 2")
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 8
#########################
tid = tid + 1
tasktext = "Keep exploring. For instance explore how the infection rate of double-infected hosts, _b~12~_ affects the competition between the pathogens, and how that depends on other parameter settings. If you want, think of 2 real-world pathogens that might compete (e.g. influenza and RSV) and see if you can find parameter values that approximate those pathogens and their possible interactions. Multi-pathogen models can easily get complicated, and here we left out a good bit of possible interactions. If you want to dig deeper, you can get yourself the code for this app (see _Further Resources_) and modify it to allow for further processes and interactions between the pathogens."
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}

References



Try the DSAIDE package in your browser

Any scripts or data that you put into this service are public.

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