Torsten: A Pharmacokinetic/Pharmacodynamic Model Library for Stan
User Manual
(Torsten Version 0.84, Stan version 2.17.1)

Table of Contents

Development team

Bill Gillespie

billg@metrumrg.com, Metrum Research Group, LLC

Yi Zhang

yiz@metrumrg.com, Metrum Research Group, LLC

Charles Margossian

charles.margossian@columbia.edu, Columbia University, Department of Statistics(formerly Metrum Research Group, LLC)

Acknowledgements

Institutions

We thank Metrum Research Group, Columbia University, and AstraZeneca.

Funding

This work was funded in part by the following organizations:

Office of Naval Research (ONR) contract N00014-16-P-2039

provided as part of the Small Business Technology Transfer (STTR) program. The content of the information presented in this document does not necessarily reflect the position or policy of the Government and no official endorsement should be inferred.

Bill & Melinda Gates Foundation.

Individuals

We thank the Stan Development Team for giving us guidance on how to create new Stan functions and adding features to Stan's core language that facilitate building ODE-based models.

We also thank Kyle Baron and Hunter Ford for helpful advice on coding in C++ and using GitHub, Curtis Johnston for reviewing the User Manual, and Yaming Su for using Torsten and giving us feedback.

1 Introduction

Stan is an open source probabilistic programing language designed primarily to do Bayesian data analysis carpenter17_stan. Several of its features make it a powerful tool to specify and fit complex models. Notably, its language is extremely flexible and its No U-Turn Sampler (NUTS), an adaptative Hamiltonian Monte Carlo algorithm, has proven more efficient than commonly used Monte Carlo Markov Chains (MCMC) samplers for complex high dimensional problems hoffman_no-u-turn_2011. Our goal is to harness these innovative features and make Stan a better software for pharmacometrics modeling. Our efforts are twofold:

  1. We contribute to the development of new mathematical tools, such as functions that support differential equations based models, and implement them directly into Stan's core language.
  2. We develop Torsten, an extension with specialized pharmacometrics functions.

Throughout the process, we work very closely with the Stan Development Team. We have benefited immensely from their mentorship, advice, and feedback. Just like Stan, Torsten is an open source project that fosters collaborative work. Interested in contributing? Shoot us an e-mail and we will help you help us(billg@metrumrg.com)!

Torsten is licensed under the BSD 3-clause license.

1.1 Installing Torsten

Installation files are available on GitHub

We are working with Stan development team to create a system to add and share Stan packages. In the mean time, the current repo contains forked version of Stan with Torsten. The latest version of Torsten (v0.84) is compatible with Stan v2.17.1. Torsten is agnostic to which Stan interface you use. Here we provide command line and R interfaces.

First, download the project. The root path of the project is your torsten_path. Set the envionment variable TORSTEN_PATH to torsten_path. In bash this is

export TORSTEN_PATH=torsten_path

1.1.1 Command line version

There is no need to install command line version. To build a Stan model with model_name in model_path using command line, in torsten_path/cmdstan, do

make model_path/model_name

1.1.2 R version

The R version is based on rstan, the Stan's interface for R. To install R version of Torsten, at torsten_path, in R

source('install.R')

Please ensure the R toolchain includes a C++ compiler with C++11 support. In particular, R 3.3.0 and later is recommended as it contains toolchain based on gcc 4.9.3. On Windows platform, such a toolchain can be found in Rtools33 and later.

Please ensure CXXFLAGS in .R/Makevars constains flag -std=c++11.

For more information of installation troubleshooting, please consult rstan wiki.

1.2 Overview

Torsten is a collection of Stan functions to facilitate analysis of pharmacometric data using Stan. The current version includes:

  • Specific linear compartment models:
    • One compartment model with first order absorption.
    • Two compartment model with elimination from and first order absorption into central compartment
  • General linear compartment model described by a system of first-order \underline{linear} Ordinary Differential Equations (ODEs).
  • General compartment model described by a system of first order ODEs.
  • Mix compartment model with PK forcing function described by a linear one or two compartment model.

The models and data format are based on NONMEM\textregistered{}1/NMTRAN/PREDPP conventions including:

  • Recursive calculation of model predictions
    • This permits piecewise constant covariate values
  • Bolus or constant rate inputs into any compartment
  • Handles single dose and multiple dose histories
  • Handles steady state dosing histories
    • Note: The infusion time must be shorter than the inter-dose interval.
  • Implemented NMTRAN data items include: TIME, EVID, CMT, AMT, RATE, ADDL, II, SS

In general, all real variables may be passed as Stan parameters. A few exceptions apply /to functions which use a numerical integrator/(i.e. the general and the mix compartment models). The below listed cases present technical difficulties, which we expect to overcome in Torsten's next release:

  • The RATE and TIME arguments must be fixed
  • In the case of a multiple truncated infusion rate dosing regimen:
    • The bioavailability (F) and the amount (AMT) must be fixed.

This library provides Stan language functions that calculate amounts in each compartment, given an event schedule and an ODE system.

1.3 Implementation details

  • Stan version 2.17.1
  • All functions are programmed in C++ and are compatible with the Stan math automatic differentiation library carpenter15_stan_math_librar
  • One and two compartment models: hand-coded analytical solutions
  • General linear compartment models with semi-analytical solutions using the built-in matrix exponential function
  • General compartment models with numerical solutions using built-in ODE integrators in Stan. The tuning parameters of the solver are adjustable. The steady state solution is calculated using a numerical algebraic solver.
  • Mix compartment model: the PK forcing function is solved analytically and the forced ODE system is solved numerically.

1.4 Development plans

Our current plans for future development of Torsten include the following:

  • Build a system to easily share packages of Stan functions (written in C++ or in the Stan language)
  • Allow numerical methods to handle RATE, AMT, TIME, and the bioavailability fraction (F) as parameters in all cases.
  • Optimize Matrix exponential functions
    • Function for the action of Matrix Exponential on a vector
    • Hand-coded gradients
    • Special algorithm for matrices with special properties
  • Fix issue that arises when computing the adjoint of the lag time parameter (in a dosing compartment) evaluated at \(t_{\text{lag}} = 0\).
  • Extend formal tests
    • We want more C++ Google unit tests to address cases users may encounter
    • Comparison with simulations from the R package mrgsolve and the software NONMEM\textregistered{}
    • Recruit non-developer users to conduct beta testing

1.5 Updates since Torsten 0.83

  • Torsten is now up to date with Stan version 2.17.1.
  • Add piecewise linear interpolation function.
  • Add univariate integral functions.
  • Minor revisions to User Manual.

2 Using Torsten

The reader should have a basic understanding of how Stan works before reading this chapter. There are excellent resources online to get started with Stan

In this section we go through the different functions Torsten adds to Stan. It will be helpful to apply these functions to a simple example. We have uploaded code and data on

We use the following example to demonstrate use of several functions in the Torsten library.

Two compartment model

We model drug absorption in a single patient and simulate plasma drug concentrations:

  • Multiple Doses: 1250 mg, every 12 hours, for a total of 15 doses
  • PK: plasma concentrations of parent drug (\(c\))
  • PK measured at 0.083, 0.167, 0.25, 0.5, 0.75, 1, 1.5, 2, 4, 6, 8, 10 and 12 hours after 1st, 2nd, and 15th dose. In addition, the PK is measured every 12 hours throughout the trial.

The plasma concentration (\(c\)) are simulated according to the following equations:

\begin{eqnarray*} \log\left(c\right) &\sim& N\left(\log\left(\widehat{c}\right),\sigma^2\right) \\ \widehat{c} &=& f_{2cpt}\left(t, CL, Q, V_2, V_3, k_a\right) \\ \left(CL, Q, V_2, V_3, ka\right) &=& \left(5\ {\rm L/h}, 8\ {\rm L/h}, 20\ {\rm L}, 70\ {\rm L}, 1.2\ {\rm h^{-1}} \right) \\ \sigma^2 &=& 0.01 \end{eqnarray*}

and the drug concentration is given by \(c = y_2/V_2\).

where the mass of drug in the central compartment (\(y_2\)) is obtained by solving the system of ordinary differential equations (ODEs):

\begin{align*} y_1' &= -k_a y_1 \\ y_2' &= k_a y_1 - \left(\frac{CL}{V_2} + \frac{Q}{V_2}\right) y_2 + \frac{Q}{V_3} y_3 \\ y_3' &= \frac{Q}{V_2} y_2 - \frac{Q}{V_3} y_3 \end{align*}

The data are generated using the R package mrgsolve Baron000.

2.1 One Compartment Model

real[] = PKModelOneCpt(real[] time, real[] amt, real[] ate,
                       real[] ii, int[] evid, int[] cmt,
                       real[] addl, int[] ss, real[] theta,
                       real[] biovar, real[] tlag)

Parameters in theta : \(CL\), \(V_2\), \(k_a\).

The event arguments. time, amt, rate, ii, evid, cmt, addl, and ss, describe the event schedule of the clinical trial. All arrays have the same length, which corresponds to the number of events.

The model arguments, theta contains the ODE parameters, biovar the bioavailability fraction in each compartment, and tlag the lag time in each compartment. The model arguments may be either one or two dimensional arrays. If they are one dimensional arrays, the parameters are constant for all events. If they are two dimensional arrays then each row contains the parameters for the interval [ time[i-1], time[i] ]. The number of rows should equal the number of events.

Setting \(ka\) to 0 eliminates the first-order absorption.

cptModels.png

Figure 1: One and two compartment models with first order absorption implemented in Torsten.

2.2 Two Compartment Model

real[] = PKModelTwoCpt(real[] time, real[] amt, real[] ate,
                       real[] ii, int[] evid, int[] cmt,
                       real[] addl, int[] ss, real[] theta,
                       real[] biovar, real[] tlag)

Parameters in theta : \(CL\), \(Q\), \(V_2\), \(V_3\), \(k_a\)

See section 1 for function description.

Code example below shows Stan example of using TwoCptModel to fit Two compartment model. Three MCMC chains of 2000 iterations were simulated. The first 1000 iteration of each chain were discarded. Thus 1000 MCMC samples per chain were used for the subsequent analyses.

data{
  int<lower = 1> nt;  // number of events
  int<lower = 1> nObs;  // number of observation
  int<lower = 1> iObs[nObs];  // index of observation

  // NONMEM data
  int<lower = 1> cmt[nt];
  int evid[nt];
  int addl[nt];
  int ss[nt];
  real amt[nt];
  real time[nt];
  real rate[nt];
  real ii[nt];

  vector<lower = 0>[nObs] cObs;  // observed concentration (Dependent Variable)
}

transformed data{
  vector[nObs] logCObs = log(cObs);
  int nTheta = 5;  // number of ODE parameters in Two Compartment Model
  int nCmt = 3;  // number of compartments in model

  // Since we're not trying to evaluate the bio-variability (F) and 
  // the lag times, we declare them as data.
  real biovar[nCmt];
  real tlag[nCmt];

  biovar[1] = 1;
  biovar[2] = 1;
  biovar[3] = 1;

  tlag[1] = 0;
  tlag[2] = 0;
  tlag[3] = 0;

}

parameters{
  real<lower = 0> CL;
  real<lower = 0> Q;
  real<lower = 0> V1;
  real<lower = 0> V2;
  real<lower = 0> ka;
  real<lower = 0> sigma;

}

transformed parameters{
  real theta[nTheta];  // ODE parameters
  vector<lower = 0>[nt] cHat;
  vector<lower = 0>[nObs] cHatObs;
  matrix<lower = 0>[nt, nCmt] x; 

  theta[1] = CL;
  theta[2] = Q;
  theta[3] = V1;
  theta[4] = V2;
  theta[5] = ka;

  // PKModelTwoCpt takes in the NONMEM data, followed by the parameter
  // arrays abd returns a matrix with the predicted amount in each 
  // compartment at each event.
  x = PKModelTwoCpt(time, amt, rate, ii, evid, cmt, addl, ss,
                   theta, biovar, tlag);

  cHat = col(x, 2) ./ V1; // we're interested in the amount in the second compartment 

  cHatObs = cHat[iObs]; // predictions for observed data recors

The MCMC history plots(Figure 2) suggest that the 3 chains have converged to common distributions for all of the key model parameters. The fit to the plasma concentration data (Figure 5) are in close agreement with the data, which is not surprising since the fitted model is identical to the one used to simulate the data. Similarly the parameter estimates summarized in Table 1 and Figure 4 are consistent with the values used for simulation.

Table 1: Summary of the MCMC simulations of the marginal posterior distribu- tions of the model parameters
  mean semean sd 2.5% 25% 50% 75% 97.5% neff Rhat
CL 4.823 0.002 0.092 4.647 4.762 4.823 4.883 5.012 2392.155 1.00
Q 7.596 0.013 0.586 6.479 7.201 7.594 7.977 8.785 1923.939 1.00
V1 21.073 0.069 2.573 16.017 19.352 21.046 22.817 26.097 1385.883 1.00
V2 76.365 0.105 5.611 65.805 72.623 76.172 79.916 87.971 2862.184 1.00
ka 1.231 0.004 0.177 0.907 1.107 1.221 1.344 1.599 1581.825 1.00
sigma 0.109 0.000 0.012 0.089 0.100 0.108 0.116 0.134 2560.112 1.00

../example-models/R/deliv/figure/TwoCptModel/TwoCptModelPlots001.pdf

../example-models/R/deliv/figure/TwoCptModel/TwoCptModelPlots002.pdf

../example-models/R/deliv/figure/TwoCptModel/TwoCptModelPlots003.pdf

../example-models/R/deliv/figure/TwoCptModel/TwoCptModelPlots006.pdf

2.3 General Linear ODE Model Function

real[] = linOdeModel(real[] time, real[] amt, real[] rate,
                     real[] ii, int[] evid, int[] cmt,
                     real[] addl, int[] ss,
                     matrix K, real[] biovar, real[] tlag)

A general linear ODE model refers to a model that may be described in terms of a system of first order linear differential equations with (piecewise) constant coefficients, i.e., a differential equation of the form:

\begin{equation} y^\prime\left(t\right) = Ky\left(t\right) \end{equation}

where \(K\) is a matrix. For example \(K\) for a two compartment model (equation \eqref{eq:TwoCpt}) with first order absorption is:

\begin{equation} K = \left[\begin{array}{ccc} -k_a & 0 & 0 \\ k_a & -\left(k_{10} + k_{12}\right) & k_{21} \\ 0 & k_{12} & -k_{21} \end{array}\right] \end{equation}

where \(k_{10}=CL/V_2\), \(k_{12}=Q/V_2\), and \(k_{21}=Q/V_3\).

System parameters is in K, with K being the constant rate matrix is the same for all events or an array of constant rate matrices. The length of the array is the number of events and each element corresponds to the matrix at the interval [ time[i-1], time[i] ]. Note that K contains all the ODE parameters, so we no longer need theta.

The following Stan example illustrate the use of General linear model for fitting a two compartment model with first order absorption.

// LinTwoCptModelExample.stan
// Run two compartment model using matrix exponential solution
// Heavily anotated to help new users

data{
  int<lower = 1> nt;  // number of events
  int<lower = 1> nObs;  // number of observations
  int<lower = 1> iObs[nObs];  // index of observation

  // NONMEM data
  int<lower = 1> cmt[nt];
  int evid[nt];
  int addl[nt];
  int ss[nt];
  real amt[nt];
  real time[nt];
  real rate[nt];
  real ii[nt];

  vector<lower = 0>[nObs] cObs;  // observed concentration (dependent variable)
}

transformed data{
  vector[nObs] logCObs = log(cObs);
  int nCmt = 3;
  real biovar[nCmt];
  real tlag[nCmt];

  for (i in 1:nCmt) {
    biovar[i] = 1;
    tlag[i] = 0;
  }
}

parameters{
  real<lower = 0> CL;
  real<lower = 0> Q;
  real<lower = 0> V1;
  real<lower = 0> V2;
  real<lower = 0> ka;
  real<lower = 0> sigma;

}

transformed parameters{
  matrix[3, 3] K;
  real k10 = CL / V1;
  real k12 = Q / V1;
  real k21 = Q / V2;
  vector<lower = 0>[nt] cHat;
  vector<lower = 0>[nObs] cHatObs;
  matrix<lower = 0>[nt, 3] x;

  K = rep_matrix(0, 3, 3);

  K[1, 1] = -ka;
  K[2, 1] = ka;
  K[2, 2] = -(k10 + k12);
  K[2, 3] = k21;
  K[3, 2] = k12;
  K[3, 3] = -k21;

  // linModel takes in the constant rate matrix, the object theta which
  // contains the biovariability fraction and the lag time of each compartment,
  // and the NONMEM data.
  x = linOdeModel(time, amt, rate, ii, evid, cmt, addl, ss,
                  K, biovar, tlag);

2.4 General ODE Model Function

real[] model_name(ODE_system, int nCmt,
                  real[] time, real[] amt, real[] rate, real[] ii,
                  int[] evid, int[] cmt, real[] addl, int[] ss,
                  real[] theta, real[] biovar, real[] tlag,                      
                  real rel_tol, real abs_tol, int max_step)

Torsten may be used to fit models described by a system of user-specified first-order ODEs, i.e., differential equations of the form:

\begin{equation*} y'(t) = f(t, y(t)) \end{equation*}

In the case where the rate vector \(R\) is non-zero, this equation becomes:

\begin{equation*} y'(t) = f(t, y(t)) + R \end{equation*}

ODE_system specifies \(f(t, y(t))\), which the user defines inside the functions block (see section 19.2 of the Stan reference manual for details and Figure~\ref{GenTwoCptModelCode} for an example). The user does NOT include the rates in their definition of \(f\). Torsten automatically corrects the derivatives when the rates are non-zero.

nCmt is the number of compartments (or, equivalently, the number of ODEs) in the model. rel_tol, abs_tol, and max_step are tuning parameters for the ODE integrator: respectively the relative tolerance, the absolute tolerance, and the maximum number of steps.

The options for model_name are:

  • generalOdeModel_rk45
  • generalOdeModel_bdf

They respectively call the built-in Runge-Kutta 4th/5th order (rk45) integrator, recommended for non-stiff ODEs, and the Backward Differentiation (BDF) integrator, recommended for stiff ODEs. Which value to use for the tuning parameters depends on the integrator and the specifics of the ODE system. Reducing the tolerance parameters and increasing the number of steps make for a more robust integrator but can significantly slow down the algorithm. The following can be used as a starting point:

  • rel_tol = 1e-6
  • abs_tol = 1e-6
  • max_step = 1e+6

for rk45 integrator and

  • rel_tol = 1e-10
  • abs_tol = 1e-10
  • max_step = 1e+8

for the bdf integrator 2. Users should be prepared to adjust these values. For additional information, see the Stan User's Manual stan_team_2017.

A few notable restrictions apply to :

  • rate and time cannot be passed as parameters.
  • In the case of a multiple truncated infusion rate dosing regimen:
    • The bioavailability biovar) and the amount amt) cannot be passed as parameters.

These restrictions also apply to mixOdeCpt_* functions, discussed in the next section.

// GenTwoCptModelExample.stan
// Run two compartment model using numerical solution
// Heavily anotated to help new users

functions{

  // define ODE system for two compartmnt model
  real[] twoCptModelODE(real t,
                        real[] x,
                        real[] parms,
                        real[] rate,  // in this example, rate is treated as data
                        int[] dummy){

    // Parameters
    real CL = parms[1];
    real Q = parms[2];
    real V1 = parms[3];
    real V2 = parms[4];
    real ka = parms[5];

    // Re-parametrization
    real k10 = CL / V1;
    real k12 = Q / V1;
    real k21 = Q / V2;

    // Return object (derivative)
    real y[3];  // 1 element per compartment of
                // the model

    // PK component of the ODE system
    y[1] = -ka*x[1];
    y[2] = ka*x[1] - (k10 + k12)*x[2] + k21*x[3];
    y[3] = k12*x[2] - k21*x[3];

    return y;
  }

}
data{
  int<lower = 1> nt;  // number of events
  int<lower = 1> nObs;  // number of observations
  int<lower = 1> iObs[nObs];  // index of observation

  // NONMEM data
  int<lower = 1> cmt[nt];
  int evid[nt];
  int addl[nt];
  int ss[nt];
  real amt[nt];
  real time[nt];
  real rate[nt];
  real ii[nt];

  vector<lower = 0>[nObs] cObs;  // observed concentration (dependent variable)
}

transformed data{
  vector[nObs] logCObs = log(cObs);
  int nTheta = 5;   // number of parameters
  int nCmt = 3;   // number of compartments
  real biovar[nCmt];
  real tlag[nCmt];

  for (i in 1:nCmt) {
    biovar[i] = 1;
    tlag[i] = 0;
  }
}

parameters{
  real<lower = 0> CL;
  real<lower = 0> Q;
  real<lower = 0> V1;
  real<lower = 0> V2;
  real<lower = 0> ka;
  real<lower = 0> sigma;
}

transformed parameters{
  real theta[nTheta];
  vector<lower = 0>[nt] cHat;
  vector<lower = 0>[nObs] cHatObs;
  matrix<lower = 0>[nt, 3] x; 

  theta[1] = CL;
  theta[2] = Q;
  theta[3] = V1;
  theta[4] = V2;
  theta[5] = ka;

  // generalCptModel takes in the ODE system, the number of compartment 
  // (here we have a two compartment model with first order absorption, so
  // three compartments), the parameters matrix, the NONEM data, and the tuning
  // parameters (relative tolerance, absolute tolerance, and maximum number of steps)
  // of the ODE integrator. The user can choose between the bdf and the rk45 integrator.
  // Returns a matrix with the predicted amount in each compartment 
  // at each event.

//  x = generalOdeModel_bdf(twoCptModelODE, 3,
//                          time, amt, rate, ii, evid, cmt, addl, ss,
//                          theta, biovar, tlag,
//                          1e-8, 1e-8, 1e8);

   x = generalOdeModel_rk45(twoCptModelODE, 3,
                           time, amt, rate, ii, evid, cmt, addl, ss,
                           theta, biovar, tlag,
                           1e-6, 1e-6, 1e6);

  cHat = col(x, 2) ./ V1;

  for(i in 1:nObs){
    cHatObs[i] = cHat[iObs[i]];  // predictions for observed data records
  }
}

2.5 Mixed ODE Model Function

real[] model_name(reduced_ODE_system, int nOde,
                  real[] time, real[] amt, real[] rate,
                  real[] ii, int[] evid, int[] cmt, real[]
                  addl, int[] ss,
                  real[] theta, real[] biovar, real[] tlag,
                  real rel_tol, real abs_tol, real max_step)

In certain cases, an ODE system can be divided in two subsystems:

\begin{align*} y_1^\prime &= f_1(t, y_1) \\ y_2^\prime &= f_2(t, y_1, y_2) \end{align*}

where \(y_1\), \(y_2\), \(f_1\), and \(f_2\) are vector-valued functions, and \(y_1^\prime\) is independent of \(y_2\). This structure arises in PK/PD models, where \(y_1\) describes a forcing PK function and \(y_2\) the PD effects. If \(y_1\) has an analytical solution, we can construct a \textit{mixed solver}, which analytically solves \(y_1\) and numerically integrates \(y_2\). This approach leads to an appreciable gain in computational efficiency. In the example of a Friberg-Karlsson semi-mechanistic model, we observe an average speedup of \(\sim 47 \pm 18 \%\) when using the mix solver in lieu of the numerical integrator. Torsten supports the mixed solver for cases where \(y_1\) solves the ODEs for a One or Two Compartment model with a first-order absorption.

The reduced_ODE_system specifies the system we numerically solve (\(y_2\) in the above discussion, also called the reduced system and nOde the number of equations in the \underline{reduced} system. The function that defines a reduced system has an almost identical signature to that used for a full system, but takes one additional argument: \(y_1\), the PKstates, i.e. solution to the PK ODEs.

real[] reducedODE(real t,  // time
                  real[] y,  // reduced state
                  real[] y1,  // PK states}'
                  real[] theta,  // parameters
                  real[] x_r,  // data (real)
                  int[] x_int)  // data (integer)

The options for modelName are:

  • mixOde1CptModel_rk45
  • mixOde1CptModel_bdf
  • mixOde2CptModel_rk45
  • mixOde2CptModel_bdf

These four functions correspond to all the permutations we can obtain when using a forcing One or Two Compartment function, and the Runge-Kutta 4th/5th order (rk45) or Backward Differentiation (BDF) integration method. The mixed ODE functions can be used to compute the steady state solutions supported by the general ODE model functions.

Restrictions regarding which arguments may be passed as parameters for generalOdeModel_* also apply to mixOdeCptModel_*.

We cannot apply the mixed solver to the Two Compartment example we have been using so far. Instead, we will consider the model which motivated the implementation of the method in the first place.

2.6 Friberg-Karlsson Semi-Mechanistic Model

In this second example, we add to our two Compartment model a PD effect, described by a system of nonlinear ODEs friberg_mechanistic_2003.

Friberg-Karlsson Model

Neutropenia is observed in patients receiving an ME-2 drug. Our goal is to model the relation between neutrophil counts and drug exposure. Using a feedback mechanism, the body maintains the number of neutrophils at a baseline value (Figure~\ref{FK}). While in the patient's blood, the drug impedes the production of neutrophils. As a result, the neutrophil count goes down. After the drug clears out, the feedback mechanism kicks in and brings the neutrophil count back to baseline.

\begin{align} \log(ANC_i) \sim& N(\log(Circ), \sigma^2_{ANC}) \\ Circ =& f_{FK}(MTT, Circ_{0}, \alpha, \gamma, c) \\ (MTT, Circ_{0}, \alpha, \gamma, ktr) =& (125, 5.0, 3 \times 10^{-4}, 0.17) \\ \sigma^2_{ANC} =& 0.001 \end{align}

where \(c\) is the drug concentration in the blood we get from the Two Compartment model, and \(Circ\) is obtained by solving the following system of nonlinear ODEs:

\begin{eqnarray} \begin{aligned} y_\mathrm{prol}' &= k_\mathrm{prol} y_\mathrm{prol} (1 - E_\mathrm{drug})\left(\frac{Circ_0}{y_\mathrm{circ}}\right)^\gamma - k_\mathrm{tr}y_\mathrm{prol} \\ y_\mathrm{trans1}' &= k_\mathrm{tr} y_\mathrm{prol} - k_\mathrm{tr} y_\mathrm{trans1} \\ y_\mathrm{trans2}' &= k_\mathrm{tr} y_\mathrm{trans1} - k_\mathrm{tr} y_\mathrm{trans2} \\ y_\mathrm{trans3}' &= k_\mathrm{tr} y_\mathrm{trans2} - k_\mathrm{tr} y_\mathrm{trans3} \\ y_\mathrm{circ}' &= k_\mathrm{tr} y_\mathrm{trans3} - k_\mathrm{tr} y_\mathrm{circ} \end{aligned} \label{eq:FK} \end{eqnarray}

where \(E_{drug} = \alpha c\).

The ODEs specifying the Two Compartment Model (equation~\ref{eq:TwoCpt}) do not depend on the PD ODEs (equation~\ref{eq:FK}) and can be solved analytically by Torsten. We therefore specify our model using a mixed solver function. We do not expect our system to be stiff and use the Runge-Kutta 4th/5th order integrator (Figures \ref{FK_reduced} and \ref{FK_mix2}).

functions{
  real[] FK_ODE(real t,
                real[] y,
                real[] y_pk,
                real[] theta,
                real[] rdummy,
                int[] idummy){
    /* PK variables */
    real VC = theta[3];

    /* PD variable */
    real mtt      = theta[6];
    real circ0    = theta[7];
    real alpha    = theta[8];
    real gamma    = theta[9];
    real ktr      = 4.0 / mtt;
    real prol     = y[1] + circ0;
    real transit1 = y[2] + circ0;
    real transit2 = y[3] + circ0;
    real transit3 = y[4] + circ0;
    real circ     = fmax(machine_precision(), y[5] + circ0);
    real conc     = y_pk[2] / VC;
    real EDrug    = alpha * conc;

    real dydt[5];

    dydt[1] = ktr * prol * ((1 - EDrug) * ((circ0 / circ)^gamma) - 1);
    dydt[2] = ktr * (prol - transit1);
    dydt[3] = ktr * (transit1 - transit2);
    dydt[4] = ktr * (transit2 - transit3);
    dydt[5] = ktr * (transit3 - circ);

    return dydt;
  }
}

data{
  int<lower = 1> nt;
  int<lower = 1> nObsPK;
  int<lower = 1> nObsPD;
  int<lower = 1> iObsPK[nObsPK];
  int<lower = 1> iObsPD[nObsPD];
  real<lower = 0> amt[nt];
  int<lower = 1> cmt[nt];
  int<lower = 0> evid[nt];
  real<lower = 0> time[nt];
  real<lower = 0> ii[nt];
  int<lower = 0> addl[nt];
  int<lower = 0> ss[nt];
  real rate[nt];
  vector<lower = 0>[nObsPK] cObs;
  vector<lower = 0>[nObsPD] neutObs;

  real<lower = 0> circ0Prior;
  real<lower = 0> circ0PriorCV;
  real<lower = 0> mttPrior;
  real<lower = 0> mttPriorCV;
  real<lower = 0> gammaPrior;
  real<lower = 0> gammaPriorCV;
  real<lower = 0> alphaPrior;
  real<lower = 0> alphaPriorCV;
}

transformed data{
  int nOde = 5;
  vector[nObsPK] logCObs;
  vector[nObsPD] logNeutObs;
//  int idummy[0];
//  real rdummy[0];

  int nTheta;
  int nIIV;

  int n;                        /* ODE dimension */
  real rtol;
  real atol;
  int max_step;

  n = 8;
  rtol = 1e-8;
  atol = 1e-8;
  max_step = 100000;

  logCObs = log(cObs);
  logNeutObs = log(neutObs);

  nIIV = 7; // parameters with IIV
  nTheta = 9; // number of parameters
}

parameters{

  real<lower = 0> CL;
  real<lower = 0> Q;
  real<lower = 0> VC;
  real<lower = 0> VP;
  real<lower = 0> ka;
  real<lower = 0> mtt;
  real<lower = 0> circ0;
  real<lower = 0> alpha;
  real<lower = 0> gamma;
  real<lower = 0> sigma;
  real<lower = 0> sigmaNeut;

  // IIV parameters
  cholesky_factor_corr[nIIV] L;
  vector<lower = 0>[nIIV] omega;
}

transformed parameters{
  vector[nt] cHat;
  vector<lower = 0>[nObsPK] cHatObs;
  vector[nt] neutHat;
  vector<lower = 0>[nObsPD] neutHatObs;
  real<lower = 0> theta[nTheta];
  matrix[nt, nOde + 3] x;
  real biovar[nTheta];
  real tlag[nTheta];

  for (i in 1:nTheta) {
    biovar[i] = 1.0;
    tlag[i] = 0.0;
  }

  theta[1] = CL;
  theta[2] = Q;
  theta[3] = VC;
  theta[4] = VP;
  theta[5] = ka;
  theta[6] = mtt;
  theta[7] = circ0;
  theta[8] = alpha;
  theta[9] = gamma;

  x = mixOde2CptModel_rk45(FK_ODE, nOde, time, amt, rate, ii, evid, cmt, addl, ss, theta, biovar, tlag, rtol, atol, max_step);

  cHat = col(x, 2) / VC;
  neutHat = col(x, 8) + circ0;

  for(i in 1:nObsPK) cHatObs[i]    = cHat[iObsPK[i]];
  for(i in 1:nObsPD) neutHatObs[i] = neutHat[iObsPD[i]];

}

model {

  // Priors
  CL    ~ normal(0, 20);
  Q     ~ normal(0, 20);
  VC    ~ normal(0, 100);
  VP    ~ normal(0, 1000);
  ka    ~ normal(0, 5);
  sigma ~ cauchy(0, 1);

  mtt       ~ lognormal(log(mttPrior), mttPriorCV);
  circ0     ~ lognormal(log(circ0Prior), circ0PriorCV);
  alpha     ~ lognormal(log(alphaPrior), alphaPriorCV);
  gamma     ~ lognormal(log(gammaPrior), gammaPriorCV);
  sigmaNeut ~ cauchy(0, 1);

  // Parameters for Matt's trick
  L ~ lkj_corr_cholesky(1);
  omega ~ cauchy(0, 1);

  // observed data likelihood
  logCObs ~ normal(log(cObs), sigma);
  logNeutObs ~ normal(log(neutObs), sigmaNeut);
}

\appendix \printindex \backmatter

\bibliography{torsten}

Footnotes:

1
NONMEM\textregistered{} is licensed and distributed by ICON Development Solutions.
2
These are the default tuning parameters the integrators. Torsten functions do not have a default values for these parameters. The user must explicitly pass the tuning parameters to generalOdeModel_*() .

Author: Yi Zhang

Created: 2018-03-22 Thu 00:01

Validate