Matlab → R, lost in translation: why ode45 yields totally different results?

Question

I am neither very familliar with coding nor stackoverflow, I am trying to replicate this following work: https://doi.org/10.1007/978-1-0716-0191-4_12 using the data provided in the book chapter.

I do not have access to MATLAB, I tried to run the code in a trial version but there was too many errors. Thus, I decided to go with R.

Here is my code (issue is described at the end of the post):

# 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)

I have to mention the following warning after runing the optim function:

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
...

As mentioned in the title, the results I get are totally different from those in the document. While I should get the following result:

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

I obtained this:

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 

I tried to increase the maxsteps to get rid of the warnings without success.

Out of curiosity, I also tried to run the simulation with the Estimated_Parameters I should get and... surprise! The simulation result is worse than the original. That's why I suspect the solver for now, but I assume a mistake with the code is more likely. Tell me if more details are needed, thanks!

Edit: experimental data table ↓

time	Xv	Xd	glc	gln
0	302000000	20000000	37.42	7.59
12	396000000	25700000	36.81	7.09
24	621000000	43600000	34.00	6.40
36	1020000000	59300000	35.02	5.50
48	1470000000	100000000	31.73	4.64
60	2120000000	160000000	25.58	3.65
72	3520000000	222000000	23.05	2.49
84	5610000000	329000000	19.70	1.41
96	7510000000	503000000	16.81	0.33
108	10000000000	613000000	16.23	0.11
120	11300000000	766000000	12.71	0.09
132	12700000000	866000000	8.63	0.11
144	12900000000	1000000000	4.36	0.15
156	8260000000	6080000000	2.14	0.17
168	6350000000	8860000000	1.39	0.16
180	4640000000	9710000000	0.49	0.16