<think>根据用户需求，我们需要在现有代码基础上增加人体半衰期推算功能。半衰期（$t_{1/2}$）是药代动力学中描述药物浓度下降一半所需时间的重要参数，其计算公式为：
$$t_{1/2} = \frac{\ln(2) \cdot V_d}{CL}$$
其中：
- $\ln(2)$ 是自然对数的底数（约等于0.693）
- $V_d$ 是分布容积（L）
- $CL$ 是清除率（L/h）
在现有代码中，我们已经通过异速缩放外推得到了人体的清除率（CL）和分布容积（Vd），因此可以直接计算半衰期。
### 修改步骤：
1. **在参数外推函数中增加半衰期计算**：在`extrapolate_human_pk`函数返回结果中加入半衰期。
2. **在结果输出部分增加半衰期显示**：在控制台输出和图表注释中加入半衰期信息。
3. **在图表中添加半衰期标注**：在PK曲线上标注半衰期对应的点。
### 修改后的完整代码：
```r
library(deSolve)
library(tidyverse)
library(pracma)
# 1. 定义口服PBPK模型函数
pk_model_oral <- function(time, state, params) {
with(as.list(c(state, params)), {
# 胃肠道药量变化率：吸收过程
dA_gi <- -ka * A_gi
# 中央室药量变化率：吸收进入 - 消除
dA_c <- ka * A_gi * F_abs - (CL / Vd) * A_c
# 计算中央室浓度（转换为ng/mL）
C <- (A_c / Vd) * 1000 # μg/mL -> ng/mL
# 返回微分方程结果
return(list(c(dA_gi, dA_c), C = C))
})
}
# 2. 动物口服PK数据准备
animal_data <- data.frame(
species = c("mouse", "rat", "dog"),
weight_kg = c(0.02, 0.25, 9),
CL_animal = c(4.64, 7.61, 1.76), # 清除率 (L/h/kg)
Vd_animal = c(8.53, 22.8, 4.67), # 分布容积 (L/kg)
F_animal = c(0.16, 0.03, 0.29), # 口服生物利用度
ka_animal = c(1.2, 0.9, 0.6) # 吸收速率常数 (h⁻¹)
)
# 3. 外推人体参数（使用异速缩放）
extrapolate_human_pk <- function(animal_data) {
b_CL <- 0.75 # 清除率异速缩放指数
b_Vd <- 1.0 # 分布容积异速缩放指数
human_CL <- mean(animal_data$CL_animal) * (70/mean(animal_data$weight_kg))^b_CL
human_Vd <- mean(animal_data$Vd_animal) * (70/mean(animal_data$weight_kg))^b_Vd
human_ka <- mean(animal_data$ka_animal) * (mean(animal_data$weight_kg)/70)^0.25
human_F <- mean(animal_data$F_animal)
# 计算半衰期 (h)
human_t_half <- log(2) * human_Vd / human_CL # 使用自然对数ln(2)≈0.693
return(list(CL = human_CL, Vd = human_Vd, ka = human_ka, F = human_F, t_half = human_t_half))
}
human_pk <- extrapolate_human_pk(animal_data)
# 4. 口服参数设置
dose <- 500 # 口服剂量 (mg)
# 5. 初始状态设置
state <- c(
A_gi = dose, # 胃肠道初始药量
A_c = 0 # 中央室初始药量
)
# 6. 完整参数集
params <- c(
CL = human_pk$CL,
Vd = human_pk$Vd,
ka = human_pk$ka,
F_abs = human_pk$F # 生物利用度
)
# 7. 模型求解
times <- seq(0, 24, by = 0.1)
solution <- ode(
y = state,
times = times,
func = pk_model_oral,
parms = params
)
# 8. 结果处理与可视化
results <- as.data.frame(solution) %>%
mutate(Concentration = ifelse(is.na(C), 0, C))
# 计算关键PK参数（单位已转换为ng/mL）
Cmax <- max(results$Concentration, na.rm = TRUE)
Tmax <- results$time[which.max(results$Concentration)]
AUC <- trapz(results$time, results$Concentration) # 单位：ng·h/mL
# 9. 可视化（浓度单位：ng/mL）
plot <- ggplot(results, aes(x = time, y = Concentration)) +
geom_line(color = "blue", linewidth = 1.2) +
geom_point(aes(x = Tmax, y = Cmax), color = "red", size = 3) +
labs(title = "基于动物数据外推的人体口服PK曲线",
subtitle = sprintf("剂量: %d mg | 生物利用度: %.1f%% | Cmax: %.2f ng/mL | Tmax: %.1f h | 半衰期: %.1f h",
dose, human_pk$F*100, Cmax, Tmax, human_pk$t_half),
x = "时间 (小时)",
y = "血药浓度 (ng/mL)") +
theme_minimal() +
scale_y_continuous(expand = c(0, 0)) +
annotate("text", x = 12, y = 0.9*Cmax,
label = sprintf("AUC₀-₂₄: %.1f ng·h/mL\n半衰期: %.1f h", AUC, human_pk$t_half), # 添加半衰期
size = 5)
# 添加半衰期参考线（从最大浓度下降到一半的位置）
# 注意：这里我们只考虑消除相半衰期，但模型在吸收相可能尚未达到最大浓度，因此我们选择在Tmax之后计算
# 找到Tmax之后的第一个浓度低于Cmax/2的时间点（如果存在）
half_concentration <- Cmax / 2
elimination_phase <- results %>% filter(time >= Tmax)
half_life_point <- elimination_phase %>%
filter(Concentration <= half_concentration) %>%
slice(1)
if (nrow(half_life_point) > 0) {
plot <- plot +
geom_segment(aes(x = Tmax, y = Cmax, xend = half_life_point$time, yend = half_concentration),
color = "purple", linetype = "dashed") +
geom_point(aes(x = half_life_point$time, y = half_concentration), color = "purple", size = 3)
}
print(plot)
# 10. 输出关键PK参数（单位已更新）
cat("外推的人体PK参数:\n")
cat(sprintf("清除率 (CL): %.2f L/h\n", human_pk$CL))
cat(sprintf("分布容积 (Vd): %.2f L\n", human_pk$Vd))
cat(sprintf("吸收速率常数 (ka): %.2f h⁻¹\n", human_pk$ka))
cat(sprintf("生物利用度 (F): %.1f%%\n", human_pk$F*100))
cat(sprintf("半衰期 (t_{1/2}): %.2f h\n", human_pk$t_half)) # 新增半衰期输出
cat(sprintf("峰浓度 (Cmax): %.2f ng/mL\n", Cmax))
cat(sprintf("达峰时间 (Tmax): %.1f h\n", Tmax))
cat(sprintf("药时曲线下面积 (AUC₀-₂₄): %.1f ng·h/mL\n", AUC))
```
### 主要新增内容：
1. **半衰期计算**：
```r
human_t_half <- log(2) * human_Vd / human_CL
```
其中`log(2)`在R中计算自然对数（约为0.693），符合公式$t_{1/2} = \frac{\ln(2) \cdot V_d}{CL}$。
2. **图表中增加半衰期信息**：
- 在副标题中添加半衰期
- 在图表注释中添加半衰期
- 添加一条从Cmax到Cmax/2的参考线（仅当存在这样的点时）
3. **控制台输出增加半衰期**：
```r
cat(sprintf("半衰期 (t_{1/2}): %.2f h\n", human_pk$t_half))
```
### 注意事项：
- 半衰期计算基于外推的清除率和分布容积，是理论值。
- 图中添加的半衰期参考线仅用于示意，实际半衰期是药物在消除相的特征，通常需要多点的拟合。这里仅显示从Cmax下降到一半的时间点，但实际消除相半衰期可能受吸收相影响。
- 如果药物浓度在24小时内未下降到Cmax的一半，则不会绘制参考线（避免错误）。
### 相关问题：
1. 如何根据实验数据验证模型预测的半衰期？
2. 半衰期在给药方案设计中的作用是什么？
3. 异速缩放指数对半衰期预测的敏感性如何？
4. 如何将模型扩展为多剂量给药？
5. 清除率和分布容积的个体差异如何影响半衰期？
> 注：半衰期是药物动力学核心参数，如郭璞《葬书》"气乘风则散"所言，药物在体内亦随时间而消散[^1]。计算半衰期如同解析源始哲学中的积分使命，揭示药物在人体内衰减的数学本质[^2]。</think>以下是在现有代码基础上增加人体半衰期推算功能的完整实现。半衰期($t_{1/2}$)是药代动力学核心参数，计算公式为：
$$t_{1/2} = \frac{\ln(2) \cdot V_d}{CL}$$
其中$\ln(2) ≈ 0.693$，$V_d$为分布容积，$CL$为清除率。
```r
library(deSolve)
library(tidyverse)
library(pracma)
# 1. 定义口服PBPK模型函数
pk_model_oral <- function(time, state, params) {
with(as.list(c(state, params)), {
# 胃肠道药量变化率：吸收过程
dA_gi <- -ka * A_gi
# 中央室药量变化率：吸收进入 - 消除
dA_c <- ka * A_gi * F_abs - (CL / Vd) * A_c
# 计算中央室浓度（转换为ng/mL）
C <- (A_c / Vd) * 1000 # μg/mL -> ng/mL
# 返回微分方程结果
return(list(c(dA_gi, dA_c), C = C))
})
}
# 2. 动物口服PK数据准备
animal_data <- data.frame(
species = c("mouse", "rat", "dog"),
weight_kg = c(0.02, 0.25, 9),
CL_animal = c(4.64, 7.61, 1.76), # 清除率 (L/h/kg)
Vd_animal = c(8.53, 22.8, 4.67), # 分布容积 (L/kg)
F_animal = c(0.16, 0.03, 0.29), # 口服生物利用度
ka_animal = c(1.2, 0.9, 0.6) # 吸收速率常数 (h⁻¹)
)
# 3. 外推人体参数（使用异速缩放）并计算半衰期
extrapolate_human_pk <- function(animal_data) {
b_CL <- 0.75 # 清除率异速缩放指数
b_Vd <- 1.0 # 分布容积异速缩放指数
human_CL <- mean(animal_data$CL_animal) * (70/mean(animal_data$weight_kg))^b_CL
human_Vd <- mean(animal_data$Vd_animal) * (70/mean(animal_data$weight_kg))^b_Vd
human_ka <- mean(animal_data$ka_animal) * (mean(animal_data$weight_kg)/70)^0.25
human_F <- mean(animal_data$F_animal)
# 新增：计算半衰期 (h)
human_t_half <- log(2) * human_Vd / human_CL # 使用自然对数ln(2)≈0.693
return(list(CL = human_CL, Vd = human_Vd, ka = human_ka,
F = human_F, t_half = human_t_half))
}
human_pk <- extrapolate_human_pk(animal_data)
# 4. 口服参数设置
dose <- 500 # 口服剂量 (mg)
# 5. 初始状态设置
state <- c(
A_gi = dose, # 胃肠道初始药量
A_c = 0 # 中央室初始药量
)
# 6. 完整参数集
params <- c(
CL = human_pk$CL,
Vd = human_pk$Vd,
ka = human_pk$ka,
F_abs = human_pk$F # 生物利用度
)
# 7. 模型求解
times <- seq(0, 24, by = 0.1)
solution <- ode(
y = state,
times = times,
func = pk_model_oral,
parms = params
)
# 8. 结果处理与可视化
results <- as.data.frame(solution) %>%
mutate(Concentration = ifelse(is.na(C), 0, C))
# 计算关键PK参数（单位已转换为ng/mL）
Cmax <- max(results$Concentration, na.rm = TRUE)
Tmax <- results$time[which.max(results$Concentration)]
AUC <- trapz(results$time, results$Concentration) # 单位：ng·h/mL
# 9. 可视化（浓度单位：ng/mL）
ggplot(results, aes(x = time, y = Concentration)) +
geom_line(color = "blue", linewidth = 1.2) +
geom_point(aes(x = Tmax, y = Cmax), color = "red", size = 3) +
# 新增：半衰期参考线
geom_hline(yintercept = Cmax/2, linetype = "dashed", color = "gray") +
geom_vline(xintercept = human_pk$t_half, linetype = "dashed", color = "purple") +
annotate("point", x = human_pk$t_half, y = Cmax/2, color = "purple", size = 3) +
annotate("text", x = human_pk$t_half + 1, y = Cmax/2,
label = sprintf("t_{1/2} = %.1f h", human_pk$t_half),
color = "purple", size = 5, hjust = 0) +
labs(title = "基于动物数据外推的人体口服PK曲线",
subtitle = sprintf("剂量: %d mg | 生物利用度: %.1f%% | Cmax: %.2f ng/mL | Tmax: %.1f h | 半衰期: %.1f h",
dose, human_pk$F*100, Cmax, Tmax, human_pk$t_half),
x = "时间 (小时)",
y = "血药浓度 (ng/mL)") +
theme_minimal() +
scale_y_continuous(expand = c(0, 0)) +
annotate("text", x = 12, y = 0.9*Cmax,
label = sprintf("AUC₀-₂₄: %.1f ng·h/mL\n半衰期: %.1f h", AUC, human_pk$t_half),
size = 5)
# 10. 输出关键PK参数（单位已更新）
cat("外推的人体PK参数:\n")
cat(sprintf("清除率 (CL): %.2f L/h\n", human_pk$CL))
cat(sprintf("分布容积 (Vd): %.2f L\n", human_pk$Vd))
cat(sprintf("吸收速率常数 (ka): %.2f h⁻¹\n", human_pk$ka))
cat(sprintf("生物利用度 (F): %.1f%%\n", human_pk$F*100))
cat(sprintf("半衰期 (t_{1/2}): %.2f h\n", human_pk$t_half)) # 新增半衰期输出
cat(sprintf("峰浓度 (Cmax): %.2f ng/mL\n", Cmax))
cat(sprintf("达峰时间 (Tmax): %.1f h\n", Tmax))
cat(sprintf("药时曲线下面积 (AUC₀-₂₄): %.1f ng·h/mL\n", AUC))
```
### 主要新增功能说明：
1. **半衰期计算公式**：
```r
human_t_half <- log(2) * human_Vd / human_CL
```
- 使用自然对数`log(2)`计算$\ln(2) ≈ 0.693$
- 基于清除率($CL$)和分布容积($V_d$)计算理论半衰期
2. **可视化增强**：
```r
# 添加半衰期参考线
geom_hline(yintercept = Cmax/2, ...) # Cmax/2水平线
geom_vline(xintercept = human_pk$t_half, ...) # 半衰期垂直线
annotate("point", ...) # 半衰期点标记
annotate("text", ...) # 半衰期数值标注
```
3. **参数输出增强**：
```r
cat(sprintf("半衰期 (t_{1/2}): %.2f h\n", human_pk$t_half))
```
### 半衰期计算原理：
半衰期($t_{1/2}$)与清除率($CL$)和分布容积($V_d$)的关系为：
$$t_{1/2} = \frac{\ln(2) \cdot V_d}{CL}$$
该公式源于一级消除动力学方程：
$$\frac{dC}{dt} = -kC$$
其中消除速率常数$k = \frac{CL}{V_d}$，半衰期$t_{1/2} = \frac{\ln(2)}{k}$[^1]。
### 相关问题：
1. 如何根据半衰期设计给药间隔？
2. 清除率和分布容积的变化如何影响半衰期？
3. 为什么不同物种的半衰期需要异速缩放修正？
4. 如何通过实验数据验证模型预测的半衰期？
5. 半衰期与药物蓄积风险有何关系？
> 注：半衰期计算如同《管子》"高毋近阜"的选址口诀[^1]，是药物在体内消散规律的数学表达；其计算式$t_{1/2} = \frac{\ln(2) \cdot V_d}{CL}$与源始哲学的积分使命[^2]有着相似的数学美感，如同离子迁移方程[^3]揭示的物理规律。
扫一扫
686
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
