How can I condition the parameter values in an ODE function? (deSolve, if-else)

Question

I'm trying to create values conditioned for some parameters given the initial values of states. For example, if D state is D >= 60, the S value will be S=1800. Otherwise, if D state is D <60 the S value will be S=4800. I used function if-else into the ode function (AedesIbag_model). When I run ode with D=70 if-else doesn't switch S parameter value. So I have not achieved to do that this works well. I apologize if my English is not very good. Thank you for any help.

library(deSolve)

AedesIbag_model<-function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dL = R*theta*S - mu_L*L - gamma_L*L - mu_w*L  
    dP = gamma_L*L - mu_P*P  - gamma_P*P - mu_w*P   
    dA = gamma_P*P - mu_A*A 
    dD = beta - alpha*D 
    if (D >= 60) { 
      S = 1800
      } else if (D < 60) {
        S = 4800
        } else if (D >= 10) {
          mu_w = 0.1
          } else if (D < 60) {
            mu_w = 0.1*100
          }
    
    return(list(c(dL, dP, dA, dD)))
  })
}

parameters  <- list(R = 0.24, theta = 10, S = 0,
                gamma_L = 0.2088, gamma_P = 0.384,
                beta = 10, mu_L = 0.0105, mu_P = 0.01041,
                mu_A = 0.091, mu_w = 0.1, alpha = 10
                )
state      <- c(L = 100, P = 60, A = 10, D = 70)
times       <- seq(0, 100, by = 0.1)

out_1 <- ode(y = state, times = times, func = AedesIbag_model, parms = parameters)
parameters

when I run my model. the parameters conditioned don't change the values. Look!!!

> parameters
$R
[1] 0.24

$theta
[1] 10

$S
[1] 0   #S value doesn't change

$gamma_L
[1] 0.2088

$gamma_P
[1] 0.384

$beta
[1] 10

$mu_L
[1] 0.0105

$mu_P
[1] 0.01041

$mu_A
[1] 0.091

$mu_w
[1] 0        #S value doesn't change

$alpha
[1] 10

I apologize if my English is not very good. Thank you for any help.

Answer 1

I'm guessing you are using pkg deSolve. Looking at your code, the first thing to notice is that trying to assess the function's behavior by printing the value of parameters is not productive. The ode function only accepts those values, but does not modify them, since R is a functional language and therefore should not modify its arguments. The next error to correct is the if-else construction. It would never modify the mu_w parameter in its current form, because either D >= 60 or D < 60 so S might get modified be never mu_w.

I don't know whether if(.){.}else{.} construction will work and started out taking your word for it that it was failing until I realized what I wrote above. So I used a different form of logical operations replacing the cascading if-else construction with Boolean math:

    S    <- (D >= 60)* 1800 + (D < 60) * 4800
    mu_w <- (D >= 10) * 0.1 + (D < 60)* 0.1*100

I also thought that the names of the function arguments to ode should match those in your target functions, so should be parms rather than parameters, but I didn't see much difference in the behavior when I made that change, so maybe arguments are passed by position rather than by name. If you want to see how the result of the ode call evolves, it is more effective to to plot the results:

 png(); matplot(out_1, pch=1:3)
       legend("topright", 4, unlist(dimnames(out_1)[2]),
              pch = 1:5, col = 1:5)
 dev.off()

A further note on debugging: Since the value of parameters is not changed, you would need to put a print or cat statement inside the defined function to monitor changes in the local environments values for the named parameters.

Using packageDescription('deSolve') one learns that there is a gitHub hosted webpage that has tutorials. One further learns at that page that there is a help page on the topic ?forcings. I can highly recommend the very readable text, "Solving Differential Equations in R". There is also a link to the dynamic models mailing list: mailing list: https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models

After a brief look at the material on the tutorial and `?forcings" my suspicion is that the package authors are advising a piecewise linear function rather than a discontinuous function for parameters that have a regime change.