# 加载必要的包
library(deSolve) # 用于求解微分方程
library(httk) # 获取人体生理参数
library(tidyverse) # 数据处理和绘图
library(minpack.lm) # 非线性最小二乘拟合
# 定义化合物参数
compound_params <- list(
MW = 581.29, # 分子量 (g/mol)
logP = 1.59, # 油水分配系数
pKa = 8.8, # 碱基解离常数
Clint_human = 45.55, # 人体肝微粒体固有清除率 (μL/min/mg protein)
fup_human = 0.0818, # 血浆游离药物分数 (8.18%)
Rbp_human = 0.954, # 全血/血浆药物浓度比
Papp = 0.45 # Caco-2渗透性 (10-6 cm/s)
)
# 溶解度数据 (mg/mL)
solubility <- data.frame(
pH = c(1.4, 2.5, 4.0, 5.0, 6.5, 7.4),
value = c(0.08557, 0.06794, 0.07188, 0.07629, 0.08165, 0.07515)
)
# 1. 估算人体吸收速率常数 (Ka)
calculate_ka <- function(Papp) {
# 经验公式: Ka = 0.1 * Papp * 3600 * 100 / 10000
ka <- 0.1 * Papp * 3600 * 100 / 10000
return(ka)
}
Ka_human <- calculate_ka(compound_params$Papp)
cat(sprintf("估算的人体吸收速率常数 Ka = %.3f 1/h\n", Ka_human))
# 2. 计算pH依赖性溶解度
fit_solubility <- function(sol_data) {
# Henderson-Hasselbalch方程
hh_model <- function(pH, S0, pKa) {
S0 * (1 + 10^(pKa - pH))
}
# 非线性拟合
fit <- nlsLM(value ~ hh_model(pH, S0, pKa),
data = sol_data,
start = list(S0 = 0.08, pKa = 8.8))
coefs <- coef(fit)
return(list(S0 = coefs[["S0"]], pKa = coefs[["pKa"]]))
}
sol_fit <- fit_solubility(solubility)
cat(sprintf("拟合的固有溶解度 S0 = %.5f mg/mL, pKa = %.2f\n", sol_fit$S0, sol_fit$pKa))
# 3. 计算组织分配系数 (Kp)
calculate_kp <- function(logP, pKa) {
# 组织成分数据 (简化版)
tissue_data <- data.frame(
Tissue = c("adipose", "bone", "brain", "gut", "heart", "kidney", "liver", "lung", "muscle", "skin"),
f_water = c(0.15, 0.40, 0.78, 0.77, 0.76, 0.77, 0.75, 0.80, 0.76, 0.70),
f_lipids = c(0.80, 0.25, 0.10, 0.05, 0.03, 0.03, 0.05, 0.03, 0.05, 0.10),
f_proteins = c(0.05, 0.35, 0.12, 0.18, 0.21, 0.20, 0.20, 0.17, 0.19, 0.20)
)
# 计算各组织Kp
kp_values <- sapply(1:nrow(tissue_data), function(i) {
f_water <- tissue_data$f_water[i]
f_lipid <- tissue_data$f_lipids[i]
f_protein <- tissue_data$f_proteins[i]
P <- 10^logP
Kp_neutral <- f_water + 0.3 * f_protein + f_lipid * P
pH <- 7.4
fraction_ionized <- 1 / (1 + 10^(pKa - pH))
Kp_ion <- f_water + 0.3 * f_protein + f_lipid * (P / 100)
Kp <- fraction_ionized * Kp_ion + (1 - fraction_ionized) * Kp_neutral
# 特殊组织调整
if (tissue_data$Tissue[i] == "liver") Kp <- Kp * 1.2
if (tissue_data$Tissue[i] == "kidney") Kp <- Kp * 1.1
return(Kp)
})
names(kp_values) <- tissue_data$Tissue
return(kp_values)
}
Kp_human <- calculate_kp(compound_params$logP, compound_params$pKa)
cat("预测的人体组织-血浆分配系数 (Kp):\n")
print(Kp_human)
# 4. 准备人体生理参数 (70kg标准人体)
get_physio_params <- function(bw = 70) {
# 器官体积 (L)
organ_volumes <- c(
gut = 1.5,
liver = 1.8,
kidney = 0.3,
adipose = 15.0,
muscle = 29.0,
heart = 0.3,
brain = 1.4,
lung = 0.8
) * bw / 70 # 按体重缩放
# 血流速率 (L/h)
Q_total <- 6.5 * 60 # 心输出量 (L/min) 转换为 L/h
organ_flows <- c(
gut = 0.15 * Q_total,
liver = 0.25 * Q_total,
kidney = 0.2 * Q_total,
adipose = 0.05 * Q_total,
muscle = 0.17 * Q_total,
heart = 0.04 * Q_total,
brain = 0.12 * Q_total,
lung = 0.1 * Q_total
)
# 其他参数
other_params <- list(
BW = bw,
Q_total = Q_total,
V_venous = 3.5 * bw / 70,
V_arterial = 1.5 * bw / 70
)
return(list(
volumes = organ_volumes,
flows = organ_flows,
other = other_params
))
}
physio_human <- get_physio_params(70)
# 5. 构建PBPK模型函数 - 修复导数返回问题
pbpk_model <- function(time, state, parameters) {
# 提取状态变量
A_gi <- state["A_gi"]
A_portal <- state["A_portal"]
A_gut <- state["A_gut"]
A_liver <- state["A_liver"]
A_kidney <- state["A_kidney"]
A_adipose <- state["A_adipose"]
A_muscle <- state["A_muscle"]
A_heart <- state["A_heart"]
A_brain <- state["A_brain"]
A_lung <- state["A_lung"]
A_venous <- state["A_venous"]
A_arterial <- state["A_arterial"]
# 提取参数
Ka <- parameters["Ka"]
fup <- parameters["fup"]
Clint <- parameters["Clint"]
# 组织体积
V_gut <- parameters["V_gut"]
V_liver <- parameters["V_liver"]
V_kidney <- parameters["V_kidney"]
V_adipose <- parameters["V_adipose"]
V_muscle <- parameters["V_muscle"]
V_heart <- parameters["V_heart"]
V_brain <- parameters["V_brain"]
V_lung <- parameters["V_lung"]
V_venous <- parameters["V_venous"]
V_arterial <- parameters["V_arterial"]
# 血流参数
Q_total <- parameters["Q_total"]
Q_gut <- parameters["Q_gut"]
Q_liver <- parameters["Q_liver"]
Q_hepatic_art <- parameters["Q_hepatic_art"]
Q_kidney <- parameters["Q_kidney"]
Q_adipose <- parameters["Q_adipose"]
Q_muscle <- parameters["Q_muscle"]
Q_heart <- parameters["Q_heart"]
Q_brain <- parameters["Q_brain"]
Q_lung <- parameters["Q_lung"]
# Kp值
Kp_gut <- parameters["Kp_gut"]
Kp_liver <- parameters["Kp_liver"]
Kp_kidney <- parameters["Kp_kidney"]
Kp_adipose <- parameters["Kp_adipose"]
Kp_muscle <- parameters["Kp_muscle"]
Kp_heart <- parameters["Kp_heart"]
Kp_brain <- parameters["Kp_brain"]
Kp_lung <- parameters["Kp_lung"]
# 清除率参数
CL_renal <- parameters["CL_renal"]
# 计算组织浓度
C_liver <- A_liver / V_liver
C_gut <- A_gut / V_gut
C_kidney <- A_kidney / V_kidney
C_adipose <- A_adipose / V_adipose
C_muscle <- A_muscle / V_muscle
C_heart <- A_heart / V_heart
C_brain <- A_brain / V_brain
C_lung <- A_lung / V_lung
C_venous <- A_venous / V_venous
C_arterial <- A_arterial / V_arterial
# 门静脉浓度使用动脉浓度近似
C_portal <- C_arterial
# 计算肝脏清除率
Clint_liver <- Clint * 60 * 1e-6 * 45 * V_liver * 1e3
CLint_free <- Clint_liver * fup
CL_hepatic <- Q_liver * (CLint_free / (Q_liver + CLint_free))
# 微分方程系统
dA_gi <- -Ka * A_gi
dA_portal <- Ka * A_gi - Q_gut * C_portal
dA_liver <- Q_gut * C_portal + Q_hepatic_art * C_arterial -
(Q_gut + Q_hepatic_art) * C_liver / Kp_liver -
CL_hepatic * C_liver * fup
dA_gut <- Q_gut * (C_arterial - C_gut / Kp_gut)
dA_kidney <- Q_kidney * (C_arterial - C_kidney / Kp_kidney) -
CL_renal * C_kidney * fup
dA_adipose <- Q_adipose * (C_arterial - C_adipose / Kp_adipose)
dA_muscle <- Q_muscle * (C_arterial - C_muscle / Kp_muscle)
dA_heart <- Q_heart * (C_arterial - C_heart / Kp_heart)
dA_brain <- Q_brain * (C_arterial - C_brain / Kp_brain)
dA_lung <- Q_lung * (C_arterial - C_lung / Kp_lung)
dA_venous <- (Q_liver * C_liver / Kp_liver +
Q_kidney * C_kidney / Kp_kidney +
Q_adipose * C_adipose / Kp_adipose +
Q_muscle * C_muscle / Kp_muscle +
Q_heart * C_heart / Kp_heart +
Q_brain * C_brain / Kp_brain +
Q_lung * C_lung / Kp_lung) -
Q_total * C_venous
dA_arterial <- Q_total * (C_venous - C_arterial)
# 返回导数向量 - 关键修复:确保返回12个导数值
return(list(c(
dA_gi,
dA_portal,
dA_gut,
dA_liver,
dA_kidney,
dA_adipose,
dA_muscle,
dA_heart,
dA_brain,
dA_lung,
dA_venous,
dA_arterial
)))
}
# 6. 参数设置
parameters <- c(
Ka = Ka_human,
fup = compound_params$fup_human,
Clint = compound_params$Clint_human,
# Kp值
Kp_gut = Kp_human["gut"],
Kp_liver = Kp_human["liver"],
Kp_kidney = Kp_human["kidney"],
Kp_adipose = Kp_human["adipose"],
Kp_muscle = Kp_human["muscle"],
Kp_heart = Kp_human["heart"],
Kp_brain = Kp_human["brain"],
Kp_lung = Kp_human["lung"],
# 生理参数
V_gi = 0.25, # 胃肠道体积 (L)
V_gut = physio_human$volumes["gut"],
V_liver = physio_human$volumes["liver"],
V_kidney = physio_human$volumes["kidney"],
V_adipose = physio_human$volumes["adipose"],
V_muscle = physio_human$volumes["muscle"],
V_heart = physio_human$volumes["heart"],
V_brain = physio_human$volumes["brain"],
V_lung = physio_human$volumes["lung"],
V_venous = physio_human$other$V_venous,
V_arterial = physio_human$other$V_arterial,
# 血流参数
Q_total = physio_human$other$Q_total,
Q_gut = physio_human$flows["gut"],
Q_liver = physio_human$flows["liver"],
Q_hepatic_art = physio_human$flows["liver"] - physio_human$flows["gut"],
Q_kidney = physio_human$flows["kidney"],
Q_adipose = physio_human$flows["adipose"],
Q_muscle = physio_human$flows["muscle"],
Q_heart = physio_human$flows["heart"],
Q_brain = physio_human$flows["brain"],
Q_lung = physio_human$flows["lung"],
# 清除率参数
CL_renal = 1.5 # 肾脏清除率 (L/h)
)
# 7. 初始状态 (口服给药 100 mg)
state <- c(
A_gi = 100, # 胃肠道药物量 (mg)
A_portal = 0, # 门静脉药物量
A_gut = 0, # 肠道组织
A_liver = 0, # 肝脏
A_kidney = 0, # 肾脏
A_adipose = 0, # 脂肪组织
A_muscle = 0, # 肌肉组织
A_heart = 0, # 心脏
A_brain = 0, # 脑
A_lung = 0, # 肺
A_venous = 0, # 静脉血池
A_arterial = 0 # 动脉血池
)
# 8. 模拟时间范围 (0-24小时)
times <- seq(0, 24, by = 0.1)
# 9. 运行模型
out <- ode(y = state, times = times, func = pbpk_model, parms = parameters)
# 10. 结果处理
results <- as.data.frame(out) %>%
mutate(Concentration_plasma = A_arterial / (parameters["V_arterial"] * compound_params$Rbp_human))
# 11. 计算药代动力学参数
pk_analysis <- function(conc, time) {
# AUC计算 (梯形法)
auc <- sum(diff(time) * (head(conc, -1) + tail(conc, -1)))/2
# 终末消除速率常数
terminal_phase <- tail(data.frame(time, conc), 20)
fit <- lm(log(conc) ~ time, data = terminal_phase)
kel <- -coef(fit)[2]
t_half <- log(2)/kel
# 清除率 (基于静脉给药AUC)
dose_iv <- 100 # 假设静脉给药剂量
auc_iv <- auc * 1.2 # 假设静脉AUC比口服低
CL <- dose_iv / auc_iv
# 分布容积
Vd <- CL / kel
# 生物利用度
F_oral <- auc / auc_iv * (dose_iv / 100) # 口服剂量100mg
return(list(AUC = auc, t_half = t_half, CL = CL, Vd = Vd, F_abs = F_oral))
}
pk_params <- pk_analysis(results$Concentration_plasma, results$time)
# 12. 结果输出
cat("\n预测的人体药代动力学参数:\n")
cat(sprintf("吸收速率常数 Ka = %.3f 1/h\n", Ka_human))
cat(sprintf("清除率 CL = %.2f L/h\n", pk_params$CL))
cat(sprintf("表观分布容积 Vd = %.2f L\n", pk_params$Vd))
cat(sprintf("消除半衰期 t1/2 = %.2f h\n", pk_params$t_half))
cat(sprintf("口服生物利用度 F = %.1f%%\n", pk_params$F_abs * 100))
# 13. 绘图
ggplot(results, aes(x = time, y = Concentration_plasma)) +
geom_line(color = "blue", linewidth = 1) +
labs(title = "PBPK模型模拟 - 血浆药物浓度时间曲线",
x = "时间 (小时)",
y = "血浆浓度 (mg/L)") +
theme_minimal()
# 14. 组织分布可视化
tissue_concentrations <- results %>%
select(time, A_liver, A_kidney, A_muscle, A_adipose) %>%
pivot_longer(cols = -time, names_to = "Tissue", values_to = "Amount") %>%
mutate(Concentration = case_when(
Tissue == "A_liver" ~ Amount / parameters["V_liver"],
Tissue == "A_kidney" ~ Amount / parameters["V_kidney"],
Tissue == "A_muscle" ~ Amount / parameters["V_muscle"],
Tissue == "A_adipose" ~ Amount / parameters["V_adipose"]
))
ggplot(tissue_concentrations, aes(x = time, y = Concentration, color = Tissue)) +
geom_line(linewidth = 1) +
labs(title = "组织药物浓度时间曲线",
x = "时间 (小时)",
y = "组织浓度 (mg/L)") +
scale_color_manual(values = c("red", "green", "blue", "purple"),
labels = c("肝脏", "肾脏", "肌肉", "脂肪")) +
theme_minimal()
# 加载必要的包
library(deSolve) # 用于求解微分方程
library(httk) # 获取人体生理参数
library(tidyverse) # 数据处理和绘图
library(minpack.lm) # 非线性最小二乘拟合
# 定义化合物参数
compound_params <- list(
MW = 581.29, # 分子量 (g/mol)
logP = 1.59, # 油水分配系数
pKa = 8.8, # 碱基解离常数
Clint_human = 45.55, # 人体肝微粒体固有清除率 (μL/min/mg protein)
fup_human = 0.0818, # 血浆游离药物分数 (8.18%)
Rbp_human = 0.954, # 全血/血浆药物浓度比
Papp = 0.45 # Caco-2渗透性 (10-6 cm/s)
)
# 溶解度数据 (mg/mL)
solubility <- data.frame(
pH = c(1.4, 2.5, 4.0, 5.0, 6.5, 7.4),
value = c(0.08557, 0.06794, 0.07188, 0.07629, 0.08165, 0.07515)
)
# 1. 估算人体吸收速率常数 (Ka)
calculate_ka <- function(Papp) {
# 经验公式: Ka = 0.1 * Papp * 3600 * 100 / 10000
ka <- 0.1 * Papp * 3600 * 100 / 10000
return(ka)
}
Ka_human <- calculate_ka(compound_params$Papp)
cat(sprintf("估算的人体吸收速率常数 Ka = %.3f 1/h\n", Ka_human))
# 2. 计算pH依赖性溶解度
fit_solubility <- function(sol_data) {
# Henderson-Hasselbalch方程
hh_model <- function(pH, S0, pKa) {
S0 * (1 + 10^(pKa - pH))
}
# 非线性拟合
fit <- nlsLM(value ~ hh_model(pH, S0, pKa),
data = sol_data,
start = list(S0 = 0.08, pKa = 8.8))
coefs <- coef(fit)
return(list(S0 = coefs[["S0"]], pKa = coefs[["pKa"]]))
}
sol_fit <- fit_solubility(solubility)
cat(sprintf("拟合的固有溶解度 S0 = %.5f mg/mL, pKa = %.2f\n", sol_fit$S0, sol_fit$pKa))
# 3. 计算组织分配系数 (Kp)
calculate_kp <- function(logP, pKa) {
# 组织成分数据 (简化版)
tissue_data <- data.frame(
Tissue = c("adipose", "bone", "brain", "gut", "heart", "kidney", "liver", "lung", "muscle", "skin"),
f_water = c(0.15, 0.40, 0.78, 0.77, 0.76, 0.77, 0.75, 0.80, 0.76, 0.70),
f_lipids = c(0.80, 0.25, 0.10, 0.05, 0.03, 0.03, 0.05, 0.03, 0.05, 0.10),
f_proteins = c(0.05, 0.35, 0.12, 0.18, 0.21, 0.20, 0.20, 0.17, 0.19, 0.20)
)
# 计算各组织Kp
kp_values <- sapply(1:nrow(tissue_data), function(i) {
f_water <- tissue_data$f_water[i]
f_lipid <- tissue_data$f_lipids[i]
f_protein <- tissue_data$f_proteins[i]
P <- 10^logP
Kp_neutral <- f_water + 0.3 * f_protein + f_lipid * P
pH <- 7.4
fraction_ionized <- 1 / (1 + 10^(pKa - pH))
Kp_ion <- f_water + 0.3 * f_protein + f_lipid * (P / 100)
Kp <- fraction_ionized * Kp_ion + (1 - fraction_ionized) * Kp_neutral
# 特殊组织调整
if (tissue_data$Tissue[i] == "liver") Kp <- Kp * 1.2
if (tissue_data$Tissue[i] == "kidney") Kp <- Kp * 1.1
return(Kp)
})
names(kp_values) <- tissue_data$Tissue
return(kp_values)
}
Kp_human <- calculate_kp(compound_params$logP, compound_params$pKa)
cat("预测的人体组织-血浆分配系数 (Kp):\n")
print(Kp_human)
# 4. 准备人体生理参数 (70kg标准人体)
get_physio_params <- function(bw = 70) {
# 器官体积 (L)
organ_volumes <- c(
gut = 1.5,
liver = 1.8,
kidney = 0.3,
adipose = 15.0,
muscle = 29.0,
heart = 0.3,
brain = 1.4,
lung = 0.8
) * bw / 70 # 按体重缩放
# 血流速率 (L/h)
Q_total <- 6.5 * 60 # 心输出量 (L/min) 转换为 L/h
organ_flows <- c(
gut = 0.15 * Q_total,
liver = 0.25 * Q_total,
kidney = 0.2 * Q_total,
adipose = 0.05 * Q_total,
muscle = 0.17 * Q_total,
heart = 0.04 * Q_total,
brain = 0.12 * Q_total,
lung = 0.1 * Q_total
)
# 其他参数
other_params <- list(
BW = bw,
Q_total = Q_total,
V_venous = 3.5 * bw / 70,
V_arterial = 1.5 * bw / 70
)
return(list(
volumes = organ_volumes,
flows = organ_flows,
other = other_params
))
}
physio_human <- get_physio_params(70)
# 5. 构建PBPK模型函数 - 修复导数返回问题
pbpk_model <- function(time, state, parameters) {
# 提取状态变量
A_gi <- state["A_gi"]
A_portal <- state["A_portal"]
A_gut <- state["A_gut"]
A_liver <- state["A_liver"]
A_kidney <- state["A_kidney"]
A_adipose <- state["A_adipose"]
A_muscle <- state["A_muscle"]
A_heart <- state["A_heart"]
A_brain <- state["A_brain"]
A_lung <- state["A_lung"]
A_venous <- state["A_venous"]
A_arterial <- state["A_arterial"]
# 提取参数
Ka <- parameters["Ka"]
fup <- parameters["fup"]
Clint <- parameters["Clint"]
# 组织体积
V_gut <- parameters["V_gut"]
V_liver <- parameters["V_liver"]
V_kidney <- parameters["V_kidney"]
V_adipose <- parameters["V_adipose"]
V_muscle <- parameters["V_muscle"]
V_heart <- parameters["V_heart"]
V_brain <- parameters["V_brain"]
V_lung <- parameters["V_lung"]
V_venous <- parameters["V_venous"]
V_arterial <- parameters["V_arterial"]
# 血流参数
Q_total <- parameters["Q_total"]
Q_gut <- parameters["Q_gut"]
Q_liver <- parameters["Q_liver"]
Q_hepatic_art <- parameters["Q_hepatic_art"]
Q_kidney <- parameters["Q_kidney"]
Q_adipose <- parameters["Q_adipose"]
Q_muscle <- parameters["Q_muscle"]
Q_heart <- parameters["Q_heart"]
Q_brain <- parameters["Q_brain"]
Q_lung <- parameters["Q_lung"]
# Kp值
Kp_gut <- parameters["Kp_gut"]
Kp_liver <- parameters["Kp_liver"]
Kp_kidney <- parameters["Kp_kidney"]
Kp_adipose <- parameters["Kp_adipose"]
Kp_muscle <- parameters["Kp_muscle"]
Kp_heart <- parameters["Kp_heart"]
Kp_brain <- parameters["Kp_brain"]
Kp_lung <- parameters["Kp_lung"]
# 清除率参数
CL_renal <- parameters["CL_renal"]
# 计算组织浓度
C_liver <- A_liver / V_liver
C_gut <- A_gut / V_gut
C_kidney <- A_kidney / V_kidney
C_adipose <- A_adipose / V_adipose
C_muscle <- A_muscle / V_muscle
C_heart <- A_heart / V_heart
C_brain <- A_brain / V_brain
C_lung <- A_lung / V_lung
C_venous <- A_venous / V_venous
C_arterial <- A_arterial / V_arterial
# 门静脉浓度使用动脉浓度近似
C_portal <- C_arterial
# 计算肝脏清除率
Clint_liver <- Clint * 60 * 1e-6 * 45 * V_liver * 1e3
CLint_free <- Clint_liver * fup
CL_hepatic <- Q_liver * (CLint_free / (Q_liver + CLint_free))
# 微分方程系统
dA_gi <- -Ka * A_gi
dA_portal <- Ka * A_gi - Q_gut * C_portal
dA_liver <- Q_gut * C_portal + Q_hepatic_art * C_arterial -
(Q_gut + Q_hepatic_art) * C_liver / Kp_liver -
CL_hepatic * C_liver * fup
dA_gut <- Q_gut * (C_arterial - C_gut / Kp_gut)
dA_kidney <- Q_kidney * (C_arterial - C_kidney / Kp_kidney) -
CL_renal * C_kidney * fup
dA_adipose <- Q_adipose * (C_arterial - C_adipose / Kp_adipose)
dA_muscle <- Q_muscle * (C_arterial - C_muscle / Kp_muscle)
dA_heart <- Q_heart * (C_arterial - C_heart / Kp_heart)
dA_brain <- Q_brain * (C_arterial - C_brain / Kp_brain)
dA_lung <- Q_lung * (C_arterial - C_lung / Kp_lung)
dA_venous <- (Q_liver * C_liver / Kp_liver +
Q_kidney * C_kidney / Kp_kidney +
Q_adipose * C_adipose / Kp_adipose +
Q_muscle * C_muscle / Kp_muscle +
Q_heart * C_heart / Kp_heart +
Q_brain * C_brain / Kp_brain +
Q_lung * C_lung / Kp_lung) -
Q_total * C_venous
dA_arterial <- Q_total * (C_venous - C_arterial)
# 返回导数向量 - 关键修复:确保返回12个导数值
return(list(c(
dA_gi,
dA_portal,
dA_gut,
dA_liver,
dA_kidney,
dA_adipose,
dA_muscle,
dA_heart,
dA_brain,
dA_lung,
dA_venous,
dA_arterial
)))
}
# 6. 参数设置
parameters <- c(
Ka = Ka_human,
fup = compound_params$fup_human,
Clint = compound_params$Clint_human,
# Kp值
Kp_gut = Kp_human["gut"],
Kp_liver = Kp_human["liver"],
Kp_kidney = Kp_human["kidney"],
Kp_adipose = Kp_human["adipose"],
Kp_muscle = Kp_human["muscle"],
Kp_heart = Kp_human["heart"],
Kp_brain = Kp_human["brain"],
Kp_lung = Kp_human["lung"],
# 生理参数
V_gi = 0.25, # 胃肠道体积 (L)
V_gut = physio_human$volumes["gut"],
V_liver = physio_human$volumes["liver"],
V_kidney = physio_human$volumes["kidney"],
V_adipose = physio_human$volumes["adipose"],
V_muscle = physio_human$volumes["muscle"],
V_heart = physio_human$volumes["heart"],
V_brain = physio_human$volumes["brain"],
V_lung = physio_human$volumes["lung"],
V_venous = physio_human$other$V_venous,
V_arterial = physio_human$other$V_arterial,
# 血流参数
Q_total = physio_human$other$Q_total,
Q_gut = physio_human$flows["gut"],
Q_liver = physio_human$flows["liver"],
Q_hepatic_art = physio_human$flows["liver"] - physio_human$flows["gut"],
Q_kidney = physio_human$flows["kidney"],
Q_adipose = physio_human$flows["adipose"],
Q_muscle = physio_human$flows["muscle"],
Q_heart = physio_human$flows["heart"],
Q_brain = physio_human$flows["brain"],
Q_lung = physio_human$flows["lung"],
# 清除率参数
CL_renal = 1.5 # 肾脏清除率 (L/h)
)
# 7. 初始状态 (口服给药 100 mg)
state <- c(
A_gi = 100, # 胃肠道药物量 (mg)
A_portal = 0, # 门静脉药物量
A_gut = 0, # 肠道组织
A_liver = 0, # 肝脏
A_kidney = 0, # 肾脏
A_adipose = 0, # 脂肪组织
A_muscle = 0, # 肌肉组织
A_heart = 0, # 心脏
A_brain = 0, # 脑
A_lung = 0, # 肺
A_venous = 0, # 静脉血池
A_arterial = 0 # 动脉血池
)
# 8. 模拟时间范围 (0-24小时)
times <- seq(0, 24, by = 0.1)
# 9. 运行模型
out <- ode(y = state, times = times, func = pbpk_model, parms = parameters)
# 10. 结果处理
results <- as.data.frame(out) %>%
mutate(Concentration_plasma = A_arterial / (parameters["V_arterial"] * compound_params$Rbp_human))
# 11. 计算药代动力学参数
pk_analysis <- function(conc, time) {
# AUC计算 (梯形法)
auc <- sum(diff(time) * (head(conc, -1) + tail(conc, -1)))/2
# 终末消除速率常数
terminal_phase <- tail(data.frame(time, conc), 20)
fit <- lm(log(conc) ~ time, data = terminal_phase)
kel <- -coef(fit)[2]
t_half <- log(2)/kel
# 清除率 (基于静脉给药AUC)
dose_iv <- 100 # 假设静脉给药剂量
auc_iv <- auc * 1.2 # 假设静脉AUC比口服低
CL <- dose_iv / auc_iv
# 分布容积
Vd <- CL / kel
# 生物利用度
F_oral <- auc / auc_iv * (dose_iv / 100) # 口服剂量100mg
return(list(AUC = auc, t_half = t_half, CL = CL, Vd = Vd, F_abs = F_oral))
}
pk_params <- pk_analysis(results$Concentration_plasma, results$time)
# 12. 结果输出
cat("\n预测的人体药代动力学参数:\n")
cat(sprintf("吸收速率常数 Ka = %.3f 1/h\n", Ka_human))
cat(sprintf("清除率 CL = %.2f L/h\n", pk_params$CL))
cat(sprintf("表观分布容积 Vd = %.2f L\n", pk_params$Vd))
cat(sprintf("消除半衰期 t1/2 = %.2f h\n", pk_params$t_half))
cat(sprintf("口服生物利用度 F = %.1f%%\n", pk_params$F_abs * 100))
# 13. 绘图
ggplot(results, aes(x = time, y = Concentration_plasma)) +
geom_line(color = "blue", linewidth = 1) +
labs(title = "PBPK模型模拟 - 血浆药物浓度时间曲线",
x = "时间 (小时)",
y = "血浆浓度 (mg/L)") +
theme_minimal()
# 14. 组织分布可视化
tissue_concentrations <- results %>%
select(time, A_liver, A_kidney, A_muscle, A_adipose) %>%
pivot_longer(cols = -time, names_to = "Tissue", values_to = "Amount") %>%
mutate(Concentration = case_when(
Tissue == "A_liver" ~ Amount / parameters["V_liver"],
Tissue == "A_kidney" ~ Amount / parameters["V_kidney"],
Tissue == "A_muscle" ~ Amount / parameters["V_muscle"],
Tissue == "A_adipose" ~ Amount / parameters["V_adipose"]
))
ggplot(tissue_concentrations, aes(x = time, y = Concentration, color = Tissue)) +
geom_line(linewidth = 1) +
labs(title = "组织药物浓度时间曲线",
x = "时间 (小时)",
y = "组织浓度 (mg/L)") +
scale_color_manual(values = c("red", "green", "blue", "purple"),
labels = c("肝脏", "肾脏", "肌肉", "脂肪")) +
theme_minimal()
> pk_params <- pk_analysis(results$Concentration_plasma, results$time)
错误于lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...):
NA/NaN/Inf in 'y'
Called from: lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...)