knitr::opts_chunk$set( collapse = TRUE, comment = "#>", echo = TRUE, message = FALSE, warning = FALSE, fig.align="center", fig.height= 4, fig.width = 6 )
Target-controlled infusion (TCI) systems calculate infusion rates required to reach target concentrations or effects within a patient. Where pharmacokinetic (PK) and pharmacodynamic (PD) models describe a time course of concentrations or effects, respectively, associated with a series of doses, TCI algorithms calculate the inverse relationship: what doses must be administered to achieve target responses?
tci package implements TCI algorithms for PK and PK-PD models for drugs described by compartmental models and administered via intravenous infusion. The package provides closed-form solutions for one, two, or three compartment mammillary models (i.e., all peripheral compartments are joined to a central compartment), as well as a three-compartment model with an adjoining effect-site. PK model code is based on solutions published by @Abuhelwa2015 and models are implemented in C++ via
Rcpp. TCI algorithms for plasma- and effect-site targeting are implemented based on work by @Jacobs1990 and @Shafer1992, respectively. Users can specify alternative PK models or TCI algorithms. See the
custom vignette for further details.
library(tci) library(ggplot2) # ggplot for plotting library(gridExtra) # arrangeGrob to arrange plots library(reshape2) # melt function
old <- theme_set(theme_bw()) ggplot <- function(...) ggplot2::ggplot(...) + scale_color_brewer(palette="Pastel1") # scale_color_manual(values = c("black","steelblue","seagreen"))
tci package is built around S3 classes
poppkmod, created with the functions
pkmod objects serve as containers for 1) functions implementing the structural PK model (e.g., a 1-compartment model with first-order elimination) and the PD model, if applicable, 2) the parameters for the respective functions, and 3) initial concentrations, and 4) information relevant for simulating observations or implementing TCI control, such as the compartment number associated with observations or with an effect-site.
poppkmod are wrapper objects that contain one or more
pkmod objects associated with published population PK models: the Marsh, Schnider, and Eleveld models for propofol, and the Minto, Kim, and Eleveld models for remifentanil.
poppkmod objects have associated
simulate methods that can be used to predict concentrations and simulate observations (PK or PD) given an infusion schedule. Infusion schedules, in turn, are created either manually via
inf_manual or by applying a TCI algorithm to reach designated targets via
pkmod objects are additionally equipped with a
update method that allows for model components (e.g., parameter values, initial concentrations) to be easily modified. Both
simulate methods pass additional arguments via the ellipses argument,
update.pkmod to readily allow for prediction or simulation under different conditions.
Examples in this vignette will focus on illustrating the lower-level functions of
tci applied to
pkmod objects. See the vignette on population PK models for illustration of higher-level functions and applications to population PK models for propofol and remifentanil.
Equations implementing 1-,2-,3-compartment and 3-compartment-effect structural PK models are included in the
tci package. The function
pkmod will automatically infer the correct structure based on the parameter names.
# 1-compartment model (mod1cpt <- pkmod(pars_pk = c(cl = 10, v = 15))) # 3-compartment model with effect site (mod3ecpt <- pkmod(pars_pk = c(cl = 10, q2 = 2, q3 =20, v = 15, v2 = 30, v3 = 50, ke0 = 1.2)))
Acceptable parameter names can be viewed by calling
list_parnms(). Less-commonly used parameters, such as clearance from a peripheral compartment, are also permissible.
# acceptable parameter names list_parnms()
pkmod objects can be modified through an
update.pkmod method. Perhaps most usefully, this allows for partial modifications to PK-PD parameters. For example, the effect-site equilibrium constant can be easily updated.
update(mod3ecpt, pars_pk = c(ke0 = 0.9), init = c(1,0.2,0.3,1))
Most functions in the
tci package pass additional arguments to
update.pkmod allowing for easy modification of
pkmod objects as needed.
An infusion schedule is required to for
simulate methods. This schedule should be a matrix with column labels "begin", "end", and "infrt", indicating infusion begin times, end times, and infusion rates. It can be created directly by the user, or outputted by the
inf_tci functions. In the former function, the user specifies infusion start times, durations, and infusion rates.
# single infusion (single_inf <- inf_manual(inf_tms = 0, duration = 0.5, inf_rate = 100)) # multiple infusions (multi_inf <- inf_manual(inf_tms = c(0,3,6), duration = c(1,0.5,0.25), inf_rate = 100))
Typically, however, the
inf_tci will be used to calculate infusion rates required to reach specified targets.
inf_tci requires 1) a set of target concentrations (or PD response values) and corresponding times at which the target is set, and 2) a
pkmod object. It has "plasma" and "effect" settings, implementing the Jacobs and Shafer algorithms, respectively. Custom algorithms can be specified through the
custom_alg argument. See the vignette on custom models and algorithms for more details.
# plasma targeting for one-compartment model inf_1cpt <- inf_tci(target_vals = c(2,3,4,4), target_tms = c(0,2,3,10), pkmod = mod1cpt, type = "plasma") head(inf_1cpt) # effect-site targeting for three-compartment effect site model inf_3ecpt <- inf_tci(target_vals = c(2,3,4,4), target_tms = c(0,2,3,10), pkmod = mod3ecpt, type = "effect") head(inf_3ecpt)
By default, plasma- and effect-targeting algorithms are updated in increments of 1/6. If a PK model elimination parameters have units of minutes (as do commonly used models for the anesthetic propofol), this will correspond to updating TCI targets at 10-second intervals. If elimination rates are in different units, such as hours, then the TCI update frequency should be modified by the argument
The infusion schedule can be applied to the
pkmod object using
simulate.pkmod methods to predict concentrations or simulate observations, respectively. Using the three-compartment model as illustration
# prediction/observation times tms_pred <- seq(0,10,0.01) tms_obs <- c(0.5,1,2,4,6,10) pre <- predict(mod3ecpt, inf = inf_3ecpt, tms = tms_pred) obs <- simulate(mod3ecpt, seed = 1, inf = inf_3ecpt, tms = tms_obs, sigma_mult = 0.2) # plot results dat <- data.frame(time = tms_pred, `plasma (3 cmpt)` = pre[,"c1"], `effect (ke0=1.2)` = pre[,"c4"], check.names = FALSE) datm <- melt(dat, id = "time") dat_obs <- data.frame(time = tms_obs, con = obs, variable = "plasma (3 cmpt)") p <- ggplot(datm, aes(x = time, y = value, color = variable)) + geom_line() + geom_point(data = dat_obs, aes(x = time, y = con)) + xlab("Minutes") + ylab("Concentration (mg/L)") p
pkmod object used in the predict and simulate methods does not need to be the same as the one used to calculate the infusion schedule. This permits the user to evaluate the effect of model misspecification either 1) by passing different parameter values to
simulate.pkmod, or 2) by using a different
To illustrate the parameter misspecification, we can evaluate predictions with a new effect-site equilibrium constant.
# evaluate with different ke0 parameter pre_misspec <- predict(mod3ecpt, inf = inf_3ecpt, tms = tms_pred, pars_pk = c(ke0 = 0.8)) dat_misspec <- data.frame(pre_misspec, variable = "effect (ke0=0.8)", time = tms_pred) p + geom_line(data = dat_misspec, aes(x = time, y = c4, color = variable))
To illustrate structural model misspecification, we can consider the case where PK are described by a one-compartment model, but infusions were calculated according to a three-compartment model.
# predicted concentrations pre_1cpt <- predict(mod1cpt, inf = inf_3ecpt, tms = tms_pred) dat_1cpt <- data.frame(pre_1cpt, variable = "plasma (1 cmpt)", time = tms_pred) # simulated observations obs_1cpt <- simulate(mod1cpt, seed = 1, inf = inf_3ecpt, tms = tms_obs, sigma_mult = 0.2) p + geom_line(data = dat_1cpt, aes(x = time, y = c1, color = variable)) + geom_point(data = data.frame(time = tms_obs, con = obs_1cpt, variable = "plasma (1 cmpt)"), aes(x = time, y = con), inherit.aes = FALSE, color = "green4")
All of the functions in
tci can be extended to include pharmacodynamic (PD) models. Unlike PK models, the equations describing PD models are typically invertible, allowing one to readily calculate the target effect-site concentration associated with a desired effect. The user, therefore, supplies to a
pkmod functions implementing the PD response (i.e., compute response from concentrations), and its inverse (i.e., concentrations from a response), as well as associated parameter values.
Four-parameter E-max models are commonly used to describe PD responses and are implemented in
tci. E-max models describe a response in terms of its minimum and maximum values,
e0, respectively, the concentration associated with 50\% effect,
c50, and the slope of the dose-response curve at c50,
gamma. In anesthesia, the Bispectral Index (BIS) is a commonly used measurement of a patient's depth of hypnosis and is often described by an E-max model. BIS is derived from EEG measurements and calibrated to vary between BIS=100, indicating a fully-alert state, and BIS=0, in which little brain activity is registered. BIS values between 40 and 60 typically indicate that a patient is sufficiently sedated for general anesthesia.
modpd <- update(mod3ecpt, pdfn = emax, pdinv = emax_inv, pars_pd = c(e0 = 100, emx = 100, c50 = 3.5, gamma = 2.2))
PD targets are passed along with the updated
inf_tci, which will assume values are PD values (unless overridden by the
ignore_pd argument of
inf_pd <- inf_tci(target_vals = c(70,60,50,50), target_tms = c(0,2,3,10), pkmod = modpd, type = "effect")
We can then similarly use
simulate.pkmod methods to predict and simulate PD responses. BIS measurements may be collected at a rate of one observation per 10-20 seconds, depending on the BIS device settings.
# predict responses pre_pd <- predict(modpd, inf = inf_pd, tms = tms_pred) # pd observations: 10 sec = 1/6 min tms_pd_obs <- seq(1/6,10,1/6) # simulate responses with additive error and parameter misspecification obs_pd <- simulate(modpd, seed = 1, inf = inf_pd, tms = tms_pd_obs, sigma_add = 5, pars_pk = c(ke0 = 0.7), pars_pd = c(c50 = 3, gamma = 1.8)) # plot results dat_pd <- data.frame(time = tms_pred, `plasma (3 cmpt)` = pre_pd[,"c1"], `effect (ke0=1.2)` = pre_pd[,"c4"], BIS = pre_pd[,"pdresp"], check.names = FALSE) dat_pdm <- melt(dat_pd, id = "time") dat_pdm$type <- as.factor(ifelse(dat_pdm$variable == "BIS", "PD","PK")) dat_pd_obs <- data.frame(time = tms_pd_obs, BIS = obs_pd, type = factor("PD"), variable = "BIS") levels(dat_pdm$type) <- levels(dat_pd_obs$type) <- c("Bispectral Index", "Concentration (mg/L)") ggplot(dat_pdm, aes(x = time, y = value, color = variable)) + facet_wrap(type~., scales = "free", nrow = 2) + geom_line() + geom_point(data = dat_pd_obs, aes(x = time, y = BIS)) + xlab("Minutes") + ylab("")
Simulations with potential model misspecification are most easily implemented using the function
simulate_tci which can be used for both
poppkmod classes. Required arguments to
simulate_tci are 1) a prior PK model (
pkmod_prior) that is used to calculate infusion rates and may be updated throughout the simulation if update times are provided, 2) a true PK model (
pkmod_true) that is used to simulate observations, 3) TCI target values, 4) TCI target times, and 5) times to simulate observations. If update times are specified then Bayesian updates will be performed to update parameters based on the (simulated) data available at each time. Data processing delays can be incorporated through the argument
To illustrate open-loop control, we simulate PK responses from a three-compartment model at times 1, 2, 3, 4, 8, and 12 over a 24 hour period in which effect-site targeting is used and the target concentration is raised from 2 mg/L to 4 mg/L.
mod_true <- update(mod3ecpt, pars_pk = c(cl = 20, q2 = 1.5, ke0 = 1.8)) sim_ol <- simulate_tci(pkmod_prior = mod3ecpt, pkmod_true = mod_true, target_vals = c(2,3,4,4), target_tms = c(0,2,3,24), obs_tms = c(1,2,3,4,8,12), seed = 1) ggplot(melt(sim_ol$resp, id.vars = c("time","type"))) + geom_line(aes(x = time, y = value, color = variable)) + geom_point(data = sim_ol$obs, aes(x = time, y = obs)) + facet_wrap(~type) + labs(x = "Hours", y = "Concentration (mg/L)")
Closed-loop control is implemented by specifying a set of update times. For model parameters to be updated,
pkmod_prior must have an "Omega" matrix specifying the variability in each parameter. This matrix is used as the prior variance-covariance matrix in the updates, while the prior model parameters are used as the prior point estimates.
Using the example above, we simulate samples drawn at 1, 2, 4, and 8 hours, with a processing time of 4 hours for each sample.
mod3ecpt <- update(mod3ecpt, sigma_mult = 0.2, Omega = matrix(diag(c(1.2,0.6,1.5,0.05)), 4,4, dimnames = list(NULL, c("cl","q2","v","ke0")))) sim_cl <- simulate_tci(pkmod_prior = mod3ecpt, pkmod_true = mod_true, target_vals = c(2,3,4,4), target_tms = c(0,2,3,24), obs_tms = c(1,2,3,4,8,12), update_tms = c(6,12,16), delay = 0, seed = 1) ggplot(melt(sim_cl$resp, id.vars = c("time","type"))) + geom_line(aes(x = time, y = value, color = variable)) + geom_point(data = sim_cl$obs, aes(x = time, y = obs)) + facet_wrap(~type) + labs(x = "Hours", y = "Concentration (mg/L)")
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.