Function definitions

Like in Template 1, we define a function that returns a named vector with initial state values:

InitLogistic5 <- function(N = 5) {
  # Gather function arguments in a named vector
  y <- c(N = N)
  
  # Return the named vector
  return(y)
}

Like in Template 1, we define a function that returns a named vector with parameter values:

ParmsLogistic5 <- function(r = 0.3, # intrinsic population growth rate
                           K = 100  # carrying capacity
                           ) {
  # Gather function arguments in a named vector
  y <- c(r = r,
         K = K)
  
  # Return the named vector
  return(y)
}

Like in Template 1, we define a function that returns a list with the rates of change of the state variables:

RatesLogistic5 <- function(t, y, parms) {
  # use the with() function to be able to access the parameters and states easily
  # remember that the with() function is with(data, expr), where
  # - 'data' is ideally a list with named elements
  #   thus: 'y' and 'parms' are combined using function c and converted using function as.list
  # - 'expr' is an expression (i.e. code) that is evaluated
  #   this can span multiple lines of code when embraced with curly brackets {}
  with(
    as.list(c(y, parms)),{
      
      # Optional: forcing functions: get value of the driving variables at time t (if any)
      # see template 2
      
      # Optional: auxiliary equations

      # Rate equations
      dNdt <- r*N*(1-N/K)
      
      # Gather all rates of change in a vector
      # - the rates should be in the same order as the states (as specified in 'y')
      # - it can be a named vector, but does not need to be
      RATES <- c(dNdt = dNdt)
      
      # Optional: get in/out flow used to compute mass balances (or set MB <- NULL)
      # see template 4
      MB <- NULL

      # Optional: gather auxiliary variables that should be returned (or set AUX <- NULL)
      # - this should be a named vector or list!
      # here we pass the rate of change of N with respect to time back as auxiliary variable
      AUX <- RATES 
      
      # Return result as a list
      # - the first element is a vector with the rates of change (in the same order as 'y')
      # - all other elements are (optional) extra output, which should be named
      outList <- list(c(RATES, # the rates of change of the state variables
                        MB),   # the rates of change of the mass balance terms if MB is not NULL
                      AUX) # optional additional output per time step
      return(outList)
    }
  )
}

Sensitivity analysis

After defining the 3 above-mentioned functions, we we can now solve the model using the ode function when we have loaded the deSolve package, like we did in Template 1.

Here, instead of solving the model with the set parameter values, we now actually want to vary the parameters from their default values (just a tiny bit, say 0.1%), both decreasing and increasing it. We are doing this one at a time: so while varying the value of one parameter, we leave the values of all the others at their default values. Then, we want to compute what the difference of the state variable $N$ is (say after 5 as well as 50 time units), and how this relates to the difference in the parameter values.

### Settings

inits  <- InitLogistic5()           # initial states
params <- ParmsLogistic5()          # parameters
times  <- seq(from=0, to=50, by=1)  # times at which we request output
method <- "ode45"                   # integration method

sensiParms <- c("r","K")            # the parameters to compute sensitivity indices for
                                    # this can be a subset of all parameters

stateName <- "N"                    # the state for which we want the sensitivity index

changeFraction <- 0.001             # the fraction by which each parameter is changed (here thus: 0.1%)
                                    # (thus increased / decreased)

evalTime <- c(5, 50)                # time(s) for which we want to compute the sensitivity/elasticity indices

### Perform sensitivity analysis

# Create data.frame to hold the output for each state in "sensiParms" and time in "evalTime" (stateDiff)
stateDiff <- data.frame(time = evalTime)
stateDiff[,stateName]  <- NA # Column for the value of the state given the unchanged parameters
stateDiff[,sensiParms] <- NA # Column(s) for the difference in state values given changes in parameters

# Create empty vectors to hold the parameter differences (parmDiff)
parmDiff <- params[sensiParms] # this makes of copy of the named parameter vector
parmDiff[]  <- NA # This sets all elements within the vector to value NA

# Update the times vector with the value of evalTime (so that any evalTime is possible)
times <- sort(unique(c(times, evalTime)))

# In a 'for'-loop with iterator "i" (which values specfied by 'sensiParms'):
# - create 2 copies of 'params': called 'paramsLo' and 'paramsHi'
# - reduce the value of the i-th parameter in paramsLo with 'changeFraction'
# - increase the value of the i-th parameter in paramsHi with 'changeFraction'
# - solve the ODE model for both parameter sets (they are identical except for the i-th parameter)
# - get the values of the state after 'evalTime' time units and store the difference in 'stateDiff'
# - store the difference between the elevated and reduced parameter value in 'parmDiff'
for(i in sensiParms) {
  # Create copies of params
  paramsLo <- params
  paramsHi <- params
  
  # Reduce/increase the value of the i-th parameter
  paramsLo[i] <- (1 - changeFraction) * params[i]
  paramsHi[i] <- (1 + changeFraction) * params[i]
  
  # Solve the ODE model for both sets of parameters
  states_lo <- ode(y = inits,
                   times = times,
                   parms = paramsLo,
                   func = RatesLogistic5,
                   method = method)
  states_hi <- ode(y = inits,
                   times = times,
                   parms = paramsHi,
                   func = RatesLogistic5,
                   method = method)
  
  # Retrieve the values of the state variable at evalTime time units
  subsetLo <- subset(as.data.frame(states_lo), time %in% evalTime)[[stateName]]
  subsetHi <- subset(as.data.frame(states_hi), time %in% evalTime)[[stateName]]
  
  # Compute the differences and store in stateDiff
  stateDiff[,i] <- subsetHi - subsetLo
  
  # Store the difference between the elevated and reduced parameter value in the i-th element of parmDiff
  parmDiff[i] <- paramsHi[i] - paramsLo[i]
}

# Add the value of the state variable given the (unchanged!) parameter values
states <- ode(y = inits,
              times = times,
              parms = params,
              func = RatesLogistic5,
              method = method)
stateDiff[,stateName] <- subset(as.data.frame(states), time %in% evalTime)[[stateName]]

When we have run this code, we have two objects that contain the information needed to compute the sensitivity index: the change in the values of the state variables at the requested times (stored in object stateDiff), and the change in the parameter values (stored in the object parmDiff).

stateDiff

##   time        N            r          K
## 1    5 19.08590 4.632963e-02 0.00565984
## 2   50 99.99942 1.743696e-05 0.19999761

parmDiff

##     r     K 
## 6e-04 2e-01

Now we can calculate the sensitivity index, let’s add them to an object called SI:

SI <- stateDiff
for(i in sensiParms) {
  SI[,paste("SI_d",i,sep="_")] <- stateDiff[,i] / parmDiff[i]
}
SI

##   time        N            r          K      SI_d_r    SI_d_K
## 1    5 19.08590 4.632963e-02 0.00565984 77.21604667 0.0282992
## 2   50 99.99942 1.743696e-05 0.19999761  0.02906161 0.9999881

We see that object SI has multiple columns:

• time: the time at which we wanted to computed the sensitivity/elasticity index (specified in evalTime);
• state (name as specified in stateName, here “N”): value of the state variable with parameters at default values;
• parameters listed in sensiParms (here: “r” and “K”): the difference in the state when changing each parameter one at a time;
• parameters listed in sensiParms with prefix “SI_d_” (thus here “SI_d_r” and “SI_d_K”): the sensitivity index for that parameter and for the specified state at the specified time.

After 50 time steps, we see that the state variable $N$ is most sensitive to variation in the parameter $K$, and not very sensitive to parameter $r$.

However, the picture becomes very different when we look at the change in the value of the state variable after only 5 time units! Namely, the system’s state $N$ is very sensitive to parameter $r$ when evaluating it after 5 time units! It should make sense to you: $K$ is a parameter specifying the carrying capacity, thus the equilibrium value of state $N$, whereas parameter $r$ specifies the intrinsic growth rate, thus the rate with which population $N$ is approaching its carrying capacity. Thus, at long time-frames, the system should be sensitive to parameter $K$, whereas at short time-scale, the systems should be sensitive to parameter $r$.

Downloadable R script

An R script with the code chunks of the above example can be downloaded here.