DSAIRM Solutions: Chronic Virus and IR Model
2022-05-18
This is a walk-through for the Chronic Virus and IR Model app in DSAIRM.
If you haven’t already loaded the DSAIRM package, it needs to be loaded.
library(DSAIRM)
library(ggplot2) #we'll use it for plottingTask 1
Let’s start by considering a chronic viral infection in the absence of the immune response. Set the initial conditions to 105 uninfected cells, no infected cells, and 10 virus particles. 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). Set the birth rate of uninfected cells to 104, and the life span od the uninfected cells to 10 days. Set the starting values for the immune response variables F and T to 0, and also set the two immune response induction rates to 0. With those choices, the values for the other immune response related parameters do not matter. If you want, you can also set them all to 0. Run the simulation for 200 days. You should get one major rise in virus, followed by a few minor peaks and then settling down to a steady state where the virus levels are at around 147500.
Record
Number of uninfected cells at the end of the simulation
Number of infected cells at the end of the simulation
As you have seen before, we can reproduce what you get by running the models through the GUI with a few simple lines of code. The Further Resources section of the app tells you the name of the simulation function for each app that you need to call. For this one, it is simulate_chronicvirusir_ode().
Here is the function call for the first task.
#call the simulator function with the specified settings
res <- simulate_chronicvirusir_ode(U = 1e+05, I = 0, V = 10, F = 0, T = 0,
n = 1e4, dU = 0.1,
dI = 1, dV = 4, b = 1e-06, p = 100, g = 1,
rF = 0, dF = 0, kF = 0,
rT = 0, dT = 0, kT = 0,
tstart = 0,
tfinal = 200,
dt = 0.1
)The following lines of code then generate the figure and text you see in the GUI.
generate_ggplot(list(res))
generate_text(list(res))Minimum / Maximum / Final value of U: 24698.06 / 1e+05 / 40404.04
Minimum / Maximum / Final value of I: 0.00 / 21121.91 / 5959.60
Minimum / Maximum / Final value of V: 4.76 / 517886.89 / 147500.00
Minimum / Maximum / Final value of F: 0.00 / 0.00 / 0.00
Minimum / Maximum / Final value of T: 0.00 / 0.00 / 0.00
This pulls out the numbers you are asked to report.
Ufinal <- round(tail(res$ts$U,1),0)
Ifinal <- round(tail(res$ts$I,1),0)Answer for Task 1
Number of uninfected cells at the end of the simulation: 40404
Number of infected cells at the end of the simulation: 5960
Task 2
As you might expect, and as discussed in the Basic Virus app, changing the parameters associated with the virus and cell kinetics can change the results. How a certain parameter impacts the steady/chronic level is not always intuitively obvious. Give that a quick try. Double the rate at which uninfected cells are produced, i.e. set n = 2x104. One might expect that this should increase the number of uninfected cells at steady state. Run the simulation to determine that this is NOT the case, that instead the number of infected cells changes. Revisit the Basic Virus app if it is unclear to you why this happens.
Record
Number of uninfected cells at the end of the simulation
Number of infected cells at the end of the simulation
Now with twice the rate of uninfected cell production.
#call the simulator function with the specified settings
res <- simulate_chronicvirusir_ode(U = 1e+05, I = 0, V = 10, F = 0, T = 0,
n = 2e4, dU = 0.1,
dI = 1, dV = 4, b = 1e-06, p = 100, g = 1,
rF = 0, dF = 0, kF = 0,
rT = 0, dT = 0, kT = 0,
tstart = 0,
tfinal = 200,
dt = 0.1
)generate_ggplot(list(res))
generate_text(list(res))Minimum / Maximum / Final value of U: 19220.59 / 147790.44 / 40404.04
Minimum / Maximum / Final value of I: 0.00 / 57096.36 / 15959.60
Minimum / Maximum / Final value of V: 4.80 / 1380273.92 / 4e+05
Minimum / Maximum / Final value of F: 0.00 / 0.00 / 0.00
Minimum / Maximum / Final value of T: 0.00 / 0.00 / 0.00
The number of uninfected cells does not change, but the infected cells do change.
This again pulls out the numbers you are asked to report.
Ufinal <- round(tail(res$ts$U,1),0)
Ifinal <- round(tail(res$ts$I,1),0)Answer for Task 2
Number of uninfected cells at the end of the simulation: 40404
Number of infected cells at the end of the simulation: 15960
Task 3
Now turn on the immune response. Set everything as you did for the first task. Then change the starting value for T to 1, set rT = 10-6, dT = 0.01 and kT = 0. Keep the innate response turned off. Run the simulation for 200 days. You might need to change the plot to a log y-axis so you can see all variables. You will see that - in contrast to the Acute Virus model, the adaptive response is triggered. However, since we set its action to 0, it does not yet impact the virus and cell kinetics. Let’s change that. Set kT = 10-5. Run the simulation. You should see the adaptive response slowly coming in and clearing infected cells, which leads to a reduction in virus load and a recovery of the uninfected cells.
Record
Number of uninfected cells at the end of the simulation, non-zero kT
Number of infected cells at the end of the simulation, non-zero kT
Now we’ll turn on the immune response and run again. First with immune action at zero.
#call the simulator function with the specified settings
res <- simulate_chronicvirusir_ode(U = 1e+05, I = 0, V = 10, F = 0, T = 1,
n = 1e4, dU = 0.1,
dI = 1, dV = 4, b = 1e-06, p = 100, g = 1,
rF = 0, dF = 0, kF = 0,
rT = 1e-6, dT = 0.01, kT = 0,
tstart = 0,
tfinal = 200,
dt = 0.1
)Here is what we get from this.
generate_ggplot(list(res))
generate_text(list(res))Minimum / Maximum / Final value of U: 24698.06 / 1e+05 / 40404.04
Minimum / Maximum / Final value of I: 0.00 / 21121.91 / 5959.60
Minimum / Maximum / Final value of V: 4.76 / 517886.89 / 147500.00
Minimum / Maximum / Final value of F: 0.00 / 0.00 / 0.00
Minimum / Maximum / Final value of T: 0.93 / 6.1e+11 / 6.1e+11
We see the adaptive response kick in, but it doesn’t act on the rest of the system, so no impact there.
Let’s re-do, now with the adaptive response killing infected cells.
#call the simulator function with the specified settings
res <- simulate_chronicvirusir_ode(U = 1e+05, I = 0, V = 10, F = 0, T = 1,
n = 1e4, dU = 0.1,
dI = 1, dV = 4, b = 1e-06, p = 100, g = 1,
rF = 0, dF = 0, kF = 0,
rT = 1e-6, dT = 0.01, kT = 1e-5,
tstart = 0,
tfinal = 200,
dt = 0.1
)generate_ggplot(list(res))
generate_text(list(res))Minimum / Maximum / Final value of U: 24699.07 / 1e+05 / 90538.25
Minimum / Maximum / Final value of I: 0.00 / 21121.33 / 419.75
Minimum / Maximum / Final value of V: 4.76 / 517872.45 / 10264.46
Minimum / Maximum / Final value of F: 0.00 / 0.00 / 0.00
Minimum / Maximum / Final value of T: 0.93 / 121517.96 / 121517.96
Now we do see an impact.
This again pulls out the numbers you are asked to report.
Ufinal <- round(tail(res$ts$U,1),0)
Ifinal <- round(tail(res$ts$I,1),0)Answer for Task 3
Number of uninfected cells at the end of the simulation, non-zero kT: 90538
Number of infected cells at the end of the simulation, non-zero kT: 420
Task 4
Explore how different strengths of the adaptive response impact the results. Keep all settings at the values you just had. Then run the model for kT = 10-2, 10-3, 10-4,… all the way to 10-12 (that’s 11 different values). Each time, run for 200 days and record final number of infected cells and virus. Then make two plots, both with the kT values on the x-axis and the first with infected cell on the y-axis, the second with virus numbers on the y-axis. From these plots, determine (approximately) the values of kT for which infected cell numbers are at 500, and the value at which virus load is at 20,000. Ideally, you would do this task writing a bit of R code to run a loop over the values and draw the plots. But of course you can also do it manually.
Record
- True or False: The value of kT at which virus load is 20,000 is greater than 1.5 x 10-10
This is where the approach using code instead of the GUI is really useful. We can fully automate this task by looping over the different parameter values, with these lines of code.
#vector of kt values
ktvec = 10^-seq(2,12,by=1)
#those will hold the outcomes of interest
Ifinalvec = rep(0,length(ktvec))
Vfinalvec = rep(0,length(ktvec))
for (n in 1:length(ktvec))
{
#run simulation for each kt value
res <- simulate_chronicvirusir_ode(U = 1e+05, I = 0, V = 10, F = 0, T = 1,
n = 1e4, dU = 0.1,
dI = 1, dV = 4, b = 1e-06, p = 100, g = 1,
rF = 0, dF = 0, kF = 0,
rT = 1e-6, dT = 0.01, kT = ktvec[n],
tstart = 0,
tfinal = 200,
dt = 0.1
)
#record final number of infected cells and virus for each kt value
Ifinalvec[n] <- tail(res$ts$I,1)
Vfinalvec[n] <- tail(res$ts$V,1)
}
plot(ktvec,Ifinalvec, type = "b",log = "x")
plot(ktvec,Vfinalvec, type = "b",log = "x")
We can do a zoomed-in plot to better see the area of interest.
plot(ktvec,Vfinalvec, type = "b", xlim = c(1e-10,1e-9), ylim = c(15000,25000))
abline(a=20000,b=0,lty = "dashed")
From this, we can guess that at a virus load of 20,000, the value for kT is greater than 1.5 x 10-10. We can’t be fully certain since we didn’t try a lot of values. We could of course just run the simulation at that value. Or we could get an overall nicer look at the curves by running the simulation for more kT values. That’s tedious by hand, but very easy with code, so let’s do it. The only change is how many values we chose, the rest of the code is exactly the same:
#many more kt values
ktvec = 10^-seq(2,12,length=100)
#those will hold the outcomes of interest
Ifinalvec = rep(0,length(ktvec))
Vfinalvec = rep(0,length(ktvec))
for (n in 1:length(ktvec))
{
#run simulation for each kt value
res <- simulate_chronicvirusir_ode(U = 1e+05, I = 0, V = 10, F = 0, T = 1,
n = 1e4, dU = 0.1,
dI = 1, dV = 4, b = 1e-06, p = 100, g = 1,
rF = 0, dF = 0, kF = 0,
rT = 1e-6, dT = 0.01, kT = ktvec[n],
tstart = 0,
tfinal = 200,
dt = 0.1
)
#record final number of infected cells and virus for each kt value
Ifinalvec[n] <- tail(res$ts$I,1)
Vfinalvec[n] <- tail(res$ts$V,1)
}
plot(ktvec,Ifinalvec, type = "b",log = "x")
plot(ktvec,Vfinalvec, type = "b",log = "x")
plot(ktvec,Vfinalvec, type = "b", xlim = c(1e-10,5e-10), ylim = c(15000,25000))
abline(a=20000,b=0,lty = "dashed")
Answer for Task 4
True or False: The value of kT at which virus load is 20,000 is greater than 1.5 x 10-10: TRUE
Task 5
Now let’s add the innate response to the model. Keep everything as you just had it. That means U = 105, I = 0, V = 10, F = 0, T = 1, n = 104, dU = 0.1, dI = 1, dV = 4, b = 10-6, p = 100, g = 1, rT = 10-6, dT = 0.01 and kT = 10-5. The set F = 1, rF = 10, dF = 1. Start with kF = 0 and run the simulation for 200 days. You see the inate response coming up, but since it doesn’t have an effect yet, the virus and cell kinetics doesn’t change from the results you got above in the absence of the innate response. Now re-run with kF = 10-5. You should see that the peak virus load is lower. Maybe surprisingly, the virus load at the end is actually higher. The reason is that the reduced virus load leads to a slower induction of the adaptive response, and the weaker adaptive response in turn does not clear infected cells as well (by day 200, if you run it longer the adaptive response will do its job).
Record
Peak virus load with kF = 0
Final virus load with with kF = 0
Peak virus load with kF = 10-5
Final virus load with with kF = 10-5
Now with the innate response present as well, first it’s not active
#call the simulator function with the specified settings
res <- simulate_chronicvirusir_ode(U = 1e+05, I = 0, V = 10, F = 1, T = 1,
n = 1e4, dU = 0.1,
dI = 1, dV = 4, b = 1e-06, p = 100, g = 1,
rF = 10, dF = 1, kF = 0,
rT = 1e-6, dT = 0.01, kT = 1e-5,
tstart = 0,
tfinal = 200,
dt = 0.1
)generate_ggplot(list(res))
generate_text(list(res))Minimum / Maximum / Final value of U: 24699.07 / 1e+05 / 90538.25
Minimum / Maximum / Final value of I: 0.00 / 21121.33 / 419.75
Minimum / Maximum / Final value of V: 4.76 / 517872.45 / 10264.46
Minimum / Maximum / Final value of F: 0.95 / 190129.97 / 4202.67
Minimum / Maximum / Final value of T: 0.93 / 121517.96 / 121517.96
Here are the values we want
Vmax1 <- round(max(res$ts$V),0)
Vfinal1 <- round(tail(res$ts$V,1),0)Now with the innate response present and active
#call the simulator function with the specified settings
res <- simulate_chronicvirusir_ode(U = 1e+05, I = 0, V = 10, F = 1, T = 1,
n = 1e4, dU = 0.1,
dI = 1, dV = 4, b = 1e-06, p = 100, g = 1,
rF = 10, dF = 1, kF = 1e-5,
rT = 1e-6, dT = 0.01, kT = 1e-5,
tstart = 0,
tfinal = 200,
dt = 0.1
)generate_ggplot(list(res))
generate_text(list(res))Minimum / Maximum / Final value of U: 57302.53 / 1e+05 / 72758.35
Minimum / Maximum / Final value of I: 0.00 / 8854.72 / 1507.54
Minimum / Maximum / Final value of V: 4.76 / 124480.69 / 32227.72
Minimum / Maximum / Final value of F: 0.95 / 83471.75 / 15447.63
Minimum / Maximum / Final value of T: 0.93 / 57941.92 / 57941.92
Again pulling out the values we want
Vmax2 <- round(max(res$ts$V),0)
Vfinal2 <- round(tail(res$ts$V,1),0)Answer for Task 5
Peak virus load with kF = 0: 517872
Final virus load with with kF = 0: 10264
Peak virus load with kF = 10-5: 124481
Final virus load with with kF = 10-5: 32228
Task 6
We finish with one more systematic exploration. Keep all settings as in the previous task, but change simulation time to 300 days. Run the model for these kF values: 10-8, 10-7, … , 10-3 (6 values). For each value, record peak and final virus load. Repeat all of these runs, but now with rT = 10-7. Make a plot with the kF values on the x-axis and the peak virus values on the y-axis, for both levels of rT. Then make a second plot, now with final virus load on the y-axis.
Record
- True or False: The final virus load decreases monotonically as kF increases, for both values of rT
This is easily done using code. There are fancier ways of doing things, but a conceptual easy one is to have neste loops. We’ll run an outer loop over the two values for rT and then for each value of that parameter, we run an inner loop over the kF values. Here is the code:
#rt values
rtvec = c(1e-6,1e-7)
#kt values
kfvec = 10^-seq(8,3,by = -1)
#this time we'll store everything in a data frame
allres = data.frame(kfval = double(), rtval = double(), Vpeak = double(), Vfinal = double())
ct = 1 # a counter for easy indexing
#outer loop
for (i in 1:length(rtvec))
{
#inner loop
for (n in 1:length(ktvec))
{
#run simulation for each kt and rt value
res <- simulate_chronicvirusir_ode(U = 1e+05, I = 0, V = 10, F = 1, T = 1,
n = 1e4, dU = 0.1,
dI = 1, dV = 4, b = 1e-06, p = 100, g = 1,
rF = 10, dF = 1, kF = kfvec[n],
rT = rtvec[i], dT = 0.01, kT = 1e-5,
tstart = 0,
tfinal = 300,
dt = 0.1
)
#record final number of infected cells and virus for each kt value
Vfinal <- tail(res$ts$V,1)
Vpeak <- max(res$ts$V)
allres[ct, ] = c(kfvec[n],rtvec[i],Vpeak,Vfinal)
ct = ct + 1 #increase counter
} #end inner loop
} #end outer loop
#make a plot with ggplot
p1 <- ggplot(allres) +
geom_line(aes(x=kfval,y=Vpeak,color=as.factor(rtval))) +
scale_x_log10()
plot(p1)
p2 <- ggplot(allres) +
geom_line(aes(x=kfval,y=Vfinal,color=as.factor(rtval))) +
scale_x_log10()
plot(p2)
The curves suggest that the virus peak decreases monotonically as kF increases, and in fact the 2 curves seem to be on top of each other. But the final virus load does not and has an intermediate peak for rT = 10-6.
This is an example where a few data points don’t give the best idea of the shape of the curve. Fortunately, we can easily run this again for a lot more values for kF.
#rt values
rtvec = c(1e-6,1e-7)
#kt values
kfvec = 10^-seq(8,3,length = 100)
#this time we'll store everything in a data frame
allres = data.frame(kfval = double(), rtval = double(), Vpeak = double(), Vfinal = double())
ct = 1 # a counter for easy indexing
#outer loop
for (i in 1:length(rtvec))
{
#inner loop
for (n in 1:length(ktvec))
{
#run simulation for each kt and rt value
res <- simulate_chronicvirusir_ode(U = 1e+05, I = 0, V = 10, F = 1, T = 1,
n = 1e4, dU = 0.1,
dI = 1, dV = 4, b = 1e-06, p = 100, g = 1,
rF = 10, dF = 1, kF = kfvec[n],
rT = rtvec[i], dT = 0.01, kT = 1e-5,
tstart = 0,
tfinal = 300,
dt = 0.1
)
#record final number of infected cells and virus for each kt value
Vfinal <- tail(res$ts$V,1)
Vpeak <- max(res$ts$V)
allres[ct, ] = c(kfvec[n],rtvec[i],Vpeak,Vfinal)
ct = ct + 1 #increase counter
} #end inner loop
} #end outer loop
#make a plot with ggplot
p1 <- ggplot(allres) +
geom_line(aes(x=kfval,y=Vpeak,color=as.factor(rtval))) +
scale_x_log10()
plot(p1)
p2 <- ggplot(allres) +
geom_line(aes(x=kfval,y=Vfinal,color=as.factor(rtval))) +
scale_x_log10()
plot(p2)
These curves look better and confirm that the final virus load does not decrease monotonically for one of the two rT values. To understand why this pattern happens, it might be worth actually looking at the time-series to try and understand what exactly is going on. I’m not doing that here, but I’m sure by now you know how to do that. One approach is to save all results returned as res for each run in one big list or data frame, then plot (some of) them.
Answer for Task 6
True or False: The final virus load decreases monotonically as kF increases, for both values of rT: FALSE
Task 7
Continue to explore how different strengths of innate and adaptive response kinetics impact the virus and infected cell dynamics, both during the acute and the chronic phase.
Record
- Nothing
Nothing to report for this task, just an open exploration.
Answer for Task 7
Nothing