lsoda error when simulating steady-state continuous infusion with reduced bio-availability

[rfaelens]

When simulating an ODE-type model with steady-state continuous infusion, mrgsolve stops simulating due to an LSODA error.

[lsoda] tout too close to t to start integration

[mrgsolve] lsoda returned with negative istate: -3
  illegal input detected (see printed message).

Error: simulation terminated.

Minimal reproducible example:

library(mrgsolve)
library(tidyverse)

m1 <- mrgsolve::mcode("myModel", '
$PARAM 
CL = 1.1, V = 20
BIOAV = 1

$CMT CENT  

$MAIN
F_CENT = BIOAV;

$ODE
dxdt_CENT =  -CL/V*CENT;
')
dataset <- tibble(ID = 1, time = 0, evid = 1, amt = 100, cmt = 1, 
     ss = 1, ii = 10, rate = 4.12, BIOAV = 0.412)
m1 %>% 
  data_set(dataset) %>% 
  mrgsim_df()

[kylebaron]

Thanks for the report, @rfaelens; I confirmed the behavior.

The critical detail in the input is that BIOAV changes the duration of the infusion to the dosing interval ii. Will look into how these records are sorted or otherwise processed.

[kylebaron]

This requires reduced bioav; doesn't happen when infusion duration is same as dosing interval with bioav==1

library(mrgsolve)
#> 
#> Attaching package: 'mrgsolve'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(tidyverse)

m1 <- mrgsolve::mcode("myModel", '
$PARAM 
CL = 1.1, V = 20
BIOAV = 1

$CMT CENT  

$MAIN
F_CENT = BIOAV;

$ODE
dxdt_CENT =  -CL/V*CENT;
')
#> Building myModel ...
#> done.

dataset <- tibble(ID = 1, time = 0, evid = 1, amt = 100, cmt = 1, 
                  ss = 1, ii = 10, rate = 10, BIOAV = 1)
m1 %>% 
  data_set(dataset) %>% 
  mrgsim() %>% plot()

^{Created on 2022-03-22 by the reprex package (v2.0.1)}

[kylebaron]

The root cause was directly comparing t and tout, which came out different when adjacent records should have been at the same time.

• Refactor how that comparison was done
• I also write some code to just run the (true) zero-order constant infusion routine whenever ii was about equal to the infusion duration; this runs the system to steady state for infusion with rate; then we have to dose one more intermittent infusion and then turn it off; this worked, but the code would need some additional refactoring to make sure that terminal dose happened "next" .... just passing more information through so the records could be properly sorted.

[kylebaron]

Example

library(mrgsolve)
#> 
#> Attaching package: 'mrgsolve'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(tidyverse)

m1 <- mrgsolve::mcode("myModel", '
$PARAM 
CL = 1.1, V = 20
BIOAV = 1

$CMT CENT  

$MAIN
F_CENT = BIOAV;

$ODE
dxdt_CENT =  -CL/V*CENT;
')
#> Building myModel ...
#> done.

dataset <- tibble(ID = 1, time = 0, evid = 1, amt = 100, cmt = 1, 
                  ss = 1, ii = 10, rate = 4.120, BIOAV = 0.4120)
m1 %>% 
  data_set(dataset) %>% 
  mrgsim() %>% plot()

^{Created on 2022-03-22 by the reprex package (v2.0.1)}

[kylebaron]

@rfaelens Here's a workaround in the meantime; this will get fixed in the next release. Kyle

library(mrgsolve)
#> 
#> Attaching package: 'mrgsolve'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(tidyverse)

m1 <- mrgsolve::mcode("myModel", '
$PARAM 
CL = 1.1, V = 20
BIOAV = 1

$CMT CENT  

$MAIN
F_CENT = BIOAV;

$ODE
dxdt_CENT =  -CL/V*CENT;
')
#> Building myModel ...
#> done.

data <- c(
  ev(amt = 0,   ss = 1,  rate = 4.12, BIOAV = 1), 
  ev(amt = 100, ii = 10, rate = 4.12, BIOAV = 0.412)
)

m1 %>% 
  data_set(data) %>% 
  mrgsim() %>% plot()

^{Created on 2022-03-23 by the reprex package (v2.0.1)}