---
title: "(7) Incorporating Events into a Simulation"
author: "Dustin Kapraun"
date: "`r Sys.Date()`"
output: html_vignette
vignette: >
  %\VignetteIndexEntry{(7) Incorporating Events into a Simulation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE,
  warning = FALSE,
  collapse = TRUE,
  comment = "#>"
)
```

## 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](https://CRAN.R-project.org/package=deSolve/deSolve.pdf). Here we will demonstrate how to perform a simulation that involves a few simple events.

## A Classical Pharmacokinetic Model

For this demonstration, we will use the classical pharmacokinetic (PK) model illustrated in the figure below.
```{r, echo=FALSE, fig.cap="A classical PK model", out.width = '100%'}
knitr::include_graphics("pk1.png")
```

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](https://www.gnu.org/software/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.
```{r}
library(MCSimMod)
```

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.
```{r, results='hide'}
# 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.
```{r, results='hide'}
# Load the model.
pk1_mod$loadModel()
```

## 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.)
```{r, results='hide'}
# 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.
```{r}
pk1_mod$parms
pk1_mod$Y0
```

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.
```{r}
events_df <- data.frame(
  var = "A0", time = seq(from = 0, to = 48, by = 12), value = 50,
  method = "add"
)
head(events_df)
```

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.
```{r, results='hide'}
# 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.
```{r, echo=FALSE, results='asis'}
library(knitr)
kable(out[1199:1203, ])
```
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.
```{r, fig.dim=c(6, 4), fig.align='center'}
# Plot simulation results.
plot(out[, "time"], out[, "C"],
  type = "l", lty = 1, lwd = 2,
  xlab = "Time (h)", ylab = "Concentration (mg/L)"
)
```
