R/calc_mc_css.R

In httk: High-Throughput Toxicokinetics

#' Distribution of chemical steady state concentration with uncertainty and variability
#' 
#' @description
#' For a given chemical and fixed dose rate this function determines a 
#' distribution of steady-state concentrations reflecting measurment uncertainty
#' an population variability. Uncertainty and variability are simulated via the
#' Monte Carlo method -- many sets of model parameters are drawn according to
#' probability distributions described in 
#' Ring et al. (2017) (\doi{10.1016/j.envint.2017.06.004})
#' for human variability and Wambaugh et al. (2019) 
#' (\doi{10.1093/toxsci/kfz205}) 
#' for measurement uncertainty. Monte Carlo samples are generated by the 
#' function \code{\link{create_mc_samples}}. To allow rapid application of the 
#' Monte Carlo 
#' method we make use of analytical solutions for the steady-state concentration
#' for a particular model via a given route (when available) as opposed to
#' solving the model numerically (that is, using differential equations).
#' For each sample of the Monte Carlo method (as specified by argument
#' \code{samples}) the parameters for the analytical solution are varied. An
#' ensemble of steady-state predictions are produced, though by default only
#' the quantiles specified by argument \code{which.quantile} are provided. If 
#' the full set of predicted values are desired use set the argument 
#' \code{return.samples} to \code{TRUE}.

#'
#' @details
#' The chemical-specific steady-state concentration for a dose rate of 
#' 1 mg/kg body weight/day can be used for in \emph{in vitro-in vivo} 
#' extrapolation (IVIVE) of
#' a bioactive \emph{in vitro} concentration by dividing the \emph{in vitro}
#' concentration
#' by the steady-state concentration to predict a dose rate (mg/kg/day) that
#' would produce that concentration in plasma. Using quantiles of the 
#' distribution (such as the upper 95th percentile) allow incorporation of 
#' uncertainty and variability into IVIVE.
#' 
#' Reverse Dosimetry Toxicodynamic IVIVE
#' 
#' \if{html}{\figure{ivive.png}{options: width="60\%" alt="Reverse Dosimetry Toxicodynamic IVIVE"}}
#' \if{latex}{\figure{ivive.pdf}{options: width=12cm alt="Reverse Dosimetry Toxicodynamic IVIVE"}}
#' 
#' Figure from Breen et al. (2021) (\doi{10.1080/17425255.2021.1935867}) shows 
#' reverse dosimetry toxicodynamic IVIVE. Equivalent external dose is
#' determined by solving the TK model in reverse by deriving the external dose
#' (that is, TK model input) that produces a specified internal concentration 
#' (that is, TK model output). Reverse dosimetry and IVIVE using HTTK relies on 
#' the linearity of the models. We calculate a scaling factor to relate 
#' \emph{in vitro} 
#' concentrations (uM) to administered equivalent doses (AED). The scaling 
#' factor is the inverse of the steady state plasma concentration (Css) 
#' predicted for a 1 mg/kg/day exposure dose rate. We use Monte Carlo to 
#' simulate variability and propagate uncertainty; for example, to calculate an 
#' upper 95th percentile Css,95 for individuals who get higher plasma 
#' concentrations from the same exposure.
#'
#' The Monte Carlo methods used here were recently updated and described by
#' Breen et al. (submitted).
#'
#' 
#' httk-pop is used only for humans. For non-human species biological 
#' variability is simulated by drawing parameters from uncorellated log-normal
#' distributions. 
#' 
#' Chemical-specific httk data are available primarily for human and for a few
#' hundred chemicals in rat. All in silico predictions are for human. Thus, 
#' when species is specified as rabbit, dog, or mouse, the user can choose
#' to set the argument \code{default.to.human} to \code{TRUE} so that this 
#' function uses the
#' appropriate physiological data (volumes and flows) but substitutes human
#' fraction unbound, partition coefficients, and intrinsic hepatic clearance.
#' 
#' If the argument \code{tissue} is used, the steady-state concentration in that
#' tissue, if available, is provided. If that tissue is included in the model
#' used (specified by arguement \code{model}) then the actual tissue 
#' concentration is
#' provided. Otherwise, the tissue-specific partition coefficient is used to
#' estimate the concentration from the plasma. 
#' 
#' The six sets of plausible IVIVE
#' assumptions identified by Honda et al. (2019) 
#' (\doi{10.1371/journal.pone.0217564}) 
#' are: \tabular{lrrrr}{
#' \tab \emph{in vivo} Conc. \tab Metabolic Clearance \tab Bioactive Chemical
#' Conc. \tab TK Statistic Used* \cr Honda1 \tab Veinous (Plasma) \tab
#' Restrictive \tab Free \tab Mean Conc. \cr Honda2 \tab Veinous \tab
#' Restrictive \tab Free \tab Max Conc. \cr Honda3 \tab Veinous \tab
#' Non-restrictive \tab Total \tab Mean Conc. \cr Honda4 \tab Veinous \tab
#' Non-restrictive \tab Total \tab Max Conc. \cr Honda5 \tab Target Tissue \tab
#' Non-restrictive \tab Total \tab Mean Conc. \cr Honda6 \tab Target Tissue
#' \tab Non-restrictive \tab Total \tab Max Conc. \cr } *Assumption is
#' currently ignored because analytical steady-state solutions are currently
#' used by this function.
#' 
#' @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 which.quantile Which quantile from Monte Carlo simulation is
#' requested. Can be a vector.
#'
#' @param species Species desired (either "Rat", "Rabbit", "Dog", "Mouse", or
#' default "Human").  Species must be set to "Human" to run httkpop model.
#'
#' @param daily.dose Total daily dose, mg/kg BW.
#' 
#' @param suppress.messages Whether or not to suppress output message.
#'
#' @param model Model used in calculation,'gas_pbtk' for the gas pbtk model, 
#' '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 population generator and sampler from
#' httkpop.  This is overwrites censored.params and vary.params and is only for
#' human physiology.  Species must also be set to 'Human'.
#'
#' @param invitrouv Logical to indicate whether to include in vitro parameters
#' 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 sublists. 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 concentration. Default is of NULL
#' typically gives whole body plasma concentration.
#'
#' @param concentration Desired concentration type: 'blood','tissue', or default 
#' 'plasma'. In the case that the concentration is for plasma, selecting "blood"
#' will use the blood:plasma ratio to estimate blood concentration. In the case
#' that the argument 'tissue' specifies a particular tissue of the body, 
#' concentration defaults to 'tissue' -- that is, the concentration in the 
#' If cocentration is set to 'blood' or 'plasma' and 'tissue' specifies a
#' specific tissue then the value returned is for the plasma or blood in that
#' specific tissue.
#'
#' @param output.units Plasma concentration units, either uM or default mg/L.
#'
#' @param invitro.mc.arg.list List of additional parameters passed to 
#' \code{\link{invitro_mc}}
#'
#' @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 parameterize.arg.list A list of arguments to be passed to the model
#' parameterization function (that is, parameterize_MODEL) corresponding to
#' argument "model". (Defaults to NULL.)  
#'
#' @param calc.analytic.css.arg.list Additional parameters passed to 
#' \code{\link{calc_analytic_css}}.
#'
#' @author Caroline Ring, Robert Pearce, John Wambaugh, Miyuki Breen
#'
#' @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. 
#' 
#' Honda, Gregory S., et al. "Using the Concordance of In Vitro and 
#' In Vivo Data to Evaluate Extrapolation Assumptions." 2019. PLoS ONE 14(5): e0217564.
#'                                                       
#' Rowland, Malcolm, Leslie Z. Benet, and Garry G. Graham. "Clearance concepts in 
#' pharmacokinetics." Journal of pharmacokinetics and biopharmaceutics 1.2 (1973): 123-136.
#'
#' @keywords Monte-Carlo Steady-State
#'
#' @seealso \code{\link{calc_analytic_css}}
#'
#' @seealso \code{\link{create_mc_samples}}
#'
#' @return
#' Quantiles (specified by which.quantile) of the distribution of plasma
#' steady-stae concentration (Css) from the Monte Carlo simulation
#'
#' @examples
#' \donttest{
#' # Set the number of samples (NSAMP) low for rapid testing, increase NSAMP 
#' # for more stable results. Default value is 1000:
#' NSAMP = 10
#'
#' # Basic in vitro - in vivo extrapolation with httk, convert 3 uM in vitro
#' # concentration of chemical with CAS 2451-62-9 to mg/kg/day:
#' set.seed(1234)
#' 3/calc_mc_css(chem.cas="2451-62-9", samples=NSAMP, output.units="uM")
#' # The significant digits should give the same answer as:
#' set.seed(1234)
#' calc_mc_oral_equiv(chem.cas="2451-62-9", conc=3, samples=NSAMP)  
#'
#'  set.seed(1234)
#'  calc_mc_css(chem.name='Bisphenol A', output.units='uM',
#'              samples=NSAMP, return.samples=TRUE)
#' 
#'  set.seed(1234)
#'  calc_mc_css(chem.name='Bisphenol A', output.units='uM',
#'              samples=NSAMP,
#'              httkpop.generate.arg.list=list(method='vi'))
#'                           
#'  # The following example should result in an error since we do not 
#'  # estimate tissue partitioning with '3compartmentss'.                         
#'  set.seed(1234)
#'  try(calc_mc_css(chem.name='2,4-d', which.quantile=.9,
#'              samples=NSAMP,
#'              httkpop=FALSE, tissue='heart'))
#'  
#' # But heart will work with PBTK, even though it's lumped since we estimate
#' # a partition coefficient before lumping:
#'  set.seed(1234)
#'  calc_mc_css(chem.name='2,4-d', model='pbtk',
#'              samples=NSAMP,
#'              which.quantile=.9, httkpop=FALSE, tissue='heart')
#' 
#'  set.seed(1234)
#'  calc_mc_css(chem.cas = "80-05-7", which.quantile = 0.5,
#'              output.units = "uM", samples = NSAMP,
#'              httkpop.generate.arg.list=list(method='vi', gendernum=NULL, 
#'              agelim_years=NULL, agelim_months=NULL,
#'              weight_category = c("Underweight","Normal","Overweight","Obese")))
#' 
#'  params <- parameterize_pbtk(chem.cas="80-05-7")
#'  set.seed(1234)
#'  calc_mc_css(parameters=params,model="pbtk", samples=NSAMP)
#'
#'  set.seed(1234)
#'  # Standard HTTK Monte Carlo 
#'  calc_mc_css(chem.cas="90-43-7", model="pbtk", samples=NSAMP)
#'  set.seed(1234)
#'  # HTTK Monte Carlo with no measurement uncertainty (pre v1.10.0):
#'  calc_mc_css(chem.cas="90-43-7",
#'  model="pbtk",
#'  samples=NSAMP,
#'  invitro.mc.arg.list = list(
#'    adjusted.Funbound.plasma = TRUE,
#'    poormetab = TRUE, 
#'    fup.censored.dist = FALSE, 
#'    fup.lod = 0.01, 
#'    fup.meas.cv = 0.0, 
#'    clint.meas.cv = 0.0, 
#'    fup.pop.cv = 0.3, 
#'    clint.pop.cv = 0.3))
#'
#'  # HTTK Monte Carlo with no HTTK-Pop physiological variability):
#'  set.seed(1234)
#'  calc_mc_css(chem.cas="90-43-7",model="pbtk",samples=NSAMP,httkpop=FALSE)
#'
#'  # HTTK Monte Carlo with no in vitro uncertainty and variability):
#'  set.seed(1234)
#'  calc_mc_css(chem.cas="90-43-7",model="pbtk",samples=NSAMP,invitrouv=FALSE)
#'
#'  # HTTK Monte Carlo with no HTTK-Pop and no in vitro uncertainty and variability):
#'  set.seed(1234)
#'  calc_mc_css(chem.cas="90-43-7" ,model="pbtk",
#'              samples=NSAMP, httkpop=FALSE, invitrouv=FALSE)
#'
#'  # Should be the same as the mean result:
#'  calc_analytic_css(chem.cas="90-43-7",model="pbtk",output.units="mg/L")
#'
#'  # HTTK Monte Carlo using basic Monte Carlo sampler:
#'  set.seed(1234)
#'  calc_mc_css(chem.cas="90-43-7",
#'              model="pbtk",
#'              samples=NSAMP,
#'              httkpop=FALSE,
#'              invitrouv=FALSE,
#'              vary.params=list(Pow=0.3))
#' }
#'
#' @import stats
#' @importFrom purrr compact 
#' @export calc_mc_css
calc_mc_css <- function(chem.cas = NULL,
                        chem.name = NULL,
                        dtxsid = NULL,
                        parameters = NULL,
                        samples = 1000,
                        which.quantile = 0.95,
                        species = "Human",
                        daily.dose = 1,
                        suppress.messages = FALSE,
                        model = '3compartmentss',
                        httkpop = TRUE,
                        invitrouv = TRUE,
                        calcrb2p = TRUE,
                        censored.params = list(),
                        vary.params = list(),
                        return.samples = FALSE,
                        tissue = NULL,
                        concentration = "plasma",
                        output.units = "mg/L",
                        invitro.mc.arg.list = NULL,
                        httkpop.generate.arg.list = 
                          list(method = "direct resampling"),
                        convert.httkpop.arg.list = NULL,
                        parameterize.arg.list = NULL,
                        calc.analytic.css.arg.list = NULL
                        ) 
{
# 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.")
  
  #Appease R CMD check --as-cran variable binding:
  Css <- NULL
  
  # 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=", ")))
  } 
  parameterize_function <- model.list[[model]]$parameterize.func

  # Error handling for tissue argument:
  if (!is.null(tissue))
  {
    if (is.null(model.list[[model]]$alltissues))
    {
      stop(paste("Tissues are not available for model", model))
    }
    if (!(tissue %in% model.list[[model]]$alltissues))
    {
      stop(paste("Tissue", tissue, "not available for model", model))
    }
  }  
  
  # Error handling for concentration arugment:
  if (!(concentration %in% c("blood","tissue","plasma")))
  {
    stop("Concentration must be one of blood, tissue, or plasma")
  }
    
#
#
# CREATE A TABLE OF PARAMETER VALUES WHERE EACH ROW IS A SEPARATE SET OF 
# VALUES FOR WHICH Css SHOULD BE CALCULATED
#
#
  parameter.dt <- do.call(create_mc_samples,
# we use purrr::compact to drop NULL values from arguments list:
                          args=purrr::compact(c(list(
                              chem.cas=chem.cas,
                              chem.name=chem.name,
                              dtxsid = dtxsid,
                              parameters=parameters,
                              samples=samples,
                              species=species,
                              suppress.messages=suppress.messages,
                              model=model,
                              httkpop=httkpop,
                              invitrouv=invitrouv,
                              calcrb2p=calcrb2p,
                              censored.params=censored.params,
                              vary.params=vary.params,
                              return.samples=FALSE,
                              invitro.mc.arg.list=invitro.mc.arg.list,
                              httkpop.generate.arg.list=httkpop.generate.arg.list,
                              convert.httkpop.arg.list=convert.httkpop.arg.list,
                              parameterize.arg.list=parameterize.arg.list))))

#
# HERE LIES THE ACTUAL MONTE CARLO STEP:
#
# Calculate CSS for each row in the parameter matrix (each row corresponds to
# a different individual):
#

  parameter.dt[,Css:= do.call(calc_analytic_css,
# we use purrr::compact to drop NULL values from arguments list:
                              args=purrr::compact(c(list(parameters=.SD,
                              model=model,
                              suppress.messages=TRUE,
                              chem.cas=chem.cas,
                              chem.name=chem.name,
                              dtxsid=dtxsid,
                              tissue=tissue,
                              concentration=concentration,
                              output.units=output.units,
                              daily.dose=daily.dose,
                              clint.pvalue.threshold=
                                parameterize.arg.list$clint.pvalue.threshold),
                              calc.analytic.css.arg.list)))]

  css.list <- parameter.dt$Css 
  
  if (!return.samples) out <- quantile(css.list,which.quantile,na.rm=TRUE)  
  else out <- css.list  

  if (!suppress.messages)
  {
    if (!return.samples)
    {
      if (is.null(tissue)) cat(paste(toupper(substr(species,1,1)),
        substr(species,2,nchar(species)),sep=''),
        "plasma concentration returned in",
        output.units,
        "units for",
        which.quantile,
        "quantile.\n") 
        else cat(paste(toupper(substr(species,1,1)),
          substr(species,2,nchar(species)),sep=''),
          tissue,
          "concentration returned in",
          output.units,
          "units for",
          which.quantile,
          "quantile.\n")    
    } else {
      if (is.null(tissue)) cat(paste(toupper(substr(species,1,1)),
        substr(species,2,nchar(species)),sep=''),
        "plasma concentration returned in",
        output.units,
        "units.\n") 
      else cat(paste(toupper(substr(species,1,1)),
        substr(species,2,nchar(species)),sep=''),
        tissue,
        "concentration returned in",
        output.units,
        "units.\n") 
    }
  }
 
# Cannot guarantee arbitrary precision:
  out <- set_httk_precision(out)
  
  return(out)
}