---

This article is part of the monograph Application of Technology to Chemical Mixture Research.

Address correspondence to R.S.H. Yang, Center for Environmental Toxicology and Technology, Colorado State University, Foothills Campus, Fort Collins, CO 80523 USA. Telephone: (970) 491-5652. Fax: (970) 491-8304. E-mail: raymond.yang@colostate.edu

The research work and related concept development on chemical mixtures were supported in part by Superfund Basic Research Program Project grant P42 ES05949 from the National Institute of Environmental Health Sciences (NIEHS), cooperative agreement U61/ATU881475 from the Agency for Toxic Substances and Disease Registry, NIEHS Quantitative Toxicology training grant 1 T32 ES 07321, and research grant RO1 ES-09655 from NIEHS. Without such generous support for biomedical research, this work could never have been possible.

Received 18 December 2001; accepted 27 September 2002.

One question frequently asked in the area of toxicology of chemical mixtures is, "How does one deal with a complex chemical mixture?" Here, when we discuss complex chemical mixtures, we are referring to something such as the smoke from the burning oil fields in Kuwait during the Persian Gulf War. That smoke not only contains hundreds or even thousands of chemicals but also has the characteristics of changing composition with time. When we say "deal with," we are referring to a systematic way of deducing the composition of the complex chemical mixture as well as the effect(s) from exposure to such a complex chemical mixture. In other words, we are implying some sort of predictive capability. For many years even the eternal optimists have not been able to provide a reasonable answer to this question, though we believe intuitively that there might be a solution to this complex problem.

Soon after the last mixture conference at Colorado State University (1), we felt, for the first time since our involvement with chemical mixture research, that there is hope in dealing with complex chemical mixtures. The approach we are advocating is integrated computer modeling of complex biologic processes. In this article we begin the discussion with some background information, present some recent successes in computer modeling of relatively simple chemical mixtures (i.e., fewer than five chemicals), and then conclude the discussion by introducing reaction network modeling and its integration into physiologically based pharmacokinetic/pharmacodynamic (PBPK/PD) modeling.

Experimental and Computational Approaches

Chemical Mixtures and Multiple Stressors

Just as we cannot ignore the scientific issues on chemical mixtures because they are complex, we should not be looking at chemicals alone when we are interested in the global issue of public health. The Gulf War syndrome taught us to look beyond the chemicals into the area of multiple stressors (2,3). Thus, the smoke in the burning oil field is but one piece of the puzzle in the overall assessment of Gulf War syndrome (2,3).

A more detailed discussion on multiple stressors was given elsewhere (3). Briefly, any chemical, physical, or biologic insult on the body is a form of stress; therefore, multiple stressors include chemicals, drugs, and physical and biologic agents. However, in the context of the Gulf War syndrome, multiple stressors may also include environmental hardship (e.g., extreme heat, poor resting conditions, poor food or water intake, heavy and nonbreathable equipment and clothing, insect or other pests), occupational hazard (e.g., dangerous tasks; injuries from work; exposure to fuels, burning oil fields, possible nerve gases, radioactive residues), and psychologic stress (e.g., threat of death and injuries, fear of exposure to chemical and biologic warfare agents, being away from home, poor living conditions).

If chemical mixtures are already sufficiently complicated, would not the addition of multiple stressors render the situation impossible? Indeed it might. However, our thinking is that although the number of combinations of chemicals or stressors is infinite, the number of biologic processes is finite. Therefore, in considering an integrated computer modeling approach, we must work on the finite biologic processes rather than the infinite combinations of chemicals and stressors.

Simple Chemical Mixtures: Interaction Thresholds

We first demonstrate one utility of PBPK modeling by estimating the threshold point for toxicologic interactions in the low occupational exposure region. Co-exposure to multiple chemicals may significantly affect the pharmacokinetics of one or more mixture components and alter the target tissue dose of the toxic moiety. In 1996 we introduced the idea of interaction thresholds as the minimal level of change in tissue dosimetry associated with a significant health effect (4). When two or more interactive chemicals are studied together, theoretically there could be infinite interaction thresholds. This is because, in the case of a binary or higher-order chemical mixtures, if we vary the concentrations of two or more chemicals, we would get, theoretically, an infinite number of interaction thresholds. However, if we specify certain occupational or environmental exposure concentrations for all the other chemicals in the mixture except one, we may obtain an interaction threshold for that set of exposure conditions. This is important because human risk from exposure to multiple chemicals may not always obey the rule of additivity. In a 2001 publication from our laboratory (5), Dobrev et al. estimated the interaction thresholds of three common volatile organic solvents, trichloroethylene (TCE), tetrachloroethylene (perchloroethylene, PERC), and 1,1,1-trichloroethane (methyl chloroform, MC), under different dosing conditions. Briefly, an interactive PBPK model was built where PERC and MC are competitive inhibitors for TCE. The model was developed and validated by gas uptake pharmacokinetic studies in Fischer 344 (F344) (Harlan Sprague Dawley, Indianapolis, IN, USA) rats at relatively high doses of single chemicals, binary mixtures, and the ternary mixture. Using computer simulation to extrapolate from high to low concentrations, we investigated the toxicologic interactions at occupational exposure levels, specifically at around threshold limit value/time-weighted average (TLV/TWA). Because long-term toxicity and carcinogenicity of these three solvents are clearly associated with their metabolism, and TCE is the most extensively metabolized among them, we focused our study on changes in internal TCE dose measures related to the mixture co-exposure. Using a 10% elevation in parent compound blood level as a criterion for significant interaction, we estimated interaction thresholds with two of the three chemicals remaining at constant concentrations. Thus, we fixed the TCE and PERC exposure concentrations in the gas uptake pharmacokinetic studies to 50 and 25 ppm, respectively, their TLVs/TWAs, and estimated the interaction threshold by varying the exposure concentration of MC to 0, 100, 130, 175, 250, or 350 ppm (350 ppm is the TLV for MC). This latter work is based on computer simulations using the interactive PBPK model; thus, it is experimentation in silico. Dobrev et al. (5) reported that under the above exposure conditions (i.e., TCE and PERC at their TLVs), the interaction threshold for the ternary mixture is 50, 130, and 25 ppm for TCE, MC, and PERC, respectively. If one wishes to use a higher criterion than 10% elevation in blood level for interaction threshold, Dobrev et al. (5) provide possible interaction thresholds for 17% (50, 250, 25 ppm for TCE, MC, PERC) and 22% (50, 350, 25 ppm for TCE, MC, PERC) elevations in blood level of TCE. This work has recently been extended, in silico, to human exposure to this three-chemical mixture and the estimation of interaction thresholds for humans (6).

Simple Chemical Mixtures: Mixture Formula

In the derivation of an occupational exposure limit (OEL) of chemical mixture exposure, the general approach is to first determine if those chemicals in the mixture cause similar toxic responses and then to implement the "mixture formula" (Equation 1) to assess if the exposure might be problematic.

[1]

The mixture formula, also referred to as the "unity calculation" (denoted E_m), calculates the ratio of worker exposure to OEL for each chemical in the mixture. If the sum of these ratios exceeds unity (1.0), an overexposure is suggested. The potential problems with this approach are that it assumes additivity, and pharmacokinetics and pharmacodynamics are not taken into consideration. To illustrate these points, Dennison et al. (7,8) first modified the E_m by incorporating PBPK modeling and came up with a new pharmacokinetically based E_m, the E_mPK, shown in Equation 2:

[2]

where C_i is the concentration of the chemical in the target tissue (obtained through PBPK modeling) either during an exposure to mixtures or to the OEL for the single chemical. As with the E_m formula, the E_mPK is a summation of ratios of the doses for the chemical in a mixture to those of the single chemicals. Using an established PBPK model for alkyl-substituted benzenes published by Tardiff et al. (9), Dennison et al. (7,8) took into consideration both pharmacokinetics and the interactive enzyme inhibition among the component chemicals in the mixture. Dennison et al. (7,8) then gave a few case studies, one of which is reproduced in Table 1, to illustrate the differences between the conventional E_m approach versus the E_mPK approach.

As shown in Table 1, Dennison et al. (7,8) provided five hypothetical chemical mixtures by varying the concentrations of the three component chemicals, toluene, ethylbenzene, and xylenes. The four columns of data on the left are for the E_m approach and the four columns of data on the right are derived from the E_mPK approach. Mixtures 2-5 are all at allowable exposure concentrations under the present E_m approach because the E_m values derived are all under unity. However, when we take pharmacokinetics (tissue dosimetry from PBPK modeling) and pharmacodynamics (critical effects) into consideration, we see an immediate problem, which is explained below.

Critical effects for toluene and ethylbenzene, which are used to set the OELs, include depression of the central nervous system. Because it is not clear if xylene also causes such a critical effect, under the E_m approach, xylene is not considered in the derivation of E_m. This results in the E_m values being less than unity in mixtures 2-5. However, as xylene interferes with the metabolism of toluene and benzene thereby increasing their tissue dosimetry, the E_mPK derived in mixtures 2-5 went over unity by about 1.2- to 3-fold (Table 1). Thus, under the E_mPK approach, these mixtures would have been unallowable. This exercise demonstrates some inadequacies in current risk assessment methodology for chemical mixtures.

A Ray of Hope for Complex Chemical Mixtures: Reaction Network Modeling

Reaction network modeling and its application to petroleum engineering. Reaction network modeling has been used very successfully in the area of petroleum and chemical engineering for very complex chemical processes involving hundreds and even thousands of chemicals. It had never been applied to biology until the interdisciplinary collaborative effort between Colorado State University and Rutgers University. To appreciate the potential of this modeling approach, it is helpful to understand its historical development.

Reaction network modeling in the fields of chemical and petroleum engineering has progressed greatly over the past 25 years, including developments in the areas of group contribution methods (10), graph theory (11,12), Monte Carlo techniques (13-15), and quantum chemistry (16,17).

Intense activity in molecular reaction engineering has generated several approaches to the creation of molecular reaction models by computer (18-23). The essential idea is to input some representation of reactant structure and chemical reaction. Algorithms and grammar for representing and determining species connectivity, uniqueness, and relationships (the reaction network) create an output that is, in some form, the controlling kinetic equation (governing ordinary differential equations) for the reaction model.

Developing a fundamental kinetic scheme requires modeling of the chemistry at the mechanistic level. This, in turn, leads to a large number of species, reactions in the governing network, and associated rate constants. Therefore, although a mechanistic model incorporating detailed molecular representations and fundamental kinetic data is needed, the inherent complexity and size of such a model is a deterrent to its development. This motivates a simplification of the system through the organization of species into related families of compounds and the automation of not only the solution but also the construction of the model.

In an attempt to describe hydrocarbon mixture properties, early models relied on the techniques of lumping (24), where a relatively small number of lumps were used to describe the mixture. In these coarsely lumped kinetic models, thousands of individual constituents in a complex feedstock were grouped into broad but measurable categories of compound classes or boiling range, with simplified reaction networks between the lumps.

More recently, Quann and Jaffe (23) developed a more fine-grained lumping approach they named structure-oriented lumping (SOL). SOL was developed in response to the need for incorporating molecular detail in petroleum chemistry to predict product compositions and properties. It is a group contribution method describing the structure of molecules, facilitating both molecular property estimation and a description of process chemistry. It was also designed to be consistent with current limitations of analytic capability to determine molecular detail. The concept central to the SOL approach is that any molecule can be described and represented by a set of certain structural features or groups. The SOL method organizes this set as a vector, with the elements of the vector representing the number of specific structural features sufficient to construct any molecule. Different molecules with the same set of structural groups, i.e., certain isomers, are lumped and represented by the same vector. The structure vector provides a framework to enable rule-based generation of reaction networks and rate equations involving thousands of components and many thousands of reactions.

An even finer-grain methodology was developed by Broadbelt and co-workers (16,25,26), who used concepts of graph theory to represent species connectivity. They also made use of computational quantum chemistry (CQC) and linear free-energy relationships (LFERs) (27) to automate the process of determining reaction rate constants. CQC was applied to determine the optimal conformation and molecular properties, such as the electron affinity, electron density, bond order, and heat of formation, associated with each of reactant and product structures (17,28). LFERs were then used to give a correlation of the rate or equilibrium constants with a property of a molecule or intermediate for a family of reactions. Thus, the CQC calculations ultimately provided an estimate of rate or equilibrium coefficients for reactants that had not been studied experimentally.

This general framework, which Klein (29) refers to as the kinetic modeler's toolbox (KMT), allows for the convenient construction and solution of even highly complex chemical reaction networks. Using this approach, various investigators have been able to encode complex hydrocarbon mixtures and create rule sets for a wide variety of reactions within the mixture. Results from model simulations (2,30) have shown good agreement with experimental observations in tracking the evolution of thousands of molecular components, then predicting mixture properties such as normal boiling points, specific gravity, and narrow boiling-cut yields. Joshi and co-workers (22) made use of the KMT for the analysis of gas oil catalytic cracking and were able to derive optimized parameter values (activation energies and frequency factors) and lumped fractions that were in good agreement with experimental results reported in the literature. Along similar lines, Mizan and Klein (31) found good agreement in terms of product yields and yield profiles between simulation results and experimental data in the reaction network modeling of
n-hexadecane hydroisomerization.

Reaction network modeling is a tool for predicting the amounts of reactants, intermediates, and products as a function of time for a series of coupled chemical reactions (potentially numbering in the tens of thousands of reactions for some systems). It is usually a mathematic and symbolic formulation suitable for solution on the computer. A reaction network model builder is a tool for the computer generation of a reaction network model. The model builder thus can be used not only to solve the kinetic equations of interest but also to generate the reaction mechanisms, rate constants, and reaction equations themselves.

Essentially, the model builder works as follows:

• The concentrations of the species to be reacted or metabolized are input to the model builder.
  
  
• For each species in turn, the model builder performs a test against each of a set of reaction rules to determine whether the species is susceptible to a particular chemical reaction.
  
  
• If it is not susceptible to any reactions, no further action is taken on this species.
  
  
• If it is susceptible, a transformation of the species into one or more product species is performed based on the particular chemical reaction.
  
  
• Each of these product species then undergoes the same susceptibility tests and a similar transformation sequence. This leads to a linking of all reactants with intermediate and ultimately with final products. This linking forms the structure of the chemical reaction network.
  
  
• After the reaction network is established, the rate constants for the reactions are retrieved or computed.
  
  
• The coupled differential equations governing the reaction kinetics for the network are then formulated by the model builder.

Finally, the kinetic equations, i.e., the model equations, are solved numerically, leading to the concentrations of all species as a function of time. For those interested in a more specific, detailed description of reaction network modeling, the article by Klein et al. in this monograph (32) should be consulted.

Reaction network modeling and its application to biomedical research. It is important to note that metabolism and toxic mechanisms of chemicals and chemical mixtures often involve complex reactions in which the species associated with one reaction are constituents of many other reactions. Interactions among chemicals are common, and this interplay among reaction pathways is the primary reason for toxicologic interaction. This interdependent, coupled set of biochemical reactions can be regarded as a reaction network and can occur for even relatively simple systems.

In the collaborative research effort between our laboratory (the Quantitative and Computational Toxicology Group) and Klein's group at Rutgers University, we intend to apply reaction network modeling to the metabolic pathways of complex mixtures in biologic systems. The approach and the modeling software being developed is "BioMOL," where "bio" represents biologic and "MOL" is the acronym for molecule-oriented lumping. The framework of KMT has served as the starting point to build BioMOL.

Figure 1. Major metabolic activation pathways for B[a]P in biologic systems. BPDE, 10-N+2-dG, 10-(deoxyguanosin-N+2-yl)-7,8,9- trihydroxy-7,8,9,10-tetrahydrobenzo[a]pyrene.

Our first application is reaction network modeling of benzo[a]pyrene (B[a]P). The metabolic pathway of B[a]P in a biologic system is used to build the first model by BioMOL. We use B[a]P because a) B[a]P is a human carcinogen; b) its metabolic pathway is extremely complex but extensively studied; and c) the reactions B[a]P undergoes, e.g., epoxidation and hydrolysis, are also common for other xenobiotics, especially polycyclic aromatic hydrocarbons. B[a]P is metabolically activated to its ultimate carcinogen, B[a]P-7,8-dihydrodiol-9,10-epoxide, via a series of reactions catalyzed by cytochrome P450 and epoxide hydrolase (Figure 1) (33-35). In addition to the formation of B[a]P-7,8-dihydrodiol-9,10-epoxide, B[a]P metabolism constitutes a complex reaction network in biologic systems. Briefly, B[a]P is first oxidized at several aromatic bonds to form arene oxides via reactions catalyzed by cytochrome P450. Once arene oxides are formed, they can undergo the following three reactions: a) rearrangement to phenols spontaneously through an NIH shift (36); b) hydrolysis catalyzed by epoxide hydrolase to form trans-dihydrodiols; and c) conjugation with glutathione, followed by excretion. Dihydrodiols might again form the corresponding epoxides in a reaction catalyzed by cytochrome P450 or undergo sulfation and glucuronidation. Bay-region dihydrodiol epoxides, e.g., B[a]P-7,8-dihydrodiol-9,10-epoxide, are unusually reactive to biologic macromolecules. The reactivity of these compounds is related to their resistance to enzymatic detoxification (37). The non-bay-region dihydrodiol epoxides may hydrolyze spontaneously to tetraols or conjugate with glutathione (33). The phenols may, in turn, be oxidized to quinones or undergo conjugation (34). On the other hand, a one-electron oxidation pathway may be responsible for the formation of 6-hydroxy-B[a]P and subsequent metabolites 1,6-, 3,6-, and 6,12-quinones (38).

Generation of B[a]P reaction networks. To generate the reaction network at molecule level, we use graph theory to convert chemical structures and reaction rules into computer code. The atomic connectivity of the molecule is represented by the bond-electron matrix. In bond-electron matrices, the off-diagonal elements denote the bond order between two atoms, and the diagonal elements represent the unpaired electron (e.g., free radicals). The structure changes of molecules caused by chemical reactions can also be described by graph theory, i.e., the reaction matrix. Most reactions involve the change of connectivity between only few atoms. Therefore, reaction matrices can be used as a compact representation of bond formation and breakage caused by chemical reactions. The addition of the reaction matrix to the reduced reactant bond-electron matrix, which consists of only the atoms involved in the reaction, forms a new matrix representing the product of the reaction (i.e., product matrix). The epoxidation of an aromatic bond is used as an example to describe the matrix operation of a chemical reaction (Figure 2). The product is then checked for its uniqueness to ensure the molecule was not previously created by other reactions. If the product is unique, it is added to the unreacted species list and could become the candidate reactant for other reactions.

Figure 2. Matrix representation of epoxidation reaction in aromatic bond.

Table 2

The reaction rule, determined by modeler's understanding the fundamental chemistry and biochemistry of the reaction, plays a central role in the model. Following the reaction rules, BioMOL can search the eligible site of reaction throughout the structure of possible reactants and create the corresponding product by means of the matrix operator. The reaction sites and matrices for major metabolic pathways of B[a]P are summarized in Table 2. Any double bond in B[a]P is eligible for the epoxidation reaction. However, certain restrictions should be applied in the reaction rule to eliminate unrealistic products. For instance, arene oxides are unstable and not likely to stay long enough to undergo the second consecutive epoxidation reaction. Instead, arene oxides are more likely to be a substrate for hydrolysis, NIH shift, and glutathione conjugation (Table 2). To form the ultimate carcinogenic metabolite B[a]P-7,8-dihydrodiol-9,10-epoxide, B[a]P should encounter two distinct epoxidation reactions. In the BioMOL model, these two steps use the same reaction rule and reaction matrix operator. The only difference is that the rule and matrix are applied to distinct reactants, namely, B[a]P and B[a]P-dihydrodiol. Therefore, a complex reaction network may start from few reactions with few reaction rules.

Estimation of B[a]P metabolic reaction rate constants. By applying the given reaction rule, the BioMOL model generates all the possible reactions and corresponding products. Some of the products may never be observed because the reaction rates are either too low for their formation or too fast for their subsequent metabolism. BioMOL will use quantitative structure/reactivity correlations (QSRCs) to estimate the reaction rate constants (k_i). QSRCs are semiempiric methods in which sterically similar reactions are lumped into reaction families. The idea of this correlation is expressed by the following equation:

[3]

where i is a component in the reaction family j. RI represents the reactivity index. The parameters a and b are calibrated by experimental data (39). The candidates for the reactivity index of B[a]P biotransformation are heat of reaction, 1⁄4-electron density, and bond order. Semiempiric quantum chemistry software like MOPAC 2000 (Schrodinger, Portland, OR, USA), which is integrated into the BioMOL model, can calculate these reactivity indices.

Figure 3. Classification of rate constants for B[a]P metabolic activation pathways. EH, epoxide hydrolase.

To calculate the reaction rate based on QSRCs, the rate constants for B[a]P metabolic activation have been divided into four groups, enzymatic, nonenzymatic, association, and dissociation (Figure 3). The QSRCs for enzymatic rate constants rely on the understanding of the chemical mechanism. For instance, the enzymatic epoxidation of double bonds catalyzed by cytochrome P450 is likely to start with the formation of a charge-transfer complex between the ferryl oxygen on P450 and 1⁄4 bond on the substrate (40). Therefore, the rate constants of the enzymatic epoxidation at distinct positions of BaP are expected to correlate with bond order of the reaction site. Loew et al. (41) have shown a qualitative correlation between the observed B[a]P metabolites and 1⁄4-bond reactivity calculated by quantum chemistry. The relative abundance of the two possible phenol products resulting from an NIH shift also correlated with the electron densities of the carbon atoms attached to oxygen (41). On the other hand, the nonenzymatic reaction rate constant may correlate with the heat of reaction. The real challenge for QSRCs is the search of reactivity indices for association and dissociation constants. The binding of substrates to enzymes is related to the three-dimensional structure of both the substrate and the enzyme, which is not fully understood. Our laboratory is currently developing the QSRCs from both analytic experiments and literature data. After reaction families are defined, the kinetics of the reaction network become accessible by using QSRC-derived parameters and reactivity indices to represent thousands of rate constants.

In summary, BioMOL is a computer-assisted modeling tool and a promising approach to handling extremely complex biologic systems and their related biochemical reactions and reaction networks. The BioMOL model has the potential for predicting a variety of reaction networks in biologic systems, as it is molecule based. The determination of reaction rules and reaction families should be supported by a solid understanding of reaction chemistry. Fundamentals of enzymatic kinetics will help QSRC estimation of reaction rate constants.

Future Perspective Second-Generation PBPK/PD Models

One of the areas of interest in our laboratory is combination chemotherapies. We are interested in any type of combination chemotherapy including cancer-, AIDS-, antibiotic-combination chemotherapies, and others. Thus, we will use cancer-combination chemotherapy as an example to illustrate our thinking in future directions. As reported in the literature (42), the potential target sites for cancer-combination chemotherapies are in basic biologic processes including purine/pyrimidine metabolic pathways leading to DNA synthesis, RNA production, and protein/enzyme and microtubule production (42). These biochemical pathways are already an immensely complex system without the addition of chemotherapeutic agents. Experimentally, it is almost impossible to study such a complex system simultaneously. Therefore, it is essential to build a simulation platform that can be used to globally examine these biologic pathways in normal and malignant cells with or without chemotherapeutic agents. A variety of terms such as "virtual cells," "virtual laboratories," and "in silico experimentation" have appeared in the recent literature. Our intent, as an interdisciplinary team of scientists and engineers, is to conduct experiments on computers, once a validated PBPK/PD model is available for these biologic pathways. With advances in computational technology, the complexity of these biologic pathways and interactions will not be a limiting factor in our understanding of the biologic foundation of combination chemotherapies.

From a modeling perspective, we believe cancer cells represent parameter changes in certain specific biologic processes inherent in normal cells. Similarly, the introduction of chemotherapeutic drugs into the cells also represents perturbations in biologic processes in normal cells. Thus, theoretically, we are proposing one PBPK/PD simulation model for the basal biologic pathways in normal cells; those same pathways in cancer cells or under combination chemotherapies are merely quantitative variations (i.e., parameter changes of certain processes) of the same model.

Our thinking goes beyond the traditional PBPK modeling. We plan to incorporate the concept of reaction network modeling by incorporating the essential biologic pathways involved. For instance, if we are studying breast cancer-combination chemotherapy, we will be adding scores of biochemical reactions into the breast tissue compartment. We will also identify specific genomic and proteomic changes in relation to the biochemical reactions to elucidate mechanistic bases for combination chemotherapies. Thus, pharmacodynamic interaction(s) will be incorporated into the PBPK modeling to transform the model to a PBPK/PD model. In that sense, we are interested in developing second-generation PBPK/PD models.

Figure 4. Enzymatic kinetics of the inhibition of FdUMP on dTMP synthesis by 5-FU. Abbreviations: DHF, dihydrofolate; dTMP, deoxythymidine monophosphate; dUMP, deoxyuridine monophosphate; FdUMP, 5-fluorodeoxyuridylate; MTHF, N5,N10-methylenetetrahydrofolate; PPi, pyrophosphate; PRPP, 5-phosphoribosyl-a-pyrophosphate. Equations were derived based on the mechanisms proposed by Hardy et al. (43) and Voet and Voet (44).

Figure 4 provides an example of such second-generation PBPK/PD modeling. The enzymatic pathway in Figure 4 is merely a small portion of the reaction network of purine/pyrimidine metabolic pathways involved in the cancer-combination chemotherapies. Here, we summarize the mechanism of action of 5-fluorouracil (5-FU) on thymidylate synthase (TS), an important enzyme in DNA synthetic pathway, and provide the enzyme kinetics involved in TS inhibition. The kinetic equations will be incorporated into the breast tissue compartment in the PBPK model for 5-FU. The reaction rate constants for the enzymatic processes may be obtained in three different ways: a) by mining the literature for quantitative data (i.e., from step 1); b) by calculation via semiempiric quantum chemical methods based on known enzyme-substrate molecular interactions or from QSRCs (32); and c) by acquisition through laboratory experimentation with relevant systems. Through such integration of individual reactions into the PBPK/PD modeling process, we will be going through a "de-lumping" process analogous to that which occurred in the field of chemical engineering in the last 30 years or so. In our case we may consider that the transformation of classic compartmental pharmacokinetics into PBPK modeling represents the first stage of de-lumping--that is, going from 2 or 3 compartments to 5-10 or 15 compartments. The second-generation PBPK/PD models will mean that, in the critical compartment, further de-lumping will carry PBPK/PD modeling down to the molecular mechanism level. If we fully utilize the power of computational technology, such second-generation PBPK/PD modeling will have the potential to handle very complex biologic systems, thereby handling exposure to complex chemical, biologic, and physical agents.