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
#'   \insertCite{ring2017identifying;textual}{httk}
#'   (\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
#'   \insertCite{wambaugh2019assessing;textual}{httk}
#'   (\doi{10.1093/toxsci/kfz205}). Finally, tissue-specific partition
#'   coefficients are predicted for each draw using the
#'   \insertCite{schmitt2008general;textual}{httk}
#'   (\doi{10.1016/j.tiv.2007.09.010}) method as calibrated to \emph{in vivo}
#'   data by \insertCite{pearce2017evaluation;textual}{httk}
#'   (\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 \insertCite{breen2022simulating;textual}{httk}.
#'
#'   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.
#'
#'   NOTE: This function calculates oral bioavailability parameters based on
#'   Caco2 data. If you are comparing the output of this function to the output
#'   of a model-parameterization function `parameterize_MODEL()`, you will need
#'   to to ensure that `parameterize_MODEL()` also uses Caco2 data. The built-in
#'   model parameterization functions default to using *in vivo* oral
#'   bioavailability data when available. To force them to use Caco2 data
#'   instead, please pass the following argument to `parameterize_MODEL()`:
#'   `Caco2.options = list(overwrite.invivo = TRUE`.
#'
#' @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
#'   \insertCite{ring2017identifying;textual}{httk} "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.args.list Additional parameters passed to the
#'   parameterize_* function for the model.
#'
#' @param Caco2.options Arguments describing how to handle Caco2 absorption data
#'   that are passed to \code{\link{invitro_mc}} and the parameterize_[MODEL]
#'   functions. See \code{\link{get_fbio}} for further details.
#'
#' @param adjusted.Funbound.plasma Deprecated argument -- use
#'   parameterize.args.list
#'
#' @param adjusted.Clint Deprecated argument -- use parameterize.args.list
#'
#' @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
#'
#' \insertAllCited{}
#'
#' @keywords Monte-Carlo
#'
#' @examples
#'
#' \donttest{
#' # We can use the Monte Carlo functions by passing a table
#' # where each row represents a different Monte Carlo draw of parameters:
#' p <- create_mc_samples(chem.cas="80-05-7")
#'
#' # Use data.table for steady-state plasma concentration (Css) Monte Carlo:
#' calc_mc_css(parameters=p)
#'
#' # Using the same table gives the same answer:
#' calc_mc_css(parameters=p)
#'
#' # Use Css for 1 mg/kg/day for simple reverse toxicokinetics
#' # in Vitro-In Vivo Extrapolation to convert 15 uM to mg/kg/day:
#' 15/calc_mc_css(parameters=p, output.units="uM")
#'
#' # Can do the same with calc_mc_oral_equiv:
#' calc_mc_oral_equiv(15, parameters=p)
#'
#' #Generate a population using the virtual-individuals method,
#' #including 80 females and 20 males,
#' #including only ages 20-65,
#' #including only Mexican American and
#' #Non-Hispanic Black individuals,
#' #including only non-obese individuals
#' set.seed(42)
#' mypop <- httkpop_generate(method = 'virtual individuals',
#'                           gendernum=list(Female=80,
#'                           Male=20),
#'                           agelim_years=c(20,65),
#'                           reths=c('Mexican American',
#'                           'Non-Hispanic Black'),
#'                           weight_category=c('Underweight',
#'                           'Normal',
#'                           'Overweight'))
#' # Including a httkpop.dt argument will overwrite the number of sample and
#' # the httkpop on/off logical switch:
#' samps1 <- create_mc_samples(chem.name="bisphenola",
#'                            httkpop=FALSE,
#'                            httkpop.dt=mypop)
#' samps2 <- create_mc_samples(chem.name="bisphenola",
#'                            httkpop.dt=mypop)
#' # But we can turn httkpop off altogether if desired:
#' samps3 <- create_mc_samples(chem.name="bisphenola",
#'                            httkpop=FALSE)
#' }
#'
#' @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,
                        adjusted.Funbound.plasma=NA,
                        adjusted.Clint=NA,
                        httkpop.generate.arg.list=
                          list(method = "direct resampling"),
                        convert.httkpop.arg.list=NULL,
                        propagate.invitrouv.arg.list=NULL,
                        parameterize.args.list =NULL,
                        Caco2.options=NULL)
{

  ## Setting up binding for Global Variables ##
  Fabs <- Fgut <- Fabsgut <- NULL
  unadjusted.Funbound.plasma <- Funbound.plasma.adjustment  <- NULL
  Fhep.assay.correction <- Clint <- NULL
  ####
  
#
#
# ERROR CHECKING AND INITIALIZATION:
#
#

  # These arguments are deprecated, provide warning for now:
  if (!is.na(adjusted.Funbound.plasma))
  {
    warning("Argument adjusted.Funbound.plasma to calc_mc_samples is deprecated. Set this using parameterize.args.list")
    parameterize.args.list[["adjusted.Funbound.plasma"]] <- adjusted.Funbound.plasma 
  }
  if (is.null(parameterize.args.list[["adjusted.Funbound.plasma"]]))
    parameterize.args.list[["adjusted.Funbound.plasma"]] <- TRUE
  if (!is.na(adjusted.Clint))
  {
    warning("Argument adjusted.Clint to calc_mc_samples is deprecated. Set this using parameterize.args.list")
    parameterize.args.list[["adjusted.Clint"]] <- adjusted.Funbound.plasma 
  }
  if (is.null(parameterize.args.list[["adjusted.Clint"]]))
    parameterize.args.list[["adjusted.Clint"]] <- TRUE
    
 # If user supplied httkpot.dt make sure our MC simulation matches:
 if (!is.null(httkpop.dt))
 {
   if (samples != dim(httkpop.dt)[1])
   {
     samples <- dim(httkpop.dt)[1]
     warning(paste0("samples set to ",
                    dim(httkpop.dt)[1],
                    " to match number of rows in httkpop.dt supplied as argument"
                    ))
   }
   if (!httkpop)
   {
     httkpop <- TRUE
     warning("httkpop set to TRUE because httkpop.dt supplied as argument.")
   }
 }

# 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=", ")))
  }
  if (is.null(model.list[[model]]$monte.carlo))
    stop(paste("Model", model, "does not yet work with Monte Carlo sampling"))
  if (!model.list[[model]]$monte.carlo)
    stop(paste("Model", model, "does not yet work with Monte Carlo sampling"))
    
  # 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
  tmp_tissuelist <- as.list(
    setdiff(
      model.list[[model]]$alltissues,
      "red blood cells")
  )
  names(tmp_tissuelist) <-   setdiff(
    model.list[[model]]$alltissues,
    "red blood cells")
  parameterize.args.list<- purrr::compact(c(list(chem.cas=chem.cas,
                                             chem.name=chem.name,
                                             dtxsid=dtxsid,
                                             species=species,
                                        parameters=parameters,
                                        suppress.messages=suppress.messages,
                                        tissuelist = tmp_tissuelist),
                                        parameterize.args.list))
  if (!is.null(Caco2.options)) parameterize.args.list[["Caco2.options"]] <- Caco2.options
  

  mean.args.list <- parameterize.args.list
  mean.args.list[["adjusted.Funbound.plasma"]] <- FALSE # We want the unadjusted in vitro measured value for the mean
  mean.args.list[["adjusted.Clint"]] <- FALSE # We want the unadjusted in vitro measured value for the mean

  # 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:
    mean.args.list <- mean.args.list[names(mean.args.list) %in% 
                                             methods::formalArgs(paramfun)]
    parameters.mean <- do.call(getFromNamespace(paramfun, "httk"),
                         args=purrr::compact(mean.args.list))
  } 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 <- mean.args.list[which(
    names(mean.args.list) %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=purrr::compact(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 = parameters.mean,
                     censored.params=censored.params,
                     cv.params=vary.params,
                     samples=samples,
                     suppress.messages=suppress.messages)
  
  if(('Funbound.plasma' %in% names(vary.params)) &
     !('unadjusted.Funbound.plasma' %in% names(vary.params))){
    parameters.dt[, unadjusted.Funbound.plasma := Funbound.plasma]
  }
                      
#
#
# 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,
                       physiology.dt=physiology.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) 
  {
    
    arglist <- c(
      list(
        parameters.dt=parameters.dt,
        samples=samples
        ),
#        adjusted.Clint = adjusted.Clint,
#        adjusted.Funbound.plasma = adjusted.Funbound.plasma),
      Caco2.options[names(Caco2.options)[
        !names(Caco2.options) %in%
          # invitro_mc doesn't make use of these arguments because we've already called
          # parameterize_[MODEL]:
          c("Caco2.Pab.default", "overwrite.invivo")]],
      invitro.mc.arg.list
    )
    
    parameters.dt <- do.call(invitro_mc,
                       args=purrr::compact(
                         arglist
                       )
    )
  } else {
  # Clint and fup are adjusted within invitro_mc, if not called we need to adjust them here:
  #
  # Check to see if we are adjusting for differences between in vitro and 
  # physiological lipid partitioning (Pearce, 2017):
    if (parameterize.args.list[["adjusted.Funbound.plasma"]])
    {
      if (!"unadjusted.Funbound.plasma" %in% colnames(parameters.dt))
      {
        parameters.dt[, unadjusted.Funbound.plasma:=Funbound.plasma]
      }
      if (all(c("Pow","pKa_Donor","pKa_Accept") %in% names(parameters.dt)) | 
          ("Dow74" %in% names(parameters.dt)))
      {
        # need a default value:
        if (is.null(parameterize.args.list$default.to.human))
          parameterize.args.list[["default.to.human"]] <- FALSE
        # put the unadjusted fup where calc_fup_correction will look for it:
        parameters.dt[, Funbound.plasma := unadjusted.Funbound.plasma]
        parameters.dt[, Funbound.plasma.adjustment :=
                        do.call(calc_fup_correction,
                                args = purrr::compact(
                                  list(
                                    parameters = parameters.dt,
                                    species = species,
                                    default.to.human = parameterize.args.list$default.to.human
                                  )
                                ))]
        
        parameters.dt[, Funbound.plasma := 
          apply_fup_adjustment(fup = Funbound.plasma,
                      fup.correction = Funbound.plasma.adjustment)]
      } else stop("Missing phys-chem parameters in invitro_mc for calc_fup_correction.") 
    } else {
      parameters.dt[, Funbound.plasma.adjustment:=1]
    }
  
  # Correct for fraction of chemical unbound in in vitro hepatocyte assay:
    if (parameterize.args.list[["adjusted.Clint"]])
    {
      parameters.dt[, Clint := apply_clint_adjustment(
                                 Clint = Clint,
                                 Fu_hep=Fhep.assay.correction,
                                 suppress.messages=TRUE)]
    }
  }
 

# 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_args <- purrr::compact(
      list(chem.cas = chem.cas,
           chem.name = chem.name,
           dtxsid = dtxsid,
           species = species,
           default.to.human = parameterize.args.list$default.to.human)
    )
    
    Rb2p.invivo <- do.call(what = get_rblood2plasma,
                           args = rb2p_args)
    
  } 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=purrr::compact(list(
               parameters=parameters.dt,
               args.schmitt,
               tissues=model.list[[model]]$alltissues,
               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 PC to be the same for all individuals since we 
# don't have the phys-chem to recalculate. But we need the adjusted values:
    parameterize.args.list$adjusted.Funbound.plasma <- TRUE
    parameterize.args.list$suppress.messages=TRUE
    if (is.null(parameters))
    {
      adj.parameters.mean <- do.call(getFromNamespace(paramfun, "httk"),
                           args=purrr::compact(parameterize.args.list))
    } else adj.parameters.mean <- parameters
    
    
    parameters.dt[,Krbc2pu:=do.call(calc_krbc2pu,
                                    args = purrr::compact(
                                      list(
                                        Rb2p = adj.parameters.mean$Rblood2plasma,
                                        Funbound.plasma = adj.parameters.mean$Funbound.plasma,
                                        species = species,
                                        default.to.human = parameterize.args.list$default.to.human    
                                      )
                                    )
   )]
  }

# If the model uses partion coefficients we need to lump each individual
# separately in case rest of body organ volumes or PCs vary:
# But do not lump tissues twice (if 'Vrestc' and 'Krest2pu' both already exists in parameters.dt)
  if (model.list[[model]]$calcpc & !all(c("Vrestc", "Krest2pu") %in% names(parameters.dt)))
  {
     lumptissues <- lump_tissues(
       Ktissue2pu.in = 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:=do.call(
        calc_krbc2pu,
        args = purrr::compact(
          list(Rb2p = Rb2p.invivo,
                                            Funbound.plasma = Funbound.plasma,
                                            hematocrit = hematocrit,
                                            species = species,
                                            default.to.human = parameterize.args.list$default.to.human)
          ))]
    } 
# Calculate Rblood2plasma based on hematocrit, Krbc2plasma, and Funbound.plasma. 
# This is the ratio of chemical in blood vs. in plasma.
    parameters.dt[,Rblood2plasma := do.call(calc_rblood2plasma,
                                            args = purrr::compact(
                                              list(
                                                hematocrit=hematocrit,
                                                Krbc2pu=Krbc2pu,
                                                Funbound.plasma=Funbound.plasma,
                                                species = species,
                                                default.to.human = parameterize.args.list$default.to.human,
                                                # We can set this to TRUE because the value in Funbound.plasma is either adjusted
                                                # or not adjusted already:
                                                adjusted.Funbound.plasma=TRUE)
                                            ))]
                             
                                      
    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,
# We can set this to TRUE because the value in Funbound.plasma is either adjusted
# or not adjusted already:                           
adjusted.Funbound.plasma=TRUE,
                            suppress.messages=suppress.messages)]
    }
  }

# Calculate fraction of oral dose surviving first-pass hepatic metabolism before
# reaching systemic blood (if needed):
  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 <- do.call(calc_hep_clearance,
                args = purrr::compact(
                  list(
                    parameters=parameters.dt,
                    species = species,
          hepatic.model='unscaled',
          restrictive.clearance=parameterize.args.list[["restrictive.clearance"]],
          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(
      species = species,
      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),
      default.to.human = parameterize.args.list$default.to.human,
      restrictive.clearance=parameterize.args.list[["restrictive.clearance"]])))]
  }
  
  # If Caco2.options given use those, otherwise use defaults:
  if ("keepit100" %in% names(Caco2.options)) keepit100 <- Caco2.options[["keepit100"]]
  else keepit100 <- FALSE
  if ("Caco2.Fgut" %in% names(Caco2.options)) Caco2.Fgut <- Caco2.options[["Caco2.Fgut"]]
  else Caco2.Fgut <- TRUE
  if ("Caco2.Fabs" %in% names(Caco2.options)) Caco2.Fabs <- Caco2.options[["Caco2.Fabs"]]
  else Caco2.Fabs <- TRUE
  #
  # Propagate any parameter changes into oral bioavailability:
  # (Note that we do not need to include Fhep since that is already accounted for
  # by first pass metabolism)
  #
  if (!keepit100) 
  {
    bioavail <- calc_fbio.oral(parameters = parameters.dt,
                               species = species) 
    if (Caco2.Fabs) parameters.dt[,Fabs:=
                                   bioavail$fabs.oral]
    if (Caco2.Fgut) parameters.dt[,Fgut:=
                                   bioavail$fgut.oral]
    parameters.dt[,Fabsgut:=Fabs*Fgut]
  }
  else {
    parameters.dt[,Fabs:=1]
    parameters.dt[,Fgut:=1]
    parameters.dt[,Fabsgut:=1]
  }
#
# Do any model-specific uncertainty propagation
#
  propagateuvfun <- model.list[[model]]$propagateuv.func
  if (!is.null(propagateuvfun))
    parameters.dt <- do.call(propagateuvfun, args=purrr::compact(c(list(
                       parameters.dt=parameters.dt,
                       species = species),
                       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])
}