Planteamiento del problema

Consideramos una cascada parente-metabolito relevante para la farmacocinética/farmacodinamia (PK/PD):

donde \(A\) es el parent administrado (o profármaco), \(B\) es la moiety activa que impulsa el efecto farmacodinámico (PD), y \(C\) es un producto inactivo posteriormente eliminado. La eliminación de \(C\) es irrelevante para la evolución temporal de \(B\) y se ignora en todo el desarrollo.

Vinculación clínica/PD. Supongamos que el efecto farmacodinámico está directa e instantáneamente impulsado por la concentración del metabolito activo, \(E(t)\propto B(t)\). En consecuencia, el tiempo de efecto máximo coincide con el tiempo de concentración máxima del metabolito activo, denotado \(t_\text{peak}\).

Datos conocidos

• Semivida de eliminación del metabolito activo \(B\): \(T_{B\to C}>0\), que da la tasa de eliminación de primer orden \(k_e=\ln 2/T_{B\to C}\).
• Tiempo observado hasta la concentración máxima del metabolito activo (y del efecto): \(t_\text{peak}>0\).
• Condición inicial que refleja un bolus IV (dosis normalizada): \(A(0)=1,\;B(0)=0\).

El objetivo es calcular la semivida de formación del parent al metabolito activo: \(T_{A\to B}=\ln 2/k_f\) basándose únicamente en los parámetros conocidos.

Modelo farmacocinético

Se trata de una cascada lineal unidireccional de primer orden sin retroalimentación:

• Formación de \(B\) a partir de \(A\) a una tasa \(k_f\) (unidades: 1/tiempo). Farmacológicamente, \(k_f\) captura la aparición neta de la moiety activa a partir de su precursor (p. ej., activación metabólica), asumiendo que la distribución es rápida en comparación con la formación/eliminación para este análisis simplificado.
• Eliminación de \(B\) a una tasa \(k_e\) (unidades: 1/tiempo), que está determinada por la semivida medida \(T_{B\to C}\).

Las EDOs de balance de masa son

con \(A(0)=1\) y \(B(0)=0\). Estas implican

Sustituyendo en la ecuación de \(B\) se obtiene una EDO lineal con entrada conocida \(k_f A(t)\). Resolviendo (mediante factor integrante) se obtiene, para \(k_f\neq k_e\),

Interpretación. \(B(t)\) es la diferencia de dos exponenciales: inicialmente aumenta cuando la formación domina, y luego disminuye cuando la eliminación domina. El pico ocurre donde el flujo de entrada y salida instantáneo de \(B\) son iguales.

2) Tiempo de pico del metabolito activo (y del efecto)

Diferenciamos \(B(t)\) e igualamos a cero:

Tomando logaritmos naturales:

Esta fórmula es dimensionalmente consistente (tiempo en ambos lados) y válida para \(k_f\neq k_e\).

Caso especial \(k_f=k_e\). Usando la regla de l’Hôpital o resolviendo la EDO degenerada, se obtiene

es decir, cuando la formación y la eliminación son igual de rápidas, el pico ocurre en una vida media promedio y \(T_{A\to B}=T_{B\to C}\).

3) Inversión de la fórmula del tiempo de pico para obtener \(k_f\)

Nuestro objetivo es calcular \(k_f\) (y por tanto \(T_{A\to B}\)) a partir de los valores conocidos de \(k_e\) y el \(t_\text{peak}\) observado. Definimos el cociente adimensional

La relación del tiempo de pico se convierte en

Esta ecuación trascendente tiene una solución en forma cerrada usando la función \(W\) de Lambert. Procedemos como sigue:

Por definición de la función \(W\) de Lambert (\(W(z)\,e^{W(z)}=z\)), identificamos

Recuperamos \(k_f=r\,k_e\), y luego convertimos a la semivida deseada

¿Qué es la función \(W\) de Lambert? Es la inversa multivaluada de \(x\mapsto x e^{x}\). Muchas bibliotecas científicas la implementan como lambertw. En nuestro caso, el argumento \(z=-a e^{-a}\) se encuentra en \((-e^{-1},0)\), donde existen dos ramas reales: la rama principal \(W_0\) y la rama inferior \(W_{-1}\).

Selección de rama (¿por qué dos \(W\) posibles, y cuál usar?)

La aplicación \(r\mapsto a=\frac{\ln r}{r-1}\) es monótona decreciente para \(r>0, r\neq 1\), con límites

Por lo tanto:

• Si \(r<1\) (formación más lenta que la eliminación, $k_f<k_e$), entonces $a>1$ (el pico ocurre después de una vida media promedio de \(B\)).
• Si \(r>1\) (formación más rápida que la eliminación, $k_f>k_e$), entonces $a<1$ (el pico ocurre antes de una vida media promedio de \(B\)).
• Si \(r=1\), entonces \(a=1\) (el caso especial anterior).

Para \(z=-a e^{-a}\in(-e^{-1},0)\), las dos ramas reales se comportan así:

• \(W_0(z)\in(-1,0)\) \(\Rightarrow\) da \(r\in(0,1)\).
• \(W_{-1}(z)\in(-\infty,-1]\) \(\Rightarrow\) da \(r\in[1,\infty)\).

Regla (en términos de observables):

Solución final en forma cerrada (lista para usar)

Sea

Entonces

con la rama \(b\) elegida mediante la regla anterior (comparar \(a\) con 1).

Comprobación de unidades. \(a\) y \(z\) son adimensionales. El numerador tiene unidades de tiempo, \(W_b(z)\) es adimensional, por lo que \(T_{A\to B}\) tiene unidades de tiempo (como se requiere).

Guía práctica

1. Datos que necesitas.

   • \(T_{B\to C}\): la semivida terminal del metabolito activo \(B\) (a partir de datos de concentración-tiempo).
   • \(t_\text{peak}\): el tiempo de efecto máximo (o \(B\) máximo si está disponible). Si se usa el efecto, asegúrate de que el efecto es directamente proporcional a \(B\) sin retardo (sin histéresis).

2. Calcula el \(a\) adimensional.

   \[ a=\frac{(\ln 2)\,t_\text{peak}}{T_{B\to C}}. \]

3. Elige la rama.

   • Si \(a>1\): usa \(W_0\) (formación más lenta que la eliminación).
   • Si \(a<1\): usa \(W_{-1}\) (formación más rápida que la eliminación).
   • Si \(a=1\): \(T_{A\to B}=T_{B\to C}\).

4. Evalúa \(T_{A\to B}\).

   • Calcula \(z=-a e^{-a}\), luego evalúa \(W_b(z)\) en tu biblioteca matemática; finalmente calcula \(T_{A\to B}=-(\ln 2)\,t_\text{peak}/W_b(z)\).

5. Comprobaciones de coherencia.

   • Si \(t_\text{peak}\) es muy temprano en relación con \(T_{B\to C}\) (\(a\) pequeño), deberías obtener \(T_{A\to B}\ll T_{B\to C}\).
   • Si \(t_\text{peak}\) es muy tardío, deberías obtener \(T_{A\to B}\gg T_{B\to C}\).

Casos límite y supuestos

• Tasas iguales \(k_f=k_e\): se maneja como arriba con \(t_\text{peak}=1/k_e\).
• Alcance del modelo. La derivación asume:
  • Cinética de primer orden (lineal) para formación y eliminación.
  • Un único compartimento bien mezclado para \(B\) (sin retardos distributivos).
  • El efecto PD es proporcional a \(B\) sin retardo (sin compartimento de efecto / tolerancia). Las violaciones (p. ej., metabolismo saturable, distribución multicompartmental, respuesta indirecta) desplazarán \(t_\text{peak}\) y romperán la correspondencia.
• Ruido de medición. El \(t_\text{peak}\) estimado a partir de un muestreo disperso puede estar sesgado. Suaviza o ajusta un modelo a la evolución temporal de \(B\) (o del efecto) para estimar el verdadero tiempo de pico.

Código Python para calcular la semivida de absorción

#!/usr/bin/env python3
# Calcular T_{A->B} a partir de T_{B->C} y t_peak usando Lambert W (NumPy + SciPy).
import numpy as np
from scipy.special import lambertw

def formation_half_life(T_BC, t_peak):
    """
    Devuelve T_{A->B} (semivida de formación) en la cadena A->B->C, dados:
      - T_BC   : semivida de B->C (mismas unidades de tiempo que la salida deseada)
      - t_peak : tiempo de pico de B (y del efecto máximo si E ∝ B)

    Forma cerrada:
        T_AB = - (ln 2) * t_peak / W_b( -a * e^{-a} ),
      donde a = (ln 2) * t_peak / T_BC
    Regla de rama (solución real):
        si a > 1 -> usar W_0   (formación más lenta que la eliminación)
        si a < 1 -> usar W_{-1} (formación más rápida que la eliminación)
        si a == 1 -> T_AB = T_BC
    """
    T_BC_arr, t_peak_arr = np.broadcast_arrays(np.asarray(T_BC, float),
                                               np.asarray(t_peak, float))

    if np.any(T_BC_arr <= 0) or np.any(t_peak_arr <= 0):
        raise ValueError("T_BC y t_peak deben ser positivos.")

    ln2 = np.log(2.0)
    a = ln2 * t_peak_arr / T_BC_arr
    z = -a * np.exp(-a)

    T_AB = np.empty_like(a, dtype=float)

    # Caso especial: a == 1 -> T_AB == T_BC
    mask_eq = np.isclose(a, 1.0, rtol=1e-12, atol=0.0)
    if np.any(mask_eq):
        T_AB[mask_eq] = T_BC_arr[mask_eq]

    # a > 1 -> rama principal W_0 -> k_f < k_e -> T_AB > T_BC
    mask_gt = a > 1.0
    if np.any(mask_gt):
        W0 = lambertw(z[mask_gt], 0)         # tipo complejo
        T_AB[mask_gt] = -(ln2 * t_peak_arr[mask_gt]) / W0.real

    # a < 1 -> rama inferior W_{-1} -> k_f > k_e -> T_AB < T_BC
    mask_lt = a < 1.0
    if np.any(mask_lt):
        Wm1 = lambertw(z[mask_lt], -1)       # tipo complejo
        T_AB[mask_lt] = -(ln2 * t_peak_arr[mask_lt]) / Wm1.real

    # Devolver escalar si las entradas eran escalares
    if np.isscalar(T_BC) and np.isscalar(t_peak):
        return float(T_AB.item())
    return T_AB


def formation_rate_constant(T_BC, t_peak):
    """Conveniencia: devuelve k_f = ln(2) / T_{A->B} (mismo comportamiento de broadcasting)."""
    T_AB = formation_half_life(T_BC, t_peak)
    return np.log(2.0) / T_AB


if __name__ == "__main__":
    # Ejemplos
    # 1) a == 1 -> T_AB == T_BC
    T_BC = 10.0
    t_peak = T_BC / np.log(2.0)
    print("Example a=1:", formation_half_life(T_BC, t_peak))

    # 2) formación más lenta (a > 1) -> T_AB > T_BC
    print("Example a>1:", formation_half_life(8.0, 20.0))

    # 3) formación más rápida (a < 1) -> T_AB < T_BC
    print("Example a<1:", formation_half_life(12.0, 2.0))

    # 4) entradas vectorizadas
    T_BC_vec = np.array([10.0, 8.0, 12.0])
    t_peak_vec = np.array([10.0/np.log(2.0), 20.0, 2.0])
    print("Vectorized:", formation_half_life(T_BC_vec, t_peak_vec))

Ejemplo de uso: Absorción de quetiapina

Para este ejemplo, usaremos la absorción de quetiapina oral (metabolito A) a quetiapina (metabolito B).

Los datos farmacológicos experimentales según NCBI son:

• aproximadamente de 1 a 2 horas => usaremos el promedio 1.5h para \(T_{A\to B}\)
• La semivida reportada (t1/2) para quetiapina es de aproximadamente 7 horas. para \(T_{B\to C}\)

Como se evidencia en la siguiente gráfica, la concentración plasmática máxima de la literatura coincide con la simulación analítica basada en la semivida de absorción calculada.

Salida:

T_BC = 6.0 h, Tmax = 1.5 h
Computed uptake half-life T_AB ≈ 0.3424 h (20.5 min)
kf = 2.0246 1/h, ke = 0.1155 1/h

#!/usr/bin/env python3
import numpy as np, math
import matplotlib.pyplot as plt
plt.style.use("ggplot")

from scipy.special import lambertw

ln2 = math.log(2.0)
T_BC_hours = 6.0    # semivida de eliminación de quetiapina (B -> eliminación)
t_peak_hours = 1.5  # Tmax (tiempo hasta la concentración plasmática máxima)

def compute_T_AB(T_BC, t_peak):
    """Devuelve la semivida de absorción T_AB a partir de T_BC y t_peak.
       Usa Lambert W si SciPy está disponible; en caso contrario bisección sobre a=ln r/(r-1)."""
    if T_BC <= 0 or t_peak <= 0:
        raise ValueError("T_BC y t_peak deben ser positivos.")
    a = ln2 * t_peak / T_BC
    if abs(a - 1.0) <= 1e-12:
        return T_BC  # caso de tasas iguales

    # Calcular usando Lambert W (SciPy presente)
    z = -a * math.exp(-a)
    branch = 0 if a > 1 else -1
    W = lambertw(z, branch)
    return float(-(ln2 * t_peak) / W.real)

# Calcular semividas y tasas
T_AB_hours = compute_T_AB(T_BC_hours, t_peak_hours)
ke = ln2 / T_BC_hours
kf = ln2 / T_AB_hours

print(f"T_BC = {T_BC_hours} h, Tmax = {t_peak_hours} h")
print(f"Computed uptake half-life T_AB ≈ {T_AB_hours:.4f} h ({T_AB_hours*60:.1f} min)")
print(f"kf = {kf:.4f} 1/h, ke = {ke:.4f} 1/h")

# Rejilla temporal y cantidades normalizadas (A(0)=1)
t = np.linspace(0.0, 24.0, 1000)  # horas
A = np.exp(-kf * t)
B = (kf / (ke - kf)) * (np.exp(-kf * t) - np.exp(-ke * t))  # forma Bateman normalizada

# Gráfica
plt.figure(figsize=(8,5))
plt.plot(t, A, label="A(t): quetiapina oral (normalizada)")
plt.plot(t, B, label="B(t): quetiapina (normalizada)")
plt.axvline(t_peak_hours, linestyle="--", label=f"Tmax de literatura ≈ {t_peak_hours:.2f} h")
plt.xlabel("Tiempo (horas)")
plt.ylabel("Cantidad / concentración normalizada (u.a.)")
plt.title("Quetiapina — semivida de absorción derivada y curvas PK normalizadas")
plt.legend()
plt.tight_layout()
plt.savefig("Quetiapine pharmacokinetics.svg")