deSolve to denim

Migrate deSolve code to denim

Original code in deSolve

The model used for demonstrating the process of migrating code from deSolve to denim is as followed

# --- Model definition in deSolve
transition_func <- function(t, state, param){
  with(as.list( c(state, param) ), {
      dS = -beta*S*I/N
      dI1 = beta*S*I/N - rate*I1
      dI2 = rate*I1 - rate*I2
      dI =  dI1 + dI2
      dR = rate*I2
      list(c(dS, dI, dI1, dI2, dR))
  })
}

# ---- Model configuration 
parameters <- c(beta = 0.3, rate = 1/3, N = 1000) 
initialValues <- c(S = 999, I = 1, I1 = 1, I2=0, R=0)

# ---- Run simulation
times <- seq(0, 100) # simulation duration
ode_mod <- ode(y = initialValues, times = times, parms = parameters, func = transition_func)

# --- show output
ode_mod <- as.data.frame(ode_mod)
head(ode_mod[-1, c("time", "S", "I", "R")])

##   time        S        I          R
## 2    1 998.6561 1.294182 0.04969647
## 3    2 998.2239 1.594547 0.18152864
## 4    3 997.6985 1.921416 0.38010527
## 5    4 997.0695 2.290735 0.63971843
## 6    5 996.3225 2.716551 0.96093694
## 7    6 995.4387 3.212577 1.34869352

Model definition

Unlike deSolve where transitions between compartments are defined by a system of ODEs, transitions in denim must be defined by (i) a distribution of dwell-time (either parametric or non-parametric) or (ii) a mathematical expression.

User must first identify the transitions that best describe the ones in their deSolve model.

Identify distribution

# --- Model definition in deSolve
transition_func <- function(t, state, param){
  with(as.list( c(state, param) ), {
      
      # For S -> I transition, since it involves parameters (beta, N), 
      # the best transition to describe this is using a mathematical formula
      dS = -beta*S*I/N
      
      # For I -> R transition, linear chain trick is applied --> implies Erlang distributed dwell time
      # Hence, we can use d_gamma from denim
      dI1 = beta*S*I/N - rate*I1
      dI2 = rate*I1 - rate*I2
      dI =  dI1 + dI2
      dR = rate*I2
      list(c(dS, dI, dI1, dI2, dR))
  })
}

With the transitions identified, user can then define the model in denim.

When using denim DSL, the model structure is given as a set of key-value pairs where

key is a string showing the transition direction between compartments
value is the built-in distribution function that describe the transition

# --- Model definition in denim
transitions <- denim_dsl({
  S -> I = beta * S * I/N * timeStep
  # shape is 2 from number of I sub compartments
  I -> R = d_gamma(rate = 1/3, shape = 2) 
})

Note that when converting an ODE to a math expression, we must multiply the original expression with time step duration (variable timeStep in the model definition code). This is because we are essentially using Euler’s method to estimate the solution for the ODE.

Model configurations

Similar to deSolve, denim also ask users to provide the initial values and any additional parameters in the form of named vectors.

For the example deSolve code, while users can use the initalValues from the deSolve code as is (denim will ignore unused I1, I2 compartments as these sub-compartments will be automatically computed internally), it is recommended to remove redundant compartments (in this example, I1 and I2).

For parameters, since rate is already defined in the distribution functions, users only need to keep beta and N from the initial parameters vector. We do not need to specify the value for timeStep variable as this is a special variable in denim and will be defined later on.

# remove I1, I2 compartments
denim_initialValues <- c(S = 999, I = 1, R=0)
denim_parameters <- c(beta = 0.3, N = 1000)

Initialization of sub-compartments: when there are multiple sub-compartments (e.g., compartment I consist of I1 and I2 sub-compartments), the initial population is always assigned to the first sub-compartment. In our example, since I = 1, denim will assign I1 = 1 and I2 = 0.

There is also an option to distribute initial value across sub-compartments based on the specified distribution. To do this, simply set dist_init parameter of distribution function to TRUE.

transitions <- denim_dsl({
  S -> I = beta * S * I/N * timeStep
  I -> R = d_gamma(rate = 1/3, shape = 2, dist_init = TRUE) 
})

However, for comparison purposes, we will keep this option FALSE for the remaining of this demonstration.

Simulation

Lastly, users need to define the simulation duration and time step for denim to run. Unlike deSolve which takes a time sequence, denim only require the simulation duration and time step.

Since denim is a discrete time model, time step must be set to a small value for the result to closely follow that of deSolve (in this example, 0.01).

Note that the value for parameter timeStep is also the value for variable timeStep in the mathematical equation.

mod <- sim(transitions = transitions,
             initialValues = denim_initialValues, 
             parameters = denim_parameters,
             simulationDuration = 100,
             timeStep = 0.01)

head(mod[mod$Time %% 1 == 0, ])