双宿主SI模型在R语言中实现μ、α参数动态调整的技术问询

问题背景 #

研究无恢复舱室的双宿主（雄性&雌性）流行病学SI模型，需同时调整μ（取值0-1）和α（取值0-1）两个参数，实现：

• 生成不同参数组合下的动态模拟结果
• 绘制3D图或相图
• 判断种群是否消亡

已有基础代码框架如下：

# 加载所需包
library(deSolve) 
library(ggplot2)

# 定义双宿主SI模型
KModel <- function(time, state, params){
  with(as.list(c(state, params)),{
    N <- SF+IF+SM+IM
    dSF <- r*(SF+alpha*IF)-r*N*SF-BFM*(SF*IM)/N
    dIF <- (BFM*(SF*IM)/N)-r*N*IF-mu*IF
    dSM <- r*(SF+alpha*IF)-r*N*SM-BMF*(SM*IF)/N
    dIM <- (BMF*(SM*IF)/N)-r*N*IM-mu*IM
    
    return(list(c(dSF, dIF, dSM, dIM)))
  })
}

# 初始参数与状态设置
r = 0.2
BFM = 1.2
BMF = 1
mu = 0
alpha = 0

params<-c(r,BFM,BMF,mu,alpha)
initial_state<-c(SF=0.49 ,IF=0.01, SM=0.49,IM=0.01)
times<-0:60

# 数值求解
out1<-ode(y=initial_state, times=times, func=KModel, parms=params, method="ode23")
out<-as.data.frame(out1)
plot(out1)

实现方案 #

1. 批量遍历μ和α参数组合

生成参数网格，批量运行模型并存储所有结果：

# 生成μ和α的参数网格（步长可自定义，示例用0.1）
mu_vals <- seq(0, 1, by = 0.1)
alpha_vals <- seq(0, 1, by = 0.1)
param_grid <- expand.grid(mu = mu_vals, alpha = alpha_vals)

# 定义单组参数的模型运行函数
run_model <- function(mu, alpha) {
  params <- c(r = r, BFM = BFM, BMF = BMF, mu = mu, alpha = alpha)
  out <- ode(y = initial_state, times = times, func = KModel, parms = params, method = "ode23")
  as.data.frame(out) |>
    mutate(mu = mu, alpha = alpha) # 标记当前参数组合
}

# 批量运行所有参数组合（用purrr简化循环）
library(purrr)
all_results <- pmap_dfr(param_grid, run_model)

2. 动态结果可视化

用分面图展示不同μ和α组合下的种群动态变化：

ggplot(all_results, aes(x = time)) +
  geom_line(aes(y = SF, color = "易感雌性")) +
  geom_line(aes(y = IF, color = "感染雌性")) +
  geom_line(aes(y = SM, color = "易感雄性")) +
  geom_line(aes(y = IM, color = "感染雄性")) +
  facet_grid(mu ~ alpha) + # 按μ行、α列分面展示
  labs(x = "时间", y = "种群比例", color = "种群类型") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

3. 3D参数扫描图（直观判断种群消亡）

以μ和α为坐标轴，最终总种群数量为高度，绘制交互式3D图：

# 提取每个参数组合的最终时刻总种群数量
final_results <- all_results |>
  filter(time == max(times)) |>
  mutate(total_pop = SF + IF + SM + IM)

# 用plotly绘制交互式3D曲面图
library(plotly)
plot_ly(final_results, x = ~mu, y = ~alpha, z = ~total_pop, 
        type = "surface", colors = "viridis") |>
  layout(scene = list(xaxis = list(title = "μ"), 
                      yaxis = list(title = "α"), 
                      zaxis = list(title = "最终总种群比例")))

4. 相图分析（单参数组合下的相空间）

选择特定μ和α值，绘制感染雌性与感染雄性的相轨迹：

# 指定目标参数组合
target_mu <- 0.3
target_alpha <- 0.5
params_target <- c(r = r, BFM = BFM, BMF = BMF, mu = target_mu, alpha = target_alpha)
out_target <- ode(y = initial_state, times = times, func = KModel, parms = params_target, method = "ode23")
out_target_df <- as.data.frame(out_target)

# 绘制IF vs IM的相图
ggplot(out_target_df, aes(x = IF, y = IM)) +
  geom_path(color = "blue", linewidth = 1) +
  geom_point(aes(x = IF[1], y = IM[1]), color = "red", size = 3) +
  labs(x = "感染雌性比例", y = "感染雄性比例", 
       title = paste0("相图 (μ=", target_mu, ", α=", target_alpha, ")")) +
  theme_bw()

5. 种群消亡判断与可视化

设定消亡阈值，筛选并可视化导致种群消亡的参数区域：

# 设定种群消亡阈值（示例：总种群比例<0.01）
extinction_threshold <- 0.01

# 筛选消亡参数组合
extinction_params <- final_results |>
  filter(total_pop < extinction_threshold) |>
  select(mu, alpha, total_pop)

# 输出消亡参数
cat("导致种群消亡的μ和α组合：\n")
print(extinction_params)

# 热图可视化消亡区域
ggplot(final_results, aes(x = mu, y = alpha, fill = total_pop < extinction_threshold)) +
  geom_tile() +
  scale_fill_manual(values = c("FALSE" = "white", "TRUE" = "darkred"), 
                    labels = c("种群存活", "种群消亡")) +
  labs(x = "μ", y = "α", fill = "种群状态") +
  theme_bw()

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