PK models in R and in Julia

1.1 Pharmacokinetic models with an analytic solution

Pharmacokinetics is the study of the absorption and elimination of drugs and their metabolites in the body. As described in the wikipedia article there are several parameters, such rate constants, volumes of distribution and clearances, that we wish to estimate from pharmacokinetic data, which usually consists of measurements of blood serum concentrations following certain doses at certain times. My particular interest is in analysis of population pharmacokinetic data using nonlinear mixed-effects models and in the design of experiments for population pharmacokinetic studies.

For many of the models based on linear elimination kinetics the concentration at time t after a dose d has an analytic solution. As part of the PFIM software developed in France Mentre's lab at Inserm, Julie Bertrand created a catalogue of these expressions for 1, 2 or 3 compartment models after oral or bolus or infusion administration of a single dose or a fixed number of doses or in steady state after multiple doses. Julie and I created the PKPDmodels package for R to provide these expressions as functions that also evaluate the gradient with respect to the parameters and, optionally, the Hessian.

At the recent PODE: Population Optimum Design of Experiments conference organized by France, the wonderful Emmanuelle Comets mentioned that in her saemix package for R she only needs the predicted concentration values. She asked if it was slowing things down to evaluate the gradient at every function evaluation. I decided to benchmark.

1.2 Model function and benchmarks in R

One of the (regrettably few) published examples of population pharmacokinetic data is from a study of Theophylline kinetics. These data are available as the Theoph data set in R. The experimental subjects were given a single oral dose of the drug. Assuming a 1-compartment model with linear elimination kinetics the expression for the model is

> library(PKPDmodels)
> PKexpr("oral", "sd")
~(dose/V) * (ka/(ka - k)) * (exp(-k * t) - exp(-ka * t))

in terms of the parameters V, the volume of distribution, ka, the absorption rate constant, and k, the elimination rate constant.

It turns out that the combination of V and k is not a particularly good choice of parameters as their estimates usually end up being highly correlated. A better choice is V and Cl, the clearance, which is defined as Cl = kV. We can achieve this parameterization by providing a list of transformation expressions from the new parameters to the old. Nothing sophisticated is going on, just a substitution of expressions in another expression. While we are at it we could switch to the logarithms of pharmacokinetic parameters which provide a log-likelihood function that is closer to being quadratic.

> PKexpr("oral", "sd", list(ka ~ exp(lka), V ~ exp(lV), k ~ exp(lCl)/V))
~(dose/exp(lV)) * (exp(lka)/(exp(lka) - exp(lCl)/V)) * (exp(-(exp(lCl)/V) * 
    t) - exp(-exp(lka) * t))

A call to PKmod with the same arguments produces a byte-compiled function that evaluates both the function and the gradient.

> (mfg <- PKmod("oral", "sd", list(ka ~ exp(lka), V ~ exp(lV), k ~ exp(lCl - lV))))
function (dose, t, lV, lka, lCl) 
{
    .expr1 <- exp(lV)
    .expr2 <- dose/.expr1
    .expr3 <- exp(lka)
    .expr5 <- exp(lCl - lV)
    .expr6 <- .expr3 - .expr5
    .expr7 <- .expr3/.expr6
    .expr8 <- .expr2 * .expr7
    .expr11 <- exp(-.expr5 * t)
    .expr14 <- exp(-.expr3 * t)
    .expr15 <- .expr11 - .expr14
    .expr19 <- .expr8 * (.expr11 * (.expr5 * t))
    .expr21 <- .expr6^2
    .expr23 <- .expr2 * (.expr3 * .expr5/.expr21)
    .value <- .expr8 * .expr15
    .grad <- array(0, c(length(.value), 3L), list(NULL, c("lV", 
        "lka", "lCl")))
    .grad[, "lV"] <- .expr19 - (.expr23 + dose * .expr1/.expr1^2 * 
        .expr7) * .expr15
    .grad[, "lka"] <- .expr2 * (.expr7 - .expr3 * .expr3/.expr21) * 
        .expr15 + .expr8 * (.expr14 * (.expr3 * t))
    .grad[, "lCl"] <- .expr23 * .expr15 - .expr19
    attr(.value, "gradient") <- .grad
    .value
}
<bytecode: 0x2bd7c88>

Because the symbolic differentiation in R to produce the gradient also performs common subexpression elimination this function produces a cleaner evaluation of the model function and a comparatively compact expression of the gradient.

For comparison I use two types of evaluations of the model function itself, one with a nested function call and one with the transformations performed internally.

> mf0 <- cmpfun(function(dose, t, V, ka, k) (dose/V) * (ka/(ka - k)) * (exp(-k * t) - exp(-ka * t)))
> mf1 <- cmpfun(function(dose, t, lV, lka, lCl) mf0(dose, t, exp(lV), exp(lka), exp(lCl - lV)))
> mf2 <- cmpfun(function(dose, t, lV, lka, lCl) {V <- exp(lV); ka <- exp(lka); k <- exp(lCl - lV); (dose/V) * (ka/(ka - k)) * (exp(-k * t) - exp(-ka * t))})

Now would be a good time to check that I haven't made transcription mistakes

> Dose <- Theoph$Dose
> str(Dose)
 num [1:132] 4.02 4.02 4.02 4.02 4.02 4.02 4.02 4.02 4.02 4.02 ...
> Time <- Theoph$Time
> lV <- rep.int(-1, length(Dose))
> lka <- rep.int(0.6, length(Dose))
> lCl <- rep.int(-4, length(Dose))
> all.equal(mf1(Dose, Time, lV, lka, lCl), mf2(Dose, Time, lV, lka, lCl))
[1] TRUE
> all.equal(mf1(Dose, Time, lV, lka, lCl), as.vector(mfg(Dose, Time, lV, lka, lCl)))
[1] TRUE

Finally we benchmark

> library(rbenchmark)
> cols <- c("test", "elapsed", "relative", "user.self", "sys.self")
> benchmark(mf1(Dose, Time, lV, lka, lCl), mf2(Dose, Time, lV, lka, lCl), mfg(Dose, Time, lV, lka, lCl), replications=1000L, columns=cols, order="elapsed")
                           test elapsed relative user.self sys.self
2 mf2(Dose, Time, lV, lka, lCl)   0.081 1.000000     0.080        0
1 mf1(Dose, Time, lV, lka, lCl)   0.084 1.037037     0.084        0
3 mfg(Dose, Time, lV, lka, lCl)   0.219 2.703704     0.220        0

So the conclusion is that there is very little difference between the nested function calls and the explicit assignment of transformed values and that it takes 2.7 times as long to evaluate both the fitted values and the gradient than does evaluation of only the fitted values. However, the actual execution times for 1000 function evaluations are sufficiently small that this operation should not be a bottleneck.

1.3 Translation into Julia

I decided to translate this code into Julia, partially to gain experience with the language. I have said that the most valuable character trait for a programmer is unbounded pessimism because you spend your working time thinking, "now what can go wrong here". Because of this tendency I was careful to distinguish scalar and vector arguments and to test for correct sizes of vector arguments, etc. A bit of testing showed me that this was unnecessary. Write it the simple way and let the functions that you call check for numeric arguments and vector-scalar operations and correct vector-vector operations, which they will do.

The only changes I needed to make were to replace '*' by '.*' and '/' by './'. The "dot" forms of these operators perform componentwise vector-vector operations. Because they are also defined for scalar-scalar operations and scalar-vector operations there is no harm in using them throughout.

The result is anticlimactic. The function definitions are

mf0(dose, t, V, ka, k) = (dose ./ V) .* (ka ./ (ka - k)) .* (exp(-k .* t) - exp(-ka .* t))
mf1(dose, t, lV, lka, lCl) = mf0(dose, t, exp(lV), exp(lka), exp(lCl - lV))

function mfg(dose, t, lV, lka, lCl)
    V      = exp(lV)
    expr2  = dose ./ V
    ka     = exp(lka)
    Cl     = exp(lCl)
    k      = Cl ./ V
    expr6  = ka - k
    expr7  = ka ./ expr6
    expr8  = expr2 .* expr7
    expr11 = exp(-k .* t)
    expr14 = exp(-ka .* t)
    expr15 = expr11 - expr14
    expr18 = V .* V
    expr19 = Cl .* V ./ expr18
    expr24 = expr6 .* expr6
    expr8 .* expr15, hcat(expr8 .* (expr11 .* (expr19 .* t)) - (expr2 .* (ka .* expr19 ./ expr24) + dose .* V ./ expr18 .* expr7) .* expr15,
                          expr2 .* (expr7 - ka .* ka ./ expr24) .* expr15 + expr8 .* (expr14 .* (ka .* t)),
                          expr2 .* (ka .* k ./ expr24) .* expr15 - expr8 .* (expr11 .* (k .* t)))
end

and timing the same calculations on the same machine produced

elapsed time: 0.08809614181518555 seconds
elapsed time: 0.17201519012451172 seconds

which are very close to the timings for the byte-compiled functions in R.