Three-dimensional spatially resolved geometrical and functional models of human liver tissue reveal new aspects of NAFLD progression.

Early disease diagnosis is key to the effective treatment of diseases. Histopathological analysis of human biopsies is the gold standard to diagnose tissue alterations. However, this approach has low resolution and overlooks 3D (three-dimensional) structural changes resulting from functional alterations. Here, we applied multiphoton imaging, 3D digital reconstructions and computational simulations to generate spatially resolved geometrical and functional models of human liver tissue at different stages of non-alcoholic fatty liver disease (NAFLD). We identified a set of morphometric cellular and tissue parameters correlated with disease progression, and discover profound topological defects in the 3D bile canalicular (BC) network. Personalized biliary fluid dynamic simulations predicted an increased pericentral biliary pressure and micro-cholestasis, consistent with elevated cholestatic biomarkers in patients' sera. Our spatially resolved models of human liver tissue can contribute to high-definition medicine by identifying quantitative multiparametric cellular and tissue signatures to define disease progression and provide new insights into NAFLD pathophysiology.

Introduction

High definition medicine is emerging as an integrated approach to profile and restore the health of an individual using a pipeline of multi-parametric analytical and therapeutic technologies[1]. It relies on large data sets, e.g. genomics, metabolomics, as well as imaging and computational modelling approaches to identify functional and structural abnormalities in organs and tissues associated with a disease state. Histology remains the method of choice to characterize pathological alterations of tissue structure[2]. However, this technique has several disadvantages, e.g. it is subjective (depends on the pathologist’s skills), is often semi-quantitative and provides only two-dimensional (2D) information[3]. In recent years, an increasing number of studies have highlighted the importance of considering three-dimensional (3D) information for the histopathological examination of tissues[4,5]. The liver is a pertinent example of complex 3D tissue organization[6]. It consists of functional units, the liver lobuli[7], containing two intertwined networks, the sinusoids for blood flow and the bile canaliculi (BC) for bile secretion and flux[7]. Sinusoids and BC run antiparallel along the central vein (CV)-portal vein (PV) axis. The hepatocytes are the major parenchymal cells and display a peculiar and unique type of cell polarity distinct from that of simple epithelia[8]. Whereas in simple epithelia the apical surfaces of all cells face the lumen of the organ, hepatocytes are sandwiched between the sinusoidal endothelial cells and share the apical surface with multiple neighbouring hepatocytes to form a 3D BC network[9]. Such an architecture makes it difficult to grasp the 3D organization of cells, BC and sinusoidal networks, and overall tissue structure from 2D histological sections. Recent advances in optical clearing and multi-photon microscopy allow imaging thick sections of tissues such that 3D information can be captured[10,11]. The image stacks can be processed to generate 3D digital reconstructions, i.e. geometrical models[6], with sub-cellular resolution. The geometric models provide a detailed quantitative description of the different cells and micro-structures forming the tissue and can be used to generate predictive models of tissue function e.g. biliary fluid dynamics[12], thus gaining novel insights into liver tissue organization and function. Thereby, geometrical models can be used to improve our understanding of liver (patho)biology.

Non-alcoholic fatty liver disease (NAFLD), defined as an accumulation of triglycerides and lipid droplets (LD) in the liver in absence of alcohol intake, is rising to become the most common chronic liver disease worldwide[13]. NAFLD includes a spectrum of liver diseases, ranging from simple steatosis to non-alcoholic steatohepatitis (NASH)[13]. Whereas steatosis is considered as a “non-progressive” status of the disease, NASH has the potential to progress to more severe stages, such as cirrhosis and hepatocellular carcinoma, leading eventually to liver failure and transplantation[3,13]. Thus, the understanding of the transition from steatosis (STEA) to early NASH (eNASH)[14], as a disease-defining moment for NAFLD prognosis is key to a deeper understanding of disease pathophysiology. Liver biopsy and histological inspection of thin tissue slices (<10 μm) reveal morphometric features such as hepatocyte ballooning and LD content, and constitute the current gold standard for the diagnosis of steatosis and NASH[2,13,15]. However, due to the limitations of histological analysis, alterations in 3D tissue structures such as BC[3,13], have been overlooked. In this study, we generated 3D spatially resolved geometrical and functional models of human liver tissue for different stages of NAFLD to contribute to a high definition medical diagnosis of disease establishment and progression.

Results

Three-dimensional geometrical models of human liver tissue

To quantitatively characterize the transition from simple STEA to eNASH, we stained, imaged and digitally reconstructed human liver tissue in 2D and 3D from biopsies of 25 patients classified into four groups: normal control (NC, n = 6), healthy obese (HO, n = 4), steatosis (STEA, n = 8) and eNASH (n = 7). The demographic, clinical and histological details of the patients are summarized in Table 1 and Supplementary Table 1. We focused on lobule size, cell and nuclear morphology, LD, and tissue features such as BC and sinusoidal networks. For this we tested 27 antibodies combinatorially (Supplementary Table 2). 2D analysis of full liver slices showed an increase in median lobule radius from 513±27 μm in NC to 581±19 μm in eNASH (Extended Data 1a-b), in agreement with previous reports[16]. Our 3D liver tissue reconstruction and analysis pipeline include the following steps. First, 100 μm liver sections were antigen retrieved using citric acid buffer (CAAR), stained for BC (CD13), sinusoids (fibronectin), nuclei (DAPI), LD (BODIPY) and cell borders (LDLR), optically cleared and imaged at high resolution using multiphoton microscopy (Extended Data 1c-f). Because this protocol did not provide a good cell border staining of the pericentral hepatocytes in STEA and eNASH, for cell-based measurements (see below Fig. 3), we used an alternative protocol without antigen retrieval enabling the staining of sinusoids (fibronectin), nuclei (DAPI), LD (BODIPY) and cell borders (phalloidin). All images cover one complete CV-PV axis within a liver lobule. The axis was oriented according to the direction of bile flow[12]. Finally, we reconstructed the tissue using our open-source Motion Tracking software, as described[6] (Fig. 1 and Supplementary Video 1-2). The generation of geometrical models constituted the basis for the quantitative and structural characterization of the liver tissue in NAFLD biopsies.

Table 1

Summary of analysed samples.

The number of patients in each phenotypic category is provided together with demographic and histologic characteristics and blood parameters. All numeric traits are shown as the median with the range provided in parenthesis. BMI, body mass index; NAS, NAFLD activity score.

	Control (NC)	Healthy obese (HO)	Steatosis (STEA)	early NASH (eNASH)
N	6	4	8	7
Male %	33	25	63	14
Age, years	69 (54-85)	36.5 (29-68)	42 (34-51)	51 (39-58)
BMI	25 (21-27)	40.5 (32-45)	45.5 (39-60)	51 (45-75)
NAS-Score	0 (0-1)	0.5 (0-1)	3 (1-3)	4 (3-5)
NAS-Fat	0 (0-1)	0.5 (0-1)	2.5 (1-3)	2 (2-3)
NAS-Balloning	0	0	0 (0-1)	1 (0-1)
NAS-Inflammation	0	0	0	1 (1-1)
Fat content %	0 (0-8)	2.5 (0-5)	60 (25-80)	60 (40-80)
Fibrosis	0	0	0	0 (0-1)
GGT (U/l)	27.5 (18-60)	19 (17-20)	57 (21-112)	77 (28-139)
AP (U/l)	70 (63-118)	55 (50-59)	76 (60-90)	86 (59-125)
Bilirubin (μmol/l)	7 (0-14)	5 (3-6)	6.5 (5-8)	8 (6-18)
ALT (U/l)	25.5 (15-51)	25 (15-27)	34 (23-55)	42 (18-71)
AST (U/l)	21 (11-44)	15 (10-25)	61 (25-80)	57 (19-76)

Extended Data Fig. 1

Immunofluorescence of human liver tissue.

a, Human liver sections were stained for glutathione synthetase (GS) to visualize CV and DAPI. Scale bar, 1000 μm. Representative images from NC = 4 samples and eNASH = 5 samples. b, 2D analysis of liver lobule radius represented by box-plots (median values as red lines, 25th and 75th percentiles as blue bottom and top edges of the boxes, extreme data points by whiskers). NC = 4 samples, HO = 4 samples, STEA = 7 samples, eNASH = 5 samples. One-sided hypothesis test. *p-values < 0.05, **p-values < 0.01, ***p-values < 0.001. c-f, Liver sections (~100 um thick) were stained for bile canaliculi (CD13), sinusoids (fibronectin), nucleus (DAPI), lipid droplets (BODIPY) and cell border (LDLR), optically cleared with SeeDB and imaged at high resolution using multiphoton microscopy (0.3 μm x 0.3 μm x 0.3 μm per voxel). Orthogonal view of NC (c), HO (d), STEA (e) and eNASH (f). Scale bar, 50 μm. Representative images from NC = 5 samples, HO = 3 samples, STEA = 4 samples, eNASH = 4 samples.

Figure 3

Cell based analysis of NAFLD.

Quantification of the number of hepatocytes per tissue volume unit (a), number of hepatocytes per 100 μm lobule section (b) and cell volume (c) along the liver lobule and the overall average. d, Cell volume distribution. For the population analysis, the hepatocytes from all the groups were pulled together and the populations were defined based on their volume distribution. By fitting the volume distribution with two log-normal distributions, the volume values defining three populations’ boundaries were identified: small (< 5800 μm3), medium (5800 – 11000 μm3) and large (> 11000 μm3) (d). The percentage of cellular volume occupied by the different populations is shown in (e, f, and g). Percentage of the cell volume occupied by LD: distribution (h) and statistics along the CV-PV axis and overall (i). Hepatocytes with percentage of LD volume lower than 0.001% are not presented in the distributions, which were normalized such that their integrals are equal to 1000 (h). 11278 cells from 16 reconstructions (NC = 5 samples, HO = 3 samples, STEA = 4 samples, eNASH = 4 samples) were analysed. Spatially-resolved quantification represented by median ± MAD per region and overall quantifications by box-plots (median values as red lines, 25th and 75th percentiles as blue bottom and top edges of the boxes, extreme data points by whiskers). One-tailed hypothesis test. *p-values < 0.05, **p-values < 0.01, ***p-values < 0.001. j, Representative cells reconstructed in 3D and selected from region 3 and 8. Apical, basal and lateral surface are shown in green, magenta and grey, respectively. LD are shown in red. Scale bar, 10 μm. 16 reconstructions (NC = 5 samples, HO = 3 samples, STEA = 4 samples, eNASH = 4 samples) were repeated independently with similar results.

Fig. 1

3D reconstruction and quantitative analysis of human liver morphology.

Human liver sections obtained by biopsy (~100 um thick) were stained for bile canaliculi (CD13), sinusoids (fibronectin), nucleus (DAPI), lipid droplets (BODIPY) and cell border (LDLR), optically cleared with SeeDB and imaged at high resolution using multiphoton microscopy (0.3 μm x 0.3 μm x 0.3 μm per voxel). For each sample, the central vein (light blue), portal vein (orange), bile canaliculus (green), sinusoids (magenta), lipid droplets (red), nuclei (random colours) and hepatocytes (random colours) were reconstructed. a, Normal control (NC). b, Healthy obese (HO). c, Steatosis (STEA). d, Early NASH (eNASH). Scale bar 30μm. Representative reconstructions from NC = 5 samples, HO = 3 samples, STEA = 4 samples, eNASH = 4 samples.

Nuclear-based analysis of NAFLD progression

We first quantified properties of hepatocytes nuclei, such as cell nuclearity and ploidy, since hepatocytes are heterogenous[6,17]. We quantified nuclear vacuolation/glycogenation given it is a common histological characteristic in NAFLD[18]. Finally, we measured nuclear texture homogeneity, a feature associated with various pathological conditions[19,20], methylation and acetylation status[21] and, more recently, transcriptional activity[22]. We neither observed differences in the proportion of mono/binuclear cells nor in the ploidy between the groups (Extended Data 2b-c). The average values of several nuclear features showed only modest variations (Extended Data 2). However, many functional and morphological features of the liver change along the CV-PV axis, such as metabolic zonation[23], ploidy[6], BC[12] and lipid composition[24]. Therefore, to account for potential morphological heterogeneities along the CV-PV axis, we computationally divided it into ten equidistant regions (Extended Data 2d). We found major differences in nuclear elongation around the CV and increased vacuolated nuclei around the PV (Extended Data 2a,e,f) as previously reported[25]. The biological significance of nuclear vacuolation is unclear, but it correlates with hepatocyte senescence[26]. Moreover, we identified spatially distributed differences in nuclear homogeneity as disease progresses (Extended Data 2g-j). Therefore, our analysis reveals that the combined spatially distributed values of nuclear vacuolation and texture homogeneity could contribute to differentiate patient groups.

Extended Data Fig. 2

Morphometric features of the nuclei.

a, Representative IF images of fixed human liver tissue sections stained with DAPI. Shown is a single-plane covering an entire CV-PV axis. Arrowhead indicates some examples of vacuolated nuclei. Representative images from NC = 5 samples, HO = 3 samples, STEA = 4 samples, eNASH = 4 samples. Quantitative characterization of hepatocytes nuclei with respect to the proportion of mono/binuclear cells (b) and ploidy (c). Only the four major populations (i.e. 1x2n, 1x4n, 2x2n and 2x4n), which account for >90% of the hepatocytes, are shown. d, Definition of the regions within the liver lobule. The CV-PV axis was divided in 10 equidistant regions. Regions 1 and 10 are adjacent to the CV and PV, respectively. Quantitative characterization of hepatocytes nuclear elongation (e) and texture based on their DAPI intensity (see Methods for details): nuclear vacuolation (f), homogeneity (Angular Second Moment) (g), local homogeneity (Inverse Difference Moment) (h), Contrast (i) and Entropy (j). NC = 5 samples, HO = 3 samples, STEA = 4 samples, eNASH = 4 samples. Spatially-resolved quantification represented by median ± MAD per region and overall quantifications by box-plots (median values as red lines, 25th and 75th percentiles as blue bottom and top edges of the boxes, extreme data points by whiskers). One-tailed hypothesis test. *p-values < 0.05, **p-values < 0.01, ***p-values < 0.001.

Morphometric parameters of LD correlate with disease progression

The finding that quantitative, spatially-resolved analysis of nuclear parameters can reveal changes that are not evident in an average estimate prompted us to re-evaluate the morphometric characterization of LD. Even though LD are the hallmark of NAFLD, a detailed quantitative description of size and spatial localization within the liver lobule has not been achieved yet. Contrary to traditional histology[3], immunostaining of thick tissue sections preserved most LD (Fig 2a). In agreement with the histopathological description of NAFLD[3,13], a major increase in LD between the second and fifth regions was observed in STEA and eNASH (Fig. 2b). However, the LD occupy a higher volume of the tissue in eNASH than STEA, consistent with[27]. It is known that the LD can present massive differences in size[28]. The LD in the human liver samples ranged from ~1 μm3 to 80,000 μm3 (Fig. 2c). We performed a population analysis based on their volume distributions (Fig. 2c). We defined three sub-populations of LD, small (< 8 μm3), medium (8 –1000 μm3) and large (> 1000 μm3)(Fig. 2d-f). Whereas small LD were evenly distributed along the CV-PV in all conditions, large LD were highly enriched towards the pericentral zone in STEA and eNASH, occupying up to 25% of the tissue volume (Fig. 2f). Medium LD occupy a higher volume of the tissue in the periportal area in eNASH than STEA (Fig. 2e).

Fig. 2

Quantitative characterization of LD along the CV-PV axis.

a, Representative IF images of fixed human liver tissue sections stained with BODIPY. Shown is a maximum projection of a 60 μm z-stack covering an entire CV-PV axis. NC = 5 samples, HO = 3 samples, STEA = 4 samples, eNASH = 4 samples were repeated independently with similar results. b, Quantification of the percentage of tissue volume occupied by LD along the CV-PV axis and overall values (i.e. over the whole CV-PV axis). c, LD volume distribution for each disease condition. The inset shows the difference of the normalized LD volume distribution for all conditions (HO+STEA+eNASH) and the one from the NC (red curve). By fitting this distribution with two log-normal distributions, we defined three LD populations: small (< 8 μm3), medium (8 – 1000 μm3) and large (> 1000 μm3). 122538 LD from 16 reconstructions (NC = 5 samples, HO = 3 samples, STEA = 4 samples, eNASH = 4 samples) were analysed. Each volume distribution was normalized such that their integrals are equal to 10000. Quantification of the percentage of tissue volume occupied by the LD along the CV-PV for (d) small, (e) medium and (f) large LD. NC = 5 samples, HO = 3 samples, STEA = 4 samples, eNASH = 4 samples. Spatially-resolved quantification represented by median ± MAD per region and overall quantifications by box-plots (median values as red lines, 25th and 75th percentiles as blue bottom and top edges of the boxes, extreme data points by whiskers). One-tailed hypothesis test. *p-values < 0.05, **p-values < 0.01, ***p-values < 0.001.

The increased LD size during disease progression is such that some LD in the pericentral region were even larger than an average hepatocyte. This leads to global changes in tissue structure. To quantify them, we measured the spatial distribution of cell density, number of hepatocytes per lobule cross-section, cell volume and percentage of cell volume occupied by LD. We found ~50% reduction in the number of hepatocytes located between the CV and the middle zone in STEA and eNASH (Fig. 3a). Next, we estimated the number of hepatocytes per lobule section using the density of hepatocytes (Fig. 3a) and the radius of the liver lobule in 2D (Extended Data 1b). We observed a decrease in the number of hepatocytes in HO and STEA, which partially reverted to normal in eNASH (Fig. 3b). The reduction in the number of hepatocytes could be due to apoptosis, a phenomenon reported in NAFLD[29]. This reduction was compensated by a massive increase in cell volume (Fig. 3c). Hepatocytes were twice larger than the average size (Fig. 3c), reaching values up to ~100,000 μm3 for STEA and eNASH (twenty times larger than a small hepatocyte) (Fig. 3d). A population analysis of the hepatocytes based on volume revealed a characteristic distribution of different cell populations along the liver lobule (Fig. 3d-g). STEA and eNASH were characterized predominantly by small and large hepatocytes which are anti-correlated along CV-PV axis (Fig. 3e-g). Even though cell density and cell volume were practically indistinguishable between STEA and eNASH, we observed a remarkable phenotype regarding the fraction of cell volume occupied by LD (Fig. 3h-j). In eNASH, hepatocytes accumulated LD even in the periportal region (Fig. 3i,j and Supplementary Video 3), suggesting that LD accumulation progressively extends to the PV as the disease progresses.

Altogether, these data reveal profound quantitative morphological disparities in cell size and LD content along the CV-PV axis. Specifically, the percentage of cell volume occupied by LD in the PV and CV areas holds the potential to discriminate between STEA and eNASH.

Alterations in apical protein trafficking

The massive presence of LD that occupy a large portion of the cytoplasm raises the question of whether apical transport of proteins in hepatocytes is affected. We analysed the localization of four apical proteins, aminopeptidase N (CD13), bile salt export pump (BSEP), multidrug resistant-associated protein-2 (MRP2) and dipeptidylpeptidase-4 (DPPIV). CD13, BSEP and MRP2 were correctly localized to the apical membrane in all conditions (Extended Data 1c-f, 3a-b). DPPIV was enriched on the apical membrane with a small fraction on the basal membrane in NC and HO (Extended Data 3c). Strikingly, DPPIV was redistributed to the lateral membrane in pericentral hepatocytes, showing already a clear trend in STEA that becomes statistically significant in eNASH, while retaining its normal localization in the periportal zone (Extended Data 3c-e). Considering that DPPIV follows the transcytotic route to the apical surface[30] whereas BSEP and MRP2 do not[31,32], the mislocalization of DPPIV suggests a possible disruption of this transport route in pericentral hepatocytes. This supports previous findings regarding the misregulation of membrane protein trafficking in NAFLD[33] and prompted us to evaluate whether the integrity and 3D organization of the BC could be affected during disease progression.

Extended Data Fig. 3

Mislocalization of DPPIV in pericentral hepatocytes in STEA and eNASH.

Representative confocal microscopy images of human liver sections stained for the apical markers BSEP (a), MRP2 (b) and DPPIV (c). Merged images of the apical markers, phalloidin and DAPI are shown in the right panels. Arrowhead indicates the lateral membrane. Scale bar, 10μm. NC = 3 samples, STEA = 4 samples, eNASH = 4 samples were repeated independently with similar results. d-e, Large field images of a single-plane of liver tissue stained with DPPIV (d). Scale bar, 50μm. Apical, basal and lateral membrane of the hepatocytes were segmented based on BSEP (not shown), DPPIV and phalloidin (not shown) in an area covering a radius of 125 μm around the CV and PV. DPPIV intensity was quantified and normalised to the area covered by the different sub-domains (e). NC = 3 samples, STEA = 4 samples, eNASH = 4 samples. Quantifications by box-plots (median values as red lines, 25th and 75th percentiles as blue bottom and top edges of the boxes, extreme data points by whiskers). One-tailed hypothesis test. *p-values < 0.05, **p-values < 0.01, ***p-values < 0.001.

Bile canaliculi network shows geometrical and topological defects in NAFLD

We carried out a geometrical and topological characterization of BC and sinusoidal networks. Even though there is a tendency towards a reduction in volume fraction and total length of the sinusoidal network in STEA and eNASH (Extended Data 4b,e), no major defects in sinusoidal microanatomy were detected (Extended Data 4a,c-d,f-h). In contrast, we detected profound quantitative alterations in the architecture of the BC network between patient groups. Contrary to the packed and homogeneous appearance in NC and HO, the BC in STEA and eNASH displayed clear morphological defects which were more pronounced in the pericentral region (Fig. 4a, b). A detailed analysis revealed a sustained increase in BC radius in STEA and eNASH especially pronounced toward the periportal region (Fig. 4d). In addition, in both STEA and eNASH, we observed a strong reduction in the total length and branches crossing regions of the BC towards the pericentral zone (Fig. 4c and i). Other geometrical properties such as volume fraction and junction density were unaffected (Fig. 4c and e). These changes are significant as we found very little variability in the geometric and topological features of the BC network among liver lobules of the same patient and among lobules within the same group of patients (Extended Data 5a-d).

Extended Data Fig. 4

Structural and topological characterization of the sinusoidal network.

a, Representative IF images of fixed human liver tissue sections stained with fibronectin after CAAR. Shown is a maximum projection of a 30 μm z-stack covering an entire CV-PV axis. Representative images from NC = 5 samples, HO = 3 samples, STEA = 5 samples, eNASH = 3 samples Quantification of the tissue volume fraction occupied by the sinusoids (b), radius (c), number of junctions (d), total length per unit tissue volume (e), fraction of connected network (f) connectivity density (g) and branches crossing regions (h) for the sinusoidal network along the CV-PV axis and overall. NC = 5 samples, HO = 3 samples, STEA = 5 samples, eNASH = 3 samples. Spatially-resolved quantification represented by median ± MAD per region and overall quantifications by box-plots (median values as red lines, 25th and 75th percentiles as blue bottom and top edges of the boxes, extreme data points by whiskers). Two-tailed hypothesis test. *p-values < 0.05, **p-values < 0.01, ***p-values < 0.001.

Figure 4

Structural and topological defects of bile canaliculi revealed by spatial 3D analysis.

a, Representative IF images of fixed human liver tissue sections stained with CD13 after citric acid antigen retrieval. Shown is a maximum projection of a 60 μm z-stack covering an entire CV-PV axis. b, Inset showing 3D representation of the BC highlighted in a. NC = 6 samples, HO = 4 samples, STEA = 8 samples, eNASH = 7 samples were repeated independently with similar results. Quantification of the volume fraction of tissue occupied by BC (c), radius (d), number of junctions (e), total length per volume (f), fraction of connected network (g) connectivity density (h) and branches crossing regions (i) of the BC network along the CV-PV axis and overall (See Methods for details). NC = 6 samples, HO = 4 samples, STEA = 8 samples, eNASH = 7 samples. Spatially-resolved quantification represented by median ± MAD per region and overall quantifications by box-plots (median values as red lines, 25th and 75th percentiles as blue bottom and top edges of the boxes, extreme data points by whiskers). Two-tailed hypothesis test. *p-values < 0.05, **p-values < 0.01, ***p-values < 0.001. j, Dependency of the predictive classification accuracy (10-k folds) on the number of parameters used by the classifier. The k-fold validation was performed 50 times. Whereas the blue dots represent the mean value, the error bars show the sem. The predictive accuracy is defined as the complement of the cross-validation loss of the model. k, 2D representation of the training set. Each sample is represented by its two Principal Components (explaining 80.8% of the point variability. Filled shapes show samples that were wrongly classified. Filling colour indicates the pathologist classification. 52 reconstructed BC networks were used for the classification analysis from healthy tissue (23 images of NC and HO), STEA (15 images) and eNASH (14 images). l, Confusion matrix of a 10k-fold prediction of the classifier showing the 'true' classes versus the predicted ones.

Extended Data Fig. 5

Geometric and topological variability of the BC network among liver lobules.

BC network was reconstructed from three CV-PV axes from different lobules for each patient. NC = 4 patients, HO = 3 patients STEA = 3 patients, eNASH = 3 patients. Quantification of the tissue volume fraction occupied by the BC, radius, total length per unit tissue volume and fraction of connected network (a-d) along the CV-PV axis and overall. Spatially-resolved quantification represented by median ± MAD per region and overall quantifications by box-plots (median values as red lines, 25th and 75th percentiles as blue bottom and top edges of the boxes, extreme data points by whiskers). Two-tailed hypothesis test. *p-values < 0.05, **p-values < 0.01, ***p-values < 0.001.

To investigate the topological properties of the BC network, we performed an analysis of network connectivity. Surprisingly, we found a pronounced decay in the connectivity in STEA and eNASH towards the pericentral region (Fig. 4a, b, g and h, Supplementary Video 4). One possibility is that the alterations in BC connectivity are only apparent: As the volume of hepatocytes is increased (Fig. 3g), the sectioning may cut both the BC and sinusoidal networks and, therefore, give only the impression of an altered connectivity. However, this is not the case for the sinusoidal network, which indeed does not appear disrupted (Extended Data 4a,f,g), supporting the idea that the BC network is specifically affected. Thus, our data point at specific geometrical and topological alterations in the pericentral BC network in both STEA and eNASH.

These results suggest that an unbiased classification of patients based on the quantitative morphometric and topological features of the BC network is feasible. To test for this, we analysed the prognostic power of the BC network parameters to classify patients in different disease stages using a machine learning framework. We used Support vector machine (SVM)[34], a supervised classifier successfully applied to several classification problems of disease conditions[35]. We analysed 52 reconstructed BC networks from healthy tissue (NC and HO), STEA and eNASH. We used 5 features (volume fraction, radius, length per volume, fraction of connected network and branches crossing regions) at 10 regions along the CV-PV axis (Figure 4 c-d, f-g), resulting in a set of 49 parameters. Remarkably, our analysis showed that a subset of just 7 parameters, which are only measurable in 3D, gives the maximum predictive accuracy (86.96±0.45 %, mean±sem) (Figure 4j). Whereas almost no errors were found in the discrimination of healthy tissue, few discrepancies were found between STEA and eNASH (Fig. 4k-l). These results suggest that a relatively small set of parameters describing the alterations of BC network may be sufficient to discriminate between disease stages.

Personalised model of bile flow predicts increase in bile pressure in the pericentral region

The observed alterations of BC network architecture are likely to have consequences for liver tissue function, particularly for bile flow. However, it is not yet possible to measure bile flow in the BC of human liver. We recently developed a computational model of bile fluid dynamics, validated its quantitative predictions in mouse models and demonstrated that bile velocity and bile pressure distributions along the liver lobule strongly depend on BC geometry[12]. Here, we extended this model[12] in a spatially heterogenous fashion (Fig. 5a) to handle the extreme inhomogeneity of BC density in STEA and eNASH. Briefly, the refined model is based on conservation of mass for water and osmolytes and Darcy’s Law for laminar flow. Since we obtained morphometric data from individual patients, we can now aim at developing personalised models, i.e. parameterized by integrating previously reported values (viscosity[36], permeability[37,38] and osmolyte secretion rate[39]) with individual geometrical and topological BC measurements (Fig. 4c,d,g, and Extended Data 6). No free parameters remained and, hence, no parameter fitting was needed. Next, we applied this model to predict bile velocity, flux and pressure across the liver lobule for individual patients from all four groups.

Figure 5

Individual-based model prediction of bile pressure p and bile fluid flux profiles based on measured bile canalicular geometries.

a, Abstraction of liver lobule by cylinder symmetry with radial coordinate ρ. The mechanistic model considers secretion of osmolytes (green) and osmotic water influx (blue) in a porous medium with ρ-dependent properties (see supplemental model description). Darcy’s law is assumed with a proportionality constant K(ρ) depending on viscosity μ, tortuosity τ, bile canalicular volume fraction εBC, bile canalicular radius rBC. All geometric parameters have been measured per patient. b-e, Model prediction for bile fluid pressure (solid line, left axis) and fluid flux (defined as average velocity times 2πρ) (dashed line, right axis) profiles for individual patients (colour) and for disease groups. f, Scatter plot of measured GGT levels versus predicted pericentral (region 1) bile fluid pressure from individual patients from all groups reveals a statistically significant positive correlation. One-sided hypothesis test. P-values and Spearman correlation coefficient are indicated in the plot. NC = 6 samples, HO = 4 samples, STEA = 8 samples, eNASH = 7 samples.

Extended Data Fig. 6

Estimates for a fraction of free lumen in total volume of a bile canaliculus.

a, Representative images of bile canaliculi for NC and eNASH liver tissue samples, used for making the estimates. Microvilli are well preserved. A red dashed line indicates lumen of a bile canaliculus. TJ, tight junction. NC = 3 samples, HO = 3 samples, STEA = 3 samples, eNASH = 3 samples. Scalebar, 500 nm. b, Estimation of fraction of free lumen by stereological point counting (the Cavalieri estimator). For each set of samples and each region (central / portal vein) a minimum of five EM images was used. NC = 3 samples, HO = 3 samples, STEA = 3 samples, eNASH = 3 samples, median ± MAD.

The model predicts similar bile velocities (5-10 μm/sec) and flux (10-35 μm2/sec) for the different groups in the periportal area (Fig. 5b-e). Comparable velocities have been reported in mouse[12]. However, the predicted pressure in the pericentral area differed significantly between the patient groups. In the NC and HO groups this pressure was relatively close to the periportal pressure (1508 ± 243 Pa) (Fig. 5b-c). In the STEA and eNASH groups, the model predicted a significant increase of bile pressure towards the pericentral region (2501 ± 1197 Pa). None of the pericentral pressures in the NC and HO patients exceeded 2000 Pa, but this was the case for 8 patients (53%) in the STEA and eNASH groups. In four of the STEA and eNASH patients (27%) the pericentral pressure exceeded 3000 Pa and in two patients (13%) even 4000 Pa (Fig. 5d-e). Therefore, our model predicts an increase in pericentral bile pressure in STEA and eNASH conditions, spanning different levels of severity, depending on the BC geometry and topology of individual patients.

We next set out to validate the model predictions. As it is impossible to measure bile flow and pressure at this resolution in the human liver, we considered possible consequences of changes in bile pressure. Increased bile pressure is a hallmark of cholestasis[40,41]. Therefore, as readout of increased bile pressure, we analysed the most commonly used cholestatic and liver damage biomarkers in serum, including bilirubin, gamma glutamyl transpeptidase (GGT), alkaline phosphatase (AP), total and primary bile acids (BAs), aspartate aminotransferase (AST) and alanine aminotransferase (ALT). To increase the statistical power, we analysed additional sera samples for the different groups (NC = 22, HO = 27, STEA = 31 and eNASH = 24 samples) (Supplementary Table 1). Whereas GGT was elevated in STEA and eNASH, and primary BAs only in eNASH, we did not detect significant changes in the levels of bilirubin, AP and total BAs between the groups (Extended Data 7).

Extended Data Fig. 7

Profile of serum cholestatic and liver injury biomarkers as well as bile acids during disease progression.

The levels of bilirubin (a), GGT (b), AP (c), AST (d), ALT (e), BA precursors (cholesterol, 7α-hydroxycholesterol and 27-hydroxycholesterol) (f), individual (CA, CDCA) and total primary BAs (g), individual (DCA, LCA, UDCA) and total secondary BAs (h), total BAs (i) and ratio secondary to primary BAs (j) were measured in the serum of the patients and represented in box-plots (median values as red lines, 25th and 75th percentiles as blue bottom and top edges of the boxes, extreme data points by whiskers). NC = 22 samples, HO = 27 samples, STEA = 31 samples, eNASH = 24 samples. One-tailed hypothesis test. *p-values < 0.05, **p-values < 0.01, ***p-values < 0.001.

Finally, we analysed the correlation between the predicted pericentral bile pressure and the biomarkers for individual patients (Extended Data 8a-f). Strikingly, we found a strong correlation for GGT and a significant correlation for AP (Fig. 5f and Extended Data 8b). There was a significant correlation also for primary BAs (however with one outlier h7252). AST and ALT, biomarkers of hepatocellular liver damage, also showed elevated levels in STEA and eNASH (Extended Data 7d-e) and a significant correlation to bile pressure (Extended Data 8e-f). This is very likely due to the reported correlation between AST and ALT levels, and fat content[16]. Altogether, our model predicts a significant degree of pericentral micro-cholestasis as a new component of the NAFLD pathophysiology.

Extended Data Fig. 8

Scatter plots and regression analysis of measured liver biomarkers and bile acids.

bilirubin (a), AP (b), total BAs (c), primary BAs (d), AST (e), ALT (f) and ratio secondary to primary BAs (g) measured in the serum versus the model-derived pericentral pressure in individual patients from all groups. Arrow indicates an outlier for primary BAs (h7252). NC = 6 samples, HO = 4 samples, STEA = 8 samples, eNASH = 7 samples. P-values and Spearman correlation coefficient are indicated in the plot.

Discussion

High definition medicine provides a novel approach to understand human health of individuals with unprecedented precision[1]. One of its pillars is the combination of image analysis and computational modelling to uncover tissue alterations at different structural and functional levels. During the last years, there has been an urge to gain a better understanding of NAFLD establishment and progression due to its growing impact on public health[42]. A lot of attention has been mostly drawn to signalling pathways, microbiome, metabolism, genetic risk factors, BAs, etc[13,42]. However, a major challenge is to understand how the molecular alterations detected are expression of the organ dysfunction, manifested as morphological and functional alterations of cells and tissue architecture. The classical histological analysis has provided insights into fundamental aspects of NAFLD. However, a quantitative description of the 3D tissue morphology is indispensable, particularly for the liver. Here, we used high resolution multiphoton microscopy and 3D digital reconstructions to generate a comparative dataset of structural changes of human liver tissue from NC, HO, STEA and eNASH patients. We identified a set of spatially distributed morphological alterations such as a characteristic size distribution of LD, nuclear texture homogeneity and BC morphology, that could be used as tissue biomarkers to resolve different stages of NAFLD progression. Although several morphological defects could be inferred from 2D images, topological characterisation of the BC and sinusoidal networks can only be extracted from a 3D reconstruction. Indeed, our 3D digital reconstruction of human liver tissue provided the first evidence that BC integrity is disrupted in NAFLD, bringing BC integrity and the mechanisms involved in its maintenance and homeostasis (polarized transport, bile flow, BAs turnover, etc.) into focus for future studies. Quantitative 3D features of BC architecture are therefore candidate parameters for unbiased classification of tissue samples.

The 3D spatially-resolved quantitative analysis of human liver samples revealed a set of unknown morphological features, ranging from the sub-cellular to tissue level, that seem to be perturbed during NAFLD progression. First, we detected changes in nuclear texture in pericentral hepatocytes, which have been reported in several diseases[20,43], and may reflect changes in transcriptional activity[21,22]. Second, although LD accumulation is a characteristic feature of NAFLD, our analysis revealed quantitative changes in their size distribution, with the medium-sized LD enriched in the periportal region in eNASH. In healthy conditions, the LD number and size are accurately regulated[44] and changes in their distribution point at cell-specific alterations in the mechanisms regulating LD biogenesis and catabolism. Third, and most striking, we observed alterations of the apical plasma membrane of hepatocytes and of the BC network. The large pericentral hepatocytes showed mislocalization of DPPIV, pointing towards a dysregulation in apical protein trafficking[33]. Interestingly, not all apical proteins were missorted, suggesting that trafficking defects could be pathway- (transcytosis) and/or cargo-specific. Transcytosis defects may be due to physical constraints caused by the large LD in the swollen hepatocytes. In addition, the uneven size of hepatocytes in the pericentral zone may introduce perturbations in cell packing leading to disruption of BC connectivity. The unaltered architecture of the sinusoidal network makes it unlikely that the reduction in BC connectivity is due to the sectioning of the enlarged hepatocytes within the thickness of the tissue slice, since this should have then affected both BC and sinusoidal networks in a similar fashion.

The weaker BC network connectivity and dilated BC radius are partly compensatory and counteracting effects. Whereas the first hinders bile flow, the latter eases it. Therefore, without a mathematical model it would have been impossible to predict the combined effect. Based on the geometrical and topological information extracted from the BC, our computational personalised model predicted high bile pressure in the pericentral area in STEA and eNASH patients. This suggests that STEA and eNASH livers are affected by a pericentral micro-cholestasis, a prediction that was supported by the detection of cholestatic biomarkers in serum. Cholestasis is reflected by a range of biomarkers. The strong correlation of systemic GGT levels to pericentral biliary pressure may reflect localized apical membrane stress, that does not lead to a bilirubin excretion problem at the organ level. Therefore, GGT may be an indicator of more advanced tissue micro-alterations that lead to bile pressure increase.

In recent years, a lot of research has been devoted to the role of BAs and BAs receptor FXR, in NAFLD[45-47]. However, there is currently no explanation for the alterations in BAs composition in blood, the decreased ratio of secondary/primary BAs observed by us (Extended Data 7f-j) and others[45,47,48], and whether it correlates with changes in tissue morphology[47]. Our data shed new light on this problem. The altered BC microanatomy and consequent pericentral micro-cholestasis may hamper BAs secretion into BC, as apical pumps (BSEP, MRP2) have to operate against elevated luminal BAs concentrations. This would lead to back-flux of primary BAs into the blood, reducing their availability for conversion into secondary BAs by the intestinal microbiota (Extended Data 7j and 8g), thus contributing to the changes in BAs composition observed in NAFLD[45,47,48]. Bile flow is essential for normal liver function. Bile accumulation, due to its detergent-like properties, can cause liver damage[49,50] and bile pressure can affect metabolism[51]. Indeed, the accumulation of LD[52] and BAs[53] could induce oxidative stress and trigger apoptosis[29], which is consistent with the reduction in pericentral hepatocytes observed in STEA and eNASH. The occurrence of pericentral micro-cholestasis is a new piece in the NAFLD pathophysiology puzzle that contributes to clarify some aspects of the disease so far without explanation, e.g. increase of GGT levels[54], BAs in serum[55], upregulation of MRP3 in NASH[56] and the beneficial effect of UDCA treatment in NAFLD[57], all signs of ongoing cholestasis[41,58].

The combination of experimental data with computational models of tissues has proven successful in elucidating pathogenetic mechanisms using animal models[12,59]. However, animal models very often fail to mimic human diseases, including NAFLD[13]. In this study, the geometrical models of liver tissue from human biopsies combined with spatially heterogeneous computational simulations revealed new aspects of NAFLD pathology. A firm separation of STEA from NASH will require further integration between the molecular analysis and 3D morphological tissue parameters. This approach may help to identify new biomarkers for early disease diagnosis and predict the functional status of the tissue with potential applications in high-definition medicine[1,60].

Methods

Human liver samples

Liver samples were obtained intraoperatively in patients in whom an intraoperative liver biopsy was indicated on clinical grounds such as exclusion of liver malignancy during major oncologic surgery or assessment of liver histology during bariatric surgery. Standardized histopathologic assessment was performed by a single pathologist (C. R.) based on the NAFLD activity score, the NASH-CRN[61]. The samples were divided into 4 groups: normal control (NC), healthy obese (HO), steatosis (STEA) and early NASH (eNASH) based on the following criteria: NC samples showed steatosis ≤5%, no inflammation, no ballooning and no fibrosis and were obtained during liver resections or biopsy in non-hepatobiliary malignancy in patients with a BMI<30kg/m2. Livers were either free of metastases or a distance of at least 2 cm to the next metastasis observed. Samples for the HO, STEA and eNASH categories were obtained in patients undergoing bariatric surgery with a BMI > 30 and subgroups defined by liver histology as follows: HO samples showed a normal liver histology such as the NC samples. STEA samples had >5% fat, no inflammation, fibrosis ≤1 and ballooning ≤1. eNASH samples were characterized by >5% fat, an inflammation grade of at least 1, fibrosis ≤1 and ballooning ≤1. All patients were of white Caucasian descent. None of the individuals underwent preoperative chemotherapy and liver histology demonstrated absence of both cirrhosis and malignancy. Biopsy specimens were fixed immediately in 4% paraformaldehyde 0.1%Tween-20/PBS for 2-3 days. Patients with evidence of viral hepatitis, hemochromatosis, or alcohol consumption greater than 20 g/day for women and 30 g/day for men were excluded. All patients provided written informed consent. The study protocol was approved by the institutional review board (Ethikkommission der Medizinischen Fakultät der Universität Kiel, D425/07, A111/99) before study commencement.

Immunolabeling, optical clearing and imaging

For the immunostainings where heat-induced epitope retrieval was required, liver slices were transferred to an eppendorf tube with citrate buffer pH 6.0 (Sigma-Aldrich Cat# C9999) and heated for 30 min at 80ºC. Next, immunolabeling (Supplementary Table 02) and optical clearing were performed as described previously[8].

Liver samples were imaged (0.3 μm voxel size) in an inverted multiphoton laser-scanning microscope (Zeiss LSM 780 NLO) using a 63x 1.3 numerical aperture glycerol immersion objective (Zeiss). DAPI and phalloidin A647 were excited at 790nm using a Chameleon Ti-Sapphire 2-photon laser and detected with non-descanned detectors (NDD). Alexa Fluor 488, 555 and 594 were excited with 488, 561 and 594 laser lines and detected with Gallium arsenide phosphide (GaAsp) detectors.

3D tissue reconstruction

The different components of live tissue (i.e. BC and sinusoidal networks, nuclei, lipid droplets and hepatocytes) were reconstructed from high-resolution (voxel size: 0.3 x 0.3 x 0.3 μm) fluorescent image stacks (~100 μm depth) of fixed liver tissue stained by specific antibodies and/or small fluorescent molecules (i.e. CD13, fibronectin, DAPI, BODIPY, and LDLR). To cover entire CV-PV axes, tiles of 2x2 or 3x1 image stacks were stitched using the Image Stitching plug-in of Fiji[62]. In general, one CV-PV axis was imaged per patient. Variability among lobules and within patients was tested by comparing the geometric and topological BC parameters used for the bile flow model, from three distant CV-PV axes from different parts of the biopsies for 13 patients covering the four groups. Then, all images were reconstructed using the software MotionTracking (http://motiontracking.mpi-cbg.de) as described in[6]. Briefly, all stitched images were desoised using the Bayesian foreground/background discrimination (BFBD) de-noising algorithm proposed in[6]. Next, the tubular structures (BC and sinusoidal networks) as well as nuclei were segmented using Maximum entropy local thresholding algorithm. Artefacts generated by the segmentation (holes and tiny isolated object) were removed by standard morphological operations (opening/closing). The triangulation mesh of the segmented surfaces was generated by the cube marching algorithm and tuned using an active mesh approach. In the case of tubular structures, representations of the skeleton, refereed as central lines, of the networks were generated. Central lines are represented as 3D graphs. Finally, the shape of the cell surface (i.e. hepatocytes) was found using an active mash expansion from the reconstructed nuclei. For details, refer to[6].

Morphological spatial analysis

To account for the variability of different morphological parameters along the liver lobule, the CV-PV axis was computationally divided in 10 equidistance zone. The zones were defined in terms of the distances to the closest CV and PV, as follows:

For i = 1:N, where N = 10, is the number of zones, ⌊⊣⌋ is the floor function, d and d are the distances to the closest CV and PV respectively. The average value of the different morphological parameters was calculated in every region.

Nuclear vacuolation

To determine whether nuclei are vacuolated, we calculated the DAPI mean intensity in the middle of the nucleus (inside a sphere of 1.2 μm located at the centre of the nucleus) as well as in the inner surface of it (over a layer of thickness 0.6 μm). If the inner DAPI intensity was 3 times lower than the DAPI intensity in the surface, the nucleus was defined as vacuolated.

Nuclear texture

We analysed the nuclear texture based on the Haralick texture features[63] measured over images of DAPI. We examined the four features related to the texture homogeneity/heterogeneity of the image, i.e. homogeneity (Angular Second Moment), local homogeneity (Inverse Difference Moment), Contrast and Entropy. To avoid boundary artefacts, we analysed only the intensity of DAPI inside a cube located in the centre of the nucleus and of length equal to half of the nucleus radius. All features were extracted from the normalized grey-level co-occurrence matrix (256 grey-levels) averaged over the 13 directional co-occurrence matrices (in 3D) at distances from 1 to 3 pixels. Vacuolated nuclei were excluded from the analysis.

Network connectivity

Two approaches were used to quantify the connectivity of the tubular networks (i.e. BC and sinusoidal networks). The first one, also referred as ‘fraction of connected network’, is based on the Central lines of the network and defined as the ratio between the length of the largest connected graph of the network and the total length of the network. The second one is based on the segmented 3D image of network and calculates the connectivity density using the Euler characteristic of the network (maximum number of branches that can be removed before the network is separated in two parts). The former calculation was performed in Fiji using the plugin BoneJ[64].

Lobule radius

Liver slices were stained with DAPI and GS antibody to facilitate the visual discrimination between CV and PV. Then, the sections were imaged with a Plan-Apochromat 10x/0.45 M27 objective at 2μm pixel size. The distance between the CV to the closest portal tract were measured in each 2D image. The lobule radius is reported as the median of all the measurements for each 2D image.

Hepatocytes per lobule

In order to estimate the number of hepatocytes per lobule section, we used the 3D information of the hepatocyte’s density δ and an approximation of the lobule section V, assuming lobule to be hexagonal prism, as follows: Where N, is the number of hepatocytes per lobule section, δ is the numerical density of hepatocytes (# of hepatocytes per tissue volume), R is the radius of the lobule and H is the height of the lobule section. We used the estimated radius on the lobule (in 2D) for R, and H = 100μm.

Optimal set of parameters for cluster analysis

In order to find the most relevant parameters for the classification, we used a greedy approach to obtain an optimal subset (OS) of parameters as described in[65]. Briefly, the values of individual parameters were systematically added to OS based on the classifier predictive accuracy The multi-class classification was performed using an error-correcting output codes (ECOC) model[66] for SVM, and the predictive accuracy was measured using 10-fold cross-validation.

Quantification of DPPIV localization

100 μm liver slices were stained with DAPI, Phalloidin, BSEP and DPPIV antibody. 20μm under the tissue surface, single plane images were acquired with an LCI Plan-Neofluar 63x/1.3 objective, 0.3 μm pixel size and using the same microscope settings for all images. The cell borders and apical domains were segmented using a threshold-based segmentation algorithm[6]. Since, no specific staining for the basal membrane was present, the basal domains were segmented manually based on phalloidin staining and the morphology. To analyse changes in the localization of DIPPV, the mean intensity of DPPIV was measured at each domain (i.e. apical, basal and lateral) and in an area covering a radius of covering 125 μm of tissue around each vein.

Electron microscopy

For the ultrastructural analysis of bile canaliculi, we used 100-μm thick vibratome sections of human liver biopsy tissue originally fixed with 4% PFA for several days and then stored in PBS. Before embedding, vibratome sections were fixed with 1% glutaraldehyde in 200 mM HEPES at least overnight and then cut into small pieces. Next, tissue was post-fixed with 1% osmium tetroxide prepared in 1.5% potassium ferricyanide, for 1 hr at room temperature. After washing with water, tissue was incubated with 1% tannic acid dissolved in 100 mM HEPES, pH 7.0 for 20 min, followed by incubation with 1% disodium sulfate for 5 min, and then by several washes with water. After that tissue was incubated with 2% aqueous uranyl acetate for 2 hrs at room temperature and protected from light. A graded ethanol series was used for dehydration: 70-80-90-96%, each 10 min, followed by absolute ethanol, 4x 15 min. Tissue was progressively infiltrated with epon over 24 hrs. Finally, tissue pieces were flat embedded between two teflon-coated glass slides and heat polymerized overnight. Embedded tissue pieces were remounted for longitudinal sectioning. Periportal or pericentral regions were selected on 1-μm sections, stained with methylene blue-azur II and examined in a light microsope. Then 70-80-nm thin sections were cut. These were stained with 0.4% lead citrate for 1 min and imaged in the Tecnai T12 transmission electron microscope (ThermoFisher), equipped with an axial Tietz CCD camera (TVIPS). Images of bile canaliculi were taken at 6,800x maginification by systematic and random screening. Bile canaliculi with poorly preserved microvilli were not included in the analysis.

To estimate a fraction of free lumen (i.e., lumen not occupied by microvilli) in bile canaliculi, we applied stereological point counting (the Cavalieri estimator), using Fiji[67] software. A test grid of vertical and horizontal lines (area per point 5000 nm2) was laid over images. Cross points over total lumen profile area and over free lumen profile area were separately counted. For each set of samples and each region (central / portal vein) a minimum of five EM images was used, and counts obtained from individual images were summed up to obtain total counts. The analysis was performed in a blinded way.

Bile flow model in the hepatic lobule

The disease phenotype of zonal cell swelling entails weaker BC network connectivity but also dilated BC radii, potentially indicative of increased BC fluid pressure. These are partly compensatory and counteracting effects as the first hinders flow but the latter eases flow. Therefore, without the help of a mathematical model it would be difficult to predict the combined effect on liver function (here bile flow) by all the observed alterations of the BC micro-anatomy (which are only attainable through a 3D imaging approach). Our mechanistic model integrates 6 quantitative data sets on the (1) BC volume fraction (Fig. 4c), (2) BC radius r(ρ) (Fig. 4d), (3) intra-canalicular cross-section fraction (1-α(ρ)) occupied by microvilli, (4) fraction of connected BC (excluding dead-end branches that do not carry flow) (Fig. 4g), (5) canaliculi tortuosity τ, (6) apical surface density A(ρ) and then predicts (i) higher mechanical stress on pericentral hepatocytes through elevated BC fluid pressure, (ii) lower overall bile flux and (iii) slightly higher canalicular bile acid concentrations for patients with fatty liver disease. When multiple samples were available for a patient, the median of the simulated values was determined (Extended Data 5).

We consider a porous medium in a three-dimensional domain Ω with an osmolyte influx density and a fluid influx density leading to a lumped concentration profile for osmolytes throughout the bulk of the domain. The bulk velocity field of the fluid is given in most generality by Darcy’s law according to which velocity is proportional, with a resistance coefficient to the pressure gradient This approach is similar to previous models of blood circulation through the sinusoidal network[68,69]. At steady-state, this system is described by conservation of mass for the osmolytes as well as the fluid and the connection between pressure gradient and velocity as follows:

To first approximation, the liver lobule can be assumed to have cylindrical symmetry. Therefore, we take Ω as a cylindrical domain with radial symmetry and radius L, approximating a single prismatic liver lobule. Hence, we need to consider only the radial variable which will be denoted by ρ. Following [12], we consider the fluid influx j as proportional to the difference in osmotic pressure (RTc − RTc0) and hydrostatic pressure p. Here c0 parametrises the osmotic background pressure in the cytosol of surrounding hepatocytes. Osmolyte diffusion with diffusion constant D along BC can be neglected (setting D = 0) for uniform osmolyte influx density and parameters chosen here. Below, we will quantify the small contribution of molecular diffusion to transport on the length scale of paths through the BC network. Using the divergence in cylinder coordinates with the radial velocity component w and the apostrophe denoting the derivative with respect to ρ, we obtain the following equations: where A(ρ) denotes the spatial density of hepatocyte apical membrane as osmotically active surface and κ denotes the membrane permeability for water. Note that w is the bulk velocity from which we can calculate the canaliculi velocity as wBC = w/ε where ε is the bulk porosity.

These are coupled first-order ordinary differential equations (ODEs) with the mixed boundary conditions w(ρ0) = 0 and p(L) = p where ρ0 is the radius of the central vein and p is the hydrostatic pressure in the bile duct. Since some model parameters such as A(ρ) and K(ρ) will below be derived directly from the quantification of image data, these parameters are given as heterogeneous profiles and the ODEs therefore possess heterogeneous coefficients.

For the special case of constant parameters and c ≈ c0 cand the equations can be solved as follows: As we have nonconstant parameters we had to solve the equations numerically and used this approximate solution for verification of the numerical solver.

The equation for the concentration profile can be readily integrated as follows

The integral expression represents the cumulative osmolyte source flux inside radius ρ and will be abbreviated as G(ρ). Using this equation, we can eliminate the concentration term in the equation for w′ as follows: or rewritten

Together with the equation for p′, these are coupled first-order ordinary differential equations (ODEs) with the mixed boundary conditions w(ρ0) = 0 and p(L) = p where ρ0 is the radius of the central vein and p is the hydrostatic pressure in the bile duct. Since some model parameters such as A(ρ) and K(ρ) will below be derived directly from the quantification of image data, these parameters are given as heterogeneous profiles and the ODEs therefore possess heterogeneous coefficients.

We can nondimensionalize ρ as follows: which modifies the equations as follows (introducing and immediately dropping the bars): Note that there is no singularity at ρ = ρ0 even though w approaches zero and appears in the denominator because we have, after applying L’Hôspital’s rule: and hence a finite limit for ρ → ρ0.

The equations are solved numerically using the shooting method. For this method, we set the initial values to w(ρ0) = 0 (the proper condition) and p(ρ0) = p* with some arbitrarily chosen value p*. We solve the ODEs for w and p with an implicit integration scheme till ρ = 1. Next, we iteratively updated p* until p(1) = p. To avoid any division by w = 0 at ρ = ρ0 in the first integration step, we make use of the above equation for the limit ρ → ρ0. This resultes in an algebraic equation for w′(ρ0): which has the solution with

The equations were integrated using an implicit Euler method (for the freely available source code see https://github.com/MichaelKuecken/bileflow).

For the determination of K, we consider a network of tubes with radii r(ρ) (Fig. 4d) and tortuosity τ that cause a bulk porosity of ε(ρ). The porosity ε(ρ) is determined from patient data as the product of the BC volume fraction (Fig. 4c) and the BC connectivity profile (Fig. 4g), amounting to pruning of branches that do not relay the flow. Where BC connectivity was quantified as 0, we set it to 0.01 assuming that at least some remote connection may exist and thereby potentially underestimating the pressure profile. We take into account that a significant portion of the BC lumen is taken up by microvilli. From our EM measurements, we know the fraction of free to total lumen α(ρ), therefore the effective porosity and radius are given by εBC(ρ) = α(ρ) ε(ρ). To determine the effective BC radius we consider two extreme cases. In the first scenario all microvilli are considered to be pushed to the canaliculi walls by the flow and we have rBC(ρ) = reffr(ρ) where is approximately 0.4 to 0.6 here. In the other scenario all microvilli behave like a porous medium inside each BC and there is no free lumen in the middle of the canaliculi. In this case porous media theory[70] predicts a free radius of α(ρ)rvilli/(1 − α(ρ)) where rvilli∼6 · 10−8m is the radius of the microvilli[71]. In this case we obtain an effective radius of rBC(ρ) = reffr(ρ) where reff is approximately 0.03. In reality the actual reff will be between the two extreme cases. We use in our simulations reff = 0.2 which is slightly smaller than the corresponding value of 0.344 computed for the mouse model[12]. This is justified by the higher microvilli density observed in humans (Extended Data Fig. 5).

The tortuosity of a curve (τ) is defined as the ratio between its length and the distance between its endpoints. In our system there are two contributions to the tortuosity. At first, there is a contribution τBC from the fact that the canaliculi are not straight. This number is determined from individual patient data (τBC ranges from 1.5 to 2.3). The second contribution τM originates from the flow around the microvilli in the canaliculi. This contribution is estimated as τM = 1.6. The overall tortuosity is then determined as τ = τBCτM.

We also take into account that the geometry of the liver lobule is not simply cylindrical but that the portal bile ducts have approximately a hexagonal arrangement. Therefore, there is a convergent flow of bile towards the bile dcts in the periportal area[12]. We model this effect as an additional factor 1/f(ρ) in the determination of K(ρ) where f(ρ) is defined as follows

With these quantities K(ρ) is given[72] in terms of patient’s data as

We use the universal gas constant R = 8.314 J/mol/K and T = 293 K. The values for A(ρ) and L are of the order of 5x104x 1/m and 4x10-4 m, respectively, but the exact values are taken from our measurements for each individual patient sample.

The central vein radius was measured to be 2x10-5m in our samples. The osmolyte loading of hepatocytes is chosen similar to that of blood and bile, c=300 mmol/l. We determine the approximate osmolyte flux g as follows: about 500 ml of bile with an osmolarity of 300 mmol/l[39] drain every day from a human liver of 1 l in size. This means that 300 mmol/l*0.5l = 150 mmol of osmolytes are secreted through the apical membranes of hepatocytes in a fraction 1/6 of a day, resembling the modulation of bile turnover by dietary cycles. This leads to g = 150 mmol/10-3m3/(4*3600s) = 0.01 mol/m3/s. The viscosity of bile was chosen as μ = 9.2x10-4 Pa s which is close to the viscosity of water. The membrane permeability κ = p/(16*μ) = 6.97*10-12 m/(16*μ) = 4.7x10-10 m/(Pa s) has been calculated from Table S2 in[12]. We consider an outer pressure p = 1000 Pa (estimated according to measurements reported in[73,74]. All these parameter values are assumed to be identical across all samples whereas disease could affect them as shown by three-fold higher median bile viscosity in patients with cholesterol gallstones[75]. Our choices of identically low viscosity and identically low portal pressure boundary value p represent a conservative limit for the disease cases and together with identical bile volume and cytoplasmic osmolarity let us focus on the important role of the altered micro-anatomy of BC.

Therewith all parameter values have been determined from previously reported data and geometrical measurements in this work such that no parameter fitting is necessary. Using the shooting method described above, we predict the spatial profiles at steady state for bile pressure and flux as shown in Fig. 5. The obtained solutions are qualitatively in agreement with the asymptotic states in numerical simulations of the time-dependent problem along a single tube which we studied using the software Morpheus[38,76]. In extreme cases, altered BC geometry can raise the fluid pressure to a point comparable with the osmotic driving and then trigger a so-called wet-tip to dry-tip transition where flow stalls near the closed pericentral branch ends[38,77]. In addition to the osmotic driving, contractility of the apical cortex of hepatocytes adds temporally fluctuating pressure modulations that accelerate bile flow by 50% on average as observed in mouse[12]. Here, we calculated bile velocities purely based on the osmotic driving. Taking contractility into account for humans as well would approximately double the bile velocity prediction but leave the temporally averaged pressures unchanged. Further, we calculated the relative contribution of diffusion to the canalicular transport flux given the concentration profiles obtained under the assumption of dominating advective transport. We have found that diffusion only contributes a negligible additional flux of ~10-4 relative to the advective flux, verifying the self-consistency of the initial assumption. Given our accurate 3D reconstructions of the BC network in health and disease, future computational work can address the osmotically-driven flow problem in the actual spatially-resolved canalicular network, replacing the porous medium approach used here by Computational Fluid Dynamics (CFD) methods as recently applied to blood flow through the sinusoidal network[78,79]. However, note that for bile, the fluid sources are adaptive (feedback loops via osmolyte concentration and fluid pressure) and distributed over the entire surfaces of all BC lumens.

Bile acids and serum markers

Individual serum bile acids and their precursors were analysed either by gas-liquid chromatography mass spectrometry (GLC-MS) or liquid chromatography (HPLC) as described[80].

Statistical analysis

For the spatially-resolved quantifications (i.e. statistics along the CV-PV axis), the median ± MAD (Median absolute deviation) values per region are shown.

Global quantifications (i.e liver enzymes, bile acids and overall tissue parameters) are shown using box-plots. Whereas median values were shown as red lines, the 25th and 75th percentiles were represented by the blue bottom and top edges of the boxes, respectively. The whiskers were extended up to the most extreme data points that are not considered outliers. The statistical significance analysis was performed using: 1) Kruskal-Wallis (non-parametric method for testing whether samples originate from the same distribution), one-way ANOVA on ranks, to compare the data within conditions. 2) Wilcoxon rank sum test to compare data between conditions (test all conditions against the normal control). The results of the correlation analysis are now reported as Spearman correlation coefficients with Wilcoxon rank sum test values. One-sided hypothesis testing was used if priori knowledge (previous studies) suggested one side hypothesis otherwise two-sided was applied, in both cases assuming unequal variances. Median, MAD, number of samples and p-values for all figures are provided in Supplementary Table 3. All statistical analysis was performed using MATLAB2018b.