pharmacokinetics

name	pharmacokinetics
description	Pharmacokinetic and pharmacodynamic analysis in R, including NCA, compartmental modeling, and bioequivalence.

Pharmacokinetics and Pharmacodynamics in R

Overview

Comprehensive pharmacokinetic (PK) and pharmacodynamic (PD) modeling in R covering non-compartmental analysis (NCA), compartmental PK modeling, population PK with nonlinear mixed effects, bioequivalence assessment, PK/PD modeling, and drug-drug interaction evaluation.

Non-Compartmental Analysis (NCA)

Using PKNCA Package

library(PKNCA)

# Prepare concentration data
conc_data <- data.frame(
  subject = rep(1:3, each = 8),
  time = rep(c(0, 0.5, 1, 2, 4, 8, 12, 24), 3),
  concentration = c(
    0, 450, 380, 290, 180, 85, 40, 12,
    0, 520, 410, 320, 200, 95, 45, 15,
    0, 480, 395, 305, 190, 90, 42, 13
  )
)

# Prepare dose data
dose_data <- data.frame(
  subject = 1:3,
  time = 0,
  dose = 100  # mg
)

# Create PKNCA objects
conc_obj <- PKNCAconc(concentration ~ time | subject, data = conc_data)
dose_obj <- PKNCAdose(dose ~ time | subject, data = dose_data)

# Combine into analysis object
data_obj <- PKNCAdata(conc_obj, dose_obj)

# Run NCA
results <- pk.nca(data_obj)

# View results
summary(results)
as.data.frame(results)

Custom NCA Parameters

library(PKNCA)

# Specify intervals and parameters
intervals <- data.frame(
  start = 0,
  end = 24,
  cmax = TRUE,
  tmax = TRUE,
  auclast = TRUE,
  aucinf.obs = TRUE,
  half.life = TRUE,
  cl.obs = TRUE,
  vss.obs = TRUE,
  mrt.last = TRUE
)

# Create analysis with custom intervals
data_obj <- PKNCAdata(
  conc_obj,
  dose_obj,
  intervals = intervals
)

results <- pk.nca(data_obj)

# Extract specific parameters
pk_summary <- summary(results)

# Individual subject parameters
individual_params <- as.data.frame(results) |>
  tidyr::pivot_wider(
    id_cols = subject,
    names_from = PPTESTCD,
    values_from = PPORRES
  )

Manual NCA Calculations

# Basic NCA calculations for a single subject
calculate_nca <- function(time, conc, dose) {

  # Cmax and Tmax
  cmax <- max(conc)
  tmax <- time[which.max(conc)]

  # AUC by linear trapezoidal rule
  auc_last <- sum(diff(time) * (head(conc, -1) + tail(conc, -1)) / 2)

  # Terminal phase (last 3 points for lambda_z)
  n_terminal <- 3
  terminal_idx <- (length(time) - n_terminal + 1):length(time)
  terminal_time <- time[terminal_idx]
  terminal_conc <- conc[terminal_idx]

  # Lambda_z from log-linear regression
  fit <- lm(log(terminal_conc) ~ terminal_time)
  lambda_z <- -coef(fit)[2]

  # Half-life
  half_life <- log(2) / lambda_z

  # AUC extrapolated to infinity
  auc_extrap <- tail(conc, 1) / lambda_z
  auc_inf <- auc_last + auc_extrap

  # Clearance (CL/F for oral)
  cl_f <- dose / auc_inf

  # Volume of distribution (Vd/F)
  vd_f <- cl_f / lambda_z

  # Mean residence time
  # AUMC (first moment)
  aumc <- sum(diff(time) * (head(time * conc, -1) + tail(time * conc, -1)) / 2)
  aumc_extrap <- tail(conc, 1) / lambda_z^2 + tail(time * conc, 1) / lambda_z
  aumc_inf <- aumc + aumc_extrap
  mrt <- aumc_inf / auc_inf

  return(data.frame(
    Cmax = cmax,
    Tmax = tmax,
    AUClast = auc_last,
    AUCinf = auc_inf,
    Lambda_z = lambda_z,
    Half_life = half_life,
    CL_F = cl_f,
    Vd_F = vd_f,
    MRT = mrt
  ))
}

Concentration-Time Visualization

Spaghetti Plots

library(ggplot2)

# Individual profiles (spaghetti plot)
ggplot(conc_data, aes(x = time, y = concentration, group = subject)) +
  geom_line(alpha = 0.5) +
  geom_point(alpha = 0.5) +
  stat_summary(aes(group = 1), fun = mean, geom = "line",
               color = "red", size = 1.5) +
  stat_summary(aes(group = 1), fun = mean, geom = "point",
               color = "red", size = 2) +
  labs(x = "Time (hours)", y = "Concentration (ng/mL)",
       title = "Concentration-Time Profile") +
  theme_bw()

# Semi-log plot
ggplot(conc_data, aes(x = time, y = concentration, group = subject)) +
  geom_line(alpha = 0.5) +
  geom_point(alpha = 0.5) +
  scale_y_log10() +
  labs(x = "Time (hours)", y = "Concentration (ng/mL) - Log Scale",
       title = "Semi-Logarithmic Concentration-Time Profile") +
  theme_bw()

Summary Plots with Error Bars

library(ggplot2)
library(dplyr)

# Calculate summary statistics
conc_summary <- conc_data |>
  group_by(time) |>
  summarise(
    mean_conc = mean(concentration),
    sd_conc = sd(concentration),
    geom_mean = exp(mean(log(concentration + 0.001))),
    n = n(),
    se_conc = sd_conc / sqrt(n)
  )

# Mean + SD plot
ggplot(conc_summary, aes(x = time, y = mean_conc)) +
  geom_line(size = 1) +
  geom_point(size = 2) +
  geom_errorbar(aes(ymin = mean_conc - sd_conc,
                    ymax = mean_conc + sd_conc),
                width = 0.5) +
  labs(x = "Time (hours)",
       y = "Concentration (ng/mL)",
       title = "Mean (± SD) Concentration-Time Profile") +
  theme_bw()

# Geometric mean plot
ggplot(conc_summary, aes(x = time, y = geom_mean)) +
  geom_line(size = 1) +
  geom_point(size = 2) +
  scale_y_log10() +
  labs(x = "Time (hours)",
       y = "Geometric Mean Concentration (ng/mL)",
       title = "Geometric Mean Concentration Profile") +
  theme_bw()

Compartmental PK with mrgsolve

One-Compartment Model

library(mrgsolve)

# Define one-compartment model with first-order absorption
code_1cmt <- '
$PARAM CL = 10, V = 100, KA = 1.5

$CMT GUT CENT

$ODE
dxdt_GUT = -KA * GUT;
dxdt_CENT = KA * GUT - (CL/V) * CENT;

$TABLE
double CP = CENT / V;

$CAPTURE CP
'

# Compile model
mod_1cmt <- mcode("pk1cmt", code_1cmt)

# Single dose simulation
out <- mod_1cmt |>
  ev(amt = 100, cmt = 1) |>  # 100 mg oral dose
  mrgsim(end = 24, delta = 0.1)

# Plot
plot(out, CP ~ time)

Two-Compartment Model

library(mrgsolve)

# Two-compartment model with central and peripheral
code_2cmt <- '
$PARAM CL = 10, V1 = 50, V2 = 100, Q = 15, KA = 1.2

$CMT GUT CENT PERIPH

$ODE
double K10 = CL / V1;
double K12 = Q / V1;
double K21 = Q / V2;

dxdt_GUT = -KA * GUT;
dxdt_CENT = KA * GUT - K10 * CENT - K12 * CENT + K21 * PERIPH;
dxdt_PERIPH = K12 * CENT - K21 * PERIPH;

$TABLE
double CP = CENT / V1;

$CAPTURE CP
'

mod_2cmt <- mcode("pk2cmt", code_2cmt)

# Simulate IV bolus
out_iv <- mod_2cmt |>
  ev(amt = 100, cmt = 2) |>  # IV bolus to central
  mrgsim(end = 48, delta = 0.1)

plot(out_iv, CP ~ time, log = TRUE)

Multiple Dosing

library(mrgsolve)

# Multiple dose regimen
dosing_regimen <- ev(
  amt = 100,
  ii = 12,      # Dosing interval (hours)
  addl = 6,     # Additional doses
  cmt = 1       # Oral
)

out_multi <- mod_1cmt |>
  ev(dosing_regimen) |>
  mrgsim(end = 96, delta = 0.1)

plot(out_multi, CP ~ time)

# Steady-state check
out_ss <- mod_1cmt |>
  ev(amt = 100, ii = 12, addl = 100, cmt = 1, ss = 1) |>
  mrgsim(end = 24, delta = 0.1)

Population PK with nlmixr2

One-Compartment PopPK Model

library(nlmixr2)

# Define population PK model
one_cmt_pop <- function() {
  ini({
    # Fixed effects (thetas)
    tka <- log(1)       # Log Ka
    tcl <- log(10)      # Log CL
    tv <- log(100)      # Log V

    # Random effects (omegas)
    eta.ka ~ 0.6        # IIV on Ka
    eta.cl ~ 0.3        # IIV on CL
    eta.v ~ 0.1         # IIV on V

    # Residual error
    add.err <- 0.1      # Additive
    prop.err <- 0.1     # Proportional
  })
  model({
    # Individual parameters
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)

    # Differential equations
    d/dt(depot) <- -ka * depot
    d/dt(central) <- ka * depot - cl/v * central

    # Concentration
    cp <- central / v

    # Combined error model
    cp ~ add(add.err) + prop(prop.err)
  })
}

# Prepare data (nlmixr2 format)
pk_data <- data.frame(
  ID = rep(1:10, each = 8),
  TIME = rep(c(0, 0.5, 1, 2, 4, 8, 12, 24), 10),
  DV = rnorm(80, 100, 20),  # Observed concentrations
  AMT = c(rep(c(100, rep(0, 7)), 10)),
  EVID = c(rep(c(1, rep(0, 7)), 10)),
  CMT = 1
)

# Fit model using SAEM
fit <- nlmixr2(one_cmt_pop, pk_data, est = "saem",
               control = saemControl(print = 50))

# Summary
summary(fit)

# Parameter estimates
fit$parFixed      # Fixed effects
fit$omega         # Random effects variance
fit$sigma         # Residual error

# Individual predictions
fit$ipred

Covariate Model

library(nlmixr2)

# Model with covariates (weight, age)
cov_model <- function() {
  ini({
    tka <- log(1)
    tcl <- log(10)
    tv <- log(100)

    # Covariate effects
    cl_wt <- 0.75       # Allometric on CL
    v_wt <- 1           # Allometric on V

    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1

    add.err <- 0.1
    prop.err <- 0.1
  })
  model({
    # Covariate-adjusted parameters
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl) * (WT/70)^cl_wt
    v <- exp(tv + eta.v) * (WT/70)^v_wt

    d/dt(depot) <- -ka * depot
    d/dt(central) <- ka * depot - cl/v * central

    cp <- central / v
    cp ~ add(add.err) + prop(prop.err)
  })
}

# Fit with covariates
fit_cov <- nlmixr2(cov_model, pk_data_with_cov, est = "saem")

Model Diagnostics

library(nlmixr2)
library(xpose.nlmixr2)

# Create xpose database
xpdb <- xpose_data_nlmixr2(fit)

# Goodness-of-fit plots
dv_vs_pred(xpdb)           # DV vs PRED
dv_vs_ipred(xpdb)          # DV vs IPRED
res_vs_pred(xpdb)          # CWRES vs PRED
res_vs_idv(xpdb)           # CWRES vs TIME

# Individual plots
ind_plots(xpdb, nrow = 3, ncol = 3)

# VPC (Visual Predictive Check)
vpc_result <- vpc(fit, n = 500)
plot(vpc_result)

Bioequivalence Analysis

Using BE Package

library(BE)

# 2x2 crossover BE study data
be_data <- data.frame(
  subj = rep(1:24, each = 2),
  seq = rep(rep(c("TR", "RT"), each = 12), each = 2),
  prd = rep(c(1, 2), 24),