About Me

MS in Safety, Health and Environmental Engineering @ National United University

PhD in Bioenvironmental Systems Engineering @ National Taiwan University

Content

1 Basic Pharmacokinetics Concepts

2 Pharmacokinetic Models

3 Computational Toolkits

4 Hands-on Exercise

Pharmacokinetics (PK) / Toxiccokinetics (TK)

Pharmacokinetics / Pharmacodynamics

Pharmacokinetics profile

https://en.wikipedia.org/wiki/Pharmacokinetics

Absorption

Movement of chemical into the body

Bioavailability ia a key factor for substance becomes available to the target tissue after administration.

Exposure routes

Different exposure route can related to different PK profile

Absorption & Distribution

https://doi.org/10.1183/09031936.00074905

Distribution

Metabolism/transformation

• Enzyme systems, etc. covered in separate class.

  • Phase I – conversion to more polar metabolite, e.g., oxidation by cytochrome P450
  • Phase II – conjugation with an endogenous substrate to increase solubility, e.g., GSH conjugation, glucuronidation

• Key issues relevant to risk assessment

  • What are the metabolites?
  • In what tissues does metabolism occur?
  • What are the rates of metabolism?
  • What is known about interspecies and intraspecies differences?

• Metabolism process

  • First order kinetics (linear; constant)
  • Michaelis-Menten Kinetics (Non-linear; concentration dependent)

Metabolism - Acetaminophen

Trichloroethylene GSH conjugation pathway

Excretion

Excretion is the removal from the body

Excretion and elimination are often confused

• Excretion is reserved for exiting the body
• Elimination is the disappearance of a chemical, including both excretion and metabolism

Elimination kinetics

Impact of TK on risk assessment

Why Pharmacokinetic Modeling

• Pharmaceutical research

• Drug development

• Health risk assessment

  .  
  .  
  .  
  

High Throughput Risk Prioritization

According to Toxic Substances Control Act (TSCA), to address thousands of chemicals, we need to use new approach methodologies (NAMs) to prioritize the existing and new chemicals for testing.

source: USEPA

Pharmacokinetic Model

What information we want to know through "prediction"?

Mathematical method

How to conduct simulation?

Intravenous administration

A_0 <- 10
t <- c(0,1,2,4,8,10)
k_e <- 0.5
A_t <- A_0*exp(-k_e*t) 
plot(t, A_t, type = "b")

Oral administration

Empirical compartmental models

Example pf empirical compartmental models

Luo, Y.S., Hsieh, N.H., Soldatow, V.Y., Chiu, W.A. and Rusyn, I., 2018. Comparative Analysis of Metabolism of Trichloroethylene and Tetrachloroethylene among Mouse Tissues and Strains. Toxicology, 409, pp.33-43.

Physiological-Based Pharmacokinetic (PBPK)

Mathematically transcribing anatomical, physiological, physical, and chemical descriptions of the phenomena involved in the complex ADME processes.

Why PBPK models?

• There are some processes that empirical models have difficulty simulating – such as inhalation exposures.

• We know a lot about the physiology/anatomy of the organisms we study and want to protect.

  • Organ/tissue sizes
  • Connections between organs/tissues
  • Blood flow rates
  • Ventilation rates
  • Metabolizing enzymes
  • Etc…

By developing a model, we can incorporate that information to help us make inferences beyond the experimental data at hand.

Building a PBPK model

• Decide on what compartments are needed

  • Target vs. non-target tissues
  • Tissues with in vivo vs. in vitro vs. no data
  • Exposure routes

• Specify the (interrelated) equations for each compartment

• Specify parameter values

  • Physiological/anatomical parameters
  • Chemical-specific parameters

• Set up the inputs (doses, exposure concentrations, etc.) and outputs (blood, tissue, exhaled breath concentrations).

• Solve using numerical differential equation solver and compare with data

Iteration may be needed for refining model and/or conducting additional experiments

Compartmental equations

• Flow in and out of “storage” compartments involve only blood flow

• Key parameters/variables

  $Q_i$ = tissue blood flow

  $P_i$ = tissue-blood partition coefficient

  $V_i$ = volume of tissue

  $C_a$ = arterial blood concentration

  $C{V}_i$ = venous blood concentration leaving tissue

  $A_i$ = amount in tissue

  $C_i$ = concentration in tissue

Basic transport equations

The model equations follow the principles of mass transport, fluid dynamics, and biochemistry in order to simulate the fate of a substance in the body.

rate in $=Q_i{C}_{art}$ & rate out $=Q_i{C}_{VT}=Q_i\cdot C_i/P_i$ (Unit: mass/time)

Change in amount = rate in – rate out

$A_i$ is amount of chemical, $Q_i$ is blood flow, $C_{art}$ incoming arterial blood concentration, $P_i$ the tissue over blood partition coefficient, and $V_i$ the volume of compartment $i$.

Different exposure routes

Computational toolkits

Multiple software packages available to solve differential equations

All of them incorporate methods for “optimizing” the model fit by adjusting parameter values.

Computational toolkits

Advantages of MCSim and R

"Free and Open Source Software" under GNU General Public License

• The freedom to run the program as you wish, for any purpose (freedom 0).
• The freedom to study how the program works, and change it so it does your computing as you wish (freedom 1). Access to the source code is a precondition for this.
• The freedom to redistribute copies so you can help others (freedom 2).
• The freedom to distribute copies of your modified versions to others (freedom 3). By doing this you can give the whole community a chance to benefit from your changes. Access to the source code is a precondition for this.

Advantages

"Parallel computing"

- Run multiple MCMC chains with multiple CPUs

- High performance (cloud) computing

source

Advantages - Reproducible research

Reproducibility and Replicability in Science (2019) http://nap.edu/25303

Conducting data analysis through R

• Implement a wide variety of statistical and graphical techniques

• Comprehensive research workflow and toolkits (e.g., RMarkdown)

• Integration with low-level language (e.g., C, C++, Fortran)

• Highly extensible through the use of user-submitted packages

• Webapp development (e.g., http://webpopix.org/shiny/ShinyExamples.html)

Image source

RStudio

Free and open-source integrated development environment

Powerful and user friendly programming interface

Designed to make it easy to write scripts

Easy to view and interact with the objects

R project with version control (e.g., git)

Support cloud computing https://rstudio.cloud/

GNU MCSim

GNU MCSim Overview

The GNU MCSim consists in two pieces, a model generator and a simulation engine:

The model generator, "mod"

• Created to facilitate structural model definition and maintenance, while keeping execution time short. You can code your model using a simplified syntax and use mod to translate it to C (model.c).

The simulation engine, "sim"

• A set of routines which are linked to your model during compilation to produce executable program (mcsim.model). After that, you can run simulations of your model under a variety of conditions, specify an associated statistical model, and perform simulations.

GNU MCSim Workflow

Types of Simulation

Simple simulation

• Straight simulations (set parameter values and initial conditions).

Used to: Model testing when building the model (e.g., mass balance)

Monte Carlo simulations

• Perform repeated (stochastic) simulations across a randomly sampled region of the model parameter space.

Used to: Check possible simulation (under given parameter distributions) results before model calibration

SetPoints simulation

• Solves the model for a series of specified parameter sets. You can create these parameter sets yourself or use the output of a previous Monte Carlo or MCMC simulation.

Used to: Posterior analysis, Local/global sensitivity analysis

Types of Simulation

Markov-chain Monte Carlo (MCMC) simulation

• Performs a series of simulations along a Markov chain in the model parameter space.
• They can be used to obtain the Bayesian posterior distribution of the model parameters, given a statistical model, prior parameter distributions and data for which a likelihood function can be computed.
• GNU MCSim can handle hierarchical statistical models as well.

MCSim-related R packages

httk: R Package for High-Throughput Toxicokinetics

MCSim (Bois and Maszle 1997) was used for converting the model equations into C code, which is used with deSolve (Soetaert et al. 2016) in solving each system of equations.

GNU MCSim model code C code deSolve package Prediction

MCSim-related R packages

pksensi: R Package for Global Sensitivity Analysis in Pharmacokinetic Modeling

Nan-Hung Hsieh, Brad Reisfeld, Weihsueh A. Chiu

pksensi implements the global sensitivity analysis workflow to investigate the parameter uncertainty and sensitivity in pharmacokinetic (PK) models, especially the physiologically based pharmacokinetic (PBPK) model with multivariate outputs.

Two types of model solver

solve_fun()
GNU MCSim model code C code deSolve package Prediction

solve_mcsim()
GNU MCSim model code Prediction

Note: solve_mcsim() is faster than solve_fun()

Disadvantages in MCSim and R

Difficult learning curve (command line interface-based)

Requires coding/programming skill

Requires installation of extra program or package

Requires "debugging"

Git

- Manage your code (models, inputs, R script)

- Transfer your file (through GitHub or GitLab)

- Collaboration work

GitHub

High Performance Research Computing

The efficient way to conduct parallel computing

Linux

The open source Unix-like operating systems for high performance research computing

Summary

• Basic pharmacokinetics

  • A,D,M,E each may affect internal dose, and differences between species or among individuals.

• Role of pharmacokinetics in risk assessment

  • Exposure / concentration prediction

• Pharmacokinetic modelling

  • Information from physiology and chemistry are common ways to characterize the concentration-time relationship.

• Computational toolkits

  • The open soure computational tools (e.g., R & GNU MCSim) are powerful and can help us conduct the model simulation and data analysis

Hands on Exercise

Task 1. Exploratory analysis of PK data (code: https://rpubs.com/Nanhung/SRP19_1)

• Learn how to use R to conduct basic analysis for PK data

Task 2. PK model development (code: https://rpubs.com/Nanhung/SRP19_2)

• Learn how to use R and MCSim to build a model

Task 3. Parameter setting and model simulation (code https://rpubs.com/Nanhung/SRP19_3)

• Understand the parameter seting in PK model and conduct the simulation

Task 4. PBPK model development (code https://rpubs.com/Nanhung/SRP19_4)

• Instead of PK model, we need to know how to build a PBPK model

Task 5. Application of PBPK model (code https://rpubs.com/Nanhung/SRP19_5)

• Here, we have a well-built PBPK model and its parameters, let's apply the model in the exposure assessment

Hands on Exercise

Task 1: Exploratory analysis of PK data

• Now, we have a theophylline PK dataset. The purpose of this exercise is to develop the simple PK model and use it to describe the PK of theophylline. First, look into the theophylline dataset. The Theoph data frame has 132 rows and 5 columns of data from an experiment on the pharmacokinetics of theophylline. Then, find the Cmax and Tmax for each individual.

head(Theoph)

##   Subject   Wt Dose Time  conc
## 1       1 79.6 4.02 0.00  0.74
## 2       1 79.6 4.02 0.25  2.84
## 3       1 79.6 4.02 0.57  6.57
## 4       1 79.6 4.02 1.12 10.50
## 5       1 79.6 4.02 2.02  9.66
## 6       1 79.6 4.02 3.82  8.58

• Plot the pharmacokinetic diagram for each individual.

Hands on Exercise

Task 2: PK model development

• Develop the non compartment model and compartment model in R and MCSim

• Non compartment model

$$C\left(t\right)=\frac{F\cdot D\cdot ka}{V(k_a-k_e)}\cdot (e^{-k_{e}t}-e^{-k_{a}t})$$

• Compartment model

$$\frac{d{A}_{gut}}{dt}=-k_a{A}_{gut}$$

$$\frac{dA}{dt}=k_a{A}_{gut}-k_e{A}_e$$ $$C=A/V$$

Hands on Exercise

Task 3: Parameter setting and simulation

• Extract the chemical information from httk package

  library(httk)
  parms <- httk::parameterize_1comp(chem.name = "theophylline")

• Use the parameters in the developed model and conduct the simulation

• Compare the difference between data and the model simulation result (Cmax, Tmax)

Hands on Exercise

Task 4: PBPK model development

• Instead of simple PK model, in this exercise we want to apply the well developed PBPK model for 1,3 butadiene.

• First, reproduce the simulation result from the published paper*.

Inputs

C_inh <- approxfun(x = c(0, 120), y=c(10,0), method = "constant", f = 0, rule = 2)

plot(C_inh(1:300), type="l", xlab = "Time (min)", ylab = "Butadiene inhaled concentration (ppm)")

Parameters and outputs

parameters <- c("BDM" = 73,            # Body mass (kg)
                "Height" = 1.6,        # Body height (m)
                "Age" = 40,            # in years
                "Sex" = 1,             # code 1 is male, 2 is female
                "Flow_pul" = 5,        # Pulmonary ventilation rate (L/min)
                "Pct_Deadspace" = 0.7, # Fraction of pulmonary deadspace
                "Vent_Perf" = 1.14,    # Ventilation over perfusion ratio
                "Pct_LBDM_wp" = 0.2,   # wp tissue as fraction of lean mass
                "Pct_Flow_fat" = 0.1,  # Fraction of cardiac output to fat
                "Pct_Flow_pp" = 0.35,  # ~ to pp
                "PC_art" = 2,          # Blood/air partition coefficient
                "PC_fat" = 22,         # Fat/blood ~
                "PC_wp" = 0.8,         # wp/blood ~
                "PC_pp" = 0.8,         # pp/blood ~
                "Kmetwp" = 0.25)       # Rate constant for metabolism

y <- c("Q_fat" = 0, # Quantity of butadiene in fat (mg)
       "Q_wp" = 0,  # ~ in well-perfused (mg)
       "Q_pp" = 0,  # ~ in poorly-perfused (mg)
       "Q_met" = 0) # ~ metabolized (mg)

Model structure

# Define the model equations
bd.model <- function(t, y, parameters) {
  with (as.list(y), {
    with (as.list(parameters), {
      # Define constants
      # Calculate flow and volumes
      # Calculate the tissue, blood, and air
      # Differentials for quantities
      # The function bd.model must return at least the derivatives
      list(c(dQ_fat, dQ_wp, dQ_pp, dQ_met), # derivatives
           c("C_ven" = C_ven, "C_art" = C_art)) # extra outputs
    }) # end with parameters
  }) # end with y
} # end bd.model

Constants

Known constants

# Known constants
Height = 1.6                            # use to calculate fraction of body fat
Age = 40                                # use to calculate fraction of body fat
Sex = 1                                 # use to calculate fraction of body fat
MW_bu = 54.0914                         # butadiene molecular weight (in grams)

Conversions from/to ppm

ppm_per_mM = 24450                       # ppm to mM under normal conditions
ppm_per_mg_per_l = ppm_per_mM / MW_bu
mg_per_l_per_ppm = 1 / ppm_per_mg_per_l

Flows and volumes

Air and blood flow

# Calculate Flow_alv from total pulmonary flow
Flow_alv = Flow_pul * (1 - Pct_Deadspace)
# Calculate total blood flow from Flow_alv and the V/P ratio
Flow_tot = Flow_alv / Vent_Perf
# Calculate fraction of body fat
Pct_BDM_fat = (1.2 * BDM / (Height * Height) - 10.8 *(2 - Sex) + 0.23 * Age - 5.4) * 0.01
# Calculate actual blood flows from total flow and percent flows
Flow_fat = Pct_Flow_fat * Flow_tot
Flow_pp = Pct_Flow_pp * Flow_tot
Flow_wp = Flow_tot * (1 - Pct_Flow_pp - Pct_Flow_fat)

Volumes

# Actual volumes, 10% of body mass (bones…) get no butadiene
Eff_V_fat = Pct_BDM_fat * BDM
Eff_V_wp = Pct_LBDM_wp * BDM * (1 - Pct_BDM_fat)
Eff_V_pp = 0.9 * BDM - Eff_V_fat - Eff_V_wp

Concentrations - tissues & blood

# Calculate the concentrations
C_fat = Q_fat / Eff_V_fat
C_wp = Q_wp / Eff_V_wp
C_pp = Q_pp / Eff_V_pp
# Venous blood concentrations at the organ exit
Cout_fat = C_fat / PC_fat
Cout_wp = C_wp / PC_wp
Cout_pp = C_pp / PC_pp
# Sum of Flow * Concentration for all compartments
dQ_ven = Flow_fat * Cout_fat + Flow_wp * Cout_wp + Flow_pp * Cout_pp
C_inh.current = C_inh(t)     # to avoid calling C_inh() twice
# Arterial blood concentration
# Convert input given in ppm to mg/l to match other units
C_art = (Flow_alv * C_inh.current * mg_per_l_per_ppm + dQ_ven) / (Flow_tot + Flow_alv / PC_art)
# Venous blood concentration (mg/L)
C_ven = dQ_ven / Flow_tot

Concentrations - air

# Alveolar air concentration (mg/L)
C_alv = C_art / PC_art
# Exhaled air concentration (ppm)
if (C_alv <= 0) {
  C_exh = 10E-30    # avoid round off errors
} else {
  C_exh = (1 - Pct_Deadspace) * C_alv * ppm_per_mg_per_l + Pct_Deadspace * C_inh.current
}

Differentials for quantities

# Quantity metabolized in liver (included in well-perfused)
dQmet_wp = Kmetwp * Q_wp
# Differentials for quantities
dQ_fat = Flow_fat * (C_art - Cout_fat)
dQ_wp = Flow_wp * (C_art - Cout_wp) - dQmet_wp
dQ_pp = Flow_pp * (C_art - Cout_pp)
dQ_met = dQmet_wp

Outputs

# Define the computation output times
t <- seq(from=0, to=1440, by=10)
# Solve ODE
library(deSolve)
out <- ode(times=t, func=bd.model, 
           y=y, parms=parameters)
head(out)

##      time      Q_fat       Q_wp       Q_pp      Q_met       C_ven      C_art
## [1,]    0 0.00000000 0.00000000 0.00000000 0.00000000 0.000000000 0.01606403
## [2,]   10 0.02293618 0.03724892 0.07427798 0.06645654 0.003338318 0.01819035
## [3,]   20 0.04722954 0.04026245 0.14189019 0.16431358 0.004379098 0.01885327
## [4,]   30 0.07221315 0.04152080 0.20176415 0.26661643 0.005210310 0.01938270
## [5,]   40 0.09777386 0.04256838 0.25471138 0.37175647 0.005941936 0.01984870
## [6,]   50 0.12383371 0.04349410 0.30153555 0.47935418 0.006589697 0.02026129

plot(out)

Hands on Exercise

Task 5: Apply the PBPK model

• Based on the OSHA permissible exposure limit (PEL), the time-weighted average is 1 ppm. If a worker is under the exposure for a long time, what is the estimated blood concentration.