Workflow of pharmacokinetic modeling (model construction + data manipulation + visualization)
Nan-Hung Hsieh
2019/12/09 (update: 2019-12-06)
0 Prerequisites
Import the R packages
library(simuloR)
library(httk)
library(deSolve)
library(tidyverse)
theme_set(theme_bw())Task 1: Exploratory data analysis
- 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
Theophdata 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)## Grouped Data: conc ~ Time | Subject
## 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.58tail(Theoph)## Grouped Data: conc ~ Time | Subject
## Subject Wt Dose Time conc
## 127 12 60.5 5.3 3.52 9.75
## 128 12 60.5 5.3 5.07 8.57
## 129 12 60.5 5.3 7.07 6.59
## 130 12 60.5 5.3 9.03 6.11
## 131 12 60.5 5.3 12.05 4.57
## 132 12 60.5 5.3 24.15 1.17subset(Theoph, Subject == 1)## Grouped Data: conc ~ Time | Subject
## 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
## 7 1 79.6 4.02 5.10 8.36
## 8 1 79.6 4.02 7.03 7.47
## 9 1 79.6 4.02 9.05 6.89
## 10 1 79.6 4.02 12.12 5.94
## 11 1 79.6 4.02 24.37 3.28subset(Theoph, Subject == 1) %>% summary()## Subject Wt Dose Time conc
## 1:11 Min. :79.6 Min. :4.02 Min. : 0.000 Min. : 0.740
## 1st Qu.:79.6 1st Qu.:4.02 1st Qu.: 0.845 1st Qu.: 4.610
## Median :79.6 Median :4.02 Median : 3.820 Median : 6.890
## Mean :79.6 Mean :4.02 Mean : 5.950 Mean : 6.439
## 3rd Qu.:79.6 3rd Qu.:4.02 3rd Qu.: 8.040 3rd Qu.: 8.470
## Max. :79.6 Max. :4.02 Max. :24.370 Max. :10.500s1 <- subset(Theoph, Subject == 1)
plot(s1$Time, s1$conc)
plot(Theoph$Time, Theoph$conc, col=Theoph$Subject)
par(mfrow=c(3,4))
for(i in 1:12){
dat <- subset(Theoph, Subject == i)
plot(dat$Time, dat$conc, main = i)
}
Theoph %>% ggplot(aes(x = Time, y =conc, color = Subject)) +
geom_point()
Theoph %>% ggplot(aes(x = Time, y =conc, color = Subject)) +
geom_point() + geom_line() +
labs(x = "Time (hr)", y = "Concentration (mg/L)")
Theoph %>% ggplot(aes(x = Time, y =conc)) +
geom_point() + geom_line() +
facet_wrap(~Subject) +
labs(x = "Time (hr)", y = "Concentration (mg/L)")
class(Theoph$Subject)## [1] "ordered" "factor"Theoph$Subject## [1] 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 3
## [24] 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 5 5
## [47] 5 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 7 7 7
## [70] 7 7 7 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 9 9 9 9
## [93] 9 9 9 9 9 9 9 10 10 10 10 10 10 10 10 10 10 10 11 11 11 11 11
## [116] 11 11 11 11 11 11 12 12 12 12 12 12 12 12 12 12 12
## Levels: 6 < 7 < 8 < 11 < 3 < 2 < 4 < 9 < 12 < 10 < 1 < 5Theoph$Subject <- factor(Theoph$Subject,
level = c("1", "2", "3", "4", "5", "6",
"7", "8", "9", "10", "11","12"))Theoph$Subject## [1] 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 3
## [24] 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 5 5
## [47] 5 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6 6 6 7 7 7
## [70] 7 7 7 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 9 9 9 9
## [93] 9 9 9 9 9 9 9 10 10 10 10 10 10 10 10 10 10 10 11 11 11 11 11
## [116] 11 11 11 11 11 11 12 12 12 12 12 12 12 12 12 12 12
## Levels: 1 < 2 < 3 < 4 < 5 < 6 < 7 < 8 < 9 < 10 < 11 < 12Theoph %>% ggplot(aes(x = Time, y =conc)) +
geom_point() + geom_line() +
facet_wrap(~Subject) +
labs(x = "Time (hr)", y = "Concentration (mg/L)")
s1.Cmax <- max(s1$conc)
s1.Cmax## [1] 10.5which(s1$conc == s1.Cmax)## [1] 4t.pt <- which(s1$conc == s1.Cmax)
T_max <- s1$Time[t.pt]Cmax_pop <- Theoph %>% group_by(Subject) %>% summarise(C_max = max(conc))
Cmax_pop## # A tibble: 12 x 2
## Subject C_max
## <ord> <dbl>
## 1 1 10.5
## 2 2 8.33
## 3 3 8.2
## 4 4 8.6
## 5 5 11.4
## 6 6 6.44
## 7 7 7.09
## 8 8 7.56
## 9 9 9.03
## 10 10 10.2
## 11 11 8
## 12 12 9.75hist(Cmax_pop$C_max)
boxplot(Cmax_pop$C_max)
Theoph %>% group_by(Subject) %>% summarise(C_max = max(conc), Wt=mean(Wt))## # A tibble: 12 x 3
## Subject C_max Wt
## <ord> <dbl> <dbl>
## 1 1 10.5 79.6
## 2 2 8.33 72.4
## 3 3 8.2 70.5
## 4 4 8.6 72.7
## 5 5 11.4 54.6
## 6 6 6.44 80
## 7 7 7.09 64.6
## 8 8 7.56 70.5
## 9 9 9.03 86.4
## 10 10 10.2 58.2
## 11 11 8 65
## 12 12 9.75 60.5Theoph %>% group_by(Subject) %>%
summarise(C_max = max(conc), Wt=mean(Wt)) %>%
ggplot(aes(x=Wt, y=C_max)) + geom_point()
Theoph %>% group_by(Subject) %>%
summarise(C_max = max(conc), Wt=mean(Wt)) %>%
ggplot(aes(x=Wt, y=C_max)) + geom_point() + geom_smooth(method = "lm")
df <- Theoph %>% group_by(Subject) %>% summarize(C_max = max(conc))
df$T_max <- NA
df## # A tibble: 12 x 3
## Subject C_max T_max
## <ord> <dbl> <lgl>
## 1 1 10.5 NA
## 2 2 8.33 NA
## 3 3 8.2 NA
## 4 4 8.6 NA
## 5 5 11.4 NA
## 6 6 6.44 NA
## 7 7 7.09 NA
## 8 8 7.56 NA
## 9 9 9.03 NA
## 10 10 10.2 NA
## 11 11 8 NA
## 12 12 9.75 NAdat <- subset(Theoph, Subject==1)
C_max.pt <- which(dat$conc == df$C_max[1])
df$T_max[1] <- dat$Time[C_max.pt]
df## # A tibble: 12 x 3
## Subject C_max T_max
## <ord> <dbl> <dbl>
## 1 1 10.5 1.12
## 2 2 8.33 NA
## 3 3 8.2 NA
## 4 4 8.6 NA
## 5 5 11.4 NA
## 6 6 6.44 NA
## 7 7 7.09 NA
## 8 8 7.56 NA
## 9 9 9.03 NA
## 10 10 10.2 NA
## 11 11 8 NA
## 12 12 9.75 NAfor(i in 2:12){
dat <- subset(Theoph, Subject==i)
C_max.pt <- which(dat$conc == df$C_max[i])
df$T_max[i] <- dat$Time[C_max.pt]
}
df## # A tibble: 12 x 3
## Subject C_max T_max
## <ord> <dbl> <dbl>
## 1 1 10.5 1.12
## 2 2 8.33 1.92
## 3 3 8.2 1.02
## 4 4 8.6 1.07
## 5 5 11.4 1
## 6 6 6.44 1.15
## 7 7 7.09 3.48
## 8 8 7.56 2.02
## 9 9 9.03 0.63
## 10 10 10.2 3.55
## 11 11 8 0.98
## 12 12 9.75 3.52Task 2: Model development
- Develop the non compartment model and compartment model in R and MCSim
The MCSim Syntax of the model description file
# Model description file (this is a comment)
<Global variable specifications>
States = {
<state variables for the model, such as quantity>
}
Outputs = {
<output variables, such as concentration>
}
Inputs = {
<input variables, such as exposure dose>
}
Initialize {
<Equations for initializing or scaling model parameters>
}
Dynamics {
<Equations for computing derivatives of the state variables>
}
CalcOutputs {
<Equations for computing output variables>
}
End. # mandatory ending keywordGeneral syntax
Variable assignments
<variable-name> = <constant-value-or-expression> ;Colon conditional assignments
<variable-name> = (<test> ? <value-if-true> : <value-if-false>);For example
Adjusted_param = (Input_var > 0.0 ? Param * 1.1 : Param);Syntax of (simulation) input-file
For the basic simulation
# Input-file (text after # are comments)
<Global assignments and specifications>
Simulation {
<Specifications for first simulation>
}
Simulation {
<Specifications for second simulation>
}
# Unlimited number of simulation specifications
End. # Mandatory End keyword. Everything after this line is ignored1.1 Define the state variables States{}
States are variables for which a first-order differential equation is defined in the Dynamics{} section. Remember, when naming the variables, be sure to let all variable names begin with a capital letter followed by meaningful lower case subscripts. In States{} section, each variable needs to separate by comma ,. In addition, remember to provide the unit that is used for the state variables.
1.2 Define the outputs variables Outputs{}
Outputs are dependent model variables (obtainable at any time as analytical functions of the states, inputs or parameters) that do not have dynamics. They must receive assignments in either the Dynamics or CalcOutputs sections. When constructing the model, we need to have an output variable for checking mass balance. Then, assign the variable that we are interested in or have an experiment to make the comparison, such as the blood concentration.
1.3 Define the input variable(s) Inputs{}
Inputs are variables independent of the other variables, but eventually varying with time (for example an exposure concentration to a chemical). GNU MCSim has various functions to define the different exposure types as:
PerDose(): # specifies a periodic input of constant
PerDose(<magnitude>, <period>, <initial-time>, <exposure-time>);PerExp(): # specifies a periodic exponential input.
PerExp(<magnitude>, <period>, <initial-time>, <decay-constant>);PerTransit(): #models a delayed input mechanism (e.g., gut absorption)
PerTransit(<magnitude>, <period>, <initial-time-in-period>, <decay-constant>, <number-of-input-compartments>);NDoses(): # specifies a number of stepwise inputs of variable magnitude and their starting times
NDoses(<n>, <list-of-magnitudes>, <list-of-initial-times>);Spikes(): # specifies a number of instantaneous inputs of variable magnitude and their exact times of occurrence.
Spikes(<n>, <list-of-magnitudes>, <list-of-times>);
1.4 Parameters of the model
Here we need to define the “default” value of input parameters. Unlike previous variables that need to use curly brackets to group these variables, the input parameters can be put in any places without the restriction. In addition, don’t forget to put a semicolon (;) after the definition of each parameter. If the input parameter doesn’t have a default value, just put a semicolon after the name of the parameter or used an arbitrary number.
The PBPK parameters always include the parameters that can describe absorption, distribution (e.g., partition coefficient), metabolism (e.g., Michaelis-Menten constant), and elimination (e.g., rate Constant). These parameters can be classified to chemical-specific (e.g., molecular weight) and anatomy (e.g., body weight, blood flow, tissue weight/volume) parameter. In addition, the unit conversion factor can define in this part.
The parameters of the PBPK model might include:
Unit conversion factor
Exposure modeling parameters
Chemical-specific parameters
Physiological parameters
Pharmacokinetic parameters
Scale parameters (no default value)
1.5 Define the parameter initialization and scaling Initialize{}
This section is used to scale specific model parameters and resolve dependencies between parameters. Generally, the scaling involves a change of units or conversion from percentage to actual units.
1.6 Define the dynamics of the simulation Dynamics{}
The equations given in this section will be called by the integration routines at each integration step. It includes specification of differential equations.
The derivative of a state variable is defined using the dt() operator, as shown here:
dt(state-variable) '=' constant-value-or-expression ';'
1.7 Define the output calculations CalcOutputs{}
In this section, the equations are given in this section will be called by the simulation program at each output time specified by a Print() or PrintStep() statement. In this way, output computations are done efficiently, only when values are to be saved. Here is the example to assign the variables in CalcOutputs{} section. Again, don’t forget the semicolon and put the keyword End. at the end of model.
Same as model file, we highly recommend to use comments (#) to annotate your code in the beginning. This comment might include:
- Brief description of this input file
- What model file that is used to conduct the simulation
- Others…
1.8 Setting integrator Integrate()
GNU MCSim provides three types of integrator that include Lsodes, Cvodes, and Euler. Here, we used Lsodes (Livermore Solver for Ordinary Differential Equation) as a primary solver, which was originated from the SLAC Fortran library. The Lsodes function uses adaptive numerical methods to advance a solution to a system of ordinary differential equations one time-step, given values for the variables Y and X. It solves the initial value problem for stiff or non-stiff systems of ordinary differential equations (ODE).
The syntax for Lsodes is: Integrate(Lsodes, <rtol>, <atol>, <method>);
where <rtol> is a scalar specifying the relative error tolerance for each integration step. The scalar <atol> specifies the absolute error tolerance parameter. Those tolerances are used for all state variables. The estimated local error for a state variable y is controlled so as to be roughly less (in magnitude) than rtol*|y| + atol. Thus the local error test passes if, for each state variable, either the absolute error is less than <rtol>. Set <rtol> to zero for pure absolute error control, and set <atol> to zero for pure relative error control. Caution: actual (global) errors may exceed these local tolerances, so choose them conservatively. Decreasing the tolerance leads to more accurate results but at the cost of significant increase in time taken.
The <method> flag should be 0 (zero) for non-stiff differential systems and 1 or 2 for stiff systems. You should try flag 0 or 1 and select the fastest for equal accuracy of output unless insight from your system leads you to choose one of them a priori. If you specify Lsodes with parameters 1e-5, 1e-7 and 1.
1.9 Setting paramerer(s)
For the simple simulation, we can assign the given parameter values in the input file to replace the default values in the model file. The parameter can be set in the simulation{} section as a local setting or outside the section as a global setting. Same as the parameter setting in the model file, use = to define the given value of parameters and put ; after the definition.
Input functions
These functions can use to different exposure types
- PerDose(): # specifies a periodic input of constant
PerDose(<magnitude>, <period>, <initial-time>, <exposure-time>);
- PerExp(): # specifies a periodic exponential input.
PerExp(<magnitude>, <period>, <initial-time>, <decay-constant>);
- PerTransit(): models a delayed input mechanism
PerTransit(<magnitude>, <period>, <initial-time-in-period>,
<decay-constant>, <number-of-input-compartments>);
- NDoses(): specifies a number of stepwise inputs of variable magnitude and their starting times
NDoses(<n>, <list-of-magnitudes>, <list-of-initial-times>);
- Spikes(): specifies a number of instantaneous inputs of variable magnitude and their exact times of occurrence.
Spikes(<n>, <list-of-magnitudes>, <list-of-times>);
1.10 Setting simulation simulation{}
After the global specifications, we can define the local parameters. In addition, we need to define the output through Print() or PrintStep() functions. The arguments of Print() and PrintStep() are the comma-separated list of variable names (at least one and up to MAX_PRINT_VARS, which is set to 10 by default) and followed by given time-points (limited by the available memory at run time) as:
Print(<identifier1>, <identifier2>, ..., <time1>, <time2>, ...);
PrintStep(<identifier1>, <identifier2>, ..., <start-time>, <end-time>, <time-step>);
Note: The output variables cannot receive assignments in simulation input files.
- Non compartment model
\[C(t) = \frac{F\cdot D \cdot ka}{V(k_a - k_e)}\cdot(\mathrm{e}^{-k_et}-\mathrm{e}^{-k_a{t}})\]
- Compartment model
\[\frac{dA_{gut}}{dt}=-k_aA_{gut}\]
\[\frac{dA}{dt}=k_aA_{gut}-k_eA_e\] \[C=A/V\]
- Non compartment model (R)
ffpk <- function(D, Fa, ka, ke, V, t){
MW <- 180.167
A <- (D * Fa * ka)/(V * (ka - ke))
C <- A * exp(-ke * t) - A * exp(-ka * t)
#C_uM <- C/MW*1000
return(C)
}- Non compartment model (MCSim)
Outputs = {
A,
Conc,
C_rest_uM
}
## ingested input (mg)
IngDose = 1.0;
## Molecular weight (g/mole)
MW = 100;
## Absorption fraction (-)
Fgutabs = 1.0; #
## Absorption rate constant (/h)
kgutabs = 1.0; #
## Elimination rate constants (/h)
kelim = 1.0;
## Distribution volumes (L)
Vdist = 1.0;
CalcOutputs {
A = (IngDose * Fgutabs * kgutabs)/(Vdist * (kgutabs - kelim));
Conc = A * exp(-kelim * t) - A * exp(-kgutabs * t);
Conc_uM = Conc / MW * 1000;
}
End.- Compartment model (R)
pbtk1cpt <- function(t, state, parameters) {
with(as.list(c(state, parameters)), {
dAgutlument = - kgutabs * Agutlument
dAcompartment = kgutabs * Agutlument - ke * Acompartment
dAmetabolized = ke * Acompartment
Ccompartment = Acompartment / vdist * BW;
list(c(dAgutlument, dAcompartment, dAmetabolized),
"Ccompartment" = Ccompartment)
})
}- Compartment model (MCSim)
States = {
A_rest,
A_elim
}; # Quantity in central compartment (mg)
Inputs = {
Oral_input
}; # Drug input (mg)
Outputs = {
C_rest, (mg/l)
C_rest_uM, (uM)
}; # Concentration in central compartment (uM)
# Oral input modeling
Absdose;
IngDose = 1.0; # ingested input (mg)
Fgutabs = 1.0; #
Period = 0.0; # period of the exposure/no exposure cycle (h)
Tlag = 0.0; # Absorption lagtime (h)
kgutabs = 0.1; # Intestinal absortion rate (/h)
Oral_input = PerExp (Absdose, Period, Tlag, kgutabs);
# Molecular weight (g/mole)
MW = 100;
# Distribution volumes (L)
Vdist = 1;
# Elimination rate constants (/h)
kelim = 1;
Initialize{ Absdose = IngDose * Fgutabs; }
Dynamics {
dt (A_rest) = kgutabs * Oral_input - kelim * A_rest;
dt (A_elim) = kelim * A_rest;
}
CalcOutputs {
C_rest = A_rest / Vdist;
C_rest_uM = C_rest / MW * 1000;
}
End.Task 3: Parameter setting and simulation
- Extract the chemical information from httk package
parms <- httk::parameterize_1comp(chem.name = "theophylline")## Warning in available_rblood2plasma(chem.cas = chem.cas, chem.name =
## chem.name, : Human in vivo Rblood2plasma returned.parms## $Vdist
## [1] 0.7801997
##
## $kelim
## [1] 0.2653932
##
## $Clint
## [1] 0.7476636
##
## $Clint.dist
## [1] NA
##
## $Funbound.plasma
## [1] 0.8965508
##
## $Funbound.plasma.dist
## [1] NA
##
## $Funbound.plasma.adjustment
## [1] 0.9961675
##
## $Fhep.assay.correction
## [1] 0.9564623
##
## $Pow
## [1] 0.8078951
##
## $pKa_Donor
## pKa_Donor
## "7.82"
##
## $pKa_Accept
## pKa_Accept
## NA
##
## $MA
## [1] NA
##
## $kgutabs
## [1] 2.18
##
## $Rblood2plasma
## [1] 1.33
##
## $million.cells.per.gliver
## [1] 110
##
## $liver.density
## [1] 1.05
##
## $hematocrit
## [1] 0.44
##
## $MW
## [1] 180.167
##
## $Fgutabs
## [1] 0.98
##
## $hepatic.bioavailability
## [1] 0.9328532
##
## $BW
## [1] 70- Use the parameters in the developed model and conduct the simulation
- Compare the difference between data and the model simulation result (Cmax, Tmax)
R
MW <- parms$MW
Fgutabs <- parms$Fgutabs * parms$hepatic.bioavailability
kgutabs <- parms$kgutabs
kelim <- parms$kelim
Vdist <- parms$Vdistdose1 <- subset(Theoph, Subject == 1) %>% summarise(dose = mean(Dose))
t1 <- subset(Theoph, Subject == 1) %>% select(Time)
out1 <- ffpk(D = dose1$dose, Fa = Fgutabs, ka = kgutabs, ke = kelim, V = Vdist, t=t1$Time)
out1## [1] 0.000000000 1.909156333 3.062364221 3.517523615 3.072105720
## [6] 1.944726071 1.385452594 0.830168992 0.485673409 0.215030371
## [11] 0.008328742conc1 <- subset(Theoph, Subject == 1) %>% select(conc)
plot(t1$Time, conc1$conc)
lines(t1$Time, out1)
data.frame(t1, out1)## Time out1
## 1 0.00 0.000000000
## 2 0.25 1.909156333
## 3 0.57 3.062364221
## 4 1.12 3.517523615
## 5 2.02 3.072105720
## 6 3.82 1.944726071
## 7 5.10 1.385452594
## 8 7.03 0.830168992
## 9 9.05 0.485673409
## 10 12.12 0.215030371
## 11 24.37 0.008328742MCSim
Integrate (Lsodes, 1e-8, 1e-8, 1);
MW = 180.167; # moleculor mass (g/mol)
IngDose = 4.02; # ingested dose (mg)
Period = 48; # period of the exposure/no exposure cycle (h)
Fgutabs = 0.9142; # bioavailability (-)
kgutabs = 2.18; # absortion rate (/h)
kelim = 0.2654; # elimination rate constant (/h)
Vdist = 0.7802; # volumes (L)
Simulation {
Print (C_rest 0.00 0.25 0.57 1.12 2.02 3.82 5.10 7.03 9.05 12.12 24.37);
}
End.mcsim(model = "models/pbtk1cpt.model", input = "inputs/pbtk1cpt.in")## * Creating executable program, pleas wait...## * Created executable program 'models/mcsim.pbtk1cpt.model'.## Executing...## Time Conc
## 1 0.25 1.90916000
## 2 0.57 3.06237000
## 3 1.12 3.51752000
## 4 2.02 3.07209000
## 5 3.82 1.94469000
## 6 5.10 1.38541000
## 7 7.03 0.83013500
## 8 9.05 0.48564700
## 9 12.12 0.21501400
## 10 24.37 0.00832746read.delim("sim.out")## Results.of.Simulation.1
## Time Conc
## 0.25 1.90916
## 0.57 3.06237
## 1.12 3.51752
## 2.02 3.07209
## 3.82 1.94469
## 5.1 1.38541
## 7.03 0.830135
## 9.05 0.485647
## 12.12 0.215014
## 24.37 0.00832746read.delim("sim.out", skip = 1)## Time Conc
## 1 0.25 1.90916000
## 2 0.57 3.06237000
## 3 1.12 3.51752000
## 4 2.02 3.07209000
## 5 3.82 1.94469000
## 6 5.10 1.38541000
## 7 7.03 0.83013500
## 8 9.05 0.48564700
## 9 12.12 0.21501400
## 10 24.37 0.00832746out <- mcsim(model = "models/pbtk1cpt.model", input = "inputs/pbtk1cpt.in")## * Creating executable program, pleas wait...## * Created executable program 'models/mcsim.pbtk1cpt.model'.## Executing...out## Time Conc
## 1 0.25 1.90916000
## 2 0.57 3.06237000
## 3 1.12 3.51752000
## 4 2.02 3.07209000
## 5 3.82 1.94469000
## 6 5.10 1.38541000
## 7 7.03 0.83013500
## 8 9.05 0.48564700
## 9 12.12 0.21501400
## 10 24.37 0.00832746plot(t1$Time, conc1$conc)
lines(out$Time, out$Conc)
Task 4: Developing the PBPK model
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 metabolismy <- 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.modelConstants
# 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)
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
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_wpConcentrations - 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_totConcentrations - 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# Define the model equations
bd.model = function(t, y, parameters) {
with (as.list(y), {
with (as.list(parameters), {
# Define some 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
# 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
# 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
# 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)
# 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
# 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
}
# 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
# 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# Define the computation output times
t <- seq(from=0, to=1440, by=10)
# Solve ODE
out <- ode(times=t, func=bd.model, y=y, parms=parameters)
head(out)## time Q_fat Q_wp Q_pp Q_met C_ven
## [1,] 0 0.00000000 0.00000000 0.00000000 0.00000000 0.000000000
## [2,] 10 0.02293618 0.03724892 0.07427798 0.06645654 0.003338318
## [3,] 20 0.04722954 0.04026245 0.14189019 0.16431358 0.004379098
## [4,] 30 0.07221315 0.04152080 0.20176415 0.26661643 0.005210310
## [5,] 40 0.09777386 0.04256838 0.25471138 0.37175647 0.005941936
## [6,] 50 0.12383371 0.04349410 0.30153555 0.47935418 0.006589697
## C_art
## [1,] 0.01606403
## [2,] 0.01819035
## [3,] 0.01885327
## [4,] 0.01938270
## [5,] 0.01984870
## [6,] 0.02026129plot(out)
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.
t_exposure <- 43200 # 30 days
C_inh <- approxfun(x = c(0, t_exposure), y=c(10,0), method = "constant", f = 0, rule = 2) plot(C_inh(1:t_exposure), type="l", xlab = "Time (min)", ylab = "Butadiene inhaled concentration (ppm)")
t <- seq(from=0, to=t_exposure, by=100)
out <- ode(times=t, func=bd.model, y=y, parms=parameters)tail(out, 10)## time Q_fat Q_wp Q_pp Q_met C_ven C_art
## [424,] 42300 11.24631 0.05555652 0.7245784 582.4532 0.01517190 0.025727662
## [425,] 42400 11.24632 0.05555653 0.7245785 583.8421 0.01517190 0.025727663
## [426,] 42500 11.24632 0.05555653 0.7245785 585.2310 0.01517190 0.025727664
## [427,] 42600 11.24632 0.05555653 0.7245785 586.6199 0.01517190 0.025727664
## [428,] 42700 11.24633 0.05555653 0.7245785 588.0088 0.01517190 0.025727665
## [429,] 42800 11.24633 0.05555653 0.7245785 589.3977 0.01517190 0.025727665
## [430,] 42900 11.24633 0.05555653 0.7245785 590.7866 0.01517191 0.025727666
## [431,] 43000 11.24633 0.05555653 0.7245785 592.1756 0.01517191 0.025727666
## [432,] 43100 11.24634 0.05555653 0.7245786 593.5645 0.01517191 0.025727667
## [433,] 43200 11.24634 0.05555619 0.7245784 594.9534 0.01517188 0.009663619plot(out)
Session info
sessionInfo()## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: elementary OS 5.0 Juno
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] forcats_0.4.0 stringr_1.4.0 dplyr_0.8.3 purrr_0.3.3
## [5] readr_1.3.1 tidyr_1.0.0 tibble_2.1.3 ggplot2_3.2.1
## [9] tidyverse_1.3.0 deSolve_1.25 httk_1.10.1 simuloR_0.0.1
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.3 lubridate_1.7.4 msm_1.6.7
## [4] mvtnorm_1.0-11 lattice_0.20-38 utf8_1.1.4
## [7] assertthat_0.2.1 zeallot_0.1.0 digest_0.6.23
## [10] R6_2.4.1 cellranger_1.1.0 backports_1.1.5
## [13] reprex_0.3.0 survey_3.36 evaluate_0.14
## [16] httr_1.4.1 pillar_1.4.2 rlang_0.4.2
## [19] lazyeval_0.2.2 readxl_1.3.1 rstudioapi_0.10
## [22] data.table_1.12.6 Matrix_1.2-17 rmarkdown_1.15
## [25] labeling_0.3 splines_3.6.1 munsell_0.5.0
## [28] broom_0.5.2 compiler_3.6.1 modelr_0.1.5
## [31] xfun_0.9 pkgconfig_2.0.3 htmltools_0.3.6
## [34] mitools_2.4 tidyselect_0.2.5 expm_0.999-4
## [37] fansi_0.4.0 viridisLite_0.3.0 crayon_1.3.4
## [40] dbplyr_1.4.2 withr_2.1.2 grid_3.6.1
## [43] nlme_3.1-142 jsonlite_1.6 gtable_0.3.0
## [46] lifecycle_0.1.0 DBI_1.0.0 magrittr_1.5
## [49] scales_1.1.0 cli_1.1.0 stringi_1.4.3
## [52] farver_2.0.1 fs_1.3.1 xml2_1.2.2
## [55] vctrs_0.2.0 generics_0.0.2 tools_3.6.1
## [58] glue_1.3.1 hms_0.5.2 survival_2.44-1.1
## [61] yaml_2.2.0 colorspace_1.4-1 rvest_0.3.5
## [64] knitr_1.25 haven_2.2.0Reference
Bois F.Y. 2013. Bayesian Inference. In: Reisfeld B., Mayeno A. (eds) Computational Toxicology. Methods in Molecular Biology (Methods and Protocols), vol 930. Humana Press, Totowa, NJ
Bois F.Y., Brochot C. (2016) Modeling Pharmacokinetics. In: Benfenati E. (eds) In Silico Methods for Predicting Drug Toxicity. Methods in Molecular Biology, vol 1425. Humana Press, New York, NY
Pearce, R.G., Setzer, R.W., Strope, C.L., Wambaugh, J.F. and Sipes, N.S., 2017. httk: R package for high-throughput toxicokinetics. Journal of statistical software, 79(4), p.1.
Soetaert, K.E., Petzoldt, T. and Setzer, R.W., 2010. Solving differential equations in R: package deSolve. Journal of statistical software, 33.