A Spatial Model of Hepatic Calcium Signaling and Glucose Metabolism Under Autonomic Control Reveals Functional Consequences of Varying Liver Innervation Patterns Across Species

Rapid breakdown of hepatic glycogen stores into glucose plays an important role during intense physical exercise to maintain systemic euglycemia. Hepatic glycogenolysis is governed by several different liver-intrinsic and systemic factors such as hepatic zonation, circulating catecholamines, hepatocellular calcium signaling, hepatic neuroanatomy, and the central nervous system (CNS). Of the factors regulating hepatic glycogenolysis, the extent of lobular innervation varies significantly between humans and rodents. While rodents display very few autonomic nerve terminals in the liver, nearly every hepatic layer in the human liver receives neural input. In the present study, we developed a multi-scale, multi-organ model of hepatic metabolism incorporating liver zonation, lobular scale calcium signaling, hepatic innervation, and direct and peripheral organ-mediated communication between the liver and the CNS. We evaluated the effect of each of these governing factors on the total hepatic glucose output and zonal glycogenolytic patterns within liver lobules during simulated physical exercise. Our simulations revealed that direct neuronal stimulation of the liver and an increase in circulating catecholamines increases hepatic glucose output mediated by mobilization of intracellular calcium stores and lobular scale calcium waves. Comparing simulated glycogenolysis between human-like and rodent-like hepatic innervation patterns (extensive vs. minimal) suggested that propagation of calcium transients across liver lobules acts as a compensatory mechanism to improve hepatic glucose output in sparsely innervated livers. Interestingly, our simulations suggested that catecholamine-driven glycogenolysis is reduced under portal hypertension. However, increased innervation coupled with strong intercellular communication can improve the total hepatic glucose output under portal hypertension. In summary, our modeling and simulation study reveals a complex interplay of intercellular and multi-organ interactions that can lead to differing calcium dynamics and spatial distributions of glycogenolysis at the lobular scale in the liver.

Introduction

The liver plays a critical role in mammalian physiology by regulating a wide range of metabolic processes including the balance of blood glucose levels. With its ability to control glucose uptake from and release into circulating blood by tuning the activity of appropriate enzymes, the liver is able to modulate glycogenolysis, gluconeogenesis, glycogenesis, and glycolysis – the major pathways involved in glucose metabolism in the liver – thereby maintaining short and longer-term glucose homeostasis (Nordlie et al., 1999; Radziuk and Pye, 2001; Wahren and Ekberg, 2007; Han et al., 2016). Disruption of glucose metabolism in the liver can lead to chronic disease as is well-known in the association of non-alcoholic fatty liver disease (NAFLD) with metabolic syndrome and type II diabetes (Jiang and Zhang, 2003; Parekh and Anania, 2007; Chao et al., 2019).

The mammalian liver serves as a reservoir of glycogen, which can be broken down by hepatocytes into glucose and released into the circulating blood. Hepatic glycogenolysis is primarily regulated by glucagon secreted by the pancreas. In response to increased demand of circulating glucose, such as during exercise, synthesis and secretion of glucagon increases by autonomic activation and signaling through the splanchnic nerve. Increased circulating glucagon induces rapid glycogenolysis in the liver. In addition to glucagon, glycogenolysis in the liver can be modulated by catecholamines during periods of starvation or intense physical activity (Chu et al., 1996). Release of epinephrine into the circulating blood supply is driven by the central nervous system (CNS) sensing of blood glucose levels and subsequent autonomic signaling to the adrenal glands through the activation of splanchnic nerve fibers (Livingstone et al., 1995; Burcelin et al., 2000; Marty et al., 2007; Oosterveer and Kristina, 2014). Activation of adrenergic receptors on hepatocytes by circulating epinephrine leads to glycogenolysis-inducing downstream signaling (Sherwin and Sacca, 1984). Additionally, activation of the adrenergic receptors elicits calcium spiking within hepatocytes through a Gq protein coupled receptor signaling pathway (Orellana et al., 1987; Baukal et al., 1988; Wolf et al., 1988; Gaspers et al., 2019). Interestingly, glycogenolysis is a calcium-dependent process. An increase in cytosolic calcium can enhance the rate of glycogenolysis through the activation of glycogen phosphorylase kinase (Assimacopoulos-Jeannet et al., 1977; Billah and Michell, 1979). Glycogenolysis can be affected similarly by the release of norepinephrine at the autonomic nerve endings within the liver, representing a direct mode of CNS control of hepatic glucose metabolism (Exton, 1987). Catecholamines can thus initiate glucose release and simultaneously accelerate its production to facilitate quick response of the liver to organismal metabolic demands.

Although all hepatocytes can break down glycogen, there is heterogeneity in the extent of glycogenolysis within hepatocytes (Jungermann, 1983). This heterogeneity is driven by a combination of several factors. The liver consists of repeating hexagonal units, termed lobules, each of which is organized around a central vein. Vertices of the hexagonal lobules are occupied by portal triads, containing the portal vein, bile duct, and hepatic artery. Hepatocytes are aligned along sinusoids that connect the portal region of lobules to the central vein. Blood flows through the sinusoids from the periportal to the pericentral region. The directionality of blood flow through the lobule creates a natural gradient of circulating hormonal stimulus to hepatocytes lined along the sinusoids, which translates into varying levels of activation of intracellular signaling and metabolic processes (Jungermann and Keitzmann, 1996). The inherent zonation of hepatocellular metabolic processes and cell-intrinsic calcium signaling regulators contributes to the heterogeneity of glycogenolytic activity of hepatocytes. Expression of glycogenolytic genes has been shown to be more active in the periportal region of liver lobules, with a progressive decrease toward the pericentral region (Jungermann, 1988; Jungermann and Kietzmann, 2000). In contrast, calcium “waves” in the liver lobule, driven by intercellular exchange of inositol 1,4,5-trisphosphate (IP₃), travel from the pericentral zone to the periportal zone (Amaya and Nathanson, 2013). Such intercellular calcium “waves” are responsible for integrating hormonal control of glucose output in the intact liver, as demonstrated by Gaspers et al. (2019). The opposing gradients of cell intrinsic calcium signaling parameters and circulating hormones can combine to influence zonal contributions toward total glucose output within lobules.

An additional layer of hepatic glycogenolysis control is provided by hepatic innervation. The liver is innervated by both sympathetic fibers from the splanchnic nerve and parasympathetic nerve fibers under the control of the vagus nerve (Berthoud, 2004; McCuskey, 2004; Yi et al., 2010; Jensen et al., 2013; Mizuno and Ueno, 2016). Interestingly, the extent of hepatic innervation varies significantly across species. At the liver lobule scale, very few nerve fibers make contact with hepatocytes in rodents, limited only to the periportal region. On the other hand, in humans, nearly every hepatic lobular layer receives neural input (Figure 1A; Reilly et al., 1978; Metz and Forssmann, 1980; Akiyoshi et al., 1998). As a result of autonomic stimulation, norepinephrine can be released at synapses between hepatic neurons and hepatocytes. Although norepinephrine does not play a significant role in glycogenolysis, it can increase glycogenolysis rates in hepatocytes indirectly by inducing Ca²⁺ oscillations in hepatocytes (Yamaguchi and Garceau, 1980; Jungermann et al., 1987). Innervation patterns of the liver can thus modulate zonal contributions to hepatic glucose synthesis through induction of Ca²⁺ signaling. Yet, it has been shown that both rats and humans promote glycogenolysis through similar mechanisms of activation (Amaya and Nathanson, 2013). This raises an interesting question: How are rodents able to compensate for the sparsity of direct neural input to liver parenchyma? One hypothesis postulates that the extent of hepatic innervation is inversely proportional to the amount of gap junctions across species. This hypothesis stems from the work of Forssmann and Ito (1977) and Reilly et al. (1978) in which they show that the density of gap junctions is low in the tree shrew and the guinea pig but high in rats. Similar to humans, tree shrews and guinea pigs exhibit nerve contacts with nearly all hepatocytes (Forssmann and Ito, 1977; Reilly et al., 1978; Metz and Forssmann, 1980). This potential compensatory relationship between innervation and gap junction connectivity remains unresolved.

FIGURE 1

Figure 1. Schematic representation of the multi-organ model of glucose metabolism and calcium signaling under autonomic control. (A) Hepatic innervation patterns across species. (B) Lobular scale model with varying degrees of hepatic innervation tuned to species. The lobule is represented as a repeating chain of hepatocytes and sinusoidal blood compartments that interact via cross-membrane transport. Hepatocytes can freely exchange IP3 with adjacent neighbors. Blood progresses from hepatic compartment 1–15 and exits into a larger compartment representing the rest of the body. The systemic blood compartment interacts with the pancreas, adrenal glands, adipose tissue, and central nervous system (CNS). (C) Inter and intra-hepatic network model with zonated intracellular processes and regulation. (D) Multi-organ interaction network of systemic signals. Figures B and C were adapted and modified from Ashworth et al. (2016).

Since innervation is crucial for metabolic functionality of the liver, denervation, or loss of function of hepatic nerves can affect liver function adversely causing systemic issues. Functional or physical loss of hepatic innervation can occur due to different reasons. Loss of nerve function was observed in livers of streptozotocin-induced diabetic mice (Stümpel et al., 1996). Cirrhotic livers have been reported to exhibit loss of parenchymal innervation (Lee et al., 1992). Effects of loss of hepatic innervation have been widely studied in the context of insulin sensitivity (Mikines et al., 1985; Xie and Lautt, 1996). However, the effects of hepatic denervation on liver glucose output and the effect due to lobular scale calcium dynamics remains poorly understood.

In the present study, we used a multi-scale modeling approach to explore the dynamic behavior of a complex network of intracellular signaling, intercellular, and multi-organ interactions governing calcium signaling and glycogenolysis in the liver. We developed a multi-scale, multi-organ model of central regulation of hepatic glucose metabolism, in which we incorporated intracellular glucose and glycogen metabolism, lobular scale zonation of these metabolic processes, hepatic innervation, peripheral organ-mediated hormonal signaling, and tissue scale calcium signaling elicited by systemic hormones (Figures 1B–D). We built upon previously published models of hepatic glucose metabolism, intracellular calcium signaling, and lobule-scale calcium dynamics, integrating a new model of direct and indirect control of liver metabolism by the central nervous system through catecholaminergic processes. We evaluated, via simulation, a range of scenarios including species-specific differences in liver innervation, alteration of lobule-scale gap junction connectivity, and disruption of local and global neural signals. Specifically, we examined the impact of these variations on calcium-mediated regulation of liver glucose metabolism and hepatic glucose output.

Materials and Methods

Multi-Scale Model Development

Our model builds on the previously developed multi-organ model of hepatic glucose, lipids, and triglyceride metabolism from Ashworth et al. (2016). Briefly, Ashworth et al. (2016) developed an ordinary differential equations (ODE)-based model to simulate and explore zonal features of insulin resistance and hepatic steatosis in non-alcoholic fatty liver disease under periodic feeding-fasting cycles. They modeled metabolic zonation in the liver lobules by treating the porto-central axis of the sinusoid as a repeating unit, with the blood supply and associated hepatocytes split into compartments based on their positions along the lobular axis (Figure 1B). In this scheme, each sinusoidal compartment interacts with local hepatocytes, each of which is affected by hormone signals and exchanges metabolites within its own blood compartment. Blood enters the periportal compartment and flows through to the pericentral compartment where it deposits into a larger body compartment. The original Ashworth et al. (2016) model captures the extensive metabolic and hormone signaling network present in hepatocytes. A detailed description of the original model and the parametric sensitivity analysis can be found in Ashworth et al. (2016). In the present study, we extended the Ashworth et al. (2016) model to include intracellular calcium signaling and intercellular IP3 signaling models into the compartmental model of the liver lobule (Figure 1C). In addition, we included adrenal glands and the central nervous system as part of the body compartments external to the liver (Figure 1D). Our expanded multi-scale model includes additional features and modifications necessary for studying autonomic control of zonated hepatic glucose metabolism, as detailed below.

Modeling Porto-Central Axis of Human Liver

We expanded the 8-compartment Ashworth et al. (2016) model to consist of 15 compartments within a hepatic lobule to account for the extent of the porto-central axis within a human liver more accurately. In the expanded hepatic lobule model, the first and last four compartments correspond to the periportal and pericentral regions, respectively, and the middle seven compartments correspond to the midlobular region. All 15 compartments exchange molecules with their immediate neighboring blood compartments while the first and fifteenth compartments are in direct contact with systemic blood and exchange molecules within the body compartment. Molecular exchange between all compartments is modeled as mass exchange terms dependent on blood flow rates, according to the following equations:

$\begin{array}{ll}(d{M}_{(i=1:N)})/dt=bf·(M_{(i-1)}-M_i)& (1)\end{array}$

$\begin{array}{ll}(d{M}_0)/dt=(bf·(M_N-M_0))/s& (2)\end{array}$

where M_i denotes the concentration of species M in blood compartment i (i = 1 …. 15), M₀ represents the concentration of species M in systemic blood, and s represents the body to hepatic compartment ratio: s = 5 · N (s = 75 for the N = 15 compartments used). The parameter bf represents blood flow through the compartments, which is taken to be bf = 0.15·n s⁻¹ [bf = 2.25s⁻¹ for the 15 compartments used]. The index i = 1 corresponds to the periportal (PP) sinusoidal, and associated hepatocyte compartment, while the index i = 15 corresponds to pericentral (PC) compartment.

Modeling Zonated Metabolic Processes in a Liver Lobule

Molecular exchange between adjacent sinusoidal blood compartments is modeled as mass transfer kinetics. The parameters derived from experimental data of various species were used when considering molecular exchange between a sinusoidal blood compartment and associated hepatocyte. The hepatocyte metabolic network contains both zonated and non-zonated processes consisting of glucose, lipid, fatty acid, and phosphate molecules with their relevant intermediates. Since the oxygen concentration in the sinusoidal blood decreases from the PP to PC region due to uptake by hepatocytes, there exists a gradient of metabolic processes following a specific zonated profile (Jungermann and Kietzmann, 1997). The Ashworth et al. model considers zonation of hepatocyte functions to be dependent on the hepatocyte-specific oxygen concentration along the porto-central axis. The metabolic reactions with higher flux in the periportal zone in which the oxygen concentration is the highest are free fatty acid (FFA) transport between sinusoidal blood compartment and hepatocyte, conversion of FFA into palmitate, ADP to ATP conversion, pyruvate to glucose-3-phosphate (G3P) conversion, G3P to glucose-6-phosphate (G6P) conversion, G6P to glucose conversion, G6P to glycogen conversion, and glycogen to G6P conversion. The metabolic reactions with higher flux in the pericentral zone in which the oxygen concentration is the lowest are glucose to G6P conversion, G3P to pyruvate conversion, and Acetyl-Co-A to FFA conversion (Figure 1C).

Modeling Systemic Insulin and Glucagon Secretion

The model accounts for systemic insulin and glucagon secretion by the pancreas. Both hormones are metabolized during their passage through the sinusoidal compartments from PP to PC, leading to a zonated concentration profile. Rates of glucagon and insulin secretion in the model are dependent on systemic glucose concentration and the systemic glucose reference value (G_ref). If the systemic glucose concentration is less than the reference value, glucagon is secreted and if the systemic glucose concentration is greater than the reference value, insulin is secreted. The rate of glucagon (Glu) and insulin (Ins) secretion are given by the following equations and the details of parameterization can be found in Ashworth et al. (2016):

$\begin{array}{l}dGlu/dt=1/{\tau }_{glu}·\text{ln}{(G_{ref}/G_B)}^{(n_{glu})}/[{(k_m^{glu})}^{(n_{glu})}\\ +\text{ ln}{(G_{ref}/G_B)}^{(n_{glu})}],if{G}_B<G_{ref}& (3)\end{array}$

$\begin{array}{l}dIns/dt=1/{\tau }_{ins}·\text{ln}{(G_{ref}/G_B)}^{(n_{ins})}/[{(k_m^{ins})}^{(n_{ins})}\\ +\text{ ln}{(G_{ref}/G_B)}^{(n_{ins})}],if{G}_B>G_{ref}& (4)\end{array}$

where G_ref = 2 a.u., τ_glu = 1.3 s, $k_m^{glu}=1.3nM$, n_glu = 2, τ_ins = 0.9 s, $k_m^{ins}=0.9nM$, n_ins = 2.

Both triglyceride synthesis and lipolysis from the gut and adipose tissue are modeled as single lumped equations with rates modeled as Hill-like equations. Additionally, the model includes inputs of blood glucose and free fatty acids (FFAs) to simulate a daily feeding cycle. However, we assume that there is no external blood glucose or FFA input, for example, through dietary intake during the course of simulation.

Modeling CNS Activation by Systemic Glucose

We expanded the above model to include autonomic control of zonated hepatic glucose metabolism. Activation of the CNS was a key component in the expanded model, to account for central autonomic control of multiple organs for coordinating the response to a systemic glucose demand. CNS activation was modeled using systemic glucose-dependent saturation kinetics, given by the following equation:

$\begin{array}{l}CNSAct=(G_{symp}-G_B)/(K_m^{symp}\\ +(G_{symp}-G_B)),if{G}_B<G_{symp}& (5)\end{array}$

where $K_m^{symp}=500\mu M$ and G_symp = 4500 μM.

The formulation considers CNS activation as occurring only when blood glucose levels decrease below 4,500 μM, as normal blood glucose levels range between 3,500 and 5,500 μM (Güemes et al., 2015). For simplicity, we assume that blood glucose sensing occurs only at the hypothalamus, which orchestrates all peripheral hormonal signals and metabolic processes (Livingstone et al., 1995; Marty et al., 2007).

Modeling Calcium Regulation of Glucose Homeostasis

Calcium dynamics have been shown to regulate glucose dynamics in hepatocytes, as an increase in Ca²⁺ results in an increase in glucose output and glycogenolysis (Billah and Michell, 1979). This regulatory effect of Ca²⁺ is mediated by the activation of glycogen phosphorylase kinase (Assimacopoulos-Jeannet et al., 1977). Our model assumes basal activation of glycogen phosphorylase kinase and a minimum additional cytosolic Ca²⁺ amplitude-dependent activation threshold (CaI^GPK = 0.44 μM). We account for the subsequent glycogen phosphorylase activation implicitly in the glycogenolysis as follows:

$\begin{array}{ll}{v}_i^{GPK}=\{\begin{array}{l}\frac{Ca{I}_i}{K_{CaI}^{GPK}+Ca{I}_i};Ca{I}_i>Ca{I}^{GPK}\\ v_0;Ca{I}_i<Ca{I}^{GPK}\end{array}& (6)\end{array}$

$\begin{array}{l}dGlyc/dt=-(v_{brk}·K_{Phos}^{max}\\ ·{(Glyc)}^{(n_{brk})})/({(Glyc)}^{(n_{brk})}\\ +{(K_m^{Glyc})}^{(n_{brk})})·Phos/(K_m^{Phos}\\ +Phos)·(1+v_{CaI}^{max}·v_i^{GPK})& (7)\end{array}$

$\begin{array}{ll}{K}_{Phos}^{max}=(Glu+k_{LP})/(Ins+k_{IP})& (8)\end{array}$

where $K_{CaI}^{GPK}=3\mu M$, v₀ = 0.05, $K_m^{Glyc}=100mM$, $K_m^{Phos}=4000\mu M$, n_brk = 4, $v_{brk}=5\mu M·s^{-1}$, k_LP = 45 pM, k_IP = 26.66 pM. $v_{CaI}^{max}=3$ is the maximum increase in glycogen phosphorylase induced by Ca²⁺, with a Michaelis-Menten constant of $K_m^{CaI}=0.6\mu M$. Phos represents the intracellular phosphate concentration. The glycogen phosphorylase kinase activity parameter is a dimensionless quantity that modulates the rate of glycogenolysis within each hepatocyte.

Modeling Calcium Signaling Across the Porto-Central Axis

The present expanded model couples the intra-hepatocyte metabolic model developed by Ashworth et al. (2016) with a Ca²⁺ spiking model (Verma et al., 2016) based on the IP₃-Ca²⁺ cross coupling proposed by Meyer and Stryer (1991). Briefly, the Ca²⁺ spiking model accounts for extracellular stimulus-induced activation of the intracellular Ca²⁺ signaling cascade. Ca²⁺ spiking occurs due to the positive feedback by IP₃-induced Ca²⁺ release and Ca²⁺-stimulated IP₃ formation. We model the induction of Ca²⁺ signals by extracellular catecholamine (epinephrine and norepinephrine) stimulus through α₁-adrenergic receptor activation, however, similar Ca²⁺ dynamics can be induced by other extracellular stimuli such as vasopressin (Assimacopoulos-Jeannet et al., 1977; Schöfl et al., 1993; Verma et al., 2018). The equations for intracellular Ca²⁺ fluxes and IP₃ receptor dynamics are retained from Meyer and Stryer (1991) with the exception of IP₃ synthesis, which is related to the hormone receptor dynamics instead of constant receptor activation in the original Meyer and Stryer (1991) model.

The present multicellular model of Ca²⁺ signal propagation includes gap junction-mediated intracellular communication for a chain of N connected hepatocytes. In Equations 9-13 below, the subscript i represents the cell index for hepatocyte specific species and parameters. For each adjacent cell pair, IP₃ exchange is modeled by a mass exchange term with the mass transfer parameter G_ij. Although Ca²⁺ can normally pass through the gap junctions, for the simplicity of our model and the lower diffusion constant of Ca²⁺ compared with IP₃ (560 vs. 280 um²/s respectively), the model only assumes IP₃ exchange between adjacent hepatocytes. Additionally, considering the short intracellular diffusion time compared with typical spike temporal width (~30 s) and inter-spike intervals (~100 s), the IP₃ concentration is assumed to be uniform throughout the cytosol. We approximate IP₃ exchange between adjacent hepatocytes as a mass exchange term owing to the small latency period compared with typical spike widths and inter-spike intervals.

Binding of extracellular catecholamines to cognate receptors on the cell membrane leads to activation of the intracellular Ca²⁺ signaling cascade. While receptor internalization is possible, the model assumes that the total receptor number remains constant. The Ca²⁺ spiking state of every hepatocyte is characterized by four different variables: cytosolic Ca²⁺ concentration (CaI), cytosolic IP₃ (IP3), ratio of free (unbound) to total IP₃ receptors (g) on the endoplasmic reticulum surface, and the ratio of free (unbound) to total stimulus receptors on the cell surface (r). In Equation 9, R_f denotes the cell surface density of free α₁-adrenergic receptors, whereas R_total represents the total cell surface density of α₁-adrenergic receptors. The rate of increase in free receptor density is assumed to be proportional to the number of bound receptors, with rate constant k_r. The density of free receptors is assumed to decrease through two different modes of binding: non-specific binding of α₁-adrenergic receptors is represented by k_d · R_f and binding of free α₁-adrenergic receptors to catecholamines is assumed to follow mass action kinetics with rate constant k_Hr. Dividing Equation 9 by R_total yields the following equation for the balance of free α₁-adrenergic receptor fraction, as shown in Equation 10. Here, r represents the fraction of free (unbound) α₁-adrenergic receptors on the cell surface.

$\begin{array}{l}(d{R}_{(f,i)})/dt=k_{ri}·(R_{total}-R_{(f,i)})\\ -k_d·R_{(f,i)}-k_{Hr}·H_i·R_{(f,i)}& (9)\end{array}$

$\begin{array}{l}(d{r}_i)/dt=k_{(r_i)}·(1-r_i)-k_d·r_i\\ -k_{Hr}·H_i·r_i& (10)\end{array}$

where $k_{(r_i)}\in \left[0.5,1\right]s^{-1}$, $k_d=0.34s^{-1}$, $k_{Hr}=1n{M}^{-1}·s^{-1}$. Parameterization of H_i (nM) is described in further detail in Equations (16) and (17).

We assume saturation kinetics for catecholamine induced IP₃ synthesis in the cytosol as expressed in Equation (11). The positive feedback between cytosolic Ca²⁺ (CaI) and IP₃ is expressed as a multiplicative first order Hill function. Metabolism of IP₃ in the cytosolic domain can occur through two different routes: conversion into IP₄ by a 3-kinase and conversion into IP₂, IP₁, and, eventually, inositol by a phosphatase (Sims and Allbritton, 1998). In the current model, the rate of IP₃ decrease is assumed to be single mass action decay, which is expressed as follows:

$\begin{array}{l}(dIP{3}_i)/dt=(k_{(IP{3}_i)}·H·r_i)/(k_{cat}\\ +r_i)·(Ca{I}_i)/(Ca{I}_i+k_3)-D·IP{3}_i\\ -G_{ij}·(IP{3}_i-IP{3}_j)& (11)\end{array}$

where $k_{(IP{3}_i)}\in \left[0.5,0.6\right]s^{-1}$, $G_{ij}=\left\{0,2.5,5\right\}{s}^{-1}$, k_cat = 0.45, k₃ = 1 μM, D = 0.8 s⁻¹.

The cytosolic Ca²⁺ concentration is modeled such that the rate of increase in cytosolic Ca²⁺ depends on the fraction of bound IP₃R (1 − g). Additionally, there is positive feedback coupling between IP₃ and cytosolic Ca²⁺. The parameter L represents constitutive leakage from the cellular Ca²⁺ stores into the cytosol, while CaT represents the concentration of Ca²⁺ in the cellular stores. A second order Hill function describes the initiation of active transport mechanisms pumping Ca²⁺ back into ER stores at the expense of ATP, when there is an increase in cytosolic Ca²⁺. This process is modeled as given in Equation 12.

$\begin{array}{l}(dCa{I}_i)/dt=(1-g)·((A·IP{3}_i^4)/({(k_1+IP{3}_i)}^4)\\ +L)·(Ca{T}_i-Ca{I}_i)-(B·Ca{I}_i^2)/(k_2^2\\ +Ca{I}_i^2)& (12)\end{array}$

where A = 0.20 μM · s⁻¹, k₁ = 1 μM, L = 0.00015 μM · s⁻¹, B = 0.082 μM · s⁻¹, k₂ = 0.15 μM.

The balance of free (unbound) to total IP₃ receptors is represented by Equation 13, where g is the fraction of free (unbound) IP₃R on the ER surface. The rate of inactivation of IP₃R by cytosolic Ca²⁺, E, is proportional to the fourth power of cytosolic Ca²⁺ concentration. The model also assumes a constitutive reactivation rate of IP₃R, denoted by F.

$\begin{array}{ll}(d{g}_i)/dt=E·Ca{I}_i^4·(1-g_i)-F& (13)\end{array}$

where E = 1 μM⁻⁴ · s⁻¹, F = 0.01 s⁻¹.

Detailed description of parameters and initial conditions for the lobular scale Ca²⁺ portion of the model can be found in Tables 1, 2, respectively. We explicitly consider the zonation profiles of the following parameters: α₁-adrenergic receptors exhibit a decreasing PC to PP gradient (Halpern et al., 2017; Verma et al., 2018). Therefore, the intracellular Ca²⁺ signaling parameters k_r and k_IP3 are initialized as decreasing gradients from the PC to PP region (k_r ∈ [0.5, 0.6]; k_IP3 ∈ [0.5, 1]). The parameter gradients across the 15 hepatic layers can be found in Table 3.

TABLE 1

Table 1. Nominal values/ranges used to initialize hepatocyte intracellular parameters in the main model of lobular scale Ca²⁺ dynamics.

TABLE 2

Table 2. Initial conditions for all hepatocytes in the model of lobular scale Ca²⁺ dynamics.

TABLE 3

Table 3. k_r and k_IP3 parameter gradients along the lobule applicable to the main model.

Modeling CNS Control of Circulating Catecholamines

Circulating catecholamines play an important role in systemic glucose homeostasis by inducing Ca²⁺ spiking in hepatocytes leading to increased glycogenolysis and glucose output from the liver. We assume a CNS-mediated activation of adrenal glands leading to release of epinephrine and norepinephrine into the systemic blood supply. We use the following equations for rates of increase in systemic epinephrine and norepinephrine:

$\begin{array}{ll}dEpn/dt=v_{Epn}(1+X_{Epn}*CNSAct)& (14)\end{array}$

$\begin{array}{ll}dNEpn/dt=v_{NEpn}(1+X_{NEpn}*CNSAct)& (15)\end{array}$

where $v_{Epn}=0.5nM·s^{-1}$; $v_{NEpn}=1.5nM·s^{-1}$; X_Epn = 6; X_NEpn = 14. X_Epn and X_NEpn describe the stimulatory factors due to CNS activation. v_NEpn > v_Epn is in accordance with the observation that norepinephrine levels rise faster and attain higher values during exercise (Christensen et al., 1979; Horton et al., 1998). Catecholamines released in the blood stream then pass through the sinusoidal blood compartments and induce intracellular Ca²⁺ spiking. We assume that catecholamines are cleared only in the liver during their passage through sinusoids following mass action kinetics.

Modeling Catecholamine-Mediated Activation of Signaling in a Liver Lobule

The zonated, short range Ca²⁺ mobilizing effect of norepinephrine was approximated by additional stimulus only to the innervated hepatocytes as shown below:

$\begin{array}{l}{H}_i=Ep{n}_i+NEp{n}_i+K_{NEp{n}_{stim}}·CNSAct;\\ i\in innervatedhepatocytes& (16)\end{array}$

$\begin{array}{ll}{H}_i=Ep{n}_i+NEp{n}_i;i\in non-innervatedhepatocytes& (17)\end{array}$

where Epn_i and NEpn_i represent the concentrations of the catecholamines in the i^th sinusoidal blood compartment and H_i represents the Ca²⁺ mobilizing stimulus to the i^th hepatocyte compartment. The value of K_NEpnstim was set to 2 nM in all simulations with innervation.

Modeling Variations in Blood Flow

Given a normal blood flow rate, hepatic portal venous gas (HVPG), or the blood pressure at the hepatic portal vein, is roughly 5 mmHg. In cases of severe portal hypertension, HVPG may increase up to 3-fold, reaching 15 mm Hg, thereby altering normal blood flow rate (Kumar et al., 2008; Lake-Bakaar, 2015). We used Poiseuille's law as denoted in Equation 18 below to describe the blood pressure differential across a liver sinusoid between PP and PC compartments. The proportional relationship between blood pressure and blood flow rate allows for calculation of blood flow rate during severe portal hypertension.

$\begin{array}{ll}\Delta p=8\mu LQ/(\pi R^2)& (18)\end{array}$

where Δp = blood pressure difference across the lobule, μ = viscosity, L = length of blood vessel, Q = blood flow rate, R = sinusoidal radius. For simplicity, we assume there is no change in viscosity of the blood, sinusoidal length, or diameter due to an increase in portal hypertension. Therefore, increased blood flow during portal hypertension can be simulated by increasing the blood flow parameter (bf = 6.75 s⁻¹ for a 15-compartment lobule).

Model Calibration and Assumptions

Metabolism and Zonation

In the current model we have preserved the dynamics of glucose and glycogen metabolism, zonation of related parameters, and systemic glucagon and insulin dynamics published in the original model. We assume that the increase in glucagon modeled in the original work includes the effect of autonomic activation implicitly.

Intracellular and Lobular Scale Ca²⁺ Signaling

Under a hormone stimulus, hepatocytes typically exhibit inter-spike intervals of 100 s (Robb-Gaspers and Thomas, 1995; Verma et al., 2018). Based on this rough estimate of inter-spike interval for catecholamine induced Ca²⁺ spikes in hepatocytes, nominal parameter values and initial conditions used in simulations are listed in Tables 1, 2, respectively. The spatial zonation of hormone receptor expression and IP₃ synthesis governing calcium dynamics are modeled as graded levels of k_r and k_IP3 increasing from the PP to PC region, with corresponding values listed in Table 3.

It must, however, be noted that cytosolic Ca²⁺ spiking exhibits wide ranges of inter-spike intervals, durations, and latency periods due to the complex intracellular Ca²⁺ signaling machinery and cell-to-cell heterogeneity. Using single values for model parameters that describe cytosolic Ca²⁺ spikes does not capture this heterogeneity exhibited by individual cells. Additionally, change in extracellular stimulus strength alters the inter-spike interval.

Upon hormonal stimulation, Ca²⁺ waves have typically been observed to propagate at rates between 15 and 25 μm/s (Robb-Gaspers and Thomas, 1995; Gaspers and Thomas, 2005). We tuned the model parameters to match this approximate Ca²⁺ propagation speed. Additionally, we considered Ca²⁺ oscillations induced only by catecholamines.

Effect of Catecholamines on Glycogenolysis

We model the effect of catecholamines on hepatic glycogenolysis as mediated explicitly via the potentiation of intracellular Ca²⁺ signaling and a downstream increase in glycogen phosphorylase kinase activity. We calibrated the exercise-induced increase in systemic glucose, epinephrine, and norepinephrine to reflect typical fold changes in human subjects (Horton et al., 1998).

Innervation

We assume different spatial innervation patterns in rodent vs. human livers. In the rodent case, we assume innervation and norepinephrine release at the four hepatocyte layers nearest the PP region. In contrast, in the human case, we assume all hepatocyte layers are innervated.

Simulation Methodology—Effects of Exercise on Glucose Homeostasis in Liver

We analyzed the simulated model output for different Ca²⁺ mobilizing stimuli and gap junction conductivity profiles during exercise. Exercise was modeled as a constantly elevated rate of glucose and FFA uptake by the muscles over a 1 h interval after an initial 500 s interval without exercise. We assumed a 3.5-times elevated rate of glucose uptake during the simulated period of exercise. We did not consider any food intake during the simulated interval. Our model does not consider the hemodynamic changes induced in the liver by CNS activation (Gardemann et al., 1993).

Although connexins are expressed heterogeneously between the PC and PP regions, we assume uniform gap junction conductivity between all adjacent hepatocyte pairs. The strength of gap junction connectivity on the system was investigated by setting G_ij ∈ {0, 2.5, 5}, where G_ij = 0 represents uncoupled gap junctions, G_ij = 2.5 represents weak gap junction connectivity, and G_ij = 5 represents strong gap junction connectivity. Species-specific autonomic control was simulated by altering the number of hepatic layers that were innervated. In simulating human-like extensive innervation, all hepatic layers were innervated and followed the structure of Equation 16. Rodent-like minimal innervation simulations had innervation of the first four periportal hepatocytes, while the following eleven were denervated. Therefore, the first four layers followed the structure of Equation 16, and the following eleven layers followed the structure of Equation 17. Finally, complete denervation of the liver was simulated by setting parameter K_(NEpnStim) = 0 in Equation 16, thereby following the structure of Equation 17 for all hepatic layers.

All simulations were performed using XPPAUT (ver. 8.0, Ermentrout, 2002), which is freely available for download at http://www.math.pitt.edu/~bard/xpp/xpp.html, and an XPP-MATLAB interface (ver. 070626, https://github.com/robclewley/xpp-matlab), available for download from GitHub (https://github.com/Daniel-Baugh-Institute/xpp-matlab). Each simulation ran for ~60 s of wall clock time on a 56-core server with Intel Xeon CPU E5-2697. The model was simulated as a simultaneous system of Ordinary Differential Equations (ODEs).

Model Alternative

In addition to the main model discussed in this manuscript, we simulated a model alternative. The alternative model differs from the main model in the following aspects: (1) Systemic catecholamine synthesis rates are not influenced by CNS feedback, representing a case of denervation of the adrenal medulla, and (2) CNS feedback has an additive effect on glucagon synthesis, representing a case where alternative pathways of glucagon synthesis are activated in the pancreas due to CNS feedback. This yields an alternate set of equations for the glucagon, epinephrine, and norepinephrine balance to replace above Equations (3), (14), and (15) with the below:

$\begin{array}{l}dGlu/dt=1/{\tau }_{glu}·{(\text{ln}(G_{ref}/G_B)}^{(n_{glu})}\\ /[{(k_m^{gln})}^{(n_{glu})}+\text{ln}({(G_{ref}/G_B)}^{(n_{glu})})]\\ +(v_{glu}·CNSAct)/(k_m^{glu}+CNSAct),if{G}_B\\ <G_{ref}& (19)\end{array}$

$\begin{array}{ll}dEpn/dt=v_{Epn}& (20)\end{array}$

$\begin{array}{ll}dNEpn/dt=v_{NEpn}& (21)\end{array}$

where $v_{glu}=0.1nM{s}^{-1}“and{}^{″}{k}_m^{glu}=0.1$. All other alternative model parameters are equivalent to the main model parameters.

In the alternative implementation of the model, Gij values were considered as 1 and 2 for weak and strong gap junction connectivity, respectively. The list of altered parameters is included in Table 4. The intracellular Ca²⁺ signaling parameters k_r and k_IP3 are initialized as decreasing gradients from the PC to PP region (k_r ∈ [1, 2]; k_IP3 ∈ [1, 2]) in the alternate model as these values led to physiologically relevant results during calibration of the model (Table 5). We discuss additional key differences in parameterization between the two models given the case of human-like innervation in Supplementary Figure 1.

TABLE 4

Table 4. Nominal values/ranges used to initialize hepatocyte intracellular parameters in the alternative model are the same as the main model, with the exception of:

TABLE 5

Table 5. k_r and k_IP3 parameter gradients along the lobule applicable to the alternative model.

Model Reproducibility and Credible Practice

All simulation results and figures presented in the current work were repeated independently by a member of our research group not involved in the original modeling and simulation study. The reproduced results were in agreement with the results presented here. Parameter values in the code were cross-checked for consistency against the values provided in Tables 1–3 of the manuscript. For the purpose of assessing reproducibility, all model files were downloaded from GitHub (https://github.com/Daniel-Baugh-Institute/SpatialLiverModel; ver. 2, 2021), and are also available as a supplement to the present manuscript. Self-assessment of the Ten Simple Rules for Credible Practice in Modeling and Simulation in Healthcare was performed and is included in the supplement (Supplementary File 3, Erdemir et al., 2020).

Results

We developed a multi-scale multi-organ model of autonomic control of hepatic glucose metabolism through regulation of systemic, intercellular, and intracellular signaling (Figure 1). Built on previous efforts to account for metabolic zonation and blood flow along porto-central axis within a liver lobule (Ashworth et al., 2016), and intra and intercellular calcium signaling along a liver sinusoid (Verma et al., 2016, 2018), the present model combines these earlier models with a new model of hepatic innervation and regulation of systemic hormonal levels to account for CNS control of organismal scale signaling and hepatic calcium dynamics and glucose metabolism (Figure 1). We simulated the multi-scale model to explore hepatic glucose output and systemic glucose recovery under a range of biological conditions, including increased intensity and duration of exercise, and denervation of the adrenal gland and liver. By accounting explicitly for the extent of neural innervation along a liver sinusoid in the model, we were able to assess the impact of variations in hepatic innervation across species. The simulations and results are discussed in the following sections.

Modeling Metabolic Response to Synaptic and Hormonal Signals in a Liver With Human-Like Extensive Hepatic Innervation

As a qualitative validation of the model, we simulated the effect of no exercise on the model to compare simulated catecholamine and glucagon levels against those during exercise. In agreement with levels reported in Horton et al. (1998), there was roughly a 2-fold increase in systemic epinephrine and norepinephrine concentrations (Figure 2A), while the systemic glucagon concentration (Figure 2B) agreed with the results from the previously developed Ashworth et al. (2016) model. Further validation of the model was conducted by simulating the effect of systemic and locally released hormones on hepatic glycogenolysis. Our simulations showed that hepatic glycogenolysis is primarily driven by glucagon. Systemic catecholamines and norepinephrine released intra-hepatically result in a relatively small increase (13.27% for strong gap junction connectivity) in hepatic glycogenolysis mediated by the initiation of Ca²⁺ response in hepatocytes (Figure 2C).

FIGURE 2

Figure 2. Model calibration. (A) Systemic catecholamine concentrations at t = 4,000 s, when a pseudo steady-state is reached, show roughly a 2-fold increase with exercise, consistent with the values published in Horton et al. (1998). (B) Systemic glucagon concentration at t = 4,000 s, when a pseudo steady state is reached, also show roughly a 2-fold increase with exercise. (C) Hepatic glycogenolysis only increased by 13.27% for strong gap junction connectivity when systemic catecholamines are released intra-hepatically.

We simulated the multi-scale multi-organ model in which all layers of the liver lobule received neural inputs, mimicking that of human hepatic innervation (Figure 3A). We considered three distinct scenarios: (1) no gap junctions between hepatocytes, hence no exchange of IP₃ (G_ij = 0); (2) weak exchange of IP₃ via gap junctions (G_ij = 2.5); and (3) strong connectivity promoting high IP₃ flux between adjacent hepatocyte layers (G_ij = 5). Our simulation results show that the total lobular glucose output was not substantially different across the three scenarios, however, a system with strong gap junction connectivity resulted in the largest glucose output, closely followed by a system with weak gap junction connectivity (Figure 3B). Glucose production within the first and second hepatic layers closest to the PP region were most sensitive to variation in gap junction connectivity (Figure 3C). This can be explained by induction of robust intracellular Ca²⁺ response in PP hepatocytes due to gap junction-mediated IP₃ transfer from adjacent cells. Our simulations further indicated that the periportal zonation of glycogenolysis is preserved in highly innervated livers, as evidenced by the gradient of glucose production shown in Figure 3C. Intrahepatic glucose production was highest in the periportal region, decreasing monotonically toward the pericentral region, irrespective of weak or strong gap junction connectivity (Figure 3C). An increase in gap junction connectivity caused a small increase in the contribution of PP hepatocytes toward total hepatic glucose output (Figure 3D).

FIGURE 3

Figure 3. Model of human-like extensive innervation. (A) Schematic representation of the model with human-like extensive innervation to all 15 hepatic layers. (B) Total hepatic glucose produced increases with strength of gap junction connectivity. (C) Hepatic glucose produced across each layer monotonically decreased from the periportal (PP) to pericentral (PC) region for weak and strong gap junction connectivity. (D) A small increase in PP glucose output contributes toward the higher total glucose output as gap junction connectivity increases. (E) Calcium spiking along the lobule is unsynchronized when there is no gap junction coupling. Weak and strong gap junction connectivity show synchronized calcium wave propagation from the PC to PP region. Calcium wave velocity increases from a system of weak to strong gap junction connectivity.

These results suggest that when all layers of the liver lobule are innervated, gap junction mediated IP₃ exchange between hepatocytes has little to no effect on the overall hepatic glucose response to neural and systemic catecholamine signals. This result is consistent with the existing knowledge that hepatic glycogenolysis is regulated primarily by glucagon. However, the intracellular calcium signals elicited in the three scenarios were significantly different. In the simulation with uncoupled hepatocytes, intracellular Ca²⁺ dynamics are unsynchronized, with the Ca²⁺ spikes for each hepatocyte showing its own intrinsic frequency (Figure 3E, left). Gap junction coupling leads to synchronized waves that start in the PC region and propagate toward the PP region, consistent with Ca²⁺ dynamics reported in recent literature (Figure 3E, middle and right; Verma et al., 2018; Gaspers et al., 2019). While the period of Ca²⁺ waves is very similar between the two coupled systems, stronger gap junction connectivity resulted in a shorter lag time (defined as the time between calcium spike maxima of the extreme PC and PP hepatocytes) and a larger Ca²⁺ wave velocity compared with the system with strong gap junction connectivity. However, the system with weak connectivity still displays increased glucose production relative to the uncoupled case, calcium dynamics, and wave velocities similar to those seen experimentally (Robb-Gaspers and Thomas, 1995; Gaspers and Thomas, 2005). Additionally, the amplitude of Ca²⁺ spikes in hepatocytes in the PP region was lower compared with that in mid-lobular and PC regions for both weak and strongly coupled gap junctions (Supplementary Figure 2).

Modeling Metabolic Response to Synaptic and Hormonal Signals in a Liver With Rodent-Like Minimal Hepatic Innervation

Next, we simulated the model innervating only the four hepatocyte layers residing closest to the portal triad (Figure 4A). This scheme of minimal neural inputs in the liver parenchyma serves as a surrogate for the extent of hepatic innervation generally observed in mice and rats (Reilly et al., 1978; Metz and Forssmann, 1980; Akiyoshi et al., 1998). In this case, our simulation results show lobular scale Ca²⁺ dynamics different from those observed following human-like innervation. Ca²⁺ waves in periportally innervated lobules start in the PP region and propagate toward the PC region despite the simulated cell-intrinsic gradients of intracellular Ca²⁺ signaling parameters favoring a PC to PP flow of Ca²⁺ waves (Figure 4B). A combination of direct neuronal stimulus, gradient of circulating catecholamines along the porto-central axis, and spatial gradients of cell-intrinsic Ca²⁺ signaling parameters were responsible for the shift in wave start site between the two cases.

FIGURE 4

Figure 4. Model of rodent-like minimal innervation. (A) Schematic representation of the model with rodent-like minimal innervation to the first 4 periportal layers. (B) Calcium waves begin in the PC region despite gap junction connectivity. (C) Total hepatic glucose production increases with strength of gap junction connectivity. (D) Hepatic glucose production decreases from the PP to PC region with the exception of uncoupled gap junctions. (E) Despite gap junction connectivity, there is an increase in total hepatic glucose production from the case of rodent-like minimal innervation to human-like extensive innervation.

Similar to the case with human-like innervation, our simulations indicate increased total hepatic glucose output with increasing gap junction connectivity (Figure 4C). Despite the PP to PC Ca²⁺ wave across the lobule, the PP zonation of glycogenolysis is preserved (Figure 4D). There is consistently higher hepatic glucose output in the human-like innervation case vs. the rodent-like innervation case (Figure 4E) irrespective of gap junction connectivity due to more innervation and consequently higher catecholamine release at additional neural synapses. Our simulations show that the hepatic glucose output in the case of rodent-like innervation exhibiting strong gap junction connectivity is similar to that of human-like innervation, indicating the compensatory role played by strong gap junction connectivity in poorly innervated livers (Figure 4E).

Next, we asked the question: can the PC to PP Ca²⁺ wave propagation in lobules exhibiting minimal innervation be restored by modulating systemic catecholamine levels? To simulate this scenario, we increased catecholamine synthesis rates by adrenal glands by 20% (Figures 5A,B), as baseline circulating catecholamine concentrations are higher in rodents, as reported in de Champlain et al. (1976). Our simulations revealed that with high gap junction connectivity, Ca²⁺ waves travel from PC to PP (Figure 5C). In addition, a higher circulating stimulus led to an increase in Ca²⁺ wave propagation speed and higher hepatic glucose output (Figures 5D,E). Thus, our simulations revealed a compensatory mechanism for poorly innervated livers, mediated by stimulation of extra-hepatic processes, that can mimic Ca²⁺ dynamics and increase total glucose output to levels near those observed in highly innervated livers (Figure 5E).

FIGURE 5

Figure 5. Model of rodent-like minimal innervation with slightly amplified catecholamine secretion (1.2X). (A) Schematic representation of model with rodent-like minimal innervation to the first 4 periportal layers. (B) Network diagram of systemic signals in which the adrenal glands increase catecholamine secretion by 20%. (C) Calcium waves now propagate from the PC to PP regions following a 1.2-fold increase in circulating catecholamine secretion. (D) Strong gap junction connectivity in a system with increased catecholamine secretion results in PC to PP calcium waves with higher velocities compared with that of a system with normal catecholamine secretion, in which waves propagate from the PC to PP region with lower velocities. (E) Despite gap junction connectivity, there is an increase in total hepatic glucose output from a system of normal catecholamine secretion to one with a 20% increase in catecholamine secretion.

Modeling Metabolic Response to Hormonal Signals in a Denervated Liver

Next, we explored the overall hepatic glucose response to hormonal signals in the absence of hepatic innervation (Figures 6A,B). We were interested in examining whether gap junction coupling between adjacent hepatocyte layers would provide a mechanism for local tissue scale coordination of calcium signaling in the absence of top-down neural inputs to each hepatocyte layer. We simulated the model to explore the effect of change in gap junction coupling on hepatic glucose output in the case of liver denervation. The denervated state of the liver was represented in our model by excluding the CNS activation dependent additional stimulus received by hepatocytes (K_Stim = 0 in Equation 16).

FIGURE 6

Figure 6. Model of a denervated liver. (A) Schematic representation of a denervated liver. (B) Network diagram showing no direct CNS input to the liver. (C) Hepatic glucose production is reduced in the denervated liver as compared with the innervated liver. (D) Calcium spiking is unsynchronized along the lobule when Gij=0. When gap junction connectivity is weak, calcium waves exhibit “islands” of synchronized Ca²⁺ response. In the strongly coupled system, calcium waves begin in the PP region and propagate toward to PC region. (E) The zonation of hepatic glycogenolysis is consistent with system of innervation: hepatic glucose output decreases from the PP to PC region. (F) Total glucose output is highest when gap junction connectivity is strong.

Our simulations reveal that total glucose output is diminished in the completely denervated liver as compared with the human-like and rodent-like innervated models (Figure 6C). In denervated livers, Ca²⁺ waves start in the PP region and propagate toward the PC region (Figure 6D). In this case, the Ca²⁺ waves follow the gradient of circulating catecholamines across the liver lobule, which is the sole model species regulating Ca²⁺ response in the absence of innervation. An interesting case arises for low gap junction conductivity where a superposition of lobular scale and localized Ca²⁺ responses result in complex signaling patterns (Figure 6D, t = 2,950). We have previously reported “islands” of synchronized Ca²⁺ response within lobular scale Ca²⁺ waves in murine livers (Verma et al., 2018). Our simulations indicate that low gap junction connectivity in combination with extracellular stimulus profiles could drive complex Ca²⁺ signaling patterns across liver lobules. The simulated zonation of glycogenolysis in denervated livers is consistent with that in innervated livers, with PP hepatocytes exhibiting the highest levels of glycogenolysis (Figure 6E).

In denervated livers, strong gap junction connectivity can increase hepatic glucose output (Figure 6G). The increase in hepatic glucose output observed in the simulations can be attributed to lobular scale Ca²⁺ signaling. In our model, PP hepatocytes show the highest rate of glycogenolysis but the lowest Ca²⁺ spiking amplitudes in response to a mobilizing stimulus. Strong gap junction connectivity results in a more robust Ca²⁺ spiking in the PP hepatocytes which further increases glycogenolysis rates.

In the absence of liver innervation, our model simulations show a PP to PC flow of Ca²⁺ signal in the lobule. However, ex vivo experiments have revealed that Ca²⁺ waves propagate from the PC to PP region (Verma et al., 2018). Increasing the base catecholamine synthesis rate by 1.4-fold from the adrenal glands as represented schematically in Figures 7A,B, and restored the experimentally observed directionality of Ca²⁺ signal propagation in our simulations (Figure 7C). In addition, there is an increase in hepatic glucose output compared with simulations of the denervated liver under normal circulating catecholamine levels (Figure 7E). Total glucose output is consistent across simulations of liver denervation with a 1.4-fold increase in rate of catecholamine secretion, rodent-like innervation with a 20% increase in catecholamine secretion, and human-like innervation with baseline rates of catecholamine secretion, all despite varying gap junction connectivity (Figure 7F). Comparing the results in Figures 6, 7 with that of Figures 4, 5, it appears that strong gap junctional connectivity is needed as the hepatic innervation diminishes in order to maintain the calcium spatial patterns elicited by the integrated effect of circulating and synaptically released catecholamines in the liver.

FIGURE 7

Figure 7. Model of a denervated liver with amplified adrenal catecholamine secretion (1.4X). (A) Schematic representation of a denervated liver. (B) Network diagram of systemic signals with the adrenal glands increasing catecholamine secretion by 1.4-fold. (C) Hepatic glucose production decreases from the PP to PC region, consistent with simulations from the denervated liver with normal catecholamine secretion. (D) Calcium waves now propagate from the PC to PP region given weak and strong gap junction connectivity and increased catecholamine secretion. Calcium wave velocities increase with strength of gap junction connectivity. (E) There is an increase in total glucose output from a denervated system with normal catecholamine secretion to one with increased catecholamine section (1.4X). (F) Total glucose output is consistent across scenarios of denervation with a 1.4-fold increase in catecholamine secretion, rodent-like innervation with a 1.2-fold increase in catecholamine secretion, and human-like innervation with normal catecholamine secretion rates.

Modeling the Impact of Portal Hypertension Mediated Increase in Sinusoidal Blood Flow

Next, we investigated the effect of elevated blood flow rates on hepatic glucose output driven by catecholamines during periods of increased systemic glucose demand (Figure 8A). We used Poiseuille's law to calculate the blood flow rate through liver lobules corresponding to a 3-fold increase in portal pressure (see Methods). An increase in blood flow serves as a proxy for simulating hepatic Ca²⁺ and glycogenolytic dynamics in cases with increase in portal hypertension. Consistent with clinical findings (Joly et al., 1967; Gaudin et al., 1991), we simulated equal rates of catecholamine secretion from the adrenal glands for normotensive and portal hypertensive cases and focused our analysis on the interplay of innervation and blood flow to modulate hepatic glycogenolysis in the two cases.

FIGURE 8

Figure 8. Model of an innervated liver with hypertensive-like blood flow. (A) Schematic representation of a liver with hypertensive-like blood flow and innervation tuned to species. (B) Total hepatic glucose production is higher in a system of rodent-like innervation with normotensive conditions compared with hypertensive conditions, despite gap junction connectivity. (C) Total hepatic glucose production is higher in a system of human-like innervation with normotensive conditions compared with hypertensive conditions, despite gap junction connectivity. (D) The hepatic norepinephrine concentration decreases in the PP region and increases in the PC region during hypertension relative to normotensive conditions. (E) The hepatic epinephrine concentration also decreases in the PP region and increases in the PC region during hypertension relative to normotensive conditions. (F) Given minimal rodent-like innervation, PC (PP) hepatic layers show higher (lower) glucose output in the hypertensive case as compared with the normotensive case. (G) In agreement with rodent-like innervation, extensive human-like innervation simulations also show PC (PP) hepatic layers with higher (lower) glucose output in the hypertensive case as compared with the normotensive case. (H) Total glucose output was highest when gap junction connectivity was strong in the case of minimal rodent-like innervation. (I) Total glucose output was highest when gap junction connectivity was weak in the case of extensive human-like innervation.

Simulation results indicate lower hepatic glucose output with increased rate of blood flow, irrespective of the extent of innervation (rodent-like minimal or human-like extensive) and the gap junction connectivity, suggesting that the capacity of the liver to break down intracellular glycogen reserves into glucose is reduced in the case of higher sinusoidal blood flow (Figures 8B,C). The reduction in hepatic glucose output is consistent with reduced catecholamine concentration profiles in the portal compartment (Figures 8D,E). Systemic catecholamines also experience a decrease in concentration during portal hypertensive cases as the glucose synthesized by hepatocytes enter the systemic compartment faster during periods of increased blood flow. This is reflected by a lower drop in systemic glucose during portal hypertension, as displayed in Supplementary Figure 3A. In addition, the simulated hemodynamic changes lead to a redistribution of circulating catecholamines across the lobule resulting in altered glycolysis rates throughout the lobule. Increase in blood flow rates results in a shallower gradient of circulating catecholamines (Supplementary Figures 3B,C). The effect of catecholamine redistribution through the lobule is reflected in the simulated zonal contribution of hepatocytes toward total hepatic glucose output. In rodent-like and human-like innervation schemes, hepatocytes closer to the PC (PP) region showed higher (lower) glucose output in the hypertensive case compared with the normotensive case, as indicated by the PP hepatocytes lying below and PC hepatocytes lying above the diagonal line (Figures 8F,G).

Interestingly, our simulations show that intercellular communication and extent of innervation played an important role in hepatic glucose response to hormonal/neural signals in the case of elevated sinusoidal blood flow. For rodent-like innervation, hepatic glucose output was highest for strong gap junction connectivity in normotensive as well as hypertensive cases (Figure 8H). However, for human-like innervation with portal hypertension, hepatic glucose output is highest when gap junction connectivity is weak (Figure 8I). Our simulations thus show that in addition to zonation and Ca²⁺ signaling, blood flow could add as an additional factor regulating lobular patterns of glycogenolysis.

Due to redistribution of circulating catecholamines discussed previously, Ca²⁺ response is dampened in the portal hypertensive case. Increased intercellular communication (i.e., stronger gap junction connectivity) led to an improved intracellular Ca²⁺ response in the PP hepatocytes (Supplementary Figures 4A,B) and consequently, to a net increase in hepatic glucose output compared with the weakly coupled system given rodent-like minimal innervation (Figure 8H). In a system with human-like extensive innervation, however, weak gap junction connectivity resulted in the highest total glucose output due to a higher overall Ca²⁺ signal resulting in minimal increases in midlobular and PC glucose production relative to the system with strong gap junction connectivity (Figure 8I; Supplementary Figures 4C,D).

Discussion

In this work, we presented a systematic analysis of the regulation of Ca²⁺ signal propagation and hepatic glucose production by zonation, multi-organ interactions, and cell-cell communication in different modes of hepatic innervation. Starting with a previously developed model of hepatic glucose metabolism (Ashworth et al., 2016), we added the effects of CNS activation on the liver directly through innervation and indirectly through increased stimulation of the adrenal glands to release systemic hormones. We expanded the model further to include 15 hepatic layers, and incorporated calcium signaling to capture the effect of gap-junction mediated Ca²⁺ responses on hepatic glycogenolysis. Our simulations show that hepatic glycogenolysis is regulated primarily by glucagon and that catecholamine induced Ca²⁺ dynamics result in only small increases. In all simulated cases, there is a PP to PC zonation of glycogenolysis within liver lobules. However, lobular scale Ca²⁺ dynamics are different despite the initialization of the cell-intrinsic Ca²⁺ signaling parameters being the same across all simulated cases. Variability in spatio-temporal Ca²⁺ dynamics in liver lobules can cause small perturbations in glycogenolysis rates in different hepatocyte layers.

Our simulations demonstrate that the complex interplay of intracellular and multi-organ interactions can override the default modes of Ca²⁺ dynamics and spatial distribution of glycogenolysis in liver lobules. It has been experimentally shown that cell intrinsic Ca²⁺ signaling components are overexpressed in the PC region compared with the PP region. This leads to propagation of Ca²⁺ waves from the PC region to the PP region in gap junction coupled hepatocytes. However, our simulations indicate that depending on the extent of innervation and circulating catecholamines, Ca²⁺ waves can start at different points in lobules showing canonical zonation patterns. Although glycogenolysis is consistently higher in the PP region in all the simulated cases, innervation, gap junction connectivity and circulating catecholamines can influence total hepatic glucose production across liver lobules differentially. Additionally, our simulations show that increased gap junction connectivity resulted in increased hepatic glucose synthesis.

We note that across many simulated cases, the difference in magnitude of total hepatic glucose output is small. One of the key factors likely responsible for the small differences is initialization of cell-intrinsic metabolic and Ca²⁺ signaling parameters. Cell-intrinsic parameters are held constant across all simulations, with the same values used in Ashworth et al. (2016) and Verma et al. (2016), for ease of interpretation of simulation results. However, we expect there to be substantial inter-subject variability, and possibly distinct zonation profiles across liver lobules in these parameters. Nevertheless, our analysis provides qualitative insights into the trends of shifts in intra-hepatic processes that can be expected when modulating key processes. High-resolution proteomic data acquired for model calibration using novel high throughput techniques such as matrix assisted laser desorption ionization (MALDI) will enable more precise quantification of these differences.

Gap junction-mediated synchronization of Ca²⁺ signal propagation has been investigated experimentally in the context of cholestasis, chronic ethanol adaptation, and, more recently, hepatic metabolism (Kruglov et al., 2011; Bartlett et al., 2017; Gaspers et al., 2019, 2021). Our simulations provide novel predictions regarding catecholamine-induced hepatic glycogenolysis in portal hypertension and the role of gap junctions in hypertensive human subjects. Our simulations indicate lower glycogenolysis in highly innervated livers during periods of increased systemic glucose demand under hypertension. Under normal blood flow rates, increase in gap junction coupling leads to increased hepatic glucose output in highly innervated livers. In contrast, hepatic glucose output from highly innervated livers with portal hypertensive blood flow is the highest under weak gap junction coupling. These results suggest that hepatic blood flow is an additional regulatory component that can modulate hepatic glycogenolysis.

Our simulations show that stimulation of adrenal glands by the CNS and the resulting increase of systemic catecholamines can reverse the suppression of hepatic glycogenolysis caused by loss of liver innervation. This is an interesting case in which the site of catecholamine release (adrenal glands) is distinct from the organ of interest (the liver) and is mediated by a third remote entity (CNS), capable of reversing glycogenolytic effects. Novel targeted therapies, rooted in peripheral organ stimulation of the heart (Hanna et al., 2021), gut (Ma et al., 2019), and stomach (Tan et al., 2021) have gained increasing attention recently. The present results using multi-scale modeling of complex interactions between multiple organs across different spatial and temporal scales can support future studies to develop and test neurostimulation therapies. This work and similar efforts could provide a computational proof-of-concept platform for assessing feasibility of such therapies, in addition to providing useful mechanistic insight into the underlying phenomena.

In this research, we combined a model of hepatic energy metabolism (Ashworth et al., 2016) with a model of lobular scale Ca²⁺ signaling (Verma et al., 2016) and added multi-organ interactions to explore the zonal patterns of hepatic glycogenolysis across liver lobules. Similar strategies could be used to explore interesting extensions of our study. For example, the integrated model considers a one-dimensional liver lobule wherein each hepatocyte has either one or two adjacent hepatocytes. Extension of the model to 2-D, 2.5-D, or 3-D cases would be capable of describing a more realistic case with multiple cell-cell interactions between hepatocytes and their effect on overall hepatic glycogenolysis. Availability of specialized applications such as CompuCell3D (Swat et al., 2012) could facilitate building and simulations of spatially accurate computational models. The scope of the model can be further increased by integrating larger models of hepatic physiology such as HepatoNet (Gille et al., 2010). The model can also be extended for computational investigation of disease pathologies such as the interaction of endoplasmic reticulum (ER) stress and hepatic insulin resistance. The mechanisms of obesity associated ER stress and its role in hepatic insulin resistance are well-documented (Kim et al., 2015; Villalobos-Labra et al., 2019). Additionally, ER stress could impact intracellular Ca²⁺ signaling adversely (Krebs et al., 2015). Inclusion of these relevant mechanisms into the current model could provide further insight into the lobular scale Ca²⁺ dynamics and energy metabolism in liver lobules. The extent to which insulin resistance can be reversed by therapeutic strategies targeting intra-hepatic and systemic signals and processes can then be elucidated.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

AV, AM, and RV designed the study. AV and AM developed the model. AM and RN performed simulations. JH, BO, and RV supervised simulations. AV, AM, JH, BO, and RV interpreted the results. AV, AM, RV, and BO edited the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was financially supported by the National Institute of Biomedical Imaging and Bioengineering U01 EB023224 and National Institute on Alcohol Abuse and Alcoholism R01 AA018873. The funding sponsors had no role in the design of the study, in the collection, analyses, or interpretation of data, in the writing of the manuscript and in the decision to publish the results.

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.

Acknowledgments

We thank Ms. Alison Moss for independently reproducing the computational modeling results presented in this study.

Supplementary Material

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

The model file in XPP format can be found in Supplementary File 1 (ModelFile.ode). The alternate model in XPP format can be found in Supplementary File 2 (AlternateModelFile.ode). The MATLAB file used for running the simulations on the XPP-MATLAB interface can be found in Supplementary File 3 (RunOdeFiles.m). The document describing the self-assessed conformance of our model to the state of the art Ten Simple Rules for Credible Practice in Modeling and Simulation in Healthcare can be found in Supplementary File 4.

References

Akiyoshi, H., Gonda, T., and Terada, T. (1998). A comparative histochemical and immunohistochemical study of aminergic, cholinergic and peptidergic innervation in rat, hamster, guinea pig, dog and human livers. Liver 18, 352–359. doi: 10.1111/j.1600-0676.1998.tb00817.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Amaya, M. J., and Nathanson, M. H. (2013). Calcium signaling in the liver. Compr. Physiol. 3, 515–539. doi: 10.1002/cphy.c120013

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashworth, W. B., Davies, N. A., and Bogle, I. D. (2016). A computational model of hepatic energy metabolism: understanding zonated damage and steatosis in NAFLD. PLoS Comput. Biol. 12:e1005105. doi: 10.1371/journal.pcbi.1005105

PubMed Abstract | CrossRef Full Text | Google Scholar

Assimacopoulos-Jeannet, F. D., Blackmore, P. F., and Exton, J. H. (1977). Studies on alpha-adrenergic activation of hepatic glucose output. Studies on role of calcium in alpha-adrenergic activation of phosphorylase. J. Biol. Chem. 252, 2662–2669. doi: 10.1016/S0021-9258(17)40509-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartlett, P. J., Antony, A. N., Agarwal, A., Hilly, M., Prince, V. L., Combettes, L., et al. (2017). Chronic alcohol feeding potentiates hormone-induced calcium signalling in hepatocytes. J. Physiol. 595, 3143–3164. doi: 10.1113/JP273891

PubMed Abstract | CrossRef Full Text | Google Scholar

Baukal, A. J., Balla, T., Hunyady, L., Hausdorff, W., Guillemette, G., and Catt, K. J. (1988). Angiotensin II and guanine nucleotides stimulate formation of inositol 1,4,5-trisphosphate and its metabolites in permeabilized adrenal glomerulosa cells. J. Biol. Chem. 263, 6087–6092. doi: 10.1016/S0021-9258(18)68753-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Berthoud, H.R. (2004). Anatomy and function of sensory hepatic nerves. Anat. Rec. 280A, 827–835. doi: 10.1002/ar.a.20088

PubMed Abstract | CrossRef Full Text | Google Scholar

Billah, M. M., and Michell, R. H. (1979). Phosphatidylinositol metabolism in rat hepatocytes stimulated by glycogenolytic hormones. Effects of angiotensin, vasopressin, adrenaline, ionophore A23187 and calcium-ion deprivation. Biochem. J. 182, 661–668. doi: 10.1042/bj1820661

PubMed Abstract | CrossRef Full Text | Google Scholar

Burcelin, R., Dolci, W., and Thorens, B. (2000). Portal glucose infusion in the mouse induces hypoglycemia: evidence that the hepatoportal glucose sensor stimulates glucose utilization. Diabetes 49, 1635–1642. doi: 10.2337/diabetes.49.10.1635

PubMed Abstract | CrossRef Full Text | Google Scholar

Chao, H.-W., Chao, S.-W., Lin, H., Ku, H.-C., and Cheng, C.-F. (2019). Homeostasis of glucose and lipid in non-alcoholic fatty liver disease. Int. J. Mol. Sci. 20:298. doi: 10.3390/ijms20020298

PubMed Abstract | CrossRef Full Text | Google Scholar

Christensen, N. J., Galbo, H., Hansen, J. F., Hesse, B., Richter, E. A., and Trap-Jensen, J. (1979). Catecholamines and exercise. Diabetes 28, 58–62. doi: 10.2337/diab.28.1.S58

PubMed Abstract | CrossRef Full Text | Google Scholar

Chu, C. A., Sindelar, D. K., Neal, D. W., and Cherrington, A. D. (1996). Direct effects of catecholamines on hepatic glucose production in conscious dog are due to glycogenolysis. Am. J. Physiol. Endocrinol. Metab. 271, E127–E37. doi: 10.1152/ajpendo.1996.271.1.E127

PubMed Abstract | CrossRef Full Text | Google Scholar

de Champlain, J., Farley, L., Cousineau, D., and van Ameringen, M. R. (1976). Circulating catecholamine levels in human and experimental hypertension. Circ. Res. 38, 109–114. doi: 10.1161/01.RES.38.2.109

PubMed Abstract | CrossRef Full Text | Google Scholar

Erdemir, A., Mulugeta, L., Ku, J. P., Drach, A., Horner, M., Morrison, T. M., et al. (2020). Credible practice of modeling and simulation in healthcare: ten rules from a multidisciplinary perspective. J. Transl. Med. 18:369. doi: 10.1186/s12967-020-02540-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Ermentrout, B. (2002). Simulating, analyzing, and animating dynamical systems: a guide to Xppaut for researchers and students. Soc. Ind. Appl. Math. doi: 10.1137/1.9780898718195

CrossRef Full Text | Google Scholar

Exton, J. H. (1987). Mechanisms of hormonal regulation of hepatic glucose metabolism. Diab. Metab. Rev. 3, 163–183. doi: 10.1002/dmr.5610030108

PubMed Abstract | CrossRef Full Text | Google Scholar

Forssmann, W. G., and Ito, S. (1977). Hepatocyte innervation in primates. J. Cell Biol. 74, 299–313. doi: 10.1083/jcb.74.1.299

PubMed Abstract | CrossRef Full Text | Google Scholar

Gardemann, A., Püschel, G. P., and Jungermann, K. (1993). Nervous control of liver metabolism and hemodynamics. EJB Rev. 143–155. doi: 10.1007/978-3-642-78046-2_12

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaspers, L. D., Pierobon, N., and Thomas, A. P. (2019). Intercellular calcium waves integrate hormonal control of glucose output in the intact liver. J. Physiol. 597, 2867–2885. doi: 10.1113/JP277650

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaspers, L. D., and Thomas, A. P. (2005). Calcium signaling in liver. Cell Calcium 38, 329–342. doi: 10.1016/j.ceca.2005.06.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaspers, L. D., Thomas, A. P., Hoek, J. B., and Bartlett, P. J. (2021). Ethanol disrupts hormone-induced calcium signaling in liver. Function 2:zqab002. doi: 10.1093/function/zqab002

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaudin, C., Braillon, A., Poo, J. L., Kleber, G., Moreau, R., and Lebrec, D. (1991). Plasma catecholamines in patients with presinusoidal portal hypertension: comparison with cirrhotic patients and nonportal hypertensive subjects. Hepatology 13, 913–916. doi: 10.1002/hep.1840130518

PubMed Abstract | CrossRef Full Text | Google Scholar

Gille, C., Bölling, C., Hoppe, A., Bulik, S., Hoffmann, S., Hübner, K., et al. (2010). HepatoNet1: a comprehensive metabolic reconstruction of the human hepatocyte for the analysis of liver physiology. Mol. Syst. Biol. 6:411. doi: 10.1038/msb.2010.62

PubMed Abstract | CrossRef Full Text | Google Scholar

Güemes, M., Rahman, S. A., and Hussain, K. (2015). What is a normal blood glucose? Arch. Dis. Childh. 101, 569–574. doi: 10.1136/archdischild-2015-308336

PubMed Abstract | CrossRef Full Text | Google Scholar

Halpern, K. B., Shenhav, R., Matcovitch-Natan, O., Tóth, 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

Han, H.-S., Kang, G., Kim, J. S., Choi, B. H., and Koo, S.-H. (2016). Regulation of glucose metabolism from a liver-centric perspective. Exp. Mol. Med. 48:e218. doi: 10.1038/emm.2015.122

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanna, P., Dacey, M. J., Brennan, J., Moss, A., Robbins, S., Achanta, S., et al. (2021). Innervation and neuronal control of the mammalian sinoatrial node a comprehensive atlas. Circ. Res. 128, 1279–1296. doi: 10.1161/CIRCRESAHA.120.318458

PubMed Abstract | CrossRef Full Text | Google Scholar

Horton, T. J., Pagliassotti, M. J., Hobbs, K., and Hill, J. O. (1998). Fuel metabolism in men and women during and after long-duration exercise. J. Appl. Physiol. 85, 1823–1832. doi: 10.1152/jappl.1998.85.5.1823

PubMed Abstract | CrossRef Full Text | Google Scholar

Jensen, K. J., Alpini, G., and Glaser, S. (2013). hepatic nervous system and neurobiology of the liver. Compreh. Physiol. 655–665. doi: 10.1002/cphy.c120018

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, G., and Zhang, B. B. (2003). Glucagon and regulation of glucose metabolism. Am. J. Physiol. Endocrinol. Metab. 284, E671–E678. doi: 10.1152/ajpendo.00492.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Joly, J. G., Leduc, J., Bernier, J., Lavoie, P., and Viallet, A. (1967). Catecholamine levels in portal, hepatic, and systemic venous blood in portal hypertension. Lancet 2, 121–123. doi: 10.1016/S0140-6736(67)92962-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Jungermann, K. (1983). Functional significance of hepatocyte heterogeneity for glycolysis and gluconeogenesis. Pharmacol. Biochem. Behav. 18, 409–414. doi: 10.1016/0091-3057(83)90208-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Jungermann, K. (1988). Metabolic Zonation of liver parenchyma. Semin. Liver Dis. 8, 329–341. doi: 10.1055/s-2008-1040554

PubMed Abstract | CrossRef Full Text | Google Scholar

Jungermann, K., Gardemann, A., Beuers, U., Ballé, C., Sannemann, J., Beckh, K., et al. (1987). Regulation of liver metabolism by the hepatic nerves. Adv. Enzyme Regul. 26, 63–88. doi: 10.1016/0065-2571(87)90006-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Jungermann, K., and Keitzmann, T. (1996). Zonation of parenchymal and nonparenchymal metabolism in liver. Annu. Rev. Nutr. 16, 179–203. doi: 10.1146/annurev.nu.16.070196.001143

PubMed Abstract | CrossRef Full Text | Google Scholar

Jungermann, K., and Kietzmann, T. (1997). Role of oxygen in the zonation of carbohydrate metabolism and gene expression in liver. Kidney Int. 51, 402–412. doi: 10.1038/ki.1997.53

PubMed Abstract | CrossRef Full Text | Google Scholar

Jungermann, K., and Kietzmann, T. (2000). Oxygen: modulator of metabolic zonation and disease of the liver. Hepatology 31, 255–260. doi: 10.1002/hep.510310201

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, O. K., Jun, W., and Lee, J. (2015). Mechanism of ER stress and inflammation for hepatic insulin resistance in obesity. Ann. Nutr. Metab. 67, 218–227. doi: 10.1159/000440905

PubMed Abstract | CrossRef Full Text | Google Scholar

Krebs, J., Agellon, L. B., and Michalak, M. (2015). Ca(2+) homeostasis and endoplasmic reticulum (ER) stress: an integrated view of calcium signaling. Biochem. Biophys. Res. Commun. 460, 114–121. doi: 10.1016/j.bbrc.2015.02.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Kruglov, E. A., Gautam, S., Guerra, M. T., and Nathanson, M. H. (2011). Type 2 inositol 1,4,5-trisphosphate receptor modulates bile salt export pump activity in rat hepatocytes. Hepatology 54, 1790–1799. doi: 10.1002/hep.24548

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, A., Sharma, P., and Sarin, S. K. (2008). Hepatic venous pressure gradient measurement: time to learn! Indian J. Gastroenterol. 27, 74–80.

PubMed Abstract | Google Scholar

Lake-Bakaar, G. (2015). Hagen-Poiseuille's law: the link between cirrhosis, liver stiffness, portal hypertension and hepatic decompensation. World J. Hepatol. 7:28. doi: 10.4254/wjh.v7.i1.28

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, J. A., Ahmed, Q., Hines, J. E., and Burt, A. D. (1992). Disappearance of hepatic parenchymal nerves in human liver cirrhosis. Gut 33, 87–91. doi: 10.1136/gut.33.1.87

PubMed Abstract | CrossRef Full Text | Google Scholar

Livingstone, C., Lyall, H., and Gould, G. W. (1995). Hypothalamic GLUT 4 expression: a glucose- and insulin-sensing mechanism? Mol. Cell. Endocrinol. 107, 67–70. doi: 10.1016/0303-7207(94)03423-Q

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, Q., Xing, C., Long, W., Wang, H. Y., Liu, Q., and Wang, R. F. (2019). Impact of microbiota on central nervous system and neurological diseases: the gut-brain axis. J. Neuroinflamm. 16:53. doi: 10.1186/s12974-019-1434-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Marty, N., Dallaporta, M., and Thorens, B. (2007). Brain Glucose Sensing, Counterregulation, and Energy Homeostasis. Physiology. 22, 241–251. doi: 10.1152/physiol.00010.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

McCuskey, R. S. (2004). Anatomy of efferent hepatic nerves. Anat. Rec. 280A, 821–826. doi: 10.1002/ar.a.20087

PubMed Abstract | CrossRef Full Text | Google Scholar

Metz, W., and Forssmann, W. G. (1980). Innervation of the liver in guinea pig and rat. Anat. Embryol. 160, 239–252. doi: 10.1007/BF00305105

PubMed Abstract | CrossRef Full Text | Google Scholar

Meyer, T., and Stryer, L. (1991). Calcium spiking. Annu. Rev. Biophys. Biophys. Chem. 20, 153–174. doi: 10.1146/annurev.bb.20.060191.001101

PubMed Abstract | CrossRef Full Text | Google Scholar

Mikines, K. J., Sonne, B., Richter, E. A., Christensen, N. J., and Galbo, H. (1985). Glucose turnover during insulin-induced hypoglycemia in liver-denervated rats. Am. J. Physiol. 248, E327–E332. doi: 10.1152/ajpendo.1985.248.3.E327

PubMed Abstract | CrossRef Full Text | Google Scholar

Mizuno, K., and Ueno, Y. (2016). Autonomic nervous system and the liver. Hepatol. Res. 47, 160–165. doi: 10.1111/hepr.12760

PubMed Abstract | CrossRef Full Text | Google Scholar

Nordlie, R. C., Foster, J. D., and Lange, A. J. (1999). regulation of glucose production by the liver. Annu. Rev. Nutr. 19, 379–406. doi: 10.1146/annurev.nutr.19.1.379

PubMed Abstract | CrossRef Full Text | Google Scholar

Oosterveer, M. H., and Kristina, S. (2014). Hepatic glucose sensing and integrative pathways in the liver. Cell. Mol. Life Sci. 71, 1453–1467. doi: 10.1007/s00018-013-1505-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Orellana, S., Solski, P. A., and Brown, J. H. (1987). Guanosine 5'-O-(thiotriphosphate)-dependent inositol trisphosphate formation in membranes is inhibited by phorbol ester and protein kinase C. J. Biol. Chem. 262, 1638–1643. doi: 10.1016/S0021-9258(19)75684-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Parekh, S., and Anania, F. A. (2007). Abnormal lipid and glucose metabolism in obesity: implications for nonalcoholic fatty liver disease. Gastroenterology. 132, 2191–2207. doi: 10.1053/j.gastro.2007.03.055

PubMed Abstract | CrossRef Full Text | Google Scholar

Radziuk, J., and Pye, S. (2001). Hepatic glucose uptake, gluconeogenesis and the regulation of glycogen synthesis. Diab. Metab. Res. Rev. 17, 250–272. doi: 10.1002/dmrr.217

PubMed Abstract | CrossRef Full Text | Google Scholar

Reilly, F. D., McCuskey, P. A., and McCuskey, R. S. (1978). Intrahepatic distribution of nerves in the rat. Anat. Rec. 191, 55–67. doi: 10.1002/ar.1091910106

PubMed Abstract | CrossRef Full Text | Google Scholar

Robb-Gaspers, L. D., and Thomas, A. P. (1995). Coordination of Ca2+ signaling by intercellular propagation of Ca2+ waves in the intact liver. J. Biol. Chem. 270, 8102–8107. doi: 10.1074/jbc.270.14.8102

PubMed Abstract | CrossRef Full Text | Google Scholar

Schöfl, C., Brabant, G., Hesch, R. D., von zur Mühlen, A., Cobbold, P. H., and Cuthbertson, K. S. (1993). Temporal patterns of alpha 1-receptor stimulation regulate amplitude and frequency of calcium transients. Am. J. Physiol. 265, C1030–C1036. doi: 10.1152/ajpcell.1993.265.4.C1030

PubMed Abstract | CrossRef Full Text | Google Scholar

Sherwin, R. S., and Sacca, L. (1984). Effect of epinephrine on glucose metabolism in humans: contribution of the liver. Am. J. Physiol. Endocrinol. Metab. 247, E157–E165. doi: 10.1152/ajpendo.1984.247.2.E157

PubMed Abstract | CrossRef Full Text | Google Scholar

Sims, C. E., and Allbritton, N. L. (1998). Metabolism of inositol 1,4,5-trisphosphate and inositol 1,3,4,5-tetrakisphosphate by the oocytes of Xenopus laevis. J. Biol. Chem. 273, 4052–4058. doi: 10.1074/jbc.273.7.4052

PubMed Abstract | CrossRef Full Text | Google Scholar

Stümpel, F., Kucera, T., Bazotte, R., and Püschel, G. P. (1996). Loss of regulation by sympathetic hepatic nerves of liver metabolism and haemodynamics in chronically streptozotocin-diabetic rats. Diabetologia 39, 161–165. doi: 10.1007/BF00403958

PubMed Abstract | CrossRef Full Text | Google Scholar

Swat, M. H., Thomas, G. L., Belmonte, J. M., Shirinifard, A., Hmeljak, D., and Glazier, J. A. (2012). Multi-scale modeling of tissues using CompuCell3D. Methods Cell Biol. 110, 325–366. doi: 10.1016/B978-0-12-388403-9.00013-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Tan, Z. T., Ward, M., Phillips, R. J., Zhang, X., Jaffey, D. M., Chesney, L., et al. (2021). Stomach region stimulated determines effects on duodenal motility in rats. Am. J. Physiol. Regul. Integr. Comp. Physiol. 320, R331–R341. doi: 10.1152/ajpregu.00111.2020

PubMed Abstract | CrossRef Full Text | Google Scholar

Verma, A., Antony, A. N., Ogunnaike, B. A., Hoek, J. B., and Vadigepalli, R. (2018). Causality analysis and cell network modeling of spatial calcium signaling patterns in liver lobules. Front. Physiol. 9:1377. doi: 10.3389/fphys.2018.01377

PubMed Abstract | CrossRef Full Text | Google Scholar

Verma, A., Makadia, H., Hoek, J. B., Ogunnaike, B. A., and Vadigepalli, R. (2016). Computational modeling of spatiotemporal Ca2+signal propagation along hepatocyte cords. IEEE Trans. Biomed. Eng. 63, 2047–2055. doi: 10.1109/TBME.2016.2550045

PubMed Abstract | CrossRef Full Text | Google Scholar

Villalobos-Labra, R., Subiabre, M., Toledo, F., Pardo, F., and Sobrevia, L. (2019). Endoplasmic reticulum stress and development of insulin resistance in adipose, skeletal, liver, and foetoplacental tissue in diabesity. Mol. Aspects Med. 66, 49–61. doi: 10.1016/j.mam.2018.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Wahren, J., and Ekberg, K. (2007). Splanchnic regulation of glucose production. Annu. Rev. Nutr. 27, 329–345. doi: 10.1146/annurev.nutr.27.061406.093806

PubMed Abstract | CrossRef Full Text | Google Scholar

Wolf, B. A., Florholmen, J., Turk, J., and McDaniel, M. L. (1988). Studies of the Ca2+ requirements for glucose- and carbachol-induced augmentation of inositol trisphosphate and inositol tetrakisphosphate accumulation in digitonin-permeabilized islets. Evidence for a glucose recognition site in insulin secretion. J. Biol. Chem. 263, 3565–3575. doi: 10.1016/S0021-9258(18)68962-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, H., and Lautt, W. W. (1996). Insulin resistance caused by hepatic cholinergic interruption and reversed by acetylcholine administration. Am. J. Physiol. 271, E587–E592. doi: 10.1152/ajpendo.1996.271.3.E587

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamaguchi, N., and Garceau, D. (1980). Correlations between hemodynamic parameters of the liver and norepinephrine release upon hepatic nerve stimulation in the dog. Can. J. Physiol. Pharmacol. 58, 1347–1355. doi: 10.1139/y80-204

PubMed Abstract | CrossRef Full Text | Google Scholar

Yi, C.-X., la Fleur, S. E., Fliers, E., and Kalsbeek, A. (2010). The role of the autonomic nervous liver innervation in the control of energy metabolism. Biochim. Biophys. Acta 1802, 416–431. doi: 10.1016/j.bbadis.2010.01.006

PubMed Abstract | CrossRef Full Text | Google Scholar