DSAIRM Solutions: Basic Virus Model
2022-05-17
This is a walk-through for the Basic Virus Model app in DSAIRM.
If you haven’t already loaded the DSAIRM package, it needs to be loaded.
library(DSAIRM)Task 1
Let’s start by considering an acute viral infection. Set the initial conditions to 105 uninfected cells, no infected cells, and 10 virus particles. We make the assumption that on the timescale of an acute infection (several days), the processes of natural, uninfected cell turnover are so slow that they can be ignored. Set values for the uninfected cell birth and death rates to reflect this assumption. Assume also that infected cells have an average life-span of 1 day, and virus has a life-span of 6 hours (remember that the inverse of the lifespan is the rate of death, and make sure you convert to the right units). Set that the virus production by an infected cell is 100 virions per day and that the rate at which new cells become infected is 10-6. Assume there is no need to do any unit adjustment/conversion (i.e. the value of that parameter is 1).
Run the simulation for 50 days; produce plots with and without log-scales. You should get a single, acute infection with virus and infected cells rising at first and then declining. At the end you should be left with around 11069 uninfected cells and no infected cells or virus.
Record
- Number of infected cells at the peak of infection
You likely performed the tasks by changing settings in the graphical interface. As you know from the package tutorial, it is also possible to call the underlying simulation functions directly. This is the approach we take for this and all other solutions. The Further Resources section of the app tells you that for this app, the relevant model function is called simulate_basicvirus_ode().
We can run the underlying simulation model to get results for specific parameter choices by calling the simulation function with the specified settings, which can be done with this line of code.
#call the simulator function with the specified settings
sim_result <- simulate_basicvirus_ode(U = 1e+05,
I = 0,
V = 10,
n = 0,
dU = 0,
dI = 1,
dV = 4,
b = 1e-06,
p = 100,
g = 1,
tstart = 0,
tfinal = 50,
dt = 0.1
)The result that is returned is saved in the object sim_result. This is a list object, and it contains an object called ts which is the time series data from the simulation. That is, ts contains the values for each of the variables (U, I, and V, here) and the time. All simulator functions return a list, the elements contained in that list can differ between simulators. The help file for each simulator function describes what is returned.
We can plot the result obtained from the simulator using basic R plotting commands. The following lines of code produce a figure that looks similar to the one you get through the graphical interface.
#plot using the base R plotting functions; all further plotting will be done using the functions included in the package
plot(sim_result$ts[,"time"], sim_result$ts[,"U"], type = 'l', col='blue',
xlab='Time', ylab='Numbers', ylim=c(min(sim_result$ts)-0.1*min(sim_result$ts),max(sim_result$ts)+0.1*max(sim_result$ts)))
lines(sim_result$ts[,"time"], sim_result$ts[,"I"], type = 'l', col='red')
lines(sim_result$ts[,"time"], sim_result$ts[,"V"], type = 'l', col='green')
legend("topright", legend=c('U','I', 'V'), lwd=2, col=c('blue','red', 'green'))
To match the figure produced by the GUI, we will instead use the built-in plotting function generate_ggplot that comes with DSIARM. A similar function, generate_text, can be used to produce the information printed below each plot in the user interface. These functions require input to be in a specific form: a vector of lists. These lists are the time series ts from our simulation results (or sometimes lists with other structures). Here, there is only one figure, therefore the vector is of length 1. To create multiple plots, we would add another list to the vector for each plot we wanted to make.
For this outcome, we can simply apply the list command to the result returned from the simulator, and then send it to the generate_ functions. The following lines of code then generate the figure and text you see in the GUI.
generate_ggplot(list(sim_result))
generate_text(list(sim_result))Minimum / Maximum / Final value of U: 11068.65 / 1e+05 / 11068.65
Minimum / Maximum / Final value of I: 0.00 / 18305.62 / 1.5e-06
Minimum / Maximum / Final value of V: 4.6e-05 / 448313.13 / 4.6e-05
For the remainder of the solutions, we will use the built-in plot and text generation functions.
We can reproduce the above plot this time showing the y-axis on the log 10 scale using the following code.
generate_ggplot(list(list(ts = sim_result$ts, yscale = "log10")))
generate_text(list(sim_result))Minimum / Maximum / Final value of U: 11068.65 / 1e+05 / 11068.65
Minimum / Maximum / Final value of I: 0.00 / 18305.62 / 1.5e-06
Minimum / Maximum / Final value of V: 4.6e-05 / 448313.13 / 4.6e-05
Note that the simulation result time series is nested inside two list() functions. The argument for the generate_ggplot() function takes a list of each figure to create (i.e., the first, outermost list()) and each figure element is also a list object (i.e., the second, innermost list()) containing the time series dataframe and, optionally, appropriate metadata (e.g., a yscale argument). You could also append the metadata to the original sim_result object, but that would override the defaults necessitating another specification to return the output plot to the original scale.
The text output provides the answer to the task question. We can also pull it out directly from the sim_result object, without using the generate_text function. We want the number of infected cells at the peak of infection which is the maximum number in the I column reported in the ts dataframe. To get this value we use the the max() function. We then round to the nearest integer and print the result. The code below shows how the value is obtained.
Imax=max(sim_result$ts$I)
Imax_rounded=round(Imax,0)
print(Imax_rounded)## [1] 18306Answer for Task 1
Number of infected cells at the peak of infection: 18306
Task 2
Slowly increase the virus death rate in increments of 0.5. Contemplate what you expect to see, then run the simulation to compare. Keep increasing until you get essentially no more infection (i.e., < 1% initial uninfected cells become infected). You will have to adjust the simulation time for that, too.
Record
Virus death rate at which no infection occurs
Virus lifespan (in days) corresponding to the death rate at which no infection occurs
For this task, we will run several simulations each varying the virus death rate until we get essentially no infection. To do so, we will use a for loop to iterate the simulation function over a sequence of values for the virus death rate parameter and output each simulation result to a new list object using the following code:
virus_death_rates <- seq(4, 12, by = 0.5)
sim_results <- list()
for(i in 1:length(virus_death_rates)){
sim_results[[i]] <- simulate_basicvirus_ode(U = 1e+05,
I = 0,
V = 10,
n = 0,
dU = 0,
dI = 1,
dV = virus_death_rates[i],
b = 1e-06,
p = 100,
g = 1,
tstart = 0,
tfinal = 500, # increased sim time
dt = 0.1
)
} Now with the simulations run for each of the virus death rate values, we can assess the extent of the infection for each parameter value. One way to accomplish this is to plot the number of uninfected cells at the end of the simulations against the virus death rate. First, we will extract the U value in the last row of the ts dataframe using the sapply() and tail() functions. Then, we will create a simple plot() to show the results.
Ufinals <- sapply(sim_results, function(x){tail(x$ts$U,1)})
plot(x = virus_death_rates, y = Ufinals, pch = 16)
From the above plot, we can note an inflection point when the virus death rate parameter is 10. To identify this inflection point programmatically, we need to think a bit differently. One way is to set a threshold for where we define essentially no infection. Suppose we say that if less than 1% of cells are infected, then there was essentially no infection. Then, we can determine the minimum virus death rate corresponding to when that condition is met, as in the following code:
Uinitial <- 100000
threshold <- Uinitial * 0.01
# min() function wrapped around subset virus_death_rates sequence
# subset virus_death_rates using [] and which() functions
# logical operator < to evaluate the condition
inflection <- min(virus_death_rates[which(Uinitial - Ufinals < threshold)])
print(inflection)## [1] 10Now that we know the virus death rate of the inflection point, we can convert this to the average virus lifespan via the following code:
# inverse of death rate is lifespan
lifespan <- 1 / inflectionAnswer for Task 2
Virus death rate at which no infection occurs: 10
Virus lifespan (in days) corresponding to the death rate at which no infection occurs: 0.1
Task 3
Set the virus death rate back to what it was in task 1. Now change the virus production rate (in increments of 10) until you reach the value at which the virus does not cause any infection. You can also repeat this process for the infected cell death rate and the infection rate.
Record
- Virus production rate at which the virus does not cause any infection
For this task, we will again run several simulations each varying the virus production rate until we get essentially no infection via a for loop. Then, we extract the values for the uninfected cells at the end of the simulation to assess the inflection point (both visually and programmatically).
virus_production_rates <- seq(20, 100, by = 10)
sim_results <- list()
for(i in 1:length(virus_production_rates)){
sim_results[[i]] <- simulate_basicvirus_ode(U = 1e+05,
I = 0,
V = 10,
n = 0,
dU = 0,
dI = 1,
dV = 4,
b = 1e-06,
p = virus_production_rates[i],
g = 1,
tstart = 0,
tfinal = 250, # increased sim time
dt = 0.1
)
}
Ufinals <- sapply(sim_results, function(x){tail(x$ts$U,1)})
plot(x = virus_production_rates, y = Ufinals, pch = 16)
Uinitial <- 100000
threshold <- Uinitial * 0.01
# max() function wrapped around subset virus_death_rates sequence
# subset virus_death_rates using [] and which() functions
# logical operator < to evaluate the condition
inflection <- max(virus_production_rates[which(Uinitial - Ufinals < threshold)])
print(inflection)## [1] 40From the above plot and the calculations, we note the inflection occurs when the virus production rate parameter is 40.
Answer for Task 3
Virus production rate at which the virus does not cause any infection: 40
Task 4
A well-studied quantity in infectious disease epidemiology is the basic reproductive number (R0), which determines if a pathogen can cause an outbreak at the population level. In the absence of any interventions or other changes, the basic reproductive number also determines the size of the outbreak. An equivalent R0 can be defined for a within-host model to determine if you get an infection or not. For this virus model (with no births and deaths of uninfected cells, i.e. n=dU=0), R0 = bpU0/(dV dI).
Plug numbers for the parameters from your simulations in task 2 and 3 into the equation for R0 to figure out what value R0 needs to be for there (not) to be an infection. Figure out the threshold value for R0 at which you go from no infection to having an infection.
To learn more about R0, see e.g. (Heffernan, Smith, and Wahl 2005; Roberts 2007; Beauchemin et al. 2008). Some of those references describe R0 in the context of infectious disease epidemiology, but if you replace humans/hosts with cells, the same concepts apply at the within-host level.
Record
- Threshold value for R0 at which you go from no infection to having an infection
In its most basic form, R is a calculator. So, we can plug in the values from the previous tasks to calculate the value for R0 given.
# parameter values
U = 1e+05
dI = 1
dV = 4
b = 1e-06
p = 100
# original
R0 <- b*p*U/(dV*dI)
print(R0)## [1] 2.5# using other virus death rate
R0 <- b*p*U/(10*dI)
print(R0)## [1] 1# using other virus production rate
R0 <- b*40*U/(dV*dI)
print(R0)## [1] 1The basic reproduction number refers to the average number of secondary infections generated from a single infection in a completely susceptible population. When an infectious disease (or agent) has an R0 > 1, the infection can spread. When an infectious disease (or agent) has an R0 < 1, the infection wanes. So, R0 = 1 is the threshold or tipping point at which you go from no infection to having an infection.
Answer for Task 4
Threshold value for R0 at which you go from no infection to having an infection: 1
Task 5
Without birth/production of new uninfected cells, the most you can get is a single acute infection (or no infection at all). To convince yourself that it is impossible to produce a chronic infection, play around with the model, try all kinds of parameter values (but keep n=0). Production of new uninfected cells is an example of resource replenishment. This is needed to allow a steady state/chronic infection, and this concept applies in general (e.g. on the population level, new susceptible individuals need to be created either through birth or through losing immunity). Let’s explore the model with uninfected cell production, i.e. resource replenishment.
We start by focusing on the dynamics of uninfected cells only. To that end, set the number of initial infected cells and virus particles to 0. Keep the number of uninfected cells at 105, set birth and death of uninfected cells to zero. Run the simulation. Nothing should happen, uninfected cells should stay at their starting value. Now, play around with birth rate and death rate of uninfected cells and see how that affects the dynamics. The number of uninfected cells once the system has settled down only depends on the birth and death rate, not the starting conditions. Confirm this by trying different values for U0 while keeping birth and death rate at some fixed values. One can write down an equation for uninfected cells at steady state as a function of birth and death rate, i.e. \(U_s = f(n,d_U)\), where f() is the mathematical symbol for some function. In this case, it is a very simple function. Based on your explorations of different values for birth and death rate and the resulting values of Us, figure out this equation.
To test your solution, set birth rate to 20000, the initial uninfected cells to 105, and initial virus to 0, and discover whether the value for the death rate keeps the number of uninfected cells unchanged at 105.
Record
- Value for uninfected cell death rate at steady state
One way we can go about solving this is to examine the given model equation for the change in the uninfected cell population:
\[\dot{U} = n - d_UU - bUV\]
Since we are eliminating the infection dynamics, i.e., \(V=0\), and we are interested in a steady state population, i.e., \(\dot{U}=0\), we can drop those terms and rearrange the equation to:
\[U_s = \frac{n}{d_U}\] which states that the steady state population size, \(U_s\), is equal to the ratio of birth rate to death rate.
Now, if we are given values for two of the three terms in the equation, as in the task, we can again rearrange to solve for the third:
\[d_U = \frac{n}{U_s}.\]
dU <- 20000 / 1e5Answer for Task 5
Value for uninfected cell death rate at steady state: 0.2
Task 6
Now we’ll explore an infection in the presence of uninfected cell birth and death. Set all parameters as in task 1. Set birth and death as described at the end of the previous task. Run the simulation. You should get an initial large increase in virus load, which then settles down and reaches a steady state of around 295000. Similarly, the variables U and I settle down to steady state values.
Record
Number of uninfected cells at steady state
Number of infected cells at steady state
Now with the steady state values input, we can rerun our simulation with the infection dynamics.
sim_result <- simulate_basicvirus_ode(U = 1e+05,
I = 0,
V = 10,
n = 20000,
dU = 0.2,
dI = 1,
dV = 4,
b = 1e-06,
p = 100,
g = 1,
tstart = 0,
tfinal = 50,
dt = 0.1
)
generate_ggplot(list(sim_result))
generate_text(list(sim_result))Minimum / Maximum / Final value of U: 30223.57 / 1e+05 / 40401.72
Minimum / Maximum / Final value of I: 0.00 / 24148.06 / 11916.83
Minimum / Maximum / Final value of V: 4.76 / 593132.68 / 294946.13
And we can extract the values for the numbers of uninfected and infected cells at steady state.
Ufinal <- tail(sim_result$ts$U,1)
Ufinal_rounded <- round(Ufinal, 0)
print(Ufinal_rounded)## [1] 40402Ifinal <- tail(sim_result$ts$I,1)
Ifinal_rounded <- round(Ifinal, 0)
print(Ifinal_rounded)## [1] 11917Answer for Task 6
Number of uninfected cells at steady state: 40402
Number of infected cells at steady state: 11917
Task 7
Investigate how the steady state values for U, I and V depend on the parameters b, p, dV and dI. You might need to increase the simulation time to ensure the system has settled down to its steady state. Once the system has settled down, there are no more changes in the numbers for each compartment. Mathematically, that means that the left side of the differential equations becomes 0, and they turn into the following algebraic equations: 0 = n - dU U - bUV, 0 = bUV - dI, 0 = pI - dV V - gb UV. One can solve those equations for each of the compartments to get a mathematical expression of what U, I and V are at steady state. If your algebra is not too rusty, try to do this. You should find that \(U_s = d_I d_V/(b(p - d_I g))\), \(V_s = (bnp-bd_Ign-d_Id_Ud_V)/bd_Id_V\) and and equation with the same numerator but slightly different denominator for Is (the subscript s denotes that these are the steady-state values of the variables). If you haven’t done algebra in a while, or if you find doing math by hand too tedious, modern computer software often helps. R cannot solve such equations analytically, but other software packages can. The main ones used for analytic math are Mathematica and Maple. Both are powerful and expensive. If you only need to solve simple equations occasionally, there is Maxima, which is free. You can download it and enter the equations above and it will solve it for you. Note that once you go beyond 4-5 variables, the steady state equations are usually very complicated, often so much so that they are not useful anymore. And once you go beyond 5 variables, in most cases your software will struggle to give you something meaningful. Fortunately, while it is less quick and elegant, you can always simulate your model and see what (if any) steady state it reaches.
Once you found the steady state equations, either by hand or with the help of some computer software, check that your equations agree with the simulations. Plug the values for the parameters into each of the equations and see if the steady state values Us, Is and Vs you computed with the equations is the same as you get as steady state value from the simulation. If that’s not the case, it means your equations aren’t right yet. It is useful to note that while the total numbers for each variable do not change at steady state, this is a dynamic equilibrium. There are still constantly cells and virus being produced and destroyed, it just so happens that the production and destruction mechanisms are equally strong and thus the overall numbers do not change.
Record
- TRUE or FALSE: There is only a single, unique combination of the parameters b, p, dV, and dI that can produce a steady state.
We’ve shown in the previous task one set of values which can generate a steady state in our model. To show simply that there exist multiple combinations of parameters that can each produce a steady state, we can sub in zero values for the parameters. Technically, these are steady state solutions since nothing is changing throughout the simulation. Furthermore, there exist many other non-zero combinations of parameter values which can produce a steady state infection. With increasing model complexity, this set of values expands as well which can lead to issues of model identifiability.
Answer for Task 7
TRUE or FALSE: There is only a single, unique combination of the parameters b, p, dV, and dI that can produce a steady state.: FALSE
Task 8
Continue to explore the model. Even though it’s a fairly simple model, you can get interesting dynamics from it, such as acute infections and chronic infections. Contemplate what specific pathogens this model could represent. Also note that this model does not contain an immune response. The interactions between cells and virus are enough to produce patterns of infection dynamics that broadly agree with patterns we can see for real infections. This of course does not mean the immune response is not important. But it does illustrate that if all we have is (noisy) virus kinetics data, we are likely able to capture that dynamics with many different types of models, including a simple one like this that is likely not too realistic for any given pathogen.
Record
- Nothing
This is an open-ended exploration, with nothing to report.
Answer for Task 8
Nothing