A computational investigation of the ventilation structure and maximum rate of metabolism for a physiologically based pharmacokinetic (PBPK) model of inhaled xylene

Physiologically based pharmacokinetic (PBPK) models are systems of ordinary differential equations that estimate internal doses following exposure to toxicants. Most PBPK models use standard equations to describe inhalation and concentrations 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 can result in similar compartmental predictions when used with different equations describing inhalation. Simulations are performed using values based on existing literature. Simulated data is also used to determine specific 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.


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 [ [48]. Most PBPK model equations use the same basic structure and assumptions, such as the assumption that compartments are well-mixed and that the transfer of some chemicals is at equilibrium. 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, V max (mg/h). Yokley [46] considered three structures for modeling ventilation, described as "equilibrium," "non-equilibrium," and "hybrid" models. The results were very similar for the equilibrium and hybrid models, and therefore only two structures will be used to model ventilation in the current study. The equilibrium model uses the standard quotient for the concentration of toxicant 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.

A. Model Background
In [46], a PBPK model for xylene exposure in humans from [39] was used with different  [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 compartment 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 Appendix.
The amount of xylene within each compartment depicted in Figure 1 is modeled by a separate differential equation. Additional equations are added to the model to incorporate the ventilation structure. For each model equation, the compartmental concentration is defined as where C xyl j represents the concentration of the chemical, xylene, in compartment j (mg/L); A xyl j represents the amount of the chemical, xylene, in compartment j (mg); and V j 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 compartment j (excluding the liver) takes the form where Q j represents the rate of blood flow to compartment j (L/h); C xyl art is the concentration of the chemcial, xylene, in the arterial blood (mg/L); and P xyl j 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 metabolized 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 where V xyl max is the maximum rate of metabolism for xylene (mg/h) and K xyl M 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 represent the ventilation structures, which are outlined below.
The three ventilation structures used in [46] were classified as "equilibrium," "nonequilibrium," 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: where Q c represents the cardiac flow rate, defined as the sum of the compartmental flow rates (L/h); Q p is the pulmonary flow or alveolar ventilation (L/h); and P bl:air is the blood/air partition coefficient (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 which are described more fully in [46]. The "hybrid" case involves using the "non-equilibrium" equations (6)-(8) with the alveolar partition coefficient, P i alv , set to be 1. Because the results from [46] were very similar between the "equilibrium" and "hybrid" cases, the "hybrid" scenario is omitted 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 is a sufficient conversion. Physiological parameters used in the PBPK model are presented in Tables I and II. As described in [46], the lung:air partition coefficient (P xyl alv:air ) was obtained from [42] and was used with the blood:air partition coefficient from [39] in order to find P xyl alv (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 xyl max = V xyl maxc BW 0.75 K xyl m = 0.2, and cardiac output and alveolar ventilation are assumed to follow allometric scaling with scaling factors used from the source model

C. V maxc values
Metabolism of xylene is described using Michaelis-Menten kinetics in the model [39] which involves a nonlinear term with the parameter V xyl max representing the maximum rate of the enzymatic process (see Equation (3)). The value of V xyl max is scaled to body weight as follows: Different values of V xyl maxc have been used in PBPK modeling and other research as is shown in Table III. Parameter values have been rescaled for  [47] EXCEPT FOR Vr , V bl , V alv , Vart, AND Vven. Vr WAS ALTERED TO HAVE ALL VOLUMES SUM TO 100%. V bl V alv , Vart, AND Vven ARE BASED ON [5].
Volume Value consistency to represent V xyl maxc 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 xyl maxc values used in the simulations. V xyl maxc 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 xyl maxc , 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 xyl maxci 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  [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], [ 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 concentration 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 xyl maxc 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 xyl maxc values could produce similar results with the two different models, the following procedure is employed. Simulated data is generated from each model using a beginning V xyl maxc 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 xyl maxc , minimizing a least squares cost function where x i represents the state of the differential equation model (using a particular value of V xyl maxc ) at time t i corresponding to data point d i . Optimization is performed using the fminsearchbnd function [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 exercise for this study, the maximum model outputs are generated for both the equilibrium and nonequilibrium models and compared to determine if there is any overlap. Figure 2 depicts the relationship between maximum predicted xylene amount in the exhaled air and values of V xyl maxc for both the equlibrium and non-equlibrium models. It is evident from this graph that there is little to no overlap in the maximum predicted xylene amount for both models.  The maximum amounts generated from the nonequilibrium model are consistently higher than those of the equilibrium model. The significant difference between the model outputs would suggest 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 between values of V xyl maxc with the maximum concentration of xylene in the venous blood and maxi-  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 either the equilibrium model or the non-equilibrium model with different values for V xyl maxc . A similar relationship for the area below the curve estimates for concentation of xylene in the liver is presented in Figure 5.  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 efficiently 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 exhaled air using simulated data from the equilibrium model are contained in Figure 6, and examples using 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

Noise
Best Fit to "Equilibrium" Best Fit to "Non-equilibrium" air concentration is presented in Table IV.
In Figure 6, it can be observed that the equilibrium 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 equilibrium model produces a least squares cost of 2.1612e −6 and the non-equilibrium model produces 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 xyl maxc in the range of 3-7. Examples of graphical results are contained  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 equilibrium and non-equilibrium models to venous blood concentration data that was generated by the equilibrium model with a V xyl maxc 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 xyl maxc 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 nonequilibrium 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

Noise
Best Fit to "Equilibrium" Best Fit to "Non-equilibrium" non-equilibrium data. The non-equilibrium data in these graphs were generated with a V xyl maxc value of 6. The best fit of the equilibrium model to the data with 2% noise is found with V xyl maxc equal to 3.3163 and has a least squares cost of 0.0201. A V xyl maxc 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 xyl maxc 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 exposure 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 Figure 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 xyl maxc values could be used to make similar predictions for the concentration of xylene in venous blood with the different ventilation structures, but different V xyl maxc values are not found to produce similar model output for the amount in exhaled air. Specifically, using V xyl maxc around 3 in the "equilibrium" model produces similar blood concentration results as using V xyl maxc 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 xyl maxc 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 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 calculated using the "equilibrium" model with V xyl maxc 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 problematic for use in optimizing metabolic parameters. Additionally, the model used in the current study is a more simplistic PBPK model, and more investigation 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 University Funding for Research and Development for support for this project.