Inferring parameters from an ODE model

ode_model_inference.R

# Systems Biology - Université Lyon 1 
# Fall 2021
#
# Inferring parameters from a ODE model 
# --------------------------------------
#
# This script shows how to fit an ODE model to 
# a dataset.
#
# The dataset is the first subject of the Indometh dataset
# contained in the R distribution
# 
# The model is a two-compartment model for the concentration
# of indomethacin, a nonsteroidal anti-inflammatory drug (NSAD) 
# in blood. The two compartments are blood (plasma concentration) 
# and tissues (peripheral tissues). Drug can diffuse between plasma
# and peripheral tissues, and can possibly be degraded.
# The resulting ODE mode is (assuming equal volumes in each compartments)
# 
#   dys/dt <- - ks*ys - kst*ys + kts*yt # equation for y_s (blood)
#   dyt/dt <- - kt*yt + kst*ys - kts*yt # equation for y_t (tissues)
#
# The model has two dynamical variables and four ODE parameter. In addition
# the initial drug concentration `ys0` is unknown, making a total of five
# model parameter to estimate.
# Two helper functions are defined: `rhs` and `odemodel`. rhs evaluates the 
# derivatives dys/dt and dyt/dt, while odemodel numerically solve the ODEs 
# 
# The R function `nls` is used to infer model parameters by minimizing the 
# sum of square of the errors.
#
# to run the code in Rstudio, press the Source button on the upper right corner of this
# window. The first try does not work, and results in an error. This is because the 
# model has too many parameters and `nls` cannot choose. Comment out the four lines 
# followinf the 'First try to fit the parameters' and run the code again, this should
# now work!


require(deSolve) # we need this package to solve the ODE system

# get first subject from dataset Indometh
Indo.1 <- Indometh[Indometh$Subject == 1, ]

# define the system of ODEs
rhs <- function(time, y, ode_parameters) {
  ############################################################################
  # rhs: right-hand-side of a system of ordinary differential equations (ODE)
  # rhs evaluates the derivative of the system
  # 
  #   dy/dt = f(t,y,ode_parameters)
  #
  # Arguments:
  # 
  # time    independent variable (time)
  # y       a named list or numeric vector of dynamical variables
  # ode_parameters  a named list or numeric vector of parameters
  #
  # Value:
  # 
  # A numerical list of derivatives dy1/dt, dy2/dt, etc
  ############################################################################
  with(as.list(c(y,ode_parameters)), {  
  
  ### begin user-defined equations - make changes below
  dysdt <- - ks*ys - kst*ys + kts*yt # equation for y_s (blood)
  dytdt <- - kt*yt + kst*ys - kts*yt # equation for y_t (tissues)
	return(list(c(dysdt, dytdt)))
	### end user-defined equations - make changes above
  
  })
}

odemodel <- function(times, ys0, ks, kt, kst, kts) {
  ###############################################################
  # odemodel: numerical solution to a system of ODE
  #
  # Arguments:
  # 
  # times   time points at which the solution should be computed
  # ...     variable number of model parameters
  #
  # Value:
  # 
  # a vector of the first component of the solution (y1)
  ###############################################################
  timespan <- sort(times) # if data points are not in increasing times, sort them
  addedzero <- FALSE
  ode_params <- abs(c(ks = ks, kt = kt, kst = kst, kts = kts)) # named vector of parameters
  if (timespan[1] > 0)          # if first time point is > 0, add t0 = 0
  {
    timespan <- c(0, times)
    addedzero <- TRUE
  }
  dynamical_variables <- c(ys = ys0, yt = 0) # named vector of dynamical variables
  expmodel.sol <- ode(dynamical_variables, timespan, rhs, ode_params) # find the solution of the ODE at data point
  if (addedzero)
  {
    return( expmodel.sol[-1,2])
  }
  else
  {
    return(expmodel.sol[,2])
  }
}

# plot the data
with(Indo.1, plot(time,conc))  # plot Indo.1 data: concentration vs time

# First try to fit the parameters
indometh.ode <- nls(
		conc ~ odemodel(time, ys0, ks, kt, kst, kts),  # model
   data = Indo.1,                                 # dataset
   start = list(ys0 = 5, ks = 0.5, kt = 0.5, kst = 1.1, kts = 0.1))  # initial guess for the parameters

# This doesn't work! The 'singular gradient' in the error message:
#
#  Error in nls(conc ~ odemodel(time, ys0, ks, kt, kst, kts), data = Indo.1,  : 
#                 singular gradient
#
# indicates that we have too many parameters. Most likely this is between ks and kt, so we 
# try again without ks (we set ks = 0)

# Second try 
ks <- 0 # fix ks to 0
indometh.ode <- nls(
  conc ~ odemodel(time, ys0, ks, kt, kst, kts),  # model
  data = Indo.1,                                 # dataset
  start = list(ys0 = 5, kt = 0.5, kst = 1.1, kts = 0.1))  # initial guess, without ks


# It works! summary of the results
summary(indometh.ode)

# plot the prediction
lines(seq(0,8, by = 0.1),  predict(indometh.ode, list(time = seq(0,8, by = 0.1))))