clc: Simulate closed-loop control

In tci: Target Controlled Infusion (TCI)

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)
## Not run: 
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)

## End(Not run)

tci documentation built on Aug. 15, 2022, 9:09 a.m.