Pharmacokinetics and ordinary differential equations
Since pharmacokinetics (PK) always involves time \(t\), we make it explicit in the population data formulation from the last chapter. To simplify presentation we only consider homoscedastic residual error and omit the subscripts: $$ y = f(x, \phi, t) + \epsilon. $$ In general the \(x\), \(y\), \(\phi\), and \(\epsilon\) are vectors.
For example, let us consider \(x=D\), \(D\) being the single bolus dose amount, and \(\theta=(V, k)\) with \(V\) for the volumne of distribution and \(k\) the elminination coefficient, in one-compartment models. For the plasma concentration observation \(y\), we have $$ y=\frac{D}{V}e^{-kt} +\epsilon, $$ thus the canonical form of \(f\) is $$ F=f(D, (V, k), t)=\frac{D}{V}e^{-kt} $$
Here \(F\) denotes the model prediction, the return of \(f\) that depends on dose, subject parameters, and time.
PRED routine
In general the user must provide function \(f()\), called the PRED routine in NONMEM. It takes THETA and each data record (encoding covariates and time), computes the prediction, and stores it in the reserved variable F.
We can specify any PRED by following the above signature. 1 For example, we can use PRED to describe a logistic \(f()\): $$ f(x, \theta, t) = f(t_0, (\alpha, \beta, \gamma), t) = \frac{\gamma}{1+\alpha e^{-\beta(t-t_0)}}, $$ and use it in a disease progression model to describe clinical outcomes that follow an S‑shaped time course.
In a NONMEM population model, we often encounter two types of predictions.
- Evaluation of \(f(x_{ij}, \phi_i, t)\) in which \(\phi_i\) is modeled by \(\mu(z_i, \theta, \eta_i)\). NONMEM calls this value individual prediction and denotes it IPRED.
- Evaluation of \(f(x_{ij}, \phi_i, t)\) at \(\eta_i=0\): \(\phi_i=\mu(z_i, \theta, \eta_i=0)\). That is, the prediction does not incorporate the inter-individual variation, by assuming identical prediction for subjects with same covariates. NONMEM calls this value PRED (not to be confused with the PRED routine). PRED is among NONMEM's default output.
Ordinary differential equations
Ordinary differential equations (ODEs) are often used in mechanistic PK models. In fact, the one-compartment model from earlier is based on an ODE model: our \(f(D, (V,k), t)\) is the solution of $$ \frac{df}{dt} = -kf(t), \quad f(t=0) = \frac{D}{V}. $$ In general, since the PK ODEs are constructed based on mass conservation, we would like to formulate them in terms of drug amount \(A(t)\), and write the above ODE as $$ \frac{dA}{dt} = -kA(t), \quad A(t=0) = A, $$ and recover the concentration as \(c(t) = A(t)/V\). This conversion between \(c(t)\) and \(A(t)\) is an example of scaling, with \(V\) being the scaling parameter here. Scaling is often used when the observation is concentration while the PK model is based on drug (or other substance) amount. Again, NONMEM has a reserved variable A corresponding to ODE states \(A\) in the model, S for scale parameters, and a control record $DES (Differential Equation Systems) to specifiy the ODEs.
Let us consider the two-compartment model with first-order absorption as another example. The ODE systsem is
$$ \frac{dA_0}{dt} = -k_aA_0(t),\\ \frac{dA_1}{dt} = k_aA_0(t) - (k_{12}+k_e)A_1(t) + k_{21}A_2(t),\\ \frac{dA_2}{dt} = k_{12}A_1(t)-k_{21}A_2(t). $$ Here \(A_0\), \(A_1\), \(A_2\) are the drug amount in depot (gut), central, and peripheral compartments, respectively, and the system has four coefficients: \(k_a\), \(k_{12}\), \(k_{21}\), \(k_e\), to quantify the absorption, the transfer between the central and peripheral compartments, as well as the elimination.
In general, the system can be formulated as
\begin{equation} \frac{dA(t)}{dt}=g(A(t), x, p, t). \end{equation}Here both left-hand-side (LHS) and right-hand-size (RHS) are vectors of dimention \(n\): \(A=[A_1, \dots, A_n]^T\), \(g=[g_1, \dots, g_n]^T\).
As the dimention of \(A()\) and the complexity of RHS increases to account for more complicated models it quickly becomes intractable to encode it in PRED, most ODEs do not have even close form to begin with. NONMEM tackles this problem by allowing user to specify the RHS \(g()\) directly, so that it can solve the ODEs numerically. In fact, NONMEM provides a library of commonly-used PK models that we need only to specify the parameters of \(g()\). Since these models describe how \(A(t)\) advances in time, NONMEM calls them ADVAN routines. We select ADVAN in the $SUBROUTINE control record.
General linear ODE model
A special case of the ODE models is when \(g()\) is linear function of A
\begin{equation} \frac{dA(t)}{dt}=K^TA(t). \end{equation}Here \(K\) is a \(n\times n\) matrix that could depend on \(x\). \(K(i, j)\) quantifies the first-order distribution of drug from compartment number i to compartment number j. For example, the above two-compartment model with first-order absorption can be cast as (unspecified elements have value 0): $$ A=[A_0, A_1, A_2]^T,\\ K_{12}=k_a,\quad K_{22}=-k_e,\quad K_{23}=k_{21},\quad K_{32}=k_{12} $$ NONMEM provides special subroutine ADVAN5 and ADVAN7 for this type of systems, and uses the above numbering scheme for matrix \(K\).
Parameterization
Another way to write the above two-compartment models is to use PK parameters such as clearance (\(\text{CL}\) and \(\text{Q}\)) and volumn of distribution (\(V_1\) and \(V_2\)). They are related to the \(k\)'s as (\(k_a\) remains the same) $$ k_e=\frac{\text{CL}}{V_1},\quad k_{12}=\frac{\text{Q}}{V_1},\quad k_{21}=\frac{\text{Q}}{V_2}. $$ NONMEM provides TRANS (for "transformation") routines for this reparameterization of the built-in PK models. We make selection of TRANS routine along with ADVAN routine, in the $SUBROUTINE control record.
Dose
It is helpful to see doses in PK studies from the pespective of ODEs, as dosing can be represented as either setting the initial condition of the ODEs, or amending the RHS function \(g()\). Specifically, consider a bolus dose \(D\) administored to the depot/gut compartment \(A_0\) at \(t=t_0\). This is equivalent to setting the initial condition $$ A_0(t_0) = D. $$ If instead \(D\) is administored in the central compartment, we have $$ A_1(t_0) = D. $$
As to an infusion with rate \(r\) given between time \(t_1\) and \(t_2\), it is equivalent to adding \(r\) to the RHS \(g()\). Then the two-compartment model with an infusion to the central compartment becomes (the other compartments remain the same) $$ \frac{dA_1}{dt} = k_aA_0(t) - (k_{12}+k_e)A_1(t) + k_{21}A_2(t) + r,\quad t \text{ in } (t_1, t_2). $$ during the infusion, and $$ \frac{dA_1}{dt} = k_aA_0(t) - (k_{12}+k_e)A_1(t) + k_{21}A_2(t),\quad t > t_2. $$ after the infusion.
A function signature specifies the function's input (arguments) and output (return). NONMEM's implementation of PRED actually involves additional arguments, but they are not essential here. In NONMEM Guides function signature is called function preface.