Tutorial 1 - Walk-through of Working Models

0 Overview

In this tutorial we want to know,

How to write Physiologically based pharmacokinetic (PBPK) model?

How to write input file and conduct simulation?

How to analyze GNU MCSim output with R?

Here, we’ll use tetrachloroethylene (PERC) and ethylbenzene (EB) as examples to illustrate how to perform straight simulation in GNU MCSim and analyze the output result in R.

1 Create PBPK model

This section illustrates how to build a PBPK model in GNU MCSim. In the beginning, use comments (#) to annotate your code. Here is the example that you can use to comment on your model file.

Brief description of this model
What is the version number
When you create this model
Summarize the unit used in this model

Then, follow the next seven steps to build the PBPK model.

1.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.

2 Create input/simulation file

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…

2.1 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 , or the relative 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 2 you should provide the Jacobian of your differential system. The good starting point for and is about 1e-6. The default integration method in Lsodes with parameters 1e-5, 1e-7 and 1.

2.2 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.

2.3 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.

After finish the model- and input-file, we can use mcsim() to run the simulation. This is the example of use mcsim() to solve Tetrachloroethylene (PERC) PBPK model

mcsim(model = "perc.model.R", input = "perc.in.R", dir = "modeling/perc")

## * Created executable program 'mcsim.perc.model.R.exe'.

## Execute: ./mcsim.perc.model.R.exe modeling/perc/perc.in.R

##       Time   C_exh_ug     C_ven
## 1    239.9 282.580000 1.9999800
## 2    245.0  97.737100 1.8346600
## 3    270.0  78.504900 1.4736500
## 4    360.0  45.480900 0.8537420
## 5   1320.0   4.268720 0.0801300
## 6   2760.0   2.913030 0.0546818
## 7   4260.0   2.261710 0.0424556
## 8   5700.0   1.774920 0.0333177
## 9   8580.0   1.093980 0.0205355
## 10 10020.0   0.859231 0.0161290

2.4 Working with the R package deSolve

In addition to use MCSim to solve ODE, there is an alternative way to solve ODE. The integrators provided by deSolve are improved implementations of the lsode family of integrators used by GNU MCSim. They provide a few more options than GNU MCSim. such as the functions that solve initial value problems of a system of partial differential equations (‘PDE’), of differential algebraic equations (‘DAE’), and of delay differential equations. They provide an interface to the FORTRAN functions ‘lsoda’, ‘lsodar’, ‘lsode’, ‘lsodes’ of the ‘ODEPACK’ collection, to the FORTRAN functions ‘dvode’, ‘zvode’ and ‘daspk’ and a C-implementation of solvers of the ‘Runge-Kutta’ family with fixed or variable time steps. However, GNU MCSim is the fastest option in intensive simulation.

3 Analyze GNU MCSim output with R

GNU MCSim doesn’t include a comprehensive tool to analyze and visualize the simulation result. The best way to do further data analysis and visualization is to use R. This instruction will provide some basic R functions that can use to learn and understand how to analyze GNU MCSim output with R.

3.1 Loading output

The simulation output is written under tab-delimited text file with default name sim.out. It is located in the directory. You can check your file through files tab on the bottom-right panel in RStudio. Or use file.exists() to check the file as,

file.exists("sim.out")

## [1] TRUE

If you can find your output file in your working directory, then you can us R function read.delim() to read it. Since the first line of the output file is used to describe the number of simulation, it can’t be used in the formal data structure. Therefore, we need to assign skip = 1 to ignore this line. Then run the following code and check your result:

read.delim(file = "sim.out", skip = 1)

##       Time   C_exh_ug     C_ven
## 1    239.9 282.580000 1.9999800
## 2    245.0  97.737100 1.8346600
## 3    270.0  78.504900 1.4736500
## 4    360.0  45.480900 0.8537420
## 5   1320.0   4.268720 0.0801300
## 6   2760.0   2.913030 0.0546818
## 7   4260.0   2.261710 0.0424556
## 8   5700.0   1.774920 0.0333177
## 9   8580.0   1.093980 0.0205355
## 10 10020.0   0.859231 0.0161290

Through this function, we can only read the content in the output file, but we can’t use the data unless we assign it to a variable. Therefore, we need to assign a variable for the output to help us analyze it. The assigned variables will be temporarily stored in memory.

3.2 Assign variables

Use the assignment operator <- to create new variables.

out <- read.delim(file = "sim.out", skip = 1)

Once a variable is created, we can use the variable name to refer to the simulation data it was assigned. To view the variable in RStudio, you can interact with it in the top-right panel. Or, use View() to check it:

View(out)

3.3 Manipulating output

Now, our data is loaded into R. We can start doing further analysis with them. First, let’s check what type it is:

class(out)

## [1] "data.frame"

The output tells us that it’s a “data frame”. This type of data structure constructed by the specific number of column (variable) and row (observation). We can see the dimensions of the data frame with the function dim():

dim(out)

## [1] 10  3

Or use nrow() and ncol() to check the information:

nrow(out)

## [1] 10

ncol(out)

## [1] 3

Also, the names() can be used to get and set the names of variables

names(out)

## [1] "Time"     "C_exh_ug" "C_ven"

Some basic functions can be used to summarize the result, such as mean(), min(), max(). Here is an example to find the maximum concentration in blood:

max(out$C_ven)

## [1] 1.99998

Here, the $ is used to assign the specific variable that was stored in the object. Also, you can quickly check all maximum or minimum concentration through apply:

apply(out, MARGIN = 2, FUN = max)

##        Time    C_exh_ug       C_ven 
## 10020.00000   282.58000     1.99998

You can also use which() to find the value as:

which(out$C_ven == max(out$C_ven))

## [1] 1

which(out$C_ven == min(out$C_ven))

## [1] 10

But, the which() only give us the location, not time point. Therefore, we need to use this location to find the corresponding time point by:

which_row <- which(out$C_ven == max(out$C_ven))
out$Time[which_row]

## [1] 239.9

The following code can let you set the variable name:

names(out) <- c("Time", "C_air", "C_ven")

Check again,

names(out)

## [1] "Time"  "C_air" "C_ven"

If you only want to change the specific variable, you can use:

names(out)[3] <- "C_blood"
names(out)

## [1] "Time"    "C_air"   "C_blood"

The unit conversion is the cruical part in data manipulation as well. Here we want to transfer the unit from mg/l to mol/l,

MW <- 165.822 # g/mol
out$C_blood_mol <- out$C_blood / MW * 1000 # mg/l -> mol/l
out

##       Time      C_air   C_blood C_blood_mol
## 1    239.9 282.580000 1.9999800 12.06100517
## 2    245.0  97.737100 1.8346600 11.06403252
## 3    270.0  78.504900 1.4736500  8.88693901
## 4    360.0  45.480900 0.8537420  5.14854483
## 5   1320.0   4.268720 0.0801300  0.48322900
## 6   2760.0   2.913030 0.0546818  0.32976203
## 7   4260.0   2.261710 0.0424556  0.25603117
## 8   5700.0   1.774920 0.0333177  0.20092449
## 9   8580.0   1.093980 0.0205355  0.12384062
## 10 10020.0   0.859231 0.0161290  0.09726695

Finally, if you prefer to use Excel, you can use write.csv() to transfer your object to Comma-separated values (csv) file, then you can use Excel to open it and do additional analysis.

write.csv(out, file = "sim.csv", row.names = F)

3.4 Plotting

Visualizing the output data is an important work to communicate the simulation result. After the data manipulation, we can use plot() to create the PK diagram as:

plot(x = out$Time, y = out$C_blood)

Usually, the basic assignment cannot provide high quality and comprehensive figure that can further use in our report or published paper. Therefore, we need to use additional assignments such as xlab, ylab, main to add the description or modification.

plot(x = out$Time, y= out$C_blood, 
     type = "l", col = "red", xlab = "Time post exposure (min)", ylab = "Blood concentration (mg/L)",
     log = "x", main = "PERC PBPK modeling")
legend("topright", legend = "Prediction", col="red", lty=1)

Here is an example to plot the two simulation results (concentration of exhaled air and blood) in the same plot. However, we can find that the range of y-axis doesn’t cover all simulation points.

plot(x = out$Time, y= out$C_air, type = "l", col = "red", 
     xlab = "Time post exposure (min)", ylab = "Con. (mg/L in blood & ug/l in exhaled air)", log = "xy",
     main = "PERC PBPK modeling")
lines(x = out$Time, y= out$C_blood, col = "blue")
legend("topright", legend = c("Conc. air", "Conc. blood"), col=c("red","blue"), lty=1)

To solve this problem, we can customize the range x-axis with the range() to define the range of x-axis as:

y_rng <- range(out$C_air, out$C_blood)
plot(x = out$Time, y= out$C_air, type = "l", col = "red",
     xlab = "Time post exposure (min)", ylab = "Con. (mg/L in blood & ug/l in exhaled air)",
     log = "xy", main = "PERC PBPK modeling", ylim = y_rng)
lines(x = out$Time, y= out$C_blood, col = "blue")
legend("topright", legend = c("Conc. air", "Conc. blood"), col=c("red","blue"), lty=1)

In addition, the relationship of concentration in blood and air can further be visualized through this way,

plot(x = out$C_blood, y= out$C_air, log = "xy", 
     type = "n", xlab = "Blood (mg/l)" , ylab = "Exhaled air (ug/l)", 
     main = "Phase plane (Unit of time: hr)")
text(out$C_blood, y= out$C_air, label = round(out$Time/60, 1), col="blue")
lines(x = out$C_blood, y= out$C_air, lty = "dotted")

Finally, you can use pdf() to export the plotting result to portable document format (pdf) as:

pdf("perc_pk.pdf")
y_rng <- range(out$C_air, out$C_blood)
plot(x = out$Time, y= out$C_air, type = "l", col = "red",
     xlab = "Time post exposure (min)", ylab = "Con. (mg/L in blood & ug/l in exhaled air)",
     log = "xy", main = "PERC PBPK modeling", ylim = y_rng)
lines(x = out$Time, y= out$C_blood, col = "blue")
legend("topright", legend = c("Conc. air", "Conc. blood"), col=c("red","blue"), lty=1)
dev.off()

4 One-compartment model

4.1 Modeling and simulation

This is an example of one-compartment model. Here we need to use model and input files of one.model.R and one.in.R to coduct the simulation. First, let’s check the output

out.1cpt <- mcsim(model = "one.model.R", input = "one.in.R")

## * Created executable program 'mcsim.one.model.R.exe'.

## Execute: ./mcsim.one.model.R.exe modeling/one.in.R

out.1cpt

##                       Time  Oral_input A_central  A_elim A_total C_central
## 1                        0        6000         0       0       0         0
## 2                        1     1635.19   4109.99 254.823 4364.81   68.4998
## 3                        2     445.641   4838.97 715.387 5554.36   80.6495
## 4                        3     121.451   4683.75  1194.8 5878.55   78.0624
## 5                        4     33.0994   4321.22 1645.68  5966.9   72.0204
## 6                        5     9.02064   3932.68  2058.3 5990.98   65.5446
## 7                        6     2.45841   3564.61 2432.93 5997.54   59.4102
## 8                        7    0.669995   3227.08 2772.25 5999.33   53.7846
## 9                        8    0.182595   2920.44 3079.38 5999.82    48.674
## 10                       9   0.0497629   2642.65  3357.3 5999.95   44.0441
## 11                      10    0.013562    2391.2 3608.78 5999.99   39.8534
## 12                      11  0.00369607   2163.66 3836.34    6000    36.061
## 13                      12   0.0010073   1957.76 4042.24    6000   32.6294
## 14                      13  0.00027452   1771.46 4228.54    6000   29.5243
## 15                      14 7.48155e-05   1602.88 4397.12    6000   26.7147
## 16                      15 2.03896e-05   1450.35 4549.65    6000   24.1724
## 17                      16 5.55682e-06   1312.33 4687.67    6000   21.8721
## 18                      17 1.51441e-06   1187.44 4812.56    6000   19.7907
## 19                      18 4.12725e-07   1074.44 4925.56    6000   17.9074
## 20                      19 1.12481e-07   972.196  5027.8    6000   16.2033
## 21                      20 3.06545e-08   879.679 5120.32    6000   14.6613
## 22                      21 8.35434e-09   795.967 5204.03    6000   13.2661
## 23                      22 2.27682e-09   720.221 5279.78    6000   12.0037
## 24                      23 6.20506e-10   651.682 5348.32    6000   10.8614
## 25                      24 1.69108e-10   589.667 5410.33    6000   9.82778
## 26 Results of Simulation 2                                                
## 27                    Time  Oral_input A_central  A_elim A_total C_central
## 28                       0       12000         0       0       0         0
## 29                       1     3270.38   8219.97 509.645 8729.62    102.75
## 30                       2     891.283   9677.94 1430.77 11108.7   120.974
## 31                       3     242.903   9367.49 2389.61 11757.1   117.094
## 32                       4     66.1988   8642.45 3291.36 11933.8   108.031
## 33                       5     18.0413   7865.35  4116.6   11982   98.3169
## 34                       6     4.91682   7129.22 4865.86 11995.1   89.1153
## 35                       7     1.33999   6454.16  5544.5 11998.7    80.677
## 36                       8     0.36519   5840.88 6158.75 11999.6    73.011
## 37                       9   0.0995258    5285.3  6714.6 11999.9   66.0662
## 38                      10    0.027124    4782.4 7217.57   12000     59.78
## 39                      11  0.00739214   4327.32 7672.68   12000   54.0915
## 40                      12  0.00201459   3915.52 8084.48   12000    48.944
## 41                      13 0.000549041   3542.91 8457.09   12000   44.2864
## 42                      14 0.000149631   3205.76 8794.24   12000    40.072
## 43                      15 4.07792e-05   2900.69 9099.31   12000   36.2587
## 44                      16 1.11136e-05   2624.65 9375.35   12000   32.8082
## 45                      17 3.02882e-06   2374.89 9625.11   12000   29.6861
## 46                      18 8.25449e-07   2148.89 9851.11   12000   26.8611
## 47                      19 2.24961e-07   1944.39 10055.6   12000   24.3049
## 48                      20 6.13091e-08   1759.36 10240.6   12000    21.992
## 49                      21 1.67087e-08   1591.93 10408.1   12000   19.8992
## 50                      22 4.55364e-09   1440.44 10559.6   12000   18.0055
## 51                      23 1.24101e-09   1303.36 10696.6   12000   16.2921
## 52                      24 3.38215e-10   1179.33 10820.7   12000   14.7417

Since the default function cannot load the output file with multi simulation result under ideal format to perform further analysis, we can use our built-in function readsims() to re-analyze the output object.

exp_1 <- readsims(out.1cpt, exp = 1)
exp_2 <- readsims(out.1cpt, exp = 2)
head(exp_1)

##   Time Oral_input A_central   A_elim A_total C_central
## 1    0 6000.00000      0.00    0.000    0.00    0.0000
## 2    1 1635.19000   4109.99  254.823 4364.81   68.4998
## 3    2  445.64100   4838.97  715.387 5554.36   80.6495
## 4    3  121.45100   4683.75 1194.800 5878.55   78.0624
## 5    4   33.09940   4321.22 1645.680 5966.90   72.0204
## 6    5    9.02064   3932.68 2058.300 5990.98   65.5446

head(exp_2)

##   Time Oral_input A_central   A_elim  A_total C_central
## 1    0 12000.0000      0.00    0.000     0.00    0.0000
## 2    1  3270.3800   8219.97  509.645  8729.62  102.7500
## 3    2   891.2830   9677.94 1430.770 11108.70  120.9740
## 4    3   242.9030   9367.49 2389.610 11757.10  117.0940
## 5    4    66.1988   8642.45 3291.360 11933.80  108.0310
## 6    5    18.0413   7865.35 4116.600 11982.00   98.3169

Plot the result of oral_input

plot(exp_2$Oral_input, type = "l", col = 2, xlab = "Time (h)", ylab = "Dose (mg)")
lines(exp_1$Oral_input, type = "l")

and C_central.

plot(exp_2$Time, exp_2$C_central, type = "l", col = "red", xlab = "Time (hr)", ylab = "Concentratiotn (mg/L)")
lines(exp_1$Time, exp_1$C_central)

4.2 Exercise

In perc.in.R, add the second simulation group that has an exposure concentration of 144 ppm. Besides, add Pct_metabolized in both simulation sections.

out42 <- mcsim(model = "perc.model.R", input = "perc2.in.R", dir = "modeling/perc")

## Execute: ./mcsim.perc.model.R.exe modeling/perc/perc2.in.R

Load the data to two objects,

sim_1 <- readsims(out42, exp = 1)
sim_2 <- readsims(out42, exp = 2)
sim_1

##       Time   C_exh_ug     C_ven Pct_metabolized
## 1    239.9 282.580000 1.9999800        0.239037
## 2    245.0  97.737100 1.8346600        0.244274
## 3    270.0  78.504900 1.4736500        0.269434
## 4    360.0  45.480900 0.8537420        0.356377
## 5   1320.0   4.268720 0.0801300        0.916256
## 6   2760.0   2.913030 0.0546818        1.353910
## 7   4260.0   2.261710 0.0424556        1.722850
## 8   5700.0   1.774920 0.0333177        2.012090
## 9   8580.0   1.093980 0.0205355        2.437360
## 10 10020.0   0.859231 0.0161290        2.590030

sim_2

##       Time  C_exh_ug     C_ven Pct_metabolized
## 1    239.9 567.10000 4.0363700        0.124866
## 2    245.0 197.42600 3.7059600        0.127567
## 3    270.0 158.99600 2.9845900        0.140670
## 4    360.0  93.00910 1.7459100        0.186878
## 5   1320.0   9.24358 0.1735150        0.555255
## 6   2760.0   6.25813 0.1174740        0.917613
## 7   4260.0   4.86477 0.0913187        1.237990
## 8   5700.0   3.81868 0.0716821        1.499000
## 9   8580.0   2.35107 0.0441330        1.900070
## 10 10020.0   1.84498 0.0346329        2.049430

Plot time vs exhaled air

plot(sim_2$Time, sim_2$C_exh_ug, type = "l", log = "y", col = "red", 
     xlab = "Time (hr)", ylab = "Exhaled air (ug)")
lines(sim_1$Time, sim_1$C_exh_ug)
legend("topright", legend = c("144 ppm", "72 ppm"), col=c(2, 1), lty=1)

Plot time vs blood concentration

plot(sim_2$Time, sim_2$C_ven, type = "l", log = "y", col = "red", 
     xlab = "Time (hr)", ylab = "Blood concentration (mg/L)")
lines(sim_1$Time, sim_1$C_ven)
legend("topright", legend = c("144 ppm", "72 ppm"), col=c(2, 1), lty=1)

Plot time vs Percent metabolized

plot(sim_1$Time, sim_1$Pct_metabolized, type = "l", log = "y", col = "red", 
     xlab = "Time (hr)", ylab = "Percent metabolized")
lines(sim_2$Time, sim_2$Pct_metabolized)
legend("bottomright", legend = c("144 ppm", "72 ppm"), col=c(2, 1), lty=1)

5 Example: Ethylbenzene (EB) PBPK model

In the EB-PBPK example, we generated three exposure scenarios in the input file, which were (1) 8-hr continuous exposure to 100 ppm EB, (2) 8-hr continuous exposure to 1 mg/m³ EB, and (3) 4-hr continuous exposure to 100 ppm EB. For the first and second scenarios, we’ll compare the EB concentration in arterial blood (Cart), pulmonary (Cvipu), and venous blood (Cvtot). The third case will be used to compare with the experimental data from the published study (Tardif et al., 1996).

5.1 Modeling and simulation

Same as PERC-PBPK case, we firstly use mcsim() to simulate the pharmacokinetics under the default parameter setting for the three given scenarios.

out <- mcsim(model = "EB.model.R", input = "EB.in.R", dir = "modeling/EB")

## * Created executable program 'mcsim.EB.model.R.exe'.

## Execute: ./mcsim.EB.model.R.exe modeling/EB/EB.in.R

Use readsims() to re-analyze the output object:

exp_1 <- readsims(out, exp = 1) 
exp_2 <- readsims(out, exp = 2) 
exp_3 <- readsims(out, exp = 3) 
nrow(exp_1)

## [1] 481

nrow(exp_2)

## [1] 481

nrow(exp_3)

## [1] 61

According to the above information, we can find that we have a relative “long” output format than the previous case for PERC. Therefore, we used head() and tail() to check the output and prevent the long report. The assignment n can use any integer to change the length of data.

head(exp_1)

##   Time        Cart       Cvipu       Cvtot
## 1  0.0 0.00000e+00 4.29882e-06 0.00000e+00
## 2  0.1 7.26594e-06 7.65351e-06 3.43935e-06
## 3  0.2 7.87140e-06 8.28883e-06 4.09070e-06
## 4  0.3 8.22359e-06 8.65849e-06 4.46969e-06
## 5  0.4 8.45244e-06 8.89881e-06 4.71608e-06
## 6  0.5 8.61830e-06 9.07306e-06 4.89472e-06

head(exp_2)

##   Time        Cart       Cvipu       Cvtot
## 1  0.0 0.00000e+00 9.88729e-12 0.00000e+00
## 2  0.1 1.55497e-11 1.63675e-11 6.64377e-12
## 3  0.2 1.69320e-11 1.77011e-11 8.01098e-12
## 4  0.3 1.74813e-11 1.83864e-11 8.71360e-12
## 5  0.4 1.79426e-11 1.88599e-11 9.19905e-12
## 6  0.5 1.82155e-11 1.91822e-11 9.52949e-12

tail(exp_3, n = 10)

##    Time       Cvtot
## 52  5.1 5.12186e-07
## 53  5.2 4.58776e-07
## 54  5.3 4.10999e-07
## 55  5.4 3.68233e-07
## 56  5.5 3.29935e-07
## 57  5.6 2.95634e-07
## 58  5.7 2.64907e-07
## 59  5.8 2.37379e-07
## 60  5.9 2.12717e-07
## 61  6.0 1.90621e-07

Same as the previous case, use apply() to find the maximum concentration in scenario 1 and 2.

apply(exp_1 , MARGIN = 2, FUN = max)

##        Time        Cart       Cvipu       Cvtot 
## 4.80000e+01 9.78695e-06 1.03014e-05 6.15409e-06

apply(exp_2 , MARGIN = 2, FUN = max)

##        Time        Cart       Cvipu       Cvtot 
## 4.80000e+01 2.02248e-11 2.13077e-11 1.17087e-11

For case 3, we need to further transfer the original unit (mol/L) to mg/L that can compare with the result in the published paper.

MW <- 106.16 # g/mol
exp_3$Cvtot_mg <- exp_3$Cvtot * MW * 1000 # mol/L -> mg/L 
apply(exp_3 , MARGIN = 2, FUN = max)

##         Time        Cvtot     Cvtot_mg 
## 6.000000e+00 6.123350e-06 6.500548e-01

After the data manipulation, we can further visualize our simulation result to compare the EB concentration sin arterial blood (Cart), pulmonary (Cvipu), and venous blood (Cvtot) under 8-hr continuous exposure for

100 ppm, and

plot(exp_1$Time, exp_1$Cvipu, col = "green", type = "l", main = "Exposure: 100 ppm", 
     xlab = "Time (hr)", ylab = "Concentration (mol/L)")
lines(exp_1$Time, exp_1$Cart, col = "red", type = "l")
lines(exp_1$Time, exp_1$Cvtot, col = "blue", type = "l")
legend("topright", legend = c("Cvipu", "Cart", "Cvtot"), col=c("green", "red","blue"), lty=1)

1 mg/m³

plot(exp_2$Time, exp_2$Cvipu, col = "green", type = "l", main = "Exposure: 1 mg/m3", 
     xlab = "Time (hr)", ylab = "Concentration (mol/L)")
lines(exp_2$Time, exp_2$Cart, col = "red", type = "l")
lines(exp_2$Time, exp_2$Cvtot, col = "blue", type = "l")
legend("topright", legend = c("Cvipu", "Cart", "Cvtot"), col=c("green", "red","blue"), lty=1)

Finally, we want to compare the simulation result with experiment data as:

data_x <- c(4.02, 4.5, 5, 5.5, 6)
data_y <- c(1.93, 1.29, 0.87, 0.55, 0.38)
sd <- c(0.15, 0.15, 0.24, 0.11, 0.05)
sd.up <- data_y + sd
sd.dn <- data_y - sd
plot(data_x, data_y, main = "Exposure: 100 ppm", xlab = "Time (hr)", ylab = "Concentration (mg/L)",
     xlim = c(0, 6), ylim = c(0, 2.5))
lines(exp_3$Time, exp_3$Cvtot_mg)
arrows(data_x, sd.dn, data_x, sd.up, code=3, length=0.02, angle=90)
legend("topright", legend = c("Data", "Model"), lty = c(NA, 1), pch = c(1, NA))

According to our simulation and experiment data, we found that the current parameter setting can’t provide the suitable prediction with the experiment data. Therefore, we might need to calibrate our parameter setting to perfom further modeling and simulation.

5.2 Exercise

5.2.1 Run EB-PBPK model under 100 ppm exposure for 4 hours and plot the time-course of blood concentration from 0 to 6 hour.

Modeling

out <- mcsim(model = "EB.model.R", input = "EB_exercise_1.in.R", dir = "modeling/EB")

## Execute: ./mcsim.EB.model.R.exe modeling/EB/EB_exercise_1.in.R

head(out, 10)

##    Time       Cvtot
## 1  0.00 0.00000e+00
## 2  0.01 1.94370e-06
## 3  0.02 2.44038e-06
## 4  0.03 2.65216e-06
## 5  0.04 2.79988e-06
## 6  0.05 2.92721e-06
## 7  0.06 3.04388e-06
## 8  0.07 3.15254e-06
## 9  0.08 3.25449e-06
## 10 0.09 3.34993e-06

Data manipulate
Visualize

plot(out$Time, out$Cvtot, type = "l", col = "red", 
     xlab = "Time (hr)", ylab = "Concentration (mol/L)")

5.2.2 Estimate the steady-state of arterial and venous blood concentrations associated with EB exposures (0.1 ppm to 1000 ppm).

Modeling

out <- mcsim(model = "EB.model.R", input = "EB_exercise_2.in.R", dir = "modeling/EB")

## Warning in readLines(paste0(dir, "/", input)): incomplete final line found
## on 'modeling/EB/EB_exercise_2.in.R'

## Execute: ./mcsim.EB.model.R.exe modeling/EB/EB_exercise_2.in.R

tail(out, 10)

##     Time     Cart       Cvtot
## 385   87 0.000782 0.000766287
## 386   88 0.000782 0.000766287
## 387   89 0.000782 0.000766287
## 388   90 0.000782 0.000766287
## 389   91 0.000782 0.000766287
## 390   92 0.000782 0.000766287
## 391   93 0.000782 0.000766287
## 392   94 0.000782 0.000766287
## 393   95 0.000782 0.000766287
## 394   96 0.000782 0.000766287

Data manipulate

for(i in 1:4){
  exp <- readsims(out, exp = i) %>% tail(1)
  exp$Conc <- 10^(i-2)
  if (i == 1) OUT <- exp else OUT <- rbind(OUT, exp)
}
OUT

##     Time        Cart       Cvtot  Conc
## 97    96 8.80149e-08 5.09942e-08   0.1
## 971   96 8.87348e-07 5.17646e-07   1.0
## 972   96 9.78738e-06 6.15455e-06  10.0
## 973   96 7.82000e-04 7.66287e-04 100.0

Visualize

plot(OUT$Conc, OUT$Cart, type = "l", log = "xy", col = "red", 
     xlab = "Exposure (ppm)", ylab = "Concentration (mol/L)")
lines(OUT$Conc, OUT$Cvtot, type = "l", col = "blue")
legend("bottomright", legend = c("Arterial", "Venous"), col=c("red","blue"), lty=1)

5.2.3 Construct the relationships for the estimated inhalation exposure level and the fraction of EB metabolized after 8-hr continuous exposure. In addition, estimate the percentage metabolized from liver, lung, and richly perfused tissue.

Modeling

out <- mcsim(model = "EB_v2.model.R", input = "EB_exercise_3.in.R", dir = "modeling/EB")

## * Created executable program 'mcsim.EB_v2.model.R.exe'.

## Warning in readLines(paste0(dir, "/", input)): incomplete final line found
## on 'modeling/EB/EB_exercise_3.in.R'

## Execute: ./mcsim.EB_v2.model.R.exe modeling/EB/EB_exercise_3.in.R

tail(out, 10)

##    Time         Ain     Amet_Rl    Amet_Rlu   Amet_Rvrg        Amet
## 65  3.5 0.000164746 1.98351e-05 3.48149e-05 5.19131e-05 0.000106563
## 66    4  0.00018524  2.2714e-05 4.03045e-05 5.95745e-05 0.000122593
## 67  4.5 0.000205245 2.55956e-05 4.58298e-05 6.72523e-05 0.000138678
## 68    5 0.000224819 2.84794e-05 5.13836e-05 7.49431e-05 0.000154806
## 69  5.5 0.000244014  3.1365e-05 5.69604e-05 8.26445e-05  0.00017097
## 70    6 0.000262875  3.4252e-05 6.25561e-05 9.03547e-05 0.000187163
## 71  6.5 0.000281442 3.71403e-05 6.81674e-05 9.80721e-05  0.00020338
## 72    7 0.000299749 4.00296e-05 7.37918e-05 0.000105796 0.000219617
## 73  7.5 0.000317828 4.29197e-05 7.94273e-05 0.000113524 0.000235871
## 74    8 0.000335705 4.58106e-05 8.50721e-05 0.000121257  0.00025214

Data manipulate

for(i in 1:4){
  exp <- readsims(out, exp = i) %>% tail(1)
  exp$Conc <- 10^(i-1)
  if (i == 1) OUT <- exp else OUT <- rbind(OUT, exp)
}
OUT$`% met` <- OUT$Amet / OUT$Ain
OUT$`% met_Rl` <- OUT$Amet_Rl / OUT$Ain
OUT$`% met_Rlu` <- OUT$Amet_Rlu / OUT$Ain
OUT$`% met_Rvrg` <- OUT$Amet_Rvrg / OUT$Ain
OUT

##     Time         Ain     Amet_Rl    Amet_Rlu   Amet_Rvrg        Amet Conc
## 17     8 4.53412e-07 1.48943e-07 5.00317e-08 2.44063e-07 4.43038e-07    1
## 171    8 4.53218e-06 1.48523e-06 5.03517e-07 2.43888e-06 4.42763e-06   10
## 172    8 4.50796e-05 1.42393e-05 5.43677e-06 2.42451e-05 4.39211e-05  100
## 173    8 3.35705e-04 4.58106e-05 8.50721e-05 1.21257e-04 2.52140e-04 1000
##         % met  % met_Rl % met_Rlu % met_Rvrg
## 17  0.9771201 0.3284937 0.1103449  0.5382809
## 171 0.9769316 0.3277076 0.1110982  0.5381251
## 172 0.9743010 0.3158701 0.1206038  0.5378286
## 173 0.7510761 0.1364609 0.2534133  0.3612011

Visualize

plot(OUT$Conc, OUT$`% met`, type = "l", log = "x", col = "red", ylim = c(0,1), lwd = 2,
     xlab = "Exposure (ppm)", ylab = "Percent metabolized")
lines(OUT$Conc, OUT$`% met_Rl`)
lines(OUT$Conc, OUT$`% met_Rlu`)
lines(OUT$Conc, OUT$`% met_Rvrg`)
text_x <- head(OUT$Conc, 2)[2]
text(x = text_x, y = head(OUT$`% met_Rl`, 1)+0.05, label = "Liver")
text(x = text_x, y = head(OUT$`% met_Rlu`, 1)+0.05, label = "Lung")
text(x = text_x, y = head(OUT$`% met_Rvrg`, 1)+0.05, label = "Richly perfused tissue")

5.2.4 Add additional exposure routes include oral ingestion in the EB-PBPK model and estimate the Cmax and Tmax after received a single gavage dose of 180 mg/kg.

Modeling

out <- mcsim(model = "EB_v2.model.R", input = "EB_exercise_4.in.R", dir = "modeling/EB")

## Execute: ./mcsim.EB_v2.model.R.exe modeling/EB/EB_exercise_4.in.R

exp <- readsims(out, exp = 1)
head(exp)

##   Time   Dmgkg   Cvtot
## 1  0.0 180.000       0
## 2  0.1 162.871 2640660
## 3  0.2 147.372 3360960
## 4  0.3 133.347 3788410
## 5  0.4 120.658 4053310
## 6  0.5 109.176 4212760

Data manipulate

C_max <- max(exp$Cvtot)
T_max <- exp$Time[which(exp$Cvtot == C_max)]
C_max

## [1] 4358230

T_max

## [1] 0.8

Visualize

plot(exp$Time, exp$Cvtot, type = "l", col = "red",
      xlab = "Time (hr)", ylab = "Concentration (mol/L)")
abline(v = T_max, lty = "dotted")
abline(h = C_max, lty = "dotted")

Reference

Bois, F.Y., Gelman A., Jiang J., Maszle D.R., Zeise L., Alexeef G. (1996). Population toxicokinetics of tetrachloroethylene.
Chiu, WA and Bois, Y. (2006). Revisiting the population toxicokinetics of tetrachloroethylene.
Gelman A., Bois, F.Y., Jiang J. (1996). Physiological Pharmacokinetic Analysis Using Population Modeling and Informative Prior Distributions.
GNU MCSim user manual
Soetaert, K., et al. (2018). deSolve: Solvers for Initial Value Problems of Differential Equations (‘ODE’, ‘DAE’, ‘DDE’).
Tardif, R., Charest-Tardif, G., and Brodeur, J. (1996). Comparison of the influence of binary mixtures versus a ternary mixture of inhaled aromatic hydrocarbons on their blood kinetics in the rat.

Tutorial 1 - Walk-through of Working Models

Nan-Hung Hsieh

2019/04/25 (update: 2019-05-29)

0 Overview

1 Create PBPK model

1.1 Define the state variables `States{}`

1.2 Define the outputs variables `Outputs{}`

1.3 Define the input variable(s) `Inputs{}`

1.4 Parameters of the model

1.5 Define the parameter initialization and scaling `Initialize{}`

1.6 Define the dynamics of the simulation `Dynamics{}`

1.7 Define the output calculations `CalcOutputs{}`

2 Create input/simulation file

2.1 Setting integrator `Integrate()`

2.2 Setting paramerer(s)

2.3 Setting simulation `simulation{}`

2.4 Working with the R package deSolve

3 Analyze GNU MCSim output with R

3.1 Loading output

3.2 Assign variables

3.3 Manipulating output

3.4 Plotting

4 One-compartment model

4.1 Modeling and simulation

4.2 Exercise

5 Example: Ethylbenzene (EB) PBPK model

5.1 Modeling and simulation

5.2 Exercise

5.2.1 Run EB-PBPK model under 100 ppm exposure for 4 hours and plot the time-course of blood concentration from 0 to 6 hour.

5.2.2 Estimate the steady-state of arterial and venous blood concentrations associated with EB exposures (0.1 ppm to 1000 ppm).

5.2.3 Construct the relationships for the estimated inhalation exposure level and the fraction of EB metabolized after 8-hr continuous exposure. In addition, estimate the percentage metabolized from liver, lung, and richly perfused tissue.

5.2.4 Add additional exposure routes include oral ingestion in the EB-PBPK model and estimate the Cmax and Tmax after received a single gavage dose of 180 mg/kg.

Reference

Tutorial 1 - Walk-through of Working Models

Nan-Hung Hsieh

2019/04/25 (update: 2019-05-29)

0 Overview

1 Create PBPK model

1.1 Define the state variables States{}

1.2 Define the outputs variables Outputs{}

1.3 Define the input variable(s) Inputs{}

1.4 Parameters of the model

1.5 Define the parameter initialization and scaling Initialize{}

1.6 Define the dynamics of the simulation Dynamics{}

1.7 Define the output calculations CalcOutputs{}

2 Create input/simulation file

2.1 Setting integrator Integrate()

2.2 Setting paramerer(s)

2.3 Setting simulation simulation{}

2.4 Working with the R package deSolve

3 Analyze GNU MCSim output with R

3.1 Loading output

3.2 Assign variables

3.3 Manipulating output

3.4 Plotting

4 One-compartment model

4.1 Modeling and simulation

4.2 Exercise

5 Example: Ethylbenzene (EB) PBPK model

5.1 Modeling and simulation

5.2 Exercise

5.2.1 Run EB-PBPK model under 100 ppm exposure for 4 hours and plot the time-course of blood concentration from 0 to 6 hour.

5.2.2 Estimate the steady-state of arterial and venous blood concentrations associated with EB exposures (0.1 ppm to 1000 ppm).

5.2.3 Construct the relationships for the estimated inhalation exposure level and the fraction of EB metabolized after 8-hr continuous exposure. In addition, estimate the percentage metabolized from liver, lung, and richly perfused tissue.

5.2.4 Add additional exposure routes include oral ingestion in the EB-PBPK model and estimate the Cmax and Tmax after received a single gavage dose of 180 mg/kg.

Reference

1.1 Define the state variables `States{}`

1.2 Define the outputs variables `Outputs{}`

1.3 Define the input variable(s) `Inputs{}`

1.5 Define the parameter initialization and scaling `Initialize{}`

1.6 Define the dynamics of the simulation `Dynamics{}`

1.7 Define the output calculations `CalcOutputs{}`

2.1 Setting integrator `Integrate()`

2.3 Setting simulation `simulation{}`