How to set time-varying variable in the model block?

[OmarAshkar]

Hello,

I am trying to define a time-varying variable, D1, in my model with if-else as follows:

mod <- function(){
  
  ini({
    Ts_init = 5500
    TR1_init = 0 
    k_g1 = 0.03
    k_s1 = 0.1
    k_M1 = 0.05
    k_M2 = 0.03
  })
  
  model({
    if(t < 18){
      D1 = 1
      D2 = 0
    } else{
      D1 = 0
      D2 = 0
    }

    
    Ts(0) = Ts_init
    d/dt(Ts) = k_g1*Ts - k_s1*D1*Ts - k_M1*D1*Ts + k_M2*(1-D1)*TR1_init   
  })
}

ev <- et(0:240)

x <- rxSolve(mod, ev)

plot(x)

However, this is not correct. and I don't know why D1 and D2 are showing in the plot.

I tried removing the if else, and It worked as I expected the model to behave.

mod <- function(){
  
  ini({
    Ts_init = 5500
    TR1_init = 0 
    k_g1 = 0.03
    k_s1 = 0.1
    k_M1 = 0.05
    k_M2 = 0.03
  })
  
  model({

    D1 = 1
    D2 = 0
    
    Ts(0) = Ts_init
    d/dt(Ts) = k_g1*Ts - k_s1*D1*Ts - k_M1*D1*Ts + k_M2*(1-D1)*TR1_init
    
    
  })
}

ev <- et(0:240)

x <- rxSolve(mod, ev)

plot(x)

But still, I want this behavior where I use something like the if-else and give me that second simulation.

I would appreciate clarifying that behavior.

Thank you!

[mattfidler]

Hi @OmarAshkar

I'm unsure what you are trying to model, but the model output is correct as written. You can use plot(x, Ts) to focus on the Ts endpoint. You can select whatever item you wish to plot in the output this way. The D1 and D2 are both plotted because they are output in the solve.

To show the change I have focused on time 0-48 and only shown the Ts endpoint as follows:

library(rxode2)
#> Warning: package 'rxode2' was built under R version 4.3.3
#> rxode2 2.1.3 using 4 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`

mod <- function(){

  ini({
    Ts_init = 5500
    TR1_init = 0
    k_g1 = 0.03
    k_s1 = 0.1
    k_M1 = 0.05
    k_M2 = 0.03
  })

  model({
    if(t < 18){
      D1 = 1
      D2 = 0
    } else{
      D1 = 0
      D2 = 0
    }

    Ts(0) = Ts_init
    d/dt(Ts) = k_g1*Ts - k_s1*D1*Ts - k_M1*D1*Ts + k_M2*(1-D1)*TR1_init
  })
}

ev <- et(0:48)

x <- rxSolve(mod, ev)
#> using C compiler: 'gcc.exe (GCC) 12.3.0'

plot(x, Ts)

^{Created on 2024-06-29 with reprex v2.1.0}

[mattfidler]

The growth starts when D1 equals zero because it becomes a first order growth with k_g1*Ts. If this is not intended, you need to adjust your ODEs.