基于deSolve的PBPK双关联微分方程组求解效率优化问询

问题背景 #

构建了包含**金属离子体内代谢（系统1）与螯合剂代谢（系统2）**的PBPK模型：

• 系统1模拟金属离子在体内的药代动力学过程
• 系统2模拟螯合剂的代谢，螯合剂可结合金属离子加速其清除，且系统2的房室螯合剂浓度会动态更新系统1对应房室的转移速率

当前实现的性能瓶颈 #

使用R语言deSolve库，通过循环单步求解系统1的微分方程组，每次迭代都根据系统2的数据更新转移速率。由于模型规模较大，这种单步循环的方式带来了巨大的调用与初始化开销，导致仿真耗时显著。

核心优化方案 #

1. 替换单步循环为分段连续求解 #

当前代码每一步都单独调用radau求解，会重复触发求解器的初始化与参数解析，是性能损耗的核心来源。优化思路是：将系统2中参数变化的时间点作为分段断点，把整个仿真时间拆分为多个参数恒定的子区间，在每个子区间内一次性求解ODE，而非单步迭代。

示例代码：

# 获取系统2的参数变化时间点
change_times <- model2$obs
# 合并并排序所有关键时间点（初始时间、原采样点、参数变化点）
all_times <- sort(unique(c(0, interval, change_times)))

# 初始化结果矩阵
dim <- length(model1$compnames)
decay_in_comps <- matrix(nrow = length(all_times), ncol = dim)
colnames(decay_in_comps) <- model1$compnames
decay_in_comps[1, ] <- input_parms  # 初始状态
current_state <- input_parms

# 分段求解ODE
for (i in 1:(length(all_times)-1)) {
  t_start <- all_times[i]
  t_end <- all_times[i+1]
  
  # 匹配当前区间的系统2参数
  row_idx <- which(model2$obs == t_start)
  if (length(row_idx) > 0) {
    model1$tcoeff["Blood","Dblood"] <- as.numeric(model2[row_idx,"Dblood"])*chelation_const
    model1$tcoeff["Liver1","DECFs"] <- as.numeric(model2[row_idx,"DECFs"])*0.07
    model1$tcoeff["ST0","DECFf"] <- as.numeric(model2[row_idx, "DECFf"])*chelation_const
  }
  
  # 求解当前子区间的ODE
  sol <- solve_eqn(model1$tcoeff, model1$compnames, current_state, times = c(t_start, t_end), ldec)
  
  # 保存结果并更新当前状态
  decay_in_comps[i+1, ] <- sol$solution[2, ]
  current_state <- sol$solution[2, ]
}

2. 简化ODE核心计算逻辑 #

当前dq函数每次都要从参数向量中重构转移矩阵，在大规模模型下会产生大量冗余计算。优化方式是直接将转移矩阵作为参数传递，避免重复的矩阵重构操作：

# 重构dq函数，直接接收转移矩阵
dq <- function(t, x, parms) {
  compsys <- parms$compsys
  qout <- x %*% compsys
  list(qout)
}

# 修改solve_eqn函数，提前计算转移矩阵并传递
solve_eqn <- function(transfers, compnames, initial, times, ldec ) {
  dim <- nrow(transfers)
  k <- transfers
  # 计算对角项
  for(i in 1:dim) {
    k[i,i] <- (-ldec - sum(transfers[i,]))
  }
  # 传递矩阵而非拆分的参数向量
  odesolution <- radau(func = dq, y = initial, times = times, parms = list(compsys = k))
  
  solution <- matrix(data = odesolution[,-1], nrow = length(times), byrow = FALSE)
  solution[solution < 1e-13] <- 0
  colnames(solution) <- compnames
  return(list("solution" = solution))
}

3. 优化系统2参数匹配逻辑 #

当前使用which(model2["obs"]==t)进行精确匹配，不仅效率低，还可能因浮点精度问题导致匹配失败。优化思路是预先排序系统2数据，使用findInterval快速定位时间对应的参数：

# 预先按时间排序系统2数据
model2_sorted <- model2[order(model2$obs), ]

# 在分段循环中替换参数匹配逻辑
t_current <- t_start
row_idx <- findInterval(t_current, model2_sorted$obs)
if (row_idx > 0) {
  # 取最近的有效参数值
  model1$tcoeff["Blood","Dblood"] <- model2_sorted$Dblood[row_idx] * chelation_const
  model1$tcoeff["Liver1","DECFs"] <- model2_sorted$DECFs[row_idx] * 0.07
  model1$tcoeff["ST0","DECFf"] <- model2_sorted$DECFf[row_idx] * chelation_const
}

4. 高阶优化选项 #

• 更换求解器：如果模型并非强刚性，可以尝试lsoda求解器（deSolve默认求解器），它会自动切换刚性/非刚性模式，部分场景下性能优于radau
• Rcpp重写核心逻辑：将ODE的核心计算（如矩阵乘法）用Rcpp重写，可大幅提升计算速度，尤其针对大规模模型
• 并行计算：利用deSolve的并行接口，或结合foreach包将分段求解任务并行化（需注意状态依赖，仅适合无交叉依赖的分段）

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