We’re indeed at a stage where we can start running benchmarks. We welcome suggestions, should people have problems they are interested in solving.
Our first example is a steady state model for pharmacokinetic. A patient takes a drug at a regular time interval: at the first the average drug mass in the patient’s system increases but after a few cycles it reaches a steady state. The drug mass y evolves over time according to the ODE
y_1' = - k_1 y_1 \\
y_2' = k_1 y_1 - k_2 y_2
with initial conditions (y_1^0, y_2^0). Often times, y_1^0 is the initial drug dose and y_2^0 = 0. The system is then readily solved by
y_1 = y_1^0 e^{-k_1 t} \\
y_2 = y_1^0 \frac{k_2}{k_1 - k_2} k_2(e^{-k_2t} - e^{-k_1 t}).
There is also an analytical solution when y_2^0 \neq 0, but it’s more complicated and I’m blinking on it, right now. We will need the latter solution when computing the system at steady state. Solving for the steady state means finding the initial conditions at the beginning of a cycle, such that we return to this state at the beginning of the next cycle. See https://charlesm93.github.io/files/2018-Margossian.pdf.
We measure the drug concentration, c, over time (for example over one or two cycles). Our likelihood distribution is, for example,
p(c \mid \theta) = \mathrm{Normal}\left (\frac{y_1}{V}, \sigma^2 \right),
where V is the volume of the central compartment and \theta the model parameters, here \theta = (k_0, k_1, V).
For many algorithms we care about, we need to compute
\nabla_\theta \log p(c \mid \theta).
We’ll do this with the adjoint and non-adjoint algebraic solvers.
The system can be scaled by introducing more patients and expanding the ODE. This ODE would only be pairwise coupled, but this can still help us learn about the scalability of our differentiation methods. I ran a similar experiment in the review for autodiff (https://arxiv.org/pdf/1811.05031.pdf, page 23).
Once we have this up and running, we can try a more challenging example where the ODE doesn’t admit an analytical solution and where we have more parameters and states per patient. I’m thinking a Friberg-Karlsson semi-mechanistic model – although I’m not sure patients would reach steady states in realistic applications.
Tagging some folks you might be interested in these experiments: @yizhang @wds15 @billg
Some values we can use for the simulation:
k1 = 2
k1 = 1
sigma = 1
y0 = 1000 // dose per pill / cycle
time = (0.083, 0.167, 0.25, 1, 2, 4)
ii = 5 // time between doses / cycles
y_obs = (58.1913, 104.996, 140.298, 206.566, 120.457, 23.1845)
The observed values were not simulated at steady state, but that shouldn’t really matter, since we’re primarily concerned in evaluating the derivatives, not doing a full model fit.
