Matlab转R时ode45求解结果差异过大的技术咨询

本人对编程及Stack Overflow并不熟悉，尝试复现某文献中的细胞动力学模型工作，因无MATLAB权限改用R语言实现。

编写的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: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 = 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 

实验数据 #

求助内容 #

尝试增大maxsteps仍无法消除警告；且使用文献给出的参数进行模拟时，结果反而更差。怀疑是代码存在错误或ode45求解器的问题，恳请协助排查原因。

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