基于R语言deSolve包事件功能的质量平衡校验方法问询

对于带定时注射事件的单室药物模型，质量平衡的核心逻辑是：总注入药物量 = 任意时刻血液中剩余药量 + 从初始时刻到该时刻的药物消除总量，理想状态下三者的差值（总注入 - 剩余 - 消除）应趋近于0（数值误差范围内）。以下是具体的校验步骤和代码实现：

校验步骤 #

• 计算总注入药量：统计模拟时间范围内所有注射事件的药物总量
• 计算累积消除量：通过数值积分对消除速率（b*blood）进行积分，得到从初始到每个时刻的总消除药量
• 计算平衡差值：用总注入量减去当前剩余药量和累积消除量，验证差值是否接近0

完整校验代码 #

library(deSolve)
b <- 0.6 
yini <- c(blood = 0)

pharmaco2 <- function(t, blood, p) { 
  dblood <- -b * blood 
  list(dblood) 
}

injectevents <- data.frame(var = "blood", time = 0:20, value = 40, method = "add")

times <- seq(from = 0, to = 10, by = 1/24) 

# 运行原模型
out2 <- ode(func = pharmaco2, 
            times = times, 
            y = yini, 
            parms = NULL, 
            method = "impAdams", 
            events = list(data = injectevents))

# 1. 计算模拟周期内的总注入药量
max_sim_time <- max(times)
total_injected <- sum(injectevents$value[injectevents$time <= max_sim_time])

# 2. 用梯形法计算每个时刻的累积消除量
t_vals <- out2[, "time"]
blood_vals <- out2[, "blood"]
elimination_rate <- b * blood_vals  # 实时消除速率

# 梯形法数值积分：每个时间间隔的积分值为(速率1+速率2)/2 * 时间差
cumulative_elimination <- c(0, cumsum(
  (elimination_rate[-1] + elimination_rate[-length(elimination_rate)])/2 * diff(t_vals)
))

# 3. 计算质量平衡差值
mass_balance_diff <- total_injected - blood_vals - cumulative_elimination

# 合并结果到输出数据框
out_with_balance <- cbind(out2, cumulative_elimination, mass_balance_diff)

# 查看最后几个时刻的平衡情况
print(tail(out_with_balance))

# 绘制质量平衡曲线，直观校验
plot(t_vals, mass_balance_diff, type = "l", lwd = 2, 
     xlab = "天数", ylab = "质量平衡差值", main = "单室模型质量平衡校验")
abline(h = 0, col = "red", lty = 2)

结果说明 #

• 输出的mass_balance_diff列即为质量平衡差值，正常情况下应在0附近小幅波动，波动幅度由数值求解器精度和时间步长决定（步长越小，误差越小）
• 当模拟时间远长于药物半衰期（本例半衰期约为ln(2)/0.6≈1.155天），血液中剩余药量趋近于0，累积消除量会接近总注入量，平衡差值趋近于0

内容的提问来源于stack exchange，提问作者eod