September 22, 2017, Friday
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at Boulder CU 
Search Links

MBW:PBPK Modeling of Metabolic Pathways of Bromochloromethane in Rats

From MathBio

Jump to: navigation, search


This is a summary and expansion of the paper titled "Physiologically Based Pharmacokinetic (PBPK) Modeling of Metabolic Pathways of Bromochloromethane in Rats" by Cuello et a.l [1] which can be found here.

Executive Summary

The authors set out to to study the toxicity of Bromochloromethane (BCM), a common by-product of water chlorination. Specifically, the goal of this paper was to determine how BCM gets metabolized in the body because the metabolic pathway that occurs will alter how toxic the chemical really is. In order to do this, the authors used a physiologically-based pharmocokinetic (PBPK) model, a common mathematical modeling technique used in conjunction with gas-uptake data to study the absorption, distribution, metabolism and excretion of various chemicals [2]. PBPK models are composed of a number of compartments that represent certain organs and parts of the body, and creates differential equations representing the change in the amount of the chemical in each compartment over time. Utilizing previous data from a gas-uptake experiment designed by the EPA, the authors developed a PBPK model and compared it with the true data in order to test different metabolic hypotheses. The authors believed that a new, two-binding site hypothesis in which metabolism of BCM occurs on one enzyme can similarly describe the data as the more-common two-pathway hypothesis in which metabolism responsibilities are shared by the CYP2E1 and glutathione transferase enzymes. Their results show that both hypotheses describe the experimental data in a similar manner and provide optimized metabolic parameters that match previously reported results. In order to expand on their work, I chose to improve the efficiency of the biological parameter storage and the accuracy of the metabolic parameter optimization.

History/Biological Context

The biological context of this paper is similar in nature to this article on determining metabolic constants from gas uptake data. The overall goal of most of these papers is to utilize the accuracy of PBPK simulations in order to provide some insight into the risk-assessment of a certain chemical. PBPK models are not a new contribution to the field of mathematical biology, but they also have not been around for decades [3]. These tools have not changed much structurally over the years, but instead have changed in more subtle, but very important, ways. More specifically, though the actual design of the experiment in which the gas uptake data is collected and the general form for the differential equations have not changed much, the metabolic rate equations and the selection of optimal metabolic constants has seen a huge overhaul. Instead of altering the metabolic constants by way of trial-and-error until the curves match the data, this paper focuses on mathematically finding the optimal values [4]. The inclusion of 3D error plots that I developed is believed to be one of the first attempts to graphically display the error over time, concentration, and parameter values. In addition, unlike most other PBPK models, this paper sought to incorporate the two-binding site hypothesis.

Mathematical Model

A. Compartments and Biological Parameters The compartments chosen for this paper (chamber, lungs, rapidly perfused tissue, adipose tissue, slowly perfused tissue, liver, kidney, arterial blood) were based off their relevance to the metabolism of the chemical. The specific values for the flow rates, volumes and partition coefficients were all calculated and can be seen below in Table 1:


B. Differential Equations The model is composed of six differential equations developed using the law of conservation of mass: BCM should neither be created nor destroyed. Moreover, the mass balance equations were determined by incorporating tissue volumes, blood flows, and partition coefficients that correspond to the specific male Fischer-344 rats used. However, the compartments of most interest were the "inhalation" or chamber compartment as well as the kidney and liver compartments, as they are where metabolism takes place:

{\frac  {dA_{{inh}}}{dt}}=N*QV*(C_{{exh}}-C_{{inh}})-loss*A_{{inh}}

{\frac  {dA_{{kid}}}{dt}}=q_{4}*(C_{{art}}-C_{{kid}})-Met_{k}

{\frac  {dA_{{liv}}}{dt}}=q_{5}*(C_{{art}}-C_{{liv}})-Met_{l}


N is the number of rats,

QV is the ventilation rate,

q_{i} is the flow rate for compartment i

A_{{comparment}} is the amount of BCM in compartment ,

C_{{compartment}} is the concentration given by C_{{compartment}}={\frac  {A_{{compartment}}}{V_{{compartment}}*P_{{compartment}}}} and P_{{compartment}} represents the partition coefficient for the compartment,

Met_{k} is the kidney metabolism,

Met_{l} is the liver metabolism.

Originally, there were two additional differential equations: one representing the amount of BCM in the venous compartment, and the other representing the amount of BCM in the lungs compartment. Matlab simulations demonstrated that both of these equations approached equilibrium very quickly, thereby allowing the authors to set each of these equations to zero. This led to algebraic expressions for the concentration of BCM in venous compartment and arterial compartments respectively.

C. Metabolic Hypotheses

This paper tested two different metabolic hypotheses using the basic PBPK model structure outlined above. Previous studies [5] demonstrated that metabolism by a single CYP2E1 binding site using typical Michaelis-Menten enzyme kinetics was not sufficient to describe the closed chamber data. Instead, it has been shown that a linear term was necessary to describe the closed chamber data at high concentrations. For the two-binding site hypothesis, a modified Michaelis-Menten equation that incorporates a second site for the enzyme is used to describe metabolism. Thus, The only difference between the two hypotheses' models lies in the metabolism equations outlined bel_low:


Met_{l}=.948\left({\frac  {V_{{max}}C_{{liv}}}{K_{m}+C_{{liv}}}}+k_{{gst}}C_{{liv}}V_{{liv}}\right)

Met_{k}=.052\left({\frac  {V_{{max}}C_{{kid}}}{K_{m}+C_{{kid}}}}+k_{{gst}}C_{{kid}}V_{{kid}}\right)

Two-binding site:

Met_{l}=.948\left({\frac  {V_{{max1}}C_{{liv}}+CL_{2}C_{{liv}}^{2}}{K_{{m1}}+C_{{liv}}}}\right)

Met_{k}=.052\left({\frac  {V_{{max1}}C_{{kid}}+CL_{2}C_{{kid}}^{2}}{K_{{m1}}+C_{{kid}}}}\right)


V_{{max}} is the maximum metabolic rate,

K_{m} is the Michaelis constant,

k_{{gst}} is the linear term for glutathione transferase,

V_{{max}} is the maximum metabolic rate for the first binding site,

CL_{{2}} is the clearance constant consisting of the ratio of the maximum metabolic rate and the affinity constant, {\frac  {V_{{max2}}}{K_{{m2}}}} for second binding site.


The results of metabolic parameter optimization are shown below in Table 2:


As hypothesized, the two-pathway and two-binding site models produced very similar graphs and fit the data equally well. In addition, the V_{{max}} and K_{m} values were very similar for either model. Plots of these models versus the experimental data are also seen below.



Thus, after comparing the graphs of the two models and their metabolic parameter values, the authors concluded that the two-binding site hypothesis is a plausible explanation for the metabolism of BCM.

Additional Analysis

This section covers that I did to contribute to the previously published paper. As previously noted, the authors were successful in showing that the metabolic hypotheses produce similar results and lead to almost identical metabolic constants. However, the most major downfall of their work, which can be shown via reproduction of their model, is that the process of producing the metabolic constants via optimization that line up with previous studies and literature is very sensitive. The authors use of Matlab's fminsearch command was successful in the sense that it produced the ideal values, but only after giving the optimizing program initial guesses that are fairly close to the ideal values. This is obviously not ideal.

Thus, I took on the task of finding a different way to select which parameters are most optimal. Methods such as establishing confidence region by adding noise, Matlab's SimBiology toolbox, and the common Eadie-Hofstee diagram were all considered, but in the end the process I chose was to create a 3D error surface that gives a visual demonstration of the parameter values that correspond to the lowest error. After all, the parameters associated with the lowest error values should be the ones we deem optimal. This was done by adapting the code to create and return the normalized error matrix that corresponds to the error for a given parameter value (holding the others constant) and time step. From there, all that was left was to organize this matrix by concentration and plot the the meshgrid. This was only done for the two-pathway model because it is the most biologically-accepted hypothesis. Examples of these 3D plots are seen below for K_{m}, <V_{max}</math> and k_{{gsh}} respectively:




Further inspection of these plots led to the following conclusions about the optimal values: K_{m} should be around .2-.4, V_{{max}} should be around 3.5, and K_{{gsh}} should be around 4.45. All of these values match up with what the literature says we should get, as well as the values from the paper. This confirms the authors' conclusions that the most optimal values are the ones presented, but does so more convincingly than their use of fminsearch does.

Moreover, since K_{m} and V_{{max}} are mathematically related, I decided to examine what happens when you produce error plots in which both parameters are being varied, not just one at a time. This resulted in further expansion of the authors code in which a series of 3D matrices were calculated, captured using Matlab's getFrame command, and glued together via Matlab's movie2avi command. This led to the same expected: the error value was least when both K_{m} and V_{{max}} were at their individually optimal values. Sample error movies are available upon request.


In order to assess the toxicity and risk a certain chemical poses to humans or other animals, PBPK models have been meshed together with gas inhalation data to study the total rate of metabolism as a function of the decrease in the concentration of the chemical in the chamber compartment. Because different metabolic pathways and enzymes produce different metabolites, PBPK modeling is one way to computationally examine risk assessments by extrapolating from rodent experiments to human populations [6]. As far as the authors were aware, this is the first study that made use of the two-binding site hypothesis in which metabolism of BCM involves the opening of a second binding site on a single enzyme after a certain concentration. This lies in stark contrast to the two-pathway model in which a second enzyme is responsible for higher concentration levels. The authors successfully showed that the two-binding site model produces values for the chamber concentration that are just as accurate as the scientifically-proven two-pathway model. Since these two pathways result in metabolites that correspond to different levels of toxicity, it's important to note this new hypothesis as a possible explanation for metabolism. Moreover, since the optimized metabolic parameter values match previous studies, the authors argue that their model is mathematically and biologically sound. However, when it came to finding the optimal metabolic constants, I found possible problems. Since the optimization routine was so dependent on the initial guesses for the constants, I produced 3D and 4D error plots that visually demonstrate how the error changes over time, concentration, and parameter values. Fortunately this work did reinforce the authors findings, thereby providing an additional reason why this model is to be trusted.


Due to the size and comlexity of the Matlab codes involved in this project, I have chosen not to provide them here. These codes are available upon request, along with additional 3D error plots and a sample error movie. For any additional information or clarifications, I suggest thoroughly reading the original publication


  1. W. S. Cuello, T. A. T. Janes, J. M. Jessee, et al., “Physiologically Based Pharmacokinetic (PBPK) Modeling of Metabolic Pathways of Bromochloromethane in Rats,” Journal of Toxicology, vol. 2012, Article ID 629781, 14 pages, 2012. doi:10.1155/2012/629781
  4. Andersen ME, Gargas ML, Jones RA, Jenkins LJ Jr. The use of inhalation techniques to assess the kinetic constants of 1,1-dichloroethylene metabolism. Toxicol Appl Pharmacol. 1979 Feb;47(2):395-409. PubMed PMID: 452031.
  5. M. L. Gargas, M. E. Andersen, and H. J. Clewell, “A physiologically based simulation approach for determining metabolic constants from gas uptake data,” Toxicology and Applied Pharmacology, vol. 86, no. 3, pp. 341–352, 1986.
  6. W. A. Chiu, M. S. Okino, and M. V. Evans, “Characterizing uncertainty and population variability in the toxicokinetics of trichloroethylene and metabolites in mice, rats, and humans using an updated database, physiologically based pharmacokinetic (PBPK) model, and Bayesian approach,” Toxicology and Applied Pharmacology, vol. 241, no. 1, pp. 36–60, 2009.