www.biomathforum.org/biomath/index.php/biomath ORIGINAL ARTICLE A computational investigation of the ventilation structure and maximum rate of metabolism for a physiologically based pharmacokinetic (PBPK) model of inhaled xylene Karen A. Yokley∗‡, Jaclyn Ashcraft∗, and Nicholas S. Luke† ∗Department of Mathematics and Statistics Elon University, Elon, NC 27244 †Department of Mathematics North Carolina A&T State University, Greensboro, NC 27411 ‡ Corresponding Author: kyokley@elon.edu Received: 29 August 2018, accepted: 6 January 2019, published: 4 February 2019 Abstract—Physiologically based pharmacokinetic (PBPK) models are systems of ordinary differential equations that estimate internal doses following exposure to toxicants. Most PBPK models use stan- dard equations to describe inhalation and concen- trations in blood. This study extends previous work investigating the effect of the structure of air and blood concentration equations on PBPK predictions. The current study uses an existing PBPK model of xylene to investigate if different values for the maximum rate of toxicant metabolism, V xylmax, can result in similar compartmental predictions when used with different equations describing inhalation. Simulations are performed using V xylmax values based on existing literature. Simulated data is also used to determine specific V xylmax values that result in similar predictions from different ventilation structures. Differences in ventilation equation structure may affect parameter estimates found through inverse problems, although further investigation is needed with more complicated models. Keywords-PBPK modeling, xylene I. INTRODUCTION Physiologically based pharmacokinetic (PBPK) modeling uses ordinary differential equations to describe absorption, distribution, metabolism, and excretion of toxicants following exposure. PBPK models have been developed to estimate internal doses for various toxicants in rodents [1] [9] [11] [14] [17] [18] [25] [34] [37] [40] and in humans [2] [3] [8] [10] [13] [15] [23] [26] [35] [36] [39] [38] [43] [48]. Most PBPK model equations use Copyright: c©2019 Yokley et al. This article is distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Citation: Karen A. Yokley, Jaclyn Ashcraft, Nicholas S. Luke, A computational investigation of the ventilation structure and maximum rate of metabolism for a physiologically based pharmacokinetic (PBPK) model of inhaled xylene, Biomath 8 (2019), 1901067, http://dx.doi.org/10.11145/j.biomath.2019.01.067 Page 1 of 13 http://www.biomathforum.org/biomath/index.php/biomath https://creativecommons.org/licenses/by/4.0/ https://creativecommons.org/licenses/by/4.0/ http://dx.doi.org/10.11145/j.biomath.2019.01.067 K.A. Yokley, J. Ashcraft, N.S. Luke, A computational investigation of the ventilation structure and ... the same basic structure and assumptions, such as the assumption that compartments are well-mixed and that the transfer of some chemicals is at equi- librium. However, the appropriate use of PBPK model estimates in risk assessment is dependent on multiple factors, including the biological basis of the model parameters and structure [6]. The industrial solvent xylene is a component of paints, paint thinners, and related products [42]. Workers in specific industries may be at higher risk for xylene exposure and poorly ventilated areas may amplify exposure dangers [16], and xylene may increase the effect of other chemicals when present in mixtures [4]. As in [46], the PBPK model of xylene from [39] is used for the current investigation because of the relative simplicity of the model. The purpose of this project is to expand on the previous PBPK investigation that considered different modeling approaches for ventilation [46] and to also consider variation or errors in the parameter for the maximum rate of metabolism, Vmax (mg/h). Yokley [46] considered three struc- tures for modeling ventilation, described as “equi- librium,” “non-equilibrium,” and “hybrid” models. The results were very similar for the equilibrium and hybrid models, and therefore only two struc- tures will be used to model ventilation in the current study. The equilibrium model uses the standard quotient for the concentration of toxi- cant in arterial blood, and the non-equilibrium model allows for a separate lung compartment. The equilibrium and non-equilibrium models are used to predict various xylene concentrations and amounts in the body following inhalation exposure using different values for the maximum rate of metabolism in the liver. Results of simulations with different metabolic rates are then used in inverse problems to investigate the success of optimization using different types of data. II. MODEL INVESTIGATION A. Model Background In [46], a PBPK model for xylene exposure in humans from [39] was used with different Fig. 1. The structure of the PBPK model of xylene from Tardif et al. [39] used in [46]. equations describing ventilation exposure in order to determine the effect of the structure of those equations on model output. The four compart- ment model described xylene concentration in the slowly perfused tissues, the rapidly perfused tissues, the fat, and the liver. A schematic of the overall model is presented in Figure 1. For easy reference, a summary of the notations used throughout the model is presented in the Ap- pendix. The amount of xylene within each compartment depicted in Figure 1 is modeled by a separate dif- ferential equation. Additional equations are added to the model to incorporate the ventilation struc- ture. For each model equation, the compartmental concentration is defined as C xyl j = A xyl j Vj , (1) where Cxylj represents the concentration of the Biomath 8 (2019), 1901067, http://dx.doi.org/10.11145/j.biomath.2019.01.067 Page 2 of 13 http://dx.doi.org/10.11145/j.biomath.2019.01.067 K.A. Yokley, J. Ashcraft, N.S. Luke, A computational investigation of the ventilation structure and ... chemical, xylene, in compartment j (mg/L); Axylj represents the amount of the chemical, xylene, in compartment j (mg); and Vj represents the volume of compartment j (L). Compartment j is either the slowly perfused tissue (denoted s), the rapidly perfused tissue (denoted r), the adipose tissue or fat (denoted f), or the liver (denoted l). The differential equation for the change in the amount of xylene within each internal compart- ment j (excluding the liver) takes the form dA xyl j dt = Qj ( C xyl art − C xyl j P xyl j ) , (2) where Qj represents the rate of blood flow to compartment j (L/h); Cxylart is the concentration of the chemcial, xylene, in the arterial blood (mg/L); and P xylj is the partition coefficient for xylene in tissue j (unitless). The differential equation representing the change in the amount of xylene within the liver takes a different form because xylene is metab- olized in the liver. The structure of the equation is similar to that of the other compartments, with a basic flow in minus flow out, but an additional term is subtracted to represent the metabolism. The form of the equation for the liver is dA xyl l dt =Ql ( C xyl art − C xyl l P xyl l ) − V xyl max C xyl l P xyl l K xyl m + C xyl l P xyl l , (3) where V xylmax is the maximum rate of metabolism for xylene (mg/h) and KxylM is the concentration of xylene at half saturation (mg/L). The full PBPK model is comprised of three equations of form (2) (one each for the slowly perfused, rapidly perfused, and fat compartments), equation (3), and two or more equations that repre- sent the ventilation structures, which are outlined below. The three ventilation structures used in [46] were classified as “equilibrium,” “non- equilibrium,” and “hybrid” cases. Many PBPK models use the “equilibrium” case equations for blood concentrations [14] [23] [39], which are constructed under a few assumptions including that the amount of toxicant leaving the alveolar space is equal to the amount eliminated through the blood. Equations similar to (4) and (5) below are typically used in PBPK modeling for concentrations of a particular toxicant, i, in the venous (ven) and arterial (art) blood: Ciart = QcC i ven + QpC i inh Qc + Qp P ibl:air (4) Civen = ∑ j Qj Aij P ij Vj Qc , (5) where Qc represents the cardiac flow rate, defined as the sum of the compartmental flow rates (L/h); Qp is the pulmonary flow or alveolar ventilation (L/h); and Pbl:air is the blood/air partition coeffi- cient (unitless). The “non-equilibrium” case allows for the amount of toxicant leaving the alveolar space (alv) to not equal the amount entering the alveolar space and thus the alveolar space, arterial concentration, and venous concentration are each represented with their own compartment in the model. This scenario can be modeled using the equations dAialv dt = Qp ( Ciinh − Aialv ValvP i bl:air ) +Qc ( Civen − Aialv ValvP i alv:air ) (6) dCiart dt = Qc ( Aialv ValvP i alv − Ciart ) (7) dCiven dt = 1 Vven  ∑ j Qj Aij P ij Vj − QcCiven  (8) which are described more fully in [46]. The “hy- brid” case involves using the “non-equilibrium” equations (6)-(8) with the alveolar partition coef- ficient, P ialv, set to be 1. Because the results from [46] were very similar between the “equilibrium” and “hybrid” cases, the “hybrid” scenario is omit- ted from the current investigation. B. Fixed Parameter Values All simulations use a body weight BW of 70 kg, under the assumption that 1 L = 1 kg Biomath 8 (2019), 1901067, http://dx.doi.org/10.11145/j.biomath.2019.01.067 Page 3 of 13 http://dx.doi.org/10.11145/j.biomath.2019.01.067 K.A. Yokley, J. Ashcraft, N.S. Luke, A computational investigation of the ventilation structure and ... TABLE I THE COMPARTMENTAL FLOWS AND PARTITION COEFFICIENTS USED IN SIMULATIONS OF THE XYLENE PBPK MODEL. Compartmental flow Value (% Qc) Source Fat 5 [39] Slowly perfused 25 [39] Richly perfused 44 [39] Liver 26 [39] Partition coefficient Value Source Blood:air 26.4 [39] Fat:blood 77.8 [39] Slowly perfused:blood 3.0 [39] Richly perfused:blood 4.42 [39] Liver:blood 3.02 [39] Lung:air 87.4 [42] is a sufficient conversion. Physiological parame- ters used in the PBPK model are presented in Tables I and II. As described in [46], the lung:air partition coefficient (P xylalv:air) was obtained from [42] and was used with the blood:air partition coefficient from [39] in order to find P xylalv (i.e., P xyl alv = P xyl alv:air/P xyl bl:air). Metabolic values from the PBPK source model are used as follows, V xylmax = V xyl maxcBW 0.75 Kxylm = 0.2, and cardiac output and alveolar ventilation are assumed to follow allometric scaling with scaling factors used from the source model Qc = 18.0 · BW 0.70 = Qp [39] . C. Vmaxc values Metabolism of xylene is described using Michaelis-Menten kinetics in the model [39] which involves a nonlinear term with the param- eter V xylmax representing the maximum rate of the enzymatic process (see Equation (3)). The value of V xylmax is scaled to body weight as follows: V xylmax = V xyl maxcBW 0.75 . (9) Different values of V xylmaxc have been used in PBPK modeling and other research as is shown in Ta- ble III. Parameter values have been rescaled for TABLE II THE COMPARTMENTAL VOLUMES USED IN SIMULATIONS OF THE XYLENE PBPK MODEL. ALL VALUES ARE THE SAME AS USED FOR THE “NON-AGED ADULT” IN [47] EXCEPT FOR Vr , Vbl, Valv , Vart, AND Vven. Vr WAS ALTERED TO HAVE ALL VOLUMES SUM TO 100%. Vbl Valv , Vart, AND Vven ARE BASED ON [5]. Volume Value Vl 0.026BW Vf 0.214BW Vs 0.613BW Vr 0.06BW Vbl 0.079BW Valv 0.008BW Vart 0.2Vbl Vven 0.8Vbl consistency to represent V xylmaxc as in (9). Note that the PBPK source model for the work in [46] focuses on m-xylene (as do the majority of the references in Table III). The values from [41] are estimated through conversion using different values for liver size. The values in Table III are used as a basis for V xylmaxc values used in the simulations. V xylmaxc values used in simulations or found in optimized fits to simulated output fall within a range of 1-11 mg/(h·kg) which is similar to the range of values presented in Table III. D. Investigational Methods Computational solutions are generated using MATLAB R2015b [27] with the ode15s solver. With the exception of values for V xylmaxc, all initial conditions, exposure scenarios (50 ppm xylene over 7 hours), and parameters remain unchanged from the previous investigation [46]. Solution curves are generated for values of V xyl maxc ranging from V xyl maxc0 = 1 to V xyl maxcfinal = 10, which are chosen to illustrate the trend of output as well as include the value, 8.4, used in the original xylene model [39]. A list of V xylmaxci values is generated with spacing 0.1, and computational solutions are generated for each metabolic constant in this list. The exposure duration is 7 hours, and 20 hours are used for output calculations to allow for the majority of the toxicant to be cleared from the system. Maximum output values over the Biomath 8 (2019), 1901067, http://dx.doi.org/10.11145/j.biomath.2019.01.067 Page 4 of 13 http://dx.doi.org/10.11145/j.biomath.2019.01.067 K.A. Yokley, J. Ashcraft, N.S. Luke, A computational investigation of the ventilation structure and ... TABLE III TABLE OF VARIOUS VALUES OF V xylmaxc VALUES AND REFERENCES CITING THEM. NUMERICAL VALUES HAVE BEEN RESCALED FOR CONSISTENCY. LISTED SOURCES EITHER USED THE VALUE GIVEN (NOT NECESSARILY SCALED) OR USED THE DATA FROM THE PAPER IN SOME WAY. NOTE THAT THE REFERENCE MARKED (*) CITED THE GIVEN SOURCE BUT USED 6.49. Primary Source Scaled V xylmaxc Cited By Tardif et al. 1993 [40] and 8.4 mg/(h·kg) [12], [19], [21], [22], [24], Tardif et al. 1995 [39] [31], [44], [45], [46] Tardif et al. 1997 [38] 5.5 mg/(h·kg) [26], [18]* Tassaneeyakul et al. 1996 [41] ≈11.3-13.7 mg/(h·kg) [29], [30] Price and Krishnan 2011 [33] Measured: 6.475 mg/(h·kg), [7], [32] Predicted: 5.300 mg/(h·kg) Haddad et al. 1999 [20] 6.59 mg/(h·kg) [22] 20 hour simulations are found for the amount of xylene in exhaled air, concentration of xylene in the venous blood, and amount of xylene in the liver. Output over the 20 hours is generated at each 0.05 step, and a Riemann sum is used to estimate area below the curve for the xylene concentration in the liver. Output in amount is converted to con- centration before the area is estimated. The xylene concentration in the venous blood and amount in the exhaled air are investigated because both could be physically collected, and those data could be used to estimate V xylmaxc values. The predictions for xylene in the liver are investigated because the liver is often a target organ (in general for various toxicants) and for the use of liver predictions in risk assessment. In order to ascertain how different V xylmaxc values could produce similar results with the two differ- ent models, the following procedure is employed. Simulated data is generated from each model using a beginning V xylmaxc value for one model output that conceivably could be compared to measured data. Different levels of noise are added to the data, and then both models are optimized to the simulated data over V xylmaxc, minimizing a least squares cost function J(Vmaxc) = min ∑ i ( xi(V xyl maxc) 2 − d2i ) where xi represents the state of the differential equation model (using a particular value of V xylmaxc) at time ti corresponding to data point di. Optimiza- tion is performed using the fminsearchbnd func- tion [28] from the Matlab Optimization Toolbox. Data are simulated for the concentration of xylene in exhaled air and in the venous blood for both the equilibrium and non-equilibrium models and then fit to output of each model. III. RESULTS The amount (or concentration) of a chemical within any given compartment during exposure is characterized by an increase during the interval of exposure followed by a gradual decline or clearing of the chemical after exposure has ceased. Between the interval of exposure and the clearance of the chemical, the amount of chemical in any given compartment will reach a maximum level. In order to fit chemical exposure data, it is imperative that the PBPK model under consideration can reach these maximum levels. As a premilinary ex- ercise for this study, the maximum model outputs are generated for both the equilibrium and non- equilibrium models and compared to determine if there is any overlap. Figure 2 depicts the relationship between max- imum predicted xylene amount in the exhaled air and values of V xylmaxc for both the equlibrium and non-equlibrium models. It is evident from this graph that there is little to no overlap in the max- imum predicted xylene amount for both models. Biomath 8 (2019), 1901067, http://dx.doi.org/10.11145/j.biomath.2019.01.067 Page 5 of 13 http://dx.doi.org/10.11145/j.biomath.2019.01.067 K.A. Yokley, J. Ashcraft, N.S. Luke, A computational investigation of the ventilation structure and ... 1 2 3 4 5 6 7 8 9 10 Scaling constant for Vmax 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 M a xi m u m p re d ic te d a m o u n t in e xh a le d a ir Equilibrium Non-equilibrium Fig. 2. Maximum model output for the amount of xylene in exhaled air for various V xylmaxc values. 1 2 3 4 5 6 7 8 9 10 Scaling constant for Vmax 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 M a xi m u m p re d ic te d c o n ce n tr a tio n in v e n o u s b lo o d Equilibrium Non-equilibrium Fig. 3. Maximum model output for the concentration of xylene in the venous blood for various V xylmaxc values. The maximum amounts generated from the non- equilibrium model are consistently higher than those of the equilibrium model. The significant difference between the model outputs would sug- gest that the non-equilibrium model would not be capable of fitting data generated by the equilibrium model, and the equilibrium model would not be capable of fitting data that was generated by the non-equilibrium model. Figures 3 and 4 illustrate the relationship be- tween values of V xylmaxc with the maximum concen- tration of xylene in the venous blood and maxi- 1 2 3 4 5 6 7 8 9 10 Scaling constant for Vmax 0 1 2 3 4 5 6 7 M a xi m u m p re d ic te d a m o u n t in li ve r Equilibrium Non-equilibrium Fig. 4. Maximum model output for the amount of xylene in the liver for various V xylmaxc values. 1 2 3 4 5 6 7 8 9 10 Scaling constant for Vmax 0 0.2 0.4 0.6 0.8 1 1.2 E st im a te d a re a b e lo w p re d ic te d li ve r co n c. c u rv e Equilibrium Non-equilibrium Fig. 5. Estimated area below the predicted liver concentration curve for various V xylmaxc values. mum amount of xylene in the liver, respectively. In contrast to the maximum xylene amount in the exhaled air (Figure 2), maximum values of xylene in the venous blood and liver show a considerable amount of overlap between the equilibrium and non-equilibrium models. This indicates that the same maximum level may be predicted using ei- ther the equilibrium model or the non-equilibrium model with different values for V xylmaxc. A similar relationship for the area below the curve estimates for concentation of xylene in the liver is presented in Figure 5. Biomath 8 (2019), 1901067, http://dx.doi.org/10.11145/j.biomath.2019.01.067 Page 6 of 13 http://dx.doi.org/10.11145/j.biomath.2019.01.067 K.A. Yokley, J. Ashcraft, N.S. Luke, A computational investigation of the ventilation structure and ... 0 2 4 6 8 10 12 14 16 18 20 Time (h) 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 C o n ce n tr a tio n o f xy le n e in e xh a le d a ir ( m g /L ) Equilibrium output using simulated equilibrium data (a) Equilibrium Output 0 2 4 6 8 10 12 14 16 18 20 Time (h) 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 C o n ce n tr a tio n o f xy le n e in e xh a le d a ir ( m g /L ) Non-equilibrium output using simulated equilibrium data (b) Non-Equilibrium Output Fig. 6. Model output for the concentration of xylene in exhaled air after optimization over data simulated from the equilibrium model with noise level 0.001 and V xylmaxc = 6. Figure (a) contains equilibrium model output and Figure (b) contains non- equilibrium model output. 0 2 4 6 8 10 12 14 16 18 20 Time (h) 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 C o n ce n tr a tio n o f xy le n e in e xh a le d a ir ( m g /L ) Equilibrium output using simulated non-equilibrium data (a) Equilibrium Output 0 2 4 6 8 10 12 14 16 18 20 Time (h) 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 C o n ce n tr a tio n o f xy le n e in e xh a le d a ir ( m g /L ) Non-equilibrium output using simulated non-equilibrium data (b) Non-Equilibrium Output Fig. 7. Model output for the concentration of xylene in exhaled air after optimization over data simulated from the non- equilibrium model with noise level 0.001 and V xylmaxc = 8. Figure (a) contains equilibrium model output and Figure (b) contains non-equilibrium model output. Based on these comparisons of the maximum value of xylene for each variation of the model, it seems that the equilibrium and non-equilibrium models would fit the opposite data set more ef- ficiently for the concentration of xylene in the venous blood and the amount of xylene in the liver than they would for the amount of xylene in the exhaled air. While the liver is a target organ for risk assessment, experimental data for the amount of a chemical within the liver is not easily collected; thus it is often unavailable. For these reasons, more focus is placed on predictions of blood concentrations for this study. Examples of graphical output for xylene in ex- haled air using simulated data from the equilibrium model are contained in Figure 6, and examples us- ing simulated data from the non-equilibrium model are contained in Figure 7. These figures provide examples of best fitting curves to the simulated data. A summary of simulation results for exhaled Biomath 8 (2019), 1901067, http://dx.doi.org/10.11145/j.biomath.2019.01.067 Page 7 of 13 http://dx.doi.org/10.11145/j.biomath.2019.01.067 K.A. Yokley, J. Ashcraft, N.S. Luke, A computational investigation of the ventilation structure and ... TABLE IV OPTIMAL COST USING SIMULATED DATA FOR THE CONCENTRATION OF XYLENE IN EXHALED AIR. THE V xylmaxc VALUE USED FOR SIMULATION IS LISTED IN COLUMN 1 OF THE TABLE, AND INITIAL GUESSES FOR V xylmaxc FOR THE OPTIMIZATION ROUTINE WERE 3, 6, AND 10. THE VALUE FOR J∗ THAT WAS LOWEST FOR THE THREE IS LISTED BELOW. IN CASES MARKED WITH (†) THE SAME COST WAS FOUND FOR VARIOUS V xylmaxc VALUES, ALL OF WHICH WERE CLOSE TO ZERO. Equilibrium simulated data V xylmaxc Noise Best Fit to “Equilibrium” Best Fit to “Non-equilibrium” 2 0.001 V xylmaxc = 1.9125, J ∗ = 1.6885e−6 V xylmaxc = 14.9393, J ∗ = 0.0224 4 0.001 V xylmaxc = 3.7144, J ∗ = 1.8447e−6 V xylmaxc = 14.9386, J ∗ = 0.0281 6 0.001 V xylmaxc = 5.1232, J ∗ = 2.1612e−6 V xylmaxc = 14.9383, J ∗ = 0.0299 Non-equilibrium simulated data V xylmaxc Noise Fit to “Equilibrium” Fit to “Non-equilibrium” 4 0.001 V xylmaxc = 0.0014 †, J∗ = 0.0151 V xylmaxc = 3.5138, J ∗ = 4.9275e−5 6 0.001 V xylmaxc = 0.0013 †, J∗ = 0.0123 V xylmaxc = 5.7957, J ∗ = 4.3432e−6 8 0.001 V xylmaxc = 0.0051 †, J∗ = 0.0114 V xylmaxc = 7.0188, J ∗ = 1.5291e−6 0 2 4 6 8 10 12 14 16 18 20 Time (h) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 C o n ce n tr a tio n o f xy le n e in v e n o u s b lo o d ( m g /L ) Equilibrium output using simulated equilibrium data (a) Equilibrium Output 0 2 4 6 8 10 12 14 16 18 20 Time (h) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 C o n ce n tr a tio n o f xy le n e in v e n o u s b lo o d ( m g /L ) Non-equilibrium output using simulated equilibrium data (b) Non-Equilibrium Output Fig. 8. Model output for the concentration of xylene in venous blood after optimization over data simulated from the equilibrium model with noise level 0.05 and V xylmaxc = 2. Figure (a) contains equilibrium model output and Figure (b) contains non-equilibrium model output. air concentration is presented in Table IV. In Figure 6, it can be observed that the equi- librium model (depicted in Figure 6(a)) provides a much better fit to the equilibrium data than the non-equilibrium model (depicted in Figure 6(b)) does. When fit to the equilibrium data, the equi- librium model produces a least squares cost of 2.1612e−6 and the non-equilibrium model pro- duces a cost of 0.0299. Similary, Figure 7 shows that the non-equilibrium model provides a much better fit to the non-equilibrium data than that of the equilibrium model. The non-equilibrium model yields a least squares cost of 1.5291e−6 with the non-equilbrium data, compared to the equilibrium model’s cost of 0.0114. Data are also simulated for the concentration of xylene in venous blood for both the equilibrium and non-equilibrium models. The fits to these models show some success for V xylmaxc in the range of 3-7. Examples of graphical results are contained Biomath 8 (2019), 1901067, http://dx.doi.org/10.11145/j.biomath.2019.01.067 Page 8 of 13 http://dx.doi.org/10.11145/j.biomath.2019.01.067 K.A. Yokley, J. Ashcraft, N.S. Luke, A computational investigation of the ventilation structure and ... 0 2 4 6 8 10 12 14 16 18 20 Time (h) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 C o n ce n tr a tio n o f xy le n e in v e n o u s b lo o d ( m g /L ) Equilibrium output using simulated non-equilibrium data (a) Equilibrium Output, noise=0.02 0 2 4 6 8 10 12 14 16 18 20 Time (h) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 C o n ce n tr a tio n o f xy le n e in t h e v e n o u s b lo o d ( m g /L ) Equilibrium output using simulated equilibrium data (b) Equilibrium Output, noise=0.05 0 2 4 6 8 10 12 14 16 18 20 Time (h) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 C o n ce n tr a tio n o f xy le n e in v e n o u s b lo o d ( m g /L ) Non-equilibrium output using simulated non-equilibrium data (c) Non-Equilibrium Output, noise=0.02 0 2 4 6 8 10 12 14 16 18 20 Time (h) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 C o n ce n tr a tio n o f xy le n e in v e n o u s b lo o d ( m g /L ) Non-equilibrium output using simulated non-equilibrium data (d) Non-Equilibrium Output, noise=0.05 Fig. 9. Model output for the concentration of xylene in venous blood after optimization over data simulated from the non-equilibrium model with given noise level and V xylmaxc = 6. Figures (a)-(b) contain equilibrium model output and Figures (c)-(d) contain non-equilibrium model output. in Figure 8 and Figure 9, and a summary of overall results for the concentration of xylene in venous blood is presented in Table V. Figure 8 illustrates the best fits of the equi- librium and non-equilibrium models to venous blood concentration data that was generated by the equilibrium model with a V xylmaxc value of 2. Unlike the previously presented results for the model fits to exhaled air, both models seem to provide an adequate fit to the venous blood data. The results from the equilibrium model (in Figure 8(a)) seem to capture the maximum more efficiently. The best fit for the equilibrium model is produced using a V xyl maxc value of 1.8299, with a cost function value of 0.0058. The best fit for the non-equilibrium model is produced with a V xylmaxc value of 1.8514, yielding a cost of 0.0203. The fits of the equilibrium and non-equilibrium models to non-equilibrium data are presented in Figure 9. Figures 9(a) and 9(b) display the fit of the equilibrium model to non-equilibrium data with noise levels of 2% and 5%, respectively. Results for the non-equilibrium model fitted to the non- equilibrium data with noise levels of 2% and 5% are depicted in Figures 9(c) and 9(d). A visual inspection of the graphs would suggest that both the equilibrium model and the non-equilibrium model can provide an adequate fit to the simulated Biomath 8 (2019), 1901067, http://dx.doi.org/10.11145/j.biomath.2019.01.067 Page 9 of 13 http://dx.doi.org/10.11145/j.biomath.2019.01.067 K.A. Yokley, J. Ashcraft, N.S. Luke, A computational investigation of the ventilation structure and ... TABLE V OPTIMAL COST USING SIMULATED DATA FOR THE CONCENTRATION OF XYLENE IN VENOUS BLOOD. THE V xylmaxc VALUE USED FOR SIMULATION IS LISTED IN COLUMN 1 OF THE TABLE, AND INITIAL GUESSES FOR V xylmaxc FOR THE OPTIMIZATION ROUTINE WERE 3, 6, AND 10. THE VALUE FOR J∗ THAT WAS LOWEST FOR THE THREE IS LISTED BELOW. Equilibrium simulated data V xylmaxc Noise Best Fit to “Equilibrium” Best Fit to “Non-equilibrium” 2 0.001 V xylmaxc = 1.9687, J ∗ = 2.2307e−4 V xylmaxc = 2.0752, J ∗ = 0.0247 2 0.01 V xylmaxc = 1.9596, J ∗ = 4.7682e−4 V xylmaxc = 2.0700, J ∗ = 0.0234 2 0.02 V xylmaxc = 1.9173, J ∗ = 7.6803e−4 V xylmaxc = 2.0190, J ∗ = 0.0237 2 0.05 V xylmaxc = 1.8299, J ∗ = 0.0058 V xylmaxc = 1.8514, J ∗ = 0.0203 3 0.001 V xylmaxc = 3.0000, J ∗ = 7.1896e−6 V xylmaxc = 4.8376, J ∗ = 0.0139 3 0.01 V xylmaxc = 2.9006, J ∗ = 2.3456e−4 V xylmaxc = 4.5075, J ∗ = 0.0132 3.5 0.05 V xylmaxc = 3.1610, J ∗ = 0.0058 V xylmaxc = 5.4462, J ∗ = 0.0127 Non-equilibrium simulated data V xylmaxc Noise Fit to “Equilibrium” Fit to “Non-equilibrium” 6 0.001 V xylmaxc = 3.3016, J ∗ = 0.0128 V xylmaxc = 6.0000, J ∗ = 8.5450e−6 6 0.01 V xylmaxc = 3.3067, J ∗ = 0.0154 V xylmaxc = 6.2866, J ∗ = 7.8847e−4 6 0.02 V xylmaxc = 3.3163, J ∗ = 0.0201 V xylmaxc = 6.6882, J ∗ = 0.0050 6 0.05 V xylmaxc = 3.3067, J ∗ = 0.0489 V xylmaxc = 6.2866, J ∗ = 0.0466 8 0.001 V xylmaxc = 3.6121, J ∗ = 0.0111 V xylmaxc = 10.6876, J ∗ = 8.9750e−4 8 0.01 V xylmaxc = 3.5750, J ∗ = 0.0111 V xylmaxc = 10.2790, J ∗ = 0.0017 non-equilibrium data. The non-equilibrium data in these graphs were generated with a V xylmaxc value of 6. The best fit of the equilibrium model to the data with 2% noise is found with V xylmaxc equal to 3.3163 and has a least squares cost of 0.0201. A V xylmaxc value of 3.3067 leads to the optimal fit of the equilibrium model to the non-equilibrium data with 5% noise, resulting in a cost of 0.0489. The best fits for the non-equilibrium model to the non-equilibrium data with 2% and 5% noise are produced with V xylmaxc values of 6.6882 and 6.2866 and result in cost values of 0.0050 and 0.0466, respectively. As previously stated, the amount of xylene present in the liver following an inhalation ex- posure is not easily measured and therefore not reported. For this reason, a comparison of the equilibrium and non-equilibrium model results for the amount of xylene in the liver is not conducted. Based on the maximum model output for the amount of xylene in the liver (depicted in Fig- ure 4), it is hypothesized that adequate fits to data can be found using both the equilibrium and non- equilibrium models (similar to the results reported for the venous blood concentrations above). IV. DISCUSSION AND CONCLUSIONS Simulations suggest that different V xylmaxc values could be used to make similar predictions for the concentration of xylene in venous blood with the different ventilation structures, but different V xylmaxc values are not found to produce similar model output for the amount in exhaled air. Specifically, using V xylmaxc around 3 in the “equilibrium” model produces similar blood concentration results as using V xylmaxc around 6 in the “non-equilibrium” model. As shown in Figure 3 and in Figure 9, the maximum xylene concentration is predicted to be about 0.4 mg/L by both “equilibrium” and “non-equilibrium” models when V xylmaxc is around 3 or 6, respectively. Additionally, optimization to simulated data results in best fits with similar numbers (see results in Table V for simulated V xyl maxc=3, 3.6, and 6). Hence, blood data may not be as informative as exhaled air data for identifying the maximum rate of metabolism when Biomath 8 (2019), 1901067, http://dx.doi.org/10.11145/j.biomath.2019.01.067 Page 10 of 13 http://dx.doi.org/10.11145/j.biomath.2019.01.067 K.A. Yokley, J. Ashcraft, N.S. Luke, A computational investigation of the ventilation structure and ... the exposure method is inhalation. Results may differ if a toxicant is administered in the blood directly or ingested or exposed dermally as in [36]. In cases involving exposure through methods other than inhalation, the ventilation equation structure would be expected to be less critical and so are not investigated in the current study. However, in those cases, the equations describing the entrance of toxicant into the body could be critical as well. Figures 4 and 5 contain predictions related to xylene in the liver. When these predictions are cal- culated using the “equilibrium” model with V xylmaxc around 3 and the “non-equilibrium” model with V xyl maxc around 6, however, the liver predictions are higher for the “equilibrium” model. Although the “equilibrium” model for the xylene PBPK model used here may then provide higher predictions for internal liver doses, the “equilibrium” ventilation structure may not provide higher internal dose estimates when used with more complicated PBPK models. For example, the results in Figures 4 and 5 contain predictions for the parent chemical, and concentrations of metabolites may show different dynamics. Although the PBPK models of xylene used in the current study do not focus on the excretion of toxicants, some models do make predictions of the parent chemical or metabolites in the urine. Data of excreted toxicants could also be problem- atic for use in optimizing metabolic parameters. Additionally, the model used in the current study is a more simplistic PBPK model, and more inves- tigation is needed to be able to make conclusions about toxicants with harmful metabolites or that require models with more compartments. PBPK models have been used to describe exposure to a mixture of chemicals (such as in [12] [18] [24] [33][39] [43]) which would also involve a more complicated investigation than in the current study. V. ACKNOWLEDGEMENTS The authors would like to thank Elon Univer- sity Funding for Research and Development for support for this project. APPENDIX Abbreviations: xyl xylene s slowly perfused tissues f adipose tissue or fat r rapidly perfused tissues l liver ven venous blood art arterial blood bl blood (mixed) alv alveolar space or respiratory compartment Model Notations: Qj Blood flow in/to compartment j (L/h) Qc Cardiac flow, ∑ j Qj (L/h) Qp Pulmonary flow/alveolar ventilation (L/h) Vj Volume of compartment j (L) P ij The partition coefficient for chemical i in tissue j (dimensionless) P xyl bl:air The blood/air partition coefficient for chemical xyl (dimensionless) A xyl j The amount of chemical xyl in tissue j (mg) C xyl j The concentration of chemical xyl in compartment j (mg/L) C xyl inh The concentration of chemical xyl inhaled (mg/L) V xylmax The maximum rate of metabolism for chemical xyl (mg/h) Kxylm The concentration of xyl at half saturation (mg/L) REFERENCES [1] Abbas, R, Fisher, JW (1997) A Physiologically Based Pharmacokinetic Model for Trichloroethylene and Its Metabolites, Chloral Hydrate, Trichloroacetate, Dichloroacetate, Trichloroethanol, and Trichloroethanol Biomath 8 (2019), 1901067, http://dx.doi.org/10.11145/j.biomath.2019.01.067 Page 11 of 13 http://dx.doi.org/10.11145/j.biomath.2019.01.067 K.A. Yokley, J. Ashcraft, N.S. Luke, A computational investigation of the ventilation structure and ... Glucuronide in B6C3F1 Mice, Toxicol Appl Pharm, 147:15-30. [2] Allen, BC, Fisher, JW (1993) Pharmacokinetic Mod- eling of Trichloroethylene and Trichloroacetic Acid in Humans, Risk Anal, 13(1):71-86. [3] Blancato, JN, Evans, MV, Power, FW, Caldwell, JC (2007) Development and Use of PBPK Modeling and the Impact of Metabolism on Variability in Dose Metrics for the Risk Assessment of Methyl Tertiary Butyl Ether (MTBE), Journal of Environmental Protection Science, 1:29-51. [4] Brautbar, N, Williams, J (2002) Industrial liver toxicity: Risk assessment, risk factors and mechanisms, Int J Hyg Envir Heal 205(6):479-491. [5] Brown, RP, Delp, MD, Lindstedt, SL, Rhomberg, LR, Beliles, RP (1997) Physiological Parameter Values for Physiologically Based Pharmacokinetic Models. Toxicol Ind Health 13(4):407-484. [6] Caldwell, JC, Evans, MV, Krishnan, K (2012) Cutting Edge PBPK Models and Analyses: Providing the Basis for Future Modeling Efforts and Bridges to Emerging Toxicology Paradigms, Journal of Toxicology 2012: 852384 (doi:10.1155/2012/852384). [7] Cheng, S, Bois, FY (2011) A mechanistic model- ing framework for predicting metabolic interactions in complex mixtures. Environ Health Persp 119(12):1712- 1718. [8] Clewell, HJ, Gentry, PR, Covington, TR, Gearhart, JM (2000) Development of a Physiologically Based Pharmacokinetic Model of Trichloroethylene and Its Metabolites for Use in Risk Assessment, Environ Health Persp 108(Suppl 2):283-305. [9] Cole, CE, Tran, HT, Schlosser, PM (2001) Physio- logically Based Pharmacokinetic Modeling of Benzene Metabolism in Mice Through Extrapolation from In Vitro to In Vivo, J Toxicol Env Heal A 62(6):439-465. [10] Csanády, GA, Kessler, W, Hoffmann, HD, Filser, JG (2003) A toxicokinetic model for styrene and its metabolite styrene-7,8-oxide in mouse, rat and human with special emphasis on the lung, Toxicol Lett 138(1- 2):75-102. [11] Cuello, WS, James, TAT, Jessee, JM, Venecek, MA, Sawyer, ME, Eklund, CR, Evans, MV (2012) Phys- iologically Based Pharmacokinetic (PBPK) Model- ing of Metabolic Pathways of Bromochloromethane in Rats, Journal of Toxicology 2012: 629781 (doi: 10.1155/2012/629781). [12] Dobrev, ID, Andersen, ME, Yang, RSH (2002) In silico toxicology: simulating interaction thresholds for human exposure to mixtures of trichloroethylene, tetra- chloroethylene, and 1, 1, 1-trichloroethane. Environ Health Persp 110(10):1031-1039. [13] El-Masri, HA, Portier, C (1198) Physiologically Based Pharmacokinetics Model of Primidone and Its Metabo- lites Phenobarbital and Phenylethylmalonamide in Hu- mans, Rats, and Mice, Drug Metab Dispos 26(6):585- 594. [14] Evans, MV, Crank, WD, Yang, H, Simmons, JE (1994) Applications of a sensitivity analysis to a physiologi- cally based pharmacokinetic model for carbon tetrachlo- ride in rats. Toxicol Appl Pharm 128:36-44. [15] Fisher, JW, Mahle, D, Abbas, R (1998) A Hu- man Physiologically Based Pharmacokinetic Model for Trichloroethylene and Its Metabolites, Trichloroacetic Acid and Free Trichloroethanol, Toxicol Appl Pharm 152:339-359. [16] Gagnaire, F, Langlais, C, Grossman, S, Wild, P (2007) Ototoxicity in rats exposed to ethylbenzene and to two technical xylene vapours for 13 weeks, Arch Toxicol 81(2):127-143. [17] Greenberg, MS, Burton, GA, Jr., Fisher, JW (1999) Physiologically Based Pharmacokinetic Modeling of Inhaled Trichloroethylene and Its Oxidative Metabolites in B6C3F1 Mice, Toxicol Appl Pharm, 154:264-278. [18] Haddad, S, Beliveau, M, Tardif, R, Krishnan, K (2001) A PBPK Modeling-Based Approach to Account for Interactions in the Health Risk Assessment of Chemical Mixtures, Toxicol Sci 63(1):125-131. [19] Haddad, S, Krishnan, K (1998) Physiological model- ing of toxicokinetic interactions: implications for mix- ture risk assessment. Environ Health Persp 106(Suppl 6):1377-1384. [20] Haddad, S, Tardif, Charest-Tardif, G, Krishnan, K (1999) Physiological modeling of the toxicokinetic in- teractions in a quaternary mixture of aromatic hydrocar- bons. Toxicol Appl Pharm 161(3):249-257. [21] Krishnan, K, Clewell, HJ 3rd, Andersen, ME (1994) Physiologically based pharmacokinetic analyses of sim- ple mixtures. Environ Health Persp 102(Suppl 9):151- 155. [22] Krishnan, K, Haddad, S, Béliveau, M, Tardif, R. (2002) Physiological modeling and extrapolation of pharma- cokinetic interactions from binary to more complex chemical mixtures. Environ Health Persp 110(Suppl 6):989-994. [23] Loizou, GD (2001) The application of physiologically based pharmacokinetic modelling in the analysis of oc- cupational exposure to perchloroethylene. Toxicol Lett 124:59-69. [24] Loizou, GD, McNally, K, Jones, K, Cocker, J (2015) The application of global sensitivity analysis in the development of a physiologically based pharmacokinetic model for m-xylene and ethanol co-exposure in humans. Frontiers in Pharmacology 6(135):1-19. [25] Manning, CC, Schlosser, PM, Tran, HT (2010) A Mul- ticompartment Liver-based Pharmacokinetic Model for Benzene and its Metabolites in Mice, B Math Biol, 72(3):507-540. [26] Marchand, A, Aranda-Rodriguez, R, Tardif, R, Nong, A, Haddad, S (2015) Human Inhalation Exposures to Toluene, Ethylbenzene, and M-Xylene and Physiolog- ically Based Pharmacokinetic Modeling of Exposure Biomarkers in Exhaled Air, Blood, and Urine. Toxicol Sci 144(2):414-424. Biomath 8 (2019), 1901067, http://dx.doi.org/10.11145/j.biomath.2019.01.067 Page 12 of 13 http://dx.doi.org/10.11145/j.biomath.2019.01.067 K.A. Yokley, J. Ashcraft, N.S. Luke, A computational investigation of the ventilation structure and ... [27] Matlab, http://www.mathworks.com/ [28] fminsearch.m, Matlab Central File Exchange, http://www.mathworks.com/matlabcentral/fileexchange/ 8277-fminsearchbnd--fminsearchcon/ [29] McNally, K, Cotton, R, Cocker, J, Jones, K, Bartels, M, Rick, D, Price, P, Loizou, G (2012) Reconstruction of exposure to m-Xylene from human biomonitoring data using PBPK modelling, Bayesian inference, and Markov chain Monte Carlo simulation. Journal of Toxicology 2012 (doi:10.1155/2012/760281). [30] McNally, K, Cotton, R, Loizou, GD (2011) A workflow for global sensitivity analysis of PBPK models. Fron- tiers in pharmacology 2(31):1-22. [31] Mendoza-Cantú, A, Castorena-Torres, F, Bermúdez de León, M, Cisneros, B, López-Carrillo, L, Rojas-Garcı́a, AE, Aguilar-Salinas, A, Manno, M, Albores, A (2006) Occupational toluene exposure induces cytochrome P450 2E1 mRNA expression in peripheral lymphocytes. Environ Health Persp 114(4):494-499. [32] Peyret, T, Krishnan, K (2012) Quantitative property- property relationship for screening-level prediction of intrinsic clearance of volatile organic chemicals in rats and its integration within PBPK Models to predict inhalation pharmacokinetics in humans. Journal of Tox- icology 2012 (http://dx.doi.org/10.1155/2012/286079). [33] Price, K, Krishnan, K (2011) An integrated QSAR- PBPK modelling approach for predicting the inhalation toxicokinetics of mixtures of volatile organic chemicals in the rat. SAR QSAR Environ Res 22(1-2):107-128. [34] Ramsey, JC, Andersen, ME (1984) A physiologically based description of the inhalation pharmacokinetics of styrene in rats and humans, Toxicol Appl Pharm, 73(1):159-175. [35] Sarangapani, R, Teeguarden, JG, Cruzan, G, Clewell, HJ, Andersen, ME (2002) Physiologically based phar- macokinetic modeling of styrene and styrene oxide respiratory-tract dosimetry in rodents and humans, Inhal Toxicol 14(8):789-834. [36] Sawyer, ME, Evans, MV, Wilson, CA, Beesley, LJ, Leon, LS, Eklund, CR, Croom, EL, Pegram, RA (2016) Development of a human physiologically based phar- macokinetic (PBPK) model for dermal permeability for lindane, Toxicol Lett 245: 106-109. [37] Simmons, JE, Boyes, WK, Bushnel, PJl, Raymer, JH, Limsakun, T, McDonald, A, Sey, YM, Evans, MV (2002) A Physiologically Based Pharmacokinetic Model for Trichloroethylene in the Male Long-Evans Rat, Toxicol Sci 69:3-15. [38] Tardif, R, Charest-Tardif, G, Brodeur, J, Krishnan, K (1997) Physiologically based pharmacokinetic modeling of a ternary mixture of alkyl benzenes in rats and humans. Toxicol Appl Pharm 144(1):120-134. [39] Tardif, R, Laparé, S, Charest-Tardif, G, Brodeur, J, Kr- ishnan, K (1995) Physiologically-based pharmacokinetic modeling of a mixture of toluene and xylene in humans. Risk Anal 15(3):335-342. [40] Tardif, R, Laparé, S, Krishnan, K, Brodeur, J (1993) Physiologically based modeling of the toxicokinetic interaction between toluene and m-xylene in the rat. Toxicol Appl Pharm 120(2):266-73. [41] Tassaneeyakul, W, Birkett, DJ, Edwards, JW, Veronese, ME, Tassaneeyakul, W, Tukey, RH, Miners, JO (1996) Human cytochrome P450 isoform specificity in the regioselective metabolism of toluene and o-, m-and p- xylene. J Pharmacol Exp Ther 276(1):101-108. [42] Thrall, KD, Gies, RA, Muniz, J, Woodstock, AD, Hig- gins, G (2002) Route-of-entry and brain tissue partition coefficients for common superfund contaminants. J Tox- icol Env Heal A 65:2075-2086. [43] Valcke, M, Haddad, S (2015) Assessing human variabil- ity in kinetics for exposures to multiple environmen- tal chemicals: a physiologically based pharmacokinetic modeling case study with dichloromethane, benzene, toluene, ethylbenzene, and m-xylene, J Toxicol Env Heal A 78: 409-431. [44] Yang, RS (1998) Some critical issues and concerns related to research advances on toxicology of chemi- cal mixtures. Environ Health Persp 106(Suppl 4):1059- 1063. [45] Yang, RS, Thomas, RS, Gustafson, DL, Campain, J, Benjamin, SA, Verharr, HJ, M.M. Mumtaz, MM (1998) Approaches to developing alternative and predictive toxicology based on PBPK/PD and QSAR modeling. Environ Health Persp 106(Suppl 6):1385-1393. [46] Yokley, KA (2013) Investigations on Ventilation Equa- tion Structure in Physiologically Based Pharmacokinetic (PBPK) Modeling of Inhaled Toxicants. International Journal of Pure and Applied Mathematics 82(1):95-123. [47] Yokley, KA, Evans, MV (2008) Physiological changes associated with aging result in lower internal doses of toluene and perchloroethylene in simulations using phar- macokinetic modeling. Toxicol Environ Chem 90(3):75- 492. [48] Yokley, KA, Evans, MV (2007) An Example of Model Structure Differences Using Sensitivity Analyses in Physiologically Based Pharmacokinetic Models of Trichloroethylene in Humans, B Math Biol 69(8):2591- 2625. Biomath 8 (2019), 1901067, http://dx.doi.org/10.11145/j.biomath.2019.01.067 Page 13 of 13 http://www.mathworks.com/ http://www.mathworks.com/matlabcentral/fileexchange/8277-fminsearchbnd--fminsearchcon/ http://www.mathworks.com/matlabcentral/fileexchange/8277-fminsearchbnd--fminsearchcon/ http://dx.doi.org/10.11145/j.biomath.2019.01.067 Introduction Model Investigation Model Background Fixed Parameter Values Vmaxc values Investigational Methods Results Discussion and Conclusions Acknowledgements Appendix References