Bounding posterior output of pharmacokinetic model


I would add up variances and not standard deviations…

1 Like

@wds15 You are obviously correct, it totally should have been real sigma = sqrt(square(sigma_absolute) + square(sigma_relative * mu[n])). This was based on my advice so thanks for noting and sorry for the confusion - assuming my older code is definitely OK and giving it as example is probably not a good policy going forward :-). (and I definitely hope @dbarneche checks what I propose, I definitely make a lot of mistakes).

Other than that looks roughly good, some more minor things I noticed:

Those are very wide priors, especially for sigma_relative, I think you can safely assume that a priori the relative error is not more than 1 (i.e. the error is not 100%), giving you something like sigma_relative ~ normal(0, 0.5); or even tighter. Similarly, knowing the error of your measuring device could give you better prior for the sigma_absolute.

This is also a bit weird - the constraint will only make the sampler throw errors at you when it is invalidated, and break the sampler if non-negligible posterior mass is outside the constraint. Are you sure Eff cannot physically be larger than 1? (I’ve discussed that in a previous post here, but not sure if you changed the model or you’ve seen a mistake in my reasoning).

1 Like

Thanks again @martinmodrak and @wds15! Even after adding the corrections above (i.e. adding variances rather than standard deviations, tighter priors on sigma), the peak is still slightly underestimated when compared to the model without an absolute error term only.

The crux of the problem is the efficiency equation. I came up with this efficiency reasoning myself borrowing some concepts from ecological theory (my field of research), so while I’d say it makes lots of sense to me, I admit I could be missing something.

If we start with concentration \mu = 0 at time T_0, and the dynamics are dictated by this single dose, first order absorption model, the only parameter in the equation that conceptually drives the drug absorption (as far as I understand!) is K_a. K_e and C_l are parameters related to elimination processes.

So, if the conceptual interpretation above is correct, K_a T sets the maximum possible concentration at any time point (even in the hypothetical absence of elimination processes). So, whatever the concentration value at the peak P is (or any concentration at any time point), it has to be lower than K_a T, even if mathematically higher values of Q (originally D/C_l) allow for E > 1.

Please let me know if you disagree. I’m not a pharmacologist, so I am going to reach out to more experienced folks to confirm this reasoning is correct!


I thought so for a while, but that is NOT correct:

D = 10
Ka = 1
Cl = 1
Ke = 0.1
tibble(t = seq(0,10,length.out = 100)) %>% 
  mutate(mu = D * Ka * (exp(-Ke * t) - exp(-Ka * t)) / (Cl * (Ka - Ke))) %>%
  ggplot(aes(t, mu)) + geom_line() + geom_abline(slope = Ka, color = "blue", linetype = "dashed")

The dashed line is K_a T, the curve is \mu(T)

1 Like

Hi @martinmodrak,

Thanks for this. As mentioned before, I appreciate that mathematically you can get cases where E > 1. In your plot above, for instance, you only see that because D was very high relative to C_l. But if you were to set D = 1 for example, you’d get E < 1 over the entire range of t.

I wonder two things:

1 - My proposed equation is correct, and all I am missing are the defined natural constraints of this equation that would set the values of E to always be [0,1]. I haven’t been able to find some in-depth discussion on the web about what all these parameters mean, and how (or if) they constrain one another.

2 - If there’s another efficiency equation for E that explicitly limits its value to be [0,1] regardless of the parameter space. In my mind, even if the E formula I proposed is wrong, the question is still valid, i.e. what was the efficiency of the organism in absorbing the drug up to the peak? Efficiencies by definition are never > 1. I’m not a mathematician though, so I wouldn’t know if there’s a method that allows me to derive such an equation that is independent of the parameter space. If you have ideas though, I’d love to hear them!

Thanks again!


Turns out I was using wrong formula (withoute K_e in the numerator). With the correct one, you can make E arbitrarily large even for fixed D/C_l, say:

D = 1
Ka = 1
Cl = 0.5
Ke = 10

T_p = (log(Ka) - log(Ke)) / (Ka - Ke)
P_rel = Ke * Ka * (exp(-Ke*T_p) - exp(-Ka*T_p)) / (Ka - Ke)
P = P_rel * D / Cl

tibble(t = seq(0,1,length.out = 100)) %>% 
  mutate(mu = D *Ke * Ka * (exp(-Ke * t) - exp(-Ka * t)) / (Cl * (Ka - Ke))) %>%
  ggplot(aes(t, mu)) + geom_line() + geom_abline(slope = Ka, color = "blue", linetype = "dashed")

And note that when modelling you generally shouldn’t put hard constraints based on domain knowledge (e.g that D < C_l), this should rather be represented as a soft constraint in prior. Hard constraints generally belong to unphysical/mathematically impossible situations, which E > 1 isn’t.