Problémové zadání

Zvažujeme kaskádu parent-metabolit relevantní pro farmakokinetiku/farmakodynamiku (PK/PD):

\[ A \xrightarrow{k_f} B \xrightarrow{k_e} C, \]

kde $A$ je podaný parent (nebo prodrug), $B$ je aktivní moiety, která pohání farmakodynamický (PD) efekt, a $C$ je neaktivní produkt následně vylučovaný. Odstranění $C$ je irelevantní pro časový průběh $B$ a je ignorováno.

Klinické/PD spojení. Předpokládejme, že farmakodynamický efekt je přímo a okamžitě poháněn koncentrací aktivního metabolitu, $E(t)\propto B(t)$. V důsledku toho se čas vrcholového efektu rovná času vrcholové koncentrace aktivního metabolitu, označenému $t_\text{peak}$.

Známé

Eliminační poločas aktivního metabolitu $B$: $T_{B\to C}>0$, dávající eliminační rychlost prvního řádu $k_e=\ln 2/T_{B\to C}$.
Pozorovaný čas do vrcholové koncentrace aktivního metabolitu (a efektu): $t_\text{peak}>0$.
Počáteční podmínka odrážející IV bolus (dávka normalizovaná): $A(0)=1,\;B(0)=0$.

Cílem je vypočítat formační poločas parentu na aktivní metabolit: $T_{A\to B}=\ln 2/k_f$ na základě pouze známých parametrů.

Farmakokinetický model

Toto je lineární, jednosměrná, prvního řádu kaskáda bez zpětné vazby:

Formace $B$ z $A$ rychlostí $k_f$ (jednotky: 1/čas). Farmakologicky $k_f$ zachycuje čistý výskyt aktivní moiety z jejího prekurzoru (např. metabolická aktivace), za předpokladu, že distribuce je rychlá relativně k formaci/eliminaci pro tuto zjednodušenou analýzu.
Eliminace $B$ rychlostí $k_e$ (jednotky: 1/čas), která je určena měřeným poločasem $T_{B\to C}$.

Mass-balance ODE jsou

\[ \dot A(t)=-k_f A(t), \qquad \dot B(t)=k_f A(t)-k_e B(t), \]

s $A(0)=1$ a $B(0)=0$. To implikuje

\[ A(t)=e^{-k_f t}. \]

Dosazením do $B$-rovnice dostaneme lineární ODE se známým vstupem $k_f A(t)$. Řešení (přes integrační faktor) dává, pro $k_f\neq k_e$,

\[ B(t)=\frac{k_f}{k_e-k_f}\Big(e^{-k_f t}-e^{-k_e t}\Big),\qquad t\ge 0. \]

Interpretace. $B(t)$ je rozdíl dvou exponenciál: zpočátku stoupá, když dominuje formace, pak klesá, když dominuje eliminace. Vrchol nastává tam, kde jsou okamžitý přítok a odtok $B$ stejné.

2) Vrcholový čas aktivního metabolitu (a efektu)

Derivujte $B(t)$ a nastavte na nulu:

\[ \frac{dB}{dt} =\frac{k_f}{k_e-k_f}\Big(-k_f e^{-k_f t} + k_e e^{-k_e t}\Big)=0 \quad\Longrightarrow\quad k_f e^{-k_f t}=k_e e^{-k_e t}. \]

Vezmeme-li přirozené logaritmy:

\[ \ln k_f - k_f t = \ln k_e - k_e t \;\;\Longrightarrow\;\; t_\text{peak}=\frac{\ln(k_f/k_e)}{\,k_f-k_e\,}. \]

Tento vzorec je dimenzionálně konzistentní (čas na obou stranách) a platný pro $k_f\neq k_e$.

Speciální případ $k_f=k_e$. Použitím l’Hospitalova pravidla nebo řešením degenerované ODE dostaneme

\[ B(t)=k_f t\,e^{-k_f t},\qquad t_\text{peak}=1/k_f=\frac{T_{B\to C}}{\ln 2}, \]

tj. když jsou formace a eliminace stejně rychlé, vrchol nastává při jednom průměrném čase života a $T_{A\to B}=T_{B\to C}$.

3) Inverze vzorce pro vrcholový čas pro získání $k_f$

Naším cílem je vypočítat $k_f$ (tedy $T_{A\to B}$) z známého $k_e$ a pozorovaného $t_\text{peak}$. Definujme bezrozměrný poměr

\[ r \equiv \frac{k_f}{k_e}>0,\qquad a \equiv k_e\,t_\text{peak}=\frac{(\ln 2)\,t_\text{peak}}{T_{B\to C}}>0. \]

Relace pro vrcholový čas se stává

\[ a=\frac{\ln r}{r-1}. \]

Tato transcendentální rovnice má uzavřené řešení pomocí Lambertovy $W$ funkce. Postupujte takto:

\[ \begin{aligned} a(r-1)&=\ln r &&\Longrightarrow& e^{a(r-1)}&=r\\ &&&\Longrightarrow& r\,e^{-a r}&=e^{-a}\\ &&&\Longrightarrow& (-a r)\,e^{-a r}&=(-a)\,e^{-a}. \end{aligned} \]

Z definice Lambertovy $W$ funkce ($W(z)\,e^{W(z)}=z$) identifikujeme

\[ -a r = W\!\big(-a e^{-a}\big) \quad\Longrightarrow\quad r = -\frac{1}{a}\,W\!\big(-a e^{-a}\big). \]

Získáme $k_f=r\,k_e$, pak převedeme na požadovaný poločas

\[ \boxed{\;T_{A\to B}=\frac{\ln 2}{k_f} = -\,\frac{(\ln 2)\,t_\text{peak}}{\,W_b\!\big(-a e^{-a}\big)}\;}, \qquad a=\frac{(\ln 2)\,t_\text{peak}}{T_{B\to C}}. \]

Co je Lambertova $W$? Je to vícehodnotová inverze $x\mapsto x e^{x}$. Mnoho vědeckých knihoven ji implementuje jako lambertw. V našem případě argument $z=-a e^{-a}$ leží v $(-e^{-1},0)$, kde existují dvě reálné větve: hlavní větev $W_0$ a dolní větev $W_{-1}$.

Výběr větve (proč dvě možná $W$, a které použít?)

Mapování $r\mapsto a=\frac{\ln r}{r-1}$ je monotonně klesající pro $r>0, r\neq 1$, s limitami

\[ \lim_{r\to 0^+} a=+\infty,\qquad a(1)=1,\qquad \lim_{r\to\infty} a=0^+. \]

Proto:

Pokud $r<1$ (formace pomalejší než eliminace, $k_f<k_e$), pak $a>1$ (vrchol nastává po jednom průměrném čase života $B$).
Pokud $r>1$ (formace rychlejší než eliminace, $k_f>k_e$), pak $a<1$ (vrchol nastává před jedním průměrným časem života $B$).
Pokud $r=1$, pak $a=1$ (speciální případ výše).

Pro $z=-a e^{-a}\in(-e^{-1},0)$ se dvě reálné větve chovají jako:

$W_0(z)\in(-1,0)$ $\Rightarrow$ dává $r\in(0,1)$.
$W_{-1}(z)\in(-\infty,-1]$ $\Rightarrow$ dává $r\in[1,\infty)$.

Pravidlo (v termínech pozorovatelných veličin):

\[ b=\begin{cases} 0, & a\gt 1 \quad (\text{použij } W_0 \Rightarrow k_f \lt k_e \Rightarrow T_{A\to B} \gt T_{B\to C}),\\ -1,& a\lt 1 \quad (\text{použij } W_{-1} \Rightarrow k_f \gt k_e \Rightarrow T_{A\to B} \lt T_{B\to C}),\\ \text{buď},& a=1 \quad (\text{dává } k_f=k_e). \end{cases} \]

Konečné uzavřené řešení (připravené k použití)

Nechť

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

Pak

\[ \boxed{\;T_{A\to B} = -\,\frac{(\ln 2)\,t_\text{peak}}{\,W_b(z)\,}\;}, \]

s větví $b$ vybranou pomocí pravidla výše (porovnejte $a$ s 1).

Kontrola jednotek. $a$ a $z$ jsou bezrozměrné. Čitatel má jednotky času, $W_b(z)$ je bezrozměrné, takže $T_{A\to B}$ má jednotky času (jak je vyžadováno).

Praktické pokyny

Vstupy, které potřebujete.
- $T_{B\to C}$: terminální poločas aktivního metabolitu $B$ (z koncentračně-časových dat).
- $t_\text{peak}$: čas maximálního efektu (nebo maximálního $B$, pokud je k dispozici). Pokud je použit efekt, ujistěte se, že efekt je přímo úměrný $B$ bez zpoždění (bez hysteréze).
Vypočítejte bezrozměrné $a$.
\[ a=\frac{(\ln 2)\,t_\text{peak}}{T_{B\to C}}. \]
Vyberte větev.
- Pokud $a>1$: použijte $W_0$ (formace pomalejší než eliminace).
- Pokud $a<1$: použijte $W_{-1}$ (formace rychlejší než eliminace).
- Pokud $a=1$: $T_{A\to B}=T_{B\to C}$.
Vyhodnoťte $T_{A\to B}$.
- Vypočítejte $z=-a e^{-a}$, pak vyhodnoťte $W_b(z)$ ve vaší matematické knihovně; konečně vypočítejte $T_{A\to B}=-(\ln 2)\,t_\text{peak}/W_b(z)$.
Kontroly správnosti.
- Pokud je $t_\text{peak}$ velmi brzy relativně k $T_{B\to C}$ (malé $a$), měli byste dostat $T_{A\to B}\ll T_{B\to C}$.
- Pokud je $t_\text{peak}$ velmi pozdě, měli byste dostat $T_{A\to B}\gg T_{B\to C}$.

Okrajové případy a předpoklady

Stejné rychlosti $k_f=k_e$: ošetřeno jako výše s $t_\text{peak}=1/k_e$.
Rozsah modelu. Derivace předpokládá:
- Kinetiku prvního řádu (lineární) pro formaci a eliminaci.
- Jednu dobře promíchanou kompartment pro $B$ (bez distribučních zpoždění).
- PD efekt je úměrný $B$ bez zpoždění (bez efektní kompartment / tolerance). Porušení (např. saturabilní metabolismus, multikompartment distribuce, nepřímá odpověď) posune $t_\text{peak}$ a rozbije mapování.
Měřící šum. $t_\text{peak}$ odhadnuté z řídkého vzorkování může být zkreslené. Vyhladit nebo fitovat model časového průběhu $B$ (nebo efektu) pro odhad skutečného času vrcholu.

Python kód pro výpočet poločasu absorpce

compute_formation_half_life.py

#!/usr/bin/env python3
# Compute T_{A->B} from T_{B->C} and t_peak using Lambert W (NumPy + SciPy).
import numpy as np
from scipy.special import lambertw

def formation_half_life(T_BC, t_peak):
    """
    Return T_{A->B} (formation half-life) in the A->B->C chain, given:
      - T_BC   : half-life of B->C (same time units as desired output)
      - t_peak : time of peak B (and peak effect if E ∝ B)

    Closed form:
        T_AB = - (ln 2) * t_peak / W_b( -a * e^{-a} ),
      where a = (ln 2) * t_peak / T_BC
    Branch rule (real solution):
        if a > 1 -> use W_0   (formation slower than elimination)
        if a < 1 -> use W_{-1} (formation faster than elimination)
        if 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 and t_peak must be positive.")

    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)

    # Special case: 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 -> principal branch 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)         # complex dtype
        T_AB[mask_gt] = -(ln2 * t_peak_arr[mask_gt]) / W0.real

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

    # Return scalar if inputs were scalars
    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):
    """Convenience: return k_f = ln(2) / T_{A->B} (same broadcasting behavior)."""
    T_AB = formation_half_life(T_BC, t_peak)
    return np.log(2.0) / T_AB


if __name__ == "__main__":
    # Examples
    # 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) formation slower (a > 1) -> T_AB > T_BC
    print("Example a>1:", formation_half_life(8.0, 20.0))

    # 3) formation faster (a < 1) -> T_AB < T_BC
    print("Example a<1:", formation_half_life(12.0, 2.0))

    # 4) vectorized inputs
    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))

Příklad použití: Absorpce quetiapinu

Pro tento příklad použijeme absorpci perorálního quetiapinu (metabolit A) na quetiapin (metabolit B).

Experimentální farmakologická data podle NCBI jsou:

přibližně 1 až 2 hodiny => použijeme průměr 1.5h pro $T_{A\to B}$
Reportovaný poločas (t1/2) pro quetiapin je asi 7 hodin. pro $T_{B\to C}$

Jak je zřejmé z následujícího grafu, literární vrcholová plazmatická koncentrace odpovídá analytické simulaci založené na vypočítaném poločasu absorpce.

Výstup:

quetiapine_example_output.txt

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

plot_quetiapine_pk.py

#!/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    # quetiapine elimination half-life (B -> elimination)
t_peak_hours = 1.5  # Tmax (time to peak plasma concentration)

def compute_T_AB(T_BC, t_peak):
    """Return uptake half-life T_AB from T_BC and t_peak.
       Uses Lambert W if SciPy is available; otherwise bisection on a=ln r/(r-1)."""
    if T_BC <= 0 or t_peak <= 0:
        raise ValueError("T_BC and t_peak must be positive.")
    a = ln2 * t_peak / T_BC
    if abs(a - 1.0) <= 1e-12:
        return T_BC  # equal rates case

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

# Compute half-lives and rates
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")

# Time grid and normalized amounts (A(0)=1)
t = np.linspace(0.0, 24.0, 1000)  # hours
A = np.exp(-kf * t)
B = (kf / (ke - kf)) * (np.exp(-kf * t) - np.exp(-ke * t))  # normalized Bateman form

# Plot
plt.figure(figsize=(8,5))
plt.plot(t, A, label="A(t): oral quetiapine (normalized)")
plt.plot(t, B, label="B(t): quetiapine (normalized)")
plt.axvline(t_peak_hours, linestyle="--", label=f"Literature Tmax ≈ {t_peak_hours:.2f} h")
plt.xlabel("Time (hours)")
plt.ylabel("Normalized amount / concentration (a.u.)")
plt.title("Quetiapine — derived uptake half-life and normalized PK curves")
plt.legend()
plt.tight_layout()
plt.savefig("Quetiapine pharmacokinetics.svg")

Farmakokinetika: Analytické řešení pro výpočet poločasu absorpce z času peak plazmatické koncentrace