Initial estimates for a compartmental population PK model can be obtained using babelmixr2 with the "pknca" estimation method. Also, the central compartment scaling factor can be auto-generated based on units for dosing, concentration measurement, desired volume of distribution units, and time.
You do not need to perform NCA analysis by hand; the "pknca" estimation method will perform NCA analysis using the PKNCA package automatically.
The methods used for converting NCA calculations to parameter estimates are described in the help for nlmixr2Est.pknca().
Initial model setup is the same as for any other nlmixr2 model. You must load the babelmixr2 library so that the nlmixr() function recognizes est = "pknca".
library(babelmixr2)
one.compartment <- function() {
ini({
tka <- log(1.57); label("Ka (1/hr)")
tcl <- log(2.72); label("Cl (L/hr)")
tv <- log(31.5); label("V (L)")
eta.ka ~ 0.6
eta.cl ~ 0.3
eta.v ~ 0.1
add.sd <- 0.7; label("additive residual error (mg/L)")
})
# and a model block with the error specification and model specification
model({
ka <- exp(tka + eta.ka)
cl <- exp(tcl + eta.cl)
vc <- exp(tv + eta.v)
d/dt(depot) <- -ka * depot
d/dt(center) <- ka * depot - cl / vc * center
cp <- center / vc
cp ~ add(add.sd)
})
}To use PKNCA to get initial estimates, use est = "pknca" instead of one of the other nlmixr2 estimation methods.
For unit conversions, provide the units to the control = pkncaControl() argument. Unit conversions are only supported when the units can be automatically converted; mass/volume can be converted to any other mass/volume ratio, but mass to molar or molar to mass cannot because there is not a single mass-to-molar conversion factor.
prepared <-
nlmixr2(
one.compartment,
data = theo_sd,
est = "pknca",
control = pkncaControl(concu = "ng/mL", doseu = "mg", timeu = "hr", volumeu = "L")
)
#> Loading required namespace: testthat
#> ℹ change initial estimate (0.89314878960486) and upper/lower bound (-3.50655789731998 to 3.72508541597241) of `tka`
#> → significant model change detected
#> → removed from model: '$getSplitModel'
#> ℹ change initial estimate (8.41044546236311) and upper/lower bound (5.51439905878865 to 10.899462850803) of `tcl`
#> ℹ change initial estimate (10.5377244826318) and upper/lower bound (7.94567233496473 to 13.1050053785005) of `tv`Now, you have the prepared model with updated initial estimates and the NCA results embedded. You can see the new model and the PKNCA estimates by looking at the prepared$ui (the model as interpreted by rxode2) and prepared$nca (the PKNCAresults object).
Note that in the new model, the fixed effect initial estimates are changed from their original values. The residual error and between-subject variability are unchanged.
prepared$ui
#> ── rxode2-based free-form 2-cmt ODE model ──────────────────────────────────────
#> ── Initalization: ──
#> Fixed Effects ($theta):
#> tka tcl tv add.sd
#> 0.8931488 8.4104455 10.5377245 0.7000000
#>
#> Omega ($omega):
#> eta.ka eta.cl eta.v
#> eta.ka 0.6 0.0 0.0
#> eta.cl 0.0 0.3 0.0
#> eta.v 0.0 0.0 0.1
#>
#> States ($state or $stateDf):
#> Compartment Number Compartment Name
#> 1 1 depot
#> 2 2 center
#> ── μ-referencing ($muRefTable): ──
#> theta eta level
#> 1 tka eta.ka id
#> 2 tcl eta.cl id
#> 3 tv eta.v id
#>
#> ── Model (Normalized Syntax): ──
#> function() {
#> ini({
#> tka <- c(-3.50655789731998, 0.89314878960486, 3.72508541597241)
#> label("Ka (1/hr)")
#> tcl <- c(5.51439905878865, 8.41044546236311, 10.899462850803)
#> label("Cl (L/hr)")
#> tv <- c(7.94567233496473, 10.5377244826318, 13.1050053785005)
#> label("V (L)")
#> add.sd <- c(0, 0.7)
#> label("additive residual error (mg/L)")
#> eta.ka ~ 0.6
#> eta.cl ~ 0.3
#> eta.v ~ 0.1
#> })
#> model({
#> ka <- exp(tka + eta.ka)
#> cl <- exp(tcl + eta.cl)
#> vc <- exp(tv + eta.v)
#> d/dt(depot) <- -ka * depot
#> d/dt(center) <- ka * depot - cl/vc * center
#> cp <- 1000 * center/vc
#> cp ~ add(add.sd)
#> })
#> }
knitr::knit_print(
summary(prepared$nca)
)
#> Interval Start Interval End N AUClast (hr*ng/mL) Cmax (ng/mL)
#> 0 24 12 74.6 [24.3] .
#> 0 Inf 12 . 8.65 [17.0]
#> Tmax (hr) CL (based on AUClast) (mg/(hr*ng/mL))
#> . 4.22 [23.0]
#> 1.14 [0.630, 3.55] .
#> Vss (based on AUClast) (mg/(ng/mL)) Half-life (hr) AUCinf,obs (hr*ng/mL)
#> 25.0 [18.5] . .
#> . 8.18 [2.12] 115 [28.4]
#> Cmax (dose-normalized) ((ng/mL)/mg)
#> .
#> 0.0274 [18.1]
#>
#> Caption: AUClast, Cmax, CL (based on AUClast), Vss (based on AUClast), AUCinf,obs, Cmax (dose-normalized): geometric mean and geometric coefficient of variation; Tmax: median and range; Half-life: arithmetic mean and standard deviationFrom the updated model, you can perform estimation on the new model object, as with any other model that has been created for nlmixr2:
fit <- nlmixr(prepared, data = theo_sd, est = "focei", control = list(print = 0))
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of full model...
#> ✔ done
#> → calculate jacobian
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → calculate sensitivities
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → calculate ∂(f)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → calculate ∂(R²)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in inner model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling inner model...
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> ✔ done
#> → finding duplicate expressions in FD model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in FD model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling EBE model...
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> ✔ done
#> → compiling events FD model...
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> ✔ done
#> rxode2 2.1.1 using 8 threads (see ?getRxThreads)
#> no cache: create with `rxCreateCache()`
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> done
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 5952
#> → compress parHistData in nlmixr2 object, save 5672
fit
#> ── nlmixr² FOCEi (outer: nlminb) ──
#>
#> OBJF AIC BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> FOCEi 117.2304 373.8302 394.0098 -179.9151 170.5271 9.700921
#>
#> ── Time (sec fit$time): ──
#>
#> setup optimize covariance table compress other
#> elapsed 0.020023 0.370917 0.370918 0.025 0.008 3.634142
#>
#> ── Population Parameters (fit$parFixed or fit$parFixedDf): ──
#>
#> Parameter Est. SE %RSE
#> tka Ka (1/hr) 0.473 0.22 46.6
#> tcl Cl (L/hr) 7.91 0.0655 0.827
#> tv V (L) 10.4 0.0307 0.296
#> add.sd additive residual error (mg/L) 0.694
#> Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
#> tka 1.6 (1.04, 2.47) 69.7 0.805%
#> tcl 2.73e+03 (2.4e+03, 3.11e+03) 27.2 7.47%
#> tv 3.19e+04 (3.01e+04, 3.39e+04) 16.6 18.5%
#> add.sd 0.694
#>
#> Covariance Type (fit$covMethod): r,s
#> No correlations in between subject variability (BSV) matrix
#> Full BSV covariance (fit$omega) or correlation (fit$omegaR; diagonals=SDs)
#> Distribution stats (mean/skewness/kurtosis/p-value) available in fit$shrink
#> Information about run found (fit$runInfo):
#> • gradient problems with initial estimate and covariance; see $scaleInfo
#> • last objective function was not at minimum, possible problems in optimization
#> • ETAs were reset to zero during optimization; (Can control by foceiControl(resetEtaP=.))
#> • initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=))
#> Censoring (fit$censInformation): No censoring
#> Minimization message (fit$message):
#> false convergence (8)
#> In an ODE system, false convergence may mean "useless" evaluations were performed.
#> See https://tinyurl.com/yyrrwkce
#> It could also mean the convergence is poor, check results before accepting fit
#> You may also try a good derivative free optimization:
#> nlmixr2(...,control=list(outerOpt="bobyqa"))
#>
#> ── Fit Data (object fit is a modified tibble): ──
#> # A tibble: 132 × 22
#> ID TIME DV PRED RES WRES IPRED IRES IWRES CPRED CRES CWRES
#> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0 0.74 0 0.74 1.07 0 0.74 1.07 0 0.74 1.07
#> 2 1 0.25 2.84 3.27 -0.432 -0.229 3.84 -1.00 -1.45 3.23 -0.389 -0.181
#> 3 1 0.57 6.57 5.84 0.733 0.291 6.79 -0.217 -0.312 5.78 0.791 0.281
#> # ℹ 129 more rows
#> # ℹ 10 more variables: eta.ka <dbl>, eta.cl <dbl>, eta.v <dbl>, depot <dbl>,
#> # center <dbl>, ka <dbl>, cl <dbl>, vc <dbl>, tad <dbl>, dosenum <dbl>To get the initial estimate, babelmixr2 automatically converts your modeling dataset to the format the is needed for PKNCA, and the NCA is automatically performed using all the data.
In some cases (e.g. studies with sparse data), NCA may not be feasible. In those cases, you can provide a different dataset to PKNCA compared to the full modeling dataset. Usually, the simplest method is to provide single-dose, dense-sampling, dose-ranging data (i.e. the single-ascending dose portion of the first-in-human study) to be estimated.
To do this, give your data to PKNCA using the ncaData argument to pkncaControl() as follows:
# Choose a subset of the full dataset for NCA
dNCA <- theo_sd[theo_sd$ID <= 6, ]
preparedNcaData <-
nlmixr2(
one.compartment,
data = theo_sd,
est = "pknca",
control = pkncaControl(concu = "ng/mL", doseu = "mg", timeu = "hr", volumeu = "L", ncaData = dNCA)
)
#> ℹ change initial estimate (0.929027077269762) and upper/lower bound (-3.50655789731998 to 3.32136703319919) of `tka`
#> → significant model change detected
#> → removed from model: '$getSplitModel'
#> ℹ change initial estimate (8.3955404628088) and upper/lower bound (5.85241523541802 to 10.7637056987378) of `tcl`
#> ℹ change initial estimate (10.5377244826318) and upper/lower bound (7.94370069836702 to 13.1024358787022) of `tv`
preparedNcaData$ui
#> ── rxode2-based free-form 2-cmt ODE model ──────────────────────────────────────
#> ── Initalization: ──
#> Fixed Effects ($theta):
#> tka tcl tv add.sd
#> 0.9290271 8.3955405 10.5377245 0.7000000
#>
#> Omega ($omega):
#> eta.ka eta.cl eta.v
#> eta.ka 0.6 0.0 0.0
#> eta.cl 0.0 0.3 0.0
#> eta.v 0.0 0.0 0.1
#>
#> States ($state or $stateDf):
#> Compartment Number Compartment Name
#> 1 1 depot
#> 2 2 center
#> ── μ-referencing ($muRefTable): ──
#> theta eta level
#> 1 tka eta.ka id
#> 2 tcl eta.cl id
#> 3 tv eta.v id
#>
#> ── Model (Normalized Syntax): ──
#> function() {
#> ini({
#> tka <- c(-3.50655789731998, 0.929027077269762, 3.32136703319919)
#> label("Ka (1/hr)")
#> tcl <- c(5.85241523541802, 8.3955404628088, 10.7637056987378)
#> label("Cl (L/hr)")
#> tv <- c(7.94370069836702, 10.5377244826318, 13.1024358787022)
#> label("V (L)")
#> add.sd <- c(0, 0.7)
#> label("additive residual error (mg/L)")
#> eta.ka ~ 0.6
#> eta.cl ~ 0.3
#> eta.v ~ 0.1
#> })
#> model({
#> ka <- exp(tka + eta.ka)
#> cl <- exp(tcl + eta.cl)
#> vc <- exp(tv + eta.v)
#> d/dt(depot) <- -ka * depot
#> d/dt(center) <- ka * depot - cl/vc * center
#> cp <- 1000 * center/vc
#> cp ~ add(add.sd)
#> })
#> }The initial estimates are now based on NCA calculated from the dNCA dataset rather than the full theo_sd dataset.
If you already have NCA results calculated by PKNCA with the required parameters (“tmax”, “cmax.dn”, and “cllast”), you can provide those instead using the pkncaControl(ncaResults) argument.
To update the initial estimates, the model must have parameters in the model() block with the names that are expected by est = "pknca". The expected names are:
kavcclvpvp2qq2Any of those parameter names that are found in the model() block will be automatically traced back to the initial conditions (ini() block), and the parameter values will be updated. If the parameter is estimated on the log scale, the updated parameter value will automatically be converted to the log scale.