Problem statement

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

where $A$ is the administered parent (or prodrug), $B$ is the active moiety that drives the pharmacodynamic (PD) effect, and $C$ is an inactive product subsequently cleared. Removal of $C$ is irrelevant to the time-course of $B$ and is ignored throughout.

Clinical/PD linkage. Assume the pharmacodynamic effect is directly and instantaneously driven by the active metabolite concentration, $E(t)\propto B(t)$. Consequently, the time of peak effect equals the time of peak active metabolite concentration, denoted $t_\text{peak}$.

Known

Elimination half-life of the active metabolite $B$: $T_{B\to C}>0$, giving the first-order elimination rate $k_e=\ln 2/T_{B\to C}$.
Observed time to peak active metabolite concentration (and effect): $t_\text{peak}>0$.
Initial condition reflecting an IV bolus (dose normalized): $A(0)=1,\;B(0)=0$.

The target is to compute the formation half-life of the parent to the active metabolite: $T_{A\to B}=\ln 2/k_f$ based solely on the known parameters.

Pharmacokinetic Model

This is a linear, one-way, first-order cascade with no feedback:

Formation of $B$ from $A$ at rate $k_f$ (units: 1/time). Pharmacologically, $k_f$ captures the net appearance of the active moiety from its precursor (e.g., metabolic activation), assuming distribution is fast relative to formation/elimination for this simplified analysis.
Elimination of $B$ at rate $k_e$ (units: 1/time), which is determined by the measured half-life $T_{B\to C}$.

\[ \dot A(t)=-k_f A(t), \qquad \dot B(t)=k_f A(t)-k_e B(t), \]\[ A(t)=e^{-k_f t}. \]\[ 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)$ is the difference of two exponentials: it rises initially as formation dominates, then falls as elimination dominates. The peak occurs where the instantaneous inflow and outflow of $B$ are equal.

2) Peak time of the active metabolite (and effect)

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

This formula is dimensionally consistent (time on both sides) and valid for $k_f\neq k_e$.

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

i.e., when formation and elimination are equally fast, the peak occurs at one mean lifetime and $T_{A\to B}=T_{B\to C}$.

3) Inverting the peak-time formula to obtain $k_f$

\[ 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. \]\[ a=\frac{\ln r}{r-1}. \]\[ \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} \]\[ -a r = W\!\big(-a e^{-a}\big) \quad\Longrightarrow\quad r = -\frac{1}{a}\,W\!\big(-a e^{-a}\big). \]\[ \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}}. \]

What is Lambert $W$? It is the multivalued inverse of $x\mapsto x e^{x}$. Many scientific libraries implement it as lambertw. In our case, the argument $z=-a e^{-a}$ lies in $(-e^{-1},0)$, where two real branches exist: the principal branch $W_0$ and the lower branch $W_{-1}$.

Branch selection (why two possible $W$’s, and which one to use?)

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

Therefore:

If $r<1$ (formation slower than elimination, $k_f<k_e$), then $a>1$ (peak occurs after one mean lifetime of $B$).
If $r>1$ (formation faster than elimination, $k_f>k_e$), then $a<1$ (peak occurs before one mean lifetime of $B$).
If $r=1$, then $a=1$ (the special case above).

For $z=-a e^{-a}\in(-e^{-1},0)$, the two real branches behave as:

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

Rule (in terms of observables):

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

Final closed-form solution (ready-to-use)

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

with branch $b$ chosen via the rule above (compare $a$ with 1).

Unit check. $a$ and $z$ are dimensionless. The numerator has units of time, $W_b(z)$ is dimensionless, so $T_{A\to B}$ has units of time (as required).

Practical guidance

Inputs you need.
- $T_{B\to C}$: the terminal half-life of the active metabolite $B$ (from concentration–time data).
- $t_\text{peak}$: the time of maximum effect (or maximum $B$ if available). If effect is used, ensure effect is directly proportional to $B$ without delay (no hysteresis).
\[ a=\frac{(\ln 2)\,t_\text{peak}}{T_{B\to C}}. \]
Pick the branch.
- If $a>1$: use $W_0$ (formation slower than elimination).
- If $a<1$: use $W_{-1}$ (formation faster than elimination).
- If $a=1$: $T_{A\to B}=T_{B\to C}$.
Evaluate $T_{A\to B}$.
- Compute $z=-a e^{-a}$, then evaluate $W_b(z)$ in your math library; finally compute $T_{A\to B}=-(\ln 2)\,t_\text{peak}/W_b(z)$.
Sanity checks.
- If $t_\text{peak}$ is very early relative to $T_{B\to C}$ (small $a$), you should get $T_{A\to B}\ll T_{B\to C}$.
- If $t_\text{peak}$ is very late, you should get $T_{A\to B}\gg T_{B\to C}$.

Edge cases and assumptions

Equal rates $k_f=k_e$: handled as above with $t_\text{peak}=1/k_e$.
Model scope. The derivation assumes:
- First-order (linear) kinetics for formation and elimination.
- A single well-mixed compartment for $B$ (no distributional delays).
- The PD effect is proportional to $B$ without delay (no effect compartment / tolerance).
  Violations (e.g., saturable metabolism, multi-compartment distribution, indirect response) will shift $t_\text{peak}$ and break the mapping.
Measurement noise. $t_\text{peak}$ estimated from sparse sampling may be biased. Smooth or model-fit the $B$ (or effect) time course to estimate the true peak time.

Python code to compute the uptake half-life

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

Example usecase: Uptake of quetiapine

For this example, we’ll use the uptake of oral quetiapine (metabolite A) to quetiapine (metabolite B).

The experimental pharmacological data according to NCBI is:

approximately 1 to 2 hours => we’ll use the average 1.5h for $T_{A\to B}$
The reported half-life (t1/2) for quetiapine is about 7 hours. for $T_{B\to C}$

As is evident from the following plot, the literature peak plasma concentration matches the analytical simulation based on the computed uptake half-life.

Quetiapine pharmacokinetics.svg

Output:

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

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

Pharmacokinetics: An analytical solution for computing the uptake half-life from time of peak plasma concentration