Problemstellung

Wir betrachten eine für die Pharmakokinetik/Pharmakodynamik (PK/PD) relevante Parent-Metabolit-Kaskade:

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

wobei $A$ der verabreichte Parent (oder Prodrug) ist, $B$ der aktive Wirkstoff ist, der den pharmakodynamischen (PD) Effekt antreibt, und $C$ ein inaktives Produkt ist, das anschließend ausgeschieden wird. Die Entfernung von $C$ ist für den Zeitverlauf von $B$ irrelevant und wird durchgehend ignoriert.

Klinische/PD-Verknüpfung. Angenommen, der pharmakodynamische Effekt wird direkt und augenblicklich durch die Konzentration des aktiven Metaboliten angetrieben, $E(t)\propto B(t)$. Folglich ist der Zeitpunkt der maximalen Wirkung gleich dem Zeitpunkt der maximalen aktiven Metabolitenkonzentration, bezeichnet als $t_\text{peak}$.

Bekannt

Eliminationshalbwertszeit des aktiven Metaboliten $B$: $T_{B\to C}>0$, was die Eliminationsrate erster Ordnung $k_e=\ln 2/T_{B\to C}$ ergibt.
Beobachtete Zeit bis zur maximalen aktiven Metabolitenkonzentration (und Wirkung): $t_\text{peak}>0$.
Anfangsbedingung entsprechend einem IV-Bolus (Dosis normalisiert): $A(0)=1,\;B(0)=0$.

Das Ziel ist die Berechnung der Bildungshalbwertszeit vom Parent zum aktiven Metaboliten: $T_{A\to B}=\ln 2/k_f$ basierend ausschließlich auf den bekannten Parametern.

Pharmakokinetisches Modell

Dies ist eine lineare, einseitige, erster Ordnung-Kaskade ohne Rückkopplung:

Bildung von $B$ aus $A$ mit Rate $k_f$ (Einheiten: 1/Zeit). Pharmakologisch erfasst $k_f$ das Netto-Erscheinen des aktiven Wirkstoffs aus seinem Vorläufer (z.B. metabolische Aktivierung), unter der Annahme, dass die Verteilung im Verhältnis zur Bildung/Elimination für diese vereinfachte Analyse schnell ist.
Elimination von $B$ mit Rate $k_e$ (Einheiten: 1/Zeit), die durch die gemessene Halbwertszeit $T_{B\to C}$ bestimmt wird.

Die Massenbilanz-ODEs sind

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

mit $A(0)=1$ und $B(0)=0$. Diese implizieren

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

Einsetzen in die $B$-Gleichung ergibt eine lineare ODE mit bekanntem Input $k_f A(t)$. Lösung (via Integrationsfaktor) ergibt, für $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. \]

Interpretation. $B(t)$ ist die Differenz zweier Exponentialfunktionen: Sie steigt anfangs, wenn die Bildung dominiert, und fällt dann, wenn die Elimination dominiert. Das Maximum tritt dort auf, wo der momentane Zufluss und Abfluss von $B$ gleich sind.

2) Peak-Zeit des aktiven Metaboliten (und der Wirkung)

Differenziere $B(t)$ und setze auf Null:

\[ \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}. \]

Natürliche Logarithmen nehmen:

\[ \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\,}. \]

Diese Formel ist dimensionskonsistent (Zeit auf beiden Seiten) und gültig für $k_f\neq k_e$.

Sonderfall $k_f=k_e$. Mit der Regel von l’Hospital oder durch Lösen der degenerierten ODE findet man

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

d.h. wenn Bildung und Elimination gleich schnell sind, tritt das Maximum bei einer mittleren Lebensdauer auf und $T_{A\to B}=T_{B\to C}$.

3) Invertieren der Peak-Zeit-Formel zur Bestimmung von $k_f$

Unser Ziel ist die Berechnung von $k_f$ (und damit $T_{A\to B}$) aus dem bekannten $k_e$ und dem beobachteten $t_\text{peak}$. Definiere das dimensionslose Verhältnis

\[ 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. \]

Die Peak-Zeit-Beziehung wird zu

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

Diese transzendente Gleichung hat eine geschlossene Lösung mit der Lambert-$W$-Funktion. Vorgehen wie folgt:

\[ \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} \]

Per Definition der Lambert-$W$-Funktion ($W(z)\,e^{W(z)}=z$) identifizieren wir

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

$k_f=r\,k_e$ zurückgewinnen, dann in die gewünschte Halbwertszeit umwandeln

\[ \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}}. \]

Was ist Lambert $W$? Es ist die mehrwertige Umkehrfunktion von $x\mapsto x e^{x}$. Viele wissenschaftliche Bibliotheken implementieren sie als lambertw. In unserem Fall liegt das Argument $z=-a e^{-a}$ in $(-e^{-1},0)$, wo zwei reelle Zweige existieren: der Hauptzweig $W_0$ und der untere Zweig $W_{-1}$.

Zweigauswahl (warum zwei mögliche $W$’s und welcher zu verwenden ist?)

Die Abbildung $r\mapsto a=\frac{\ln r}{r-1}$ ist monoton fallend für $r>0, r\neq 1$, mit Grenzwerten

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

Daher:

Wenn $r<1$ (Bildung langsamer als Elimination, $k_f<k_e$), dann $a>1$ (Maximum tritt nach einer mittleren Lebensdauer von $B$ auf).
Wenn $r>1$ (Bildung schneller als Elimination, $k_f>k_e$), dann $a<1$ (Maximum tritt vor einer mittleren Lebensdauer von $B$ auf).
Wenn $r=1$, dann $a=1$ (der Sonderfall oben).

Für $z=-a e^{-a}\in(-e^{-1},0)$ verhalten sich die beiden reellen Zweige wie:

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

Regel (in Bezug auf Beobachtungsgrößen):

\[ b=\begin{cases} 0, & a\gt 1 \quad (\text{verwende } W_0 \Rightarrow k_f \lt k_e \Rightarrow T_{A\to B} \gt T_{B\to C}),\\ -1,& a\lt 1 \quad (\text{verwende } W_{-1} \Rightarrow k_f \gt k_e \Rightarrow T_{A\to B} \lt T_{B\to C}),\\ \text{beliebig},& a=1 \quad (\text{ergibt } k_f=k_e). \end{cases} \]

Endgültige geschlossene Lösung (einsatzbereit)

Sei

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

Dann

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

mit Zweig $b$, gewählt über die Regel oben (vergleiche $a$ mit 1).

Einheitenprüfung. $a$ und $z$ sind dimensionslos. Der Zähler hat Einheiten von Zeit, $W_b(z)$ ist dimensionslos, also hat $T_{A\to B}$ Einheiten von Zeit (wie erforderlich).

Praktische Hinweise

Benötigte Eingaben.
- $T_{B\to C}$: die terminale Halbwertszeit des aktiven Metaboliten $B$ (aus Konzentrations-Zeit-Daten).
- $t_\text{peak}$: die Zeit der maximalen Wirkung (oder maximales $B$, falls verfügbar). Falls die Wirkung verwendet wird, stellen Sie sicher, dass die Wirkung direkt proportional zu $B$ ohne Verzögerung ist (keine Hysterese).
Berechnen Sie das dimensionslose $a$.
\[ a=\frac{(\ln 2)\,t_\text{peak}}{T_{B\to C}}. \]
Wählen Sie den Zweig.
- Wenn $a>1$: verwenden Sie $W_0$ (Bildung langsamer als Elimination).
- Wenn $a<1$: verwenden Sie $W_{-1}$ (Bildung schneller als Elimination).
- Wenn $a=1$: $T_{A\to B}=T_{B\to C}$.
Berechnen Sie $T_{A\to B}$.
- Berechnen Sie $z=-a e^{-a}$, dann evaluieren Sie $W_b(z)$ in Ihrer Mathematik-Bibliothek; schließlich berechnen Sie $T_{A\to B}=-(\ln 2)\,t_\text{peak}/W_b(z)$.
Plausibilitätsprüfungen.
- Wenn $t_\text{peak}$ sehr früh relativ zu $T_{B\to C}$ liegt (kleines $a$), sollten Sie $T_{A\to B}\ll T_{B\to C}$ erhalten.
- Wenn $t_\text{peak}$ sehr spät liegt, sollten Sie $T_{A\to B}\gg T_{B\to C}$ erhalten.

Randfälle und Annahmen

Gleiche Raten $k_f=k_e$: wie oben behandelt mit $t_\text{peak}=1/k_e$.
Modellumfang. Die Herleitung nimmt an:
- Kinetik erster Ordnung (linear) für Bildung und Elimination.
- Ein einzelnes gut-gemischtes Kompartiment für $B$ (keine Verteilungsverzögerungen).
- Der PD-Effekt ist proportional zu $B$ ohne Verzögerung (kein Effekt-Kompartiment / Toleranz). Verletzungen (z.B. sättigbare Metabolisierung, Mehrkompartiment-Verteilung, indirekte Antwort) verschieben $t_\text{peak}$ und brechen die Abbildung.
Messrauschen. $t_\text{peak}$, geschätzt aus spärlicher Probenahme, kann verzerrt sein. Glätten oder modell-anpassen Sie den $B$- (oder Wirkungs-) Zeitverlauf, um die wahre Peak-Zeit zu schätzen.

Python-Code zur Berechnung der Aufnahme-Halbwertszeit

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))

#!/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))

Beispielanwendungsfall: Aufnahme von Quetiapin

Für dieses Beispiel verwenden wir die Aufnahme von oralem Quetiapin (Metabolit A) zu Quetiapin (Metabolit B).

Die experimentellen pharmakologischen Daten gemäß NCBI sind:

approximately 1 to 2 hours => wir verwenden den Durchschnitt 1.5h für $T_{A\to B}$
The reported half-life (t1/2) for quetiapine is about 7 hours. für $T_{B\to C}$

Wie aus dem folgenden Plot ersichtlich, stimmt die literaturbekannte maximale Plasmakonzentration mit der analytischen Simulation basierend auf der berechneten Aufnahme-Halbwertszeit überein.

Ausgabe:

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

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")

#!/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")

Check out similar posts by category: Bioinformatics, Pharmacokinetics, Mathematics

If this post helped you, please consider buying me a coffee or donating via PayPal to support research & publishing of new posts on TechOverflow

Buy me a coffee

Pharmakokinetik: Eine analytische Lösung zur Berechnung der Aufnahme-Halbwertszeit aus dem Zeitpunkt der maximalen Plasmakonzentration