Pharmacokinetics

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 
  tidyr::pivot_wider(
    id_cols = subject,
    names_from = PPTESTCD,
    values_from = PPORRES
  )

Manual NCA Calculations

# Basic NCA calculations for a single subject
calculate_nca 
  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 
  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 
  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(dosing_regimen) |>
  mrgsim(end = 96, delta = 0.1)

plot(out_multi, CP ~ time)

# Steady-state check
out_ss 
  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 = 80 & ci_upper 
  ev(amt = 100, cmt = 1) |>
  mrgsim(end = 24, delta = 0.1)

# Plot PK and PD
library(ggplot2)
out_df 
  ev(amt = 100, cmt = 1) |>
  mrgsim(end = 72, delta = 0.1)

plot(out_indirect, RESP ~ time)

Drug-Drug Interaction

Competitive Inhibition

library(mrgsolve)

# DDI model with competitive inhibition
code_ddi 
  ev(amt = 100, cmt = 1) |>
  mrgsim(end = 24, delta = 0.1)

# Victim + perpetrator
out_ddi 
  ev(amt = 100, cmt = 1, time = 0) |>
  ev(amt = 200, cmt = 3, time = 0) |>
  mrgsim(end = 24, delta = 0.1)

# Calculate AUC ratio (DDI magnitude)
auc_alone 
    group_by(parameter) |>
    summarise(
      N = n(),
      Mean = mean(value),
      SD = sd(value),
      CV_percent = SD / Mean * 100,
      Median = median(value),
      Min = min(value),
      Max = max(value),
      Geom_Mean = exp(mean(log(value))),
      Geom_CV = sqrt(exp(var(log(value))) - 1) * 100
    ) |>
    mutate(across(where(is.numeric), ~round(., 3)))
}

# Format for regulatory submission
format_pk_table 
    mutate(
      `Mean (SD)` = paste0(round(Mean, 2), " (", round(SD, 2), ")"),
      `Median [Min, Max]` = paste0(round(Median, 2), " [",
                                   round(Min, 2), ", ", round(Max, 2), "]"),
      `Geom Mean (Geom CV%)` = paste0(round(Geom_Mean, 2), " (",
                                      round(Geom_CV, 1), "%)")
    ) |>
    select(parameter, N, `Mean (SD)`, `Median [Min, Max]`,
           `Geom Mean (Geom CV%)`)
}

Key Packages Summary

| Package | Purpose | |---------|---------| | PKNCA | Non-compartmental analysis | | mrgsolve | ODE-based PK/PD simulation | | nlmixr2 | Population PK modeling (NLME) | | rxode2 | ODE solver, nlmixr2 backend | | BE | Bioequivalence analysis | | PowerTOST | BE sample size and power | | pk | Classical PK calculations | | NonCompart | NCA calculations | | xpose.nlmixr2 | Diagnostic plots for nlmixr2 | | vpc | Visual predictive checks |

Best Practices

1. NCA conventions: Use linear trapezoidal up/log down for oral dosing
2. Terminal phase: Require at least 3 points and R² > 0.9 for lambda_z
3. Population PK: Start simple (1-cmt), add complexity as needed
4. Model qualification: GOF plots, VPC, bootstrap for parameter uncertainty
5. BE studies: Log-transform PK parameters, use ANOVA for crossover
6. Regulatory: Follow FDA/EMA guidance for NCA intervals and BE limits
7. Documentation: Report software versions, methods, and all assumptions
8. Covariate selection: Use stepwise approach with clinical plausibility

Source & license

This open-source skill is cataloged on AgentStack and links to its original source — we do not rehost the code.

• Author: choxos
• Source: choxos/BiostatAgent
• License: MIT

Install and usage instructions live in the source repository linked above.