Non-finite gradient in pharmacokinetic model with complex absorption

Modeling
, ,
1

Hi all,

I’m working on a pharmacokinetic model using rstan. From the initial development I’ve implemented a first-order absorption model using both the analytic and ODE solutions. Both perform fine, but the data suggest an absorption delay is required. A lag-time model performs poorly (irreperable divergences and low ESS) presumably due to numeric stability issues around the transition point.

Instead I’ve implemented a transit compartment model (https://doi.org/10.1007/s10928-007-9066-0) which in principle implements a delayed absorption without an abrupt transition. The system has bolus dosing so I used the analytic solution as in the paper. No matter what I do with it, I receive a ‘non-finite gradient at initial value’ error on attempting to sample from the model.

After some troubleshooting:

  1. there are no apparent problems with ODE performance; the model appears to be correct as I can simulate from it with Stan, and get equivalent results from deSolve. I used the same initial values. I used the print() function within the Stan code which confirms that the model predicted values appear realistic.
  2. the ‘grad_log_prob()’ function shows that indeed the mean-transit-time parameter has NaN gradient at any initial value (except very small values in which CVOde fails), but for any values I select, the log_prob apparently exists. I also fixed the ‘number of compartments’ parameter so I think there shouldn’t be an identifiability issue.
  3. don’t seem to be any programming problems; the model converges just fine if I replace the absorption model with a first-order model with otherwise exactly the same code.

Has anyone made a similar model and encountered this or have any comment about what I might be missing?

1 Like
2

@wds15 @Funko_Unko

3

Maybe something like here is happening? Modelling spatiotemporal random effects results in - #17 by martinmodrak

But you say it works with a slightly different model, so maybe not?

2 Likes
4

@Funko_Unko after much digging, I didn’t make a diagnosis but I think the source of this issue is numeric performance of the ODE.

The paper recommends exp(log()) of the transit compartment term to improve numeric stability, which I had implemented, but I wonder if this might have been causing underflow somewhere, so removed it. I rewrote the ODE system and the caller function and set the absolute tolerance much lower. So far sampling without error (albeit slowly).

Hopefully remains stable when I run the hierarchical version!

2 Likes