Skip to contents

Parameter sweeping systematically explores a grid of candidate parameter values and identifies the combination that best reproduces the observed data. In nlmixr2autoinit, this approach is used to derive initial estimates for model structures that cannot be solved analytically from one-compartment PK analysis, specifically:

  • Multi-compartment models (2-compartment and 3-compartment)
  • Nonlinear elimination models (Michaelis-Menten kinetics: VmaxV_{max} and KmK_m)

This vignette explains how to run parameter sweeping independently and how to visualise the sweep results to understand the parameter space.

Running Parameter Sweeping

We can bypass the automated base parameter calculation from nlmixr2autoinit and instead specify CLCL, VdV_d, and KaK_a directly. This is useful when we are only interested in exploring the complex model parameters such as VmaxV_{max} and KmK_m, without re-running the full estimation pipeline.

For example, sim_sens_1cmpt_mm() then performs the sweep for VmaxV_{max} and KmK_m in a one-compartment model with Michaelis-Menten elimination for IV data.

library(nlmixr2autoinit)
#> Loading required package: nlmixr2data
# Example: IV dosing scenario with automatic Vmax and Km
out1 <- sim_sens_1cmpt_mm(
  dat    = Bolus_1CPTMM,
  sim_vmax = list(mode = "manual", values = c(200, 400, 600, 800, 1000, 1200)),
  sim_km   = list(mode = "manual", values = c(5, 10, 50, 100, 200, 400)),
  sim_vd   = list(mode = "manual", values = 70),
  sim_ka   = list(mode = "manual", values = NA),
  route = "iv"
)

We can inspect the first few rows of the simulated output and see that multiple goodness-of-fit metrics are returned for each parameter combination. Definitions of these metrics are provided in ?metrics..

Here we focus on rRMSE2, the symmetric relative root mean squared error, which normalizes each predicted–observed difference by the mean of the predicted-observed pair and is used here to compare parameter sets.

head(out1)
#> # A tibble: 6 × 11
#>    Vmax    Km    Vd    Ka      APE   MAE  MAPE  RMSE rRMSE1 rRMSE2
#>   <dbl> <dbl> <dbl> <dbl>    <dbl> <dbl> <dbl> <dbl>  <dbl>  <dbl>
#> 1   200     5    70    NA 5572697.  802.  608. 1400.   87.9   73.9
#> 2   200    10    70    NA 5594331.  805.  621. 1402.   88.0   73.6
#> 3   200    50    70    NA 5753722.  828.  729. 1413.   88.8   76.8
#> 4   200   100    70    NA 5899008.  849.  821. 1425.   89.5   79.3
#> 5   200   200    70    NA 6096946.  877.  933. 1443.   90.6   81.9
#> 6   200   400    70    NA 6344339   913. 1059. 1468.   92.2   84.4
#> # ℹ 1 more variable: Cumulative.Time.Sec <dbl>
# # A tibble: 6 × 11
#    Vmax    Km    Vd    Ka      APE   MAE  MAPE  RMSE rRMSE1 rRMSE2 Cumulative.Time.Sec
#   <dbl> <dbl> <dbl> <dbl>    <dbl> <dbl> <dbl> <dbl>  <dbl>  <dbl>               <dbl>
# 1   200     5    70    NA 5572697.  802.  608. 1400.   87.9   73.9                0.33
# 2   200    10    70    NA 5594331.  805.  621. 1402.   88.0   73.6                0.64
# 3   200    50    70    NA 5753722.  828.  729. 1413.   88.8   76.8                0.99
# 4   200   100    70    NA 5899008.  849.  821. 1425.   89.5   79.3                1.29
# 5   200   200    70    NA 6096946.  877.  933. 1443.   90.6   81.9                1.6 
# 6   200   400    70    NA 6344339   913. 1059. 1468.   92.2   84.4                1.92

Visualising Sweep Results

We can visualize the sweep results as a heatmap to examine how model fit changes across combinations of VmaxV_{max} and KmK_m, with lower rRMSE2 indicating better fit. In this example, the optimal parameter combination is identified at Vmax=1000mg/hV_{max}=1000 mg/h and Km=200mg/LK_m=200 mg/L, which gives the lowest rRMSE2 among all tested values.

library(ggplot2)
best <- out1[which.min(out1$rRMSE2), ]
p1<-ggplot(out1, aes(x = factor(Vmax), y = factor(Km), fill = rRMSE2)) +
  geom_tile() +
  geom_tile(
    data = best,
    aes(x = factor(Vmax), y = factor(Km)),
    fill  = NA,
    color = "black",
    linewidth = 1.2
  ) +
  geom_text(aes(label = round(rRMSE2, 1)), size = 4.5, color = "grey20") +
  scale_fill_gradient(
    low  = "#1a9641",
    high = "#f7f7f7",
    name = "rRMSE2"
  ) +
  labs(
    x = expression(V[max]),
    y = expression(K[m])
  ) +
  theme_bw() +
  theme(
    axis.title = element_text(size = 11, face = "bold"),
    axis.text  = element_text(size = 10)
  )
print(p1)

See Also