library(deSolve)
library(tidyverse)
library(pracma)
1. 定义口服PBPK模型函数（包含多种制剂参数）
pk_model_oral <- function(time, state, params) {
with(as.list(c(state, params)), {
# 计算胃肠道不同区段的pH值（随时间变化）
if (time < gastric_emptying_time) {
# 胃部环境
current_pH <- gastric_pH
} else if (time < (gastric_emptying_time + small_intestine_transit_time)) {
# 小肠环境
current_pH <- small_intestine_pH
} else {
# 结肠环境
current_pH <- colon_pH
}
# 计算pH依赖的溶解度 F_solubility <- calculate_pH_dependent_solubility( solubility_pH1 = solubility_pH1, solubility_pH2 = solubility_pH2, solubility_pH3 = solubility_pH3, solubility_pH4 = solubility_pH4, solubility_pH5 = solubility_pH5, solubility_pH6 = solubility_pH6, current_pH = current_pH, pKa = pKa, acid_base = acid_base ) # 计算pH依赖的logD current_logD <- calculate_pH_dependent_logD( logD_neutral = logD_neutral, pKa = pKa, acid_base = acid_base, current_pH = current_pH ) # 计算pH依赖的渗透性 F_permeability <- calculate_pH_dependent_permeability( Papp_ref = Papp_ref, Papp_compound = Papp_compound, current_pH = current_pH, pKa = pKa, acid_base = acid_base ) # 计算吸收速率 absorbed <- calculate_absorption( formulation_type = formulation_type, time = time, A_gi = A_gi, ka = ka, F_abs_intrinsic = F_abs_intrinsic, F_solubility = F_solubility, F_permeability = F_permeability, F_particle_size = F_particle_size, F_formulation = F_formulation, release_start = release_start, release_end = release_end, release_rate = release_rate, delay_time = delay_time, ka_factor = ka_factor ) dA_gi <- -absorbed dA_c <- absorbed - (CL / Vd) * A_c C_plasma <- (A_c / Vd) * 1000 # μg/mL -> ng/mL (血浆浓度) # 计算全血浓度 (考虑血/浆比) C_blood <- C_plasma / blood_plasma_ratio # 计算游离药物浓度 (考虑血浆蛋白结合) C_free <- C_plasma * fu_plasma return(list(c(dA_gi, dA_c), C_plasma = C_plasma, C_blood = C_blood, C_free = C_free, current_pH = current_pH, F_solubility = F_solubility, F_permeability = F_permeability, current_logD = current_logD))
})
}
2. pH依赖性溶解度计算函数
calculate_pH_dependent_solubility <- function(solubility_pH1, solubility_pH2, solubility_pH3,
solubility_pH4, solubility_pH5, solubility_pH6,
solubility_pH7, current_pH, pKa, acid_base) {
创建溶解度-pH数据框
solubility_data <- data.frame(
pH = c(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0),
solubility = c(solubility_pH1, solubility_pH2, solubility_pH3,
solubility_pH4, solubility_pH5, solubility_pH6)
)
插值计算当前pH下的溶解度
F_solubility <- approx(solubility_data$pH, solubility_data$solubility,
xout = current_pH, rule = 2)$y
对于弱酸/弱碱，应用Henderson-Hasselbalch方程校正
if (acid_base == “acid” && !is.na(pKa)) {
# 弱酸: S = S0 * (1 + 10^(pKa - pH))
F_solubility <- solubility_pH7 * (1 + 10^(pKa - current_pH))
} else if (acid_base == “base” && !is.na(pKa)) {
# 弱碱: S = S0 * (1 + 10^(pH - pKa))
F_solubility <- solubility_pH7 * (1 + 10^(current_pH - pKa))
}
return(F_solubility)
}
3. pH依赖性logD计算函数
calculate_pH_dependent_logD <- function(logD_neutral, pKa, acid_base, current_pH) {
if (is.na(pKa)) return(logD_neutral)
if (acid_base == “acid”) {
# 弱酸: logD = logD_neutral - log(1 + 10^(pH - pKa))
logD <- logD_neutral - log10(1 + 10^(current_pH - pKa))
} else if (acid_base == “base”) {
# 弱碱: logD = logD_neutral - log(1 + 10^(pKa - pH))
logD <- logD_neutral - log10(1 + 10^(pKa - current_pH))
} else {
logD <- logD_neutral
}
return(logD)
}
4. pH依赖性渗透性计算函数
calculate_pH_dependent_permeability <- function(Papp_ref, Papp_compound, current_pH, pKa, acid_base) {
参考渗透性 (美托洛尔作为阳性对照)
ref_permeability <- 2e-4 # cm/s (高渗透性参考值)
计算基本渗透性因子
F_permeability <- Papp_compound / Papp_ref
考虑pH对分子形态的影响
if (!is.na(pKa) && !is.na(acid_base)) {
if (acid_base == “acid”) {
# 弱酸在低pH时主要以非离子形式存在，渗透性高
if (current_pH < pKa) {
F_permeability <- F_permeability * 1.5
} else {
F_permeability <- F_permeability * 0.7
}
} else if (acid_base == “base”) {
# 弱碱在高pH时主要以非离子形式存在，渗透性高
if (current_pH > pKa) {
F_permeability <- F_permeability * 1.5
} else {
F_permeability <- F_permeability * 0.7
}
}
}
限制在合理范围
F_permeability <- min(max(F_permeability, 0.1), 2.0)
return(F_permeability)
}
5. 吸收计算函数
calculate_absorption <- function(formulation_type, time, A_gi, ka,
F_abs_intrinsic, F_solubility, F_permeability,
F_particle_size, F_formulation,
release_start, release_end, release_rate,
delay_time, ka_factor) {
根据剂型计算吸收速率
if (formulation_type == “IR”) {
# 速释制剂：一级吸收
absorbed <- ka * ka_factor * A_gi * F_abs_intrinsic * F_solubility * F_permeability
} else if (formulation_type %in% c(“ER”, “SR”)) {
# 缓释制剂：零级释放（在特定时间段内恒定释放）
if (time >= release_start & time <= release_end) {
absorbed <- release_rate * F_abs_intrinsic * F_solubility * F_permeability
} else {
absorbed <- 0
}
} else if (formulation_type == “DR”) {
# 延迟释放：延迟一段时间后开始一级吸收
if (time >= delay_time) {
absorbed <- ka * ka_factor * A_gi * F_abs_intrinsic * F_solubility * F_permeability
} else {
absorbed <- 0
}
} else if (formulation_type == “Nanoparticle”) {
# 纳米制剂：增强的一级吸收
absorbed <- ka * ka_factor * 1.5 * A_gi * F_abs_intrinsic * F_solubility * F_permeability
}
应用粒径和剂型影响因子
absorbed <- absorbed * F_particle_size * F_formulation
return(absorbed)
}
6. 动物口服PK数据准备（增加种属特异性参数）
animal_data <- data.frame(
species = c(“mouse”, “rat”, “dog”, “human”),
weight_kg = c(0.02, 0.25, 9, 70),
CL_animal = c(4.64, 7.61, 1.76, 0.8), # 清除率 (L/h/kg)
Vd_animal = c(8.53, 22.8, 4.67, 6.0), # 分布容积 (L/kg)
F_animal = c(0.16, 0.03, 0.29, NA), # 口服生物利用度
blood_plasma_ratio = c(0.75, 0.85, 0.95, 1.0), # 全血/血浆药物浓度比
fu_plasma = c(0.25, 0.35, 0.15, 0.05) # 血浆蛋白游离药物分数
)
7. 外推人体参数（包含新增参数）
extrapolate_human_pk <- function(animal_data, logD_neutral, pKa, acid_base,
solubility_pH1, solubility_pH2, solubility_pH3,
solubility_pH4, solubility_pH5, solubility_pH6, solubility_pH7,
Papp_ref, Papp_compound, dose_mg,
particle_size_um = 10, formulation_type = “IR”,
human_weight = 70) {
基础PK参数计算
b_CL <- 0.75 # 清除率异速缩放指数
b_Vd <- 1.0 # 分布容积异速缩放指数
使用异速缩放方程外推人体参数
human_CL <- mean(animal_data$CL_animal[1:3]) * (human_weight/mean(animal_data$weight_kg[1:3]))^b_CL
human_Vd <- mean(animal_data$Vd_animal[1:3]) * (human_weight/mean(animal_data$weight_kg[1:3]))^b_CL
计算吸收速率常数 (ka) - 使用logD替代logP
human_ka <- 10^(0.54 * logD_neutral - 0.95) # 单位 h⁻¹
固有生物利用度（不考虑制剂因素）
F_abs_intrinsic <- mean(animal_data$F_animal[1:3], na.rm = TRUE)
=== 制剂参数影响计算 ===
1. 粒径影响因子 (F_particle_size)
if (particle_size_um <= 10) {
F_particle_size <- 1
} else if (particle_size_um <= 50) {
F_particle_size <- 0.8
} else {
F_particle_size <- 0.6
}
2. 剂型影响因子 (F_formulation)
formulation_factors <- c(
“IR” = 1.0, # 速释制剂
“ER” = 0.8, # 缓释制剂
“SR” = 0.85, # 缓释制剂
“DR” = 0.9, # 延迟释放
“Nanoparticle” = 1.2 # 纳米制剂
)
F_formulation <- formulation_factors[[formulation_type]]
计算最终生物利用度 (F) - 溶解度、渗透性在模型中动态计算
human_F <- F_abs_intrinsic * F_particle_size * F_formulation
限制在0-1范围内
human_F <- max(0, min(1, human_F))
计算半衰期 (h)
human_t_half <- log(2) * human_Vd / human_CL
种属特异性参数
human_blood_plasma_ratio <- animal_data$blood_plasma_ratio[animal_data$species == “human”]
human_fu_plasma <- animal_data$fu_plasma[animal_data$species == “human”]
return(list(
CL = human_CL,
Vd = human_Vd,
ka = human_ka,
F = human_F,
t_half = human_t_half,
logD_neutral = logD_neutral,
pKa = pKa,
acid_base = acid_base,
solubility_pH1 = solubility_pH1,
solubility_pH2 = solubility_pH2,
solubility_pH3 = solubility_pH3,
solubility_pH4 = solubility_pH4,
solubility_pH5 = solubility_pH5,
solubility_pH6 = solubility_pH6,
solubility_pH7 = solubility_pH7,
Papp_ref = Papp_ref,
Papp_compound = Papp_compound,
F_particle_size = F_particle_size,
F_formulation = F_formulation,
F_abs_intrinsic = F_abs_intrinsic,
blood_plasma_ratio = human_blood_plasma_ratio,
fu_plasma = human_fu_plasma
))
}
8. 设置化合物特性和制剂参数
logD_neutral <- 2.5 # 中性条件下的油水分配系数
pKa_value <- 4.5 # 酸碱解离平衡常数
acid_base <- “acid” # 酸性或碱性化合物(“acid"或"base”)
pH依赖性溶解度 (mg/mL)
solubility_pH1 <- 1.4 # pH 1.0
solubility_pH2 <- 2.5 # pH 2.0
solubility_pH3 <- 4.0 # pH 3.0
solubility_pH4 <- 5.0 # pH 4.0
solubility_pH5 <- 6.5 # pH 5.0
solubility_pH6 <- 7.4 # pH 6.0
渗透性参数
Papp_ref <- 20.5e-6 # 阳性对照物(美托洛尔)的Papp值 (cm/s)
Papp_compound <- 15.2e-6 # 目标化合物的Papp值 (cm/s)
制剂参数
dose <- 200 # 口服剂量 (mg)
human_weight <- 68 # 人体体重 (kg)
particle_size_um <- 15 # 粒径 (μm)
formulation_type <- “IR” # 剂型 (“IR”, “ER”, “SR”, “DR”, “Nanoparticle”)
胃肠道生理参数
gastric_pH <- 1.5 # 胃部pH
small_intestine_pH <- 6.5 # 小肠pH
colon_pH <- 7.0 # 结肠pH
gastric_emptying_time <- 0.5 # 胃排空时间 (小时)
small_intestine_transit_time <- 4 # 小肠转运时间 (小时)
9. 外推人体参数 (包含新增参数)
human_pk <- extrapolate_human_pk(
animal_data,
logD_neutral = logD_neutral,
pKa = pKa_value,
acid_base = acid_base,
solubility_pH1 = 0.08557,
solubility_pH2 = 0.06794,
solubility_pH3 = 0.07188,
solubility_pH4 = 0.07629,
solubility_pH5 = 0.08165,
solubility_pH6 = 0.07515,
Papp_ref = Papp_ref,
Papp_compound = Papp_compound,
dose_mg = dose,
particle_size_um = particle_size_um,
formulation_type = formulation_type,
human_weight = human_weight
)
10. 设置剂型特定参数
set_formulation_params <- function(formulation_type) {
params <- list()
if (formulation_type == “IR”) {
params$release_start <- 0
params$release_end <- 0
params$release_rate <- 0
params$delay_time <- 0
params$ka_factor <- 1.0
} else if (formulation_type == “ER”) {
params$release_start <- 0.5
params$release_end <- 12
params$release_rate <- dose / (params$release_end - params$release_start)
params$delay_time <- 0
params$ka_factor <- 1.0
} else if (formulation_type == “SR”) {
params$release_start <- 0.5
params$release_end <- 24
params$release_rate <- dose / (params$release_end - params$release_start)
params$delay_time <- 0
params$ka_factor <- 1.0
} else if (formulation_type == “DR”) {
params$release_start <- 0
params$release_end <- 0
params$release_rate <- 0
params$delay_time <- 1.5
params$ka_factor <- 1.0
} else if (formulation_type == “Nanoparticle”) {
params$release_start <- 0
params$release_end <- 0
params$release_rate <- 0
params$delay_time <- 0
params$ka_factor <- 1.5
}
return(params)
}
form_params <- set_formulation_params(formulation_type)
11. 初始状态设置
state <- c(
A_gi = dose, # 胃肠道初始药量 (mg)
A_c = 0 # 中央室初始药量 (mg)
)
12. 完整参数集
params <- c(
CL = human_pk$CL,
Vd = human_pk$Vd,
ka = human_pk$ka,
F_abs_intrinsic = human_pk$F_abs_intrinsic,
formulation_type = formulation_type,
release_start = form_params$release_start,
release_end = form_params$release_end,
release_rate = form_params$release_rate,
delay_time = form_params$delay_time,
ka_factor = form_params$ka_factor,
F_particle_size = human_pk$F_particle_size,
F_formulation = human_pk$F_formulation,
solubility_pH1 = human_pk$solubility_pH1,
solubility_pH2 = human_pk$solubility_pH2,
solubility_pH3 = human_pk$solubility_pH3,
solubility_pH4 = human_pk$solubility_pH4,
solubility_pH5 = human_pk$solubility_pH5,
solubility_pH6 = human_pk$solubility_pH6,
pKa = human_pk$pKa,
acid_base = human_pk$acid_base,
logD_neutral = human_pk$logD_neutral,
Papp_ref = human_pk$Papp_ref,
Papp_compound = human_pk$Papp_compound,
gastric_pH = gastric_pH,
small_intestine_pH = small_intestine_pH,
colon_pH = colon_pH,
gastric_emptying_time = gastric_emptying_time,
small_intestine_transit_time = small_intestine_transit_time,
blood_plasma_ratio = human_pk$blood_plasma_ratio,
fu_plasma = human_pk$fu_plasma
)
13. 模型求解
times <- seq(0, 24, by = 0.1)
solution <- ode(
y = state,
times = times,
func = pk_model_oral,
parms = params
)
14. 结果处理
results <- as.data.frame(solution) %>%
mutate(
Concentration_plasma = ifelse(is.na(C_plasma), 0, C_plasma),
Concentration_blood = ifelse(is.na(C_blood), 0, C_blood),
Concentration_free = ifelse(is.na(C_free), 0, C_free),
pH = ifelse(is.na(current_pH), gastric_pH, current_pH),
Solubility_factor = ifelse(is.na(F_solubility), 1, F_solubility),
Permeability_factor = ifelse(is.na(F_permeability), 1, F_permeability),
logD = ifelse(is.na(current_logD), logD_neutral, current_logD)
)
15. 可视化（多面板显示）
p1 <- ggplot(results, aes(x = time, y = Concentration_plasma)) +
geom_line(color = “blue”, linewidth = 1.2) +
labs(title = “血浆药物浓度”,
x = “时间 (小时)”,
y = “浓度 (ng/mL)”) +
theme_minimal()
p2 <- ggplot(results, aes(x = time, y = pH)) +
geom_line(color = “red”, linewidth = 1.2) +
labs(title = “胃肠道pH变化”,
x = “时间 (小时)”,
y = “pH”) +
theme_minimal()
p3 <- ggplot(results, aes(x = time)) +
geom_line(aes(y = Solubility_factor, color = “溶解度因子”), linewidth = 1.2) +
geom_line(aes(y = Permeability_factor, color = “渗透性因子”), linewidth = 1.2) +
geom_line(aes(y = logD/max(logD), color = “logD (归一化)”), linewidth = 1.2) +
scale_color_manual(values = c(“溶解度因子” = “blue”,
“渗透性因子” = “red”,
“logD (归一化)” = “green”)) +
labs(title = “吸收相关参数变化”,
x = “时间 (小时)”,
y = “因子值”) +
theme_minimal() +
theme(legend.position = “top”)
p4 <- ggplot(results, aes(x = time)) +
geom_line(aes(y = Concentration_plasma, color = “血浆总浓度”), linewidth = 1.2) +
geom_line(aes(y = Concentration_free, color = “游离药物浓度”), linewidth = 1.2) +
geom_line(aes(y = Concentration_blood, color = “全血浓度”), linewidth = 1.2) +
scale_color_manual(values = c(“血浆总浓度” = “blue”,
“游离药物浓度” = “red”,
“全血浓度” = “purple”)) +
labs(title = “不同浓度指标比较”,
x = “时间 (小时)”,
y = “浓度 (ng/mL)”) +
theme_minimal() +
theme(legend.position = “top”)
使用patchwork包组合图形
library(patchwork)
combined_plot <- (p1 + p2) / (p3 + p4) +
plot_annotation(
title = “增强型口服PBPK模型模拟结果”,
subtitle = sprintf(“化合物: %s (pKa=%.1f) | 剂型: %s | 剂量: %d mg”,
acid_base, pKa_value, formulation_type, dose),
theme = theme(plot.title = element_text(size = 16, face = “bold”))
)
print(combined_plot)
16. 输出关键PK参数
Cmax_plasma <- max(results$Concentration_plasma, na.rm = TRUE)
Tmax <- results$time[which.max(results$Concentration_plasma)]
AUC <- trapz(results$time, results$Concentration_plasma)
cat(“=== 增强型PBPK模型参数汇总 =\n")
cat(sprintf(“化合物类型: %s (pKa=%.1f)\n”, acid_base, pKa_value))
cat(sprintf(“中性logD: %.2f\n”, logD_neutral))
cat(sprintf(“吸收速率常数 (ka): %.2f h⁻¹\n”, human_pk$ka))
cat(sprintf(“清除率 (CL): %.2f L/h\n”, human_pk$CL))
cat(sprintf(“分布容积 (Vd): %.2f L\n”, human_pk$Vd))
cat(sprintf(“全血/血浆浓度比: %.2f\n”, human_pk$blood_plasma_ratio))
cat(sprintf(“血浆蛋白游离分数 (fu): %.3f\n”, human_pk$fu_plasma))
cat(sprintf(“Caco-2 Papp: %.2e cm/s (参考值: %.2e cm/s)\n”,
Papp_compound, Papp_ref))
cat(sprintf(“半衰期 (t_{1/2}): %.2f h\n”, human_pk$t_half))
cat(sprintf(“峰浓度 (Cmax): %.2f ng/mL\n”, Cmax_plasma))
cat(sprintf(“达峰时间 (Tmax): %.1f h\n”, Tmax))
cat(sprintf(“药时曲线下面积 (AUC₀-₂₄): %.1f ng·h/mL\n”, AUC))
cat("\n= pH依赖性溶解度 ===\n”)
cat(sprintf(“pH 1.0: %.2f mg/mL | pH 4.0: %.2f mg/mL\n”,
solubility_pH1, solubility_pH4))
cat(sprintf(“pH 2.0: %.2f mg/mL | pH 5.0: %.2f mg/mL\n”,
solubility_pH2, solubility_pH5))
cat(sprintf(“pH 3.0: %.2f mg/mL | pH 6.0: %.2f mg/mL\n”,
solubility_pH3, solubility_pH6))
cat(sprintf(“pH 7.0: %.2f mg/mL\n”, solubility_pH7))
错误于extrapolate_human_pk(animal_data, logD_neutral = logD_neutral, :
缺少参数"solubility_pH7",也缺失默认值，只有6个PH值的溶解度，删除"solubility_pH7"。
扫一扫
3528
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
