Introducción

Los modelos farmacocinéticos lineales están definidos por un sistema de n ecuaciones diferenciales lineales de primer grado con coeficientes constantes que expresan el balance de materia diferencial en cada compartimento. Para la administración IV-bolus, el sistema puede formularse como un sistema homogéneo, y bajo ciertas condiciones la solución es un conjunto n funciones n-exponenciales,

\[\begin{equation} x_i(t)~=~\sum_{j=1}^n~A_{ij}~e^{-\lambda_j~t} \tag{1} \end{equation}\]

donde xi(t) es la cantidad de fármaco en el compartimento i a tiempo t. Para el modelo de 2 compartimentos, las expresiones de los exponentes (λj) y de los coeficientes (Aij) son bien conocidas y pueden encontrarse en los textos de farmacocinética. Para 3 o más compartinentos la complejidad de las ecuaciones aumenta con el número de compartimentos debido a:

El método de las matrices permite calcular directamente las soluciones de los sistemas lineales homogéneos a partir del sistema de ecuaciones diferenciales; además, su implementación es muy sencilla si se dispone de una función o una subrutina para el cálculo de los autovalores y autovectores de una matriz.

En este documento se incluyen:

Modelos compartimentales lineales

El sistema de ecuaciones diferenciales que definen el modelo se deriva a partir del esquema del modelo; en este se identifican (ver figura 1)

 

Figura 1. Esquemas de dos planteamientos de los modelos compartimentales lineales

Figura 1. Esquemas de dos planteamientos de los modelos compartimentales lineales

 

En farmacocinética encontramos tres formas del balance de materia (ecuaciones 2 - 6):

En el análisis compartimental la variable de estado es la cantidad de substancia o la masa que forma cada compartimento. El balance de materia para el compartimento i es:

\[\begin{equation} \frac{d~x_i}{dt}~=~- k_{ii} x_i~+~\sum_{j=1 , j \ne i}^n~k_{ji}~x_j~+~f_i(t) \tag{2} \end{equation}\]

donde xi es la cantidad de substancia que forma el compartimento i; kij (dimensiones t-1) es la constante de velocidad de transporte o transformación desde el compartimento i al compartimento j; kii es la suma de las constantes de velocidad desde el compartimento i hacia los demás compartimentos, incluido el exterior del sistema identificado con el número 0:

\[\begin{equation} k_{ii}~=~\sum_{j=0 , j \ne i}^n~k_{ij} \tag{3} \end{equation}\]

fi(t) (dimensiones m·t-1) es la velocidad de entrada de fármaco en el sistema a través del compartimento i.

La magnitud obervada no suele ser la cantidad de substancia sino la concentración plasmática o, en el caso de radiofármacos, la actividad, razón por la que se incluye en el modelo el volumen aparente de distribución. En casos más complejos la magnitud observada es una combinación lineal de las soluciones del modelo (ver el Caso 2). Generalizando,

\[\begin{equation} \mathbf{y}(t)~=~\mathbf{B}·\mathbf{x}(t) \tag{4} \end{equation}\]

donde y(t) es el vector observado, x(t) la solución del sistema, y B la matriz de observación.

El número de filas de la matriz de observación es igual al números de magnitudes observadas, y el número de columnas igual al número de compartimentos del modelo. Si sólo se observa una magnitud (por ejemplo, la concentración plasmática), la matriz de observación se reduce a un vector fila.

En ocasiones resulta más conveniente utilizar directamente como variable de estado la concentración de fármaco en cada compartimento; de la ecuación 2 se deduce,

\[\begin{equation} V_i~\frac{d~c_i}{dt}~=~- k_{ii}~V_i~c_i~+~\sum_{j=1 , j \ne i}^n k_{ji}~V_j~c_j~+~f_i(t) \tag{5} \end{equation}\]

donde Vi es el volumen aparente de distribución del compartimento i.

En los estudios de farmacocinética de poblaciones es preferible utilizar como parámetros del modelo el aclaramiento, los volúmenes de distribución y los aclaramientos inter - compartimentales, y el balance de materia es,

\[\begin{equation} V_i~\frac{d~c_i}{dt}~=~-\left(Cl_i + \sum_{j=0}^n Q_{ij}\right)~c_i~+~\sum_{j=1,j \ne i}~Q_{ji}~c_j + f_i(t) \tag{6} \end{equation}\]

donde Cli el aclaramiento del fármaco a partir del compartimento i, y Qij el aclaramiento inter-compartimental desde el compartimento i al compartimento j. Observar que Qij = Qji. La ventaja de esta formulación es se utiliza el aclaramiento del fármaco como parámetro primario, que a su vez puede expresarse como función de las covariables del modelo (peso del individuo, sexo, etc.).

Los tres planteamientos son equivalentes desde el punto de vista matemático.

Para un modelo compartimental, el sistema de ecuaciones diferenciales puede escribirse como:

\[\begin{equation} \frac{d}{dt} \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{pmatrix} ~=~ \begin{pmatrix} -k_{11} & k_{21} & \cdots & k_{n1} \\ k_{12} & - k_{22} & \cdots & k_{n2} \\ \cdots & \cdots & \cdots & \cdots \\ k_{1n} & k_{2n} & \cdots & -k_{nn} \end{pmatrix} ~ \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{pmatrix} ~+~ \begin{pmatrix} f_1(t) \\ f_2(t) \\ \vdots \\f_n(t) \end{pmatrix} \tag{7} \end{equation}\]

Y utilizando la notación matricial,

\[\begin{equation} \bf{\dot{x}}~=~\bf{K}~\bf{x}+ \bf{f}(t) \tag{8} \end{equation}\]

La matriz K recibe el nombre de matriz del sistema.

Sistemas homogéneos

En los sistema homogéneos no existe función de entrada y por lo tanto la ecuación 4 se reduce a,

\[\begin{equation} \bf{\dot{x}}~=~\bf{K}~\bf{x} \tag{9} \end{equation}\]

Puede probarse que la solución para las condiciones iniciales x = x0 para t = 0 es

\[\begin{equation} \bf{x}(t)~=~e^{\bf{K}~t}~\bf{x_0} \tag{10} \end{equation}\]

El término exponencial tiene una interpretación similar a la expansión de la función exponencial,

\[\begin{equation} e^x~=~1 + x + \frac{x^2}{2!} + \frac{x^3}{3!}~=~\sum_{i=0}^\infty~\frac{x^i}{i!} \end{equation}\]

\[\begin{equation} e^{\mathbf{K}~t}~=~\mathbf{I} + \mathbf{K}~t + \frac{\mathbf{K}^2~t^2}{2!} + \dots + \frac{\mathbf{K}^n~t^2}{n!} + \dots \end{equation}\]

donde I es la matriz unidad. Igualmente puede probarse que

\[\begin{equation} e^{\bf{K}~t}~=~\bf{H}~e^{\bf{\lambda~I}~t}~\bf{H}^{-1} \tag{11} \end{equation}\]

λ es el vector de autovalores (valores propios) de la matriz K, es decir, las soluciones del polinomio característico de la matriz K,

\[\begin{equation} |\mathbf{K}~-~\lambda~\mathbf{I} | ~=~0 \tag{12} \end{equation}\]

y H la matriz formada por los correspondientes autovalores. El vector propio hi asociado al autovalor λi es aquel que satisface la ecuación,

\[\begin{equation} \left(\mathbf{K}~-~\lambda_i~\mathbf{I} \right)~\mathbf{h_i}~=~0 \tag{13} \end{equation}\]

Combinando las ecuaciones 10 y 11 obtenemos finalmente,

\[\begin{equation} \bf{x}(t)~=~\bf{H}~e^{\bf{\lambda~I}~t}~\bf{H}^{-1}~\bf{x_0} \tag{14} \end{equation}\]

Pueden darse tres casos:

En los modelos farmacocinéticos de uso habitual los autovalores de la matriz K son reales negativos y distintos; para ello deben cumplirse dos condiciones:

\[\begin{equation} k_{11} < k_{22} < \dots < k_{nn} \end{equation}\]

Ejemplo

En este ejemplo se muestra paso a paso la aplicación del método de las matrices para resolver el modelo de dos compartimentos mostrado en la figura 1-A. Los parámetros del modelo son: k10 = 0.1 h-1, k12 = 0.2 h-1, k20 = 0.1 h-1, k21 = 0.1 h-1, x0 = 100 mg y V = 20 l. Como sólo se observa el compartimento 1 la matriz de observación es un vector:

\[\begin{equation} \mathbf{B}~=~ \begin{pmatrix} 1/20 & 0 \end{pmatrix} \end{equation}\]

La matriz del sistema es:

\[\begin{equation} \mathbf{K}~=~ \begin{pmatrix} -(0.1+0.2) & 0.1 \\ 0.2 & -(0.1+0.1) \end{pmatrix} ~=~ \begin{pmatrix} -0.3 & 0.1 \\ 0.2 & -0.2 \end{pmatrix} \end{equation}\]

Observar que se dan las dos condiciones necesarias y suficientes para que los autovalores de la matriz sean reales negativos. El polinomio característico de la matriz del sistema es (ecuación 8):

\[\begin{align} \begin{vmatrix} -0.3 - \lambda & 0.1 \\ 0.2 & -0.2 - \lambda \end{vmatrix} ~=~(-0.3 - \lambda)~(-0.2 - \lambda) - 0.1\times 0.2~=~ \\ \lambda^2+0.5~\lambda+0.04~=~0 \end{align}\]

Las soluciones de esta ecuación de segundo grado son λ1 = - 0.4 y λ2 = - 0.1 h-1. El vector propio asociado al primer autovalor es (ecuación 13)

\[\begin{equation} \left[ \begin{pmatrix} -0.3 & 0.1 \\ 0.2 & -0.2 \end{pmatrix} ~-~ \begin{pmatrix} -0.4 & 0 \\ 0 & -0.4 \end{pmatrix} \right] ~ \begin{pmatrix} h_{11} \\ h_{12} \end{pmatrix} ~=~0 \end{equation}\]

de donde resulta el sistema lineal,

\[\begin{align} 0.1~h_{11}+0.1~h_{12}~=~0 \\ 0.2~h_{11}+0.2~h_{12}~=~0 \end{align}\]

Utilizando el método de reducción gaussiana se obtiene la ecuación del sub-espacio:

\[\begin{equation} h_{12}~=~-h_{11} \end{equation}\]

Para calcular la base se da un valor arbitrario a uno de los elementos del vector h1, por ejemplo, h11 = 1 y se calcula el segundo, resultando,

\[\begin{equation} \mathbf{h_1}~=~\begin{pmatrix} 1 \\ -1 \end{pmatrix} \end{equation}\]

El sistema lineal para el segundo autovalor λ2 = - 0.1 h-1 es,

\[\begin{align} -0.2~h_{21}+0.1~h_{22}~=~0 \\ 0.2~h_{21}-0.1~h_{22}~=~0 \end{align}\]

Multiplicando la primer ecuación por -1 y sumando ambas ecuaciones se obtiene la ecuación del subespacio,

\[\begin{equation} h_{22}~=~\frac{0.4}{0.2}~h_{21} \end{equation}\]

y haciendo h21 = 1 se obtiene el segundo vector propio:

\[\begin{equation} \mathbf{h_2}~=~\begin{pmatrix} 1 \\ 2 \end{pmatrix} \end{equation}\]

La matriz H y su inversa son:

\[\begin{equation} \mathbf{H}~=~\begin{pmatrix} 1 & 1 \\ -1 & 2 \end{pmatrix} \end{equation}\]

\[\begin{equation} \mathbf{H}^{-1}~=~\begin{pmatrix} \tfrac{2}{3} & -\tfrac{1}{3} \\ \tfrac{1}{3} & \tfrac{1}{3} \end{pmatrix} \end{equation}\]

Aplicando la ecuación 14,

\[\begin{align} \begin{pmatrix} x_1(t) \\ x_2(t) \end{pmatrix} ~=~ \begin{pmatrix} 1 & 1 \\ -1 & 2 \end{pmatrix} ~ \begin{pmatrix} e^{-0.4~t} & 0 \\ 0 & e^{-0.1~t} \end{pmatrix} ~ \begin{pmatrix} \tfrac{2}{3} & -\tfrac{1}{3} \\ \tfrac{1}{3} & \tfrac{1}{3} \end{pmatrix} ~ \begin{pmatrix} 100 \\ 0 \end{pmatrix} ~=~\\ \begin{pmatrix} 1 & 1 \\ -1 & 2 \end{pmatrix} ~ \begin{pmatrix} e^{-0.4~t} & 0 \\ 0 & e^{-0.1~t} \end{pmatrix} ~ \begin{pmatrix} \tfrac{2}{3}~100 \\ \tfrac{1}{3}~100 \end{pmatrix}~=~ \\ \begin{pmatrix} 100~\left(\frac{2}{3}~e^{-0.4~t} + \frac{1}{3}~e^{-0.1~t} \right) \\ 100~\left(- \frac{2}{3}~e^{-0.4~t} + \frac{2}{3}~e^{-0.1~t} \right) \end{pmatrix} \end{align}\]

Multiplicando la solución por la matriz de observación obtenemos la solución del ejemplo:

\[\begin{equation} c_p(t)~=~\frac{10}{3}~e^{-0.4~t} + \frac{5}{3}~e^{-0.1~t} \end{equation}\]

Normalización autovalores

Aunque no es imprescindible, es habitual normalizar los autovectores hacendo que su módulo sea igual a 1. Por ejemplo, para el primer autovector h1,

\[\begin{equation} \frac{\begin{pmatrix} 1 \\ -1 \end{pmatrix}}{\left( \begin{pmatrix}1 & -1 \end{pmatrix}~\begin{pmatrix} 1 \\ -1 \end{pmatrix} \right)^{1/2} } ~=~\begin{pmatrix}0.7072 \\ - 0.7072 \end{pmatrix} \end{equation}\]

Implementación del método en R

La función R eigen():

Función eigen()

La sintáxis de la función es:

eigen(A, symmetric, ...)

Argumentos:

  • A es una matriz real o compleja,
  • symmetric es una variable lógica para indicar que la matríz es simétrica (symmetric = TRUE) o no simétrica (symmetric = FALSE). Cuando no se incluye este argumento la función eigen() identifica el tipo de matriz utilizando la función is.symmetric().

La función eigen() devuelve un objeto tipo list con dos elementos: $values (los autovalores de la matriz) y $vectors (la matriz de autovectores).

En el código siguiente se calcula los autovalores y autovectores de la matriz del sistema del ejemplo anterior,

eigen(
  matrix(c(-0.3, 0.2, 0.1, -0.2), ncol = 2)
)
## eigen() decomposition
## $values
## [1] -0.4 -0.1
## 
## $vectors
##            [,1]       [,2]
## [1,] -0.7071068 -0.4472136
## [2,]  0.7071068 -0.8944272

Ejemplo (continuación)

Los pasos para resolver el ejemplo utilizando la función eigen() son:

Primero, definir la matriz del sistema, el vector de condiciones iniciales, la matriz de observación, y calcular los autovalores y autovectores de la matriz del sistema:

K <- matrix(c(-0.3, 0.2, 0.1, -0.2), ncol = 2)          # matriz del sistems
x0 <- matrix(c(100, 0), ncol = 1)                       # vector condiciones iniciales
B <- matrix(c(1/20, 0), ncol = 2)                       # matriz de observación
M <- eigen(K)                                           # cálculo autovalores y autovectores
M
## eigen() decomposition
## $values
## [1] -0.4 -0.1
## 
## $vectors
##            [,1]       [,2]
## [1,] -0.7071068 -0.4472136
## [2,]  0.7071068 -0.8944272

Segundo, extraer la matriz de autovectores y calcular b = H-1 x0:

H <- M$vectors
b <- solve(H) %*% x0; b
##          [,1]
## [1,] -94.2809
## [2,] -74.5356

Tercero, los coeficientes de la función poli-exponencial del compartimento i es igual al producto de los elementos del autovector hi por los elementos del vector b:

A <- cbind(H[1, ] * b , H[2, ] * b); A
##          [,1]      [,2]
## [1,] 66.66667 -66.66667
## [2,] 33.33333  66.66667

El resultado obtenido es por lo tanto:

\[\begin{align} x_1(t)~=~66.666~e^{-0.4~t} + 33.333~e^{-0.1~t} \\ x_2(t)~=~-66.666~e^{-0.4~t} + 66.666~e^{-0.1~t} \end{align}\]

Multiplicando la matriz de observación B (en este caso un vector fila) por la traspuesta de la matriz de coeficientes A obtenemos la matriz de los coeficientes de la función del tiempo de la magnitud observada, la concentración plasmática de fármaco:

B %*% t(A)
##          [,1]     [,2]
## [1,] 3.333333 1.666667

Función hlsc()

En el código siguiente se difenen dos funciones R: hlsc() y f1()

  • La función hlsc() (homogeneous linear systems coefficients) “encapsula” el método descrito. Sus argumentos son la matriz del sistema y el vector de condiciones iniciales, y devuelve el vector de autovalores y la traspuesta de la matriz de coeficientes, por lo que los coeficientes de los compartimentos están ordenados por filas.
  • La función f1() se utiliza para la simulación del modelo. Los argumentos son el vector de coeficientes del compartimento, el vector de autovalores y el vector de tiempos.
hlsc <- function(K, x0) {
  dm <- dim(K)
  n <- dim(x0)[1]
  if(dm[1] != dm[2]) stop("wrong matrix K dimmension")
  if(dm[1] != n) stop("wrong x0 length")
  M <- eigen(K)
  l <- M$values
  H <- M$vectors
  b <- solve(H) %*% x0
  A <- c()
  for(i in 1:n) A <- rbind(A, t(H[i,] * b))
  list(A=A, l=l)
}

f1 <- function(A, l, x) {
  f0 <- function(t) sum(A * exp(l*t))
  sapply(x, f0)
}

Aplicaciones

Caso 1

Schwartz y cols. (2006) propusieron un método para el cálculo de la velocidad de filtración glomerular basado en la farmacocinética del iohexol utilizando un modelo de dos compartimentos. Posteriormente, Taubert y cols. (2018) analizaron nuevamente los datos utilizando el modelo de 3 compartimentos mostrado en la figura 1-B; los valores medios de los parámetros poblacionales obtenidos por estos autores fueron los siguientes:

x0 <- matrix(c(3235, 0, 0), ncol = 1)              # dosis, mg
Cl <- 3.44                                         # aclaramiento plasmático, l/h
V1 <- 6.40 ; V2 <- 5.43 ; V3 <- 3.26               # volúmenes distribución, l
Q12 <- 3.45 ; Q13 <- 12                            # aclaramientos intercompartimentales
c(Cl = Cl, V1 = V1, V2 = V2, Q12 = Q12, Q13 = Q13)
##    Cl    V1    V2   Q12   Q13 
##  3.44  6.40  5.43  3.45 12.00

Matriz de coeficientes y autovalores (recordar que la función hlsc() devuelve la matriz de coeficientes por filas para cada compartimento):

K <- matrix(c(-(Cl+Q12+Q13)/V1 , Q12/V1 , Q13/V1,
                Q12/V2 , -Q12/V2 , 0,
                Q13/V3 , 0, -Q13/V3),
              byrow = TRUE, ncol = 3)
ej2 <- hlsc(K, x0)
ej2
## $A
##           [,1]       [,2]      [,3]
## [1,]  1406.914   857.0351  971.0511
## [2,]  -166.741 -1233.2694 1400.0104
## [3,] -2236.725  1211.4509 1025.2737
## 
## $l
## [1] -5.9963421 -1.0768888 -0.1946723

Las soluciones son por tanto:

\[\begin{align} c_1(t)~=~1406.91~e^{-5.996~t} + 857.04~e^{-1.079~t} + 971.05~e^{-0.1947~t} \\ c_2(t)~=~-166.74~e^{-5.996~t} - 1233.27~e^{-1.079~t} + 1400.01~e^{-0.1947~t} \\ c_3(t)~=~-2236.72~e^{-5.996~t} + 1211.45~e^{-1.079~t} + 1025.27~e^{-0.1947~t} \\ \end{align}\]

Simulación:

sc2 <- data.frame(
tx = seq(0.001, 12, length.out = 201))
sc2 <- within(sc2, {
              x1 <- f1(ej2$A[1, ], ej2$l, tx)
              x2 <- f1(ej2$A[2, ], ej2$l, tx)
              x3 <- f1(ej2$A[3, ], ej2$l, tx)})

plot(x1 ~ tx, data = sc2, log = "y", type = "l",
     ylab = "concentración (mg/l)", xlab = "tiempo (h)")
lines(x2 ~ tx, data = sc2, col = "red")
lines(x3 ~ tx, data = sc2, col = "blue")
legend("topright", c("c1(t)", "c2(t)", "c3(t)"), lwd = 1, col = c(1,2,3))

Caso 2

Uno y cols. desarrollaron un modelo experimental de 3 compartimentos para predecir la concentración intersticial de taloporfina sódica en el miocardio (ver figura 2). Los compartimentos 1, 2 y 3 se refieren a los espacios vascular, intersticial e intracelular.

Las magnitudes observadas fueron la concentración plasmática (c1) y la concentración en el miocardio medida por fluorescencia, igual al promedio ponderado de las concentraciones en los tres compartimentos:

\[\begin{equation} c_{myo}~=~0.080~c_1 + 0.189~c_2 + 0.731~c_3 \end{equation}\]

 

Figura 2. Modelo de tres compartimentos para el análisis de la farmacocinética de la taloporfina

 

Matriz de coeficientes y autovalores:

### constantes de velocidad en 1/h
### volúmenes aparentes de distribución en l

k10 <- 1.18 ; k12 <- 4.37 ; k21 <- 7.85
k23 <- 6.41 ; k32 = 0.75
V1 <- 0.394 ; V2 <- 0.251 ; V3 <- 0.970

x0 <- matrix(c(2.5/V1, 0, 0), ncol = 1)

K <- matrix(c(-(k10+k12) , k21*V2/V1, 0,
              k12*V1/V2, -(k21+k23), k32*V3/V2,
              0, k12*V2/V3, - k32),
            byrow = TRUE, ncol = 3)
ej3 <- hlsc(K, x0)
ej3
## $A
##            [,1]      [,2]      [,3]
## [1,]  1.2404468  4.687777 0.4169543
## [2,] -2.9298133  2.492724 0.4370893
## [3,]  0.1994398 -1.316697 1.1172571
## 
## $l
## [1] -17.3616062  -2.8907789  -0.3076149

Simulación:

par(mfrow = c(1,2))

sc3 <- data.frame(tx <- seq(0.001, 5, length.out = 201))
sc3 <- within(sc3, {
  c1 <- f1(ej3$A[1, ], ej3$l, tx)
  c2 <- f1(ej3$A[2, ], ej3$l, tx)
  c3 <- f1(ej3$A[3, ], ej3$l, tx)
})

sc3 <- within(sc3, cmyo <- 0.080*c1 + 0.189*c2 + 0.731*c3)

plot(c1 ~ tx, data = sc3, type = "l",
     xlab = "tiempo (h)", ylab = "concentración (mg/(l·kg de peso))")
lines(c2 ~ tx, data = sc3, col = "red")
lines(c3 ~ tx, data = sc3, col = "blue")
legend("topright", c("c1(t)", "c2(t)", "c3(t)"), lwd = 1, col = c(1,2,3))

plot(c1 ~ tx, data = sc3, type = "l",
     xlab = "tiempo (h)", ylab = "concentración (mg/(l·kg de peso))")
lines(cmyo ~ tx, data = sc3, col = "red")
legend("topright", c("c1(t)", "cmyo(t)"), col = c(1,2))

Comentarios

1. En este documento sólo se ha aplicado el método de las matrices a la resolución numérica de modelos farmacocinéticos y a la simulación. Otras aplicaciones son el análisis de los momentos estadísticos y análisis de sensibilidad.

2. La precisión con la que se puede calcular los autovalores depende de que la matriz del sistema esté bién condicionada. Para matrices mal condicionadas, pequeñas variaciones en los términos de las matrices dan lugar a variaciones muy grandes en los autovalores de la matriz, llegando incluso a que aparezcan autovalores complejos. Para más información sobre métodos numéricos en R ver

3. Para más información sobre la resolución numérica de sistemas de ecuaciones diferenciales en R ver