### 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.

I see that its been close to 4 years since this blog post was written, somehow missed it. Julia made some big strides and so did R in this time. Was/Is there any interest in developing a fast and flexible open source PKPD non-linear mixed effects modeling system in Julia or R that can replace NONMEM?

ReplyDeleteDo you know of any groups interested in developing or willing to collaborate on such a venture!