A methodology for solving physiologically based pharmacokinetic models without the use of simulation softwares

Toxicology Letters Toxicology Letters 85 (1996) 113-126

A methodology for solving physiologically based pharmacokinetic models without the use of simulation softwares Sami Haddad, Michael Pelekis, Kannan Krishnan* Dkpartement de mkdecine du travail et d’hygikne du milieu, Fact&C de mtidecine. UniversitPde Montrkal, 2375 Cfite Ste-Catherine, Mont&al, Quibec, Canada. H3T IA8 Received 27 November 1995; revised 19 January 1996; accepted 22 January 1996

Abstract

The objective of the: present study was to develop and validate a methodology for solving physiologically based pharmacokinetic (PBPK) models without the use of simulation software. The approach involves keying the parameter values and model equations into Microsoft Excel@ spreadsheets, and conducting simulations by solving the model equations according to Euler’s method of numerical integration. This approach was applied to simulate the pharmacokinetics of styrene in rats exposed to 80 and 600 ppm for 6 h. The simulation results were plotted along with experimental data using the regular graphic features available in Excel@,and validated by comparing them with simulation results obtained using a commercially available software (Advanced Continuous Simulation Language, ACSLe). The simulations obtained with ACSL@ and Excel@, in general, differed by < 1%. The methodology developed in the present study should lhelp informed individuals understand and solve PBPK models, without having to use ‘black-box’ kind of computer programs and simulation softwares. Keywork

PBPK modeling; Spreadsheet programs; Styrene; Physiological pharmacokinetics

1. IntroductIorl

Risk assessors are confronted with the increasing number of potentially useful mathematical models for conducting extrapolations essential for dose-response and exposure assessments. Whereas some of these models are sold as ‘ready-made’ packages, others are available to the user as programs written in a specific programming or simu* Corresponding author, Tel.: 514 343 6581; Fax: 514 343 2200.

lation language. A number of these models look like a ‘black-box’ to the analyst, who gets to see the end results but not the temporal evolution of solutions to the complex mathematical formulations constituting the basis of end results generated. The ‘internal mechanics’ of computer simulation in such cases can be visualized if the user can reconstruct the way in which the software (1) solves each equation in the model, and (2) takes the output of one equation and provides it as input to other equations of the model. Approaches allowing such reconstructions would enable a bet-

0378-4274/96/$15.00 c 1996 Elsevier Science Ireland Ltd. All rights reserved PII SO378-4274(96)03648-X

114

S. Ha&d

et al. /Toxicology

ter understanding of the manner in which the model equations are solved within simulation software to provide the end results of interest to the user. Physiologically based pharmacokinetic (PBPK) models are among those tools that are increasingly finding use in risk assessment applications. The PBPK models published so far have been written and solved using either commercially available simulation software packages (e.g., ACSLe, SCoP@, STELLAe [l-4]), or other software products that allow programming of a specific nature [5,6]). These programs/packages are easily accessed and understood by individuals who are familiar with the techniques of programming or mechanics of simulation. However, a large majority of toxicologists are not yet in this category and therefore should benefit from the availability of an altemative, simpler methodology for solving PBPK models. In this context, we chose to explore the feasibility of solving PBPK model equations using a spreadsheet program. The spreadsheet programs such as Lotus l-2-3@‘,QuattroProe and Excel’s’are among the most commonly used programs by biologists and toxicologists to organize, transform, and plot experimental data. By specifying and solving each PBPK model equation, and establishing the linkages among the various model equations within spreadsheets, it should be possible to generate simulations of the pharmacokinetics of chemicals. The objective of the present study was to develop a methodology for solving PBPK models using a spreadsheet program, namely, Excel@. 2. Approach The approach involved (1) writing physiologically based mathematical expressions of the pharmacokinetics of a volatile organic chemical (styrene), (2) keying the equations into Excel@ spreadsheets, (3) conducting simulations by solving the model equations, and (4) comparing the results with simulations obtained using ACSL@ (Advanced Continuous Simulation Language, Mitchell and Gauthier Inc., Concord, MA).

Letters 85 (19%)

113-126

2.1. Equations describing the pharmacokinetics of inhaled styrene The processes of absorption, distribution, and clearance of inhaled styrene were described mathematically, based on hypothetical interrelationships among certain critical parameters [ 11. Accordingly, the arterial blood concentration (C,) of styrene, whose pulmonary metabolism and retention are assumed to be negligible, was estimated as follows [l]:

c, = Qc G + Qp Cinh Qc+ Qdpb

where

(1)

Q, = cardiac output (l/h), QP= alveolar ventilation rate (lib), C, = venous blood concentration (mg/l), Cinh= inhaled air concentration (mg/I), Pb = blood:air pahtiOn COCfkient.

The rate of change in the amount of styrene in non-metabolizing tissues (i) (slowly perfused tissues (s), richly perfused tissues (r), fat (f)) was described as a perfusion-limited process as per the following mass balance differential equation (MBDE):

ddildt = Qi (C, - C,)

(2)

Ai = amount in tissue (mg/l), Vi = volume Of tissue (l), Pi = tissue:blood partition coefficient. The rate of change in the amount of styrene in metabolizing tissues (i.e., liver in the present example) (dAl/dt) was described as follows: dAl/dt = Ql (C, - C,,) - dAmet/dt

(4)

where (A, = amount in liver (mg), Q, = liver blood flow rate (l/h), CVl=concentration in venous blood leaving liver

=-

Al VIPI

(5)

S. Haa&d et al. /Toxicology Letters 85 (19%) 113-126

where PI = 1iver:lblood partition coefficient, Vt = volume of liver (I), dAmet/dt = rate of change in the amount metabolized (mgh) =

vmaxC”I K, + c,,

where Fg/$, m . The venous c

_

v-

(6)

V,, = maximal velocity for metabolism K,,, = Michaelis-Menten affinity constant concentration of styrene in the mixed blood (C,) was calculated as follows:

QIGI + Qr Gr + QsGs + Qf c, QC

(7)

In essence, these 7 equations represent the PBPK model for inhaled styrene [I]. The solution to these equations simulates inhalation pharmacokinetics of styrene in the species of interest. Examining these equations, it is evident that they

Fig. 1. A schematic of the PBPK model for styrene. In this model, the rat is represented as a 4-compartmental system interconnected by systemic circulation. The input for the system is the product of the inhaled concentration of styrene (C& times the alveolar ventilation rate (Q,). The resulting C, is in turn provided as the input to the tissue compartments, the effluent venous blood concentrations of which are provided as input for the calculation of mixed venous concentration (Cv). t refers to the time interval of integration. All other abbreviations are defined in Table I.

115

can be linked in some logical sequence such that the output of one equation becomes the input for another equation, and so forth (Fig. 1). 2.2. Keying the PBPK model parameters and equations into spreacisheets The ‘spreadsheet’ is a format in which a number of cells are found. Each cell within a worksheet is tagged with a column and row co-ordinate, designated by a letter and number, respectively. For example, a cell referred to as Bl signifies that this cell is placed at the junction of the second vertical line and first horizontal line in the spreadsheet. Using this basic feature of spreadsheets, the essential links to conduct PBPK model simulations can potentially be developed. A fundamental requirement for any useful simulation framework is that it should consist of features that permit: (1) the introduction of numerical values of model parameters at the beginning of each simulation exercise, and (2) the establishment of output-input links between the various equations constituting the model. These aspects can be accomplished with spreadsheet programs since they permit the use of absolute and relative reference expressions. The absolute reference expression is one which involves referring to the content of a particular cell in any other part of the spreadsheet by providing the cell-specific ‘address’ or label. This option can be used for referring to particular parameters anywhere in a spreadsheet. The relative reference expression involves referring to a cell according to its location relative to another cell where the calculation is being carried out. This option is particularly useful in situations when the output of an equation contained in one cell is to be provided as input for equation contained in another (e.g., adjacent) cell. Thus, the use of relative reference expressions in spreadsheets may be useful to facilitate loop-type calculations essential for advancing the state of a system during simulations. 2.2.1. Entry of model parameter values For solving the equations presented in Fig. 1, the numerical values of several physiological,

116

S. Ha&d

et al. /Toxicology Letters 85 (19%)

physicochemical, and biochemical parameters are required. The parameter values were entered into the spreadsheet in such a way that the numerical value of each parameter resides in a specific location, i.e., cell-specific address. Additionally, a particular cell containing a parameter value may be referred to with a combination of letters (e.g., cell C5 was identified as Qr, in the present study). In fact, in the present study, all model parameters were referred to by a combination of letters (Table 1). For example, the alveolar ventilation rate is referred to as Q,,, and whenever Qr is typed in any other cell of the spreadsheet, the numerical value found in cell C5 will be imported automatically. There is another way of referring to the parameter values in the spreadsheets. This procedure involves the use of $ sign before the column and row numbers of a cell. For example, let us say that the location of the numerical value for Q, is referred to as $C!W.Whenever SC%4is specified in an equation in any section of the spreadsheet, the numerical value of Q, found in that particular cell will be imported. However, in the present study,

113-126

the PBPK model parameters were referred to by a combination of letters as detailed in the preceding paragraph, and not by specifying the cell addresses. 2.2.2. Entry of model equations and the interconnections among them

Table 2 lists the manner in which the PBPK model equations represented in Fig. 1 were written in Excel@ spreadsheets. Aside the venous blood and gas-exchange lung compartments, 4 equations per tissue compartment were written. These equations correspond to the tracking of: (1) the rate of change in the amount of chemical in the tissue, (2) amount of chemical in tissue, (3) concentration of chemical in tissue, and (4) concentration of chemical in the venous blood leaving tissue. The form of these equations written for Excel’s’ (Table 2) contains model parameters referred to by abbreviations (listed in Table 1) and variables referred to with the use of relative reference expressions. Table 3 presents a portion of the spreadsheet depicting how the equations shown in Table 2 are

Table I List of parameters of the rat styrene PBPK model, their numerical values and location in the Excel@ spreadsheet Parameters

Cardiac output Alveolar ventilation rate Blood flow to liver Blood flow to fat Blood flow to richly perfused tissues Blood flow to slowly perfused tissues Volume of liver Volume of fat Volume of richly perfused tissues Volume of slowly perfused tissues Liver:blood partition coefftcient Fat:blood partition coeflicient Richly perfused tissues:blood partition coefftcient Slowly perfused tissues:blood partition coefticient Bloodair partition coefficient Maximal velocity of metabolism Michaelis affinity constant

Abbreviations

Numerical values b

Cell locationc

5.64 Iih 4.5 lnl 2.11 lib

c4 C5 C6 Cl C8 c9 D6 D7 D8

0.5I/h

2.35 lib 0.68 lib 0.012 I 0.027 I 0.015 I 0.219 I 2.7 50 5.7

1 40

3.6 mgIh 0.36 mgil

D9 E6 El E8 E9 FIO G6 H6

‘%e various model parameters are referred to, with the use of these abbreviations in the spreadsheet. bAll parameter estimates were obtained from Ramsey and Andersen [I]. CThecell locations provided here correspond to the column and row co-ordinates respectively, i.e., the alphabetical letters denote the columns and the arabic numerals correspond to the rows in the spreadsheet.

Q, x (Cr.,, - G,, _ ,)

C”” =

Q, x G,

C vsn = Cs,IPs

Csn = AsnK

+ Qf x Cvr” + Or x Gr,+Qs QC

Ak = dAddt,, x t + A,,, _ ,

dAddt. =

C, = Ar.IVr G, n = C,lP,

t+A,“_,

Qr X G,, - qrn _ ,)

Ar,,=dAidt,x

dAJdt, =

x Gsn

u39 =

ml

I’,, x H38 K + H38

QC

Q, x H39 + Qf x L39 + Qr x P39 + Q, x T39

439 = Q, x (D39 - T38) R39=Q39x t+R38 S39 = R39/V, T39 = S39lP,

P39 = 039/P,

039 = N39/Vr

N39 = M39 x t + N38

M39 = Q, x (D39 - P3B)

L39 = K39/P,

K39 = J39/Vf

139 = Qf x (D39 - L38)

F39 = E39 x t + F38 G39 = F394 H39 = G39/P,

Cfn = 4/r

“1”- I

E39 = Q, x (D39 - H38) -

Q, x U38 + Q, x C39 Q<+ 0-I -y-uPL

J39 = 139 x t + J38

Qr X (can - cvrn _ ,)

m

x C”l+,

D39 =

Expression in Excel@b

Am = dA,4dt, x t + A,,, _ ,

dA,Jdt,, =

C”l. = Cl”4

C,,, = AkV,

x t + A,, _ ,

QI 7 (can - CVI,, _ ,) - ;yc

A,, = dA,ldt,

dA{dt,, =

Equationsa

“All abbreviations and symbols used in the equations, except n and n - 1, are defined in Fig. I. Subscripts n and n - 1 refer to the current and previous simulation times. The difference between n and n - I in the styrene example was 0.005 h. ,‘fhe components of these equations refer either to absolute references (in the case of constant input parameters, as defined in Table 1) or to relative references (in the case of state variables).

Venous blood

Slowly perfused tissues

Richly perfused tissues

Fat

Liver

Arterial blood

Compartment

Table 2 Equations used in the calculation of tissue, arterial and venous blood concentration of styrene, and their expression in Excel’s spreadsheet

r+B38

39

Cinh

= Cinh

=

Ctnh (mgfl)

C

tQ, x C39)Y tQc + (Q#‘d

=((Q, x U38)+

= (tQc x 0) + (Qp x C38)Y tQ, + (Q,&N

C,

D

= (Q, x (D39 H38)) V’,, x H38/ W, + H38))

K,, + 0))

(V,xXO~

= (Q, x (D38 - 0)) -

dA,ldr

Liver

E

Cl

G

= (E39 x r) = F39/ + F38 v,

= (E38 x r) = F38/ v,

Al

F

= G39/ P/

= G38/ p/

C”l

H

= (Qr x (D39 L38))

= (Qf x (D38-0))

dAddr

Fat

I

= 139 x f + 538

= 138x1

Ar

J

Vf

= J39/

= J38/Vf

Cf

K

Pf

= K39/

= K38/P,

Cvf

L

=(((Qi x H39) +

(Qf x L39) + (Q, x P39) + (Q, x T39)YQc)

tQr x L38)+ tQ, x P38) + tQs x TWYQJ

= (((Q, x H38)+

C”

U

The row and column co-ordinates are designated by an arabic numeral and an alphabetical letter, respectively. The equations found in rows 38 and 39 correspond to calculations at the first and second integration intervals. For continuing the simulation, the set of equations in row 39 should be copied onto desired number of subsequent rows. In this table, column B represents the state of the system which advances during each time interval (1). Column C contains the exposure concentration at any given time during the simulation, whereas columns D and U represent calculations of arterial and venous blood concentrations of the chemical. In between, 4 columns per compartment (e.g., liver: columns E-H, fat: columns I-L, slowly perfused tissues: columns M-P, and richly perfused tissues: columns Q-T) are devoted for calculations of (I) the rate of change in the amount of chemical in tissue, (2) amount of chemical in tissue, (3) concentration of chemical in tissue, and (4) concentration of chemical in venous blood leaving tissue. All abbreviations are defined in Table I.

I

Time (h) --

38

31

36

B

Table 3 A portion of the spreadsheet depicting the entry of model equations

~8 2 y s

:

9

&

i$ 2

F Oy 5.

2

k 4

.@ %

S. Hadhd

et al. / Toxicology Letters 85 (19%)

actually entered into spreadsheets. Accordingly, in the spreadsheet, the descriptions of the 4 tissue compartments occupy 16 columns (columns E-T, in Table 3), and calculation/representation of the simulation time, exposure concentration, arterial concentration, and venous concentration occupy one column each (columns B, C, D, and U, respectively). The mixed venous blood concentration resulting from the chemical concentration in venous blood leaving the various tissue compartments described in ,a particular row is calculated in the same row (i.e., cell U38 in row 38, Table 3), and this C, then is used along with Cinh (i.e., cell C39) to calculate C, in the subsequent row (e.g., cell D39). In this structure then, according to the schematics shown in Fig. 1 and Table 3, all model equations are interconnected by specifying the proven/hypothetical input-output connections among them. 2.3. Solving the d#erential equations To advance the state of the system (i.e., calculate the temporal evolution/change in the concentration of styrene in the various compartments of the model), the set of differential equations constituting the PBPK model should be solved. In the present study, the Euler method of approximation of the solution to first order differential equations was used. For the first order differential equations such as Eq. 2, which are typically employed in PBPK models, the change in the amount of chemical for a given dt (At), corresponding to the integration interval, should be estimated. The amount at the new t (i.e., lthe solution to the differential equation) may be approximated to the amount at the beginning of the integration interval plus the product of the slope given by the differential equation and the integration step size (At): new value = old value + (slope

x

At), or

current amount in tissue = previous amount in tissue + [(dAi/dt) x At]

I13-I26

119

depicted in Table 3 (e.g., cell F39). Here, the product of integration interval (or the step size, At) and the rate of change in the amount (or the slope, d4ildt) (cell E39) is added to the old value (found in cell F38) to estimate the amount in tissue at the current time t (cell F39). One has only to repeat the calculations shown in row 39 of Table 3 for each time interval of integration until the end of the desired time length for simulation. In the present study, the time interval for integration was fixed at 0.005 h. Therefore, the set of equations in row 39 in Table 3 was copied onto desired number of subsequent rows, depending on the required length of simulation and specified integration time interval. Each line in the Excel’s’ spreadsheet then represents calculations characterizing the state of the system at every 0.005 h. In the present study, simulations were conducted for 24 h using 0.005 h as the integration interval. Therefore, in total, 4800 (= 24/0.005) lines were used up for conducting PBPK simulations. Once (1) time interval for integration is specified, (2) required number of cells chosen, (3) numerical values of model parameters provided, and (4) equations in the first and subsequent rows of the spreadsheet entered, the simulation begins. Model simulation refers to the continuous solving of the set of MBDEs for the various tissue compartments for a defined exposure condition. The solution to the set of PBPK model equations is generated every time the numerical value in a cell corresponding to exposure concentration or model parameters is changed, since these cells are specified in one or more equations appearing in the spreadsheet (Table 3). In the present work, simulations were conducted for inhaled styrene (80 and 600 ppm) in the rat, using Microsoft Excele (version 5.0) in an IBMcompatible personal computer (Model 486-DX2, 16 MB RAM). All model parameters were obtained from Ramsey and Andersen [ 11. The simulation results were plotted against time using the regular graphic features available in Excel@. 2.4. Validation of the proposed methodology

The above expressions essentially represent the Euler method of integration [7]. The implementation of this procedure for solving the MBDEs is

The C, values for a 6-h exposure of rats to 80 and 600 ppm styrene obtained with the proposed

120

S. Haddad et al. / Toxicology Letters 85 (19%)

method were compared with those obtained using ACSL@ (version 11.2.1). Simulations of styrene pharmacokinetics in ACSL@ were obtained using Euler and stiff Gear algorithms. While using Gear’s method (which is a variable order, variable step size method) in ACSLe, the communication interval (CINT), the number of integration steps per communication interval (NSTP), maximum value of integration step size (MAXT), and minimum value of integration step size (MINT) were set equal to 0.01 h, 1000, 1.0 e+‘O, and 1.0 e-lo, respectively. The integration step size for variable algorithms is calculated as follows [8]: Step size = MAX(MINT, NSTP))

MIN(MAXT,CINT/

Accordingly, in the present example, the initial integration step size would be O.Ol/lOOO (i.e., CINTYNSTP), and this would subsequently be varied/increased until it reaches the most efficient size. The upper limit of the integration step size chosen by a variable step algorithm is bounded by the CINT. While using Euler algorithm in ACSLe, the integration step size is calculated by [8]: Step size = MIN(MAXT, CINT/NSTP) In the present study, the CINT, NTSP, and MAXT were set equal to 0.005 h, 1, and 0.005 h, respectively. Setting the NSTP to 1 and MAXT to the desired integration step size (i.e., 0.005 h) results in a decoupling of the control of step size from CINT [8]. In the present study, as per the formula presented above, a fixed step size of 0.005 h was used in ACSLe when Euler algorithm was chosen. In this particular case, both the integration step size (fixed) and the CINT are equal to 0.005 h. This was done simply to keep CINT and integration step size values identical to those used in the Excel@methodology. Setting CINT equal to integration step size means that the result for each integration step size (i.e., simulated data for every 0.005 h in the present study) will be displayed/communicated. 3. Results The solution to the set of differential equations constituting styrene PBPK model obtained using

113-126

the proposed method corresponds to data on the temporal evolution of the blood and tissue concentrations of styrene in the rat. A portion of these data, calculated and displayed anew for each interval of 0.005 h until 24 h, are shown for 80 ppm styrene exposure (Fig. 2). In Fig. 2 Excel’s’simulations of C, are also compared with the corresponding experimental data [ 11. The simulation results obtained with the present method for 80 and 600 ppm exposures are compared with corresponding experimental data and simulations obtained using Euler and Gear algorithms in ACSL@ in Tables 4 and 5. These results indicate that the simulations obtained using Excel@ and ACSL@ are quite comparable, and the difference between them, in most cases, is within 1%. It should be noted that the simulation results obtained using Gear algorithm in ACSL@ correspond to those that have been previously reported by Ramsey and Andersen [ 11. The integration step size and order specified by ACSLe using Gear algorithm were obtained for the 80 and 600 ppm exposure simulations. The Gear algorithm in ACSLe used an integration interval (or step size) of -0.01 h for the simulations (not shown), and frequently performed second order integration for solving the differential equations of styrene PBPK model for the 600 ppm exposure simulations and first order integration for 80 ppm exposure simulations (not shown). 4. Discussion Simulation is the system behavior predicted for specific exposure conditions, by solving the set of differential and algebraic equations representing the quantitative interrelationships among the various elements of a model. Typically, simulation is opted when (1) the real system does not exist, (2) the real system exists but experimentation is expensive or unethical, (3) a forecasting model is required to analyze events occurring during very long periods of time or complex scenarios, and (4) the model does not have practical, analytical solutions (e.g., nonlinear systems) [9]. In the case of PBPK modeling, simulation is opted because of the last 3 reasons. Conventionally, the use of simulation or programming languages has been sought

c

1,

0 656

0.672

0.340

0340

0.346

0346

0.340

0.340

0340

0.340

0.3

0.4

06

0.7

0.6

0.9

1

05

0 706

0.696

0 665

0 640

0.620

0 594

0.557

0340

0.2

0.266 0466

Ca

0.1

(ms!L)l

0.346 0.340

0005

IClnh

0.012

03

I

I

c

,

0.441

0.005 0 005 0.005

0.001

0.001

0.449

0.432

0.005

0.422

0.410

0001

0 005

0001

0 002

0 396

0.360

0.002 0.005

0 360

0 004 0.005

0.003

0.002

0.332

Cl

0.004

1

3.6

(mwhr

0.004

Al

LIVER

Vmax

I

I

0166

0163

0.160

0 156

0.152

0 147

0.141

0.133

0.123

0 067 0.104

1 CVI I

0.36

Km (m@L

Biochemical

0.234 0.261

1

Fb

r

0.003 0 003

0.561

dAUdt

2.7

Pi

PhyJicahemical

0.257

0.262

0.266

0.270

0.272

0 274

0.274

0.271

0.262

0.133 0.236

dAfIdt

1

0 260

0 234

0.206

0.161

0.154

0 127

0 099

0 072

0.045

1

FAT 0.001

Af

9 631

6 671

7.693

6 701

5.696

4.666

3 672

2.663

1.672

0.025 0.736

CT

I

m

0 173

0.154

0.134

0.114

0.094

0.073

0.053

0.033

0.000 0.015

~vf

1 I

,

0 006

0.009

0.011

0.012

0.014

0.017

0.021

0.030

0.056

0.171

0.625

I

RICHLY dAr/dt

I

I

”

,

r

0.060

0 059

0.056

0.057

0.056

0.054

0.052

0.050

0 046

0.036

0.003

Ar

I

PERFUSED

4.006

3.949

3 063

3 607

3.710

3.613

3.467

3.322

3.056

2.413

I

TISSUES 0.206

cr

0.703

0.693

0 661

0 666

0.652

0.634

0.612

0.563

0.537

0.423

0.037

CM

Styrene 80 ppm

n

I

1

, I

0 047

0.056

0.066

0.063

0.101

0.125

0.154

0.190

0.231

0.260

0.161

dAskn

I

SLDWLY

”

I

0.140

0 135

0.126

0.121

0.112

0.101

0.067

0.070

0.049

0.024

0001

AS

I

PERFUSED

n

I

0.636

0.615

0.567

0.552

0.511

0.460

0.396

0.319

0.223

0.110

I

TISSUES 0.004

Cs

s

0 636

0615

0.567

0.552

0.511

0.460

0.396

0.319

0223

I

I

1

0.110

0004

ti

I

”

0.44!

0.43!

0 421

0.41!

0.4Ol

0.30:

0.36:

0.331

0.29!

0.231

0.04.

Cv

Fig. 2. Print out of a computer screen depicting (1) the plot of simulated (solid line) and experimental (squares, data from Ref. [I]) data on the C, of styrene during and following a 6-h exposure of rats to 80 ppm of this chemical, (2) the numerical values of the PBPK model parameters, and (3) a portion of the raw numbers corresponding to the rate of change in the amount of styrene in tissue, amount in tissue, concentration in tissue, and venous blood concentration leaving the tissues, generated during simulations between time 0.005 and 24 h. Note that only the data in 1 out of every 20 rows were chosen for display using the ‘filter’ function in Excel@.

3e

35

Time (hr)

(t)

rttewal

45

eption

5.64

2.11

38

Y

a (L/Iv) V(L)

Lung (P) Lii (I)

!+I

,

Physiological

I

r

s

Tissue

Bodv (c)

,

*

3

2

S. Ho&d

122

et 01. /Toxicology

Letters 85 (19%)

113-126

Table 4 Comparison of the simulations of arterial blood concentration of styrene in rats exposed for 6 h to 80 ppm obtained using the present methodology with corresponding experimental data and simulations obtained with ACSLo Time (h)

Arterial blood concentration (mg/l) Experimentala

Simulated Using Excel@b

I 3 5 6 6.25 7 8 10 12 14 16 18 20 24

0.633 0.709 0.7779 0.931 1.017 0.199 0.095 0.0417 0.0239 0.01689 0.0108 0.0037 0.002

0.70587 0.80193 0.84818 0.59750 0.25930 0.14413 0.09730 0.05370 0.03010 0.01709 0.00948 0.00532 0.00299 0.00094

Using ASCL@ IC

IId

0.70669 0.80276 0.84901 0.59524 0.26045 0.14458 0.09756 0.05393 0.03012 0.01691 0.00954 0.00546 0.00302 0.09064

0.70697 0.80283 0.84907 0.59816 0.25960 0.14431 0.09743 0.05377 0.03014 0.00949 0.00533 0.00299 0.00094

‘The experimental data correspond to those used by Ramsey and Andersen [l] for validation of the rat styrene model. these simulation results were obtained using the methodology described in this paper. CThesesimulations were obtained using Gear algorithm in ASCLe, and correspond to those presented previously by Ramsey and Andersen [I]. dThese simulations were obtained using Euler algorithm in ASCLe.

to solve the set of equations constituting PBPK models [l-6]. The level of complexity of the majority of PBPK models, at least those developed so far, is relatively minor compared to common engineering problems for which the use of simulation software is absolutely needed to be practical. The present study demonstrates the feasibility of solving PBPK models without the use of programming/simulation languages, by simply using the computer as an efficient ‘calculator’. The use of Excel’s’or any other spreadsheet program in the present context is sought only to facilitate the task of providing output-input linkages between the various model equations. The method of solving PBPK models presented in this paper, to begin with, was attempted without the use of computers. The present approach would permit the user not only to ‘visualize’ the mechanics of simulation but also ‘see’ actually what goes on when PBPK model simulation is

conducted. A fundamental difference between the present method and the commercially available simulation software is that neither translation of the model equations into a programming language (e.g., FORTRAN) nor their compilation into an executable (.EXE) program are required in the present method. Erroneous simulations in the present approach can result only if wrong equations and parameter values are used. The simulations obtained with the present method should be comparable to those obtained with any commercial software as long as model equations are solved in an identical manner, as demonstrated by comparison with ACSLa in the present study. The commercial softwares such as ACSL@ provide a long list of algorithms from which the user can choose the appropriate one. In the present methodology, however, the user has to know the working mechanics of the algorithmic method that he/she wants to implement. We have, in addition

S. Ha&d et al. / Toxicology Letters 85 (195%) 113-126

123

Table 5 Comparison of the simulations of arterial blood concentration of styrene in rats exposed for 6 h to 600 ppm obtained using the present methodology with corresponding experimental dam and simulations obtained with ACSLe Time (h)

Arterial blood concentration (mg!l) Experimentala

Simulated Using Excel@b

1 3 5 6 6.25 7 8 10 12 14 16 18 20 24

9.1 16.9 .20.I .25.3 16.7 19.9 13.2 4.73 1.77 0.715 10.248 0.121 0.021

11.9218 19.2791 24.9199 25.4847 18.6018 13.0011 9.36880 4.91468 2.09283 0.99757 0.53516 0.29537 0.16454 0.05155

Using ASCLe IC

IP

11.9505 19.3135 24.9591 25.5217 18.6630 13.0398 9.39983 4.93728 2.10467 1.00235 0.53787 0.29703 0. I6568 0.05195

11.9601 19.3179 24.9630 25.5285 18.6398 13.0343 9.39860 4.93849 2.10466 1.00186 0.53728 0.29651 0.16517 0.05175

‘The experimental data correspond to those used by Ramsey and Andersen [l] for validation of the rat styrene model. these simulation results were obtained using the methodology described in this paper. %ese simulations were obtained using Gear algorithm in ACSLa, and correspond to those presented previously by Ramsey and Andersen [I]. dThese simulations were obtained using Euler algorithm in ACSLe.

to Euler, implemented fourth order Runge-Kutta procedure for solving differential equations representing PBPK models in Excel@ spreadsheets (not shown). The Euler method of integration used in the present study has previously been used by Johanson and Naslund [5] to conduct PBPK model simulations. However, these authors did not address the adequacy/approlpriateness of the following aspects relating to the ‘useof Euler method for PBPK simulations: the integration step size (INTSS), order of integration. method used, and ability to deal with stiffness in the model. In the present study, the INTSS warn fixed at 0.005 h. This INTSS is 3 times lower than the smallest time constant in the model (Table 6)‘. In simulation modeling, the standard practice is to use an INTSS smaller than the smallest time constant. Whereas Gear algorithm (in ACSL@) used an INTSS slightly lower than the value of the smallest time constant (-0.01 vs. 0.015 h)!, we chose to use one-third of

the smallest time constant. As long as the INTSS is smaller than or comparable to the smallest time constant of the model, it does not make any difference as to the stability of the simulation output. The only impact is that using very small INTSS would mean wastage of computer time. In the present context, with such a small PBPK model for styrene, the computer time is not a matter of concern. Even though the INTSS used in the Excels’ method (performing Euler integration) is somewhat lower than that used when Gear algorithm was specified in ACSLo, is adequate since it does not exceed the smallest time constant in the model. The PBPK model equations generally, and specifically in the example presented in this paper, are all first order differential equations. The solution to these first order equations can be obtained using integration methods of various orders. The order of an integration method is determined by the highest degree of the polynomial representing

124

S. Ha&d

et al. /Toxicology

Table 6 Time constants for the tissue compartments of the rat styrene PBPK model Compartments

Time constants (h)a

Liver Fat Richly perfused tissues Slowly perfused tissues

0.015 2.103 0.041 0.322

VMxlated Table 1.

as pP/Q] using the parameter values shown in

the solution. It is generally believed that, for a number of modeling applications, the higher the order of the integration method, the more accurate the simulation results, the relationship being valid for up to 5 orders [lo]. Therefore, the first order Euler method used in the present study may be viewed with reservation; however, the results indicate that the difference between the simulations obtained using Gear and Euler algorithms even within ACSLs is minimal, for the conditions and scenarios examined in the present study. In other words, there does not seem to be a loss in precision for not using the more complicated Gear integration method in the present study. The error associated with the use of Euler method arises from the negligence of second and other higher order derivatives of the Taylor expansion series; the error associated with Euler integration method then is proportional to At2. Subsequently, with the use of smaller At values, the error associated with this first order integration method and, therefore, the difference between the Euler and higher/variable order integration methods (e.g., Gear), may only be nominal as shown in the present study. The ability to deal with stiff systems is a characteristic of the Gear algorithm, most commonly used in previous PBPK modeling efforts. Even though Gear algorithm has frequently been used in the past, there is at least one example indicating probable inadequacy of this algorithm for PBPK modeling of specific substances [ 111. The stiffness in models is reflected by the ratio (generally of several orders of magnitude) of the largest to the smallest time constant in a model (R). For the styrene PBPK model used in the pres-

Letters 85 (19%)

113-126

ent study, this ratio (obtained with the consideration of only the partition coefficients, blood flow rates and tissue volumes) corresponds to 183.3 (Table 6), and therefore the ability of the algorithmic method to account for stiffness is important, if not crucial. The stiffness requires the use of variable numerical values for INTSS, since the R value changes as the simulation time progresses due to compartments with the smallest time constants attaining steady-state condition gradually. In the present study, adjustment of INTSS during simulation was not done; however, the INTSS used was smaller than the smallest time constant in the model (0.005 vs. 0.015 h). The fact that the INTSS value was not changed with the progression of simulation time, in accordance with the existing stiffness state, represents only a disadvantage of taking longer time to complete a given simulation. Particularly because at steady-state conditions (at which dAi/dt = 0) larger INTSS values can be used, thus saving some time without losing the accuracy of simulation (e.g., Ref. [5]). The use of consistently smaller INTSS throughout the simulation period including the steady-state condition (which was done in the present study to keep the example simple) results in longer simulation time. The present work should not be implied as recommending Euler method over the Gear’s stiff algorithm. The present study only attempted to develop a methodology to facilitate informed individuals understand and solve PBPK models without the use of ‘black-box’ kind of simulation and programming languages. We chose to present the implementation of the present methodology using Euler algorithm only due to the simplicity and not necessarily to its adequacy. Whereas beginners can start by using this algorithm, experienced individuals can use other algorithms (e.g., Runge-Kutta second order, Gear’s) which will then permit them to ‘see’ actually what goes on when simulations are conducted using these algorithms. These aspects might look trivial to people who are knowledgeable about simulation software and equally experienced well in the area of pharmacokinetic modeling. This certainly is not the case with the large majority of toxicologists, but they are increasingly getting into this area

S. Hadahd et al. /Toxicology Letters 85 (19%) 113-126

either because of a personal interest or due to corporate needs. This is a welcoming situation, as the PBPK models appear to have useful applications, particularly for reducing some uncertainties in the risk assessment process, and for designing specific, hypothesis-driven mechanistic toxicology studies. However, the ‘black-box’ kind of commercial simulation packages should not be used by the individuals without proper training in relevant computational techniques, since such practices could result in ‘expensive’ errors relating to the use of PBPK models in health risk assessment. The methodology developed in the present study provides a feasible alternative, and should enable any lay person to examine the existing models and develop newer ones, all without having to acquire knowledge of simulation languages. The present methodology is also a valuable tool for training/teaching purposes. Institutions that provide training in PBPK modeling have conventionally devoted several hours of instruction just to help the beginners get some hands-on experience with the use of particular simulation languages [6]. Logicall!y, it would be preferable and advantageous to provide such basic training in PBPK modeling without the use of simulation languages, using a methodology similar to that presented in this article. Once people understand how the system works using the present methodology, they can then move onto using advanced techniques and specialized simulation languages offering greater flexibility, speed, and additional features (e.g., sensitivity analysis, optimization routines). The major limitations of the proposed approach relate to (1) the number of cells to be selected in the spreadsheet, and (2) the runtime required to solve more complex PBPK models. Since pharmacokinetic data generated for each INTSS during the simulation exercise are displayed in one row of cells, the number of rows to be selected depends on the INTSS and the total time for simulation. Essentially, the number of rows available in a spreadsheet is finite, and with increasing number of rows the data handling will be tedious, rendering the runtime to be very long, resulting in questionable benefit of the present methodology. During the process of data analysis, if the compu-

125

ter does not have sufficient speed and runtime memory, it will likely ‘jam’. The occurrence of such a phenomenon does not indicate the limit of the proposed methodology, but rather that of the computer used by the investigator. Ideally, for complex modeling scenarios it would be preferable to use the present methodology to ‘visualize’ the working of selected sub-sections of the model. Copies of worksheets created using Excel@’for a 4-compartmental PBPK model presented in this paper (for Macintosh or IBM-PC) are available by writing to the authors. Acknowledgements This study represents an initiative undertaken as a part of Dr. Krishnan’s research program on physiological modeling supported by grants from the Ethyl Corporation, Canadian Network of Toxicology Centres (CNTC), Fonds de la recherche en Sante du Qutbec (FRSQ), Fonds pour la formation de chercheurs et l’aide a la recherche (FCAR), and the Natural Sciences and Engineering Research Council of Canada (NSERC). References 111Ramsey,

J.C. and Andersen, M.E. (1984) A physiologically-based description of the inhalation pharmacokinetics of styrene in rats and humans. Toxicol. Appl. Phannacol. 73, 159-175. 121Kootsey, J.M., Kohn, M.C., Feezor, M.D., Mitchell, G.R. and Fletcher, P.R. (1986) SCoP: an interactive simulation control program for micro and minicomputers. Bull. Math. Biol. 48, 427-441. [31 Menxel, D.B., Wolpert, R.L., Boger, J.R. and Kootsey, J.M. (1987) Resources available for simulation in toxicology: specialized computers, generalized software and communication networks. Drinking Water Health 8, 229-254. [41 Hattis, D., White, P., Marmorstein, L. and Koch, P. (1990) Uncertainties in pharmacokinetic modeling for perchloroethylene. I. Comparison of model structure, parameters, and predictions for low dose metabolic rates for models by different authors. Risk Anal. 10,449-458. 151Johanson, G. and Naslund, P.H. (1988) Spreadsheet programming: a new approach in physiologically-based modeling of solvent toxicokinetics. Toxicol. Lett. 41, 115-127. 161Dong, M.H. (1994) Microcomputer programs for physiologically-based pharmacokinetic (PBPK) modeling. Comp. Methods Prog. Biomed. 45, 213-221.

126

S. HaaSd

et al. /Toxicology

[7] Davis, P.W. (1992) Differential Equations for Mathemat-

its, Science, and Engineering. Prentice-Hall, Englewood Cliffs, N.J. [8] Mitchell and Gauthier Associates Inc. (1993) Advanced Continuous Simulation Language (ACSL): Reference Manual. Edition 10.1, Mitchell and Gauthier, Concord, Mass. [9] Rideout, V.C. (1991) Mathematical and Computer Modeling of Physiological Systems. Prentice-Hall, N.Y.

Letters 8S (19%)

113-126

[lo] van der Bosch, P.P.J. and van der Rlaun, A.C. (1994) Modeling: Identification and Simulation of Dynamic Systems. CRC Press Inc., Boca Raton, Fla. [ll] Wada, R.D., Stanski, D.R. and Ebling, W.F. (1995) A PC-based graphical simulator for physiological pharmacokinetic models. Comp. .Methods Prog. Biomed. 46, 245-255.