0

I am using the simple code as below and I am getting an error: "Error in func(time, state, parms, ...) : unused arguments (state, parms)". It looks like this error arises whenever I use 'function(koff_WT)' function otherwise the model runs fine . Any suggestions what may be going wrong ?

library(deSolve)
kon_WT = 1e-4         
kd_WT = 0.01        
koff_WT = kon_WT*kd_WT 
R_WT = 1000       

Complex <- function(koff_WT) {
 function (t,y,parms){
  with(as.list(y,parms), {
    
    dC_WT <-  koff_WT*RL_WT -kon_WT*R_WT*C_WT 
    dRL_WT <- kon_WT*R_WT*C_WT - koff_WT*RL_WT   #nM
    dR_WT <-  koff_WT*RL_WT -kon_WT*R_WT*C_WT 
        
    return(list(c(dC_WT, dRL_WT, dR_WT)))
  })
  }
}

yini <- c(C_WT = 0.1, RL_WT = 0, R_WT= R_WT)
times <- seq(0,40000,100)      
Out <- ode(y = yini, times = times, func=Complex, parms=NULL)

Output <- data.frame(Out)
Result <- Complex(koff_WT =1) ; Result
Khushhi
  • 61
  • 4

2 Answers2

2

Complex() is a factory, i.e. a function that returns a function of the type that ode is expecting. So you need to call ode with func=Complex(...) rather than func=Complex.

Out <- ode(y = yini, times = times, 
      func=Complex(koff_WT=1), parms=NULL)
Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
1

The original code had several issues, #1 already shown by Ben Bolker (I just put it to a separate line), #2 a missing c() in the with function and #3 a too big parameter value, so that the model declined immediately. #1 and #2 are programming issues, #3 a modelling decision. Finally a little bit code cosmetics to improve readability.

library(deSolve)
kon_WT  <- 1e-4         
kd_WT   <- 0.01        
koff_WT <- kon_WT * kd_WT 
R_WT    <- 1 # 1000                  # (3) R_WT was too high

Complex <- function(koff_WT) {
  function (t, y, parms){
    with(as.list(c(y,parms)), {      # (2) note c() to combine the vectors
      
      dC_WT <-  koff_WT * RL_WT     - kon_WT  * R_WT * C_WT 
      dRL_WT <- kon_WT  * R_WT*C_WT - koff_WT * RL_WT   #nM
      dR_WT <-  koff_WT * RL_WT     - kon_WT  * R_WT * C_WT 
      
      return(list(c(dC_WT, dRL_WT, dR_WT)))
    })
  }
}

yini <- c(C_WT = 0.1, RL_WT = 0, R_WT = R_WT)
times <- seq(0, 40000, 100)

model <- Complex(koff_WT) # (1) Complex is a factory (or closure), returning a function

Out <- ode(y = yini, times = times, func = model, parms = NULL)

plot(Out)
tpetzoldt
  • 5,338
  • 2
  • 12
  • 29