Matlab → R，迷失在翻译中：为什么 ode45 会产生完全不同的结果？

Question

我既不熟悉编码也不熟悉 stackoverflow，我正在尝试复制以下工作：https://doi.org/10.1007/978-1-0716-0191-4_12 使用本书章节中提供的数据。

我无法访问 MATLAB，我尝试在试用版中运行代码，但错误太多。 因此，我决定和R一起去。

这是我的代码（问题在帖子末尾描述）：

# Loads the experimental data
library(readxl)
Data <- readxl::read_excel("exp_data.xlsx", sheet = "F1") 
Data <- data.frame(Data)

# Defines starting values from experimental data 
c0 <- c(Xv = 3.024E8, Xd = 2.01E7, Glc = 37.42, Gln = 7.59)

# Definition of model parameters 
Parameters <- c(mumax = 0.05, mudmax = 0.03, mudmin = 0.003, qGlcmax = 3E-10, qGlnmax = 8.5E-11, KsGlc = 0.03, KsGln = 0.03, KGlc = 0.19, KGln = 1) 

# Defines the model
Model <- function(t, c, Parameters) { 

# Renames

mumax <- Parameters[1] 
mudmax <- Parameters[2] 
mudmin <- Parameters[3] 
qGlcmax <- Parameters[4] 
qGlnmax <- Parameters[5] 
KsGlc <- Parameters[6] 
KsGln <- Parameters[7] 
KGlc <- Parameters[8] 
KGln <- Parameters[9] 

Xv <- c[1] # viable cell density 
Xd <- c[2] # dead cell density 
Glc <- c[3] # glucose concentration 
Gln <- c[4] # glutamine concentration 

# Representation of the kinetic relationships 

mu <- mumax*(Glc /(Glc+KsGlc))*(Gln /(Gln+KsGln)) 
mud <- mudmin + mudmax*(KsGlc/(KsGlc+Glc)) 
qglc <- qGlcmax*(Glc/(Glc+KsGlc))*(mu/(mu+mumax)+0.5) 
qgln <- qGlnmax*(Gln/(Gln+KsGln)) 

# Process model, here for batch 

dcdt <- numeric(4) 
dcdt[1] <- (mu-mud)*Xv #Xv 
dcdt[2] <- mud*Xv #Xd 
dcdt[3] <- -qglc*Xv #cglc 
dcdt[4] <- -qgln*Xv #cgln 

# Current concentration changes as output 
return(list(dcdt)) }

# Time steps to be simulated 
tspan <- 0:1:200 # [h] 

# Call of model function, solved for tspan 
library(deSolve) 
prior <- ode(y = c0, times = tspan, func = Model, parms = Parameters, method = "ode45") # prior gives a consistent result whith a cell growth even if the values are too low compared to experimental data

# Compares prior to experimental data, define objective function

tspan <- Data$time #experimental tspan up to 180h only 
Weighting <- c(100, 1, 10, 100)
Magnitude <- c(1E-9, 1E-8, 1, 1)   

objective <- function(Parameters, Weighting) { 
# Call of ode system 
c <- ode(func = Model, times = tspan, y = data_c0, parms = Parameters, method = "ode45") 
# Calculate sum of squares
Sum_of_squares <- Weighting[1]*sum((abs(c[,1]- Data[,2])*Magnitude[1])^2) + Weighting[2]*sum((abs(c[,2]-Data[,3])*Magnitude[2])^2) + Weighting[3]*sum((abs(c[,3]-Data[,4])*Magnitude[3])^2) + Weighting[4]*sum((abs(c[,4]-Data[,5])*Magnitude[4])^2) 
return(Sum_of_squares) } 

# Estimation of model parameters with Nelder-Mead algorithm 
Estimated_Parameters <- optim(par = Parameters, fn = objective, Weighting = Weighting, method = "Nelder-Mead")$par 
print(Estimated_Parameters)

运行 optim 函数后我不得不提到以下警告：

There were 50 or more warnings (use warnings() to see the first 50)

warnings()

1: In rk(y, times, func, parms, method = "ode45", ...) :
  Number of time steps 59966 exceeded maxsteps at t = 0.000199293
2: In rk(y, times, func, parms, method = "ode45", ...) :
  Number of time steps 61726 exceeded maxsteps at t = 0.00318004
3: In rk(y, times, func, parms, method = "ode45", ...) :
  Number of time steps 59978 exceeded maxsteps at t = 0.000471839
4: In rk(y, times, func, parms, method = "ode45", ...) :
  Number of time steps 59971 exceeded maxsteps at t = 0.00039857
5: In rk(y, times, func, parms, method = "ode45", ...) :
  Number of time steps 61904 exceeded maxsteps at t = 0.0050244
...

如题所述，我得到的结果与文档中的结果完全不同。 虽然我应该得到以下结果：

3.7900e-02 4.2100e-02 2.4000e-03 6.2000e-11 4.5000e-12 4.3800e-02 3.2800e-02 4.3300e-02 1.4787e+00

我得到了这个：

7.297162e-02  1.390685e-02 -1.812314e-02  3.653943e-08 -5.175317e-02  7.573135e-02  6.733384e-02  2.190989e-01  1.033314e+00 

我试图增加最大步数以消除警告但没有成功。

出于好奇，我也尝试用我应该得到的 Estimated_Parameters 来运行模拟，结果……惊喜！仿真结果比原来的差。这就是为什么我现在怀疑解算器，但我认为代码更有可能出错。告诉我是否需要更多详细信息，谢谢！