(7) Incorporating Events into a Simulation

Events

The MCSimMod package allows one to solve ordinary differential equation (ODE) initial value problems in which discrete changes in state variables may occur at specific points in time. We call these discrete changes in state variables “events,” which is the same term used by the authors of the deSolve package. The MCSimMod package uses functions from the deSolve package to solve ODE initial value problems, and these functions have a built-in mechanism for handling events. One can learn more about deSolve event handling in the deSolve package documentation. Here we will demonstrate how to perform a simulation that involves a few simple events.

A Classical Pharmacokinetic Model

A classical PK model

In this ODE model, there are three “compartments” (represented as rectangular boxes in the figure) for which amounts of substance are tracked as state variables: \(A_0\) represents the amount of substance in the “exposure compartment” (e.g., the gut in an oral dosing scenario); \(A_1\) represents the amount of substance in the “central compartment” (e.g., in the blood of the organism); and \(A_2\) represents the amount that has been cleared from the body (e.g., via metabolism or excretion). Thus, \[\begin{align} \frac{\textrm{d}}{\textrm{d}t}A_0(t) &= -k_{01} A_0(t), \\ \frac{\textrm{d}}{\textrm{d}t}A_1(t) &= k_{01} A_0(t) - k_{12} A_1(t), \textrm{ and} \\ \frac{\textrm{d}}{\textrm{d}t}A_2(t) &= k_{12} A_1(t), \\ \end{align}\] where \(k_{01}\) is the rate constant (with units of one over time) that, along with the value of \(A_0\), determines the rate of delivery of the substance to the body and \(k_{12}\) is the rate constant (with units of one over time) that, along with the value of \(A_1\), determines the rate of clearance from the body. In addition to the state variables, the model includes an “output” variable, \(C\), that represents the concentration in the central compartment. The value of this variable is computed as \[\begin{equation} C(t) = \frac{A_1(t)}{V_\textrm{d}}, \end{equation}\] where \(V_\textrm{d}\) is a parameter representing the effective volume of the central compartment, or the “volume of distribution.” The model includes another output variable, \[\begin{equation} A_\textrm{tot} = A_0 + A_1 + A_2, \end{equation}\] that represents the total amount of substance in the system. Finally, the model includes an additional state variable that represents the area under the concentration vs. time curve (AUC), as this is a quantity that is often computed in pharmacokinetic analyses. The AUC for the period from time \(0\) to time \(t\) is given by \(\textrm{AUC}(t) = \int_0^t C(\tau) \, \textrm{d}\tau\) and thus the state equation for this quantity is \[\begin{equation} \frac{\textrm{d}}{\textrm{d}t} \textrm{AUC}(t) = C(t). \end{equation}\]

In order to solve an initial value problem for the classical PK model, one needs to provide the values of the three parameters (\(k_{01}\), \(k_{12}\), and \(V_\textrm{d}\)) as well as the initial values of the four state variables (\(A_0\), \(A_1\), \(A_2\), and \(\textrm{AUC}\)). One can also provide a data structure that describes discrete changes to state variables at specific points in time. Details about how to define and use an “events” data structure are provided below.

MCSim Model Specification

We used the GNU MCSim model specification language to implement the Newton’s Law of Cooling model. The complete MCSim model specification file for this model, pk1.model, can be found in the extdata subdirectory of the MCSimMod package.

The model specification file uses the text symbols A0, A1, A2, and AUC to represent the state variables \(A_0\), \(A_1\), \(A_2\), and \(\textrm{AUC}\), respectively, and the text symbols k01, k12, and Vd to represent the parameters \(k_{01}\), \(k_{12}\), and \(V_\textrm{d}\), respectively. In addition, the text symbols A0_init, A1_init, A2_init, and AUC_init represent parameters that can be used to set (via the updateY0() method of the Model class) the initial conditions of the state variables.

Building the Model

First, we load the MCSimMod package as follows.

library(MCSimMod)
#> Loading required package: tools

Using the following commands, we create a model object (i.e., an instance of the Model class) using the model specification file pk1.model that is included in the MCSimMod package.

# Get the full name of the package directory that contains the example MCSim
# model specification file.
mod_path <- file.path(system.file(package = "MCSimMod"), "extdata")

# Create a model object using the example MCSim model specification file
# "pk1.model" included in the MCSimMod package.
pk1_mod_name <- file.path(mod_path, "pk1")
pk1_mod <- createModel(pk1_mod_name)

Once the model object is created, we can “load” the model (so that it’s ready for use in a given R session) as follows.

# Load the model.
pk1_mod$loadModel()
#> C compilation complete. Full details are available in the file C:\Users\dkapraun\AppData\Local\Temp\Rtmp4YuNTn\compiler_output.txt.
#> Hash created and saved in the file C:\Users\dkapraun\AppData\Local\Temp\Rtmp4YuNTn\pk1_model.md5.

Predicting the Blood Concentration of a Substance During Repeated Oral Dosing using Events

Suppose we want to predict the blood concentration of a substance during a 48-hour period when oral doses of 50 mg are provided every 12 hours. Let’s assume that the oral absorption rate constant (\(k_{01}\)) is 1.0 h\(^{-1}\), the clearance rate constant (\(k_{12}\)) is 0.5 h\(^{-1}\), and the volume of distribution (\(V_\textrm{d}\)) is 50 L. These are the default values of those model parameters given in the model specification file. We can change the initial condition for the state variable \(A_0\) to zero using the following commands. (Note that all other state variables have default initial conditions of zero based on the model specification file.)

# Change the values of the model parameters from their default values.
pk1_mod$updateParms(c(A0_init = 0))

# Update the initial value(s) of the state variable(s) based on the updated
# parameter value(s).
pk1_mod$updateY0()

We can verify the values of all parameters and initial conditions with the following commands.

pk1_mod$parms
#>  A0_init  A1_init  A2_init AUC_init       Vd      k01      k12 
#>      0.0      0.0      0.0      0.0     50.0      1.0      0.5
pk1_mod$Y0
#>  A0  A1  A2 AUC 
#>   0   0   0   0

To specify events, we first create an R data frame with four columns: var, which contains the names of the state variables to be modified in each event; time, which contains the times of each event; value, which contains the amounts by which the state variables are to be changed; and method, which contains descriptions of the methods by which the state variables are to be changed. As explained in the deSolve documentation, the event “method” can be “replace,” “add,” or “multiply.” We can create an appropriate data frame for our simulation (and examine it) using the following commands.

events_df <- data.frame(
  var = "A0", time = seq(from = 0, to = 48, by = 12), value = 50,
  method = "add"
)
head(events_df)
#>   var time value method
#> 1  A0    0    50    add
#> 2  A0   12    50    add
#> 3  A0   24    50    add
#> 4  A0   36    50    add
#> 5  A0   48    50    add

Finally, we can perform a simulation that provides results for the desired output times (i.e., \(t = 0, 0.01, 0.02, \ldots, 48.00\)) using the following commands.

# Define output times for simulation.
times <- seq(from = 0, to = 48, by = 0.01)

# Run simulation.
out <- pk1_mod$runModel(times, events = list(data = events_df))

Examining the Results

The final command shown above performs a model simulation (incorporating events) and stores the simulation results in a “matrix” data structure called out. There is one row for each output time, and one column for each state variable and each output variable. Five rows of this data structure are shown below. Note that the independent variable, which is \(t\) in the case of the classical PK model, is always labeled “time” in the output data structure.

Note also that each event (in this case, adding 50 to the state variable \(A_0\) every 12 hours) occurs immediately after the specified “time” of the event (e.g, at 12 hours into the simulation).

We can create a visual representation of the simulation results using the following commands.

# Plot simulation results.
plot(out[, "time"], out[, "C"],
  type = "l", lty = 1, lwd = 2,
  xlab = "Time (h)", ylab = "Concentration (mg/L)"
)