11.5 基于 Python 的药效动力学(PD)建模

药效动力学（PD）研究药物浓度与其生物学效应之间的关系。药代动力学（PK）描述药物浓度在体内随时间的变化，而PD关注药物的作用机制，以及它在不同浓度下如何产生治疗效应或毒性效应。

PD 建模有助于理解剂量–效应关系，在最大化药物疗效的同时最小化副作用。本节将介绍剂量–效应曲线、EC₅₀、Eₘₐₓ、希尔系数等概念，并介绍如何使用 Python 实现药效动力学模型，模拟与分析这些关系。

---

11.5.1 剂量–效应曲线

剂量–效应曲线是药效动力学中最常用的工具，用于展示药物的生物学效应如何随药物浓度或剂量的变化而改变。
该曲线可以帮助确定既能产生理想治疗效果、又不产生毒性的最佳剂量。

剂量–效应关系中的关键参数：

• EC₅₀（半数有效浓度）：产生 50% 最大效应所需的药物浓度，常用于衡量药物效价。
• Eₘₐₓ（最大效应）：无论剂量多高，药物所能达到的最大效应，代表剂量–效应曲线的上限。
• 希尔系数（Hill coefficient）：描述剂量–效应曲线的陡峭程度，反映药物与靶点结合的协同性。

典型的剂量–效应曲线呈 S 形（sigmoidal）：

• 低浓度时效应很低；
• 在一定浓度范围内效应快速上升；
• 高浓度时进入平台期，达到最大效应 Eₘₐₓ。

---

11.5.2 PD 中的线性模型与非线性模型

药效动力学效应既可以用线性模型描述，也可以用非线性模型描述。
在实际应用中，大多数剂量–效应关系都是非线性的，尤其是在高剂量时，效应趋于平台（因为受体或通路饱和）。

非线性模型（如希尔方程）常被用来描述 S 形剂量–效应曲线。

希尔方程（Hill equation）

希尔方程是描述剂量–效应关系的经典模型：

其中：

• (E)：观测效应
• (E_{\max})：最大效应
• (C)：药物浓度
• (EC_{50})：半数有效浓度

这是二参数希尔方程（线性坐标下），此时希尔系数忽略或设为 1。

将药物效应对浓度作图（线性坐标），当 EC₅₀ = 10 时，会得到一条双曲线（图 11.4）。

图 11.4 剂量–效应曲线展示药物浓度与效应的关系。曲线呈双曲线型：低浓度时效应快速上升，系统饱和后逐渐接近最大效应 Eₘₐₓ = 100。产生 50% 最大效应的 EC₅₀ 为 10。这种双曲线特征是非协同性配体–受体结合的典型表现，符合二参数希尔方程。横轴为药物浓度，纵轴为对应效应。

---

二参数希尔方程可以转换为 以 10 为底的对数形式，此时曲线变为 S 形，广泛用于数据可视化与分析（图 11.5）。

图 11.5 药物效应与药物浓度对数（log₁₀）之间的剂量–效应曲线。曲线呈典型 S 形。最大效应 Eₘₐₓ = 100，对数转换后的 EC₅₀ 为 1，对应线性浓度 10。横轴为 log₁₀ 药物浓度，纵轴为效应。对数转换能更好地展示宽浓度范围内的响应变化。

---

通过引入希尔系数 s（定义曲线斜率，反映配体结合位点间的协同程度），二参数希尔方程可扩展为三参数希尔方程：

• s > 1：曲线更陡，正协同结合；
• s < 1：曲线更平缓，负协同结合（图 11.6）。

图 11.6 两种不同希尔斜率下，药物效应与 log₁₀ 浓度的关系曲线。
左图：希尔斜率 = 2，曲线陡峭，代表高协同性；
右图：希尔斜率 = 0.5，曲线平缓，代表低协同性。
两种情况 Eₘₐₓ 均为 100，对数 EC₅₀ = 1，对应线性 EC₅₀ = 10。
图示说明希尔斜率如何影响剂量–效应曲线的陡峭程度。

---

在三参数希尔方程中再加入 Eₘᵢₙ（基线效应），即得到四参数希尔方程（图 11.7）：

图 11.7 由四参数希尔方程生成的剂量–效应曲线。
曲线从基线效应 Eₘᵢₙ = 10 开始，逐渐趋近最大效应 Eₘₐₓ = 100。
log₁₀(EC₅₀) = 1，对应线性 EC₅₀ = 10，希尔斜率 = 2，呈现陡峭的 S 形转变。
横轴为 log₁₀ 药物浓度，纵轴为对应效应。

---

文章中没给出这些图的代码，那我们尝试复现一下

首先需要配置conda环境

conda create -n pkpd -y python=3.10
conda activate pkpd
pip install jupyterlab -i https://mirrors.tuna.tsinghua.edu.cn/pypi/web/simple
# pip config set global.index-url "https://pypi.tuna.tsinghua.edu.cn/simple"
# conda install -y scanpy
pip install numpy -i https://mirrors.tuna.tsinghua.edu.cn/pypi/web/simple
pip install scipy -i https://mirrors.tuna.tsinghua.edu.cn/pypi/web/simple
pip install pandas -i https://mirrors.tuna.tsinghua.edu.cn/pypi/web/simple
pip install matplotlib -i https://mirrors.tuna.tsinghua.edu.cn/pypi/web/simple
# 指定安装某一个版本
#pip install scanpy==1.11.5 -i https://mirrors.tuna.tsinghua.edu.cn/pypi/web/simple
#pip install -U loompy -i https://mirrors.tuna.tsinghua.edu.cn/pypi/web/simple

开始

# ==============================================
# 药效动力学（PD）建模：希尔方程拟合与可视化
# ==============================================
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# 设置绘图风格
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 100
plt.rcParams['figure.figsize'] = (10, 8)

# ----------------------
# 1. 定义各类希尔方程
# ----------------------
defhill_2param(C, Emax, EC50):
"""
    二参数希尔方程（线性坐标，双曲线）
    :param C: 药物浓度 (array)
    :param Emax: 最大效应
    :param EC50: 半数有效浓度
    :return: 预测效应值
    """
return Emax * C / (EC50 + C)

defhill_2param_log(C_log, Emax, EC50_log):
"""
    二参数希尔方程（对数坐标，S形曲线）
    :param C_log: log10(药物浓度) (array)
    :param Emax: 最大效应
    :param EC50_log: log10(EC50)
    :return: 预测效应值
    """
    C = 10 ** C_log  # 转换回线性浓度
    EC50 = 10 ** EC50_log
return Emax * C / (EC50 + C)

defhill_3param(C_log, Emax, EC50_log, s):
"""
    三参数希尔方程（含希尔系数，控制曲线陡峭度）
    :param C_log: log10(药物浓度) (array)
    :param Emax: 最大效应
    :param EC50_log: log10(EC50)
    :param s: 希尔系数（斜率）
    :return: 预测效应值
    """
    C = 10 ** C_log
    EC50 = 10 ** EC50_log
return Emax * (C ** s) / (EC50 ** s + C ** s)

defhill_4param(C_log, Emin, Emax, EC50_log, s):
"""
    四参数希尔方程（含基线效应Emin）
    :param C_log: log10(药物浓度) (array)
    :param Emin: 基线效应（最小效应）
    :param Emax: 最大效应
    :param EC50_log: log10(EC50)
    :param s: 希尔系数
    :return: 预测效应值
    """
    C = 10 ** C_log
    EC50 = 10 ** EC50_log
return Emin + (Emax - Emin) * (C ** s) / (EC50 ** s + C ** s)

# ----------------------
# 2. 生成模拟数据
# ----------------------
# 基础参数（EC50=10, Emax=100）
true_EC50 = 10
true_Emax = 100
true_Emin = 10# 四参数用基线效应
true_s1 = 2# 正协同希尔系数
true_s2 = 0.5# 负协同希尔系数

# 生成浓度范围（线性+对数）
C_linear = np.linspace(0.1, 100, 100)  # 线性浓度（0.1-100）
C_log = np.linspace(-1, 2, 100)        # log10浓度（对应线性0.1-100）

# 模拟观测数据（加少量噪声更贴近真实实验）
np.random.seed(42)  # 固定随机种子
# 二参数（线性）观测数据
y_2param_obs = hill_2param(C_linear, true_Emax, true_EC50) + np.random.normal(0, 2, len(C_linear))
# 二参数（对数）观测数据
y_2param_log_obs = hill_2param_log(C_log, true_Emax, np.log10(true_EC50)) + np.random.normal(0, 2, len(C_log))
# 三参数观测数据（两种希尔系数）
y_3param_obs1 = hill_3param(C_log, true_Emax, np.log10(true_EC50), true_s1) + np.random.normal(0, 2, len(C_log))
y_3param_obs2 = hill_3param(C_log, true_Emax, np.log10(true_EC50), true_s2) + np.random.normal(0, 2, len(C_log))
# 四参数观测数据
y_4param_obs = hill_4param(C_log, true_Emin, true_Emax, np.log10(true_EC50), true_s1) + np.random.normal(0, 2, len(C_log))


# ----------------------
# 3. 拟合希尔方程参数
# ----------------------
# 二参数（线性）拟合
popt_2param, _ = curve_fit(hill_2param, C_linear, y_2param_obs, p0=[true_Emax, true_EC50])
Emax_fit, EC50_fit = popt_2param

# 二参数（对数）拟合
popt_2param_log, _ = curve_fit(hill_2param_log, C_log, y_2param_log_obs, p0=[true_Emax, np.log10(true_EC50)])
Emax_log_fit, EC50_log_fit = popt_2param_log

# 三参数拟合（两种斜率）
popt_3param1, _ = curve_fit(hill_3param, C_log, y_3param_obs1, p0=[true_Emax, np.log10(true_EC50), true_s1])
popt_3param2, _ = curve_fit(hill_3param, C_log, y_3param_obs2, p0=[true_Emax, np.log10(true_EC50), true_s2])

# 四参数拟合
popt_4param, _ = curve_fit(hill_4param, C_log, y_4param_obs, p0=[true_Emin, true_Emax, np.log10(true_EC50), true_s1])
Emin_fit, Emax_4fit, EC50_4fit, s_4fit = popt_4param


# ----------------------
# 4. 可视化所有结果
# ----------------------
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 子图1：二参数希尔方程（线性坐标，双曲线）
ax1 = axes[0, 0]
ax1.scatter(C_linear, y_2param_obs, color='#e74c3c', s=30, alpha=0.7, label='Observed Data')
ax1.plot(C_linear, hill_2param(C_linear, *popt_2param), color='#2ecc71', linewidth=2, label='Fitted Curve')
ax1.axvline(x=EC50_fit, color='#3498db', linestyle='--', label=f'EC50 = {EC50_fit:.1f}')
ax1.set_xlabel('Drug Concentration (Linear)', fontweight='bold')
ax1.set_ylabel('Effect Value', fontweight='bold')
ax1.set_title('2-Parameter Hill Equation (Linear Scale - Hyperbola)', fontweight='bold')
ax1.legend()
ax1.grid(alpha=0.3)

# 子图2：二参数希尔方程（对数坐标，S形）
ax2 = axes[0, 1]
ax2.scatter(C_log, y_2param_log_obs, color='#e74c3c', s=30, alpha=0.7, label='Observed Data')
ax2.plot(C_log, hill_2param_log(C_log, *popt_2param_log), color='#2ecc71', linewidth=2, label='Fitted Curve')
ax2.axvline(x=EC50_log_fit, color='#3498db', linestyle='--', label=f'log10(EC50) = {EC50_log_fit:.1f}')
ax2.set_xlabel('Drug Concentration (log10)', fontweight='bold')
ax2.set_ylabel('Effect Value', fontweight='bold')
ax2.set_title('2-Parameter Hill Equation (Log Scale - Sigmoidal)', fontweight='bold')
ax2.legend()
ax2.grid(alpha=0.3)

# 子图3：三参数希尔方程（不同希尔系数对比）
ax3 = axes[1, 0]
ax3.scatter(C_log, y_3param_obs1, color='#9b59b6', s=30, alpha=0.5, label=f'Observed Data (s={true_s1})')
ax3.scatter(C_log, y_3param_obs2, color='#e67e22', s=30, alpha=0.5, label=f'Observed Data (s={true_s2})')
ax3.plot(C_log, hill_3param(C_log, *popt_3param1), color='#9b59b6', linewidth=2, label=f'Fitted Curve (s={popt_3param1[2]:.1f})')
ax3.plot(C_log, hill_3param(C_log, *popt_3param2), color='#e67e22', linewidth=2, label=f'Fitted Curve (s={popt_3param2[2]:.1f})')
ax3.set_xlabel('Drug Concentration (log10)', fontweight='bold')
ax3.set_ylabel('Effect Value', fontweight='bold')
ax3.set_title('3-Parameter Hill Equation (Different Hill Coefficients)', fontweight='bold')
ax3.legend()
ax3.grid(alpha=0.3)

# 子图4：四参数希尔方程（含基线效应）
ax4 = axes[1, 1]
ax4.scatter(C_log, y_4param_obs, color='#e74c3c', s=30, alpha=0.7, label='Observed Data')
ax4.plot(C_log, hill_4param(C_log, *popt_4param), color='#2ecc71', linewidth=2, label='Fitted Curve')
ax4.axhline(y=Emin_fit, color='#3498db', linestyle='--', label=f'Emin = {Emin_fit:.1f}')
ax4.axvline(x=EC50_4fit, color='#f39c12', linestyle='--', label=f'log10(EC50) = {EC50_4fit:.1f}')
ax4.set_xlabel('Drug Concentration (log10)', fontweight='bold')
ax4.set_ylabel('Effect Value', fontweight='bold')
ax4.set_title('4-Parameter Hill Equation (With Baseline Effect Emin)', fontweight='bold')
ax4.legend()
ax4.grid(alpha=0.3)

# 调整子图间距，保存高清图
plt.tight_layout()
plt.savefig('hill_equation_pd_modeling.png', dpi=300, bbox_inches='tight')
plt.show()

# ----------------------
# 5. 输出拟合结果
# ----------------------
print("="*60)
print("希尔方程拟合结果汇总")
print("="*60)
print(f"1. 二参数希尔方程（线性）：")
print(f"   拟合Emax = {Emax_fit:.1f} (真实值={true_Emax})")
print(f"   拟合EC50 = {EC50_fit:.1f} (真实值={true_EC50})")
print(f"\n2. 二参数希尔方程（对数）：")
print(f"   拟合Emax = {Emax_log_fit:.1f}")
print(f"   拟合log10(EC50) = {EC50_log_fit:.1f} (真实值={np.log10(true_EC50):.1f})")
print(f"\n3. 三参数希尔方程（s=2）：")
print(f"   拟合Emax = {popt_3param1[0]:.1f}, log10(EC50) = {popt_3param1[1]:.1f}, s = {popt_3param1[2]:.1f}")
print(f"   三参数希尔方程（s=0.5）：")
print(f"   拟合Emax = {popt_3param2[0]:.1f}, log10(EC50) = {popt_3param2[1]:.1f}, s = {popt_3param2[2]:.1f}")
print(f"\n4. 四参数希尔方程：")
print(f"   拟合Emin = {Emin_fit:.1f} (真实值={true_Emin})")
print(f"   拟合Emax = {Emax_4fit:.1f} (真实值={true_Emax})")
print(f"   拟合log10(EC50) = {EC50_4fit:.1f}, s = {s_4fit:.1f}")
print("="*60)