Application of Biologically Based Lumping To Investigate the

Publication Date (Web): February 18, 2016. Copyright This article not subject to ... Published 2016 by the American Chemical Society. *Phone: 919-541-...
0 downloads 0 Views 803KB Size
Download PDF
Article pubs.acs.org/est

Application of Biologically Based Lumping To Investigate the Toxicokinetic Interactions of a Complex Gasoline Mixture Micah N. Jasper, Sheppard A. Martin, Wendy M. Oshiro, Jermaine Ford, Philip J. Bushnell, and Hisham El-Masri* National Health and Environmental Effects Research Laboratory, Office of Research and Development, U.S. Environmental Protection Agency, Research Triangle Park, North Carolina 27709, United States S Supporting Information *

ABSTRACT: People are often exposed to complex mixtures of environmental chemicals such as gasoline, tobacco smoke, water contaminants, or food additives. We developed an approach that applies chemical lumping methods to complex mixtures, in this case gasoline, based on biologically relevant parameters used in physiologically based pharmacokinetic (PBPK) modeling. Inhalation exposures were performed with rats to evaluate the performance of our PBPK model and chemical lumping method. There were 109 chemicals identified and quantified in the vapor in the chamber. The time-course toxicokinetic profiles of 10 target chemicals were also determined from blood samples collected during and following the in vivo experiments. A general PBPK model was used to compare the experimental data to the simulated values of blood concentration for 10 target chemicals with various numbers of lumps, iteratively increasing from 0 to 99. Large reductions in simulation error were gained by incorporating enzymatic chemical interactions, in comparison to simulating the individual chemicals separately. The error was further reduced by lumping the 99 nontarget chemicals. The same biologically based lumping approach can be used to simplify any complex mixture with tens, hundreds, or thousands of constituents.

■

INTRODUCTION Predictive computational and statistical methods enable the investigation of the toxicological potential of complex mixtures, while minimizing the use of animals in research. The resulting inferences are often limited to the data ranges to which they are applied. Mechanistic models are based on descriptions of biological functions (absorption, distribution, metabolism, and excretion) and interactions between coexposed chemicals and with surrounding tissues. Combining mechanistic models with statistical methods has the potential to enable extrapolations outside the range of constraining data but within quantifiable limits of uncertainty for complex mixtures.1 Physiologically based pharmacokinetic (PBPK) models are a kind of mechanistic model that have been widely used to describe toxicokinetic interactions of simple and complex mixtures.2,3 One approach that utilizes PBPK modeling and the engineering concept of chemical lumping to quantitatively investigate the toxicology of complex mixtures has been described in an approach termed a Biologically Based Lumping Methodology (BBLM) by LeFew and El-Masri.4 Chemical lumping is defined as the process of grouping chemicals based on selected quantitative properties that are relevant to the question at hand (e.g., tissue dosimetry or toxicological response). Lumped groups are treated as pseudochemicals, with each lump having physiochemical properties intended to represent the properties of the individual chemicals that constitute the lump. Considering that a complex mixture This article not subject to U.S. Copyright. Published XXXX by the American Chemical Society

consists of a large number of chemical combinations and potential metabolic and/or toxicological interactions between the individual chemicals, a BBLM can greatly simplify the complexity of the mathematical model that describes the mixture. This simplification also reduces the computational power needed to mechanistically model the associated processes that influence absorption, distribution, metabolism, and excretion (ADME). Application of BBLM requires identification of physiochemical parameters by which lumps are generated. These parameters are obtained from a general PBPK model that captures the determinants of ADME that are relevant to chemical dosimetry and potential toxicokinetic interactions within a mixture. In a BBLM, a PBPK model provides the biological description that gives the relevant chemical specific parameters (such as partition coefficients or Michaelis−Menten metabolic parameters) that are used for the generation of lumps. Figure 1 provides a diagram of the inhalation PBPK model used in this research. Gasoline is a complex mixture with frequent human exposure, and it is used as an example in this research. Large human populations are exposed to gasoline liquid and vapor,5−9 and several individual components of gasoline are also known Received: November 23, 2015 Revised: February 17, 2016 Accepted: February 18, 2016

A

DOI: 10.1021/acs.est.5b05648 Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Environmental Science & Technology

Article

In the current work, we describe an approach for the application of chemical lumping methods to a complex mixture using physiochemical parameters that are described and quantified by PBPK models. In our approach, we use clustering algorithms, specifically the k-means clustering algorithm, to simplify the gasoline complex mixture (which contained 109 identified and quantified chemicals) by grouping chemicals into lumps. The process of simplifying the mixture with lumps results in “errors” that describe any resulting divergence between simulated values and experimental data (e.g., the measured concentrations of chemicals in blood). However, because our approach is based on computational methods, the errors generated by lumping can be quantified and then used to identify the ideal number of lumps to simplify the complex mixture and reduce the error within an acceptable level of tolerance compared to experimental error.

■

MATERIAL AND METHODS Overall Computational Approach. We followed the approach illustrated in Figure 2 to simplify the complex

Figure 1. Diagram of an inhalation PBPK model illustrating the tissue compartments. CInh is the concentration of inhaled chemicals. CArt is the concentration of the chemical in arterial blood. CVBrn, CVFat, CVSlw, CVRap, CVKid, and CVLiv refer to the venous blood concentrations of the chemical in the brain, fat/adipose tissue, slowly perfused tissue, rapidly perfused tissue, the kidneys, and the liver, respectively. QTot is the total cardiac output. QPul is the pulmonary ventilation rate. QBrn, QFat, QSlw, QRap, QKid, and QLiv refer to tissue blood flow to the brain, fat/adipose tissue, slowly perfused tissue, rapidly perfused tissue, the kidneys, and the liver, respectively. AMLiv refers to the amount metabolized in the liver.

animal and human toxicants.10,11 However, toxicity studies have not been performed on every identified constituent of gasoline and cannot possibly be performed on every combination of constituents within the mixture. Furthermore, potential interactions between mixture constituents influence the overall toxicity of the mixture.2 Specifically, the toxicity of the gasoline mixture is directly related to concentrations of individual chemicals in the target tissues (e.g., the brain, or the liver, etc.). These tissue concentrations are in turn a function of circulating concentrations in blood. Therefore, any factor that modulates the ADME of these chemicals will subsequently influence the resultant toxicity, in particular the toxicokinetic interactions between chemicals.2 Several PBPK models have been developed to quantitatively estimate the impact of toxicokinetic interactions among selected sets of chemicals within gasoline.12,13,4 Considering gasoline as a complex mixture, Dennison et al.14 applied engineering concepts of chemical lumping, together with individual PBPK models for a selected set of chemicals (benzene, toluene, ethylbenzene, xylene, and hexane), to simulate clearance of vapors during gas-uptake studies with rats. The lumped component of the gasoline mixture was described using a single set of chemical parameters (such as volatility-weighted fractions) that depended on the blend of gasoline. The resulting PBPK model of Dennison et al.14 successfully accounted for the binary interactions within the mixture by incorporating competitive metabolic inhibition between the lumped component and the five target chemicals. This type of modeling approach served as one of the first examples of how the engineering concept of chemical lumping can be used to investigate toxicokinetic interactions within complex mixtures.

Figure 2. Diagram illustrating the process of selecting the number of lumps (N) using the biologically based lumping methodology.

gasoline mixture into 10 target chemicals and sets of chemical lumps. First, we examined the individual components of the gasoline mixture to identify a set of target chemicals that should be excluded from lumping with others. These target chemicals are separated from the other chemicals in the mixtures because 1) they represent a large portion of the gasoline mixture by volume, 2) their concentrations in blood were specifically determined following in vivo experiments, and/or 3) they have been previously studied.12,13 The 10 target chemicals were npentane, 2-methylpentane, toluene, n-hexane, cyclohexane, benzene, m-xylene, ethylbenzene, p-xylene, and o-xylene. Their percentage in the gasoline sample and chemical properties are displayed in Supplemental Table 1. Second, we developed a general PBPK model for inhalation exposures in rats using relevant physiological and biochemical parameters. The general model was used to 1) identify chemical parameters for use in lumping the mixture, 2) simulate chemical concentrations in blood for selected individual target chemicals after considering metabolic interactions within the mixture, and 3) calculate errors due to simplification of the complex gasoline mixture, by comparing simulations of the selected chemicals to the corresponding experimental data. Third, chemical lumping was applied using k-means clustering. This method generates B

DOI: 10.1021/acs.est.5b05648 Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Environmental Science & Technology

Article

human and rat PC data has been reported to be less than the interlaboratory error for blood:air PCs,18 and we have found that the error associated with using human PC data, in lieu of rat data, for all tissue:air PCs is less than the error generated by using quantitative estimation methods for PCs (see Supplemental Tables 2 and 3). The PC values for chemicals where experimental data was not available (rat or human) were estimated using modified equations based on chemical boiling points (BPs).19 The BP method of Perbellini et al.19 was modified to take the maximum of the tissue:air and blood:air PC for each tissue, because for chemicals with low boiling points the original linear Perbellini et al.19 method would give negative PC values (see the Supporting Information for more details). The BP values were taken from online databases.20,21 When first tested against target chemical data, the PBPK simulations overestimated the data in most cases. Several parameters (e.g., cardiac output, ventilation rate, and/or blood:air partition coefficients) could have been adjusted to provide closer fits to data. However, of all the parameters investigated, the parameter that we had the least data to support and that was sensitive to blood data (see Supplemental Figure 1) was the lung:blood partition coefficient (PLng). This parameter was optimized using a least-squares method to fit the experimental data for the ten target chemicals when simulating the interaction effects of the complex mixture (all 109 identified, quantified chemicals). The Michaelis−Menten metabolism coefficients Vmax and Km are taken from the literature where available.22−25 A quantitative structure−property relationship (QSPR) method described by Price and Krishnan24 was used to estimate metabolic parameters where data was not available. Interaction Model. Metabolism in the liver is represented by Michaelis−Menten saturable kinetic equations. Interaction between chemicals is assumed to take place via enzymatic competitive inhibition. Competitive inhibition is a form of enzyme inhibition where binding of the one chemical to the active site on the enzyme prevents binding of another chemical. Therefore, one requirement for competitive inhibition to take place is that chemicals in a mixture are metabolized by the same enzyme. It is known that many of these chemicals are metabolized by the same enzyme (e.g., CYP2E1). Many of the volatiles in the gasoline mixture have been modeled previously using competitive inhibition mechanisms because it is the most common type of inhibition for the hydrocarbons in gasoline.26,14,2,27,28,3 The rate of the amount metabolized in the liver (RAMLiv) for a chemical using competitive inhibition is described according to the following equations

various numbers of lumps that were simulated with the ten target chemicals. An evaluation of error was performed incrementally for each number of lumps, in order to determine an applicable number of lumps to represent the complex mixture. PBPK Model Development. In a PBPK model, the human or animal body is divided into several compartments representing tissues and organs. The model enables generation of quantitative estimates for the distribution of chemicals within the body and the impact of any toxicokinetic interactions between chemicals. In order to simplify the model, tissues with similar properties are often grouped together based on the perfusion characteristics (i.e., rapidly or slowly perfused tissues). Figure 1 represents the PBPK model and the organs or tissues that were modeled. The lung compartment is where the arterial blood and the inhaled air equilibrate. The lung compartment can be thought of as a tissue group representing the lung, heart, and circulatory tissue. The venous blood is considered at quasi-steady-state, because it is assumed that the chemicals rapidly equilibrate within the venous blood. Finally, metabolism in the liver is represented by Michaelis−Menten saturable kinetic equations. The differential equations for the PBPK model are provided in the Supporting Information. PBPK Model Parameters. Physiological parameters describing bodyweight (BW), tissue volumes (V), pulmonary/alveolar ventilation rate (QPul), total blood circulation rate or cardiac output (QTot), and tissue blood flow rates (Q) were obtained from the literature15 and are presented in Table 1. The cardiac output and maximal metabolic rate (Vmax) were Table 1. Physiological Parameters for the Rata parameter volumes BW VLng VBrn VKid VLiv VFat VSlw VRap flow rates QPul QTot QBrn QKid QLiv QFat QSlw QRap a

value

unit

total body weight lungs brain kidneys liver fat/adipose tissue slowly perfused tissues rapidly perfused tissues

0.25 0.0050 0.0057 0.0073 0.0366 0.07 0.5818 0.2986

kg *BW *BW *BW *BW *BW *BW *BW

pulmonary ventilation total cardiac output brain kidneys liver fat/adipose tissue slowly perfused tissues rapidly perfused tissues

60 16.5 0.02 0.141 0.174 0.07 0.159 0.580

*BW (L/h/kg) *BW0.75 (L/h/kg0.75) *QTot *QTot *QTot *QTot *QTot *QTot

RAMLivi =

15

Values taken from Brown et al.

Vmax i ∗CLivi dAMLivi = dt K m*i + CLivi

⎛ K m*i = K mi⎜⎜1 + ⎝

allometrically scaled to the body weight by a factor of 0.75 (BW0.75), and the pulmonary ventilation rate (QPul) was linearly scaled to body weight.15 The value for pulmonary ventilation is higher than the mean of the experimental data reported but well within the range of the data. Partition coefficient (PC) values for rat tissues were taken from literature where available.16,17 If experimental PC data for rats were not available, but experimental data for humans were, the human PC values were used. The difference between

∑ j

CLivj ⎞ ⎟ K ij ⎟⎠

(1)

(2)

where i is the current chemical, and j is one of all of the other chemicals in the mixture, RAMLiv is the rate of the amount metabolized in the liver, Vmax is maximum rate of metabolism in the liver, CLiv is the concentration in the liver, Km is the Michaelis−Menten affinity constant, Ki is the competitive inhibition constant, and Km* is the apparent Michaelis− Menten constant with competitive inhibition. Even though C

DOI: 10.1021/acs.est.5b05648 Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Environmental Science & Technology

Article

scales or normalizes the simulated values in relation to the average (avg), maximum (max), and minimum (min) experimental values for each time point (t), where n is the total number of time points. The RMSNE acts as an indicator of the goodness of fit and can be compared across the different target chemicals.

there may be qualitative and quantitative differences between the Ki and Km values, Dennison et al.12 show that Km and Ki are of the same order of magnitude, and Price and Krishnan24 show that Km and Ki are approximately equal so that Ki can be replaced by Km in PBPK simulations. k-Means Clustering. The k-means clustering algorithm was used to lump the nontarget chemicals into a specified number of lumps. For a specified number of clusters, k, the k-means algorithm determines centers (i.e., means) such that the sum of square distances from each point to the center of its cluster is minimized. For the purpose of this research, the clusters are the chemical lumps, the points are the chemical specific parameters for each individual chemical, and the centers become the pseudochemicals that represent each lump. The centers of the clusters returned by the k-means algorithm are used as the chemical lumps in the PBPK model. The k-means algorithm typically uses an iterative method called Lloyd’s algorithm to cluster points and converge on a local minimum for the sum of square distances of each point in a cluster to the center of the cluster.29 We used the k-means clustering function in Matlab (2014a, Mathworks) to lump the chemicals in the gasoline complex mixture. The chemicals in the gasoline mixture were clustered using the k-means algorithm based on all of the chemical specific parameters as they are used in the PBPK model, specifically the maximum rate of metabolism (Vmax), the Michaelis−Menten affinity constant (Km), the blood:air partition coefficient (PBA), and the tissue:blood partition coefficients for the lung, fat, liver, muscle, rapidly perfused tissues, slowly perfused tissues, kidney, and brain (PLng, PFat, PLiv, PMsl, PRap, PSlw, PKid, PBrn). A relative sensitivity analysis was performed to determine the impact of chemical specific parameters on the amount of chemical in the arterial blood. The results are shown in Supplemental Figure 1. It is possible that not all of the chemical specific parameters have to be used when lumping. If the number of chemical specific parameters is very large, one could select only the most sensitive parameters to cluster the chemicals lumps and then average the other nonsensitive parameters. When the chemicals are lumped, the generated lump has its own unique properties. The composite dose for each lump in the PBPK model is a sum of all of the doses for the individual chemicals in the lump, and its physiochemical properties are assigned based on the k-means clustering algorithm. In general, the physiochemical properties of a lump can be thought of as an average of the properties of the chemicals in the lump. Once each lump is defined, its properties and dose are included in the PBPK model to simulate its competitive inhibition interaction effects on the blood levels of the target chemicals. Because the PBPK model we developed is general and dynamic, it is capable of simulating any number of chemicals and/or any number of lumps. Error Calculation. In order to evaluate the impact of simplifying the complex mixture into a number of lumps we calculated the error introduced by the lumping process. This error was calculated by comparing the experimental data to the simulated values for blood concentration for each of the 10 target chemicals in blood. In this research, the error calculation is performed on one target chemical at a time. The root mean squared normalized error (RMSNE) of simulated values (xsim) to experimental data (xexp) was used as the error function. The RMSNE equation is similar to the more commonly used normalized root mean squared error (NRMSE) except that it

RMSNE =

⎛ xsim, t − xexp,avg, t ⎞2 ⎟ /n ∑ ⎜⎜ x − xexp,min, t ⎟⎠ t ⎝ exp,max, t

(3)

Animals. Adult female nonpregnant (n = 24) Long-Evans rats (Charles River Laboratories, Raleigh, NC) were maintained on a 12-h light-dark cycle (lights on at 6:00 a.m, temperature 22° ± 2 °C, and relative humidity 50 ± 10%). Rats were maintained in groups of two per polycarbonate cage (18780, Maryland Plastics, Inc., Federalsburg, MD). Each cage contained kiln-dried pine shaving bedding (Northeastern Products, Warrensburg NY) that was changed twice weekly by animal care staff. Ad libitum food (PMI 5001, Lab Diet/PMI Nutrition International, Richmond, IN) and water were available in the animal facility. Animal colonies at the US EPA (Research Triangle Park, NC) are fully accredited by the Association for Assessment and Accreditation of Laboratory Animal Care International according to NIH guidelines. Research protocols were approved by the US EPA, Office of Research and Development, National Health and Environmental Laboratory Institutional Animal Care and Use Committee. Exposure System. The vapor condensates of gasoline were generated for this project by Chevron USA, Inc. (Richmond, CA) under contract with the EPA as described by Bushnell et al.30 Rats were exposed using a nose-only chamber described by Martin et al.31 Briefly, a 52-port Cannon-style directed flow N− O exposure system was used to expose rats and enable rapid removal for sacrifice and collection of kinetic time points. Liquid vapor condensates were first expressed from pressurized storage containers into a heated J-tube, according to Bushnell et al.30 Once volatilized, the vapors were passed into the inner plenum of the exposure chamber and expressed into the breathing zones of the rats via jets within the N−O ports. Exhaled and waste vapor was removed from the breathing zone via a negative pressure vacuum flow in the outer plenum of the chamber. Vapor concentrations were measured according to Bushnell et al.30 In the few situations where chemicals coeluted together, we quantitatively separated the chemicals. The ratio of m- to pxylene was estimated from data collected by GC-FID to be 66.5% and 33.5% of the measured dose for the pair. The three other coeluted pairs (3-ethylpentane and trans-1,3-dimethylcyclopentane; 2,4-dimethylhexane and ethylcyclopentane; 3,3dimethylhexane and 1,2,4-trimethylcyclopentane) were assumed to be equally divided and split 50%−50% of the measured dose for each pair. The rats were exposed to 9000 ppmv vapor condensate for a maximum of 6 h. The percent mass by carbon (percent) and the volume fraction (ppmv) for each of the 109 identified and quantified chemicals in the gasoline vapor are listed in Supplemental Table 4. Analysis of Blood Samples. The blood and tissue samples were collected, prepared, and stored as previously described for nonpregnant rats.31 The complexities of analyzing a fuel mixture in blood necessitated that a subgroup of chemicals be D

DOI: 10.1021/acs.est.5b05648 Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Environmental Science & Technology

Article

Table 2. Chemical Constituents of the 15 Lumps Using k-Means Clustering To Group the Nontarget Chemicals parameter 10 target chemicals Lump1 Lump2 Lump3 Lump4 Lump5 Lump6 Lump7 Lump8 Lump9 Lump10 Lump11 Lump12

Lump13 Lump14 Lump15

remarks n-pentane; 2-methylpentane; toluene; n-hexane; cyclohexane; benzene; m-xylene; ethylbenzene; p-xylene; o-xylene n-butane; cyclopentane; n-octane; 1,3,5-trimethylbenzene; 1,2,3-trimethylbenzene propane; propylene 2,5-dimethylhexane; 3-methylheptane; 2,4-dimethylhexane; 2-methylheptane; 2,3-dimethylhexane; 4-methylheptane; 3,4-dimethylhexane; 3,3dimethylhexane; 2-methyl-3-ethylpentane 2-methyl-2-pentene; trans-2-hexene; trans-3-methyl-2-pentene; cis-3-methyl-2-pentene; 2-methyl-1-pentene; cis-3-hexene; 1-methylcyclopentene; cis-2hexene; 3-methyl-1-pentene; 1-hexene ctc-1,2,4-trimethylcyclopentane; 2,4,4-trimethyl-1-pentene; 1,2,4-trimethylcyclopentane; ctc-1,2,3-trimethylcyclopentane; trans-1,4-dimethylcyclohexane; cis-2-octene; trans-2-octene; isopropylcyclopentane; cis-1,2-dimethylcyclohexane; n-propylcyclopentane isobutane; trans-2-butene; cis-2-butene; isobutylene; 1-butene 2,2,4-trimethylpentane isopentane; 2-methyl-2-butene; trans-2-pentene; 2-methyl-1-butene; cis-2-pentene; 1-pentene; cyclopentene; 3-methyl-1-butene 3-methylpentane; 2,2-dimethylbutane ethanol; methanol 1,2,4-trimethylbenzene; m-ethyltoluene; p-ethyltoluene; n-propylbenzene; o-ethyltoluene; indane 2,3-dimethylbutane; 2,4-dimethylpentane; 2,3-dimethylpentane; 2-methylhexane; methylcyclohexane; cis-1,3-dimethylcyclopentane; ethylcyclopentane; 3-ethylpentane; trans-1,3-dimethylcyclopentane; 3-methyl-cis-3-hexene; 3,3-dimethylpentane; 2,2-dimethylpentane; 1,1-dimethylcyclopentane; 2,2,3trimethylbutane; trans-3-heptene; trans-3-methyl-2-hexene; 2,3-dimethyl-2-pentene; cis-4-methyl-2-hexene; 4-methyl-1-hexene; 2,4-dimethyl-1pentene n-heptane; 2,3,4-trimethylpentane; 2,2-dimethylpropane methylcyclopentane; 3-methylhexane 2,2,5-trimethylhexane; 2,3,5-trimethylhexane; 2,5-dimethylheptane; 2,2,4-trimethylhexane; 3-methyloctane; 2,4-dimethylheptane; ctt-1,2,4trimethylcyclohexane; 2,6-dimethylheptane; n-nonane; 1,1,4-trimethylcyclohexane; 2,2-dimethylheptane; 4,4-dimethylheptane; n-decane; 3methylnonane

and quantified chemicals, which constitute more than 99.7% of the gasoline mixture by volume, were separated into 10 target chemicals and 99 nontarget chemicals. The nontarget chemicals were lumped into clusters based on all chemical-specific parameters (parameters given in Supplemental Table 4). To illustrate the lumping method, we plotted three parameters (Vmax, Km, and PBA) for 5 lumps and 10 target chemicals as shown in Supplemental Figure 2. Even though three parameters are plotted, the lumps and centers were generated by the kmeans clustering algorithm using all chemical specific parameters. For each chemical lump, the “center” results in chemical specific parameters that represent the chemical cluster (shown as black crosses in the center of the lumps). For this example, the chemicals that constitute each lump are shown in Supplemental Table 6, and the chemical parameters for the five chemical lumps are shown in Supplemental Table 7. In this research, clustering was performed using varying numbers of lumps (N) and all chemical specific PBPK parameters for the 99 nontarget chemicals in the complex mixture. Then, these lumps were simulated as interacting with the 10 target chemicals using the competitive inhibition interaction model. For each number of lumps (N) an error was calculated to illustrate the divergence between simulated values and experimental data of the blood concentrations for each target chemical when interacting with the lumps and the remaining target chemicals. Optimizing the Number of Lumps. Using experimental data, errors were calculated to determine the goodness of fit of between PBPK simulations and in vivo data. The error is calculated using the NRMSE equation (eq 3). Simulations were conducted using the 10 target chemicals with varying numbers of lumps (N) until the maximum number of lumps equaled all the number of nontarget chemicals (99 lumps). The 99 lump case is equivalent to simulating all of the 109 identified chemicals in the gasoline sample. The results for the influence of the number of lumps on error are illustrated in Supplemental Figure 3. All of the individual target chemicals were simulated

tracked, rather than attempting to quantify every measurable chemical in the mixture. Approximately 20 chemicals were quantified in blood samples, of which 10 were used as target chemicals, for reasons described earlier in the Materials and Methods. For quantification, hydrocarbon standards (SigmaAldrich Inc., St. Louis, MO) were used to develop standard curves for benzene, toluene, ethylbenzene, o-xylene, m-xylene, p-xylene, n-pentane, n-hexane, 2-methylpentane, and cyclohexane. A neat-liquid gasoline standard and four separate internal standards for n-pentane, n-hexane, toluene, and benzene were also used (Sigma-Aldrich Inc., St. Louis, MO). The solution of standards described above was dissolved in cumene (4-isopropylbenzene). A “standard mix” was recreated every 14 days during the analysis to maintain the desired concentrations of these volatile hydrocarbons. Blood was analyzed using gas chromatography/mass spectrometry. In order to improve the limits of detection, dynamic headspace was used to concentrate the hydrocarbons onto an analytical trap. The limit of detection for hydrocarbons in blood was approximately 1.0 ng/g. Prior to analysis, blood samples were removed from −20 °C storage, allowed to briefly equilibrate on ice before pipetting into headspace vials. These aliquots of blood were added to 20 mL headspace vials containing deuterated standards, capped, and then immediately placed on the auto sampler. Samples were analyzed with matrix blood blanks and laboratory control duplicates as a quality control step. Control limits were ±30% for spiked samples, with a relative percent difference for duplicate samples of 20%. The experimental blood concentrations for the target chemicals are given in Supplemental Table 5.

■

RESULTS Generation of Chemical Lumps. The use of the k-means clustering algorithm gives the flexibility to specify any number of chemical lumps based on PBPK parameters, up to the total number of nontarget chemicals. For instance, the 109 identified E

DOI: 10.1021/acs.est.5b05648 Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Environmental Science & Technology

Article

as interacting both with each other and with the N lumps of nontarget chemicals. As Supplemental Figure 3 shows, increasing the number of lumps decreased the error; however, the error eventually began to plateau, and improvements due to increasing the number of lumps diminished at around 15 lumps. The chemicals that constitute each of these 15 lumps are shown in Table 2 and illustrated in Figure 3, and the chemical

alone (no interaction) (solid red line). It is also apparent that there was a minimal difference between simulating all 109 interacting chemicals (solid black line) and simulating a simplified model with 10 target chemicals and 15 lumps (dashed blue line). In fact, these simulated curves approximately superimposed each other. Furthermore, the improvement was more significant for benzene, ethylbenzene, toluene, and m-, p-, and o-xylene. The fact that lumping significantly improved the simulation error for some chemicals but not as much for others indicates that the effects of competitive inhibition are more impactful for chemicals that have a higher rate of metabolism (i.e., high Vmax and/or low Km). In several cases, however, improvements in fit to early time points and peak concentrations coincide with a decreased quality of fits to clearance data in comparison to the simulations for individual chemicals and vice versa. In each such case, the general rate of clearance was maintained, though the actual data during the clearance phase or the dosing phase were somewhat overpredicted. Finally, for some chemicals, the 10 target chemicals and 1 lump and the 10 target chemicals and 15 lumps simulations predict higher values (instead of lower values) than the simulation of all 109 chemicals; however, the error (calculated by the RMSNE error equation given previously) is still higher for these simulations.

■

DISCUSSION Experimentally, several efforts have investigated the toxicity of complex mixtures as one lumped set of chemicals.12,13 Although valuable information can be obtained from these experiments, they do not provide insight into the impact of the mixture chemical composition on any possible interactions between the mixture’s constituents. This detailed information is critical, since composition of complex mixtures can differ depending on source and the method of exposure. Quantitative and statistical methods were developed to investigate toxicological interactions of complex mixtures, in lieu of performing additional resource intensive in vivo experiments. Considering the large number of possible permutations and combinations for modeling interactions between individual components, research has been focused on simplifying complex mixtures into interacting lumps. The Biologically Based Lumping Methodology (BBLM) in this research provides a systematic approach for deciding on the number of lumps and estimating the impact of this simplification on the overall toxicological effects of the mixture, both as a whole and for any selected individual components. This approach is based on aggregating chemicals based on their biological impact (such as parameters impacting ADME) and the influence of toxicological interactions (such as enzymatic metabolism inhibition) on the ability of chemicals in the mixture to distribute to target tissues (e.g., the brain). A joint computational method for lumping (e.g., k-means or Ward clustering algorithms) and biological models to identify chemical parameters for the lumping and interaction mechanism (e.g., PBPK) were both utilized in our BBLM. One advantage of the lumping process is the simplification of complex mixtures to a level where the quantitative analysis of interactions can be effectively performed. Although increasing the number of lumps, and chemicals simulated in general, improves the accuracy of the results, the time and effort to run the simulation increases as the number of chemicals being simulated increases. Supplemental Figure 5 illustrates this increase in computational time as a function of number of

Figure 3. Plot of 15 lumps clustered from all 109 identified and quantified chemicals. Target chemicals are plotted with vertical crosses +, nontarget chemicals with dots •, the center of each lump with diagonal crosses ×. The nontarget chemicals that constitute the same lump are connected with a line. Each lump is plotted with a different color, circled, and labeled. The axes are the blood/air partition coefficient (PBA), the Michaelis−Menten metabolism term Km, and the VmaxC (Vmax is later allometrically scaled to body weight in the PBPK model; however, the unscaled VmaxC is shown in the figure.).

parameters for the 15 chemical lumps are given in Supplemental Table 8. The number of lumps to use can be chosen based on the incremental improvement in average error from the previous lump configuration. In other words, increase the number of lumps until the percentage improvement from one configuration to the next is below the predetermined threshold. It can be seen in Supplemental Figure 3 that some chemicals are more sensitive to increasing the number of lumps. This indicated that inhibition affects the metabolism for some chemicals more than others. Finally, the discrete noise in the error plots is due to the clustering algorithm, because the lumping configuration differs slightly for each number of lumps N. Comparison with Data. For each of the target chemicals, simulations were performed to evaluate the effects of chemical lumping. Supplemental Figure 4 shows the simulated arterial blood concentrations for each target chemical, plotted against experimental data using three scenarios: 1) each of the target chemicals were simulated individually (no interactions), 2) all of the 10 target chemicals were simulated interacting with 15 lumps of the remaining chemicals in the mixtures, and 3) all 109 of the identified and quantified chemicals in the mixture were simulated interacting with each other for comparison with the rest of the simulations. It is clear from Supplemental Figure 4 that the difference between the simulated values and experimental data is drastically reduced by considering the interaction with the remaining chemicals in the mixture (dashed blue line), instead of simulating each target chemical F

DOI: 10.1021/acs.est.5b05648 Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Environmental Science & Technology

Article

lumped chemicals. This computational time increase may not be substantial if all 109 chemicals are simulated for only eight simulation hours, as in this case. However, for longer simulation periods of weeks or years, as in life-stage models, the computational time may considerably increase. Additionally, the plot of computational time does not show the human effort that is required to collect data from either experiments or literature, as is needed to perform detailed investigations of complex mixtures. By determining the PK of a set of target chemicals, in combination with an appropriate number of lumps (using an acceptable tolerance for error compared to experimental error) will ultimately save time, effort, and resources. We have demonstrated that the error decreases with an increasing number of lumps. However, as expected, the improvement obtained by using more lumps gradually decreases to a point where there is little or no improvement gained by increasing the number of lumps. In this research, this point of diminishing improvement was reached at around 15 lumps. The chemicals that constitute each of the 15 lumps are shown in Table 2, and the chemical parameters for the 15 chemical lumps are shown in Supplemental Table 8. There are two large gains in simulation accuracy using interaction and lumping, compared to simulating target chemicals individually: 1) the accuracy drastically improves by simulating multiple interacting chemicals (e.g., the 10 interacting target chemicals with no lumps) instead of a single chemical with no interaction, around a 36% improvement in some cases. This is similar to the approach of Haddad et al.13 2) There is another large improvement (around 1%) by including just a single lump of the nontarget chemicals (e.g., the 10 target chemicals with 1 lump). This is similar to the approach of ref 12. The simulation accuracy will continue to improve by incorporating more lumps, but the gains are small compared to these first two leaps, and the gains continue to diminish with increasing numbers of lumps. These gains in accuracy are illustrated in Figure 4, comparing the simulation results for various lump configurations with experimental data for benzene. The improvements seen by including chemical interaction and lumping are dependent on the metabolism of the chemical itself and its corresponding Vmax and Km values. If a chemical is metabolized slowly (a low Vmax and a high Km), then almost no improvement is gained by including interaction and by including lumping. This effect is seen with n-pentane, 2methylpentane, n-hexane, and cyclohexane. However, if the chemical is metabolized quickly (a high Vmax and a low Km), then much improvement is gained by including interaction and lumps. This is seen with chemicals toluene, benzene, ethylbenzene, and m-, p-, and o-xylene. In conclusion, the Biologically Based Lumping Method (BBLM) presented in this research is one of the first examples of applying clustering algorithms and biologically based chemical lumping to a complex mixture. In BBLM, generic PBPK models provide a biological basis for considering parameters for the clustering algorithm. PBPK models provide a quantitative basis to include relevant routes of exposure, species extrapolation, and the impact of low dose concentrations on the presence or absence of interactions.32 In this example we applied BBLM to a complex gasoline mixture consisting of 109 identified and quantified chemicals. It is shown that there was a minimal difference between simulating a simplified mixture with 10 target chemicals and 15 lumps and simulating all 109 interacting chemicals. Though gasoline was

Figure 4. Plots of the mean experimental data (crosses ×) and the simulated values (lines) for the arterial blood concentrations for benzene when simulating 1) the single chemical individually with no interaction (bottom solid line), 2) 10 target chemicals interacting (no lumps), 3) the 10 target chemicals and 1 lump, 4) the 10 target chemicals and 15 lumps, and 5) all 109 identified and quantified chemicals interacting. The error bars show the maximum and minimum values of the experimental data.

specifically considered in this research, the lumping method presented can be applied to any kind of complex mixture to lump the chemicals and simplify the mixture based on biologically relevant chemical properties.

■

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.est.5b05648.

■

8 tables, 5 figures, and 11 PBPK model equations (PDF)

AUTHOR INFORMATION

Corresponding Author

*Phone: 919-541-1302. Fax: 919-541-428. E-mail: el-masri. [email protected]. Corresponding author address: US EPA/ NHEERL/ISTD/SBB, 109 T.W. Alexander Drive, B105-03, Research Triangle Park, NC 27709. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This manuscript has been reviewed in accordance with the policy of the National Health and Environmental Effects Research Laboratory, U.S. Environmental Protection Agency, and approved for publication. Approval does not signify that the contents necessarily reflect the views and policies of the Agency, nor does mention of trade names or commercial products constitute endorsement or recommendation for use. We would like to thank Drs. Marina Evans and Cecelia Tan for their valuable comments and input while reviewing this paper. G

DOI: 10.1021/acs.est.5b05648 Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Environmental Science & Technology

■

Article

L., Jr. Development of multi-route physiologically-based pharmacokinetic models for ethanol in the adult, pregnant, and neonatal rat. Inhalation Toxicol. 2012, 24 (11), 698−722. (23) Pastino, G. M.; Conolly, R. B. Application of a physiologically based pharmacokinetic model to estimate the bioavailability of ethanol in male rats: distinction between gastric and hepatic pathways of metabolic clearance. Toxicol. Sci. 2000, 55 (2), 256−265. (24) Price, K.; Krishnan, K. An integrated QSAR-PBPK modelling approach for predicting the inhalation toxicokinetics of mixtures of volatile organic chemicals in the rat. SAR and QSAR in environmental research 2011, 22 (1−2), 107−28. (25) Tephly, T. R.; Parks, R. E., Jr.; Mannering, G. J. Methanol Metabolism in the Rat. J. Pharmacol. Exp. Ther. 1964, 143, 292−300. (26) Campbell, J. L., Jr.; Fisher, J. W. A PBPK modeling assessment of the competitive metabolic interactions of JP-8 vapor with two constituents, m-xylene and ethylbenzene. Inhalation Toxicol. 2007, 19 (3), 265−73. (27) Haddad, S.; Tardif, R.; Charest-Tardif, G.; Krishnan, K. Physiological modeling of the toxicokinetic interactions in a quaternary mixture of aromatic hydrocarbons. Toxicol. Appl. Pharmacol. 1999, 161 (3), 249−57. (28) Martin, S. A.; Campbell, J. L.; Tremblay, R. T.; Fisher, J. W. Development of a physiologically based pharmacokinetic model for inhalation of jet fuels in the rat. Inhalation Toxicol. 2012, 24 (1), 1−26. (29) Lloyd, S. P. Least-Squares Quantization in Pcm. IEEE Trans. Inf. Theory 1982, 28 (2), 129−137. (30) Bushnell, P. J.; Beasley, T. E.; Evansky, P. A.; Martin, S. A.; McDaniel, K. L.; Moser, V. C.; Luebke, R. W.; Norwood, J., Jr.; Copeland, C. B.; Kleindienst, T. E.; Lonneman, W. A.; Rogers, J. M. Toxicological assessments of rats exposed prenatally to inhaled vapors of gasoline and gasoline-ethanol blends. Neurotoxicol. Teratol. 2015, 49, 19−30. (31) Martin, S. A.; Oshiro, W. M.; Evansky, P. A.; Degn, L. L.; Ledbetter, A. D.; Ford, J.; Krantz, Q. T.; LeFew, W. R.; Beasley, T. E.; El-Masri, H.; McLanahan, E. D.; Boyes, W. K.; Bushnell, P. J. Use of novel inhalation kinetic studies to refine physiologically-based pharmacokinetic models for ethanol in non-pregnant and pregnant rats. Inhalation Toxicol. 2014, 26 (10), 598−619. (32) el-Masri, H. A.; Constan, A. A.; Ramsdell, H. S.; Yang, R. S. Physiologically based pharmacodynamic modeling of an interaction threshold between trichloroethylene and 1,1-dichloroethylene in Fischer 344 rats. Toxicol. Appl. Pharmacol. 1996, 141 (1), 124−32.

REFERENCES

(1) El-Masri, H. A. Experimental and mathematical modeling methods for the investigation of toxicological interactions. Toxicol. Appl. Pharmacol. 2007, 223 (2), 148−54. (2) Haddad, S.; Beliveau, M.; Tardif, R.; Krishnan, K. A PBPK modeling-based approach to account for interactions in the health risk assessment of chemical mixtures. Toxicol. Sci. 2001, 63 (1), 125−131. (3) Yang, R. S.; El-Masri, H. A.; Thomas, R. S.; Dobrev, I. D.; Dennison, J. E., Jr.; Bae, D. S.; Campain, J. A.; Liao, K. H.; Reisfeld, B.; Andersen, M. E.; Mumtaz, M. Chemical mixture toxicology: from descriptive to mechanistic, and going on to in silico toxicology. Environ. Toxicol. Pharmacol. 2004, 18 (2), 65−81. (4) LeFew, W.; El-Masri, H. Computational estimation of errors generated by lumping of physiologically-based pharmacokinetic (PBPK) interaction models of inhaled complex chemical mixtures. Inhalation Toxicol. 2012, 24 (1), 36−46. (5) Claxton, L. D. The history, genotoxicity, and carcinogenicity of carbon-based fuels and their emissions. Part 3: diesel and gasoline. Mutat. Res., Rev. Mutat. Res. 2015, 763, 30−85. (6) Dement, J. M.; Hensley, L.; Gitelman, A. Carcinogenicity of gasoline: a review of epidemiological evidence. Ann. N. Y. Acad. Sci. 1997, 837, 53−76. (7) Hakkola, M.; Saarinen, L. Exposure of tanker drivers to gasoline and some of its components. Ann. Occup. Hyg. 1996, 40 (1), 1−10. (8) Hakkola, M. A.; Saarinen, L. H. Customer exposure to gasoline vapors during refueling at service stations. Appl. Occup. Environ. Hyg. 2000, 15 (9), 677−680. (9) Wixtrom, R. N.; Brown, S. L. Individual and population exposures to gasoline. J. Exposure Anal. Environ. Epidemiol. 1992, 2 (1), 23−78. (10) Graham, D. G.; Amarnath, V.; Valentine, W. M.; Pyle, S. J.; Anthony, D. C. Pathogenetic Studies of Hexane and Carbon-Disulfide Neurotoxicity. Crit. Rev. Toxicol. 1995, 25 (2), 91−112. (11) Laitinen, J.; Kangas, J.; Pekari, K.; Liesivuori, J. Short-Time Exposure to Benzene and Gasoline at Garages. Chemosphere 1994, 28 (1), 197−205. (12) Dennison, J. E.; Andersen, M. E.; Clewell, H. J.; Yang, R. S. Development of a physiologically based pharmacokinetic model for volatile fractions of gasoline using chemical lumping analysis. Environ. Sci. Technol. 2004, 38 (21), 5674−81. (13) Haddad, S.; Charest-Tardif, G.; Krishnan, K. Physiologically based modeling of the maximal effect of metabolic interactions on the kinetics of components of complex chemical mixtures. J. Toxicol. Environ. Health, Part A 2000, 61 (3), 209−223. (14) Dennison, J. E.; Andersen, M. E.; Yang, R. S. Characterization of the pharmacokinetics of gasoline using PBPK modeling with a complex mixtures chemical lumping approach. Inhalation Toxicol. 2003, 15 (10), 961−986. (15) Brown, R. P.; Delp, M. D.; Lindstedt, S. L.; Rhomberg, L. R.; Beliles, R. P. Physiological parameter values for physiologically based pharmacokinetic models. Toxicol. Ind. Health 1997, 13 (4), 407−84. (16) Gargas, M. L.; Burgess, R. J.; Voisard, D. E.; Cason, G. H.; Andersen, M. E. Partition coefficients of low-molecular-weight volatile chemicals in various liquids and tissues. Toxicol. Appl. Pharmacol. 1989, 98 (1), 87−99. (17) Meulenberg, C. J.; Vijverberg, H. P. Empirical relations predicting human and rat tissue:air partition coefficients of volatile organic compounds. Toxicol. Appl. Pharmacol. 2000, 165 (3), 206−16. (18) Abraham, M. H.; Ibrahim, A.; Acree, W. E., Jr. Air to blood distribution of volatile organic compounds: a linear free energy analysis. Chem. Res. Toxicol. 2005, 18 (5), 904−11. (19) Perbellini, L.; Brugnone, F.; Caretta, D.; Maranelli, G. Partition coefficients of some industrial aliphatic hydrocarbons (C5-C7) in blood and human tissues. Occup. Environ. Med. 1985, 42 (3), 162−7. (20) Pence, H. E.; Williams, A. ChemSpider: An Online Chemical Information Resource. J. Chem. Educ. 2010, 87 (11), 1123−1124. (21) USEPA EPISuite. http://www.epa.gov/opptintr/exposure/ pubs/episuite.htm (accessed June 2015). (22) Martin, S. A.; McLanahan, E. D.; El-Masri, H.; LeFew, W. R.; Bushnell, P. J.; Boyes, W. K.; Choi, K.; Clewell, H. J., 3rd; Campbell, J. H

DOI: 10.1021/acs.est.5b05648 Environ. Sci. Technol. XXXX, XXX, XXX−XXX