Simulation of the crosstalk between glucose and acetaminophen metabolism in a liver zonation model

The liver metabolizes a variety of substances that sometimes interact and regulate each other. The modeling of a single cell or a single metabolic pathway does not represent the complexity of the organ, including metabolic zonation (heterogeneity of functions) along with liver sinusoids. Here, we integrated multiple metabolic pathways into a single numerical liver zonation model, including drug and glucose metabolism. The model simulated the time-course of metabolite concentrations by the combination of dynamic simulation and metabolic flux analysis and successfully reproduced metabolic zonation and localized hepatotoxicity induced by acetaminophen (APAP). Drug metabolism was affected by nutritional status as the glucuronidation reaction rate changed. Moreover, sensitivity analysis suggested that the reported metabolic characteristics of obese adults and healthy infants in glucose metabolism could be associated with the metabolic features of those in drug metabolism. High activities of phosphoenolpyruvate carboxykinase (PEPCK) and glucose-6-phosphate phosphatase in obese adults led to increased APAP oxidation by cytochrome P450 2E1. In contrast, the high activity of glycogen synthase and low activities of PEPCK and glycogen phosphorylase in healthy infants led to low glucuronidation and high sulfation rates of APAP. In summary, this model showed the effects of glucose metabolism on drug metabolism by integrating multiple pathways into a single liver metabolic zonation model.

Introduction

The liver performs hundreds of vital functions such as metabolism of exogenous and endogenous compounds, bile production, and maintenance of blood levels of many substances, including albumin, glucose, amino acids, and vitamins (Trefts et al., 2017). Anatomically, the human liver contains approximately a million hepatic lobules with hexagonal cross-sections, and each lobule consists of a million hepatocytes and a thousand capillaries called sinusoids. The hepatic artery and the portal vein from the intestine merge at one end of the sinusoid in the periportal region of the lobule and flow into the central vein at the other end of the sinusoid in the pericentral region of the lobule (Abdel-Misih and Bloomston, 2010). Metabolic enzymes are expressed at different levels in hepatocytes along with the sinusoidal axis, which compartmentalizes hepatic functions (Halpern et al., 2017; Kietzmann, 2017; Halpern et al., 2018; Ben-Moshe et al., 2019). This spatial patterning of hepatic functions (hepatic zonation) presumably contributes to the efficient functional management of the whole organ since certain pathways antagonize each other.

Acetaminophen (APAP) is one of the most commonly used over-the-counter medications; however, its overdoses can lead to hepatotoxicity and acute liver failure (Sacks et al., 2018). APAP is intensively used as a model compound to study liver drug metabolism, and it is well known that hepatocytes in the pericentral region are initially damaged due to metabolic zonation (Anundi et al., 1993; Cunningham and Porat-Shliom, 2021). Multiple metabolic pathways are involved in APAP metabolism, and their mutual interaction and regulation add another layer of complexity to understanding the metabolic fate of APAP in the liver (Mazaleuskaya et al., 2015). In hepatocytes, APAP is eliminated through glucuronidation or sulfation or converted to the toxic intermediate metabolite N-acetyl-p-benzoquinone imine (NAPQI) by cytochrome P450 2E1 (CYP2E1), which is detoxified and eliminated through glutathione (GSH) conjugation. In adults, glucuronidation and sulfation play a major role in APAP metabolism, and approximately 5% of the substrate is converted into NAPQI (Court et al., 2001; Mutlib et al., 2006; Riches et al., 2009). Sulfation activity and GSH concentration are high in the periportal region, while glucuronidation, CYP2E1, and GSH conjugation are high in the pericentral region (Kietzmann, 2017). GSH is depleted in the pericentral region leading to site-specific accumulation of NAPQI, which covalently binds to cysteine (Cys) residues in proteins to form protein adducts, which cause hepatotoxicity (Yoon et al., 2016).

Glucuronidation accounts for the majority of APAP metabolic pathways in healthy adults (Jiang et al., 2013). Interestingly, uridine diphosphate (UDP)-glucuronic acid (UDP-GA), a substrate of glucuronide conjugation, is enzymatically produced from a glycogen metabolism metabolite constituting the glucose metabolism cascade (Adeva-Andany et al., 2016). This indicates that drug metabolism and glucose metabolism interact (Ghafoory et al., 2013; Dargue et al., 2020; Cunningham and Porat-Shliom, 2021). Glucose metabolism consisting of glycolysis, gluconeogenesis, and glycogen metabolism shows metabolic zonation in the liver (Kietzmann, 2017). Glucose uptake and glycolysis are promoted in the pericentral region, while gluconeogenesis and glucose delivery dominate in the periportal region.

The precise evaluation of complex human liver metabolism and toxicity may include in silico approaches as complementary or promising alternatives to existing in vivo and in vitro approaches. The first mathematical model for APAP metabolism (Reith et al., 2009) was followed by a multi-compartmental model (Ben-Shachar et al., 2012), a physiologically based pharmacokinetic model (Ben-Shachar et al., 2012; Jiang et al., 2013), and simplified models dealing with only major pathways (Remien et al., 2012; Reddyhoff et al., 2015). Recently, the zonation effects on APAP metabolism and the zone-specific hepatotoxicity were simulated in mathematical models (Smith et al., 2016; Franiatte et al., 2019; Kennedy et al., 2019; Means and Ho, 2019). Meanwhile, mathematical models incorporating zonal regulations of glycolysis, gluconeogenesis, and glycogen metabolisms were reported (Bulik et al., 2016; Berndt et al., 2018; Berndt and Holzhütter, 2018). However, most existing models considering metabolic zonation deal with drug metabolism or glucose metabolism individually. Thus, it is difficult to investigate crosstalk effects that impair the accuracy of evaluation under different clinical conditions.

In this study, we integrated metabolic zonation models of multiple pathways including drug and glucose metabolism into a single numerical human liver model. This model showed that changes in glucose metabolism enzyme activities affected APAP metabolism, and the crosstalk effects were dependent on liver zonation and nutritional status. These results showed the importance of incorporating multiple metabolic networks in a single zonation model.

Materials and methods

Model overview

We built a numerical liver zonation model by integrating glucose, APAP, and Cys metabolisms (Figure 1 and Supplementary Material for abbreviations). The glucose metabolism model consists of glycolysis, gluconeogenesis, and glycogen metabolism. Glucose metabolism rate equations were adopted and modified from previous studies (Bulik et al., 2016; Berndt et al., 2018; Berndt and Holzhütter, 2018). APAP metabolism utilized a previously described model (Means and Ho, 2019). Our model also contains Cys metabolism. We divided the liver into two zones for most cases and three zones for calculating glucose exchange rates and glycogen concentrations to compare with reported results (Berndt et al., 2018). The model consists of 34 variables (metabolites) and 47 reactions. We simulated metabolic zonation by changing enzyme activities depending on the zone that the cell belongs to: from the periportal-to the pericentral regions. For instance, the pericentral region has a greater expression of UTP-glucose-1-phosphate uridylyltransferase than the periportal region (Kietzmann, 2017). Thus, we used a higher Vmax value for UTP-glucose-1-phosphate uridylyltransferase in the pericentral cell. The model contains 13 zonation-dependent kinetic parameters (Supplementary Material). We assumed that the organism outside the liver was unchanged during simulations. The external glucose concentration and APAP administration rate were the boundary conditions, which we changed during simulations. We used the external glucose concentrations of 4 and 11 mM for the fasting and feeding states, respectively. We used the APAP administration rates of 0, 0.5, and 6 mM/h for no, moderate, and excessive administration, respectively. The liver zonation model was implemented in Python (ver.3.10, https://www.python.org) with NumPy library (ver. 1.21.0, https://numpy.org). The blood concentrations of glucose and APAP under different conditions, such as fasting, feeding, and moderate- and over-doses, were determined according to previous reports (Kennedy et al., 2019; Means and Ho, 2019; Lammers et al., 2020).

FIGURE 1

FIGURE 1. Network map of the constructed model. Schematic representation of the simulation model combining glucose (Glc), acetaminophen (APAP), and cysteine (Cys) metabolisms indicated by open blue, solid orange, and solid blue shapes, respectively. The definitions of the abbreviations are listed in the Supplementary Material. Glc, lactate, APAP, and Cys in the blood can be transported into the cytosol. The endoplasmic reticulum and mitochondria are shown on the background. The square boxes represent enzymes and transporters, whereas the circles represent substrates and metabolites.

Simulation method

Computer simulations of metabolic pathways often use dynamic simulation or metabolic flux analysis (MFA). The dynamic simulation reproduces the time-dependent transitions of metabolic concentrations but requires many parameters. For instance, the reversible Michaelis–Menten rate equation (commonly used for dynamic metabolic models) requires at least four parameters per enzyme reaction, leading to hundreds of parameters to construct a single metabolic model. Due to the limits of experimental data availability, optimizing many parameters is difficult. Meanwhile, MFA does not require any kinetic parameters. However, as it assumes that the entire system is at a steady state, a simple MFA cannot simulate any temporal changes in metabolite concentrations. Therefore, we used a hybrid approach in this study to compensate for the shortcomings of each method (Yugi et al., 2005). Briefly, dynamic simulation was applied to the rate-limiting reactions of the metabolic pathways, while MFA was utilized for the other reactions by assuming a pseudo-steady state at each simulation step.

In general, the dynamic behavior of a metabolic pathway is modeled by the following differential equation.

$\frac{\mathit{d}\mathit{x}}{\mathit{d}\mathit{t}}=\mathbf{S}\mathit{v}(1)$

where x is the variable vector representing metabolite concentrations (n × 1), S is the stoichiometry matrix (n × m), and v is the reaction rate vector (m × 1), which is a function of x. We divided the system into two modules: dynamic and static. The dynamic module contains key reactions that are bottlenecks of dynamic behaviors. The static module contains fast, non-bottleneck reactions. We divided v into reaction rates in the dynamic module (v_d) and the static module (v_s). We divided x into the variables only associated with dynamic module reactions (x_d) and those with one or more static module reactions (x_s). We divided S into three matrices: S_d,d, S_s,d, and S_s,s. S_d,d is the stoichiometry matrix for the dynamic module variables and dynamic module reactions, S_s,d is for the static module variables and dynamic module reactions, and S_s,s is for the static module variables and static module reactions. Our model has 13 variables and 27 reactions in the dynamic module and 21 variables and 20 reactions in the static module. The incorporation of the new symbols allows the previous equation to be written as:

$\{\begin{array}{c}\frac{\mathit{d}{\mathit{x}}_{\mathit{d}}}{\mathit{d}\mathit{t}}={\mathbf{S}}_{\mathit{d},\mathit{d}}{\mathit{v}}_{\mathit{d}}\\ \frac{\mathit{d}{\mathit{x}}_{\mathit{s}}}{\mathit{d}\mathit{t}}={\mathbf{S}}_{\mathit{s},\mathit{d}}{\mathit{v}}_{\mathit{d}}+{\mathbf{S}}_{\mathit{s},\mathit{s}}{\mathit{v}}_{\mathit{s}}\end{array}(2)$

S_d,d, S_s,d, and S_s,s values are given. v_d is a set of kinetic rate equations such as mass actions and Michaelis-Menten rate equations, i.e., v_d is a function of x_d and x_s. Thus, the right-hand side of the upper equation can be calculated.

To calculate the right-hand side of the following equation, we assume the static module is at a steady state at each time step in a simulation. Thus, the following equation becomes:

$\mathbf{0}={\mathbf{S}}_{\mathit{s},\mathit{d}}{\mathit{v}}_{\mathit{d}}+{\mathbf{S}}_{\mathit{s},\mathit{s}}{\mathit{v}}_{\mathit{s}}(3)$

In the over-determined case (including our model), there is no solution space for v_s to satisfy the previous equation. Instead, we can obtain ${\mathit{v}}_s^{pseudo}$, the most plausible alternative for v_s, by the following equation:

${\mathit{v}}_{\mathit{s}}^{\mathit{p}\mathit{s}\mathit{e}\mathit{u}\mathit{d}\mathit{o}}={{-\mathbf{S}}_{\mathit{s},\mathit{s}}^{\#}\mathbf{S}}_{\mathit{s},\mathit{d}}{\mathit{v}}_{\mathit{d}}(4)$

where ${\mathrm{S}}_{s,s}^{\#}$ is the Moore-Penrose pseudo-inverse of S_s,s. ${\mathit{v}}_s^{pseudo}$ provides the least-squares estimate of the reaction rate distribution, which minimizes ${\left|{\mathrm{S}}_{s,d}{\mathit{v}}_d+{\mathrm{S}}_{s,s}{\mathit{v}}_s\right|}^2$. This procedure equally distributes the error among the reaction rates of the static module (Yugi et al., 2005). Replacement of v_s in the second equation with ${\mathit{v}}_s^{pseudo}$ enables time evolution simulation of x_d and x_s.

Results

Comparison of the liver model with experimental data

The whole combined model was initially performed at different blood glucose levels without APAP administration. Metabolite concentrations and reaction rates are shown in Figures 2, 3, respectively. As expected, the substrate- and metabolites concentrations related to APAP metabolisms such as APAP, NAPQI, and protein adducts were zero. Furthermore, the reaction rates of APAP sulfation (v_Sulf), oxidation (v_CYP450), glucuronidation (v_UGT), and GSH conjugation (v_GSHT) were zero. In contrast, metabolite concentrations and reaction rates related to glucose metabolism changed in response to the blood glucose level and reproduced zonation patterns with the values consistent with the liver (Halpern et al., 2017, 2018; Ben-Moshe et al., 2019: (Halpern et al., 2018; Ben-Moshe et al., 2019). The pericentral region had higher reaction rates of glucokinase (v_GK), fructose-2,6-bisphosphatase (v_FBP2), aldolase (v_ALD), glyceraldehyde-3-phosphate dehydrogenase (v_GAPDH), and pyruvate kinase (v_PK) compared with the periportal region, and the positive lactate dehydrogenase reaction rate was indicative of glycolysis promotion and lactate (Lac) production (Miethke et al., 1985). The periportal region contained higher GSH concentration and greater reaction rates of fructose-1,6-bisphosphatase (v_FBP1) and phosphoenolpyruvate carboxykinase in the cytosol (PEPCK, v_PEPCK), and glucose-6-phosphate phosphatase in the endoplasmic reticulum (G6P, v_{G6P_ER}) compared with the pericentral region (Katz et al., 1977; Ma et al., 2020). Lactate dehydrogenase works in the reverse direction in the fasting state, which favors gluconeogenesis. In contrast to cytosolic PEPCK, mitochondrial PEPCK contribution to glucose metabolism is not fully understood. Nevertheless, our model predicted that the reaction rate of mitochondrial PEPCK (v_PEPCKmito) was higher in the pericentral region (Méndez-Lucas et al., 2013; Stark et al., 2014; Stark and Kibbey, 2014).

FIGURE 2

FIGURE 2. Metabolite concentrations without acetaminophen administration. Hepatocytes were in a feeding state from 0 to 12 h and in a fasting state for the rest of the simulation time. The external (blood) glucose levels were set to 4 and 11 mM for the fasting and feeding states, respectively. The solid red lines and blue dotted lines indicate the periportal and pericentral regions, respectively. The definitions of the abbreviations are listed in Supplementary Material.

FIGURE 3

FIGURE 3. Metabolite reaction rates without acetaminophen administration. Hepatocytes were in a feeding state from 0 to 12 h and in a fasting state for the rest of the simulation time. The external (blood) glucose levels were set to 4 and 11 mM for the fasting and feeding states, respectively. The solid red lines and blue dotted lines indicate the periportal and pericentral regions, respectively. The definitions of the abbreviations are listed in Supplementary Material.

The glucose exchange rate by glucose transporter 2 in response to the external (blood) glucose concentration and the change of glycogen concentration over time was investigated to further validate the glucose metabolism model (Figure 4). The glucose transporter 2reaction rate increased with increased external glucose; in other words, cellular glucose consumption increased (Figure 4A). Pericentral cells consumed glucose faster than periportal cells. Periportal cells exported 4–9 mM glucose from external glucose concentrations. The nutritional-status-dependent glycogen concentration showed that its levels were lowest in periportal cells in the fasting state (Figure 4B). However, glycogen accumulated faster in periportal cells in the feeding state (from 0 to 24 h) and reached higher concentrations than that of intermediate and pericentral cells. These glycogen concentration transitions were congruous with the reaction rates of glycogen phosphorylase and glycogen synthase (GS) shown in Figure 3. These glucose metabolism results were consistent with existing observations and experimental data (Berndt et al., 2018).

FIGURE 4

FIGURE 4. Validation of the constructed model. (A) Simulated glucose exchange rate by glucose transporter 2 in response to external (blood) glucose concentration. (B) Simulated glycogen concentrations over time. Hepatocytes were in a feeding state from 0 to 24 h and in a fasting state for the rest of the simulation time. The external (blood) glucose concentrations were set to 4 and 11 mM for the fasting and feeding states, respectively. The external lactose concentration was fixed at 1 mM in all simulations.

Cellular APAP concentration increased with its administration in the constructed model (Figure 5). After 2 h, APAP was eliminated by three APAP metabolic reactions: sulfation (v_Sulf), oxidation (v_CYP450), and glucuronidation (v_UGT) (Figure 1). The rate of APAP sulfation (v_Sulf) was higher in the periportal region than that in the pericentral region, while v_CYP450 and v_UGT showed the opposite tendencies. A prominent increase of protein adducts leading to hepatotoxicity was observed in the pericentral region due to low GSH concentration and sulfation activity. Indeed, over-administration of APAP resulted in GSH depletion. Overall, our model was consistent with the existing observations in the pericentral region (Yoon et al., 2016; Ahn et al., 2019; Means and Ho, 2019).

FIGURE 5

FIGURE 5. Simulation of acetaminophen (APAP) metabolism. APAP was administered for the initial 30 min of the simulation. The red and blue lines indicate the periportal and pericentral regions, respectively. The solid and dotted lines indicate moderate (0.5 mM/h) and excessive (6.0 mM/h) APAP administration, respectively. The definitions of the abbreviations are listed in Supplementary Material.

Finally, the contribution of each APAP metabolic reaction was investigated (Table 1). According to a previous study (Mcgill and Jaeschke, 2013), APAP is eliminated mainly through glucuronidation (50%–70%), followed by sulfation (25%–35%) and oxidation (5%–15%). In our model, glucuronidation was also the main route of APAP elimination in most cases, and the range of contribution of each metabolic reaction was consistent with previous studies (Mcgill and Jaeschke, 2013) (Table 1). Sulfation and oxidation increased, while glucuronidation decreased following excessive administration of APAP (6.0 mM/h). Increasing protein adduct production in the pericentral region under these conditions is indicative of site-specific hepatotoxicity. Taken together, our liver metabolism model was validated through reproduced contributions of main APAP metabolic reactions and localized toxicity previously observed in the liver.

TABLE 1

TABLE 1. Relative contribution of APAP metabolic pathways in healthy adults.

Prediction

It is widely known that obese and healthy individuals have different enzyme activities. For example, PEPCK and G6P activities are increased in obese individuals (Ropelle et al., 2009; Wang et al., 2011). Thus, we increased the enzyme activities (V_max) of PEPCK and G6P by 50% in our model and investigated the effects on APAP metabolism. Interestingly, increased PEPCK and G6P activities influenced APAP metabolism in only two cases out of eight (bold letters in Table 2; Figure 6): both cases were in the periportal region in the fasting state. In both cases, glucose-1-phosphate, UDP-glucose, and UDP-GA (metabolites connecting glycolysis and APAP metabolism) were depleted. UDP-GA is a substrate of glucuronidation; therefore, APAP could not be metabolized by this pathway. Decreased glucuronidation in the periportal region led to increased sulfation and oxidation by CYP2E1, suggesting the potential risk of increased hepatotoxicity in obese individuals.

TABLE 2

TABLE 2. Relative contribution of APAP metabolic pathways in obese patients.

FIGURE 6

FIGURE 6. Visual representation of the results in Tables 1, 2, and 3. Acetaminophen (APAP) was administered for the initial 30 min at 0.5 mM/h (moderate dose) or 6.0 mM/h (overdose) under either 4 mM (fasting) or 11 mM (feeding) of blood glucose. The reaction rates of APAP metabolic pathways were integrated up to 5 h, and the results are shown as relative contributions to APAP metabolism in pie charts. In each chart, the outer and inner circles represent the results in the periportal and pericentral regions, respectively. The activities of phosphoenolpyruvate carboxykinase (PEPCK) and glucose-6-phosphate phosphatase (G6P) increased in obese adults compared with healthy adults. Furthermore, the glycogen synthase activity increased, and PEPCK and glycogen phosphorylase activities decreased in healthy infants compared with healthy adults. The modified enzymatic activities are highlighted in green on the network maps. Dashed squares (red to red and blue to blue) highlight the prominent differences in relative contributions of APAP metabolic pathways among three model cases.

Infants have low PEPCK and glycogen phosphorylase activities and high GS activity compared with adults (Mitanchez, 2007). Sulfation is the major conjugation pathway in drug metabolism in children, while glucuronidation is dominant in adults (Jiang et al., 2013). We simulated the infant glucose metabolism by increasing the enzyme activities of PEPCK and glycogen phosphorylase and decreasing the GS activity using our liver metabolism model to investigate the possible mechanistic interconnection behind these two observations (Figure 6). There was decreased glucuronidation and increased sulfation in all cases in infants compared with adults (Tables 1 and 3). The impact was the most prominent in the state of feeding (Figure 6). For example, adults mainly used glucuronidation (66.3%) for APAP metabolism in the pericentral region with excessive APAP administration (Table 1). In contrast, glucuronidation decreased to 55.6% of the APAP metabolism in infants, while both sulfation and oxidation increased as compensation. This resulted in a drastic increase in protein adducts from 6.5% to 14.2%. These results showed that infants were susceptible to APAP-induced hepatotoxicity more than adults because of the reduced glucuronidation activity and increased oxidation.

TABLE 3

TABLE 3. Relative contribution of APAP metabolic pathways in healthy infants.

Discussion

In this study, we integrated multiple metabolic pathways including drug and glucose metabolism into a single numerical liver zonation model and validated it with the results and observations of previous reports in a semi-quantitative manner. Then, we investigated the possible crosstalk between those pathways. The model showed the effects of glucose metabolism on drug metabolism in different contexts, which indicated the importance of integrating multiple metabolic pathways.

The glucose metabolism model consisting of glycolysis, gluconeogenesis, and glycogen metabolism was determined by adopting and modifying previous rate equations (Bulik et al., 2016; Berndt et al., 2018; Berndt and Holzhütter, 2018). APAP metabolism was determined by modifying a previous model (Means and Ho, 2019). Our model also contained Cys metabolism. The combined model was validated by initially testing it without APAP administration at different blood glucose levels representing fasting and feeding states. Overall, the glucose metabolism metabolite concentrations and reaction rates in our model were consistent with existing observations of liver zonation (Figures 2–4) (Halpern et al., 2017; Halpern et al., 2018; Ben-Moshe et al., 2019). Drug metabolism in the constructed model was validated by monitoring metabolite concentrations and reaction rates of APAP metabolic pathways under moderate- and excessive APAP administration in fasting and feeding states (Figure 5). Besides, the relative contribution of each APAP metabolic reaction was investigated (Table 1). High reaction rates and contributions of glucuronidation and oxidation were observed in the pericentral region, while sulfation was promoted in the periportal region. Zonal characteristics and the contribution ranges of each metabolic reaction were similar to a previous study (Mcgill and Jaeschke, 2013). Besides, a prominent increase of protein adducts leading to hepatotoxicity was observed in the pericentral region due to GSH depletion. Altogether, our model reproduced the main contribution of APAP metabolic reactions and localized toxicity observed in the liver (Yoon et al., 2016; Ahn et al., 2019; Means and Ho, 2019).

The model showed that drug metabolism was affected by nutritional status (Table 1). The nutritional status changed the concentrations of UDP-glucose and UDP-GA and the GS reaction rate (Figures 2, 3). The low concentration of UDP-GA in the feeding state seemed to slow down glucuronidation. Decreased glucuronidation in the feeding state under moderate APAP administration was compensated by sulfation and oxidation followed by GSH conjugation and did not lead to increased protein adducts. Sulfation and GSH conjugation were promoted during excessive APAP administration; however, protein adducts increased in the feeding state, especially in the pericentral region. In general, drug administration during the feeding state is considered safer than that during the fasting state, mostly because the rate of drug administration through the digestive system decreases during feeding. However, our model suggested that the feeding state could have adverse effects by altering liver metabolism and result in higher hepatotoxicity, provided that other things are equal. This result may be slightly controversial; however, it at least suggests the necessity of careful interpretation of in vivo observation and demonstrates the advantage of simultaneously considering multiple metabolic pathways.

Our integration model was also useful for proposing hypothetical insights into the mechanisms of metabolic features involving glucose- and drug metabolism. Increased activities of PEPCK and G6P in glucose metabolism are reported in obese patients (Ropelle et al., 2009; Wang et al., 2011). Meanwhile, CYP2E1-mediated oxidation and sulfation in drug metabolism are higher than that of non-obese individuals (Brill et al., 2012; Van Rongen et al., 2016). We found that the contributions of CYP2E1-mediated oxidation and sulfation became higher in obese patients than that in healthy adults following increasing PEPCK- and G6P activities to reproduce glucose metabolism (Tables 1 and 2; Figure 6). This suggested that the increased activities of PEPCK and G6P could at least partly account for increased CYP2E1-mediated oxidation and sulfation in obese patients. However, there was a discrepancy in the contribution of glucuronidation between our model and reported observations (Brill et al., 2012; Van Rongen et al., 2016). Therefore, the model could be improved by considering activity changes in drug metabolism in obese patients, and the amount of glycogen stored in the liver may need to be adjusted for obese patients. On the contrary, these results suggest the difficulty of mechanism-based interpretations of in vivo observations and the complementarity of in vivo and in silico approaches. We also investigated the case of healthy infants in our model by adjusting some of the glucose metabolism enzymatic activities to reproduce healthy infant characteristics according to previous reports (Mitanchez, 2007). The adjustments affected drug metabolism, which successfully recapitulated the difference in relative contributions of glucuronidation and sulfation between adults and infants (Jiang et al., 2013). Although enzymatic activities in drug metabolism can have a larger impact on drug metabolism, our model suggested that the crosstalk between glucose and drug metabolism should not be neglected for the precise assessment of drug kinetics.

The study has some limitations. The model successfully reproduced periportal and pericentral characteristics of metabolic zonation regarding many metabolites and reaction rates in APAP and glucose metabolisms. Nonetheless, such validations were performed by comparing the present results with the results and observations of previous studies in a semi-quantitative manner, and we performed these validations using a limited number of conditions regarding glucose and APAP blood concentrations. The model contains only 13 zonation-dependent kinetic parameters (Supplementary Material), and metabolic zonation was simulated by manually changing enzyme activities depending on the zone that the cell belongs to, either the periportal or pericentral region. Conversely, the actual zonation is gradual and continuous along the sinusoidal axis without any histological borders and involves crosstalks of various metabolisms (Panday et al., 2022). This model also lacked the dynamic interaction of hepatocytes between the periportal and the pericentral regions through sinusoidal blood flow (Berndt and Holzhütter, 2018; Kennedy et al., 2019) and omitted the regulatory effects of oxygen and growth factors that are critical for maintaining metabolic zonation (Schmierer et al., 2010; Kietzmann, 2017; Scheidecker et al., 2020). Furthermore, there was a discrepancy between our model and the reported in vivo observations in drug metabolism in obese patients. To overcome these limitations, the modification of boundary conditions and parameters in drug metabolism based on actual measurements may lead to better prediction. The contribution of fatty acid metabolism is another challenge, especially for the better reproduction of energy metabolism in obese patients (Aubert et al., 2011; Berndt et al., 2021). To understand the regulatory mechanisms, the effects of hypoxia-inducible factor and coenzymes such as ATP may need to be considered as dependent variables of blood oxygen levels (Ainscow and Brand, 1999; Maher et al., 2007; Rankin et al., 2009). Considering future perspectives, the model may contribute to drug development and personalized medicine. The incorporation of multiscale regulatory networks ranging from gene expression to metabolisms will facilitate the complementary use of such models in in vivo and in vitro experiments, contributing to a better understanding and prediction of complex biological systems (Brunak et al., 2020; Waltemath et al., 2020).

To conclude, we integrated multiple metabolic zonation models into a single numerical human liver model and validated it with the results and observations of previous reports in a semi-quantitative manner. This model showed that changes in the enzyme activities of glucose metabolism affected APAP metabolism, and the crosstalk effects were dependent on liver zonation and nutritional status. These results showed the importance of incorporating multiple metabolic networks in a single zonation model.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

YS and MN conceived the study. KM, SH, and MS constructed the simulation model. KM, SH, MS, and MN wrote the manuscript and prepared the figures. MS, SY, and MN advised and assisted with the research plan and organization. KM, MS, YS, and MN secured funding. All authors contributed to drafting and revising the manuscript. All authors read and approved the submitted version.

Funding

This work was supported by the JSPS KAKENHI, Grant-in-Aid for Transformative Research Areas B (grant number: 20B205), and Grant-in-Aid for Scientific Research (A) (grant number: 22H00595).

Acknowledgments

We thank Kurebayahi at University of Tokyo for fruitful discussions.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2022.995597/full#supplementary-material

Abbreviation

APAP, Acetaminophen; CYP2E1, cytochrome P450 2E1; Cys, cysteine; G6P, glucose-6-phosphate phosphatase; GS, glycogen synthase; GSH, glutathione; MFA, metabolic flux analysis; NAPQI, N-acetyl-p-benzoquinone imine; PEPCK, phosphoenolpyruvate carboxykinase; UDP, uridine diphosphate; UDP-GA, UDP-glucuronic acid.

References

Abdel-Misih, S. R., and Bloomston, M. (2010). Liver anatomy. Surg. Clin. North Am. 90, 643–653. doi:10.1016/j.suc.2010.04.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Adeva-Andany, M. M., Pérez-Felpete, N., Fernández-Fernández, C., Donapetry-García, C., and Pazos-García, C. (2016). Liver glucose metabolism in humans. Biosci. Rep. 36, e00416. doi:10.1042/bsr20160385

PubMed Abstract | CrossRef Full Text | Google Scholar

Ahn, J., Ahn, J. H., Yoon, S., Nam, Y. S., Son, M. Y., and Oh, J. H. (2019). Human three-dimensional in vitro model of hepatic zonation to predict zonal hepatotoxicity. J. Biol. Eng. 13, 22. doi:10.1186/s13036-019-0148-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Ainscow, E. K., and Brand, M. D. (1999). Internal regulation of ATP turnover, glycolysis and oxidative phosphorylation in rat hepatocytes. Eur. J. Biochem. 266, 737–749. doi:10.1046/j.1432-1327.1999.00856.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Anundi, I., Lähteenmäki, T., Rundgren, M., Moldeus, P., and Lindros, K. O. (1993). Zonation of acetaminophen metabolism and cytochrome P450 2E1-mediated toxicity studied in isolated periportal and perivenous hepatocytes. Biochem. Pharmacol. 45, 1251–1259. doi:10.1016/0006-2952(93)90277-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Aubert, J., Begriche, K., Knockaert, L., Robin, M. A., and Fromenty, B. (2011). Increased expression of cytochrome P450 2E1 in nonalcoholic fatty liver disease: Mechanisms and pathophysiological role. Clin. Res. Hepatol. Gastroenterol. 35, 630–637. doi:10.1016/j.clinre.2011.04.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Ben-Moshe, S., Shapira, Y., Moor, A. E., Manco, R., Veg, T., Bahar Halpern, K., et al. (2019). Spatial sorting enables comprehensive characterization of liver zonation. Nat. Metab. 1, 899–911. doi:10.1038/s42255-019-0109-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Ben-Shachar, R., Chen, Y., Luo, S., Hartman, C., Reed, M., and Nijhout, H. F. (2012). The biochemistry of acetaminophen hepatotoxicity and rescue: A mathematical model. Theor. Biol. Med. Model. 9, 55. doi:10.1186/1742-4682-9-55

PubMed Abstract | CrossRef Full Text | Google Scholar

Berndt, N., and Holzhütter, H. G. (2018). Dynamic metabolic zonation of the hepatic glucose metabolism is accomplished by sinusoidal plasma gradients of nutrients and hormones. Front. Physiol. 9, 1786. doi:10.3389/fphys.2018.01786

PubMed Abstract | CrossRef Full Text | Google Scholar

Berndt, N., Horger, M. S., Bulik, S., and Holzhütter, H. G. (2018). A multiscale modelling approach to assess the impact of metabolic zonation and microperfusion on the hepatic carbohydrate metabolism. PLoS Comput. Biol. 14, e1006005. doi:10.1371/journal.pcbi.1006005

PubMed Abstract | CrossRef Full Text | Google Scholar

Berndt, N., Kolbe, E., Gajowski, R., Eckstein, J., Ott, F., Meierhofer, D., et al. (2021). Functional consequences of metabolic zonation in murine livers: Insights for an old story. Hepatology 73, 795–810. doi:10.1002/hep.31274

PubMed Abstract | CrossRef Full Text | Google Scholar

Brill, M. J., Diepstraten, J., Van Rongen, A., Van Kralingen, S., Van Den Anker, J. N., and Knibbe, C. A. (2012). Impact of obesity on drug metabolism and elimination in adults and children. Clin. Pharmacokinet. 51, 277–304. doi:10.2165/11599410-000000000-00000

PubMed Abstract | CrossRef Full Text | Google Scholar

Brunak, S., Bjerre Collin, C., Eva Ó Cathaoir, K., Golebiewski, M., Kirschner, M., Kockum, I., et al. (2020). Towards standardization guidelines for in silico approaches in personalized medicine. J. Integr. Bioinform. 17, 20200006. doi:10.1515/jib-2020-0006

PubMed Abstract | CrossRef Full Text | Google Scholar

Bulik, S., Holzhütter, H. G., and Berndt, N. (2016). The relative importance of kinetic mechanisms and variable enzyme abundances for the regulation of hepatic glucose metabolism--insights from mathematical modeling. BMC Biol. 14, 15. doi:10.1186/s12915-016-0237-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Court, M. H., Duan, S. X., Von Moltke, L. L., Greenblatt, D. J., Patten, C. J., Miners, J. O., et al. (2001). Interindividual variability in acetaminophen glucuronidation by human liver microsomes: Identification of relevant acetaminophen UDP-glucuronosyltransferase isoforms. J. Pharmacol. Exp. Ther. 299, 998

PubMed Abstract | Google Scholar

Cunningham, R. P., and Porat-Shliom, N. (2021). Liver zonation - revisiting old questions with new technologies. Front. Physiol. 12, 732929. doi:10.3389/fphys.2021.732929

PubMed Abstract | CrossRef Full Text | Google Scholar

Dargue, R., Zia, R., Lau, C., Nicholls, A. W., Dare, T. O., Lee, K., et al. (2020). Metabolism and effects on endogenous metabolism of paracetamol (acetaminophen) in a porcine model of liver failure. Toxicol. Sci. 175, 87–97. doi:10.1093/toxsci/kfaa023

PubMed Abstract | CrossRef Full Text | Google Scholar

Franiatte, S., Clarke, R., and Ho, H. (2019). A computational model for hepatotoxicity by coupling drug transport and acetaminophen metabolism equations. Int. J. Numer. Method. Biomed. Eng. 35, e3234. doi:10.1002/cnm.3234

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghafoory, S., Breitkopf-Heinlein, K., Li, Q., Scholl, C., Dooley, S., and Wölfl, S. (2013). Zonation of nitrogen and glucose metabolism gene expression upon acute liver damage in mouse. PLoS One 8, e78262. doi:10.1371/journal.pone.0078262

PubMed Abstract | CrossRef Full Text | Google Scholar

Halpern, K. B., Shenhav, R., Massalha, H., Toth, B., Egozi, A., Massasa, E. E., et al. (2018). Paired-cell sequencing enables spatial gene expression mapping of liver endothelial cells. Nat. Biotechnol. 36, 962–970. doi:10.1038/nbt.4231

PubMed Abstract | CrossRef Full Text | Google Scholar

Halpern, K. B., Shenhav, R., Matcovitch-Natan, O., Toth, B., Lemze, D., Golan, M., et al. (2017). Single-cell spatial reconstruction reveals global division of labour in the mammalian liver. Nature 542, 352–356. doi:10.1038/nature21065

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, X. L., Zhao, P., Barrett, J. S., Lesko, L. J., and Schmidt, S. (2013). Application of physiologically based pharmacokinetic modeling to predict acetaminophen metabolism and pharmacokinetics in children. CPT. Pharmacometrics Syst. Pharmacol. 2, e80. doi:10.1038/psp.2013.55

PubMed Abstract | CrossRef Full Text | Google Scholar

Katz, N., Teutsch, H. F., Sasse, D., and Jungermann, K. (1977). Heterogeneous distribution of glucose-6-phosphatase in microdissected periportal and perivenous rat liver tissue. FEBS Lett. 76, 226–230. doi:10.1016/0014-5793(77)80157-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Kennedy, R. C., Smith, A. K., Ropella, G. E. P., Mcgill, M. R., Jaeschke, H., and Hunt, C. A. (2019). Propagation of pericentral necrosis during acetaminophen-induced liver injury: Evidence for early interhepatocyte communication and information exchange. Toxicol. Sci. 169, 151–166. doi:10.1093/toxsci/kfz029

PubMed Abstract | CrossRef Full Text | Google Scholar

Kietzmann, T. (2017). Metabolic zonation of the liver: The oxygen gradient revisited. Redox Biol. 11, 622–630. doi:10.1016/j.redox.2017.01.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Lammers, L. A., Achterbergh, R., Mathôt, R. a. A., and Romijn, J. A. (2020). The effects of fasting on drug metabolism. Expert Opin. Drug Metab. Toxicol. 16, 79–85. doi:10.1080/17425255.2020.1706728

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, R., Martínez-Ramírez, A. S., Borders, T. L., Gao, F., and Sosa-Pineda, B. (2020). Metabolic and non-metabolic liver zonation is established non-synchronously and requires sinusoidal Wnts. Elife 9, e46206. doi:10.7554/eLife.46206

PubMed Abstract | CrossRef Full Text | Google Scholar

Maher, J. C., Wangpaichitr, M., Savaraj, N., Kurtoglu, M., and Lampidis, T. J. (2007). Hypoxia-inducible factor-1 confers resistance to the glycolytic inhibitor 2-deoxy-D-glucose. Mol. Cancer Ther. 6, 732–741. doi:10.1158/1535-7163.Mct-06-0407

PubMed Abstract | CrossRef Full Text | Google Scholar

Mazaleuskaya, L. L., Sangkuhl, K., Thorn, C. F., Fitzgerald, G. A., Altman, R. B., and Klein, T. E. (2015). PharmGKB summary: Pathways of acetaminophen metabolism at the therapeutic versus toxic doses. Pharmacogenet. Genomics 25, 416–426. doi:10.1097/fpc.0000000000000150

PubMed Abstract | CrossRef Full Text | Google Scholar

Mcgill, M. R., and Jaeschke, H. (2013). Metabolism and disposition of acetaminophen: Recent advances in relation to hepatotoxicity and diagnosis. Pharm. Res. 30, 2174–2187. doi:10.1007/s11095-013-1007-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Means, S. A., and Ho, H. (2019). A spatial-temporal model for zonal hepatotoxicity of acetaminophen. Drug Metab. Pharmacokinet. 34, 71–77. doi:10.1016/j.dmpk.2018.09.266

PubMed Abstract | CrossRef Full Text | Google Scholar

Méndez-Lucas, A., Duarte, J. A., Sunny, N. E., Satapati, S., He, T., Fu, X., et al. (2013). PEPCK-M expression in mouse liver potentiates, not replaces, PEPCK-C mediated gluconeogenesis. J. Hepatol. 59, 105–113. doi:10.1016/j.jhep.2013.02.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Miethke, H., Wittig, B., Nath, A., Zierz, S., and Jungermann, K. (1985). Metabolic zonation in liver of diabetic rats. Zonal distribution of phosphoenolpyruvate carboxykinase, pyruvate kinase, glucose-6-phosphatase and succinate dehydrogenase. Biol. Chem. Hoppe. Seyler. 366, 493–501. doi:10.1515/bchm3.1985.366.1.493

PubMed Abstract | CrossRef Full Text | Google Scholar

Mitanchez, D. (2007). Glucose regulation in preterm newborn infants. Horm. Res. 68, 265–271. doi:10.1159/000104174

PubMed Abstract | CrossRef Full Text | Google Scholar

Mutlib, A. E., Goosen, T. C., Bauman, J. N., Williams, J. A., Kulkarni, S., and Kostrubsky, S. (2006). Kinetics of acetaminophen glucuronidation by UDP-glucuronosyltransferases 1A1, 1A6, 1A9 and 2B15. Potential implications in acetaminophen-induced hepatotoxicity. Chem. Res. Toxicol. 19, 701–709. doi:10.1021/tx050317i

PubMed Abstract | CrossRef Full Text | Google Scholar

Panday, R., Monckton, C. P., and Khetani, S. R. (2022). The role of liver zonation in physiology, regeneration, and disease. Semin. Liver Dis. 42, 1–16. doi:10.1055/s-0041-1742279

PubMed Abstract | CrossRef Full Text | Google Scholar

Rankin, E. B., Rha, J., Selak, M. A., Unger, T. L., Keith, B., Liu, Q., et al. (2009). Hypoxia-inducible factor 2 regulates hepatic lipid metabolism. Mol. Cell. Biol. 29, 4527–4538. doi:10.1128/mcb.00200-09

PubMed Abstract | CrossRef Full Text | Google Scholar

Reddyhoff, D., Ward, J., Williams, D., Regan, S., and Webb, S. (2015). Timescale analysis of a mathematical model of acetaminophen metabolism and toxicity. J. Theor. Biol. 386, 132–146. doi:10.1016/j.jtbi.2015.08.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Reith, D., Medlicott, N. J., Kumara De Silva, R., Yang, L., Hickling, J., and Zacharias, M. (2009). Simultaneous modelling of the Michaelis-Menten kinetics of paracetamol sulphation and glucuronidation. Clin. Exp. Pharmacol. Physiol. 36, 35–42. doi:10.1111/j.1440-1681.2008.05029.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Remien, C. H., Adler, F. R., Waddoups, L., Box, T. D., and Sussman, N. L. (2012). Mathematical modeling of liver injury and dysfunction after acetaminophen overdose: Early discrimination between survival and death. Hepatology 56, 727–734. doi:10.1002/hep.25656

PubMed Abstract | CrossRef Full Text | Google Scholar

Riches, Z., Bloomer, J., Patel, A., Nolan, A., and Coughtrie, M. (2009). Assessment of cryopreserved human hepatocytes as a model system to investigate sulfation and glucuronidation and to evaluate inhibitors of drug conjugation. Xenobiotica. 39, 374–381. doi:10.1080/00498250902763440

PubMed Abstract | CrossRef Full Text | Google Scholar

Ropelle, E. R., Pauli, J. R., Cintra, D. E., Frederico, M. J., De Pinho, R. A., Velloso, L. A., et al. (2009). Acute exercise modulates the Foxo1/PGC-1alpha pathway in the liver of diet-induced obesity rats. J. Physiol. 587, 2069–2076. doi:10.1113/jphysiol.2008.164202

PubMed Abstract | CrossRef Full Text | Google Scholar

Sacks, D., Baxter, B., Campbell, B. C. V., Carpenter, J. S., Cognard, C., Dippel, D., et al. (2018). Multisociety consensus quality improvement revised consensus statement for endovascular therapy of acute ischemic stroke. Int. J. Stroke 13, 612–632. doi:10.1177/1747493018778713

PubMed Abstract | CrossRef Full Text | Google Scholar

Scheidecker, B., Shinohara, M., Sugimoto, M., Danoy, M., Nishikawa, M., and Sakai, Y. (2020). Induction of in vitro metabolic zonation in primary hepatocytes requires both near-physiological oxygen concentration and flux. Front. Bioeng. Biotechnol. 8, 524. doi:10.3389/fbioe.2020.00524

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmierer, B., Novák, B., and Schofield, C. J. (2010). Hypoxia-dependent sequestration of an oxygen sensor by a widespread structural motif can shape the hypoxic response--a predictive kinetic model. BMC Syst. Biol. 4, 139. doi:10.1186/1752-0509-4-139

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, A. K., Petersen, B. K., Ropella, G. E., Kennedy, R. C., Kaplowitz, N., Ookhtens, M., et al. (2016). Competing mechanistic hypotheses of acetaminophen-induced hepatotoxicity challenged by virtual experiments. PLoS Comput. Biol. 12, e1005253. doi:10.1371/journal.pcbi.1005253

PubMed Abstract | CrossRef Full Text | Google Scholar

Stark, R., Guebre-Egziabher, F., Zhao, X., Feriod, C., Dong, J., Alves, T. C., et al. (2014). A role for mitochondrial phosphoenolpyruvate carboxykinase (PEPCK-M) in the regulation of hepatic gluconeogenesis. J. Biol. Chem. 289, 7257–7263. doi:10.1074/jbc.C113.544759

PubMed Abstract | CrossRef Full Text | Google Scholar

Stark, R., and Kibbey, R. G. (2014). The mitochondrial isoform of phosphoenolpyruvate carboxykinase (PEPCK-M) and glucose homeostasis: Has it been overlooked? Biochim. Biophys. Acta 1840, 1313–1330. doi:10.1016/j.bbagen.2013.10.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Trefts, E., Gannon, M., and Wasserman, D. H. (2017). The liver. Curr. Biol. 27, R1147–R1151. doi:10.1016/j.cub.2017.09.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Rongen, A., Välitalo, P. a. J., Peeters, M. Y. M., Boerma, D., Huisman, F. W., Van Ramshorst, B., et al. (2016). Morbidly obese patients exhibit increased CYP2E1-mediated oxidation of acetaminophen. Clin. Pharmacokinet. 55, 833–847. doi:10.1007/s40262-015-0357-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Waltemath, D., Golebiewski, M., Blinov, M. L., Gleeson, P., Hermjakob, H., Hucka, M., et al. (2020). The first 10 years of the international coordination network for standards in systems and synthetic biology (COMBINE). J. Integr. Bioinform. 17, 20200005. doi:10.1515/jib-2020-0005

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Nakagawa, Y., Liu, L., Wang, W., Ren, X., Anghel, A., et al. (2011). Tissue-specific dysregulation of hexose-6-phosphate dehydrogenase and glucose-6-phosphate transporter production in db/db mice as a model of type 2 diabetes. Diabetologia 54, 440–450. doi:10.1007/s00125-010-1956-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoon, E., Babar, A., Choudhary, M., Kutner, M., and Pyrsopoulos, N. (2016). Acetaminophen-induced hepatotoxicity: A comprehensive update. J. Clin. Transl. Hepatol. 4, 131–142. doi:10.14218/jcth.2015.00052

PubMed Abstract | CrossRef Full Text | Google Scholar

Yugi, K., Nakayama, Y., Kinoshita, A., and Tomita, M. (2005). Hybrid dynamic/static method for large-scale simulation of metabolism. Theor. Biol. Med. Model. 2, 42. doi:10.1186/1742-4682-2-42

PubMed Abstract | CrossRef Full Text | Google Scholar