Enunciado do problema

Consideramos uma cascata precursor–metabólito relevante para a farmacocinética/farmacodinâmica (PK/PD):

onde \(A\) é o precursor administrado (ou pró-fármaco), \(B\) é a fração ativa que determina o efeito farmacodinâmico (PD) e \(C\) é um produto inativo subsequentemente eliminado. A remoção de \(C\) é irrelevante para a evolução temporal de \(B\) e é ignorada em toda a análise.

Conexão clínica/PD. Assume-se que o efeito farmacodinâmico é direta e instantaneamente determinado pela concentração do metabólito ativo, \(E(t)\propto B(t)\). Consequentemente, o tempo de pico de efeito é igual ao tempo de pico de concentração do metabólito ativo, denotado \(t_\text{peak}\).

Dados conhecidos

• Meia-vida de eliminação do metabólito ativo \(B\): \(T_{B\to C}>0\), fornecendo a taxa de eliminação de primeira ordem \(k_e=\ln 2/T_{B\to C}\).
• Tempo observado até o pico de concentração do metabólito ativo (e do efeito): \(t_\text{peak}>0\).
• Condição inicial refletindo um bolus IV (dose normalizada): \(A(0)=1,\;B(0)=0\).

O objetivo é calcular a meia-vida de formação do precursor para o metabólito ativo: \(T_{A\to B}=\ln 2/k_f\) baseando-se exclusivamente nos parâmetros conhecidos.

Modelo farmacocinético

Esta é uma cascata linear, unidirecional, de primeira ordem sem retroalimentação:

• Formação de \(B\) a partir de \(A\) à taxa \(k_f\) (unidades: 1/tempo). Farmacologicamente, \(k_f\) captura a aparição líquida da fração ativa a partir de seu precursor (ex.: ativação metabólica), assumindo que a distribuição é rápida em relação à formação/eliminação nesta análise simplificada.
• Eliminação de \(B\) à taxa \(k_e\) (unidades: 1/tempo), determinada pela meia-vida medida \(T_{B\to C}\).

As EDOs de balanço de massa são

com \(A(0)=1\) e \(B(0)=0\). Estas implicam

Substituindo na equação de \(B\) obtém-se uma EDO linear com entrada conhecida \(k_f A(t)\). A resolução (via fator integrante) fornece, para \(k_f\neq k_e\),

Interpretação. \(B(t)\) é a diferença de duas exponenciais: inicialmente cresce enquanto a formação predomina, depois decai quando a eliminação predomina. O pico ocorre onde o fluxo de entrada e de saída instantâneos de \(B\) são iguais.

2) Tempo de pico do metabólito ativo (e do efeito)

Diferencie \(B(t)\) e iguale a zero:

Aplicando logaritmos naturais:

Esta fórmula é dimensionalmente consistente (tempo em ambos os lados) e válida para \(k_f\neq k_e\).

Caso especial \(k_f=k_e\). Usando a regra de l’Hôpital ou resolvendo a EDO degenerada, obtém-se

isto é, quando a formação e a eliminação são igualmente rápidas, o pico ocorre em um tempo de vida médio e \(T_{A\to B}=T_{B\to C}\).

3) Invertendo a fórmula do tempo de pico para obter \(k_f\)

Nosso objetivo é calcular \(k_f\) (e portanto \(T_{A\to B}\)) a partir dos conhecidos \(k_e\) e do \(t_\text{peak}\) observado. Defina a razão adimensional

A relação do tempo de pico torna-se

Esta equação transcendental possui uma solução em forma fechada usando a função Lambert \(W\). Proceda da seguinte forma:

Por definição da função Lambert \(W\) (\(W(z)\,e^{W(z)}=z\)), identificamos

Recupere \(k_f=r\,k_e\), depois converta para a meia-vida desejada

O que é a função Lambert \(W\)? É a inversa multivalorada de \(x\mapsto x e^{x}\). Muitas bibliotecas científicas a implementam como lambertw. No nosso caso, o argumento \(z=-a e^{-a}\) está em \((-e^{-1},0)\), onde existem dois ramos reais: o ramo principal \(W_0\) e o ramo inferior \(W_{-1}\).

Seleção do ramo (por que dois \(W\) possíveis, e qual usar?)

O mapeamento \(r\mapsto a=\frac{\ln r}{r-1}\) é monotonamente decrescente para \(r>0, r\neq 1\), com limites

Portanto:

• Se \(r<1\) (formação mais lenta que a eliminação, $k_f<k_e$), então $a>1$ (o pico ocorre após um tempo de vida médio de \(B\)).
• Se \(r>1\) (formação mais rápida que a eliminação, $k_f>k_e$), então $a<1$ (o pico ocorre antes de um tempo de vida médio de \(B\)).
• Se \(r=1\), então \(a=1\) (o caso especial acima).

Para \(z=-a e^{-a}\in(-e^{-1},0)\), os dois ramos reais comportam-se como:

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

Regra (em termos de observáveis):

Solução final em forma fechada (pronta para uso)

Seja

Então

com o ramo \(b\) escolhido pela regra acima (compare \(a\) com 1).

Verificação de unidades. \(a\) e \(z\) são adimensionais. O numerador tem unidades de tempo, \(W_b(z)\) é adimensional, portanto \(T_{A\to B}\) tem unidades de tempo (como requerido).

Orientação prática

1. Entradas necessárias.

   • \(T_{B\to C}\): a meia-vida terminal do metabólito ativo \(B\) (a partir de dados de concentração–tempo).
   • \(t_\text{peak}\): o tempo de efeito máximo (ou máximo de \(B\), se disponível). Se o efeito for usado, garanta que o efeito seja diretamente proporcional a \(B\) sem atraso (sem histerese).

2. Calcule o \(a\) adimensional.

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

3. Escolha o ramo.

   • Se \(a>1\): use \(W_0\) (formação mais lenta que a eliminação).
   • Se \(a<1\): use \(W_{-1}\) (formação mais rápida que a eliminação).
   • Se \(a=1\): \(T_{A\to B}=T_{B\to C}\).

4. Avalie \(T_{A\to B}\).

   • Calcule \(z=-a e^{-a}\), depois avalie \(W_b(z)\) em sua biblioteca matemática; finalmente calcule \(T_{A\to B}=-(\ln 2)\,t_\text{peak}/W_b(z)\).

5. Verificações de sanidade.

   • Se \(t_\text{peak}\) for muito cedo em relação a \(T_{B\to C}\) (\(a\) pequeno), você deve obter \(T_{A\to B}\ll T_{B\to C}\).
   • Se \(t_\text{peak}\) for muito tarde, você deve obter \(T_{A\to B}\gg T_{B\to C}\).

Casos-limite e premissas

• Taxas iguais \(k_f=k_e\): tratado como acima com \(t_\text{peak}=1/k_e\).
• Escopo do modelo. A derivação assume:
  • Cinética de primeira ordem (linear) para formação e eliminação.
  • Um único compartimento bem misturado para \(B\) (sem atrasos de distribuição).
  • O efeito PD é proporcional a \(B\) sem atraso (sem compartimento de efeito / tolerância). Violações (ex.: metabolismo saturável, distribuição multicompartimental, resposta indireta) deslocarão \(t_\text{peak}\) e romperão o mapeamento.
• Ruído de medição. \(t_\text{peak}\) estimado a partir de amostragem esparsa pode ser enviesado. Alise ou ajuste por modelo a evolução temporal de \(B\) (ou do efeito) para estimar o verdadeiro tempo de pico.

Código Python para calcular a meia-vida de absorção

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

def formation_half_life(T_BC, t_peak):
    """
    Retorna T_{A->B} (meia-vida de formação) na cadeia A->B->C, dados:
      - T_BC   : meia-vida de B->C (mesmas unidades de tempo da saída desejada)
      - t_peak : tempo de pico de B (e pico de efeito se E ∝ B)

    Forma fechada:
        T_AB = - (ln 2) * t_peak / W_b( -a * e^{-a} ),
      onde a = (ln 2) * t_peak / T_BC
    Regra de ramo (solução real):
        se a > 1 -> use W_0   (formação mais lenta que eliminação)
        se a < 1 -> use W_{-1} (formação mais rápida que eliminação)
        se 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 e t_peak devem 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 -> ramo 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)         # dtype complexo
        T_AB[mask_gt] = -(ln2 * t_peak_arr[mask_gt]) / W0.real

    # a < 1 -> ramo 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)       # dtype complexo
        T_AB[mask_lt] = -(ln2 * t_peak_arr[mask_lt]) / Wm1.real

    # Retorna escalar se as entradas forem 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):
    """Conveniência: retorna k_f = ln(2) / T_{A->B} (mesmo comportamento de broadcasting)."""
    T_AB = formation_half_life(T_BC, t_peak)
    return np.log(2.0) / T_AB


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

    # 2) formação mais lenta (a > 1) -> T_AB > T_BC
    print("Exemplo a>1:", formation_half_life(8.0, 20.0))

    # 3) formação mais rápida (a < 1) -> T_AB < T_BC
    print("Exemplo a<1:", formation_half_life(12.0, 2.0))

    # 4) entradas vetorizadas
    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("Vetorizado:", formation_half_life(T_BC_vec, t_peak_vec))

Exemplo de uso: Absorção de quetiapina

Para este exemplo, usaremos a absorção de quetiapina oral (metabólito A) em quetiapina (metabólito B).

Os dados farmacológicos experimentais de acordo com a NCBI são:

• aproximadamente 1 a 2 horas => usaremos a média 1,5h para \(T_{A\to B}\)
• A meia-vida relatada (t1/2) da quetiapina é de cerca de 7 horas. para \(T_{B\to C}\)

Como é evidente no gráfico a seguir, a concentração plasmática de pico da literatura corresponde à simulação analítica baseada na meia-vida de absorção calculada.

Saída:

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    # meia-vida de eliminação da quetiapina (B -> eliminação)
t_peak_hours = 1.5  # Tmax (tempo até a concentração plasmática de pico)

def compute_T_AB(T_BC, t_peak):
    """Retorna a meia-vida de absorção T_AB a partir de T_BC e t_peak.
       Usa Lambert W se SciPy estiver disponível; caso contrário, bissecção em a=ln r/(r-1)."""
    if T_BC <= 0 or t_peak <= 0:
        raise ValueError("T_BC e t_peak devem ser positivos.")
    a = ln2 * t_peak / T_BC
    if abs(a - 1.0) <= 1e-12:
        return T_BC  # caso de taxas iguais

    # Cálculo 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 meias-vidas e taxas
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"Meia-vida de absorção calculada 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")

# Grade temporal e quantidades 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 de Bateman normalizada

# Gráfico
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 da literatura ≈ {t_peak_hours:.2f} h")
plt.xlabel("Tempo (horas)")
plt.ylabel("Quantidade / concentração normalizada (u.a.)")
plt.title("Quetiapina — meia-vida de absorção derivada e curvas PK normalizadas")
plt.legend()
plt.tight_layout()
plt.savefig("Quetiapine pharmacokinetics.svg")