Divergent transitions in PK models



I am using one compartmental PK model with 5 subjects (simulated data) and still get divergent transitions while using non-centered parametrization. Any suggestions how to get rid of them? Pairs plot shows that population k forms a funnel with sd of k (the same for V).

Thank you,

vanco1cmprestlog.stan (3.0 KB)


That population k is in the generated quantities, so that doesn’t affect sampling or divergence—the generated quantities get calculated once after the Hamiltonian simulation.

You may want to impose stronger priors than Cauchy—those can let things run off into the tails where you can run into divergences because of arithmetic issues.

If your data heavily constrains the posterior (i.e., lots of data or very narrowly distributed data), the centered parameterization works better.

Maybe @betanalpha or @wds15 will recognize the model and have some advice.


are you simulating from this model, or did you use some other software to simulate?

One thought is, if you simulated with another software and have real looking profieles, I would suggest re-checking your equations - they do not look right.

Your single dose is flooring concentration to 0 when t is less than td, and likewise your multidose is automatically negating any accumulation by setting c = 0. Vancomycin you’re going to be targeting accumulation to at least a trough of 10 mg/L / AUC/MIC ratio of 400, so the scenario your showing shouldn’t have full washout, thus my concern over not handling accumulation.

For multidose equations you should be either using the analytical solution with the dosing interval baked in, or you must incrementally save and update based on the accumulation by the next dose, which you are not doing (that I could see).

Why don’t you use torsten for this? Just use the one compartment analytical solution it provides, and is made specifically for PK applications?


A bit against the usual advise to use the non centered parametrization, i found that in quite some of my applications, the cp parametrization perfoms better. This seem to happen whenever you see the pk profiles clearly from what i have seen.



I actually simulated from this model in stan and checked with PKPDsim R
library. Since I have a first order elimination I am using superposition
of each administration. So the single dose function calculates profile
just for a single dose administration while multidose function simply
adds those profiles.

I actually thought about using torsten but the data I have is not just
plasma concentration but it also includes BUN/serum creatinine etc. So I
thought that making my one model may be simpler.

Maybe you could help with the math when k changes during the multidose
treatment - basically I am trying to model renal function improvement
during the treatment what means that vancomycin clearance changes as
well. I have a bunch of labs to play with: BUN/serum creatinine/urine
creatinine/cystatin C and bit overwhelmed how to reflect all this data
in the changing clearance. I am planning to use available correlations
but perhaps there is a better way. The goal is to construct predictive
posterior for a population.



do not use k. Distinguishing physiological changes back to a rate constant with capacity for extrapolation is a wild goose chase. Use physiological parameterizations - CL and V.

Vancomycin is 70+% cleared by glomerular filtration, and changes in clearance are tied to changes in GFR - so use EGFR as your covariate.

I don’t know if you are in peds, elderly, or just people in renal failure, but the rest of your model can reflect the population a little bit as well.

For vancomycin I would allometrically scale WT on CL/V with EGFR on CL as well either linearly or using allometry (0.75).

If you have access to nonmem, here is a neonatal model I’ve developed that is robust to 7-fold changes in serum creatinine over time

$PROBLEM   ID: 143 fold:  NOBS: 3/3
$DATA ../../mdata/id_143_nobs_3_fold_04.csv 



D1 = 1 ; 1 hour infusion durations

; constants 
NORM_WT = 2.9
NORM_PMA = 34.8

;; because log transforming all thetas, want to keep the exponent in a consistent domain


V = EXP(TVV + ETA(2))
S1 = V

( 0.1 , 0.16824 ) ; Typical Value Clearance for 2.9 kg baby with EGFR of 30
( 0.1 , 1.76117 ) ; Typical Value Volumei for 2.9 kg baby
( 0.1 , 0.77429 ) ; exponent for EGFR~clearance

0.16824 FIX; Typical Value Clearance
1.76117 FIX; Typical Value Volume
0.77429 FIX; exponent for EGFR~clearance

3.708500e-06 FIXED 
7.174766e-06 6.119516e-04 
-4.559485e-07 1.959912e-05 1.742600e-04
2.910000e-03 2.316000e-02
4.613000e-02 FIXED 
2.910000e-03 2.316000e-02
4.810000e-03 1.195230e+00
2.918000e-02 FIXED 
4.810000e-03 1.195230e+00

Y = IPRED*(1 + ERR(1)) + ERR(2)

"97 FORMAT(I12,1X,F14.0,10(1X,1PG12.5))



I’m working with @charlesm93 to translate it into torsten and stan, but both of us are tight on time.

For example, red (observed samples) fs the predicted profiles.

Hope this will help, but be careful using this parameterization exactly if you are in adults or older peds - I have not done an QA/QC outside of neonates, adults with renal insufficiency could have a different set of parameters that might be better pulled from literature.


In the moment when k changes, the elimination rate, then the superposition principle breaks down. If the changes of the renal function is slow over time then the easiest way to deal with it is to say that k is a stepwise constant function. The way to solve things is stepwise and take advantage of the fact that the odes solutions are specified by the initial. You have to make reasonable stepsizes in comparison to the change rate of k.

This works… i have fit a time changing clearwnce 2cmt model myself a while ago.


Thanks a lot for great suggestions. Yes, I assumed that k (or clearance as
suggested before) is a stepwise function. I tried to implement in stan but
it became very nasty since superposition is not anymore. Do you happen to
have the codes (for 2cmt model)? Maybe torsten handles covariates changing
over the time?

The data includes patients for all age groups so intention is to consider
separate populations (neonates/pediatric/adults/elderly). It is very
sparse for serum concentrations but there is a decent amount of labs
(creatinine clearance/BUN/albumin serum/total protein/age/weight/WBC/etc)
so thought is to implement those as covariates that influence GFR and
consequently vancomycin clearance.

The ultimate goal is to determine population predictive posterior to be
used for loading doses.


Awesome fit! If interested I would be glad to help. I think I have enough
familiarity in NONMEM to transform your model to stan - perhaps we will
have to use different priors and cholesky decomposition to handle
correlation. I am not sure about torsten since it lacks IV model. The
challenging part is handling time changing covariates.


torsten can handle IV infusions I believe, the biggest thing it does give you is an event system, such that should also (in theory) allow easier management of time-varying covariates.

Here is a repo that @charlesm93 and I started to do an example of an IV infusion. We’ve done a bit of work offline, and he has it on his computer. I’ll get him to push it up https://github.com/dpastoor/infusion-stan

Feel free to hack away at it, I would be very interested to see/test a stan implementation. I look forward to the day I can get off nonmem so I can use all the cores I want :-)


By the way when I moved to CL V parametrization and informative priors the
divergencies disappeared!


It looks like Pred1_oneCpt has the code for IV - after intelligent
replacement of ka with the rate of infusion. Unfortunately I had some make
errors when compiling torsten example so probably will proceed with plain
stan model. Having real data (at least for 5+ subjects) would be helpful -
otherwise I just use simulated data.


Stan model works fine for simulated data for 5 patients. Black - true
concentration, green dash - 95% CI, green - median, red - observed. The
model may be improved in following ways:

  • cholesky for correlated parameters
  • non-centered parametrization
  • follow torsten framework more closely

I guess I could improve the model on items 1&2 if more realistic data is
supplied. R and stan models are attached.

vanco1cmp.r (8.52 KB)

vanco1cmppoplog1.stan (2.93 KB)

vanco1cmppoprestlog1.stan (4.09 KB)


awesome! I didn’t take too close of a look - I am turning in my dissertation in a little over a week and am currently 980 pages in, and still have more writing, but I will circle back to this for sure once I turn it in.

Charles pushed up some of the changes - its still not there yet, but this also has a mixed residual error, and is using torsten

I’ll also be dockerizing torsten in August so it is a ‘one-click’ install for anyone that wants to try out that branch, so maybe you can give it a whirl then.


Good luck with thesis. It is funny that my structure is quite similar to
yours. It looks like that you have everything except EGFR. It would be
great if you could sent me some data to play with so I can improve the stan
model so eventually we can move to torsten.


simulate my dude!

overseer and PKPDmisc are personal packages, so

install_github(“dpastoor/”) to use them

If you can get those fit, I can throw you some real profiles and see how it performs, but I warn you they are nasty - there is a sample individual

ID_422_rep_10_nobs_10.csv (4.2 KB)


Actually I would like to start with the nasty data. I am already full
with the simulated data (frankly I want WT and EGFR as covariates). I
have a question about the IV start/finish time. In the csv file AMT=30
and RATE=-2. If IV starts at 0 when does it end? At 15? But the next
dosing event is at 7.44. How does it work? Maybe you can drop
id_143_nobs_3_fold_04.csv? At least I will have NONMEM’s answer to
compare with…


You may want to consider reducing page count rather than increasing it!


Attached is the updated model and results of with simulated data (WT and
EGFR). I have followed the structure you have used while leaving non
relevant items such as ii or ss out. Would like to crank data for 5+
Black line - simulated concentration
Green line - median of predicted concentration
Blue dash - 95% CI of predicted concentration
Red - concentration values used for fit

The first stan model simulates the data, the second stan model does the fit
and makes prediction

vanco1cmppop.r (8.58 KB)

vanco1cmppoplog1.stan (2.41 KB)

vanco1cmppoprestlog1.stan (3.47 KB)


Hi all,
I’m hopping in a bit late, but I’ve been following.

Maybe torsten handles covariates changing over the time?

Yes. When passing the theta argument to torsten, you may elect to pass a 2D array to reflect piece-wise constant parameters which change from one event to the other.

I am not sure about torsten since it lacks IV model.

You can use PKModelOneCpt and specify ka = 0.

It looks like Pred1_oneCpt has the code for IV - after intelligent replacement of ka with the rate of infusion. Unfortunately I had some make errors when compiling torsten example so probably will proceed with plain stan model.

There was indeed a bug in torsten. I fixed it while working on Devin’s model this Friday, but the fix is in the develop, not the released version. Using ka as the rate infusion is a neat trick, and speaks to the flexibility of compartment models.

In the csv file AMT=30 and RATE=-2. If IV starts at 0 when does it end? At 15? But the next dosing event is at 7.44. How does it work?

If rate = 2, then yes, the infusion ends at 15 hours. (I’m not sure about rate = -2). You then get overlapping infusions, since a new infusion starts before the previous one ends. I suspect Torsten should be able to handle this fine, but I haven’t tried it out yet. There is a notable problem if the system is at steady state, because we need to count the number of overlapping infusions, which introduces a discrete parameter – but SS doesn’t seem to be relevant here.

Also, how often, during clinical trials is amt / rate > ii, i.e infusion time longer than the interdose interval?