Warn on flip-flop kinetics? about pmxtools

[billdenney]

Nevermind about changing the ratio, we can just use ka < kel based on this quick simulation. When ka < kel based on flip-flip, the half-life is only about 20% higher than expected:

library(tidyverse)
#> Warning: package 'ggplot2' was built under R version 4.2.2
#> Warning: package 'dplyr' was built under R version 4.2.2
#> Warning: package 'stringr' was built under R version 4.2.2
library(nlmixr2)
#> Warning: package 'nlmixr2' was built under R version 4.2.2
#> Loading required package: nlmixr2data

## The basic model consiss of an ini block that has initial estimates
one.compartment <- function() {
  ini({
    tka <- 0.45 # Log Ka
    tcl <- 1 # Log Cl
    tv <- 3.45    # Log V
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    add.sd <- 0.7
  })
  # and a model block with the error sppecification and model specification
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    d/dt(depot) = -ka * depot
    d/dt(center) = ka * depot - cl / v * center
    cp = center / v
    cp ~ add(add.sd)
  })
}

## The fit is performed by the function nlmixr/nlmix2 specifying the model, data and estimate
fit <- nlmixr2(one.compartment, theo_sd,  est="saem", saemControl(print=0))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of saem model...
#> ✔ done
#> → finding duplicate expressions in saem model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in saem model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> ✔ done
#> rxode2 2.0.11.9000 using 8 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
#> Calculating covariance matrix
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of saem model...
#> ✔ done
#> → finding duplicate expressions in saem predOnly model 0...
#> → finding duplicate expressions in saem predOnly model 1...
#> → optimizing duplicate expressions in saem predOnly model 1...
#> → finding duplicate expressions in saem predOnly model 2...
#> ✔ done
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 5952
#> → compress phiM in nlmixr2 object, save 62360
#> → compress parHist in nlmixr2 object, save 9592
#> → compress saem0 in nlmixr2 object, save 28752
fit
#> ── nlmixr² SAEM OBJF by FOCEi approximation ──
#> 
#>  Gaussian/Laplacian Likelihoods: AIC(fit) or fit$objf etc. 
#>  FOCEi CWRES & Likelihoods: addCwres(fit) 
#> 
#> ── Time (sec fit$time): ──
#> 
#>         setup covariance saem table compress other
#> elapsed 0.002       0.02 5.32  0.05     0.03 2.718
#> 
#> ── Population Parameters (fit$parFixed or fit$parFixedDf): ──
#> 
#>        Parameter  Est.     SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
#> tka       Log Ka 0.454  0.196 43.1       1.57 (1.07, 2.31)     71.5   -0.0203% 
#> tcl       Log Cl  1.02 0.0853  8.4       2.76 (2.34, 3.26)     27.6      3.46% 
#> tv         Log V  3.45 0.0454 1.32       31.5 (28.8, 34.4)     13.4      9.89% 
#> add.sd           0.693                               0.693                     
#>  
#>   Covariance Type (fit$covMethod): linFim
#>   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 
#>   Censoring (fit$censInformation): No censoring
#> 
#> ── Fit Data (object fit is a modified tibble): ──
#> # A tibble: 132 × 19
#>   ID     TIME    DV  PRED    RES IPRED   IRES  IWRES eta.ka eta.cl   eta.v    cp
#>   <fct> <dbl> <dbl> <dbl>  <dbl> <dbl>  <dbl>  <dbl>  <dbl>  <dbl>   <dbl> <dbl>
#> 1 1      0     0.74  0     0.74   0     0.74   1.07   0.103 -0.491 -0.0820  0   
#> 2 1      0.25  2.84  3.27 -0.426  3.87 -1.03  -1.48   0.103 -0.491 -0.0820  3.87
#> 3 1      0.57  6.57  5.85  0.723  6.82 -0.246 -0.356  0.103 -0.491 -0.0820  6.82
#> # … with 129 more rows, and 7 more variables: depot <dbl>, center <dbl>,
#> #   ka <dbl>, cl <dbl>, v <dbl>, tad <dbl>, dosenum <dbl>

min(theo_sd$DV[theo_sd$DV > 0])
#> [1] 0.15
lloq <- 0.15
modcoef <- exp(coef(fit)$fixed)

time_step <- 0.1
times <- seq(0, 240, by= time_step)

fit_hl_simple <- function(conc) {
  tlast <- PKNCA::pk.calc.tlast(conc = conc, time = seq_along(conc))
  stopifnot(tlast < length(conc))
  fit <- stats::.lm.fit(x = cbind(1, seq_len(3)*time_step), y = log(conc[(tlast - 2):tlast]))
  log(2)/-fit$coefficients[2]
}

d_plot <-
  tibble(
    ka = bsd.report::seq_exp(0.01, 1, length.out = 1000),
    kel = modcoef["tcl"]/modcoef["tv"]
  ) %>%
  mutate(
    conc =
      pmap(
        .l = list(ka = ka),
        .f = pmxTools::calc_sd_1cmt_linear_oral_1,
        t = times,
        dose = 300,
        #ka = modcoef["tka"],
        CL = modcoef["tcl"],
        V = modcoef["tv"]
      ),
    conc_loq =
      pmap(
        .l = list(conc),
        .f = ~ifelse(.x > lloq, .x, 0)
      ),
    half_life=
      pmap_dbl(
        .l = list(conc = conc_loq),
        .f =  fit_hl_simple
      )
  )

half_life_model <-
  pmxTools::calc_derived_1cpt(
    ka = modcoef["tka"],
    CL = modcoef["tcl"],
    V = modcoef["tv"]
  )$thalf

ggplot(d_plot, aes(x=ka/kel, y=half_life/half_life_model)) +
  geom_line() +
  geom_hline(yintercept = 1.2, linetype = "63") +
  geom_vline(xintercept = 1, linetype = "63") +
  scale_x_log10() + scale_y_log10()

^{Created on 2023-01-12 with reprex v2.0.2}

[kestrel99]

This will go into 1.2.5

[billdenney]

Thanks!