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 |
