2

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!

Arne
  • 337
  • 1
  • 6
  • 17
  • Haven't tried it, but maybe setting to `Inf` instead of `NA` should do the trick? – tonytonov Aug 02 '17 at 15:59
  • Hey - thanks for your reply. I tried it out. Unfortunately, it had the same effect as when using "NA". The following error message cancels the algorithm as soon as the if-condition is fulfilled: "Error in optim(par = c(kA1 = P_init[1, 1], KmA = P_init[2, 1], kA2 = P_init[3, : function cannot be evaluated at initial parameters" – Arne Aug 03 '17 at 07:35
  • Maybe the problem is [somewhere else](https://stackoverflow.com/questions/15068213/error-in-optim-function-cannot-be-evaluated-at-initial-parameters) then? Try reducing the optimization problem step by step by throwing away pieces that do not affect and you'll manage to narrow it down eventually. – tonytonov Aug 03 '17 at 08:31

0 Answers0