tci (version 0.2.0)

clc: Simulate closed-loop control

Description

Simulate closed-loop control using Bayesian updates. Infusion rates are calculated using `pkmod_prior` to reach `target_vals` at `target_tms`. Data values are simulated using `pkmod_true` at `obs_tms`. `pkmod_prior` and `pkmod_true` do not need to have the same structure. Model parameters are updated at each update time using all available simulated observations. Processing delays can be added through the `delay` argument, such that observations aren't made available to the update mechanism until `update_tms >= obs_tms + delay`.

Usage

clc(
  pkmod_prior,
  pkmod_true,
  target_vals,
  target_tms,
  obs_tms,
  update_tms,
  type = c("effect", "plasma"),
  custom_alg = NULL,
  resp_bounds = NULL,
  delay = 0,
  seed = NULL
)

Arguments

pkmod_prior

`pkmod` object describing a PK/PK-PD model that is used to calculate TCI infusion rates and is updated as data are simulated and incorporated. Must have an associated Omega matrix.

pkmod_true

`pkmod` object describing the patient's "true" response. This model will be used to simulate observations.

target_vals

A vector of numeric values indicating PK or PD targets for TCI algorithm.

target_tms

A vector of numeric values indicating times at which the TCI algorithm should begin targeting each value.

obs_tms

Times at which data values should be simulated from `pkmod_true`.

update_tms

Times at which `pkmod_prior` should be updated using all available simulated observations.

type

Type of TCI algorithm to be used. Options are "plasma" and "effect". Defaults to "effect". Will be overwritten if `custom_alg` is non-null.

custom_alg

Custom TCI algorithm to overwrite default plasma- or effect-site targeting.

resp_bounds

Optional vector of two values indicating minimum and maximum values possible for the response.

delay

Optional numeric value indicating a temporal delay between when observations are simulated and when they should be made available for updating `pkmod_prior`. For example, a delay should be set to account for a processing time delay in Bispectral Index measurements or the time required to measure drug concentrations from collected samples.

seed

An integer used to initialize the random number generator.

Examples

prior_vcov <- matrix(diag(c(0.265,0.346,0.209,0.610,0.565,0.597,0.702,0.463)),
8,8, dimnames = list(NULL,c('cl','q2','q3','v','v2','v3','ke0','sigma_add')))
pkmod_prior <- pkmod(pars_pk = c(cl = 10, q2 = 2, q3 =20, v = 15, v2 = 30, v3 = 50, ke0 = 1.2),
sigma_add = 0.2, log_response = TRUE, Omega = prior_vcov)
pkmod_true  <- pkmod(pars_pk = c(cl = 16, q2 = 4, q3 =10, v = 20, v2 = 20, v3 = 80, ke0 = 0.8),
sigma_add = 0.1, log_response = TRUE)
target_vals <- c(2,3,4,3,3)
target_tms <- c(0,5,10,36,60)
obs_tms <- c(1,2,4,8,12,16,24,36,48)
update_tms <- c(5,15,25,40,50)
sim <- clc(pkmod_prior, pkmod_true, target_vals, target_tms, obs_tms, update_tms, seed = 200)
len <- 500
tms <- seq(0,60,length.out = len)
# true, prior, posterior plasma concentrations
df <- data.frame(time = rep(tms,3),
                 value = c(predict(pkmod_true, sim$inf,tms)[,1],
                 predict(pkmod_prior, sim$inf,tms)[,1],
                 predict(sim$pkmod_post, sim$inf, tms)[,1]),
                 type = c(rep("true",len),rep("prior",len),rep("posterior",len)))
library(ggplot2)
ggplot(df, aes(x = time, y = value, color = type)) +
  geom_line() +
  geom_point(data = sim$obs, aes(x = time, y = obs), inherit.aes = FALSE) +
  geom_step(data = data.frame(time = target_tms, value = target_vals),
  aes(x = time, y = value), inherit.aes = FALSE)


# PK-PD example with observation delay (30 sec)
prior_vcov <- matrix(diag(c(0.265,0.346,0.209,0.702,0.242,0.230)),6,6,
dimnames = list(NULL,c('cl','q2','q3','ke0','c50','sigma_add')))
pkpdmod_prior <- update(pkmod_prior, pars_pd = c(c50 = 2.8, gamma = 1.47, e0 = 93, emx = 93),
pdfn = emax, pdinv = emax_inv, sigma_add = 4, log_response = FALSE, Omega = prior_vcov)
pkpdmod_true <- update(pkmod_true, pars_pd = c(c50 = 3.4, gamma = 1.47, e0 = 93, emx = 93),
pdfn = emax, pdinv = emax_inv, sigma_add = 3, log_response = FALSE)
target_vals <- c(75,60,50,50)
target_tms <- c(0,3,6,10)
obs_tms <- seq(1/6,10,1/6)
update_tms <- seq(1,10,0.5)
if (FALSE) {
sim_pkpd <- clc(pkpdmod_prior, pkpdmod_true, target_vals, target_tms, obs_tms,
update_tms, seed = 201, delay = 0.5)
# plot results
tms <- seq(0,10,length.out = len)
df <- data.frame(time = rep(tms,3),
                 value = c(predict(pkpdmod_true, sim_pkpd$inf,tms)[,"pdresp"],
                 predict(pkpdmod_prior, sim_pkpd$inf,tms)[,"pdresp"],
                 predict(sim_pkpd$pkmod_post, sim_pkpd$inf, tms)[,"pdresp"]),
                 type = c(rep("true",len),rep("prior",len),rep("posterior",len)))
library(ggplot2)
ggplot(df, aes(x = time, y = value, color = type)) +
  geom_line() +
  geom_point(data = sim_pkpd$obs, aes(x = time, y = obs), inherit.aes = FALSE) +
  labs(x = "Hours", y = "Bispectral Index") + theme_bw() +
  geom_vline(xintercept = update_tms, linetype = "dotted", alpha = 0.6) +
  geom_step(data = data.frame(time = target_tms, value = target_vals),
  aes(x = time, y = value), inherit.aes = FALSE)
}