Physiologically Based Pharmacokinetic (PBPK) Modeling of Metabolic Pathways of Bromochloromethane in Rats

Abstract

Bromochloromethane (BCM) is a volatile compound and a by-product of disinfection of water by chlorination. Physiologically based pharmacokinetic (PBPK) models are used in risk assessment applications. An updated PBPK model for BCM is generated and applied to hypotheses testing calibrated using vapor uptake data. The two different metabolic hypotheses examined are (1) a two-pathway model using both CYP2E1 and glutathione transferase enzymes and (2) a two-binding site model where metabolism can occur on one enzyme, CYP2E1. Our computer simulations show that both hypotheses describe the experimental data in a similar manner. The two pathway results were comparable to previously reported values (𝑉_max=3.8 mg/hour, 𝐾_𝑚=0.35 mg/liter, and 𝑘_GST=4.7 /hour). The two binding site results were 𝑉_max1=3.7 mg/hour, 𝐾_𝑚1=0.3 mg/hour, CL₂ = 0.047 liter/hour. In addition, we explore the sensitivity of different parameters for each model using our obtained optimized values.

1. Introduction

Bromochloromethane (BCM, CH₂BrCl, CAS number 83847-49-8) is a volatile solvent historically used in portable fire extinguishers. BCM was gradually replaced by halons in the 1970s and banned in 2002 due to its ozone depletion potential. However, BCM is still used as an intermediate in the production of other solvents [1] and is rapidly absorbed via inhalation due to its volatility [2]. Similar to other halomethanes, there is a suggestion of hepato- and nephro-toxicity due to prolonged exposure [2]. Another potential source of exposure for BCM is through oral ingestion via water, since BCM is considered a water disinfection byproduct [3], and is listed within the U.S.A EPA’s Candidate Contaminant List (CCL). Additional research on brominated disinfection by-products (DBPs) has been proposed because there is evidence that brominated species can be more potent than chlorinated compounds [4]. The occurrence of brominated disinfection by-products (e.g., bromochloromethane) is an added consideration for coastal areas due to potential bromine intrusion from seawater and subsequent transformation of chlorinated by-products into brominated species [5]. Given the potential health impact of brominated chemicals, we chose to study BCM in consideration of both increased potential for exposure near coastline utilities, and on the basis of a deficiency and availability of published data. Based on a recent PubMed search, the published literature for BCM is much less than that of the well-studied bromodichloromethane (BDCM, CHBrCl₂, CAS number 75-27-4), another structurally related DBP.

Physiologically based pharmacokinetic (PBPK) models are computational tools that are increasingly used to address risk assessment issues, particularly in quantifying the relationship between measures of external exposure and internal dose [6]. PBPK models are sets of equations representing the flow of blood and toxicant into and out of organs and body tissues. Unlike classical pharmacokinetic models, PBPK models include species-specific physiological, chemical, and biochemical parameters, allowing for extrapolation to humans [7]. The interaction between modelers and biologists is crucial in determining the final representation of organs to ensure inclusion of target organs, special physicochemical properties and toxicity. An accurate mathematical description of biological processes (i.e., an accurate model) requires synthesis of mathematics, modeling, and systems biology. The expectation of such an approach is that an integrated modeling effort and its application to quantification of health risk assessment would increase our confidence in the simulated results.

For volatile halogenated methanes, metabolism is an important component studied to provide explanations of how toxicity occurs. The complex kinetics suggested for dihalomethane metabolism (e.g., BCM) consists of two components: a saturable Michaelis-Menten hyperbola superimposed with a linear pathway [9]. This linear pathway does not appear to saturate and shows a continued proportional increase at higher exposure concentrations. At lower concentrations, the contribution of this linear pathway is smaller than that of the saturable component. The question becomes what are the physiological counterparts of these two different metabolism components? What enzymes may exhibit such different characteristics? Previous PBPK modeling work with dichloromethane (DCM, CH₂Cl₂, CAS number 75-09-2) led to the hypothesis that two different enzymes are involved in its metabolism: CYP2E1, a P450 isoenzyme, and GST (GSH S-transferase). The early PBPK modeling work was followed by experiments showing that the GST pathway, although numerically small, was an important pathway toxicologically. Based on experimental evidence, the GST pathway is believed to be responsible for the carcinogenicity of these compounds [15]. The early PBPK modeling paved the way for establishing a two-pathway hypothesis which includes both P450 and glutathione transferase (GST) enzymes [9]. Based on these concepts, the proposed metabolic scheme for bromochloromethane is shown in (Figure 1). Specifically for BCM, the two-pathway hypothesis was first suggested by PBPK modeling followed by experimental confirmation of adduct formation from the GST pathway in bacterial assays [16]. Using an in vitro Salmonella assay, BCM was confirmed to be mutagenic in the presence of cytosol (GST containing fraction) relative to other structurally related halomethanes [17]. Based on structural similarity with other brominated ethanes (i.e., 1,2-dibromoethane), where plasma bromide was used to track metabolism, the existence of a small but relevant GST pathway is to be expected for BCM [18]. However, we were not able to find in vivo studies that can be used to quantify the linear pathway for BCM, other than suggestions that in vivo data is lacking [13].

Figure 1 

Schematic representation of BCM metabolism, adapted from [13].

Other than having two different enzymes describe the complex kinetics, is there another plausible explanation for the linear kinetics exhibited by the in vivo closed chamber data? If so, can PBPK modeling be used to explain the additional linear kinetics? In recent years, the realization has been made that several P450s may exhibit atypical kinetics, that is, complex kinetics that cannot be explained by the single saturable Michaelis-Menten kinetics [19]. In an attempt to better describe these increased complexities, biochemical models with two sites have been developed to describe the linear kinetics exhibited at higher concentrations [20]. These complex biochemical models have an additional binding site, which is simultaneously available to the substrate (Figure 2). Once binding to the second site occurs, a conformational change takes place that modifies the metabolism of the first site. We propose that a PBPK model for BCM can be constructed to show the plausibility of a second site model. At this point in time, we are not suggesting that one kinetic model is more realistic than the other, but are investigating the possibility of alternative pathways based on atypical kinetics and its implications for metabolism. Interestingly, one potential way to distinguish between the two mechanism is by the type of metabolites excreted. With two-sites in a single enzyme, the metabolites produced would probably be reactive products of CYP2E1 oxidation. In contrast, a two-pathway mechanisms would have glutathione containing metabolites during excretion.

Figure 2 

Comparison of single-sited and two-sited kinetics, adapted from [14].

When testing a new hypothesis using PBPK model simulations, one possible question becomes how to rank the importance of relative metabolic parameters. We believe that this question is particularly relevant when trying to distinguish between metabolic constants that are numerically small. Sensitivity analysis is a mathematical tool that allows determination of how the change in one parameter affects model predictions. When conducting sensitivity analysis, the model is run repeatedly to determine how changes in parameter values relate to changes in a prediction. Although optimization is performed before sensitivity analysis, modeling behavior insights come from the sequential sensitivity analysis. The sensitivity analysis results will also be dependent on the values obtained from the optimization. Both techniques are combined since the first is needed to find the “correct” parameter values, while sensitivity analysis looks for impact across parameters. In the present work, both techniques were used to identify potential issues with metabolic parameter quantification. Once the optimization is performed, sensitivity analysis can determine if a unique set of values exist for the optimized results. In general, sensitivity analysis has been used in combination with PBPK modeling to estimate the relative ranking of different model parameters [21], sometimes including the interaction between parameters [22]. Local sensitivity analysis is the simplest form of sensitivity analysis, which allows for estimation of the impact of each parameter individually. One potential analogy is that this type of sensitivity analysis is similar to a worst-case scenario estimate for the variability of each parameter. In addition to ease of implementation, another advantage of this analysis is its ability to sort potential identifiability issues, or to state whether the optimization provides unique estimates for the resulting values obtained for the different modeling hypotheses.

In summary, the goals of this study were (1) to use PBPK modeling in combination with BCM closed chamber data to test two potential metabolic hypotheses: a two-pathway enzyme combination of CYP2E1 and GST and a one enzyme mechanism consisting of two binding sites for CYP2E1; (2) to make use of sensitivity analysis to identify important model parameters; (3) to compare the fits across different metabolic hypotheses and suggest future experiments that lead to increased confidence in our simulations.

2. Methods

2.1. PBPK Models

2.1.1. Closed Chamber Data

Inhalation closed chamber data for BCM in rats were obtained from the published literature [9] digitized using UN-SCAN-IT 6.0 (Silk Scientific, Orem, UT, USA). Inhalation closed chamber data were collected from a sealed system, allowing for quantification of metabolism. A BCM bolus was injected into the chamber at time zero. After equilibrium between the rat and air inside the chamber, BCM decreased with time as it was distributed into tissues and metabolized by the rodent. A PBPK model was then constructed to describe the decline in air chamber concentration. The metabolic parameters were calculated using optimization techniques in combination with the inhalation data. Different metabolic hypotheses were investigated by using a PBPK model with the closed chamber-derived metabolic parameters to predict metabolism.

2.1.2. Core Model

A flow-limited PBPK model based on the structure described by Ramsey and Andersen [23] was used to describe BCM metabolism (Figure 3). A seven-compartment model (lung, blood, adipose, rapidly perfused tissue, slowly perfused tissue, liver, and kidney) was developed for BCM, with the assumption that metabolism occurred via the same enzymes in liver and kidney. All physiological parameters were obtained from the literature as shown in Table 1. Differential equations based on mass conservation principles were derived for each organ. Lung and blood compartments use a steady state assumption, simplifying the overall solutions for the model. Tissue to blood partition coefficients were calculated using experimentally-derived BCM coefficients from Gargas et al., 1986 [9] (described in the appendix). After derivation, all equations were coded into MatLab (MathWorks, version 7.11.0.584) using the equation solver ode15s.

Table 1 

Parameters for adult Fisher 344 rats and corresponding significance in the PBPK model. (Note: partition coefficients are unitless [8].)

Figure 3 

Schematic representation of PBPK model used for BCM.

2.1.3. Metabolic Hypotheses Testing

The current work tested two different metabolic hypotheses using the same basic PBPK model structure, varying only the equations for liver and kidney metabolism. Previously, Gargas et al. [9] evaluated typical Michaelis-Menten kinetics, with metabolism described by a single CYP2E1 binding site. The results of their simulations demonstrated that the addition of a linear term was necessary to describe the closed chamber data. Thus, a two-pathway model with two different enzymes, CYP2E1 exhibiting Michaelis-Menten kinetics and GST exhibiting linear kinetics, provided much better fits to the data than the model that incorporated only Michaelis-Menten kinetics [9]. Here, the two-pathway model (Michaelis-Menten and GST) described by [9] is compared to the two-binding site model described by Evans and Caldwell [24]. In the two-binding site model, a modified Michaelis-Menten equation that includes two binding sites for CYP2E1 is used to describe BCM metabolism.

2.1.4. Two-Pathway Kinetics

The two-pathway description uses Michaelis-Menten kinetics in addition to a linear term to account for GST metabolism. The equations listed below represent the additional metabolic term in both the liver (𝑀𝑒𝑡_𝑙) and the kidney (𝑀𝑒𝑡_𝑘). Liver is set to account for 94.8% of total body metabolism and the kidney 5.2% based on [25]:𝑀𝑒𝑡_𝑙=0.948𝑉^maxC_liv𝐾_𝑚+𝐶_liv+𝑘^GST𝐶_liv𝑉_liv𝑀𝑒𝑡_𝑘,=0.052𝑉^max𝐶_kid𝐾_𝑚+𝐶_kid+𝑘^GST𝐶_kid𝑉_kid,(1) where (i) 𝑉_max: maximum velocity of the reaction, mg/hour, (ii) 𝐾_𝑚: affinity constant, mg/liter, (iii) 𝑘_GST: proportionality constant for linear pathway metabolized by glutathione transferase, /hour.

2.1.5. Two-Binding Site Kinetics

Linear kinetics at higher concentrations can also be described using a dual binding site, as demonstrated for dichloromethane using closed chamber data [24]. In this case, CYP2E1 undergoes a structural change in which two separate binding sites are simultaneously available to BCM, and metabolism is achieved with one enzyme. For the dual binding site equation, the squared concentration term becomes asymptotically linear at higher exposure concentrations. This squared term is not present in the typical Michaelis-Menten equation and will lead to a different rate of ascent to the asymptote. The modified Michaelis-Menten equation reflecting two binding sites in one enzyme [20] becomes𝑀𝑒𝑡_𝑙=0.948⎛⎜⎜⎝𝑉^max1C_liv+CL₂𝐶_liv²𝑀𝑒𝑡_𝑘𝐾_𝑚1+𝐶_liv⎞⎟⎟⎠,=0.052⎛⎜⎜⎝𝑉^max1𝐶_kid+CL₂𝐶_kid²𝐾_𝑚1+𝐶_kid⎞⎟⎟⎠,(2) where (i) 𝑉_max1: maximum metabolic rate for the first binding site, mg/hour. The second binding site is accounted for by CL₂, (ii) 𝐾_𝑚1: affinity constant for the first binding site, mg/liter, (iii) CL₂: clearance constant consisting of the ratio of the maximum metabolic rate and the affinity constant, 𝑉_max2/𝐾_𝑚2, for second binding site, liter/hour.

2.2. Comparison of the Different Metabolic Models

A total of six metabolic parameters were optimized using five closed chamber data sets, each having a different starting concentration. Optimizations were performed using the natural log of the data with MatLab’s fminsearch function (MathWorks, version 7.11.0.584). For the two-pathway model, three parameters were optimized using the digitized data: 𝑉_max, 𝐾_𝑚, and 𝑘_GST. For the two-binding site model, three parameters were optimized using the digitized data: 𝑉_max1, 𝐾_𝑚1, and CL₂. The function fminsearch, available in the standard MatLab package, tries to minimize a “cost’’ as the minimum difference between data and simulation. The best overall fit minimized by fminsearch includes the entire data range (200–4000 ppm).

2.3. Physiological Value Selection

The physiological values to be used in the PBPK model have been updated from those used in the original simulations of the data in 1986. As an example, cardiac output values for F344 rats are now available [11]. The range of normalized cardiac output values for F344 rats values (QCC) is between 10–20 liters/hour/kg^0.75, due to dependence of cardiac output on age. We decided to use the mean value of 15 liters/hour/kg^0.75, since a 225 g rat is considered adult [10]. The ventilation rate for rats was calculated using the ventilation perfusion ratio, which varies between 1–4 fold due to differences in measurement methodology [22, Table 31]. Since the mean value of 15 liters/hour/kg^0.75 matched the values used in the benzene PBPK model dated from approximately the same time period as the Gargas dataset, we used the same ventilation perfusion rate as in Medinsky et al. [12] for consistency. The F344-specific tissue perfusion and volume values were adopted from [11].

2.4. Sensitivity Analysis

Sensitivity coefficients were calculated with partial derivatives for each variable of interest with respect to the model parameter being specified in the model. Sensitivity coefficients were then normalized by both the variable and model parameter. The resulting time course for each sensitivity coefficient was plotted. Sensitivity coefficients were calculated using Automated Differentiation written by Martin Fink and available at MatlabCentral(http://www.mathworks.com/matlabcentral/). Another Matlab function written by Adam Attarian (tssolve NC State University, Raleigh, NC, USA) was used to organize the partial derivatives and for the final plots. A three-dimensional surface was generated to show variation in sensitivity coefficients with time, and variability within the metabolic parameter. These tests were performed to check stability of the sensitivity coefficients when introducing variability around the optimized estimates.

3. Results

The optimization results are shown in Table 2. For all concentrations studied, the two-pathway and two-binding site models gave very similar descriptions of the closed chamber data (Figure 4). Both the 𝑉_max and 𝐾_𝑚 values are similar for the two-pathway and two-binding models. The two-pathway parameter values were as follows: 𝑉𝐾_max=3.8mg/hour,_𝑚𝑘=0.35mg/liter,_GST=4.7/hour.(3) The two-binding site hypothesis uses Michaelis-Menten parameters 𝑉_max1 and 𝐾_𝑚1 to describe the first site, and clearance CL₂, defined to be the ratio of 𝑉_max2 and 𝐾_𝑚2, for the second site. The two-binding site results were as follows: 𝑉𝐾_max1=3.7mg/hour,_𝑚1CL₂=0.3mg/liter,=0.047liters/hour.(4)

Table 2 

Optimized parameter results.

(a)

(b)

Figure 4 

PBPK modeling results using two different metabolic hypotheses. (a) shows results for the two-pathway model. (b) shows results for the two-binding site model.

When comparing the values for both metabolic hypotheses, Michaelis-Menten results for the two-binding site hypothesis are very similar to the Michaelis-Menten values for the two-pathway hypothesis. The 𝑉_max and 𝐾_𝑚 values are similar for the two pathway and two-binding models. Based on these results, both metabolic hypotheses are considered equally plausible.

Normalized sensitivity coefficients for all model parameters were calculated after optimization results were obtained. The majority of resulting sensitivity coefficients were not shown, since we are focusing on metabolic parameters: 𝑉_max, 𝐾_𝑚, 𝑘_GST, or CL₂. The authors performed sensitivity analysis using 500 ppm as the initial concentration. Each sensitivity figure consists of two panels. The first panel presents sensitivity coefficients calculated for the impact of air chamber concentration on the simulated value of 𝑉_max and the second panel presents sensitivity coefficients for liver concentration of BCM. Figure 5 presents sensitivity results for the two-pathway model. The sensitivity coefficients for 𝑉_max have the highest values. The 𝑉_max sensitivity coefficients for liver concentration are larger than those obtained using chamber air as the experimental variable. Figure 6 presents sensitivity results for the two-binding site model. The sensitivity coefficients for 𝑉_max also have the highest values, particularly when comparing liver coefficients versus chamber air coefficients. For both model hypotheses, peak 𝑉_max sensitivity occurs within the first two hours.

Figure 5 

Sensitivity coefficients for the two-pathway model at 500 ppm. Each concentration includes the analysis using air and liver concentrations.

Figure 6 

Sensitivity coefficients for the two-binding site model at 500 ppm. Each concentration includes the analysis using air and liver concentrations.

A three-dimensional plot of liver sensitivity coefficients for the two-pathway model is presented in Figure 7. The third dimension represents variability in 𝑉_max generated by repeating sensitivity analysis simulations that vary 𝑉_max between 2.5–4.5 mg/hour (to include the mean value 3.8 mg/hour). Based on the value of 𝑉_max, there is a peak in sensitivity to 𝑉_max occurring within the first 3 hours. The chamber air concentration plot shows a much shallower surface, with no obvious peaks in sensitivity towards 𝑉_max shown.

Figure 7 

Three-dimensional sensitivity analysis plot for the two-pathway model at 500 ppm.

4. Discussion

Risk assessments for volatile chemicals have used PBPK models to account for differences across routes and extrapolate from rodent experiments to human populations [26, 27]. The use of a PBPK model allows for the integration of physiological, chemical, and biochemical information with different kinetic hypotheses. Closed chamber inhalation data are a measure of total rate of metabolism represented by the decrease in the concentration of the volatile chemical inside the chamber after a bolus injection. One application of closed chamber inhalation data in combination with PBPK modeling has been the suggestion of metabolic hypotheses leading to effects in a target organ, particularly for halomethanes [9, 28]. In this context, PBPK models have been proposed as useful computational tools for risk assessments.

PBPK models have played a decisive role in the suggestion of an in vivo role for the GST pathway (two-pathway model) and its contribution to toxicity for halogenated compounds [15, 28]. In the case of DCM, the GSH conjugation pathway has a small proportionality constant (∼0.01/hour [24]), while the GST pathway for BCM is larger (∼5.3/hour [9]). Following the inclusion of a linear pathway in the PBPK model for halogens, in vitro work in bacterial systems provided the first evidence that GST was the enzyme associated with the additional linear pathway. In vitro work also suggested that as the number of bromines included in the compound increased, the resulting genotoxicity was larger. For example, Thier et al. [16] demonstrated that the genotoxicity of CH₂Cl₂ (DCM) was less than that of CH₂ClBr (BCM), which was less than that of CH₂Br₂ (dibromomethane) toxicity. Prior to these experiments, there was no experimental evidence that dihalomethanes such as BCM were metabolized via a GST pathway. The only suggestion for the existence of this second pathway came from PBPK modeling and closed chamber data.

The current PBPK model making use of a two-pathway mechanism gives a similar 𝑘_GST value to that of previous models (Table 2 versus [9, 29, 30]). (It is important to note that both PBPK models in [29, 30] were based on [9]. In addition, the 𝑉_max and 𝐾_𝑚 values are slightly different when compared to [9]: 2.57 mg/hour and 0.3 mg/liter, resp.). We attribute the difference in optimized values to differences in physiological constants used in our model. The current PBPK model made use of F344 specific values [11], which correlates to the strain used in the closed chamber experiments. It is important to note that the F344 physiological values were not available at the time when the first BCM model was published (1986). The largest differences in volume or flow for the compartments used in the two different BCM PBPK models probably lie in the lumped compartments, namely, the rapidly and slowly perfused compartments.

To our knowledge, there is no previous PBPK work that makes use of a two-binding site mechanism for BCM, making a comparison of our linear CL₂ constant with previous models not possible. An important realization of this PBPK modeling effort is that both metabolic hypotheses appeared equally plausible. In order to help convince ourselves that this was the case, the available BCM in vivo literature was studied. Early in vivo experiments with halomethanes confirmed the enzymatic role for GST for other iodinated or brominated compounds but not for BCM [31]. This paper refers to earlier work and states that specifically for BCM, a BCM-GST reaction is not present in measurable amounts. This early indication concludes that GST is not involved as a major component of BCM metabolism. Another indication that the GST pathway is very small in vivo comes from the vapor uptake experiments, where treatment with pyrazole almost completely suppressed BCM metabolism as shown in Figure 8, from [9]. Since CYP2E1 has been shown to be the major P450 contributing to its oxidative metabolism [32], the in vivo metabolic decrease observed at this high concentration argues for the predominance of CYP2E1 over other P450s. In addition, the in vitro kinetic values reported by [16] also suggest a small kinetic constant when compared to ethylene dibromide or dibromomethane. The ability to numerically estimate metabolic constants for the linear pathway using sensitivity analysis will be discussed later.

Figure 8 

BCM vapor uptake with pyrazole treatment [9].

Several P450s (2A1, 2B1, and 2E1) have shown atypical kinetics, explained by allosteric mechanisms making use of a second site at least in vitro [33]. Li et al. [33] performed simulations to confirm the relative small size of CYP2E1’s active site, stating that for this enzyme the reaction proceeds by the cooperative binding of two substrates simultaneously to form a ternary complex. The substrate that binds in the active site is usually called active substrate, while the other one that binds in the effector site, and that controls the oxidation of the active substrate is called effector substrate [33].

Effector substrates are becoming increasingly studied to offer potential explanations for drug-to-drug interactions. CYP2E1 specifically plays a diverse role in physiology, toxicity, and metabolism. For example, CYP2E1 is a P450 isoform with an important role in glycogenesis, and its physiological role is increasingly being recognized [34]. At present, this isoform is known to metabolize about 70 different compounds of varying size, including alcohols, ketones, nitrosamines, anesthetics, and even long-chain fatty acids. One of the questions recently resolved by studying human CYP2E1 structure is based on addressing how is it possible for this enzyme, with a small active site, to metabolize substrates ranging in a wide number of sizes (including fatty acids). In fact, it is the suggestion of a channel structure available as a second site, and in proximity to the active site, that becomes a solution to help explain CYP2E1’s flexibility to accommodate different sized substrates [35]. Although BCM is a small molecule, probably able to fit in the active site, the existence of a second site does describe linear kinetics at higher exposure concentrations. The second site model thus became a candidate for inclusion in a PBPK model calibrated using in vivo data for BCM. Since the second binding site in CYP2E1 is involved in fatty acid metabolism, we propose that BCM binding to this second site may impair fatty acid metabolism and may indirectly prevent normal metabolism of larger molecules that may be already impaired in chronic illnesses such as diabetes. This proposed effect is hypothesized to occur in addition to the glutathione-based toxicity based on a two-pathway hypothesis.

When using optimization to estimate metabolic model parameters, correlation can be an important issue that may impact the ability to estimate unique values for the parameters. Basic algebra reminds us that a system of equations can be uniquely solved when the number of parameters is equal to the number of equations, assuming that the estimates are independent of each other. In general, the unique estimation of 𝐾_𝑚 will depend on having an independent estimate for 𝑉_max. Since 𝐾_𝑚 is interrelated to 𝑉_max, experiments with multiple concentrations are needed, and the concentrations used must be high enough that metabolic saturation is reached.

One of the applications of local sensitivity analysis is to help illustrate the existence of these potential correlations between parameters because this affects the ability to solve for unique estimates. For more complex metabolic models, the relationship between metabolic parameters may not be obvious. Sensitivity analysis theory states that parameters are uniquely identifiable when the sensitivity coefficients cannot add to zero [36]. For example, if all sensitivity coefficients examined are positive, then their sum cannot be equal to zero. Using this reasoning, these authors examined the time course plots for 500 ppm and the sensitivity coefficients for both metabolic hypotheses. The sensitivity coefficients can be estimated for the different variables included in the model, allowing the analysis to suggest improvements in the experimental design. For this reason, liver concentration was selected as a possible experimental variable. When compared with chamber concentration, liver concentration resulted in increased sensitivity towards 𝑉_max. Since sensitivity coefficients for other metabolic parameters were nearly zero, 𝑉_max is uniquely identifiable when using closed chamber experiments. At higher concentrations, the sensitivity coefficients for both 𝑘_GST or CL₂ increase with exposure time (results not shown). This current sensitivity analysis corroborates the importance of these two additional parameters at higher concentrations.

We utilized a method using local sensitivity analysis as an initial step to help determine identifiability issues. When unique identifiability is guaranteed, the ability to have one unique solution is also guaranteed. Therefore, local sensitivity analysis can be seen as a first step towards increasing our confidence in the parameters obtained as a solution. The results of our sensitivity analysis suggest that measuring liver concentration would provide improved estimates for 𝑉_max at intermediate concentrations and improved estimates for the linear constants 𝑘_GST or CL₂ at high concentrations. To combat uncertainty in the estimates for 𝑉_max and other metabolic parameters, we implemented a novel approach by estimating the three-dimensional surface for the sensitivity coefficient changing with time, and including variability in the metabolic estimate of concern. If future experiments include liver concentration measurements, our three-dimensional analysis at 500 ppm suggests that peak information on 𝑉_max occurs before 2 hours, suggesting that experiment duration can be shortened by several hours. The application of PBPK modeling and three-dimensional sensitivity analyses can be helpful to design future experiments aimed at refining metabolic estimates.

The combination of optimization and sensitivity analysis presented in this paper also leads to different suggestions for future experimental design. In vivo experiments using the vapor uptake approach take advantage of a non-invasive approach to estimate 𝑉_max and 𝐾_𝑚 by using air chamber concentration. However, the current sensitivity analysis suggests an increased advantage to using in vitro approaches to estimate metabolic parameters. The sensitivity analysis performed consistently indicates that liver tissue measurements increase our ability to estimate 𝑉_max and 𝑘_GST for the two-pathway hypothesis. The same sensitivity analysis tools also indicate that our ability to estimate 𝑉_max1 and CL₂ for two-binding site hypothesis is similar to that for the two-pathway hypothesis. The sensitivity analysis for liver tissue concentration also suggests a peak or maximum ability to estimate 𝑉_max for an experiment taking 2 hours (instead of 6). Future experiments can include the possible combination of in vivo and in vitro vapor uptake techniques performed using intermediate concentrations just above and below metabolic saturation (500 ppm in this case).

Future in vivo research is needed to answer definitely the question as to whether BCM’s atypical metabolism can be described by a two different enzymes or by one enzyme with multiple binding sites. In this paper, we have used PBPK modeling to demonstrate the plausibility of both mechanisms in describing a linear kinetic pathway at higher concentrations. In order to determine if two separate enzymes are involved, a pharmacological agent can be used to deplete the GST pathway to determine if metabolism can proceed without GSH conjugation. As an example, such experiments have been performed for dichloromethane, using phorone as an pharmacological agent leading to GSH depletion [37]. These experiments were performed in vitro, using microsomes harvested from the GSH-depleted animals and used with vial equilibration techniques (to determine total changes in amount metabolized in vitro). Based on their results, the GST pathway was a small quantitative component of the total amount metabolized. Based on calculations provided by the current PBPK model, the GST pathway is expected to be less than 5% of the total amount metabolized. This small percent of amount metabolized difference is what is being described as either GST- or CYP2E1-dependent. In order to experimentally measure the difference in GST-mediated metabolism, a pharmacological agent such as phorone could be used for the next series of closed-chamber experiments and compared to BCM exposure without the phorone exposure. The difference between the two experiments would provide the actual difference in total metabolism observed between the two pathways.

As stated before, closed chamber inhalation experiments quantify changes in total chemical disappearance; thus, additional experiments are necessary to identify the actual P450 isoform involved in BCM metabolism. The confirmation of CYP2E1 as being the major P450 involved for in vivo BCM metabolism was obtained by using different P450 specific inhibitors and inducers to examine differences in CO production [32]. The use of inhibitors has become an important experimental tool to discern kinetic mechanisms involved in metabolism of different P450 isoforms. Pyrazole (CAS number 288-13-1), and its derivatives have been, used as CYP2E1 inhibitors by measuring total metabolic disappearance in microsomal fractions [38]. Recently, a structurally analog of pyrazole, 4-methylpyrazole or (4-MP), has been used to describe its CYP2E1 inhibition properties using a second-site kinetic mechanism [14]. These authors explained the inhibitory mechanism for 4-MP by adding a second-site to fit the experimental results using pNP (4-nitrophenol, CAS number 100-02-7) oxidation as a marker for CYP2E1 activity [39]. Future in vivo closed chamber experiments can make use of the in vitro paradigm used by [14] in that increasing pyrazole concentrations can be used to distinguish between the single-site or double-site models. These proposed experiments would add a dose response element to the single dose pyrazole inhibition experiments already performed by [9] as shown in Figure 8. A PBPK model can then be used to describe the resulting pyrazole inhibition as either containing a single or double site enzyme.

In summary, a combination of closed chamber data and PBPK modeling was used to examine different metabolic descriptions for BCM. The standard two-pathway description was compared to a two-binding site model within the same enzyme. Metabolic parameters for the different descriptions were optimized using the gas uptake data for BCM. Different metabolic parameters have different concentration ranges that determine their ability to be measured, and sensitivity analysis was used to demonstrate identifiability relationships between parameters. The benefits of these computational tools towards improved health risk assessments rely on their ability to accurately describe metabolic hypotheses.

5. MatLab Codes

There are three distinct sections to the code used to generate the results found in this paper. These sections are Data Transformations and Initialization, Compartmental Models, and Sensitivity Analysis. In addition to the brief descriptions below, each function has an extensive help section at the top of each 𝑚 file. Note that the names of the functions are italicized for clarity.

5.1. Data Transformations and Initialization

Almost all of the codes used for modeling require access to physiological parameters and the data set; the script file params writes these parameters into the working directory as well as converts the original data (in parts per million) to a form usable by the model (in mg/liter). Note that params is called before any other function as it generates a. mat file containing relevant parameters called by the other functions in the model.

5.2. Compartmental Models

The heart of the compartmental model is ode15s, a stiff ordinary differential equation (ODE) solver part of the standard MatLab package. There are several functions that build upon the output from the ODE solver. The main file, BCMmain, allows the user to choose to optimize for the metabolic parameters or to use the inputted values to plot the original data concentrations (in ppm) versus the model curves on a logarithmic scale.

BCMmain, when used strictly as a plotting tool, is self-contained (aside from calling the file generated by params and the ODE solver). In addition to a plot output, BCMmain also outputs a string containing the absolute root mean square error (rmse) and the relative rmse for the fit of each model at each concentration. The absolute error is calculated by comparing the two-norm of the vector containing the difference between the model predictions and the given data divided by the number of elements in the vector. The relative rmse error is given by the absolute error at each concentration divided by the initial concentration. This gives a sense of the percent error in the model.

When choosing to optimize for metabolic parameters, BCMmain calls on the function optthis, which uses a combination of fminsearch and a least-squares method to optimize for the parameters. The function fminsearch, available in the standard MatLab package, tries to minimize a “cost”; we define this cost to be the sum of the difference (obtained by linear-least squares) between the model predictions and data at each concentration. The equations for the PBPK model are contained in the file pbpk, which is called upon by the ODE solver and the curve-fitting functions. The output of optthis is the vector of metabolic parameters that yields the best fit (after a specified number of iterations) used to plot the model in BCMmain.

5.3. Sensitivity Analysis

There are two main functions used for the sensitivity analysis aspect of the code. The first function, plot2dsens takes in user-specified parameters and plots the sensitivity of these parameters in the chamber and the liver compartments. We choose these compartments due to their importance in obtaining data from the chamber as a noninvasive way; however data obtained from the liver may yield more insight. A second function, plot3dsens, and its associated helper function modelchoice, are designed to give a three-dimensional plot that shows the sensitivity of one parameter with respect to time as it varies over a user-input range in each compartment of interest. This analysis divides up the range of the parameter of interest into a user-specified number of subintervals and calculates the sensitivity of the parameter using the subinterval value as the parameter value in the same sense as the two-dimensional model. The sensitivity plots are then spliced together into a three-dimensional plot.

Both sensitivity functions require the use of tssolve, created by Adam Attarian (North Carolina State University, http://www4.ncsu.edu/~arattari/), which in turn requires Martin Fink’s myAD package. Both functions are available from the MathWorks website (http://www.mathworks.com).

5.4. Code Evaluation

One of the questions addressed in the present work is how to validate new code that is generated to study a chemical, such as the case for BCM. These authors approached this question on several fronts. First, a search was made for existing code using a similar chemical to BCM. The previously published DCM model [24] was selected for a comparison with BCM. The DCM code was converted to model BCM by changing chemical-specific parameters to reflect those of BCM. The simulations were performed using the BCM/DCM version. Results of these simulations were then compared with the newly created BCM code and found to be identical. Once the new BCM code was verified to be correct, the second step was to optimize for the parameters describing two-pathway kinetics and determine their similarity with previously published values (see Results). This additional step ensures that the new PBPK model description accurately depicts BCM metabolism.

Appendix

The differential equations used to describe each compartment in the PBPK model were developed using the law of conservation of mass: the amount of BCM flowing into each compartment must equal the amount of BCM leaving the compartment, whether it is through blood flow or ventilation, or through loss by metabolism or condensation. The following abbreviations can be found in Table 1:𝑑𝐴_ven𝑑𝑡=𝑞¹𝐶_rap+𝑞₂𝐶_adi+𝑞₃𝐶_slow+𝑞₄𝐶_kid+𝑞₅𝐶_liv⟹𝐶_ven−𝑄𝐶𝐶⋅𝐶_ven,(A.1)=𝑞¹𝐶_rap+𝑞₂𝐶_adi+𝑞₃𝐶_slow+𝑞₄𝐶_kid+𝑞₅𝐶_liv𝑑𝐴_lung𝑄𝐶𝐶,(A.2)𝑑𝑡=𝑄𝐶𝐶⋅𝐶^ven−𝐶_art+𝑄𝑉⋅𝐶_inh,⟹𝐶_art𝑄𝐶𝐶+𝑃𝑄𝑉_Bld,(A.3)=𝑄𝐶𝐶⋅𝐶^ven+𝑄𝑉⋅𝐶_inh𝑄𝐶𝐶+𝑄𝑉/𝑃_Bld𝐶,(A.4)_exh=𝐶^art𝑃_Bld𝑑𝐴_inh,(A.5)𝑑𝑡=𝑁⋅𝑄𝑉⋅𝐶^exh𝑑𝐴_rap−𝐶_inh−loss⋅𝐴_inh,(A.6)𝑑𝑡=𝑞¹𝑑𝐴_adi⋅𝐶_art−𝐶_rap,(A.7)𝑑𝑡=𝑞²𝑑𝐴_slow⋅𝐶_art−𝐶_adi,(A.8)𝑑𝑡=𝑞³𝑑𝐴_kid⋅𝐶_art−𝐶_slow,(A.9)𝑑𝑡=𝑞⁴𝑑𝐴_liv⋅𝐶_art−𝐶_kid−𝑀𝑒𝑡_𝐾,(A.10)𝑑𝑡=𝑞⁵⋅𝐶_art−𝐶_liv−𝑀𝑒𝑡_𝐿.(A.11)

For (A.1)–(A.11), we denote the venous concentration of BCM leaving a specific compartment comp as 𝐶_comp=𝐴_comp/(𝑉_comp⋅𝑃_comp), where 𝐴_comp is the amount of BCM in the compartment, 𝑉_comp is the volume of the compartment, and 𝑃_comp is the partition coefficient for the compartment. Note there is no associated partition coefficient for the chamber, thus the concentration in the chamber is denoted by 𝐶_inh=𝐴_inh/𝑉_cha. Partition coefficients (given in Table 1) describe the measure of the differential solubility of BCM as the chemical passes between two different phases; this unitless ratio accounts for the amount of BCM absorbed by compartmental tissue as it passes from the compartment back into the bloodstream. Previous studies have determined the tissue/air and blood/air partition coefficients for BCM [9]. To determine the tissue/blood partition coefficient for each compartment, (A.12) is used:𝑃_organ=tissue∶blood=tissue∶airblood∶air.(A.12) Several assumptions were made for the construction of each equation in the system: Initially, a differential equation, (A.1), was created to represent the amount of BCM in the venous compartment. However, simulations in MatLab demonstrated that the amount of BCM in the venous compartment approached equilibrium very quickly with respect to the length of the experiment. This allowed us to set (A.1) equal to zero to obtain an algebraic expression for the concentration of BCM in the venous compartment with respect to the concentration of BCM in the other compartments. A similar technique was used with the lung compartment (A.3), obtaining an algebraic expression for the concentration of BCM in the arterial compartment.

List of Abbreviations

Chemicals

4-MP:	4-methylpyrazole
BCM:	Bromochloromethane
BDCM:	Bromodichloromethane
DBP:	Disinfection by-products
DCM:	Dichloromethane
GST:	Glutathione transferase
PBPK:	Physiologically-based pharmacokinetic model
CCL:	Candidate Contaminant List, Office of Water, US EPA.

Equation Variables

𝑉max:	Maximum metabolic rate, mg/hour
𝑉max𝑖:	Maximum metabolic rate for binding site 𝑖 in two-binding site hypothesis, mg/hour
𝐾𝑚:	Affinity constant, mg/liter
𝐾𝑚𝑖:	Affinity constant for binding site 𝑖 in two-binding site hypothesis, mg/liter
𝑘GST:	Proportionality constant for linear pathway metabolized by GST, /hour
CL2:	clearance constant consisting of the ratio of the maximum metabolic rate and the affinity constant, 𝑉max2/𝐾𝑚2, for second binding site, liter/hour
𝐴compartment:	Amount of BCM in compartment
𝑃compartment:	Partition coefficient for compartment
𝑉compartment:	Volume of compartment, liters
𝐶compartment:	Concentration of BCM in compartment, given by 𝐶comp=𝐴comp/(𝑉comp⋅𝑃comp).

Disclaimer

This paper has been reviewed by the US Environmental Protection Agency and approved for publication. The views expressed in this paper are those of the authors and do not necessarily reflect the views or policies of the U S Environmental Protection Agency.

Acknowledgments

The authors would like to thank North Carolina State University and the US EPA for making this work possible. The following reviewers from the US EPA provided insightful comments: Drs. William Lefew, Rogelio Tornero-Velez, and Jane Ellen Simmons. Also a special thanks to the support from the NSF and NSA funding offices. The Research Experience for Undergraduates (REU) work was possible due to NSF Grant numbers DMS-0636568, DMS-0552571 and DMS-1063010 and NSA Grant number H9823-10-1-0252.