Undefined behavior in pk1 when KA=CL/V

[toddyoder]

The concentrations in the closed form pk1 model appear to be either all zeros or NaN when KA is chosen to be equal to CL/V.

library(dplyr)
library(mrgsolve)

# KA != CL/V
modlib("pk1") %>% 
  param(CL = 0.391, V = 30, KA = 0.013) %>% 
  ev(amt = 100) %>% 
  mrgsim() %>% 
  as_tibble %>% 
  pull(CENT) %>% 
  summary

# KA = CL/V
modlib("pk1") %>% 
  param(CL = 0.39, V = 30, KA = 0.013) %>% 
  ev(amt = 100) %>% 
  mrgsim() %>% 
  as_tibble %>% 
  pull(CENT) %>% 
  summary

# KA = CL/V
modlib("pk1") %>% 
  param(CL = 1, V = 20, KA = 0.05) %>% 
  ev(amt = 100) %>% 
  mrgsim() %>% 
  as_tibble %>% 
  pull(CENT) %>% 
  summary

Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.000   6.371  12.868  12.162  18.312  22.829 

Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0       0       0       0       0       0 

Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
0       0       0       0       0       0      24 

[kylebaron]

https://github.com/metrumresearchgroup/mrgsolve/blob/3d5d70189b61345478f154de30000c9288b886e8/src/odeproblem.cpp#L387

[kylebaron]

Updated behavior

library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(mrgsolve)
#> 
#> Attaching package: 'mrgsolve'
#> The following object is masked from 'package:stats':
#> 
#>     filter

# KA != CL/V
modlib("pk1") %>% 
  param(CL = 0.391, V = 30, KA = 0.013) %>% 
  ev(amt = 100) %>% 
  mrgsim() %>% 
  as_tibble %>% 
  pull(CENT) %>% 
  summary
#> Building pk1 ...
#> done.
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   0.000   6.371  12.868  12.162  18.312  22.829

# KA = CL/V
modlib("pk1") %>% 
  param(CL = 0.39, V = 30, KA = 0.013) %>% 
  ev(amt = 100) %>% 
  mrgsim() %>% 
  as_tibble %>% 
  pull(CENT) %>% 
  summary
#> Loading model from cache.
#> Error in eval(expr, envir, enclos): k10 is too close to ka for analytical solution to one-compartment model.

# KA = CL/V
modlib("pk1") %>% 
  param(CL = 1, V = 20, KA = 0.05) %>% 
  ev(amt = 100) %>% 
  mrgsim() %>% 
  as_tibble %>% 
  pull(CENT) %>% 
  summary
#> Loading model from cache.
#> Error in eval(expr, envir, enclos): k10 is too close to ka for analytical solution to one-compartment model.

^{Created on 2024-05-24 with reprex v2.0.2}