Causality Analysis and Cell Network Modeling of Spatial Calcium Signaling Patterns in Liver Lobules

Dynamics as well as localization of Ca²⁺ transients plays a vital role in liver function under homeostatic conditions, repair, and disease. In response to circulating hormonal stimuli, hepatocytes exhibit intracellular Ca²⁺ responses that propagate through liver lobules in a wave-like fashion. Although intracellular processes that control cell autonomous Ca²⁺ spiking behavior have been studied extensively, the intra- and inter-cellular signaling factors that regulate lobular scale spatial patterns and wave-like propagation of Ca²⁺ remain to be determined. To address this need, we acquired images of cytosolic Ca²⁺ transients in 1300 hepatocytes situated across several mouse liver lobules over a period of 1600 s. We analyzed this time series data using correlation network analysis, causal network analysis, and computational modeling, to characterize the spatial distribution of heterogeneity in intracellular Ca²⁺ signaling components as well as intercellular interactions that control lobular scale Ca²⁺ waves. Our causal network analysis revealed that hepatocytes are causally linked to multiple other co-localized hepatocytes, but these influences are not necessarily aligned uni-directionally along the sinusoids. Our computational model-based analysis showed that spatial gradients of intracellular Ca²⁺ signaling components as well as intercellular molecular exchange are required for lobular scale propagation of Ca²⁺ waves. Additionally, our analysis suggested that causal influences of hepatocytes on Ca²⁺ responses of multiple neighbors lead to robustness of Ca²⁺ wave propagation through liver lobules.

Introduction

The liver performs a wide variety of physiological functions, including the regulation of intermediary metabolism, lipid synthesis, bile production, and xenobiotic detoxification. Normal liver function requires both tight regulation of intracellular processes and intercellular coordination. Free Ca²⁺ in the intracellular domain participates in the regulation of such hepatocyte functions as glucose metabolism, bile secretion, proliferation, and apoptosis (Exton, 1987; McConkey and Orrenius, 1997; Canaff et al., 2001). Regulation of cytosolic Ca²⁺ is particularly important in hepatocytes, the cells responsible for the bulk of metabolic and detoxification activities in the liver. Consequently, disruption of Ca²⁺ dynamics can potentially lead to pathological conditions, such as cholestasis (Kruglov et al., 2011).

Structurally, cells in the liver are arranged in lobules, the functional units of the liver conceptualized as having roughly hexagonal cross sections delineated by the portal triad and the central veins. Circulating blood enters lobules through the portal vein and hepatic artery residing in the periportal region, and is drained into the central vein after passing through sinusoids. Hepatocytes are polarized cells arranged alongside sinusoids, with their basolateral membranes in contact with the systemic blood flow. Upon contact with Ca²⁺ mobilizing agents in the blood stream, such as ATP, hormones, or growth factors, spikes in cytosolic Ca²⁺ concentration are observed within the intracellular domains of hepatocytes (Woods et al., 1986, 1987; Serradeil-Le Gal et al., 1991). Binding of extracellular stimuli, such as hormones to receptors on the basolateral hepatocyte membranes, elicits an intracellular signaling cascade involving phospholipase C (PLC) activation, inositol triphosphate (IP₃) synthesis, IP₃ receptor (IP₃R) activation in the endoplasmic reticulum (ER) membrane, leading to a rapid efflux of Ca²⁺ from the ER into the cytosol. Once in the cytosol, Ca²⁺ can be sequestered by mitochondria, released into the extracellular region, or pumped back into the ER, thus reducing the cytosolic Ca²⁺ levels, yielding a Ca²⁺ spike. Intracellular Ca²⁺ spiking has been reported to arise primarily due to fast activation and slow inhibition of IP₃R by cytosolic Ca²⁺ operating in conjunction with active pumping of cytosolic Ca²⁺ into the ER by SERCA pumps (Atri et al., 1993; Keizer and De Young, 1993). In the intact rat liver, sustained hormone stimulation typically leads to Ca²⁺ spike trains in hepatocytes with inter-spike intervals dependent upon stimulus strength as well as intracellular signaling capacity (Robb-Gaspers and Thomas, 1995).

Heterogeneity in the expression and intracellular distribution of the Ca²⁺ signaling components as well as variability in the extracellular regulatory factors can lead to differences in characteristic Ca²⁺ spiking frequencies of individual hepatocytes. Variations in Ca²⁺ spike frequencies have been shown to lead to differential downstream gene expression and protein regulation (Dolmetsch et al., 1998; Zhu et al., 2008; Smedler and Uhlén, 2014). A coherent lobular scale response to extracellular stimuli requires that the Ca²⁺ signals in hepatocytes across the lobule be coordinated. Gap junctions are hypothesized to play a role in coordinating this response, leading to synchronization of cytosolic Ca²⁺ spikes across the liver lobule in response to G-protein coupled receptor agonists (Robb-Gaspers and Thomas, 1995; Tordjmann et al., 1997). Cytosolic Ca²⁺ and IP₃ from a hepatocyte can migrate to neighboring hepatocytes, likely responsible for inducing synchronization of Ca²⁺ spikes across the liver lobule (Sáez et al., 1989). Another mechanism of long-range coordination may involve the release of paracrine signals, such as ATP, into the extracellular milieu, which then elicits Ca²⁺ spiking in the neighboring hepatocytes by purinergic receptor activation (Schlosser et al., 1996). Exchange of molecular signals between cells through gap junctions or paracrine signaling in addition to spatially organized heterogeneity could lead to a coordinated response within a population of cells with regard to downstream processes regulated by Ca²⁺ in response to extracellular stimuli.

In liver lobules, Ca²⁺ signals commonly manifest as traveling waves (Keizer and De Young, 1993; Robb-Gaspers and Thomas, 1995; Thomas et al., 1996). Ca²⁺ waves usually start in cells located in the pericentral (PC) region of the lobule and propagate toward the periportal (PP) region (Nathanson et al., 1995; Robb-Gaspers and Thomas, 1995) upon lobule-wide stimulation by vasopressin or phenylephrine. The direction of Ca²⁺ signal propagation is opposite to the general direction of blood flow, which is from PP to PC. This observation indicates an organized spatial heterogeneity, termed liver zonation, in the Ca²⁺ signaling capacity of cells. Liver zonation has been observed in many other physiological functions in liver lobules (Gebhardt and Mecke, 1983; Jungermann, 1987; Braeuning et al., 2006; Gebhardt and Matz-Soja, 2014).

Dynamics as well as localization of Ca²⁺ transients plays a vital role in liver function under homeostatic conditions, repair, and disease (Rooney et al., 1990; Pusl and Nathanson, 2004; Gaspers and Thomas, 2008; Lagoudakis et al., 2010; Fu et al., 2012; Amaya and Nathanson, 2013; Bartlett et al., 2014; Oliveira et al., 2015). Although a wealth of information exists regarding the intracellular Ca²⁺ dynamics, understanding of lobular scale propagation of Ca²⁺ signal, its quantification as well as a clear understanding of its significance is lacking. Use of computational modeling to decipher processes that control intracellular Ca²⁺ spikes and spatial patterns of Ca²⁺ signal propagation has been a long-standing area of investigation due to the complexity of their origin and propagation. For instance, Schuster et al. (2002) present a detailed discussion of computational models developed for describing Ca²⁺ spiking, as well as propagation of Ca²⁺ waves through co-localized cells. In previous work, we used a computational model to predict that zonation of intracellular signaling components as well as gap junction-mediated IP₃ exchange between immediate neighbors are required for the propagation of a Ca²⁺ signal through a chain of connected hepatocytes (Verma et al., 2016). A recent study employing single cell RNA sequencing provided evidence in support of our predictions of lobular gradients of intracellular signaling components at the mRNA level (Halpern et al., 2017).

In this study, we combined analysis of experimentally acquired images of cytosolic Ca²⁺ dynamics in mouse liver lobules with dynamic modeling to identify spatial features of lobular scale Ca²⁺ signal propagation and putative causal linkages between adjacent hepatocytes. We imaged cytosolic Ca²⁺ levels in response to a vasopressin stimulus in a 2D optical slice of a perfused intact mouse liver to obtain a data set on Ca²⁺ transients in 1300 hepatocytes residing in different lobules, measured every 4 s over a period of 1600 s. We analyzed correlation as well as causal networks constructed using the acquired high-dimensional time series to characterize the spatial extent and directional alignment of intercellular interactions that lead to Ca²⁺ waves across liver lobules. We incorporated the causal connectivity from network analysis into an ordinary differential equation-based dynamic model of intra-and inter-cellular Ca²⁺ signaling. We utilized the model to evaluate the effect of spatial heterogeneity in the intra- and inter-cellular signaling components on spatial patterns of cytosolic Ca²⁺ signals. Our dynamic model-based analysis predicted the spatial distribution of signaling components that yield lobular scale Ca²⁺ patterns that are consistent with the experimentally observed Ca²⁺ wave propagation.

Methods

Calcium Imaging in Isolated Perfused Mouse Livers

All animal procedures used in this study were handled in accordance with mandated standards of humane care and were approved by the Thomas Jefferson University Institutional Animal Care and Use Committee. Confocal imaging of intact perfused livers was performed as previously described (Robb-Gaspers and Thomas, 1995; Bartlett et al., 2017). Briefly, livers from 8–12 weeks old male C57 BL/6 mice were perfused via the hepatic portal vein with a HEPES-buffered balance salt solution (121 mM NaCl, 25 mM HEPES, 5 mM NaHCO3, 4.7 mM KCl, 1.2 mM KH₂PO₄, 1.2 mM MgSO₄, 1.3 mM CaCl₂, 5.5 mM glucose, 0.5 mM glutamine, 3 mM lactate, 0.3 mM pyruvate, 0.2 mM bromosulfophthalein (BSP), 0.1% BSA, pH 7.4) equilibrated with 100% O₂ at 30°C. A Ca²⁺-sensitive indicator, fluo-8 AM (5 μM) was loaded into the hepatocytes in vivo by recirculating the perfusion buffer supplemented with fluo-8 AM plus 0.02% Pluronic F-127 and 2% BSA for 40–50 min. Confocal images were acquired with an EC Plan-Neofluar 10x/0.30 M27 objective using a Zeiss LSM510MP confocal microscope. Fluo-8 images (488 nm excitation, 520–600 nm emission) were captured every 4 s. Periportal and pericentral zones were identified by differential dye loading and perfusion of fluorescein-conjugated BSA.

Image Segmentation and Cytosolic Calcium Time Trace Extraction

Hepatocytes in the acquired images were segmented manually. Intensities of all pixels lying within segmented hepatocyte boundaries were added for every time slice to obtain a 400-time point cytosolic Ca²⁺ time series for all the segmented hepatocytes.

Pre-Processing Cytosolic Ca²⁺ Time Series Data

The following operations were performed on the cytosolic Ca²⁺ trace of every hepatocyte (see Supplementary Figure S1 for details):

Baseline Correction and Rescaling

The cytosolic Ca²⁺ series for each hepatocyte was detrended using an implementation of the rolling ball baseline correction algorithm contained in the baseline package (version 1.2) in R platform for statistical analysis (version 3.2.3; R Core Team, 2015) to remove low frequency components and correct for dye photobleaching during the experiment.

Low Pass Filtering and Rescaling

High frequency components in each baseline-corrected Ca²⁺ time series were removed using the smooth.fft function from the itsmr package (version 1.5) in R (version 3.2.3). This function removes frequencies corresponding to the highest nth percentile from the power spectrum of a given time series. The signal in the lowest 27.5 percentile of the frequency range in every time series was retained for subsequent analysis. The amplitudes of cytosolic Ca²⁺-dependent fluorescence signal intensity in the resultant time series were then rescaled to between 0 and 1.

Network Analysis

Undirected Correlation Network Construction and Analysis

Undirected correlation networks were constructed using pairwise Spearman rank correlation coefficient estimates of baseline corrected, low pass filtered, and rescaled time series data. Edges corresponding to correlation values < 0.75 or those that were between hepatocytes lying at a distance > 100 μm were discarded. The resultant networks were analyzed for isolated clusters, their sizes and node degrees. Analysis of the correlation network was performed using the igraph (version 1.0.1) package in R (version 3.2.3).

Transfer Entropy (TE) Based Causal Network Construction and Analysis

Transfer entropy (TE) is a measure of the directed influence between two random processes. TE from a process X to another process Y is defined as the amount of reduction in uncertainty of future values of Y by knowing the past values of X, given past values of Y. In the present study, pair-wise TE between Ca²⁺ responses of hepatocytes was estimated based on Shannon's conditional entropy, as follows:

$\begin{array}{rcll}{h}_1& =& -{\sum }_{y_{t+1},Y_t^n,X_t^m}p(y_{t+1},Y_t^n,X_t^m).logp(y_{t+1}|Y_t^n,X_t^m)& (1)\end{array}$

$\begin{array}{rcll}{h}_2& =& -{\sum }_{y_{t+1},Y_t^n,X_t^m}p(y_{t+1},Y_t^n,X_t^m).logp(y_{t+1}|Y_t^n)& (2)\end{array}$

$\begin{array}{rcll}T{E}_{X\to Y}& =& h_2-h_1& (3)\end{array}$

where x_t and y_t represent the cytosolic Ca²⁺ levels in hepatocytes X and Y, respectively, at time t; $X_t^m=\left[x_t,x_{t-1},\dots ,x_{t-m+1}\right]$ and $Y_t^n=\left[y_t,y_{t-1},\dots ,y_{t-n+1}\right]$ are past m and n values of cytosolic Ca²⁺ in respective hepatocytes X and Y; $p(y_{t+1}|Y_t^n,X_t^m)$ is the probability of the occurrence of y_t+1 given $X_t^m$ and $Y_t^n$, and $p(y_{t+1},Y_t^n,X_t^m)$ is the joint probability of occurrence of y_t+1, $X_t^m$, and $Y_t^n$. TE therefore represents the decrease in Shannon's entropy when past values of X and Y are used to predict the current value of Y compared to past values of Y alone. Information transfer is considered as occurring from X to Y if TE_X→Y> 0 (see Schreiber, 2000 for more details). In this work, m and n were taken to be 1 based on cross correlation measures (Figure S2). Additionally, a theoretical estimate of IP₃ diffusion time between two hepatocytes with diameters of ~ 25 μm is 1.1 s based on IP₃ diffusion constant values in xenopus oocyte cytosolic extracts (Allbritton et al., 1992). We therefore limited our TE analysis to a history value of 1 (4 s lag), with information flow interpreted as IP₃ exchange between neighboring hepatocytes. It must be noted that the diffusion time of IP₃ between hepatocytes in vivo might be increased due to molecular charge, gap junction channel properties and cell-intrinsic buffering. However, to our knowledge, no data exists regarding effective diffusivity of IP₃ in mouse hepatocytes in vivo.

Directed causal networks between hepatocytes based on TE were constructed using quantized Ca²⁺ time series of all hepatocytes. Pre-processed Ca²⁺ time series for each hepatocyte was quantized into high (= 1) or low (= 0) cytosolic Ca²⁺ at each time point using the 85th percentile of the cell-intrinsic intensity distribution as considered the threshold. As an additional filter to minimize the effect of noise, all high cytosolic Ca²⁺ values in the time series were changed to zero unless the value at an immediately preceding or following time point was also high, i.e., cytosolic Ca²⁺ intensity was sustained above the 85th percentile within the cell for at least 8 s.

In the absence of a good value of TE to infer cell-to-cell influence, hepatocyte-specific significance testing was employed to identify influence edges and construct TE-based causal networks. For every hepatocyte, pairwise TE values from all other hepatocytes were estimated to obtain an empirical distribution. If the TE value to the given hepatocyte from another adjacent hepatocyte was greater than the 95th percentile of its empirical TE distribution, a positive causal influence was considered from the neighbor to the hepatocyte of interest. The similarity in TE networks identified by our method based on binarized data and a continuous TE estimation method implemented in JIDT (Lizier, 2014) can be found in Figure S11. Additionally, we chose a cell-specific TE threshold instead of a global threshold to avoid inclusion of false positives (Figure S12).

Computational Modeling of Intra- and Inter-Cellular Ca²⁺ Signaling

We started with a receptor oriented, ordinary differential equation (ODE)-based model of Ca²⁺ signal propagation in a cord of hepatocytes detailed in Verma et al. (2016). Here, we consider the complex spatial features of a liver lobule by allowing the hepatocytes to be connected with more than two other hepatocytes, as was the case in the original model of Verma et al. (2016). In the present computational model, the state of every cell “i” and its interaction with a set of adjacent cells represented by the index “j” is defined by the following system of ODEs:

$\begin{array}{rcll}\frac{d{r}_i}{dt}& =& k_{r_i}(1-r_i)-k_d{r}_i-k_{Hr}H.r_i& (4)\end{array}$

$\begin{array}{rcl}\frac{dIP{3}_i}{dt}& =& \left(\frac{k_{IP{3}_i}.H.r_i}{k_{cat}+r_i}\right)\left(1-\frac{k_3}{Ca{I}_i+k_3}\right)-D\frac{IP{3}_i}{2}\end{array}$

$\begin{array}{rcl}-& {\displaystyle \sum _{j\in ad{j}_i}}{G}_{ij}(IP{3}_i-IP{3}_j)& (5)\end{array}$

$\begin{array}{rcl}\frac{dCa{I}_i}{dt}& =& (1-g_i)\left(\frac{A{(\frac{IP{3}_i}{2})}^4}{{(k_1+\frac{IP{3}_i}{2})}^4}+L\right)(Ca{T}_i-Ca{I}_i)\end{array}$

$\begin{array}{rcl}-& \frac{B.Ca{I}_i^2}{k_2^2+Ca{I}_i^2}& (6)\end{array}$

$\begin{array}{rcll}\frac{d{g}_i}{dt}& =& E.Ca{I}_i^4(1-g_i)-F& (7)\end{array}$

$\begin{array}{rcll}\frac{dCa{T}_i}{dt}& =& -\frac{dCa{I}_i}{dt}& (8)\end{array}$

Surface receptor activity (r) including non-specific binding was modeled as shown in Equation (4), where H represents the hormone level—a model parameter. The total number of receptors for each hepatocyte was assumed to be constant. Intracellular IP₃ concentration (IP3) balance was modeled using saturation kinetics for synthesis influenced by hormone binding to the receptor and cytosolic Ca²⁺ and mass action kinetics for degradation (Equation 5). IP₃ exchange between adjacent hepatocytes was modeled as a mass transfer term assuming fast kinetics, with G_ij being the mass transfer coefficient. Increase in cytosolic Ca²⁺ (CaI) in the model was regulated IP₃R (g) activation, cytosolic IP₃ levels, store Ca²⁺ content (CaT), and a constant leakage from the ER (L), whereas decrease in cytosolic Ca²⁺ caused by SERCA pump activity was modeled as a Hill function (Equation 6). IP₃R activation (g) in the model was regulated by cytosolic Ca²⁺, whereas a constitutive rate of IP₃R was considered (Equation 7). The model assumed constant total intracellular Ca²⁺ for all hepatocytes. Additional details of model development can be found in Verma et al. (2016). See Tables 1, 2 for parameter descriptions, their nominal ranges and initial values for model species. All simulations in this study were performed using Matlab (version 8.1.0.604 (R2013a) Mathworks, Natick, MA).

TABLE 1

Table 1. List of species and their initial values in the computational model.

TABLE 2

Table 2. List of nominal parameter values/ranges for the computational model.

To identify the effects of non-uniformity of gap junction conductivity between adjacent hepatocytes in our simulations, G_ij values were sampled as follows: a uniform random number r₁ ϵ [0, 1] was drawn. If r₁ exceeded a threshold value p_th (two cases considered: p_th = 0.2 or 0.5), a G_ij was sampled ϵ [0.5, 0.9]. Otherwise G_ij = 0. p_th = 0.2 or 0.5 correspond to cases where 20% or 50% G_ij values are likely to be 0 respectively.

Model Reproducibility and Comparison of Alternatives

Simulation results presented in the current work were reproduced independently using the parameter values and hepatocyte adjacency information provided as Supplementary Material with this manuscript. While the original model was implemented in Matlab as a sequentially updating model species according to their specific rate equations, the rate equations in the reproduced model were implemented as a matrix. The Matlab code for the two independent implementations is provided in the Supplementary Material). The simulation results of the two model implementations were in agreement (see Figure S3 for details).

We also considered an alternative modeling scheme, in which the store Ca²⁺ content of each hepatocyte is considered to be a constant (= 500 μM). In this alternative model, Equation (8) is excluded and Equation (6) was changed as follows:

$\begin{array}{rcll}\frac{dCa{I}_i}{dt}& =& (1-g_i)\left(\frac{A{(\frac{IP{3}_i}{2})}^4}{{(k_1+\frac{IP{3}_i}{2})}^4}+L\right)(500-Ca{I}_i)& (\mathrm{6a})\end{array}$

$\begin{array}{rc}-& \frac{B.Ca{I}_i^2}{k_2^2+Ca{I}_i^2}\end{array}$

Results

We acquired a dataset consisting of cytosolic Ca²⁺ dynamics in 1300 hepatocytes across different liver lobules over a period of 1600 s (Figure 1). The Ca²⁺ transients within the lobules were induced by a sustained vasopressin stimulus (see section Methods). At low vasopressin stimulus levels (0.1–0.5 nM), hepatocytes in intact mouse livers did not exhibit sustained cytosolic Ca²⁺ spikes (Figure S5). Vasopressin levels to which cells were exposed during the experiment were varied from 0.5 to 1 nM. The stimulus time profile is shown in Figure 1A. We used a step-wise increasing stimulus profile to identify cell-cell interactions that remain unaffected by a change in stimulus strength. Ca²⁺ response profiles for all hepatocytes in our data suggested an overarching synchronized response (Figure 1C). Cytosolic Ca²⁺ spikes as well as Ca²⁺ wave propagation through a lobular section bounded by a central vein and a portal triad are shown in Figures 1D,E, respectively. Hepatocytes generally exhibit asynchronous cytosolic Ca²⁺ spiking behavior superimposed on propagating Ca²⁺ waves.

FIGURE 1

Figure 1. (A) vasopressin stimulation profile used during the course of the experiment; (B) Left: 2D cross section of an idealized liver lobule with locations of portal triads and central vein. Blood flows from the periportal region toward the pericentral region. Right: Manually segmented hepatocytes with approximate locations of central veins and portal triads in the imaged area. Different colors represent different manually segmented hepatocytes. The hepatocytes are distributed across different liver lobules; (C) heat map of cytosolic calcium intensities across all 1,300 segmented hepatocytes over 1,600 s. The data suggests an overall synchronization of calcium spikes across the imaged area; (D) Sequential cytosolic Ca²⁺ spiking in 10 hepatocytes residing in the region between a central vein and a portal triad. All hepatocytes are in the region within the white box in (B). Propagating Ca²⁺ waves interspersed with asynchronous spiking can be seen; (E) Ca²⁺ wave propagating away from a central vein.

Correlation Network Analysis

Our data suggested an overall synchronization of intracellular Ca²⁺ dynamics across all 1300 hepatocytes that were segmented within the imaged field even though these hepatocytes were often separated by several cell layers. With a correlation-based network analysis, we sought to identify the typical spatial range within with Ca²⁺ responses of individual hepatocytes are synchronized under the experimental conditions. The correlation networks were constructed using pairwise Spearman rank correlation coefficients between Ca²⁺ traces. We used a minimum correlation coefficient cutoff (R_th) of 0.75 and maximum inter-hepatocyte distance cutoff (d_th) of 100 μm to assign network edges between two hepatocytes. The resulting network consisted of 669 hepatocytes with 565 edges between them. The node degree distribution for all hepatocytes in the network suggested that a large number of hepatocytes were synchronized with one or two other hepatocytes, for the chosen R_th and d_th values (Figure 2A). In order to identify the typical spatial extent of Ca²⁺ response synchronization among hepatocytes, we decomposed the network into isolated clusters. We found a set of 14 clusters containing more than 8 hepatocytes in each cluster (Figure 2B).

FIGURE 2

Figure 2. Cluster sizes and degree distribution for correlation-based network. (A) Node degree distribution for all nodes in the correlation network. A majority of nodes exhibit a degree < = 2; (B) cluster size distribution for all isolated clusters. Most clusters consist of 2 nodes and 1 edge. 14 clusters consisted of > = 8 hepatocytes (inset); (C) Node degree distribution for hepatocytes residing in large clusters found in the correlation network. Degree of synchronization for hepatocytes with their neighbors is frequently > = 3; (D) Average node degree distribution for large clusters. Nine out of Fourteen clusters show average degree > 2; (E) Clusters in correlation network mapped to their physical locations. Hepatocytes are represented as circles centered at their locations in the imaged field. Red stars mark the approximate locations of central veins (CV) in the imaged area. Hepatocytes belonging to a cluster have been plotted in the same color. Synchronized “islands” of hepatocytes cover only small regions of the imaged area.

We focused our analysis on clusters that consisted of at least 8 hepatocytes (herein referred to as large clusters) for further analysis. For each large cluster, we estimated node degrees for individual hepatocytes as well as the average node degree for all hepatocytes within the cluster. Node degrees of individual hepatocytes residing in large clusters ranged between 1 and 7 (Figure 2C). Eighty seven out of the 190 hepatocytes residing in the large clusters had node degrees 3 or greater. Nine out of the 14 large clusters exhibited average node degrees > 2 (Figure 2D).

Mapping of the large clusters onto their physical locations (Figure 2E) suggested the existence of “islands” of synchronized Ca²⁺ response, which were generally situated close to the central veins. The typical spatial dimension of these synchronized clusters was less than the lobular dimensions (considered to be half the typical distance between approximate locations of two central veins). It must be noted that other smaller clusters, which are not shown in the figure, may be present in the intermittent space between the larger clusters.

In summary, analysis of correlation networks constructed based on hepatocyte Ca²⁺ response showed that, despite an apparent global concurrence of Ca²⁺ peak intensities, synchronized hepatocytes formed localized clusters spanning small regions within liver lobules. Ca²⁺ responses of up to 7 hepatocytes were synchronized within these clusters. However, only a small fraction of hepatocytes was included in clusters with sizes 8 or greater (Figure 3C) suggesting that a correlation-based formulation of cell-cell interactions is insufficient to explain the observed lobular scale propagation of Ca²⁺ waves.

FIGURE 3

Figure 3. Cluster sizes and degree distribution for causal network. (A) Total node degree distribution of all hepatocytes in the network; (B) Cluster size distribution of all and large clusters (clusters with sizes > = 8, inset). Cluster sizes are much higher than those in correlation network analysis; (C) Total node degree distribution for all hepatocytes in large clusters. Hepatocytes are causally connected to up to 8 neighbors.

Causal Network Analysis

We constructed causal networks between hepatocytes to identify whether adjacent neighbor-driven intercellular interactions can account for lobular scale Ca²⁺ wave propagation. We used a cell-oriented, transfer entropy (TE) based approach to identify causal connectivity between neighboring hepatocytes (see Methods). We considered molecular exchange between hepatocytes as the basis of causal influence between spatially co-localized hepatocytes and therefore allowed causal edges to exist only between closest neighbors. In our analysis, a unidirectional alignment of causal edges would suggest an organized, wave-like information flow along hepatocytes. We analyzed the resulting causal network for average node degree, total node degree, in and out node degrees, and direction of causal edges, to identify how many neighbors typically influence a given hepatocyte and the directional orientation of cell-cell interactions.

The causal network comprised of 1,162 hepatocytes with 1,491 edges between them. The number of hepatocytes included in the causal network far exceeded that in the correlation network (669 hepatocytes with 565 edges) suggesting that causality analysis describes intercellular interactions between neighboring hepatocytes better than a correlation-based analysis. We analyzed total node degree and cluster sizes for all nodes in the graph. The total node degree distribution of the all hepatocytes in the network peaked at a value of 2, pointing to causal connections between a given hepatocyte and multiple neighbors (Figure 3A). Upon decomposing the causal network into isolated clusters, we found 19 large clusters (cluster size ≥ 8). However, unlike the large clusters in the correlation network, large clusters in the causal network consisted of a higher number of hepatocytes (up to 160 hepatocytes, Figure 3B). The total node degree distribution for all hepatocytes residing in large clusters (sizes ≥ 8) exhibited similar causal connectivity characteristics between hepatocytes and their neighbors (Figure 3C).

Figure 4A shows large isolated clusters in the causal network mapped to their physical locations within the imaged slice. The large clusters contain 929 of the 1,300 hepatocytes in the imaged area. In contrast to the correlation network, large clusters within the causal network span much larger areas of lobules compared to correlation network clusters. A visualization of the direction of causal influence between hepatocytes residing in a large cluster is shown in Figure 4B. Our analysis suggested that although hepatocytes were causally connected with a number of neighbors ranging from 1 to 8, the direction of causal influence was not consistently from the pericentral region to the periportal region.

FIGURE 4

Figure 4. Spatial mapping of large causal networks. (A) Each cluster is represented by a different color. Large clusters (sizes > = 8) in causal network include 71.4% of the segmented hepatocytes (929 out of 1,300 hepatocytes); (B) causal influence edges in a representative cluster. Green arrows represent causal edges from the pericentral to the periportal region, whereas red arrows represent causal edges from the periportal region to the pericentral region. Causal influences identified between hepatocyte pairs are not aligned in a unidirectional fashion in the expected pericentral to periportal direction.

Computational Model-Based Analysis of Spatial Ca²⁺ Wave Propagation Patterns

We next sought to determine whether a combination of the causal influence network and a dynamic model of hepatocyte Ca²⁺ response can yield propagating Ca²⁺ waves consistent with experimental observations. We started with a previously-published dynamic model of hepatocyte Ca²⁺ response (Verma et al., 2016) and extended the model to incorporate cell-cell connectivity suggested by the causal influence network (see Methods). In addition, we modified the model parameters to incorporate zonation patterns of signaling components based on results from recently published single cell RNA-seq study (Halpern et al., 2017). We mined the transcriptomic data set (Table S3 from Halpern et al., 2017) to identify zonation of mRNA expression of Ca²⁺ signaling relevant genes. Specifically, we considered the zonation patterns of arginine-vasopressin receptor 1a (Avpr1a) and Phospholipase C β-1 (Plcb1). Zonation profiles for Avpr1a and Plcb1 in the data from (Halpern et al., 2017) are shown in Figure 5A. Avpr1a expression levels and Plcb1 expression levels correspond to vasopressin receptor recycling rate (model parameter k_r), and IP₃ synthesis rate (model parameter k_IP3), respectively, in the present dynamic model. For the subsequent analysis using integrated causal network and dynamic modeling, we considered a large cluster of hepatocytes identified using the causal network analysis. Experimentally determined Ca²⁺ patterns in this cluster of hepatocytes are shown in Figure 5B. Notable features of Ca²⁺ wave propagation through the cluster were: (1) Ca²⁺ waves propagated through the cluster from the pericentral region toward the periportal region consistent with the prior expectation, and (2) Ca²⁺ waves started from multiple hepatocytes located closer to the approximate location of the central vein residing closest to the cluster. We evaluated the dynamic model of this large cluster to identify the spatial patterns of intracellular signaling components as well as gap junction connectivity patterns that are consistent with experimentally observed Ca²⁺ wave propagation. In the dynamic model, Ca²⁺ response coupling is caused by gap junction mediated IP₃ exchange, a phenomenon that has been reported previously (Tordjmann et al., 1997; Eugenín et al., 1998).

FIGURE 5

Figure 5. Ca²⁺ waves in experimental data. (A) zonation of Avpr1a (arginine—vasopressin receptor 1-a) and Plcb1 (PLC beta 1, linked to IP₃ production rate) gene expression profiles across liver lobules (data from Halpern et al., 2017). Means and standard errors for gene expression values (fraction of total cellular mRNA) of each gene across the 9 layers show decreasing expression patterns from the pericentral region (PC) toward the periportal region (PP). Layer 1 lies closest to a central vein (acinar zone 3) whereas layer 9 is the farthest away (acinar zone 1); (B) Ca²⁺ wave propagation through a cluster identified using TE-based analysis. Ca²⁺ waves are initiated by hepatocytes lying close to the bottom of each panel and propagate upward. The red star in the right panel represents the approximate location of the central vein closest to the cluster. Ca²⁺ waves start from multiple hepatocytes residing near the bottom of each panel.

We simulated the dynamic model to identify the effect of gap junctions on coupling of Ca²⁺ dynamics across hepatocytes. Simulations were performed using the spatial locations of hepatocytes for the cluster shown in Figure 5B. The connectivity structure of the causal influence network from the TE-based analysis was utilized as the adjacency matrix for cell-cell IP₃ exchange. We considered two modes of gap junction conductivity (model parameter G_ij): (1) no hepatocyte exchanges IP₃ with its neighbors, and (2) each hepatocyte exchanges IP₃ with all its neighbors. Gap junction conductivity parameter between any pair of hepatocytes, modeled as a mass transfer term assuming fast IP₃ diffusion kinetics, was set to either 0 (no IP₃ exchange) or 0.9.

In the first set of simulations, the individual hepatocyte-specific values of signaling parameters k_r and k_IP3 were sampled from uniform distributions over nominal parameter ranges listed in Verma et al. (2016) (k_r ϵ [1, 2] s⁻¹, k_IP3 ϵ [0.7, 0.9] × 10⁴ μMs⁻¹, Figure 6A). The corresponding simulation results demonstrate that lobular scale Ca²⁺ waves did not occur when the gap junction-mediated IP₃ exchange was turned off (Figure 6C). By contrast, Ca²⁺ waves propagated through the cluster, when each of the hepatocytes exchanged IP₃ with all their neighbors (Figure 6D). However, the direction of wave propagation was not necessarily consistent with the experimental observations (Figure 6B).

FIGURE 6

Figure 6. Effect of randomness of intracellular Ca²⁺ signaling parameters. Intracellular parameter values were randomly chosen from uniform distributions for each hepatocyte in the simulation (k_r ϵ [1, 2], k_IP3 ϵ [0.7, 0.9]). (A) k_r and k_IP3 values used for each hepatocyte; (B) Ca²⁺ waves in experimental data; (C) no gap junction mediated IP₃ exchange; and (D) gap junction mediated IP₃ exchange with all neighbors. With gap junctions switched off, there are no Ca²⁺ waves through the cluster (C). With gap junctions switched on (D), Ca²⁺ waves propagate, however, the direction of propagation is not consistent with that observed experimentally. The red stars show the approximate location of the closest central vein. Note that the time points shown in each case were selected to best depict Ca²⁺ dynamics and spatial propagation. The times shown in (B) are with reference to the experiment start time. Times shown in (C,D) were measured from the time when stimulus was introduced in the simulations (T = 200 s).

We next simulated the dynamic model with the values of parameters k_r and k_IP3 drawn from spatial profiles that mimicked the zonated gene expression levels observed experimentally in Halpern et al. (2017). We approximated spatial profiles for mean k_r and k_IP3 values as exponentially and linearly decreasing functions with increasing distance from central vein, respectively, confined within the nominal parameter ranges (Figure 7A; Table 1; Verma et al., 2016). Parameter values for all hepatocytes in the model were initialized based on their distance from the central vein with additive noise (see Figure S6 for a description of model parametrization). We evaluated the effect of changing gap junction conductivity, according to the two modes considered in simulations shown in Figure 6.

FIGURE 7

Figure 7. Effect of lobular gradients of intracellular Ca²⁺ signaling parameters. (A) k_r and k_IP3 values used for each hepatocyte. k_r values decay exponentially with increase in the distance of a hepatocyte from the central vein whereas k_IP3 values decrease linearly with increasing distance from the central vein; (B) Ca²⁺ waves in experimental data; (C) no gap junction mediated IP₃ exchange; and (D) gap junction mediated IP₃ exchange with all neighbors. No Ca²⁺ waves appear when gap junctions are switched off (C). With gap junctions switched on (D) Ca²⁺ waves propagate through the cluster. However, the Ca²⁺ waves start from only a few hepatocytes. A Ca²⁺ wave propagates through the cluster in ~ 30 s, as compared to ~ 50 s in the experiment. The red stars show the approximate location of the closest central vein. Note that the time points shown in each case were selected to best depict Ca²⁺ dynamics and spatial propagation. The times shown in (B) are with reference to the experiment start time. Times shown in (C,D) were measured from the time when stimulus was introduced in the simulations (T = 200 s).

The simulation results suggested that even with gradients in parameters governing receptor-mediated signaling and IP₃ synthesis, Ca²⁺ waves did not arise in the absence of molecular exchange (Figure 7C). In the simulations, propagating Ca²⁺ wave patterns consistent with experiments were observed when hepatocytes exchanged IP₃ with their neighbors (Figure 7B, D). Our simulation results differed from experiments with regards to the region where Ca²⁺ waves are initiated. In the experimental observations, Ca²⁺ waves started from hepatocytes spread out in a relatively wider area close to the central vein. This difference could be due to the fact that in our dynamic model, hepatocytes were parametrized based on their distance from a central vein approximated as a point, when in reality the pericentral hepatocyte phenotype might result from microenvironmental signals in a more diffused region surrounding the central veins, whose diameters could span a few cell layers.

We simulated our model to identify the effects of non-homogeneous gap junction conductivity by varying parameter G_ij. Heterogeneity in gap junction conductivity could account for variability in cell-cell contact areas and gap junctions themselves. We considered two modes of gap-junction non-uniformity where either 20 or 50% G_ij values were likely to be 0 to account for a fraction of hepatocyte pairs not interacting with each other. Additionally, the non-zero G_ij values in either case were randomly drawn from a uniform distribution [ϵ [0.5, 0.9]] to account for variability in gap junction conductivity and number between a pair of adjacent hepatocytes. Other cell-intrinsic parameter values were identical to those used in the simulations corresponding to Figure 7. Effect of gap junction heterogeneity on spatial patterns of Ca²⁺ signal propagation through a cluster of hepatocytes identified using the TE-based analysis are shown in Figure 8 (see Figure S7 for cell-cell connections). We observed that in our simulations Ca²⁺ waves propagated through the cluster despite 22.1% (Figure 8A) and 50% (Figure 8B) G_ij values set to 0. Consistent with expectation, Ca²⁺ waves propagated in the direction of intracellular parameter gradients in both cases. Our simulations suggested that multiplicity of hepatocyte interactions makes Ca²⁺ wave propagation robust to non-interacting hepatocyte pairs.

FIGURE 8

Figure 8. Effect of non-uniform gap junction conductivity on Ca²+ signal propagation. Normalized G_ij value histograms and Ca²⁺ wave propagation with 20% (p_th = 0.2, A) and 50% (p_th = 0.5, B) likelihood of a given IP₃ mass transfer term (G_ij) being 0. 22.1% G_ij values in (A) and 49.6% G_ij values in (B) are 0. In either case, Ca²⁺ waves propagate through the cluster. Note that the time points shown in each case were selected to best depict Ca²⁺ dynamics and spatial propagation. Times shown in (A,B) were measured from the time when stimulus was introduced in the simulations (T = 200 s).

Since the store Ca²⁺ concentration is nearly 1000 times higher than cytosolic Ca²⁺ concentration we considered an alternative model in which the store Ca²⁺ was considered to be a constant (see section Methods). In this model formulation, store Ca²⁺ could be interpreted as a driving force for influx of Ca²⁺ within the cytosol instead of being trafficked between the cytosol and the Ca²⁺ store. Ca²⁺ wave propagation characteristics in this case were similar to the case in which store Ca²⁺ was dynamic (see Figure S4).

Low vs. High Stimulus

We used a step-wise increasing vasopressin stimulus over the course of the live tissue imaging duration (Figure 1). To identify the effects of change in stimulus, we compared the induced correlation and causal networks present at different stimulus levels in our data. We divided the time course into low and high stimulus regimes based on overall increase in rates of cytosolic Ca²⁺ spikes and compared the cluster sizes and localization for correlation and causal networks (Figures S8–S10). Our analysis revealed that although there was an increase in the number of large clusters in the high stimulus regime, the size of the clusters decreased (Figure S8). In contrast, we found similar large clusters in the high stimulus regime with a moderate increase in cluster sizes present in the causal network. In either case, regions of the image spanned by the large clusters were dependent on the stimulus level (Figure S9).

Discussion

In this work, we analyzed Ca²⁺ signal propagation in a two-dimensional optical slice of a perfused and intact mouse liver at the lobular scale. We generated a large-scale data set on cytosolic Ca²⁺ responses of 1300 hepatocytes to hormonal stimuli over a period of 1600 s. We analyzed the synchronization of Ca²⁺ response of a large population of hepatocytes in intact liver using correlation analysis and TE-based causal network analysis to identify directional flow of causal influence across hepatocytes in a lobule. We employed a computational model-based analysis to identify spatial patterns of intracellular Ca²⁺ signaling components and gap junction conductivity that can yield lobular scale Ca²⁺ waves consistent with experimental observations.

Identification of functional networks within a population of colocalized cells has gained prominence in recent decades (Bullmore and Sporns, 2009; Ahrens et al., 2013; Tian et al., 2018). Correlation (Fox et al., 2005; StoŽer et al., 2013; Markovič et al., 2015) as well as causal (Lungarella and Sporns, 2006; Wollstadt et al., 2014; Seth et al., 2015) network analysis are viable strategies to identify functional connectivity in cell populations. Correlation networks are commonly used to analyze Ca²⁺ responses in a population of cells under a global stimulus. For instance, correlation networks constructed using Ca²⁺ dynamics in Islets of Langerhans exhibit stimulus-dependent synchronization characteristics when stimulated by glucose (StoŽer et al., 2013; Markovič et al., 2015). However, correlation network analysis was insufficient to explain lobular scale propagation of Ca²⁺ waves observed in our experiment. In contrast, causal network analysis of the experimental data elucidated prominent features of lobular scale Ca²⁺ wave propagation such as existence of “islands” of causally connected hepatocytes within liver lobules and lack of directional alignment of causal edges between hepatocytes from the pericentral region to the periportal region. Although causal network analysis yielded misaligned causal connections between hepatocytes residing in a cluster, it pointed toward zonation and intercellular communication as cell-level dynamics that yield lobular scale organization of Ca²⁺ response. Spatially organized heterogeneity leads to location-based differences in Ca²⁺ signaling capacity of hepatocytes. Intercellular communication results in entrainment of Ca²⁺ responses of adjacent hepatocytes which extends throughout liver lobules via local interactions to yield Ca²⁺ waves.

The computational model of Ca²⁺ dynamics used in our study is limited in its scope. Our dynamic model-based analysis, parametrized using lobular scale spatial patterns of liver gene expression, represents a specific case of a more generalized concept of functional gradients that control Ca²⁺ wave propagation in liver lobules. Functional zonation results from zonal differences in micro-RNA expression (Sekine et al., 2009; Chen and Verfaillie, 2014) as well as protein activity (Gebhardt and Mecke, 1983; Jungermann, 1987). Cytosolic Ca²⁺ spiking dynamics have previously been observed in rat hepatocytes due to activation of adrenergic (Rooney et al., 1990) and purinergic (Dixon et al., 1990) receptors. Of the wide array of cell surface receptors and extracellular stimuli that could be spatially organized in liver lobules, our model-based analysis considered zonation of Avpr1a and Plcβ1 only, and response to hormone-induced GPCR signaling cascade. Additionally, our deterministic model ignores stochasticity in cellular level phenomena. For example, we modeled gap junction mediated molecular exchange as a flux term wherein channel conductivity between a pair of hepatocytes remained constant over time. However, a probabilistic treatment of open and closed channels, possibly linked to intracellular Ca²⁺ signaling events (Török et al., 1997; Peracchia, 2004; De Vuyst et al., 2006), may capture cell-neighbor molecular interactions more accurately. Explicit consideration of a comprehensive intracellular Ca²⁺ signaling cascade with zonal information, a stochastic modeling framework, and integration of experimental data can potentially capture the complexity observed in lobular scale Ca²⁺ dynamics in the liver, such as lack of directionality of causal linkages between hepatocytes.

Although Ca²⁺ as well as IP₃ could be exchanged between neighboring hepatocytes through gap junctions and lead to Ca²⁺ efflux from intracellular stores, the effective diffusivity of IP₃ is higher than Ca²⁺ because Ca²⁺ is heavily buffered within hepatocytes (Allbritton et al., 1992). These observations suggest that IP₃ is strongly involved in coordinating Ca²⁺ responses at the lobular scale. In addition, a loss of wave-like propagation of Ca²⁺ signals has been shown upon disruption of cell-cell contacts using Ca²⁺ free buffer (Gaspers and Thomas, 2005). Intracellular Ca²⁺ mobilization could arise from paracrine ATP release and subsequent purinergic receptor activation. The relative contribution of Ca²⁺ response synchronization via gap junctions or paracrine ATP would depend on the tissue and zone-specific expression of Connexin subtypes and purinergic receptors. Disrupting cell-cell contacts between hepatocytes in perfused livers results in asynchronous Ca²⁺ spikes in hepatocytes under a vasopressin stimulus and the Ca²⁺ signals do not spread to neighbors (Gaspers and Thomas, 2005). These results suggest that the paracrine ATP release is not sufficient to drive a lobular scale Ca²⁺ signal propagation observed experimentally. That said, explicit consideration of other potential paracrine factors such as ATP will expand the scope and applicability to time scales of cell-cell interaction beyond the relatively fast timescale considered in this study.

We note that a variety of fluorophores are available for reporting intracellular calcium levels, including genetically encoded calcium reporters. For example, Fluo-8 AM and Rhod4 have been utilized for cytosolic calcium reporting in hepatocytes with small differences in Kd values and photostability (Lock et al., 2015). We have been using Fluo-8 AM with good success in previous studies (Bartlett et al., 2017) and therefore utilized this reporter for obtaining the dynamic data analyzed in the present study.

An important consideration in analyzing and modeling Ca²⁺ signal propagation within liver lobules is the three-dimensional arrangement of hepatocytes. Although the proximity of hepatocytes to either a portal triad or a central vein within a two-dimensional slice can be ascertained, information regarding the cellular adjacency and spatial localization along a third dimension is lacking in our experimental data. The lack of directional alignment of causal edges along a pericentral to periportal orientation could arise due to the presence of micro-environmental cues from other pericentral or periportal regions above or below the optical slice corresponding to the imaged area. Alternatively, multidirectional alignment of causal edges may be due to cell-autonomous Ca²⁺ responses of hepatocytes within a small region which appear independently of the global stimulus and do not propagate beyond a few cells. High spatial and temporal resolution imaging of three-dimensional tissue structure sufficient to study spatial organization of Ca²⁺ signaling in liver lobules remains a challenging problem. However, imaging techniques are constantly evolving to produce accurate three-dimensional reconstructions of tissues with high spatial resolution (Arganda-Carreras et al., 2010; Hoehme et al., 2010). Intra-vital imaging techniques for visualizing molecular dynamics in live animals (Benechet et al., 2017; Dunn and Ryan, 2017) could further augment our modeling efforts at small spatial scales. However, these methods would introduce new challenges such as lack of control over distribution of stimulus in the immediate microenvironment of hepatocytes. Application of a combination of three-dimensional reconstruction and intra-vital imaging may provide data with the high spatial and temporal resolution required for a detailed dynamic model-based accounting of Ca²⁺ signal propagation in liver lobules.

Author Contributions

AV, JH, and RV designed the study. AA conducted the experiments. AA, AN, JH, and RV analyzed experimental data and interpreted the results. AV performed causal network analysis. AV developed computational model and performed simulations. AV, BO, and RV analyzed and interpreted the network analysis and simulation results. AV and RV drafted the manuscript. AA, BO, JH, and RV edited the manuscript. JH supervised the experimental aspects. RV supervised the computational aspects.

Funding

This work was financially supported by National Institute of Biomedical Imaging and Bioengineering U01 EB023224, National Institute on Alcohol Abuse and Alcoholism R01 AA018873, and National Science Foundation EAGER 1747917. 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 Statement

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.

Acknowledgments

We thank Mr. Madhur Parihar for reproducing the computational modeling results presented in this study independently and providing the code for inclusion in the Supplementary Material.

Supplementary Material

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

References

Ahrens, M. B., Orger, M. B., Robson, D. N., Li, J. M., and Keller, P. J. (2013). Whole-brain functional imaging at cellular resolution using light-sheet microscopy. Nat. Methods 10, 413–20. doi: 10.1038/nmeth.2434

PubMed Abstract | CrossRef Full Text | Google Scholar

Allbritton, N. L., Meyer, T., and Stryer, L. (1992). Range of messenger action of calcium ion and inositol 1, 4, 5-trisphosphate. Science 258, 1812–1815.

PubMed Abstract | 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

Arganda-Carreras, I., Fernández-González, R., Mu-oz-Barrutia, A., and Ortiz-De-Solorzano, C. (2010). 3D reconstruction of histological sections: application to mammary gland tissue. Microsc. Res. Tech. 73, 1019–1029. doi: 10.1002/jemt.20829

PubMed Abstract | CrossRef Full Text | Google Scholar

Atri, A., Amundson, J., Clapham, D., and Sneyd, J. (1993). A single-pool model for intracellular calcium oscillations and waves in the Xenopus laevis oocyte. Biophys. J. 65, 1727–1739. doi: 10.1016/S0006-3495(93)81191-3

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

Bartlett, P. J., Gaspers, L. D., Pierobon, N., and Thomas, A. P. (2014). Calcium-dependent regulation of glucose homeostasis in the liver. Cell Calcium 55, 306–316. doi: 10.1016/j.ceca.2014.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Benechet, A. P., Ganzer, L., and Iannacone, M. (2017). Intravital microscopy analysis of hepatic T cell dynamics. Methods Mol. Biol. 1514, 49–61. doi: 10.1007/978-1-4939-6548-9_4

PubMed Abstract | CrossRef Full Text | Google Scholar

Braeuning, A., Ittrich, C., Köhle, C., Hailfinger, S., Bonin, M., Buchmann, A., et al. (2006). Differential gene expression in periportal and perivenous mouse hepatocytes. FEBS J. 273, 5051–5061. doi: 10.1111/j.1742-4658.2006.05503.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Bullmore, E., and Sporns, O. (2009). Complex brain networks: graph theoretical analysis of structural and functional systems. Nat. Rev. Neurosci. 10, 186–198. doi: 10.1038/nrn2575

PubMed Abstract | CrossRef Full Text | Google Scholar

Canaff, L., Petit, J. L., Kisiel, M., Watson, P. H., Gascon-Barré, M., and Hendy, G. N. (2001). Extracellular calcium-sensing receptor is expressed in rat hepatocytes coupling to intracellular calcium mobilization and stimulation of bile flow. J. Biol. Chem. 276, 4070–4079. doi: 10.1074/jbc.M009317200

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., and Verfaillie, C. M. (2014). MicroRNAs: the fine modulators of liver development and function. Liver Int. 34, 976–990. doi: 10.1111/liv.12496

PubMed Abstract | CrossRef Full Text | Google Scholar

De Vuyst, E., Decrock, E., Cabooter, L., Dubyak, G. R., Naus, C. C., Evans, W. H., et al. (2006). Intracellular calcium changes trigger connexin 32 hemichannel opening. EMBO J. 25, 34–44. doi: 10.1038/sj.emboj.7600908

PubMed Abstract | CrossRef Full Text | Google Scholar

Dixon, C. J., Woods, N. M., Cuthbertson, K. R., and Cobbold, P. H. (1990). Evidence for two Ca2⁺-mobilizing purinoceptors on rat hepatocytes. Biochem. J. 269, 499–502. doi: 10.1042/bj2690499

PubMed Abstract | CrossRef Full Text | Google Scholar

Dolmetsch, R. E., Xu, K., and Lewis, R. S. (1998). Calcium oscillations increase the efficiency and specificity of gene expression. Nature 392, 933–936.

PubMed Abstract | Google Scholar

Dunn, K. W., and Ryan, J. C. (2017). Using quantitative intravital multiphoton microscopy to dissect hepatic transport in rats. Methods 128, 40–51. doi: 10.1016/j.ymeth.2017.04.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Eugenín, E. A., González, H., Sáez, C. G., and Sáez, J. C. (1998). Gap junctional communication coordinates vasopressin-induced glycogenolysis in rat hepatocytes. Am. J. Physiol. Gastrointest. Liver Physiol. 274, G1109–G1116. doi: 10.1152/ajpgi.1998.274.6.G1109

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Fox, M. D., Snyder, A. Z., Vincent, J. L., Corbetta, M., Van Essen, D. C., and Raichle, M. E. (2005). The human brain is intrinsically organized into dynamic, anticorrelated functional networks. Proc. Natl. Acad. Sci. U. S. A. 102, 9673–9678. doi: 10.1073/pnas.0504136102

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, S., Watkins, S. M., and Hotamisligil, G. S. (2012). The role of endoplasmic reticulum in hepatic lipid homeostasis and stress signaling. Cell Metab. 15, 623–634. doi: 10.1016/j.cmet.2012.03.007

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., and Thomas, A. P. (2008). Calcium-dependent activation of mitochondrial metabolism in mammalian cells. Methods 46, 224–232. doi: 10.1016/j.ymeth.2008.09.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Gebhardt, R., and Matz-Soja, M. (2014). Liver zonation: Novel aspects of its regulation and its impact on homeostasis. World J. Gastroenterol. 20, 8491– 8504. doi: 10.3748/wjg.v20.i26.8491

PubMed Abstract | CrossRef Full Text | Google Scholar

Gebhardt, R., and Mecke, D. (1983). Heterogeneous distribution of glutamine synthetase among rat liver parenchymal cells in situ and in primary culture. EMBO J. 2, 567–570. doi: 10.1002/j.1460-2075.1983.tb01464.x

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

Hoehme, S., Brulport, M., Bauer, A., Bedawy, E., Schormann, W., Hermes, M., et al. (2010). Prediction and validation of cell alignment along microvessels as order principle to restore tissue architecture in liver regeneration. Proc. Natl Acad. Sci. U. S. A. 107, 10371–10376. doi: 10.1073/pnas.0909374107

PubMed Abstract | CrossRef Full Text | Google Scholar

Jungermann, K. (1987). Metabolic zonation of liver parenchyma: significance for the regulation of glycogen metabolism, gluconeogenesis, and glycolysis. Diabetes Metab. Res. Rev. 3, 269–293. doi: 10.1002/dmr.5610030112

PubMed Abstract | CrossRef Full Text | Google Scholar

Keizer, J., and De Young, G. (1993). Effect of voltage-gated plasma membrane Ca2+ fluxes on IP3-linked Ca2⁺ oscillations. Cell Calcium 14, 397–410. doi: 10.1016/0143-4160(93)90044-7

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

Lagoudakis, L., Garcin, I., Julien, B., Nahum, K., Gomes, D. A., Combettes, L., et al. (2010). Cytosolic calcium regulates liver regeneration in the rat. Hepatology 52, 602–611. doi: 10.1002/hep.23673

PubMed Abstract | CrossRef Full Text | Google Scholar

Lizier, J. T. (2014). JIDT: an information-theoretic toolkit for studying the dynamics of complex systems. Front. Robot. AI 1:11. doi: 10.3389/frobt.2014.00011

CrossRef Full Text | Google Scholar

Lock, J. T., Parker, I., and Smith, I. F. (2015). A comparison of fluorescent Ca2⁺ indicators for imaging local Ca2⁺ signals in cultured cells. Cell Calcium 58, 638–648. doi: 10.1016/j.ceca.2015.10.003

CrossRef Full Text | Google Scholar

Lungarella, M., and Sporns, O. (2006). Mapping information flow in sensorimotor networks. PLoS Comput. Biol. 2:e144. doi: 10.1371/journal.pcbi.0020144

PubMed Abstract | CrossRef Full Text | Google Scholar

Markovič, R., StoŽer, A., Gosak, M., Dolenšek, J., Marhl, M., and Rupnik, M. S. (2015). Progressive glucose stimulation of islet beta cells reveals a transition from segregated to integrated modular functional connectivity patterns. Sci. Rep. 5:7845. doi: 10.1038/srep07845

PubMed Abstract | CrossRef Full Text | Google Scholar

McConkey, D. J., and Orrenius, S. (1997). The role of calcium in the regulation of apoptosis. Biochem. Biophys. Res. Commun. 239, 357–366. doi: 10.1006/bbrc.1997.7409

PubMed Abstract | CrossRef Full Text | Google Scholar

Nathanson, M. H., Burgstahler, A. D., Mennone, A., Fallon, M. B., Gonzalez, C. B., and Saez, J. C. (1995). Ca2+ waves are organized among hepatocytes in the intact organ. Am. J. Physiol. Gastrointest. Liver Physiol. 269, G167–G171. doi: 10.1152/ajpgi.1995.269.1.G167

PubMed Abstract | CrossRef Full Text | Google Scholar

Oliveira, A. G., Andrade, V. A., Guimarães, E. S., Florentino, R. M., Sousa, P. A., Marques, P. E., et al. (2015). Calcium signalling from the type I inositol 1, 4, 5-trisphosphate receptor is required at early phase of liver regeneration. Liver Int. 35, 1162–1171. doi: 10.1111/liv.12587

PubMed Abstract | CrossRef Full Text | Google Scholar

Peracchia, C. (2004). Chemical gating of gap junction channels: roles of calcium, pH and calmodulin. Biochim. Biophy. Acta 1662, 61–80. doi: 10.1016/j.bbamem.2003.10.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Pusl, T., and Nathanson, M. H. (2004). The role of inositol 1, 4, 5-trisphosphate receptors in the regulation of bile secretion in health and disease. Biochem. Biophys. Res. Commun. 322, 1318–1325. doi: 10.1016/j.bbrc.2004.08.036

PubMed Abstract | CrossRef Full Text | Google Scholar

R Core Team (2015). R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing 2014.

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

Rooney, T. A., Sass, E. J., and Thomas, A. P. (1990). Agonist-induced cytosolic calcium oscillations originate from a specific locus in single hepatocytes. J. Biol. Chem. 265, 10792–10796.

PubMed Abstract | Google Scholar

Sáez, J. C., Connor, J. A., Spray, D. C., and Bennett, M. V. (1989). Hepatocyte gap junctions are permeable to the second messenger, inositol 1, 4, 5-trisphosphate, and to calcium ions. Proc. Natl Acad. Sci. U. S. A. 86, 2708–2712. doi: 10.1073/pnas.86.8.2708

PubMed Abstract | CrossRef Full Text | Google Scholar

Schlosser, S. F., Burgstahler, A. D., and Nathanson, M. H. (1996). Isolated rat hepatocytes can signal to other hepatocytes and bile duct cells by release of nucleotides. Proc. Natl Acad. Sci. U. S. A. 93, 9948–9953. doi: 10.1073/pnas.93.18.9948

PubMed Abstract | CrossRef Full Text | Google Scholar

Schreiber, T. (2000). Measuring information transfer. Phys. Rev. Lett. 85:461. doi: 10.1103/PhysRevLett.85.461

CrossRef Full Text | Google Scholar

Schuster, S., Marhl, M., and Höfer, T. (2002). Modelling of simple and complex calcium oscillations. FEBS J. 269, 1333–1355. doi: 10.1046/j.0014-2956.2001.02720.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Sekine, S., Ogawa, R., Mcmanus, M. T., Kanai, Y., and Hebrok, M. (2009). Dicer is required for proper liver zonation. J. Pathol. 219, 365–372. doi: 10.1002/path.2606

PubMed Abstract | CrossRef Full Text | Google Scholar

Serradeil-Le Gal, C., Jouneaux, C., Sanchez-Bueno, A., Raufaste, D., Roche, B., Preaux, A. M., et al. (1991). Endothelin action in rat liver. Receptors, free Ca2⁺ oscillations, and activation of glycogenolysis. J. Clin. Invest. 87, 133–138.

PubMed Abstract | Google Scholar

Seth, A. K., Barrett, A. B., and Barnett, L. (2015). Granger causality analysis in neuroscience and neuroimaging. J. Neurosci. 35, 3293–3297. doi: 10.1523/JNEUROSCI.4399-14.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Smedler, E., and Uhlén, P. (2014). Frequency decoding of calcium oscillations. Biochim. Biophy. Acta 1840, 964–969. doi: 10.1016/j.bbagen.2013.11.015

PubMed Abstract | CrossRef Full Text | Google Scholar

StoŽer, A., Gosak, M., Dolenšek, J., Perc, M., Marhl, M., Rupnik, M. S., et al. (2013). Functional connectivity in islets of Langerhans from mouse pancreas tissue slices. PLoS Comput. Biol. 9:e1002923. doi: 10.1371/journal.pcbi.1002923

PubMed Abstract | CrossRef Full Text | Google Scholar

Thomas, A. P., Bird, G., Hajnoczky, G., Robb-Gaspers, L. D., and Putney, J. W. (1996). Spatial and temporal aspects of cellular calcium signaling. FASEB J. 10, 1505–1517. doi: 10.1096/fasebj.10.13.8940296

PubMed Abstract | CrossRef Full Text | Google Scholar

Tian, L., Li, Q., Wang, C., and Yu, J. (2018). Changes in dynamic functional connections with aging. Neuroimage 172, 31–39. doi: 10.1016/j.neuroimage.2018.01.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Tordjmann, T., Berthon, B., Claret, M., and Combettes, L. (1997). Coordinated intercellular calcium waves induced by noradrenaline in rat hepatocytes: dual control by gap junction permeability and agonist. EMBO J. 16, 5398–5407. doi: 10.1093/emboj/16.17.5398

PubMed Abstract | CrossRef Full Text | Google Scholar

Török, K., Stauffer, K., and Evans, W. H. (1997). Connexin 32 of gap junctions contains two cytoplasmic calmodulin-binding domains. Biochem. J. 326, 479–483. doi: 10.1042/bj3260479

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 Ca 2+ 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

Wollstadt, P., Martínez-Zarzuela, M., Vicente, R., Díaz-Pernas, F. J., and Wibral, M. (2014). Efficient transfer entropy analysis of non-stationary neural time series. PloS ONE 9:e102833. doi: 10.1371/journal.pone.0102833

PubMed Abstract | CrossRef Full Text | Google Scholar

Woods, N. M., Cuthbertson, K. R., and Cobbold, P. H. (1986). Repetitive transient rises in cytoplasmic free calcium in hormone-stimulated hepatocytes. Nature 319, 600–602.

PubMed Abstract | Google Scholar

Woods, N. M., Cuthbertson, K. S., and Cobbold, P. H. (1987). Agonist-induced oscillations in cytoplasmic free calcium concentration in single rat hepatocytes. Cell Calcium 8, 79–100. doi: 10.1016/0143-4160(87)90038-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, L., Luo, Y., Chen, T., Chen, F., Wang, T., and Hu, Q. (2008). Ca2+ oscillation frequency regulates agonist-stimulated gene expression in vascular endothelial cells. J. Cell Sci. 121, 2511–2518. doi: 10.1242/jcs.031997

PubMed Abstract | CrossRef Full Text | Google Scholar