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?

The `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 Abuhelwa, Foster, and Upton (2015) and models
are implemented in C++ via `Rcpp`

. TCI algorithms for plasma-
and effect-site targeting are implemented based on work by Jacobs (1990) and Shafer
and Gregg (1992), 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
```

`pkmod`

and `poppkmod`

object classesThe `tci`

package is built around S3 classes
`pkmod`

and `poppkmod`

, created with the functions
`pkmod`

and `poppkmod`

, respectively.
`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.

Both `pkmod`

and `poppkmod`

objects have
associated `predict`

and `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 `inf_tci`

.
`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
`predict`

and `simulate`

methods pass additional
arguments via the ellipses argument, `...`

, to
`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
<- pkmod(pars_pk = c(cl = 10, v = 15)))
(mod1cpt #> tci pkmod object
#> See ?update.pkmod to modify or add elements
#>
#> PK model
#> 1-compartment PK model
#> PK parameters: cl = 10, v = 15
#> Initial concentrations: (0)
#> Plasma compartment: 1
#> Effect compartment: 1
#>
#> Simulation
#> Additive error SD: 0
#> Multiplicative error SD: 0
#> Logged response: FALSE
# 3-compartment model with effect site
<- pkmod(pars_pk = c(cl = 10, q2 = 2, q3 =20, v = 15, v2 = 30, v3 = 50, ke0 = 1.2)))
(mod3ecpt #> tci pkmod object
#> See ?update.pkmod to modify or add elements
#>
#> PK model
#> 4-compartment PK model
#> PK parameters: cl = 10, q2 = 2, q3 = 20, v = 15, v2 = 30, v3 = 50, ke0 = 1.2
#> Initial concentrations: (0,0,0,0)
#> Plasma compartment: 1
#> Effect compartment: 4
#>
#> Simulation
#> Additive error SD: 0
#> Multiplicative error SD: 0
#> Logged response: FALSE
```

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()
#> Acceptable names for 'pars_pk' vector (case-insensitive)
#>
#> First compartment options
#> Central volume: 'v','v1'
#> Elimination: 'cl','cl1','k10','ke'
#>
#> Second compartment options
#> Peripheral volume: 'v2'
#> Transfer: 'q','q2','k12','k21'
#> Elimination: 'cl2','k20'
#>
#> Third compartment options
#> Second peripheral volume: 'v3'
#> Transfer: 'q3','k13','k31'
#> Elimination: 'cl3','k30'
#>
#> Effect-site
#> Elimination: 'ke0'
```

Elements of `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))
#> tci pkmod object
#> See ?update.pkmod to modify or add elements
#>
#> PK model
#> 4-compartment PK model
#> PK parameters: cl = 10, q2 = 2, q3 = 20, v = 15, v2 = 30, v3 = 50, ke0 = 0.9
#> Initial concentrations: (1,0.2,0.3,1)
#> Plasma compartment: 1
#> Effect compartment: 4
#>
#> Simulation
#> Additive error SD: 0
#> Multiplicative error SD: 0
#> Logged response: FALSE
```

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 `predict`

and
`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_manual`

or
`inf_tci`

functions. In the former function, the user
specifies infusion start times, durations, and infusion rates.

```
# single infusion
<- inf_manual(inf_tms = 0, duration = 0.5, inf_rate = 100))
(single_inf #> begin end inf_rate
#> [1,] 0 0.5 100
# multiple infusions
<- inf_manual(inf_tms = c(0,3,6), duration = c(1,0.5,0.25), inf_rate = 100))
(multi_inf #> begin end inf_rate
#> [1,] 0.0 1.00 100
#> [2,] 1.0 3.00 0
#> [3,] 3.0 3.50 100
#> [4,] 3.5 6.00 0
#> [5,] 6.0 6.25 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_tci(target_vals = c(2,3,4,4), target_tms = c(0,2,3,10),
inf_1cpt pkmod = mod1cpt, type = "plasma")
head(inf_1cpt)
#> begin end inf_rate Ct c1_start c1_end
#> [1,] 0.00000 0.16667 190.1851 2 0 2
#> [2,] 0.16667 0.33333 20.0000 2 2 2
#> [3,] 0.33333 0.50000 20.0000 2 2 2
#> [4,] 0.50000 0.66667 20.0000 2 2 2
#> [5,] 0.66667 0.83333 20.0000 2 2 2
#> [6,] 0.83333 1.00000 20.0000 2 2 2
# effect-site targeting for three-compartment effect site model
<- inf_tci(target_vals = c(2,3,4,4), target_tms = c(0,2,3,10),
inf_3ecpt pkmod = mod3ecpt, type = "effect")
head(inf_3ecpt)
#> begin end inf_rate Ct c1_start c2_start c3_start c4_start
#> [1,] 0.00000 0.16667 643.6921 2 0.000000 0.00000000 0.0000000 0.0000000
#> [2,] 0.16667 0.33333 0.0000 2 6.033672 0.03532348 0.2079682 0.5966263
#> [3,] 0.33333 0.50000 0.0000 2 4.301795 0.09138590 0.5235366 1.4096671
#> [4,] 0.50000 0.66667 0.0000 2 3.136188 0.13103476 0.7267367 1.8178340
#> [5,] 0.66667 0.83333 19.3835 2 2.350071 0.15960595 0.8548449 1.9785480
#> [6,] 0.83333 1.00000 0.0000 2 2.000000 0.18174014 0.9390928 2.0109479
#> c1_end c2_end c3_end c4_end
#> [1,] 6.033672 0.03532348 0.2079682 0.5966263
#> [2,] 4.301795 0.09138590 0.5235366 1.4096671
#> [3,] 3.136188 0.13103476 0.7267367 1.8178340
#> [4,] 2.350071 0.15960595 0.8548449 1.9785480
#> [5,] 2.000000 0.18174014 0.9390928 2.0109479
#> [6,] 1.586605 0.19939667 0.9931813 1.9678507
```

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
`dtm`

.

The infusion schedule can be applied to the `pkmod`

object
using `predict.pkmod`

or `simulate.pkmod`

methods
to predict concentrations or simulate observations, respectively. Using
the three-compartment model as illustration

```
# prediction/observation times
<- seq(0,10,0.01)
tms_pred <- c(0.5,1,2,4,6,10)
tms_obs
<- predict(mod3ecpt, inf = inf_3ecpt, tms = tms_pred)
pre <- simulate(mod3ecpt, seed = 1, inf = inf_3ecpt, tms = tms_obs, sigma_mult = 0.2)
obs
# plot results
<- data.frame(time = tms_pred, `plasma (3 cmpt)` = pre[,"c1"],
dat `effect (ke0=1.2)` = pre[,"c4"],
check.names = FALSE)
<- melt(dat, id = "time")
datm <- data.frame(time = tms_obs, con = obs, variable = "plasma (3 cmpt)")
dat_obs
<- ggplot(datm, aes(x = time, y = value, color = variable)) +
p geom_line() +
geom_point(data = dat_obs, aes(x = time, y = con)) +
xlab("Minutes") + ylab("Concentration (mg/L)")
p
```

Notably, the `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 `update.pkmod`

via
`predict.pkmod`

or `simulate.pkmod`

, or 2) by
using a different `pkmod`

object.

To illustrate the parameter misspecification, we can evaluate predictions with a new effect-site equilibrium constant.

```
# evaluate with different ke0 parameter
<- predict(mod3ecpt, inf = inf_3ecpt, tms = tms_pred,
pre_misspec pars_pk = c(ke0 = 0.8))
<- data.frame(pre_misspec, variable = "effect (ke0=0.8)", time = tms_pred)
dat_misspec + geom_line(data = dat_misspec, aes(x = time, y = c4, color = variable)) p
```

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
<- predict(mod1cpt, inf = inf_3ecpt, tms = tms_pred)
pre_1cpt <- data.frame(pre_1cpt, variable = "plasma (1 cmpt)", time = tms_pred)
dat_1cpt # simulated observations
<- simulate(mod1cpt, seed = 1, inf = inf_3ecpt, tms = tms_obs, sigma_mult = 0.2)
obs_1cpt
+ geom_line(data = dat_1cpt, aes(x = time, y = c1, color = variable)) +
p 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, `emx`

and `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.

```
<- update(mod3ecpt, pdfn = emax, pdinv = emax_inv,
modpd pars_pd = c(e0 = 100, emx = 100, c50 = 3.5, gamma = 2.2))
```

PD targets are passed along with the updated `pkmod`

to
`inf_tci`

, which will assume values are PD values (unless
overridden by the `ignore_pd`

argument of
`inf_tci`

).

`<- inf_tci(target_vals = c(70,60,50,50), target_tms = c(0,2,3,10), pkmod = modpd, type = "effect") inf_pd `

We can then similarly use `predict.pkmod`

and
`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
<- predict(modpd, inf = inf_pd, tms = tms_pred)
pre_pd # pd observations: 10 sec = 1/6 min
<- seq(1/6,10,1/6)
tms_pd_obs # simulate responses with additive error and parameter misspecification
<- simulate(modpd, seed = 1, inf = inf_pd, tms = tms_pd_obs, sigma_add = 5,
obs_pd pars_pk = c(ke0 = 0.7), pars_pd = c(c50 = 3, gamma = 1.8))
# plot results
<- data.frame(time = tms_pred, `plasma (3 cmpt)` = pre_pd[,"c1"],
dat_pd `effect (ke0=1.2)` = pre_pd[,"c4"],
BIS = pre_pd[,"pdresp"],
check.names = FALSE)
<- melt(dat_pd, id = "time")
dat_pdm $type <- as.factor(ifelse(dat_pdm$variable == "BIS", "PD","PK"))
dat_pdm<- data.frame(time = tms_pd_obs, BIS = obs_pd,
dat_pd_obs 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 `pkmod`

and `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 `delay`

.

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.

```
<- update(mod3ecpt, pars_pk = c(cl = 20, q2 = 1.5, ke0 = 1.8))
mod_true <- simulate_tci(pkmod_prior = mod3ecpt,
sim_ol 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.

```
<- update(mod3ecpt, sigma_mult = 0.2,
mod3ecpt Omega = matrix(diag(c(1.2,0.6,1.5,0.05)), 4,4,
dimnames = list(NULL, c("cl","q2","v","ke0"))))
<- simulate_tci(pkmod_prior = mod3ecpt,
sim_cl 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)")
```

Abuhelwa, Ahmad Y., David J R Foster, and Richard N. Upton. 2015.
“ADVAN-style analytical solutions for common
pharmacokinetic models.” *Journal of Pharmacological
and Toxicological Methods* 73: 42–48. https://doi.org/10.1016/j.vascn.2015.03.004.

Jacobs, James R. 1990. “Algorithm for Optimal
Linear Model-Based Control with Application to Pharmacokinetic
Model-Driven Drug Delivery.” *IEEE Transactions on
Biomedical Engineering* 37 (1): 107–9. https://doi.org/10.1109/10.43622.

Shafer, Steven L., and Keith M. Gregg. 1992. “Algorithms to rapidly achieve and maintain stable drug
concentrations at the site of drug effect with a computer-controlled
infusion pump.” *Journal of Pharmacokinetics and
Biopharmaceutics* 20 (2): 147–69. https://doi.org/10.1007/BF01070999.