问题陈述

我们考虑一个与药代动力学/药效动力学（PK/PD）相关的母体-代谢物级联反应：

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

其中 $A$ 是给药的母体（或前药），$B$ 是驱动药效动力学（PD）效应的活性部分，$C$ 是随后被清除的无活性产物。$C$ 的清除与 $B$ 的时间进程无关，在全文中忽略。

临床/PD 关联。 假设药效动力学效应直接且瞬时地由活性代谢物浓度驱动，即 $E(t)\propto B(t)$。因此，峰值效应时间等于活性代谢物峰值浓度时间，记为 $t_\text{peak}$。

已知条件

活性代谢物 $B$ 的消除半衰期：$T_{B\to C}>0$，对应一级消除速率 $k_e=\ln 2/T_{B\to C}$。
观测到的活性代谢物达峰浓度（及效应）时间：$t_\text{peak}>0$。
反映静脉推注的初始条件（剂量归一化）：$A(0)=1,\;B(0)=0$。

目标是仅根据已知参数计算母体到活性代谢物的生成半衰期：$T_{A\to B}=\ln 2/k_f$。

药代动力学模型

这是一个线性、单向、一级级联反应，无反馈：

$B$ 从 $A$ 生成，速率为 $k_f$（单位：1/时间）。在药理学上，$k_f$ 表示活性部分从其前体净出现的速率（如代谢活化），假设在此简化分析中分布相对于生成/消除足够快。
$B$ 的消除速率为 $k_e$（单位：1/时间），由测量的半衰期 $T_{B\to C}$ 确定。

质量平衡 ODE 为

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

初始条件 $A(0)=1$ 和 $B(0)=0$。由此可得

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

代入 $B$ 的方程得到一个已知输入 $k_f A(t)$ 的线性 ODE。求解（通过积分因子）得，对于 $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. \]

解释。 $B(t)$ 是两个指数函数的差：初期生成占主导时上升，随后消除占主导时下降。峰值出现在 $B$ 的瞬时流入和流出相等处。

2) 活性代谢物（及效应）的达峰时间

对 $B(t)$ 求导并令其为零：

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

此公式量纲一致（两边均为时间），对 $k_f\neq k_e$ 有效。

特殊情况 $k_f=k_e$。 使用洛必达法则或求解退化 ODE，可得

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

即当生成和消除速率相等时，峰值出现在一个平均寿命处，且 $T_{A\to B}=T_{B\to C}$。

3) 反演达峰时间公式以求 $k_f$

我们的目标是从已知的 $k_e$ 和观测的 $t_\text{peak}$ 计算 $k_f$（进而求 $T_{A\to B}$）。定义无量纲比值

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

此超越方程有使用 Lambert $W$ 函数的闭式解。推导如下：

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

根据 Lambert $W$ 函数的定义（$W(z)\,e^{W(z)}=z$），我们识别出

\[ -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$，然后转换为所需的半衰期

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

什么是 Lambert $W$？ 它是 $x\mapsto x e^{x}$ 的多值反函数。许多科学库将其实现为 lambertw。在我们的情况下，参数 $z=-a e^{-a}$ 位于 $(-e^{-1},0)$ 区间内，此处存在两个实分支：主分支 $W_0$ 和下分支 $W_{-1}$。

分支选择（为什么有两个可能的 $W$，该用哪个？）

映射 $r\mapsto a=\frac{\ln r}{r-1}$ 在 $r>0, r\neq 1$ 时单调递减，极限为

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

因此：

若 $r<1$（生成慢于消除，$k_f<k_e$），则 $a>1$（峰值出现在 $B$ 的一个平均寿命之后）。
若 $r>1$（生成快于消除，$k_f>k_e$），则 $a<1$（峰值出现在 $B$ 的一个平均寿命之前）。
若 $r=1$，则 $a=1$（即上述特殊情况）。

对于 $z=-a e^{-a}\in(-e^{-1},0)$，两个实分支的行为：

$W_0(z)\in(-1,0)$ $\Rightarrow$ 得到 $r\in(0,1)$。
$W_{-1}(z)\in(-\infty,-1]$ $\Rightarrow$ 得到 $r\in[1,\infty)$。

规则（以可观测量表示）：

\[ b=\begin{cases} 0, & a\gt 1 \quad (\text{使用 } W_0 \Rightarrow k_f \lt k_e \Rightarrow T_{A\to B} \gt T_{B\to C}),\\ -1,& a\lt 1 \quad (\text{使用 } W_{-1} \Rightarrow k_f \gt k_e \Rightarrow T_{A\to B} \lt T_{B\to C}),\\ \text{均可},& a=1 \quad (\text{得到 } k_f=k_e). \end{cases} \]

最终闭式解（可直接使用）

令

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

分支 $b$ 通过上述规则选择（将 $a$ 与 1 比较）。

量纲检查。 $a$ 和 $z$ 是无量纲的。分子具有时间单位，$W_b(z)$ 无量纲，因此 $T_{A\to B}$ 具有时间单位（符合要求）。

实用指南

所需输入。
- $T_{B\to C}$：活性代谢物 $B$ 的终末半衰期（来自浓度-时间数据）。
- $t_\text{peak}$：最大效应时间（如有 $B$ 的最大值时间也可）。若使用效应，确保效应直接与 $B$ 成正比且无延迟（无滞回）。
计算无量纲 $a$。
\[ a=\frac{(\ln 2)\,t_\text{peak}}{T_{B\to C}}. \]
选择分支。
- 若 $a>1$：使用 $W_0$（生成慢于消除）。
- 若 $a<1$：使用 $W_{-1}$（生成快于消除）。
- 若 $a=1$：$T_{A\to B}=T_{B\to C}$。
计算 $T_{A\to B}$。
- 计算 $z=-a e^{-a}$，然后在数学库中计算 $W_b(z)$；最后计算 $T_{A\to B}=-(\ln 2)\,t_\text{peak}/W_b(z)$。
合理性检查。
- 若 $t_\text{peak}$ 相对于 $T_{B\to C}$ 非常早（$a$ 小），应得到 $T_{A\to B}\ll T_{B\to C}$。
- 若 $t_\text{peak}$ 非常晚，应得到 $T_{A\to B}\gg T_{B\to C}$。

边界情况与假设

等速率 $k_f=k_e$：如上处理，$t_\text{peak}=1/k_e$。
模型适用范围。 推导假设：
- 生成和消除均为一级（线性）动力学。
- $B$ 为单一充分混合的房室（无分布延迟）。
- PD 效应与 $B$ 成正比且无延迟（无效应房室/耐受性）。违反这些假设（如可饱和代谢、多房室分布、间接响应）会偏移 $t_\text{peak}$ 并破坏映射关系。
测量噪声。 稀疏采样估计的 $t_\text{peak}$ 可能有偏差。对 $B$（或效应）的时间进程进行平滑或模型拟合以估计真实达峰时间。

计算吸收半衰期的 Python 代码

compute_formation_half_life.py

#!/usr/bin/env python3
# 使用 Lambert W（NumPy + SciPy）从 T_{B->C} 和 t_peak 计算 T_{A->B}
import numpy as np
from scipy.special import lambertw

def formation_half_life(T_BC, t_peak):
    """
    返回 A->B->C 链中的 T_{A->B}（生成半衰期），给定：
      - T_BC   : B->C 的半衰期（与所需输出相同的时间单位）
      - t_peak : B 的达峰时间（若 E ∝ B，也是效应达峰时间）

    闭式解：
        T_AB = - (ln 2) * t_peak / W_b( -a * e^{-a} ),
      其中 a = (ln 2) * t_peak / T_BC
    分支规则（实数解）：
        若 a > 1 -> 使用 W_0   （生成慢于消除）
        若 a < 1 -> 使用 W_{-1}（生成快于消除）
        若 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 和 t_peak 必须为正数。")

    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)

    # 特殊情况：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 -> 主分支 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 -> 下分支 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

    # 若输入为标量则返回标量
    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):
    """便捷函数：返回 k_f = ln(2) / T_{A->B}（相同的广播行为）。"""
    T_AB = formation_half_life(T_BC, t_peak)
    return np.log(2.0) / T_AB


if __name__ == "__main__":
    # 示例
    # 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) 生成较慢 (a > 1) -> T_AB > T_BC
    print("Example a>1:", formation_half_life(8.0, 20.0))

    # 3) 生成较快 (a < 1) -> T_AB < T_BC
    print("Example a<1:", formation_half_life(12.0, 2.0))

    # 4) 向量化输入
    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
# 使用 Lambert W（NumPy + SciPy）从 T_{B->C} 和 t_peak 计算 T_{A->B}
import numpy as np
from scipy.special import lambertw

def formation_half_life(T_BC, t_peak):
    """
    返回 A->B->C 链中的 T_{A->B}（生成半衰期），给定：
      - T_BC   : B->C 的半衰期（与所需输出相同的时间单位）
      - t_peak : B 的达峰时间（若 E ∝ B，也是效应达峰时间）

    闭式解：
        T_AB = - (ln 2) * t_peak / W_b( -a * e^{-a} ),
      其中 a = (ln 2) * t_peak / T_BC
    分支规则（实数解）：
        若 a > 1 -> 使用 W_0   （生成慢于消除）
        若 a < 1 -> 使用 W_{-1}（生成快于消除）
        若 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 和 t_peak 必须为正数。")

    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)

    # 特殊情况：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 -> 主分支 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 -> 下分支 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

    # 若输入为标量则返回标量
    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):
    """便捷函数：返回 k_f = ln(2) / T_{A->B}（相同的广播行为）。"""
    T_AB = formation_half_life(T_BC, t_peak)
    return np.log(2.0) / T_AB


if __name__ == "__main__":
    # 示例
    # 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) 生成较慢 (a > 1) -> T_AB > T_BC
    print("Example a>1:", formation_half_life(8.0, 20.0))

    # 3) 生成较快 (a < 1) -> T_AB < T_BC
    print("Example a<1:", formation_half_life(12.0, 2.0))

    # 4) 向量化输入
    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))

应用示例：喹硫平的吸收

本例使用口服喹硫平（代谢物 A）到喹硫平（代谢物 B）的吸收过程。

根据 NCBI 的实验药理学数据：

约 1 至 2 小时 => 我们使用平均值 1.5h 作为 $T_{A\to B}$
喹硫平的报告半衰期（t1/2）约为 7 小时。 作为 $T_{B\to C}$

从下图可以看出，文献报道的达峰血浆浓度与基于计算吸收半衰期的解析模拟一致。

Quetiapine pharmacokinetics.svg

输出：

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    # 喹硫平消除半衰期（B -> 消除）
t_peak_hours = 1.5  # Tmax（达峰血浆浓度时间）

def compute_T_AB(T_BC, t_peak):
    """从 T_BC 和 t_peak 返回吸收半衰期 T_AB。
       若 SciPy 可用则使用 Lambert W；否则对 a=ln r/(r-1) 使用二分法。"""
    if T_BC <= 0 or t_peak <= 0:
        raise ValueError("T_BC 和 t_peak 必须为正数。")
    a = ln2 * t_peak / T_BC
    if abs(a - 1.0) <= 1e-12:
        return T_BC  # 等速率情况

    # 使用 Lambert W 计算（需要 SciPy）
    z = -a * math.exp(-a)
    branch = 0 if a > 1 else -1
    W = lambertw(z, branch)
    return float(-(ln2 * t_peak) / W.real)

# 计算半衰期和速率
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")

# 时间网格和归一化量（A(0)=1）
t = np.linspace(0.0, 24.0, 1000)  # 小时
A = np.exp(-kf * t)
B = (kf / (ke - kf)) * (np.exp(-kf * t) - np.exp(-ke * t))  # 归一化 Bateman 形式

# 绘图
plt.figure(figsize=(8,5))
plt.plot(t, A, label="A(t): 口服喹硫平（归一化）")
plt.plot(t, B, label="B(t): 喹硫平（归一化）")
plt.axvline(t_peak_hours, linestyle="--", label=f"文献 Tmax ≈ {t_peak_hours:.2f} h")
plt.xlabel("时间（小时）")
plt.ylabel("归一化量 / 浓度（a.u.）")
plt.title("喹硫平 — 推导的吸收半衰期与归一化 PK 曲线")
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    # 喹硫平消除半衰期（B -> 消除）
t_peak_hours = 1.5  # Tmax（达峰血浆浓度时间）

def compute_T_AB(T_BC, t_peak):
    """从 T_BC 和 t_peak 返回吸收半衰期 T_AB。
       若 SciPy 可用则使用 Lambert W；否则对 a=ln r/(r-1) 使用二分法。"""
    if T_BC <= 0 or t_peak <= 0:
        raise ValueError("T_BC 和 t_peak 必须为正数。")
    a = ln2 * t_peak / T_BC
    if abs(a - 1.0) <= 1e-12:
        return T_BC  # 等速率情况

    # 使用 Lambert W 计算（需要 SciPy）
    z = -a * math.exp(-a)
    branch = 0 if a > 1 else -1
    W = lambertw(z, branch)
    return float(-(ln2 * t_peak) / W.real)

# 计算半衰期和速率
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")

# 时间网格和归一化量（A(0)=1）
t = np.linspace(0.0, 24.0, 1000)  # 小时
A = np.exp(-kf * t)
B = (kf / (ke - kf)) * (np.exp(-kf * t) - np.exp(-ke * t))  # 归一化 Bateman 形式

# 绘图
plt.figure(figsize=(8,5))
plt.plot(t, A, label="A(t): 口服喹硫平（归一化）")
plt.plot(t, B, label="B(t): 喹硫平（归一化）")
plt.axvline(t_peak_hours, linestyle="--", label=f"文献 Tmax ≈ {t_peak_hours:.2f} h")
plt.xlabel("时间（小时）")
plt.ylabel("归一化量 / 浓度（a.u.）")
plt.title("喹硫平 — 推导的吸收半衰期与归一化 PK 曲线")
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

药代动力学：从达峰时间计算吸收半衰期的解析解