求助：基于MATLAB/R绘制基本再生数R0双参数等高线图

我正在练习复现某传染病模型的基本再生数（R₀）关于两个参数的等高线图，目标是复现Kifle等人论文中的附图。目前已写出单值计算R₀的MATLAB代码，但尝试绘制等高线的代码结果不正确，恳请帮忙用MATLAB或R完成正确的绘图。

单值计算R₀的MATLAB代码 #

%Parameters 
theta = 141302; %recruitment rate
mu = 0.001229; %natural death rate
tau = 0.45; %modification factor for A
zeta = 1/14; %influx from Q to S 
beta = 0.88; %transmission coefficient
alpha = 0.75214; %hospitalization rate
q = 0.31167; %influx from Q to I
eta_1 = 0.81692; %influx from E to Q
eta_2 = 0.02557; %influx from E to A
eta_3 = 1/7; %influx from E to I
delta_1 = 0.16673; %disease death rate for A
delta_2 = 0.00147; %disease death rate for I
delta_3 = 0.00038; %disease death rate for J
gamma_1 = 0.00827; %recovery rate for A
gamma_2 = 0.00787; %recovery rate for I
gamma_3 = 0.20186; %recovery rate for J

%Basic Reproduction Number
K_1 = eta_1 + eta_2 + eta_3 + mu;
K_2 = zeta + q + mu;
K_3 = gamma_1 + delta_1 + mu;
K_4 = alpha + gamma_2 + delta_2 + mu;
K_5 = gamma_3 + delta_3 + mu;

R_0 = beta*(tau*eta_2*K_2*K_4 + K_3*(eta_3*K_2 + eta_1*q))/(K_1*K_2*K_3*K_4);

出错的绘图尝试代码 #

[beta,eta_1] = meshgrid(0.1:0.001:1,0.1:0.001:1);

% 原代码存在两处问题：1. 语法错误（多了一个点和括号）；2. K_1依赖eta_1但未随网格更新
R_0 =  beta.*(tau.*eta_2.*K_2.*K_4 + K_3.*(eta_3.*K_2 + eta_1.*q).)./(K_1.*K_2.*K_3.*K_4);

%Drawing the plot
surf(beta,eta_1,R_0)
hold on
z2 = 0*beta + 1;
surf(beta,eta_1,z2,'MarkerFaceColor','red')

修正后的MATLAB代码（等高线图版本） #

% 固定参数（除了beta和eta_1）
theta = 141302;
mu = 0.001229;
tau = 0.45;
zeta = 1/14;
alpha = 0.75214;
q = 0.31167;
eta_2 = 0.02557;
eta_3 = 1/7;
delta_1 = 0.16673;
delta_2 = 0.00147;
delta_3 = 0.00038;
gamma_1 = 0.00827;
gamma_2 = 0.00787;
gamma_3 = 0.20186;

% 创建参数网格（调整步长为0.01提升效率，原0.001会生成1000*1000矩阵，资源消耗大）
[beta_grid, eta1_grid] = meshgrid(0.1:0.01:1, 0.1:0.01:1);

% 向量化计算K值和R0，注意K1依赖eta1_grid，需逐元素计算
K2 = zeta + q + mu;
K3 = gamma_1 + delta_1 + mu;
K4 = alpha + gamma_2 + delta_2 + mu;
K1 = eta1_grid + eta_2 + eta_3 + mu;

% 修正语法错误，向量化计算R0
R0_grid = beta_grid .* (tau.*eta_2.*K2.*K4 + K3.*(eta_3.*K2 + eta1_grid.*q)) ./ (K1.*K2.*K3.*K4);

% 绘制等高线图，重点标注R0=1的关键线
figure('Position',[100,100,800,600])
contourf(beta_grid, eta1_grid, R0_grid, 20); % 填充式等高线
hold on
contour(beta_grid, eta1_grid, R0_grid, [1 1], 'LineWidth',2, 'Color','red'); % 高亮R0=1的分界线
xlabel('Transmission coefficient (\beta)','FontSize',12)
ylabel('Influx from E to Q (\eta_1)','FontSize',12)
colorbar
title('Contour plot of R_0 vs \beta and \eta_1','FontSize',14)
grid on
hold off

可选的R代码实现 #

# 固定参数
theta <- 141302
mu <- 0.001229
tau <- 0.45
zeta <- 1/14
alpha <- 0.75214
q <- 0.31167
eta_2 <- 0.02557
eta_3 <- 1/7
delta_1 <- 0.16673
delta_2 <- 0.00147
delta_3 <- 0.00038
gamma_1 <- 0.00827
gamma_2 <- 0.00787
gamma_3 <- 0.20186

# 创建参数网格
beta_grid <- seq(0.1, 1, by=0.01)
eta1_grid <- seq(0.1, 1, by=0.01)
grid <- expand.grid(beta=beta_grid, eta1=eta1_grid)

# 计算K值和R0
K2 <- zeta + q + mu
K3 <- gamma_1 + delta_1 + mu
K4 <- alpha + gamma_2 + delta_2 + mu
K1 <- grid$eta1 + eta_2 + eta_3 + mu
R0 <- with(grid, beta * (tau*eta_2*K2*K4 + K3*(eta_3*K2 + eta1*q)) / (K1*K2*K3*K4))
grid$R0 <- R0

# 绘制等高线图
library(ggplot2)
ggplot(grid, aes(x=beta, y=eta1, z=R0)) +
  geom_contour_filled(bins=20) +
  geom_contour(aes(z=R0), breaks=1, color="red", linewidth=1) +
  labs(x="Transmission coefficient (β)", y="Influx from E to Q (η₁)", 
       title="Contour plot of R₀ vs β and η₁",
       fill="R₀ Value") +
  theme_minimal(base_size=12) +
  theme(plot.title = element_text(hjust=0.5))

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