Estimación de los parámetros farmacocinéticos poblacionales de disposición

Método basado en los parámetros de la función característica

library(knitr)
library(magrittr)
library(gridExtra)
library(nlme)
library(nlraa)

Introducción

Parámetros fundamentales de disposición

• En lo sucesivo se designarán parámetros fundamentales de disposición el aclaramiento plasmático (\(Cl\)), volumen aparente de distribución inicial o del compartimento central (\(V_1\)), volumen aparente de distribución en estado estacionario (\(V_{SS}\)), volumen aparente de distribución asintótico (\(V_{AS}\)) y vida media (\(t_{50}\)). Las definiciones matemáticas de los cuatro primeros se recogen en la primera columna de la tabla 2; \(E[\theta]\) es el valor esperado (media) del tiempo de residencia del fármaco en el sistema.
• Obviamente, si la farmacocinética del fármaco, en ciertas condiciones experimentales, es interpretable mediante el modelo monocompartimental \(V_1~=~V_{SS}~=~V_{AS}\).
• Los parámetros \(V_{SS}\), \(V_{AS}\) y \(t_{50}\) no son estimables directamente, sino que deben ser calculados a partir de los parámetros de la función utilizada para ajustar los datos experimentales. La estimación directa de \(Cl\) y \(V_1\) depende de la parametrización del modelo farmacocinético.

Antecedentes

• Wagner (1976) publicó las fórmulas para el cálculo de los parámetros básicos de disposición a partir de los parámetros de una función poliexponencial (coeficientes y exponentes) estimados utilizando los datos de concentración plamática de fármaco - tiempo.
• Wagner y colaboradores (1977) analizaron los resultados obtenidos en la estimación de los parámetros de disposición de 20 fármacos. Los autores:
  • Identificaron el error en el cálculo de los parámetros farmacocinéticos utilizando un modelo con administración IV-bolus cuando la administración es por perfusión intravenosa a velocidad constante (\(k_0\), dimensiones \(m~t^{-1}\)) durante un corto plazo de tiempo (\(t_0\)).
  • Proponen utilizar para el análisis de los datos experimentales la función

\[\begin{equation} c_p(t)~=~\sum_{i=1}^m~A_i^´~e^{-\lambda_i~(t-t_0)} ~~~~~~t \ge t_0 \tag{1} \end{equation}\]

y calcular los coeficientes que se hubieran obtenido tras la administración IV-bolus tras la administración de una dosis igual a \(k_0 \times t_0\) (ecuaciones 1 - 3 del artículo original; ver ecuación 8 más adelante).

Problemas a resolver en este tutorial

• Estimación de los parámetros fundamentales de diposición y la variabilidad interindividual a partir de la función característica de disposición utilizando un método bootstrap.
• No tenemos información de que se haya propuesto algún criterio para decidir si el tiempo de perfusión es suficientemente corto como para interpretar la administración por perfusión contínua como si hubiera sido una administración IV-bolus. En la sección Tiempo de perfusión y selección del modo de administración se propone como criterio \(t_0 \times \lambda_1 \le 0.05\).

Descripción de este tutorial

Modelo teórico

• Se utiliza el esquema propuesto por DiStefano y Landaw (1984a, 1984b) (ver figura 1) para interpretar la farmacocinética del fármaco.
• La disposición del fármaco se interpreta a través de la función característica.
• Las ecuaciones para los modelos farmacocinéticos lineales más habituales se obtienen por la convolución de la función de entrada por la función característica. Ver las funciones de entrada y las ecuaciones resultantes en la tabla 1.
• Los parámetros de disposición fundamentales se calculan a partir de los parámetros de la función característica (ver columna 3 de la tabla 2).

Estimación de los parámetros de disposición poblacionales

• Los parámetros de la función característica es estiman mediante regresión no lineal de efectos mixtos utilizando la función nlme::nlme().
• Los parámetros de disposición fundamentales se calculan mediante el método bootstrap implementado en nlraa::boot_nlme().

Modelo teórico

La figura 1 recoge los esquemas del modelo para el análisis farmacocinético no compartimental propuesto por DiStefano y Landaw (1984a, 1984b); las líneas de trazo contínuo representa flujo de masa; las de trazo discontínuo observación.

 

Figura 1. Esquema propuesto por DiStefano y Landaw para el análisis farmacocinético no compartimental

 

Integral de convolución

La integral de convolución es un teorema fundamental en el estudio de los sistemas lineales e invariantes en el tiempo. En esta sección se aplica al desarrollo de modelos farmacocinéticos; en lo sucesivo supondremos que la función de interés es la concentración plasmática de fármaco \(c_p(t)\).

Los modelos farmacocinéticos lineales, tanto compartimentales como no compartimentales, para distintos modos de administración, pueden derivarse a partir de la convolución de la función de entrada \(\left(I(t) \right)\) con la función característica \(\left( G(t) \right)\):

\[\begin{equation} c_p(t)~=~(I~*~G)(t)~=~\int_{\tau=0}^t I(t)~G(t- \tau)~d \tau \tag{2} \end{equation}\]

El operador \(*\) se conoce como producto de convolución.

Función de entrada

• La función de entrada expresa la velocidad de entrada de fármaco en el sistema, por lo que sus dimensiones son \(m~t^{-1}\).
• \(I(t)\) puede ser una función lineal o no lineal, contínua o definida a trozos.
• Las tres funciones de entrada de uso más frecuente en farmacocinética son las que describen la administración IV-bolus, la administración por perfusión contínua, y la absorción de primer orden. Dado que en general la fracción de dosis absorbida (\(F\)) es deconocida cuando la administración es extravasal, solo los dos primeros modos de administración son adecuados para la estimación de los parámetros de disposición.
• Se utilizan otras funciones más complejas en el diseño y evaluación de formas de cesión sostenida y en el análisis de las correlaciones in vitro - in vivo (Gillespie, 1997)

La segunda columna de la tabla 2 recoge las funciones de entrada para los tres modos de administración estudiados habitualmente en farmacocinética. \(\delta(t)\) es la función delta de Dirac, y \(U(t~-~t_0)\) la función escalón de Heaviside.

Función característica

La función característica es la respuesta del sistema cuando la entrada es un impulso unitario (adimensional). En contexto desarrollado en este tutorial las dimensiones de la función característica es \(l^{-3}\).

• Otras denominaciones de la función característica: función de ponderación (teoría de sistemas lineales); función respuesta al impulso unitario (por ejemplo, Gomeni y Bressolle-Gomeni, 2021). La transformada de Laplace de la función característica recibe el nombre de función de transferencia.
• La función característica tiene que ser necesariamente lineal e invariable en el tiempo. Por lo tanto la ecuación 2 no es aplicable a fármacos con farmacocinética no lineal o afectados por el ciclo circadiano.

Función característica en el análisis no compartimental

Aunque se han propuesto otras alternativas (van Rossum y de Bie, 1989), en farmacocinética se utiliza casi exclusivamente la función poliexponencial,

\[\begin{equation} G(t)~=~\sum_{i=1}^m~a_j~e^{- \lambda_i~t}\tag{3} \end{equation}\]

• Las dimensiones de los coeficientes a_i son los de la función característica; los exponentes λ_i tienen por dimensión t^-1.
• Por convención, los términos de la ecuación 3 se ordenan de forma que λ₁ λ₂ > … λ_n.
• Los parámetros de la ecuación 3 se estiman experimentalmente.

Función característica de los modelos compartimentales

La función característica se deriva resolviendo el sistema de ecuaciones diferenciales homogéneo asumiendo que la entrada en el sistema es un impulso unitario adimensional. Por ejemplo, para el modelo lineal de dos compartimentos y utilizando la nomenglatura habitual en los textos de farmacocinética, la función característica es:

\[\begin{equation} G(t)~=~\frac{1}{V_1}~\frac{\alpha - k_{12}}{\alpha - \beta}~e^{-\alpha~t}~+~\frac{1}{V_1}~\frac{k_{12} - \beta}{\alpha - \beta}~e^{-\beta{t}} \tag{4} \end{equation}\]

Cálculo numérico de la función característica de los modelos compartimentales

La forma más sencilla y directa de calcular numéricamente la funcion característica es el método de las matrices. Para una descripción detallada y su implementación en R ver Resolución numérica de modelos compartimentales lineales. Ver el ejemplo 2.

Equivalencia entre el análisis no compartimental y el análisis compartimental

Los modelos derivados de ambos métodos de análisis farmacocinético son equivalentes si tienen la misma función característica.

Modelos farmacocinéticos

La tabla 1 resume los modelos para la administración IV-bolus, administración IV-perfusión contínua y absorción de primer orden. Obsérvese que en estas ecuaciones no aparece de forma explícita el volumen aparente de distribución del compartimento central.

 

Tabla 1. Modelos farmacocinéticos habituales

Modo administración	\(I(t)\)	\(c_p(t)\)
IV - bolus	\(\delta(t)~X_0\)	\(c_p(t)~=~X_0 \sum_{j=1}^m a_j~e^{-\lambda_j~t}\)
Perfusión contínua IV	\(k_0 \left((1-U(t-t_0) \right)\)	\(c_p(t)~=~k_0 \left(\sum_{j=1}^m \frac{a_j}{\lambda_j} \left(1-e^{-\lambda_j~t} \right)-U(t-t_0)~\sum_{j=1}^m \frac{a_j}{\lambda_j} \left(1-e^{-\lambda_j~(t-t_0)} \right) \right)\)
Absorción orden 1	\(k_a~F~X_0~e^{-k_a~t}\)	\(c_p(t)~=~F~X_0~k_a \left(\sum_{i=1}^m \frac{a_i}{k_a-\lambda_i}e^{-\lambda_i~t} + e^{-k_a~t} \sum_{i=1}^m \frac{a_i}{\lambda_i - k_a}\right)\)

 

Administración IV-bolus

Obsérvar

\[\begin{equation} A_i~=~X_0~a_i \tag{5} \end{equation}\]

y por lo tanto las dimensiones de \(A_i\) son \(m~l^{-3}\).

Administración por perfusión contínua

La ecuación de niveles plasmáticos - tiempo tras la administración por perfusión contínua se simplifica según sea \(t \le t_0\) (fase de perfusión) o \(t > t_0\) (fase de post-perfusión):

• Fase de perfusión (\(t \le t_0\)):

\[\begin{equation} c_p(t)~=~k_0~\sum_{i=1}^m~\frac{a_i}{\lambda_i} \left(1-e^{-~\lambda_i~t} \right) \tag{6} \end{equation}\]

• Fase post-perfusión(\(t > t_0\));

\[\begin{equation} c_p(t)~=~k_0~\sum_{i=1}^m~\frac{a_i}{\lambda_i} \left(e^{\lambda_i~t_0} - 1 \right)~e^{-~\lambda_i~t} \tag{7} \end{equation}\]

Jusfificación de la ecuación 1

En la ecuación propuesta por Wagner (1977) para interpretar la curva de niveles plasmáticos tras la administración del fármaco por perfusión contínia IV (ecuación 1),

\[\begin{equation} A'_i~=~k_0~\frac{a_i}{\lambda_i}~\left(e^{\lambda_i~t_0}~-~1 \right) \tag{8} \end{equation}\]

Tiempo de perfusión y selección del modo de administración

La función \(c_p(t)\) tras la administración del fármaco por perfusión contínua puede interpretarse como una administración IV-bolus si se dan las condiciones siguientes:

El desarrollo de McLaurin de la función exponencial es

\[\begin{equation} e^x~=~\sum_{k=0}^\infty~\frac{x^k}{k!} \end{equation}\]

Si \(x \to 0\), \(e^x \approx 1+x\). Por ejemplo, si \(x=0.05\), \(e^{0.05} \approx 1.05127\). Por lo tanto, si \(\lambda_i~t_0\) son suficientemente pequeños, la función \(c_p(t)\) de la fase post-infusión puede interpretarse como la administración IV-bolus a \(t = t_0\):

\[\begin{align} c_p(t) & = k_0~\sum_{i=1}^m \left(1~+~\lambda_i~t_0~-~1 \right)~e^{-\lambda_i~(t~-~t_0)} \\ & = X_0~\sum_{i=1}^m~a_i~e^{-\lambda_i~(t~-~t_0)} \\ & = \sum_{i=1}^m~A_i~e^{-\lambda_i~(t~-~t_0)} \\ \end{align}\]

Parámetros de disposición fundamentales

La tabla 2 resume las ecuaciones para el cálculo de los parámetros básicos de disposición. Las información mostrada por columnas es la siguiente:

• las definiciones matemáticas de los parámetros básicos de disposición; \(E[\theta]\) es el valor esperado del tiempo de residencia.
• las fórmulas para el cálculo a partir de los parámetros de una función m-exponencial estimados a partir de los datos obtenidos tras la administración IV del fármaco (\(X_0\), dosis administrada),
• las fórmulas para el cálculo a partir de los parámetros de la función característica.

 

Tabla 2. Parámetros farmacocinéticos fundamentales de disposición

Definición	\(C_p(t)\)	\(G(t)\)
\(Cl~=~-~\frac{dX/dt}{c_p}\)	\(Cl~=~\frac{X_0}{\sum A_i/\lambda_i}\)	\(Cl~=~\frac{1}{\sum a_i/\lambda_i}\)
\(V_1~=~\frac{X_0}{c_p(0)}\)	\(V_1~=~\frac{X_0}{\sum A_i}\)	\(V_1~=~\frac{1}{\sum a_i}\)
\(V_{AS}~=~\frac{X(t)}{c_p(t)} \Big |_{t \to \infty}\)	\(V_{AS}~=~\frac{X_0}{\lambda_m~\sum A_i/\lambda_i}\)	\(V_{AS}~=~\frac{1}{\lambda_m~\sum~a_i/l_i}\)
\(V_{SS}~=~Cl~E[\theta]\)	\(V_{SS}~=~X_0~\frac{\sum~A_i/\lambda_i^2}{\left(\sum~A_i/\lambda_i \right)^2}\)	\(V_{SS}~=~\frac{\sum~a_i/l_i^2}{\sum~(a_i/l_i)^2}\)
\(k_{11}~=~-~\frac{d~c_p/dt}{c_p}\Big |_{t \to 0}\)	\(k_{11}~=~\frac{\sum~A_i~\lambda_i}{\sum A_i}\)	\(k_{11}~=~\frac{\sum~a_i/l_i}{\sum a_i}\)

 

Ejemplo 1

Considere el modelo de dos compartimentos con los parámetros siguientes: \(V_1 = 0.150 ~ l~kg^{-1}\), \(k_{10} = 0.06 ~ min^{-1}\), \(k12 = 0.100~min^{-1}\), \(k_{21} = 0.080~min^{-1}\). Como es habitual, se supone que el compartimento #1 es el compartimento de entrada y de observación.

En el código mostrado a continuación se muestra el cálculo de la función característica utilizando el método de las matrices:

 

# ----- model parameters
V1 <- 0.150 # l/kg weight
k10 <- 0.06 ; k12 <- 0.100; k21 <- 0.080 # 1/min

# ----- model definition
K <- matrix(c(-(k10+k12), k21, k12, - k21),
            ncol = 2, byrow = TRUE)         # system matrix
x0 <- matrix(c(1,0), ncol = 1)              # initial conditions, dimesionless
B <- matrix(c(1/V1, 0), ncol = 2)           # observation matrix/vector

# ----- characteristic function
M <- eigen(K)
l <- M$values
H <- M$vectors

b <- solve(H) %*% x0
a <- t(H[1,] * b) / V1
rbind(a, l)

##         [,1]        [,2]
##    4.6941610  1.97250570
## l -0.2179796 -0.02202041

La función característica es por lo tanto:

\[\begin{equation} G(t)~=~4.694~e^{-~0.218~t}~+~1.972~e^{-~0.0220~t} \end{equation}\]

Parámetros fundamentales de disposición:

• pk.fchar: vector de parámetros calculados a partir de la función característica (columna \(G(t)\) de la tabla 2).
• pk.mod: vector de parámetros calculados a partir de los parámetros del modelo farmacocinético.

 

pk.fchar <- c(
  Cl = 1/sum(a/-l),
  V1  = 1/(sum(a)),
  Vss = sum(a/l^2)/(sum(a/l))^2,
  Vas = 1/(l[2] * sum(a/l))
)

pk.mod <- c(
  Cl = k10*V1,
  V1 = V1,
  Vss = V1*(1+k12/k21),
  Vas = V1*k10/-l[2]
)

rbind(pk.fchar, pk.mod)

##             Cl   V1    Vss       Vas
## pk.fchar 0.009 0.15 0.3375 0.4087117
## pk.mod   0.009 0.15 0.3375 0.4087117

 

Estimación de los parámetros poblacionales de la función característica

La función nlme::nlme() (nlme) permite estimar los parámetros de modelos estadísticos no lineales para una amplia variedad de diseños experimentales. Para una introducción a la función nlme::nlme() ver Modelos de efectos mixtos con R; para una descripción detallada Pinheiro y Bates (2000).

En los estudios de farmacocinetica de poblaciones más sencillos, orientados al análisis de la variabilidad interindividual con un único nivel de agrupamiento, individuos, y sin incluir covariables, el modelo estadístico en nlme() es el siguiente:

\[\begin{array} \mathbf{y}_i~=~f(\mathbf{\phi}_i~,~\mathbf{x_i})~+~\mathbf{\epsilon}_i \\ \mathbf{\phi}_i~=~\mathbf{A}_i~\mathbf{\beta}~+~\mathbf{B}_i~\mathbf{b_i} \tag{9} \end{array}\]

donde:

• \(\mathbf{y}_i\) es el vector de observaciones de longitud \(n_i\) realizadas sobre el individuo \(i\).
• \(f()\) es la función no lineal del modelo estructural.
• \(\mathbf{\phi}_i\) es el vector de longitud \(p\) de parámetros para el individuo \(i\), función de los efectos fijos (vector \(\mathbf{\beta}\)) y de los efectos aleatorios (vector \(\mathbf{b}_i\))
• \(\mathbf{x_i}\) es el vector, o matriz si el modelo incluye covariables, de predictores para el individuo \(i\).
• \(\mathbf{\epsilon}_i\) es el vector de errores (residuos) para el individuo \(i\); su distribución depende de que el modelo se considere homocedástico o heterocedástico.
• \(\mathbf{A}_i\), matriz determinada por el diseño experimental.
• \(\mathbf{B}_i\) determina los parámetros del modelo que son efectos aleatorios o efectos fijos.

La distribución de los efectos aleatorios es:

 

\[\begin{equation} \mathbf{b}_i~ \sim ~ N(\mathbf{0}~,~\Psi) \tag{10} \end{equation}\]

En farmacocinética poblacional se ha generalizado la terminología introducida por la aplicación nonmem; por lo tanto:

• \(\Psi\), la matriz de covarianza de los efectos aleatorios equivale a matriz equivale a la matriz \(\Omega\) de nonmem
• los términos aleatorios \(b_{j,i}\) se designa por \(\eta_{ji}\) (etas) en nonmem.
• se asume que la distribución de los parámetros farmacocinéticos poblacionaes es log-normal; en `nlme´ es necesario expresar la ecuación del modelo utilizando los logarítmos (neperianos) de los parámetros.

Ejemplo 2

En el ejemplo que se desarrolla más adelante (farmacocinética del cefamandol administrado por perfusión contínua) se concluyó que con los datos disponibles, de los cuatro parámetros de la función característica, \((a_1,~ a_2,~\lambda_1,~\lambda_2)^T\), sólo los tres primeros son efectos aleatorios y que \(\lambda_2\) debe considerarse un efecto fijo. En consecuencia, el vector de parámetros para el individuo \(i\) es:

 

\[\begin{equation} \begin{pmatrix}a_{1,i} \\ a_{2,i}\\ \lambda_{1,i} \\ \lambda_{2,i} \end{pmatrix} ~=~ \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}~ \begin{pmatrix} a_1 \\ a_2 \\ l_1 \\ l_2 \end{pmatrix}~+~ \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{pmatrix} ~ \begin{pmatrix} b_{i,1} \\b_{i,2} \\ b_{i,3} \end{pmatrix} ~=~ \begin{pmatrix}a_1 + b_{1,i} \\ a_2 + b_{2,i} \\ \lambda_1 + b_{3,i} \\ \lambda_2 \end{pmatrix} \end{equation}\]

 

 

Función nlme()

La sintáxis básica es la siguiente:

nlme(model = , data = , groups = , fixed = , random = , weights = , start =)

model

Expresión R del modelo farmacocinético.

data y groups

El diseño experimental utilizado puede implementarse de dos formas:

• Utilizando como tabla de datos un objeto tipo groupedData que incluye como metadato la estructura de los datos. Para constrruir una tabla groupedData() a partir de un data.frame se utiliza la función del mismo nombre.
• Utilizar una tabla tipo data.frame y especificar la estructura de los datos utilizando el argumento groups.

fixed

Fórmula que especifica los parámetros fijos del modelo; puede expresarse de dos formas:

• fixed = p1 + p2 + ... ~ 1, donde p1, etc., son los parámetros del modelo,
• fixed = list(p1 ~ 1, p2 ~ 1, ...).

random

Especifica los términos aleatorios del modelo y la estructura de la matriz \(\Omega\). En los modelos farmacocinéticos poblacionales se asume habitualmente que la matriz \(\Omega\) (\(\Psi\)) es diagonal (los efectos aleatorios son independientes), en cuyo caso se utilizará la función pdDiag(). Ver ejemplo cefamandol más adelante.

weights

Por defecto se asume que el modelo es homocedástico; en caso contrario la heterocedasticidad puede modelarse utilizando funciones de la familia variance function classess. Ver un resumen aquí.

start

Un vector con los valores iniciales.

Ejemplo 3: cefamandol

La librería nlme incluye los datos de concentración plasmática - tiempo de cefamandole administrado por perfusión contínua (\(k_0~=~1.5~mg~kg^{-1}~min^{-1}\), \(t_0~=~10~min\)) a seis individuos (Davidian y Giltinan, 1993); se tomaron 14 muestras de plasma entre los 10 minutos (momento en que finaliza la perfusión) y los 360 minutos. Estos datos también se utilizan en el ejercicio propuesto 2 del capítulo Nonlinear mixed-effects models: basic concepts and motivating examples del libro de Pinheiro y Bates (2000), aunque no se da ninguna solución en pafrticular.

En el código siguiente se importa la tabla Cefamandole y se muestran los tres primeros registros. Observar que la tabla incluye como metadato la fórmula de agrupamiento de los datos, que se muestran en la figura siguiente.

data("Cefamandole", package = "nlme")
head(Cefamandole, 3)

## Grouped Data: conc ~ Time | Subject
##   Subject Time  conc
## 1       1   10 127.0
## 2       1   15  80.0
## 3       1   20  47.4

colrs <- c("black", "red", "blue", "darkgreen", "orange", "violet")
plot(NULL, xlim = c(0, 400), ylim = c(0.1, 250),
     log = "y", xlab = "tiempo (min)", ylab = "cp (mg/l)")
for(i in 1:6) lines(conc ~ Time, data = subset(Cefamandole, Subject == i),
                    lty = 2, type = "b", col = colrs[i], pch = 16)
legend("topright", legend = paste0("S",1:6), col = colrs,
       lty = 1, bty = "n", cex = 0.75, pch = 16)

Estimación de los parámetos del modelo

El desarrollo de un modelo farmacocinético experimental implica resolver los puntos siguientes:

• reparametrización e identificación del modelo estructural,
• identificar los parámetros efectos aleatorios y los parámetros efectos fijos,
• en el caso de heterodasticidad, seleccionar el modelo adecuado.

Para abreviar la exposición aquí nos limitaremos a comparar dos modelos. Las características comunes a ambos modelos son:

• función bi-exponencial como modelo estructural de la función característica,
• se ha reparametrizado el modelo para utilizar los logaritmos neperianos de los parámetros; la razón es que habitualmente se observa que los parámetros farmacocinéticos poblacionales tienen una distribución log-normal (ver Apéndice A más adelante).
• Los parámetros \(a_1\), \(a_2\) y \(\lambda_1\) son aleatorios; el parámetro \(\lambda_2\) es fijo.

La diferencia entre ambos modelos es:

• Modelo homocedástico: los errores (\(\mathbf{\epsilon}_i\) en la ecuación 9) se suponen independientes con distribución normal e igual varianza, \(\mathbf{\epsilon}_i \sim N(\mathbf{0}~,~\mathbf{I}_{n_i \times n_i})\).
• Modelo heterocedástico: la varianza de la respuesta se interpreta utilizando el modelo \(\sigma_y^2~=~\sigma^2~(\delta_1 + |x_{ij}|^{\delta_2})^2\) utilizando como argumento opcional respecto al modelo homocedástico weights = varConstPower(power = 0.2).

En el codigo mostrado a continuación:

• se estiman los parámetros utilizando el modelo homocedástico (cefA.nlme) y heterocedástico (cefB.nlme),
• se comparan ambos modelos mediante la función anova().

cefA.nlme <- nlme(conc ~ 1.5*((exp(la1)/exp(ll1))*(exp(exp(ll1)*10) - 1)*exp(-exp(ll1)*Time) +
                    (exp(la2)/exp(ll2))*(exp(exp(ll2)*10) - 1)*exp(-exp(ll2)*Time)),
                  data = Cefamandole,
                  fixed = la1 + la2 + ll1 + ll2 ~ 1,
                  random = pdDiag(la1 + la2 + ll1 ~ 1),
                  start = c(la1 = log(8), la1 = log(8), ll1 = log(0.1), ll2 = log(0.02)),
                  )

cefB.nlme <- update(cefA.nlme,
                    weights = varConstPower(power = 0.2))
anova(cefA.nlme, cefB.nlme)

##           Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## cefA.nlme     1  8 538.8141 558.2607 -261.4071                        
## cefB.nlme     2 10 436.5305 460.8387 -208.2653 1 vs 2 106.2836  <.0001

Análisis gráfico

Las gráficas mostradas a continuación:

• fila superior, residuos estandarizados,
• fila inferior: concentración plasmática de cefamandol en función de las concentraciones predichas a nivel poblacional y a nivel individual para el modelo cefB.nlme.

fig1a <- plot(cefA.nlme, ylim = c(-4, 4), main = "cefA.nlm2")
fig1b <- plot(cefB.nlme, ylim = c(-4, 4), main = "cefB.nlme")
fig2a <- plot(cefB.nlme, conc ~ fitted(.,0), xlim = c(0, 250), ylim = c(0, 250),
     abline = c(0, 1), main = "cefB.nlme - población")
fig2b <- plot(cefB.nlme, conc ~ fitted(.,1), xlim = c(0, 250), ylim = c(0, 250),
     abline = c(0, 1), main = "cefB.nlme - individual")
grid.arrange(fig1a, fig1b, fig2a, fig2b, nrow = 2, ncol = 2)

Resultados

El código siguiente muestra el resumen de los resultados del modelo seleccionado, cefB.nlme:

cefB.nlme

## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: conc ~ 1.5 * ((exp(la1)/exp(ll1)) * (exp(exp(ll1) * 10) - 1) *      exp(-exp(ll1) * Time) + (exp(la2)/exp(ll2)) * (exp(exp(ll2) *      10) - 1) * exp(-exp(ll2) * Time)) 
##   Data: Cefamandole 
##   Log-likelihood: -208.2653
##   Fixed: la1 + la2 + ll1 + ll2 ~ 1 
##        la1        la2        ll1        ll2 
##  2.1360194  0.1945783 -3.3083801 -4.6074789 
## 
## Random effects:
##  Formula: list(la1 ~ 1, la2 ~ 1, ll1 ~ 1)
##  Level: Subject
##  Structure: Diagonal
##             la1       la2        ll1 Residual
## StdDev: 0.25533 0.6427553 0.04691619 0.164903
## 
## Variance function:
##  Structure: Constant plus power of variance covariate
##  Formula: ~fitted(.) 
##  Parameter estimates:
##        const        power 
## 5.391697e-08 1.023145e+00 
## Number of Observations: 84
## Number of Groups: 6

Matriz de correlación efectos fijos

Una elevada correlación entre las estimadas del modelo nos indica bien una sobreparametrización o bien un diseño experimental inadecuado para estimar los parámetros del modelo. En este caso los coeficientes de correlación entre estimadas son moderadamente bajos:

print(cov2cor(vcov(cefB.nlme)), digits = 4)

##          la1      la2    ll1      ll2
## la1  1.00000 -0.01939 0.1781 -0.05664
## la2 -0.01939  1.00000 0.3212  0.43699
## ll1  0.17810  0.32122 1.0000  0.64296
## ll2 -0.05664  0.43699 0.6430  1.00000

Efectos fijos

Estimadas del los factores fijos del modelo reparametrizado (logaritmos neperianos de los parámetros de la función característica); la tabla siguiente muestra las estimadas de los efectos fijos (primera fila) y sus anti-logaritmos (\(a_1\) y \(a_2\) en \(mg~l^{-1}\); \(\lambda_1\) y \(\lambda_2\) en \(min^{-1}\)).

print(rbind(fixed.effects(cefB.nlme),
      exp(fixed.effects(cefB.nlme))), digits = 4)

##        la1    la2      ll1       ll2
## [1,] 2.136 0.1946 -3.30838 -4.607479
## [2,] 8.466 1.2148  0.03658  0.009977

De acuerdo con la teoría de la estimación de modelos no lineales, las estimadas tienen una distribución asintótica normal multivariante; por lo tanto, las antitransformadas tendrían en su caso una distribución log-normal multivariante; además, las antitransformadas no son los valores esperados (medias) sino las modas (ver Apéndice A más adelante).

Matriz omega

Matriz \(\Omega\) (en la terminología habitual en farmacocinética, equivalente a la matriz \(\Psi\) en nlme), expresada utilizando las desviaciones típicas, extraidas de la sección Random effects del objeto cefB.nlme es:

\[\begin{equation} \mathbf{\Omega} ~=~ \begin{pmatrix}0.2553^2 & 0 & 0 \\ 0 & 0.6427^2 0 \\ 0 & 0 & 0.0469^2 \end{pmatrix} \end{equation}\]

Obsérvese que la matriz \(\Omega\) es de dimensión \(3 \times 3\) porque el parámetro \(\lambda_2\) se considera fijo.

Predicciones

plot(augPred(cefB.nlme, level = 0:1), log = "y",length.out = 2)

Tiempo de perfusión y selección de la función de entrada.

¿El tiempo de perfusión utilizado, \(t_0 = 10\) min, es suficientemente corto como para interpretar la administración como si hubiera sido IV-bolus? El valor de \(\lambda_1\) y de \(t_0 \times \lambda_1\) son

pk.est <- exp(fixed.effects(cefB.nlme))
l1 <- unname(pk.est["ll1"])
c( l1 = l1, "t0*l1" = 10*l1)

##         l1      t0*l1 
## 0.03657538 0.36575375

\(t_0 \times \lambda_1 > 0.05\), por lo que no es aceptable interpretar la administración como IV-bolus.

Parámetros de disposición

Hay tres familias de métodos aplicables a la estimación de la variabilidad inter-individual de los parámetros fundamentales de disposición a partir de los parámetros poblacionales de la función característica: el método delta (Hoef, 2012), el método Monte Carlo de simulación estadística y el método bootstrap; los dos últimos forman parte del conjunto de métodos de estadística computacional. La bibliografía sobre estos temas es muy amplia; para más detalles ver por ejemplo Rizzo (2019).

El método bootstrap tiene menos restricciones respecto a la distribución estadística de los datos originales, por lo que se utiliza ferecuentemente en los estudios de farmacocinética para verificar la distribución de los parámetros poblacionales. Para más detalles ver Thai y cols. (2014). en este tutorial se utiliza la función nlraa::boot_nlme() (Miguez, 2023) que tiene por argumento un objeto nlme.

Ejemplo 3: continuación

En el código siguiente se implementa el análisis gráfico y la estimación del intervalo de confianza no parámetrico de los parámetros de disposición fundamentales del cefamandol.

Detalles:

• La tabla T5 contiene 999 muestras bootstrap generadas a partir del resultado cefB.nlme.
• Es posible que algunas de las combinaciones de individuos y residuos dentro de individuos den lugar a muestras bootstraap para las cuales la estimación de los parámetros no convergen (la función boot_nmle() emite un mensaje indicando el número de casos). Los casos en los que la convergencia ha sido exitosa se identifican en el vector lógico srow.
• La tabla pkDisp.boot contiene los valores anti-transformados logarítmicamente de los parámetros bootstrap.
• Se completa la tabla pkDisp.boot con los parámetros de disposición utilizando las fórmulas recogidas en la tabla 2.
• Los valores obtenidos de los parámetros de disposición se presentan gráficamente.
• El tiempo de ejecución es considerablemente elevado (puede llevar algunos minutos). Es recomendable utilizar el argumento cores = para activar la ejecución en paralelo de la función boot().

n.boot <- 999
t.ini <- Sys.time()
T5 <- boot_nlme(cefB.nlme, R = n.boot, cores = 4)

## Number of times model fit did not converge 80 out of 999

t.end <- Sys.time()

srow <- complete.cases(T5$t)
pkDisp.boot <- data.frame(
  a1 = exp(T5$t[srow, 1]),
  a2 = exp(T5$t[srow, 2]),
  l1 = exp(T5$t[srow, 3]),
  l2 = exp(T5$t[srow, 4])
)

pkDisp.boot <- within(pkDisp.boot, 
                     {Vas <- 1/(l2 * (a1/l1 + a2/l2))
                      Vss <- (a1/l1^2 + a2/l2^2)/((a1/l1 + a2/l2)^2)
                      V1 <- 1/(a1 + a2)
                      Cl <- 1/(a1/l1 + a2/l2)
                      })

panel.hist <- function(x, ...)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(usr[1:2], 0, 1.5) )
    h <- hist(x, plot = FALSE)
    breaks <- h$breaks; nB <- length(breaks)
    y <- h$density; y <- y/max(y)
    rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...)
}
pairs(pkDisp.boot[, 5:8], gap = 0.2, panel = panel.smooth,
      cex = 1, bg = "light blue", horOdd=TRUE,
      diag.panel = panel.hist, cex.labels = 1.5, font.labels = 2)

Tiempo de ejecución del código (minutos): 1.52.

Intervalo de confianza 95% no parámetrico

La figura anterior muestra la distribución de los 4 parámetros de disposición esperada para la población a partir de las muestras generadas por el método bootstrap.

Además:

• La elevada correlación entre los 4 parámetros de disposición. Podría entenderse la correlación entre los 3 volúmenes aparentes de distribución (\(V_1\), \(V_{AS}\) y \(V_{ES}\)), pero no entre el aclaramiento plasmático (\(Cl\)) y los volúmenes de distribución. Esta correlación es debida a que los 4 parámetros de distribución son función de los cuatro parámetros de la funcion característica (\(a_1\), \(a_2\), \(\lambda_1\), \(\lambda_2\)).
• La distribución de parámetros es aproximadamente normal. Los valores de los percentiles \(2.5\%\), \(50\%\) y $97.5%& son:

S <- apply(pkDisp.boot[, 5:8], 2, quantile, probs = c(0.025, 0.50, 0.975))
print(S, digits = 4)

##             Cl      V1    Vss    Vas
## 2.5%  0.002181 0.08082 0.1155 0.2157
## 50%   0.002797 0.10212 0.1469 0.2807
## 97.5% 0.003575 0.13094 0.1845 0.3768

La amplitud del intervalo de confianza no parámetrico respecto al percentil 50% y expresado en tanto por ciento, es una medida de la variabilidad inter-individual de los parámetros de disposición:

100*(S[3,]-S[1,])/S[2,]

##       Cl       V1      Vss      Vas 
## 49.84857 49.07337 46.98539 57.40917

Comentarios finales

El número de individuos incluidos en la tabla Cefamandole, 6, es muy pequeño para extraer conclusiones sobre la farmacocinética poblacional de este fármaco por dos razones:

• Muestra pequeña en la etapa de estimación de los parámetros poblacionales de la función característica.
• Muestra igualmente insuficiente en la estimación bootstrap de los parámetros de disposición. Obsérvese que se generaron 999 muestras bootstral (argumento R en la función boot_nlme()), mientras que el número de combinaciones con repetición de los seis individuos (\(n\)) tomados de 6 en 6 (\(r\)) es:

\[\begin{equation} \begin{pmatrix} n~+~r~-~1 \\ r \end{pmatrix} ~=~ \begin{pmatrix}11 \\ 6 \end{pmatrix}~=~462 \end{equation}\]

No obstante, los datos utilizados es un ejemplo sencillo de como implementar la metodología descrita.

Bibliografía

Davidian, M., Giltinan, D.M. (1993) Analysis of repeated measurement data using the nonlinear mixed effects models Chemometrics and Intelligent Laboratory Systems 20 1 - 24. https://doi.org/10.1016/0169-7439(93)80017-C

DiStefano, J.J., Landaw, M. (1984) Multiexponential, multicompartmental, and noncompartmental modeling. I. Methodological limitations and physiological interpretations. American Journal of Physiology 246 R651 - R654 PubMed

DiStefano, J.J., Landaw, M. (1984b) Multiexponential, multicompartmental, and noncompartmental modeling. II. Data analysis and statistical considerations. American Journal of Physiology 246 R665 - R677 PubMed

Gillespie, W.R. (1997) Convolution . based approach for in vivo - in vitro correlation modeling. In: Young, D., Devane, J.G., Butler, J. (eds) In vitro . in vivo correlations. Advances in experimental medicine and biology, vol 423, Springer, Boston, MA https://doi.org/10.1007/978-1-4684-6036-0_5

Gomeni, R., Bressolle-Gomeni, F. (2021) Modeling complex pharmacokinetics of long-acting injectable products using convolution-based models with nonparametric imput functions. The Journal of Clinical Pharmacology 61 1081 - 1095 doi: 10.1002/jcph.1842

Hoef, J.M. (2012) Who invented the delta method? The American Statistician 66 124 - 127 ResearchGate

Miguez, F. (2023) Nonlinear regression for agricultural applications CRAN.R

Pinheiro, J.C., Bates, D.M. (2000) Mixed-Effects Models in S and S-PLUS. Springer-Verlag

Rizzo, M.L. (2019) Statistical Computing with R, 2nd ed. CRC Press ISBN-10:1466553324 Thai, H.-T., Mentré, F., Holford, N.H.G., Veyrat-Follet, C., Comets, E. (2014) Evaluation of bootstrap methods for estimating uncertainty of parameters in nonlinear-mixed effects models: a simulation study in population pharmacokinetics Journal of Pharmacokinetics and Pharmacodynamics 41 15 - 33

Thai, N.-T., Mentré, F., Holdford, N.H.G., Veyrat-Follet, C., Comets, E. (2014) Evaluation of bootstrap methods for estimating uncertainty of parameters in nonlinear mixed.effects models: a simulation study in population poharmacokinetics Journal of Pharmacokinetics and Pharmacodynamics 41 15 -33 doi: 10.1007/510928-013-9343-z

van Rossum, J.M., de Bie, J.E.G.M. (1989) System dynamics in clinical pharmacokinetics. An introduction Clinical Pharmacokinetics 17 27 - 44 https://doi.org/10.2165/00003088-198917010-00003

Wagner, J.G. (1976) Linear pharmacokinetic equations allowing direct calculation of many needed pharmacokinetic parameters from the coefficients and exponents of poliexponential equation which have been fitted to the data Journal of Pharmacokinetics and Biopharmaceutics 4 443-467 https://doi.org/10.1007/BF01062831

Wagner, J. G. and coworkers (1977) Pharmacokinetic parameters estimated from intravenous data by uniform methods and some of their uses Journal of Pharmacokinetics and Biopharmaceutics 5 161 - 182 https://doi.org/10.1007/BF01066219

Apéndice

Distribución logarítmico - normal

Sea \(Y\) una variable aleatoria con distibución normal con media \(\mu\) y varianza \(\sigma^2\), lo que se expresa como \(Y~\sim~N(\mu~,~\sigma^2)\), y \(X~=~e^Y\). Se dice entonces que la variable aleatoria \(X\) tiene una distribución log-normal (también denominada logarítmco - normal y normal - logarítmica). La tabla A1 recoge las respectivas funciones de densidad, valores esperados (\(E[X]\)) y varianzas (\(var[X]\)).

Tabla A1. Funciones de densidad, valores esperados y varianzas de la distribución normal y logarítmico normal

distribución	función densidad	\(E[X]\)	\(~~moda~~\)	\(VAR[X]\)
normal	\(f(x)~=~\frac{1}{\sigma~\sqrt{2~\pi}}~exp \left(-~\frac{(x~-~\mu)^2}{2~\sigma^2} \right)\)	\(\mu\)	\(\mu\)	\(\sigma^2\) \(\sigma > 0\)
log - normal	\(f(x)~=~\frac{1}{x~\sigma~\sqrt{2~\pi}}~exp \left(-~\frac{(ln(x)~-~\mu)^2}{2~\sigma^2}\right)\)	\(e^{\mu+\sigma^2/2}\)	\(e^\mu\)	\(\left(e^{\sigma^2}~-~1 \right)~e^{2~\mu~+~\sigma^2}\) \(\sigma > 1\)

 

par(mfrow = c(1,2))
my <- 2
sy <- 0.2
mx <- exp(my + sy^2/2)
sx <- sqrt((exp(sy^2) - 1)*exp(2*my + sy^2))
plot(function(x) dnorm(x, my, sy), from = my-3*sy, to = my+3*sy,
     ylab = "densidad", ylim = c(0, 2.5), xlab = "Y",
     main = "distribución normal", cex.main = 1)
text(2.4, 2.3,
     labels = substitute(paste(mu,' = ', val, sep = ""), list(val=my)))
text(2.4, 2.1,
     labels = substitute(paste(sigma,' = ', val, sep = ""), list(val=sy)))

plot(function(x) dlnorm(x, log(mx), log(sx)), from = mx/(sx^3), to = mx*sx^3,
     ylab = "densidad", xlab = "X", ylim = c(0, 0.2),
     main = "distribución log-normal", cex.main = 1)
text(20, 0.18,
     labels = substitute(paste(mu,' = ', val, sep = ""), list(val=round(mx, 3))))
text(20, 0.16,
     labels = substitute(paste(sigma,' = ', val, sep = ""), list(val=round(sx, 3))))

Figura A1

# más información sobre expresiones matemáticas en gráficos: ?plotmath

Coeficiente de variación

\[\begin{equation} CV(X)~=~\frac{\left( \left(e^{\sigma^2}~-~1 \right)~e^{2~\mu~+~\sigma^2}\right)^{1/2}}{e^{\mu+\sigma^2/2}}~=~\sqrt{e^{\sigma^2}~-~1} \end{equation}\]

Si \(\sigma\) es pequeño, \(e^{\sigma^2}\) puede aproximarse por los dos primeros términos de la serie de McLaurin, \(e^{\sigma^2}~\approx~1~+~\sigma^2\), y por lo tanto \(CV(X)~\approx~\sigma\).