/
simulate_basicvirus_modelexploration.R
126 lines (112 loc) · 6.36 KB
/
simulate_basicvirus_modelexploration.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
#' Simulation to illustrate parameter scan of the basic virus model
#'#'
#' @description This function simulates the basic virus model ODE for a range of parameters.
#' The function returns a data frame containing the parameter that has been varied and the outcomes (see details).
#'
#' @param U : Starting value for uninfected cells : numeric
#' @param I : Starting value for infected cells : numeric
#' @param V : Starting value for virus : numeric
#' @param n : Rate of new uninfected cell replenishment : numeric
#' @param dU : Rate at which uninfected cells die : numeric
#' @param dI : Rate at which infected cells die : numeric
#' @param dV : Rate at which virus is cleared : numeric
#' @param b : Rate at which virus infects cells : numeric
#' @param p : Rate at which infected cells produce virus : numeric
#' @param g : Possible conversion factor for virus units : numeric
#' @param tstart : Start time of simulation : numeric
#' @param tfinal : Final time of simulation : numeric
#' @param dt : Times for which result is returned : numeric
#' @param samples : Number of values to run between pmin and pmax : numeric
#' @param parmin : Lower value for varied parameter : numeric
#' @param parmax : Upper value for varied parameter : numeric
#' @param samplepar : Name of parameter to be varied : character
#' @param pardist : spacing of parameter values, can be either 'lin' or 'log' : character
#' @return The function returns the output as a list,
#' list element 'dat' contains the data frame with results of interest.
#' The first column is called xvals and contains the values of the
#' parameter that has been varied as specified by 'samplepar'.
#' The remaining columns contain peak and steady state values of bacteria and immune response,
#' Upeak, Ipeak, Vpeak, Usteady, Isteady and Vsteady.
#' A final boolean variable 'steady' is returned for each simulation.
#' It is TRUE if the simulation reached steady state, otherwise FALSE.
#' @details ##this code illustrates how to do analyze a simple model.
#' A simple 3 compartment ODE model (the basic virus model introduced in the app of that name)
#' is simulated for different parameter values.
#' This function runs the model for a range of values for any one parameter,
#' while holding all other paramter values fixed.
#' The user can specify which parameter is sampled, and
#' the simulation returns for each parameter sample the peak and final value for U, I and V.
#' Also returned is the varied parameter and an indicator if steady state was reached.
#' @section Notes: The parameter dt only determines for which times the solution is returned,
#' it is not the internal time step. The latter is set automatically by the ODE solver.
#' @section Warning: This function does not perform any error checking. So if
#' you try to do something nonsensical (e.g. specify negative parameter values
#' or fractions > 1), the code will likely abort with an error message.
#' @examples
#' # To run the simulation with default parameters just call the function:
#' \dontrun{res <- simulate_basicvirus_modelexploration()}
#' # To choose parameter values other than the standard one, specify them, like such:
#' res <- simulate_basicvirus_modelexploration(samples=5, samplepar='dI', parmin=1, parmax=10)
#' # You should then use the simulation result returned from the function, like this:
#' plot(res$dat[,"xvals"],res$data[,"Vpeak"],xlab='Parameter values',ylab='Virus Peak',type='l')
#' @seealso See the shiny app documentation corresponding to this simulator
#' function for more details on this model.
#' @author Andreas Handel
#' @export
simulate_basicvirus_modelexploration <- function(U = 1e+05, I = 0, V = 1, n = 10000, dU = 0.1, dI = 1, dV = 2, b = 2e-05, p = 5, g = 1, tstart = 0, tfinal = 100, dt = 0.1, samples = 10, parmin=1, parmax=10, samplepar='p', pardist = 'lin')
{
#initialize vectors that will contain the outcomes of interest
Upeak=rep(0,samples)
Ipeak=rep(0,samples)
Vpeak=rep(0,samples)
Usteady=rep(0,samples)
Isteady=rep(0,samples)
Vsteady=rep(0,samples)
#create values for the parameter of interest to sample over
#do equal spacing in log space
if (pardist == 'lin') {parvec=seq(parmin,parmax,length=samples)}
if (pardist == 'log') {parvec=10^seq(log10(parmin),log10(parmax),length=samples)}
steady = rep(TRUE,samples) #indicates if steady state has not been reached
for (nn in 1:samples)
{
#replace value of parameter we want to vary
if (samplepar == 'n') {n = parvec[nn]}
if (samplepar == 'dU') {dU = parvec[nn]}
if (samplepar == 'dI') {dI = parvec[nn]}
if (samplepar == 'dV') {dV = parvec[nn]}
if (samplepar == 'b') {b = parvec[nn]}
if (samplepar == 'p') {p = parvec[nn]}
if (samplepar == 'g') {g = parvec[nn]}
#this runs the bacteria ODE model for each parameter sample
#all other parameters remain fixed
odeout <- simulate_basicvirus_ode(U = U, I = I, V = V,
n = n, dU = dU, dI = dI, dV = dV,
b = b, p = p, g = g,
tstart = tstart, tfinal = tfinal, dt = dt)
timeseries = odeout$ts
#get peak and steady state for variables
Upeak[nn]=max(timeseries[,"U"]);
Ipeak[nn]=max(timeseries[,"I"]);
Vpeak[nn]=max(timeseries[,"V"]);
Usteady[nn] = utils::tail(timeseries[,"U"],1)
Isteady[nn] = utils::tail(timeseries[,"I"],1)
Vsteady[nn] = utils::tail(timeseries[,"V"],1)
#a quick check to make sure the system is at steady state,
#i.e. there is no absolute change larger than 0.1
#AND the value for at the final time is not more than
#1% different than several time steps earlier
#the abs() functions are there since sometimes for small values, they can become numerically negative
vl=nrow(timeseries);
finaldiff = abs(timeseries[vl,"U"]-timeseries[vl-10,"U"])
finalrel = finaldiff/abs(timeseries[vl,"U"])
if (finaldiff>0.1 || finalrel > 1e-2)
{
steady[nn] = FALSE
}
}
#final list structure containing all results that are returned
result = list()
dat = data.frame(xvals = parvec, Upeak = Upeak, Ipeak=Ipeak, Vpeak=Vpeak, Usteady = Usteady, Isteady = Isteady, Vsteady = Vsteady, steady = steady)
result$dat = dat
return(result)
}