this is my first tryout performing a model optimization. I simulate the model using DeSolve and try to optimize using an objective function coupled to the optim-function in R. I generally apply either Nelder-Mead or L-BFGS-B.
Here is the model function:
library(deSolve)
model <- function(t, y, parms) {
with(as.list(c(y, parms)), {
y1=400*sin((2*pi)/(10^Tlog)*t)+600
dY2 = ((kA1*(Atot-y2)*y1^nA)/(KmA^nA+y1^nA)-kA2*y2)
dY3 = ((kD1*(Dtot-y3)*y1^nD)/(KmD^nD+y1^nD)-kD2*y3)
dY4 = ((kP1*(Ptot-y4)*y2^nP1)/(KmP1^nP1+y2^nP1)-(kP2*(y4)*y3^nP2)/(KmP2^nP2+y3^nP2))
dY5 = y2
dY6 = y3
dY7 = y4
list(c(dY2, dY3, dY4, dY5, dY6, dY7))
})
}
Here the optimization function - I simulate three times per optimization iteration. Note that I optimize only a subset of all parameters, one specific parameter (Tlog) is changed target-oriented for each of the three simulations:
ls <- function(par,extra) {
if(min(par)<0.01 || max(par) >10000) {
return(NA)
}
with(as.list(c(par, extra)), {
#1
parms <- c(vol=vol, kA1=kA1, nA=nA, KmA=KmA, Atot=Atot, kA2=kA2,
kD1=kD1, nD=nD, KmD=KmD, Dtot=Dtot, kD2=kD2,
kP1=kP1, nP1=nP1, KmP1=KmP1,Ptot=Ptot,kP2=kP2,nP2=nP2,KmP2=KmP2,
Tlog=-1)
yini <- c(y2 = 0, y3 = 0, y4 = 0, y5 = 0, y6 = 0, y7 = 0)
times <- seq(from = 0, to = 100, by = 0.01)
out <- as.data.frame(ode (times = times, y = yini, func = model, parms = parms, ynames = FALSE))
out[,5:7]=out[,5:7]/out[,1]
colnames(out)=c("t[s]","A[nM]","D[nM]","P[nM]","A_act[nM]","D_act[nM]","P_act[nM]")
O_1=(out$`P_act[nM]`[nrow(out)-1])
#2
parms <- c(vol=vol, kA1=kA1, nA=nA, KmA=KmA, Atot=Atot, kA2=kA2,
kD1=kD1, nD=nD, KmD=KmD, Dtot=Dtot, kD2=kD2,
kP1=kP1, nP1=nP1, KmP1=KmP1,Ptot=Ptot,kP2=kP2,nP2=nP2,KmP2=KmP2,
Tlog=1)
# y2 = A; y3 = D; y4 = P; y5 = A_Act; y6 = D_Act; y7 = P_Act
yini <- c(y2 = 0, y3 = 0, y4 = 0, y5 = 0, y6 = 0, y7 = 0)
times <- seq(from = 0, to = 1000, by = 0.01)
out <- as.data.frame(ode (times = times, y = yini, func = model, parms = parms, ynames = FALSE))
out[,5:7]=out[,5:7]/out[,1]
colnames(out)=c("t[s]","A[nM]","D[nM]","P[nM]","A_act[nM]","D_act[nM]","P_act[nM]")
O_2=(out$`P_act[nM]`[nrow(out)-1])
#3
parms <- c(vol=vol, kA1=kA1, nA=nA, KmA=KmA, Atot=Atot, kA2=kA2,
kD1=kD1, nD=nD, KmD=KmD, Dtot=Dtot, kD2=kD2,
kP1=kP1, nP1=nP1, KmP1=KmP1,Ptot=Ptot,kP2=kP2,nP2=nP2,KmP2=KmP2,
Tlog=3)
# y2 = A; y3 = D; y4 = P; y5 = A_Act; y6 = D_Act; y7 = P_Act
yini <- c(y2 = 0, y3 = 0, y4 = 0, y5 = 0, y6 = 0, y7 = 0)
times <- seq(from = 0, to = 10000, by = 0.01)
out <- as.data.frame(ode (times = times, y = yini, func = model, parms = parms, ynames = FALSE))
out[,5:7]=out[,5:7]/out[,1]
colnames(out)=c("t[s]","A[nM]","D[nM]","P[nM]","A_act[nM]","D_act[nM]","P_act[nM]")
O_3=(out$`P_act[nM]`[nrow(out)-1])
if(O_2<(max(c(O_1,O_3)))) {
return(NA)
}
Objective=1/(abs(O_2-O_1))+1/(abs(O_2-O_3))
print(c(O_1,O_2,O_3))
print(Objective)
print("####")
return(Objective)
})
}
And finally here is the call of the optimization function with the initial parameter set:
P_init=matrix(ncol=1,nrow=10)
for(i in 1:length(P_init)) {
P_init[i,1]=runif(1,0.001,10000)
}
Res=optim(par=c(kA1=P_init[1,1], KmA=P_init[2,1], kA2=P_init[3,1],
kD1=P_init[4,1], KmD=P_init[5,1],kD2=P_init[6,1],
kP1=P_init[7,1], KmP1=P_init[8,1],kP2=P_init[9,1],KmP2=P_init[10,1]),
extra=c(vol=1,nA=4,nD=4,nP1=4,nP2=4,Atot=5000,Dtot=5000,Ptot=5000),
fn=ls,control=list(trace=TRUE))
As you can see the objective function is composed of results from each of the three simulations per iteration. I implemented an if-statement implying that it should be returned NA in case a certain condition is fulfilled. Apparently this does not work. The algorithm just cancels. I could say it should return a numeric value that is undesirable as objective value but this could cause the function to believe it hit convergence when it occurs too often in a row.
Do you an idea how to force the algorihm to dismiss the current parameter set if a condition is fulfilled and skip to another one instead?
I hope this makes sense - best wishes and many thanks!