R/create_mc_samples.R

In httk: High-Throughput Toxicokinetics

#' Create a table of parameter values for Monte Carlo
#' 
#' @description
#' This is the HTTK master function for creating a data table for use with 
#' Monte Carlo methods to simulate parameter uncertainty and variabilit.
#' Each column of the output table corresponds to an HTTK model parameter and 
#' each row corresponds
#' to a different random draw (for example, different individuals when 
#' considering biological variability). This function call three different key
#' functions to simulate parameter parameter uncertainty and/or variability in
#' one of three ways. First parameters can be varied in an uncorrelated manner
#' using truncated normal distributions by the function 
#' \code{\link{monte_carlo}}. Then,
#' physiological parameters can be varied in a correlated manner according to
#' the Ring et al. (2017) (\doi{10.1016/j.envint.2017.06.004}) \emph{httk-pop} 
#' approach by the function \code{\link{httkpop_mc}}. Next, both uncertainty
#' and variability of in vitro HTTK parameters can be simulated by the function
#' \code{\link{invitro_mc}} as described by Wambaugh et al. (2019) 
#' (\doi{10.1093/toxsci/kfz205}). Finally, tissue-specific partition 
#' coefficients are predicted for each draw using the Schmitt (2008)
#' (\doi{10.1016/j.tiv.2007.09.010}) method as calibrated to \emph{in vivo}
#' data by Pearce et al. (2017) (\doi{10.1007/s10928-017-9548-7}) and 
#' implemented in \code{\link{predict_partitioning_schmitt}}.  
#'
#' @details
#' The Monte Carlo methods used here were recently updated and described by
#' Breen et al. (submitted).
#' 
#' We aim to make any function that uses chemical identifiers 
#' (name, CAS, DTXSID) 
#' also work if passed a complete vector of parameters (that is, a row from the
#' table generated by this function). This allows the
#' use of Monte Carlo to vary the parameters and therefore vary the function
#' output. Depending on the type of parameters (for example, physiological vs. 
#' in vitro measurements) we vary the parameters in different ways with 
#' different functions.
#' 
#' @param chem.cas Chemical Abstract Services Registry Number (CAS-RN) -- if
#'  parameters is not specified then the chemical must be identified by either
#'  CAS, name, or DTXISD
#' @param chem.name Chemical name (spaces and capitalization ignored) --  if
#'  parameters is not specified then the chemical must be identified by either
#'  CAS, name, or DTXISD
#' @param dtxsid EPA's DSSTox Structure ID (\url{https://comptox.epa.gov/dashboard})  
#'  -- if parameters is not specified then the chemical must be identified by 
#' either CAS, name, or DTXSIDs
#' @param parameters Parameters from the appropriate parameterization function
#' for the model indicated by argument model
#' @param samples Number of samples generated in calculating quantiles.
#' @param species Species desired (either "Rat", "Rabbit", "Dog", "Mouse", or
#' default "Human"). Species must be set to "Human" to run httkpop model. 
#' @param suppress.messages Whether or not to suppress output message.
#' @param model Model used in calculation: 'pbtk' for the multiple compartment
#' model,'3compartment' for the three compartment model, '3compartmentss' for
#' the three compartment steady state model, and '1compartment' for one
#' compartment model.  This only applies when httkpop=TRUE and species="Human",
#' otherwise '3compartmentss' is used.
#' @param httkpop Whether or not to use the Ring et al. (2017) "httkpop"
#' population generator. Species must be 'Human'.
#' @param invitrouv Logical to indicate whether to include in vitro parameters
#' such as intrinsic hepatic clearance rate and fraction unbound in plasma
#' in uncertainty and variability analysis
#' @param calcrb2p Logical determining whether or not to recalculate the 
#' chemical ratio of blood to plasma 
#' @param censored.params The parameters listed in censored.params are sampled
#' from a normal distribution that is censored for values less than the limit
#' of detection (specified separately for each parameter). This argument should
#' be a list of sub-lists. Each sublist is named for a parameter in
#' "parameters" and contains two elements: "CV" (coefficient of variation) and
#' "LOD" (limit of detection, below which parameter values are censored. New
#' values are sampled with mean equal to the value in "parameters" and standard
#' deviation equal to the mean times the CV.  Censored values are sampled on a
#' uniform distribution between 0 and the limit of detection. Not used with
#' httkpop model.
#' @param vary.params The parameters listed in vary.params are sampled from a
#' normal distribution that is truncated at zero. This argument should be a
#' list of coefficients of variation (CV) for the normal distribution. Each
#' entry in the list is named for a parameter in "parameters". New values are
#' sampled with mean equal to the value in "parameters" and standard deviation
#' equal to the mean times the CV. Not used with httkpop model.
#' @param return.samples Whether or not to return the vector containing the
#' samples from the simulation instead of the selected quantile.
#' @param tissue Desired steady state tissue conentration.
#' @param httkpop.dt A data table generated by \code{\link{httkpop_generate}}.
#' This defaults to NULL, in which case \code{\link{httkpop_generate}} is 
#' called to generate this table.
#' @param httkpop.generate.arg.list Additional parameters passed to 
#' \code{\link{httkpop_generate}}.
#' @param convert.httkpop.arg.list Additional parameters passed to the 
#' convert_httkpop_* function for the model.
#' @param invitro.mc.arg.list Additional parameters passed to 
#' \code{\link{invitro_mc}}.
#' @param propagate.invitrouv.arg.list Additional parameters passed to model's
#' associated in vitro uncertainty and variability propagation function
#' @param parameterize.arg.list Additional parameters passed to the 
#' parameterize_* function for the model.
#'
#' @return
#' A data table where each column corresponds to parameters needed for the 
#' specified model and each row represents a different Monte Carlo sample of
#' parameter values.
#'
#' @author Caroline Ring, Robert Pearce, and John Wambaugh
#'
#' @references 
#' Wambaugh, John F., et al. "Toxicokinetic triage for 
#' environmental chemicals." Toxicological Sciences 147.1 (2015): 55-67.
#'
#' Ring, Caroline L., et al. "Identifying populations sensitive to
#' environmental chemicals by simulating toxicokinetic variability."
#' Environment international 106 (2017): 105-118. 
#'
#' @keywords Monte-Carlo
#' 
#' @examples 
#' 
#' \donttest{
#' sample_set = create_mc_samples(chem.name = 'bisphenol a')
#' }
#'
#' @import stats
#' @import methods
#' @importFrom purrr compact 
#' @export create_mc_samples
create_mc_samples <- function(chem.cas=NULL,
                        chem.name=NULL,
                        dtxsid = NULL,
                        parameters=NULL,
                        samples=1000,
                        species="Human",
                        suppress.messages=FALSE,
                        model='3compartmentss',
                        httkpop=TRUE,
                        invitrouv=TRUE,
                        calcrb2p=TRUE,
                        censored.params=list(),
                        vary.params=list(),
                        return.samples=FALSE,
                        tissue=NULL,
                        httkpop.dt=NULL,
                        invitro.mc.arg.list=NULL,
                        httkpop.generate.arg.list=
                          list(method = "direct resampling"),
                        convert.httkpop.arg.list=NULL,
                        propagate.invitrouv.arg.list=NULL,
                        parameterize.arg.list=NULL)
{

#
#
# ERROR CHECKING AND INITIALIZATION:
#
#

# We need to describe the chemical to be simulated one way or another:
  if (is.null(chem.cas) & 
      is.null(chem.name) & 
      is.null(dtxsid) &
      is.null(parameters)) 
    stop('Parameters, chem.name, chem.cas, or dtxsid must be specified.')

  if (is.null(model)) stop("Model must be specified.")
# We need to know model-specific information (from modelinfo_[MODEL].R]) 
# to set up the solver:
  model <- tolower(model)
  if (!(model %in% names(model.list)))            
  {
    stop(paste("Model",model,"not available. Please select from:",
      paste(names(model.list),collapse=", ")))
  }

  # Column names for data.tables
  # Appease R CMD check --as-cran variable binding:
  variable <- Name  <- Parameter <- hematocrit <- this.chem <- Krbc2pu <- NULL
  Rblood2plasma <- Qgutf <- Funbound.plasma <- Qtotal.liverc <- NULL
  Qcardiacc <- Qliverf <- hepatic.bioavailability <- ..parameter.names <- NULL
  Species <- NULL
  
  #Depending on model, choose the function in HTTK that will return the default
  #HTTK parameters for this chemical
  paramfun <- model.list[[model]]$parameterize.func
  parameterize.args <- purrr::compact(c(list(chem.cas=chem.cas,
                                             chem.name=chem.name,
                                             dtxsid=dtxsid,
                                             species=species,
                                        parameters=parameters,
                                        adjusted.Funbound.plasma=FALSE, # We want the unadjusted in vitro measured value
                                        adjusted.Clint=FALSE, # We want the unadjusted in vitro measured value
                                        suppress.messages=suppress.messages),
                                        parameterize.arg.list))
  # Check to see if we need to call the parameterize_MODEL function:
  if (is.null(parameters))
  {
    # Make sure all the arguments are used by the parameterization function:
    parameterize.args <- parameterize.args[names(parameterize.args) %in% 
                                             methods::formalArgs(paramfun)]
    parameters.mean <- do.call(getFromNamespace(paramfun, "httk"),
                         args=parameterize.args)
  } else {
    if (!is.list(parameters)) stop(
"Argument \"parameters\" to create_mc_samples should be a list of model parameters.")
    parameters.mean <- parameters
  }
  # Pass all 'parameterize.args' arguments and the 'suppress.messages'
  # arguments to the 'parameterize_schmitt' function.
  args.schmitt <- parameterize.args[which(
    names(parameterize.args) %in% names(formals(fun = parameterize_schmitt))
    )]
  args.schmitt$suppress.messages <- TRUE
  # The Schmitt parameters are useful if we need to redo partitioning later, though
  # some models don't include partitioning so the function might fail:
  pschmitt <- try(do.call(parameterize_schmitt,
                      args = args.schmitt),
                      silent=TRUE)
  if (is(pschmitt,"try-error")) pschmitt <- NULL
  
  # But we don't want to overwrite any Schmitt params provided by the
  # argument parameters:
  pschmitt <- pschmitt[!(names(pschmitt)%in%names(parameters.mean))]
  # Add in the Schmitt parameters:
  parameters.mean <- c(parameters.mean, pschmitt)
  parameter.names <- names(parameters.mean)
    
# Make sure that parameters that monte_carlo samples won't be overwritten later:
  if (httkpop & any(c(names(censored.params),names(vary.params)) %in%
    model.list[[model]]$httkpop.params)) stop(paste("Parameters",
      paste(c(names(censored.params),names(vary.params))[
        c(names(censored.params),names(vary.params)) %in%
        model.list[[model]]$httkpop.params],collapse=", "), 
        "are specified to be sampled by monte_carlo and then overwritten by httkpop_mc."))
        
# Make sure that parameters that monte_carlo samples won't be overwritten later:
  if (invitrouv & (any(c(names(censored.params),names(vary.params)) %in%
    model.list[[model]]$invitro.params))) stop(paste("Parameters",
      paste(c(names(censored.params),names(vary.params))[
        c(names(censored.params),names(vary.params)) %in%
        model.list[[model]]$invitro.params],collapse=", "), 
        "are specified to be sampled by monte_carlo and then overwritten by invitro_mc."))
  
#
#
# MONTE CARLO STEP ONE
#
#

# Sample any parameters requested with the conventional sampler
# (Wambaugh et al., 2015):
  parameters.dt <- monte_carlo(parameters.mean,
                     censored.params=censored.params,
                     cv.params=vary.params,
                     samples=samples,
                     suppress.messages=suppress.messages)
                      
#
#
# MONTE CARLO STEP TWO
# httk-pop (Ring et al., 2017)
#
#

  if (httkpop==TRUE & tolower(species)=="human")
  {
    # we use purrr::compact to drop NULL values from arguments list:
    physiology.dt <- do.call(httkpop_mc, args=purrr::compact(c(list(
                       model=model,
                       samples=samples,
                       httkpop.dt=httkpop.dt),
                       httkpop.generate.arg.list)))
# Overwrite parameters specified by httk-pop:
    parameters.dt[,names(physiology.dt):=physiology.dt]
    
  # Convert the httk-pop parameters to appropriate model variables
    converthttkpopfun <- model.list[[model]]$convert.httkpop.func
    if (!is.null(converthttkpopfun))
  # we use purrr::compact to drop NULL values from arguments list:
      parameters.dt <- do.call(converthttkpopfun, args=purrr::compact(c(list(
                       parameters.dt=parameters.dt,
                       httkpop.dt=httkpop.dt),
                       convert.httkpop.arg.list)))
   } else {
    if(httkpop==TRUE) 
      warning('httkpop model only available for human and thus not used.\n\
Set species=\"Human\" to run httkpop model.')   
     this.tissuedata <- subset(tissue.data, tolower(Species)==tolower(species))
     these.vols <- subset(this.tissuedata,variable=="Vol (L/kg)")
     these.vols$Name <- paste("V",these.vols$Tissue,"c",sep="")
     for (this.name in these.vols$Name)
       if (!(this.name %in% names(parameters.dt)))
         parameters.dt[,eval(this.name):=subset(these.vols,Name==this.name)$value]
     these.flows <- subset(this.tissuedata,variable=="Flow (mL/min/kg^(3/4))")
     these.flows$Name <- paste("Q",these.vols$Tissue,"f",sep="")
     # We don't use httk::physiology.data because we want the user to be able to
     # edit the data, so we grab physiology data from the interactive environment
     these.flows$value <- these.flows$value/
       as.numeric(subset(physiology.data,Parameter=="Cardiac Output")[
       tolower(colnames(physiology.data))==tolower(species)])
     for (this.name in these.flows$Name)
       if (!(this.name %in% names(parameters.dt)))
         parameters.dt[,eval(this.name):=subset(these.flows,Name==this.name)$value]
     parameters.dt[, hematocrit:=
       as.numeric(subset(physiology.data,Parameter=="Hematocrit")[
       tolower(colnames(physiology.data))==tolower(species)])]
  }
#
#
# MONTE CARLO STEP THREE
# PERFORM MONTE CARLO ON THE IN VITRO PARAMETERS (Wambaugh et al., 2019)
#
#

# Next add chemical-specific Funbound.plasma and CLint values
# Just cbind them together for now
  if (invitrouv) 
  {
    parameters.dt <- do.call(invitro_mc,
                       args=c(list(
                         parameters.dt=parameters.dt,
                         samples=samples),
                         invitro.mc.arg.list))
  }

# CLEAN UP PARAMETER MATRIX (bug fix v1.10.1)
#
# Force pKa to NA_real_ so data.table doesn't replace everything with text
#  if (any(c("pKa_Donor","pKa_Accept") %in% names(parameters.dt)))
#  {
#    suppressWarnings(parameters.dt[, c("pKa_Donor","pKa_Accept") := NULL]) %>% .[, c("pKa_Donor","pKa_Accept") := NA_real_]
#  }

#
#
# MONTE CARLO STEP FOUR
# PROPAGATE ANY CHANGES IN PARAMETER VALUES:
#
#

# Do we need to calculate first pass metabolism for this model (i.e., is flow
# from the gut into the liver not included in the model)
  firstpass <- model.list[[model]]$firstpass

  if (calcrb2p & any(c(!is.null(chem.cas),
                       !is.null(chem.name),
                       !is.null(dtxsid))))
  {
    Rb2p.invivo <- get_rblood2plasma(chem.cas=chem.cas,
                                     chem.name=chem.name,
                                     dtxsid=dtxsid)
  } else Rb2p.invivo <- NA

# We need all of these parameters to recalculate values with the 
# Schmitt (2008) method:
  if (all(c("Pow","pKa_Donor","pKa_Accept") %in% colnames(parameters.dt)))
  {
  # Then check to see if we need to calculate PC's, either because we have 
  # tissue compartments:
    if (model.list[[model]]$calcpc | 
  # Or if we use the blood to plasma ratio:
      (calcrb2p & is.na(Rb2p.invivo)) | 
  # Or if we include first pass metabolism (which requires Rb2p)
      firstpass)
    { 
  #Now, predict the partitioning coefficients using Schmitt's method. The
  #result will be a list of numerical vectors, one vector for each
  #tissue-to-plasma partitioning coefficient, and one element of each vector
  #for each individual. The list element names specify which partition
  #coefficient it is, e.g. Kliver2plasma, Kgut2plasma, etc.
      PCs <- do.call(predict_partitioning_schmitt, args = list(
               parameters=parameters.dt,
               args.schmitt,
               suppress.messages=TRUE))
# Store the red blood cell to unbound plasma partition coefficient if we need
# it later:
      if (calcrb2p | firstpass) parameters.dt[, Krbc2pu:=PCs[['Krbc2pu']]]
    }
# If we can't calculate partition coefficients now (no phys-chem) but need the
# Krbc2pu partition coefficient we can back-calculate it:
  } else if ("Rblood2plasma" %in% colnames(parameters.dt) &
    (calcrb2p & is.na(Rb2p.invivo))) 
  {
# Calculate Krbc2plasma from blood:plasma ratio (if available). We use the average
# value because we want this to PC to be the same for all individuals since we 
# don't have the phys-chem to recalculate:
    parameters.dt[,Krbc2pu:=calc_krbc2pu(
      Rb2p = parameters$Rblood2plasma,
      Funbound.plasma = parameters$Funbound.plasma)]
  }

# If the model uses partion coefficients we need to lump each individual
# separately in case rest of body organ volumes or PCs vary:
  if (model.list[[model]]$calcpc)
  {
     lumptissues <- lump_tissues(
                      PCs,
                      parameters=parameters.dt,
                      tissuelist=model.list[[model]]$tissuelist,
                      species=species,
                      model=model
                      ) 
                      
     parameters.dt[, names(lumptissues):= lumptissues]
  }
  
  if (calcrb2p | firstpass)
  {
# If we have an in vivo value, then back-calculate the partition coefficient
# from the measured value:
    if (!is.na(Rb2p.invivo))
    {
# From Pearce et al. (2017):
      parameters.dt[, Krbc2pu:=calc_krbc2pu(Rb2p.invivo,
                                            parameters.mean$Funbound.plasma,
                                            parameters.mean$hematocrit)]
    } 
# Calculate Rblood2plasma based on hematocrit, Krbc2plasma, and Funbound.plasma. 
# This is the ratio of chemical in blood vs. in plasma.
    parameters.dt[,Rblood2plasma := calc_rblood2plasma(
                                      hematocrit=hematocrit,
                                      Krbc2pu=Krbc2pu,
                                      Funbound.plasma=Funbound.plasma)] 
                                      
    if (any(is.na(parameters.dt$Rblood2plasma)))
    {
      parameters.dt[is.na(Rblood2plasma),
                          Rblood2plasma := available_rblood2plasma(
                            chem.cas=chem.cas,
                            chem.name=chem.name,
                            dtxsid=dtxsid,
                            species=species,
                            adjusted.Funbound.plasma=TRUE,
                            suppress.messages=suppress.messages)]
    }
  }
  
  if (firstpass)
  {
  
# For models that don't described first pass blood flow from the gut, need the
# total liver blood flow to cacluate a hepatic bioavailability (Rowland, 1973):
    if (!("Qtotal.liverc" %in% names(parameters.dt)))
      parameters.dt[, Qtotal.liverc:=Qcardiacc*(Qgutf+Qliverf)] # L/h/(kgBW)^3/4
  
# For models that don't described first pass blood flow from the gut, need the
# unscaled hepatic clearance to cacluate a hepatic bioavailability 
# (Rowland, 1973):      
  cl <- calc_hep_clearance(parameters=parameters.dt,
          hepatic.model='unscaled',
          suppress.messages=TRUE)#L/h/kg body weight

# we use purrr::compact to drop NULL values from arguments list:
  parameters.dt[,hepatic.bioavailability := do.call(calc_hep_bioavailability,
    args=purrr::compact(list(
      parameters=list(
        Qtotal.liverc=parameters.dt$Qtotal.liverc, # L/h/kg^3/4
        Funbound.plasma=parameters.dt$Funbound.plasma,
        Clmetabolismc=cl, # L/h/kg
        Rblood2plasma=parameters.dt$Rblood2plasma,
        BW=parameters.dt$BW),
      restrictive.clearance=parameterize.arg.list[["restrictive.clearance"]])))]
  }
  
#
# Do any model-specific uncertainty propagation
#
  propagateuvfun <- model.list[[model]]$propagateuv.func
  if (!is.null(propagateuvfun))
    parameters.dt <- do.call(propagateuvfun, args=c(list(
                       parameters.dt=parameters.dt),
                       propagate.invitrouv.arg.list))
  
# set precision:
  cols <- colnames(parameters.dt)
  parameters.dt[ , (cols) := lapply(.SD, set_httk_precision), .SDcols = cols]
  
#Return only the HTTK parameters for the specified model. That is, only the
#columns whose names are in the names of the default parameter set.
  return(parameters.dt[,model.list[[model]]$param.names,with=FALSE])
}