Resilience of three-dimensional sinusoidal networks in liver tissue

Can three-dimensional, microvasculature networks still ensure blood supply if individual links fail? We address this question in the sinusoidal network, a plexus-like microvasculature network, which transports nutrient-rich blood to every hepatocyte in liver tissue, by building on recent advances in high-resolution imaging and digital reconstruction of adult mice liver tissue. We find that the topology of the three-dimensional sinusoidal network reflects its two design requirements of a space-filling network that connects all hepatocytes, while using shortest transport routes: sinusoidal networks are sub-graphs of the Delaunay graph of their set of branching points, and also contain the corresponding minimum spanning tree, both to good approximation. To overcome the spatial limitations of experimental samples and generate arbitrarily-sized networks, we developed a network generation algorithm that reproduces the statistical features of 0.3-mm-sized samples of sinusoidal networks, using multi-objective optimization for node degree and edge length distribution. Nematic order in these simulated networks implies anisotropic transport properties, characterized by an empirical linear relation between a nematic order parameter and the anisotropy of the permeability tensor. Under the assumption that all sinusoid tubes have a constant and equal flow resistance, we predict that the distribution of currents in the network is very inhomogeneous, with a small number of edges carrying a substantial part of the flow—a feature known for hierarchical networks, but unexpected for plexus-like networks. We quantify network resilience in terms of a permeability-at-risk, i.e., permeability as function of the fraction of removed edges. We find that sinusoidal networks are resilient to random removal of edges, but vulnerable to the removal of high-current edges. Our findings suggest the existence of a mechanism counteracting flow inhomogeneity to balance metabolic load on the liver.


Introduction
Leaf venation 1 , fungal mycelium 2, 3 , but also animal trails networks 4 , river deltas 5 , and even force networks in granular materials 6 , each represent natural transport networks formed by self-organization. Inside our body, the micro-vasculature forms a plexus-like, three-dimensional network of small capillaries that deliver nutrients to every cell in a tissue and removes waste products 7 .
Mathematically, transport networks are spatial networks, i.e., graphs embedded in space. It has been proposed that the presence of cycles in these graphs provides redundancy, and thereby resilience against the failure of individual links 8 . Yet, to the best of our knowledge, this concept has never been tested in three-dimensional, micro-vasculature networks. Past research addressed resilience properties almost exclusively in two-dimensional biological transport networks. It was shown that self-organization by local feedback rules can generate hierarchical networks resembling those of leaf networks, which optimize flow resistance upon removal of a single link 1 . For example, work on two-dimensional networks addressed the balance between the cost of network formation and network resilience to random failure 9 , or the cost of repair after perturbations 10 .
The restriction to two-dimensional networks in past research was largely a consequence of the considerable difficulties in imaging three-dimensional networks. Yet, topology suggests fundamental differences between two-dimensional and threedimensional spatial networks, as it imposes constraints on the distribution of cycles in the network 11 . Here, we take advantage of recent technological advances in high-resolution imaging of adult mouse liver tissue 12,13 , allowing us to study statistical geometry and resilience of three-dimensional sinusoidal networks.
The liver is the largest metabolic organ in the human body. It is responsible for storage of metabolites, secretion of digestion enzymes and detoxification of blood 7 . The liver is organized into millimeter-sized basic functional units termed lobules, see Figure. 1A. Each lobule comprises a central vein (CV) and a portal vein (PV) connected by a dense plexus, termed the sinusoidal network, which transports metabolite-rich blood from PV to CV. The sinusoidal network contacts each of the hepatocytes, the main parenchymal cell type in the liver, which take up nutrients and toxins from the blood. The liver serves as the central organ of blood detoxification, including the metabolism of medically relevant drugs. Hence, an ongoing intense research effort focuses on the pharmacological modeling of the uptake of drugs from the blood stream in the liver [14][15][16][17][18][19][20] . Understanding fluid flow in the sinusoidal network is imperative to improve these models and make functional predictions. Previous simulation studies of blood flow in liver microvasculature were either limited to large veins 16 , or considered only small spatial regions 14,19 . For the entire lobule, spanning from PV to CV, only continuum models had been used to simulate flow 14,21,22 . Finally, the sinusoidal network can become damaged upon non-lethal toxification, prompting the question of resilience properties of this network. How the permeability of the sinusoidal network responds to perturbations is not known.
Here, we analyze network geometry using a digital reconstruction of the sinusoidal network based on high-resolution image data of adult mouse liver 12,13 . We develop a network generation algorithm that reproduces statistical features of the sinusoidal network (node degree distribution, edge length distribution, mean nematic order parameter), enabling us to simulate arbitrarily sized networks from spatially restricted biological samples and, moreover, to explore in silico a design space of three-dimensional networks. While simulating random graphs with given degree distribution is a classical problem of combinatorics 23 , and popular software packages exist for common models of random spatial networks 24 , we were not aware of previous network generation algorithms that allow to prescribe both degree and edge length distribution.
Sinusoidal networks display a weak nematic alignment along the direction of flow 13,14,25 . Using our algorithm, we can systematically vary this nematic alignment in simulated networks. We empirically find a linear relationship between the anisotropic permeability of simulated networks and a nematic order parameter of the networks that quantifies their anisotropic geometry. Permeabilities allow to efficiently compute macroscopic, tissue-level flows using a continuum model 14,21,22,26 , thus providing an effective medium theory of fluid transport.
To quantify the fault tolerance of these networks, we introduce a new resilience measure, which we term permeabilityat-risk and which quantifies changes in network permeability if a given fraction of network links is removed. The resulting permeability-at-risk curves can be considered as a generalization of the bond percolation problem in the theory of random resistor networks 27,28 . We find that simulated networks with weak nematic order display a substantially increased permeability along the direction of nematic alignment. If the mean nematic order parameter equals that of sinusoidal networks, this increased permeability does not compromise network resilience as compared to isotropic simulated networks. Our minimal transport model, which assumes constant and equal flow resistance per unit length for each edge, predicts that the distribution of computed currents is very inhomogeneous in the network, with a few edges carrying most of the current. This renders these networks highly vulnerable to the removal of high-current edges, despite their resilience against random removal of edges. In the discussion, we speculate on mechanisms such as shear-dependent adaptation of the diameter of sinusoids [29][30][31] , or transient clogging by erythrocytes 32,33 , which would both affect especially high-current edges, and could homogenize the time-averaged distribution of currents in the network, thereby reducing the vulnerability of sinusoidal networks to the removal of high-current edges.

Experimental Data and Network Metrics
To analyze the statistical geometry of three-dimensional microvasculature networks, we took advantage of advances in highresolution imaging of murine liver tissue 12,13 . Based on segmented three-dimensional image data the skeleton of the hepatic sinusoidal network was computed using MotionTracking image analysis software 12, 34 , see Fig. 1B.
Next, a cleaned version of the raw network data was computed: (i) small disconnected network components not connected to the largest component were discarded, (ii) connected nodes separated by a distance smaller than a cut-off distance R c = 8 µm (outer diameter of sinusoids) were merged into a single node, (iii) in a subsequent pruning step, dead ends were removed, leaving only nodes with node degree d ≥ 2. Finally, linear-chain motifs consisting of degree-two nodes in series were replaced by a single link with weight equal to the total length of the linear chain. In rare instances, removal of a linear chain might yield triangles at the extremities of the network, which were also removed. The remaining node points are exactly the branch points of the biological network, whose positions are determined with high precision. This clean-up procedure reduces ambiguity on small network details that were difficult to resolve with current imaging techniques. It provides a faithful representation of the sinusoidal network used in all subsequent analysis, see Fig. 1C.
The skeleton of the sinusoidal network is a spatial network with set of node points P = {q i } and set of edges E that connect pairs of points. We find that the sinusoidal network is a homogeneous network with mean node degree d = 3.3 ± 0.6 and a unimodal distribution p(l) of edge lengths l with mean l = 17 ± 8 µm, see Fig. 1CD.
We characterize the relative position of node points q i in terms of their normalized radial distribution function

2/20
where ρ(r) = ∑ i δ (r − p i ) is the point density of nodes, ρ 0 = ρ the mean density (with units of an inverse volume), and r 0 the position of a 'central node'. The radial distribution function is closely related to the structure factor, which is used in condensed matter physics to describe the packing of particles and which can be considered the spatial power spectral density of ρ(r) 35 . For an ideal gas, g(r) = 1. Fig. 1F shows g(r) for the node points of the sinusoidal network. Interestingly, we observe a peak at a characteristic distance r ≈ 10 µm, indicating a characteristic distance between the branching points of this network. In fact, the resultant g(r) resembles the radial distribution function of a fluid, with short-range repulsion of fluid particles. Interestingly, the observed characteristic distance of branch points corresponds to approximately half the diameter of hepatocytes 13 , which is likely to set a characteristic mesh-size of the sinusoidal network. The ratio of branch points to number of hepatocytes is approximately 2 : 1 (i.e., n = 1643 branch points and 857 hepatocytes for sample volume shown in Fig. 1C).
The sinusoidal network is a sub-graph of its Delaunay graph Given the set P of node positions of the sinusoidal network, we construct the corresponding Delaunay graph with set of edges D(P). The Delaunay graph generalizes the familiar Delaunay triangulation in the plane to three space dimensions, and may be considered as the graph that connects 'nearest neighbors' in P. Specifically, we may assign to each node q ∈ P a neighborhood V q defined as the set of all points that are closer to q than to any other point of P \ q. Each V q is a polyhedron. Together, these polyhedra tesselate three-dimensional space, defining the so-called Voronoi tesselation of P. Now, an edge connecting nodes q 1 , q 2 ∈ P belongs to the Delaunay graph D(P) if and only if the corresponding polyhedra V q 1 and V q 2 share a common face. Correspondingly, D(P) is the dual graph of the Voronoi tesselation. The Delaunay graph exists and is unique provided the node points are in general position 36 .
Remarkably, we find that the sinusoidal network is a subgraph of the Delaunay graph to very good approximation: 99% of the edges in E are contained in D(P), see Fig. 1G. The 1% of edges contained in E but not D(P) (marked red) are longer than average (with a mean length of 57 ± 19 µm), yet contribute less than 10% to the network permeabilities computed below. Thus, the sinusoidal network can be considered a network of nearest-neighbor edges, reflecting the design requirement of a space-filling network that connects all hepatocytes in the tissue.

The sinusoidal network contains the minimum spanning tree
The set of node points P defines also a second graph, the minimum spanning tree with edges M(P). The minimum spanning tree is defined as the connected and cycle-free graph with node points P for which the sum of the edge lengths is minimal. Note that the minimum spanning tree M(P) is always a sub-graph of the Delaunay graph D(P) for any set P of points in Euclidean space. Remarkably, the sinusoidal network contains the minimum spanning tree of its node points as a sub-graph to good approximation: 90% of the edges M(P) of the minimum spanning tree belong also to E, see Fig. 1H. This finding suggests a possible optimization of the sinusoid network for shortest paths. The liver is organized into millimeter-sized lobules, each containing a central vein (CV, cyan) and several portal veins (PV, orange). Blood flows from PV to CV through the plexus-like sinusoidal network. B. Raw sinusoidal network obtained by high-resolution 3D imaging of mouse liver tissue. Shown is the network skeleton with branch points (black) and edges (blue), defining a spatial network. CV and PV shown schematically. Scale bar: 100 µm. C. Sinusoidal network N obtained from the raw network by removal of degree-2 nodes and disconnected components, used in all subsequent analysis. D, E, F. Distribution of node degree p(k), distribution of edge lengths p(l), and radial distribution function g(r) of spatial node positions for the sinusoidal network N , respectively. G. The majority of edges of the sinusoidal network are contained in the Delaunay graph corresponding to its set of node positions. Edges of the sinusoidal network not contained in the Delaunay graph are shown in red. Inset: cartoon illustrating a Delaunay graph (black) and Voronoi tessellation (gray) of a set of points. H. The sinusoidal network contains the majority of edges of the minimum spanning tree M(P) of its set of node positions P. Edges of M(P) not in E are shown in red. Inset: cartoon representation of a minimum spanning tree (green).

A network generation algorithm for spatial networks
We developed a Monte-Carlo algorithm to generate synthetic networks that faithfully reproduce the statistical features of spatially restricted samples of hepatic sinusoidal networks, using multi-objective optimization of both node degree and edge length distribution, see Fig. 2A.
Our algorithm proceeds in two steps: first, we simulate surrogate node positions P sim , followed by a second step of selecting a set of edges E sim connecting these nodes. Importantly, our analysis of sinusoidal network graphs allows us to restrict to simulated networks that are subgraphs of the Delaunay graph D(P sim ) and at the same time contain the minimum spanning tree M(P sim ), M(P sim ) ⊆ E sim ⊆ D(P sim ). This restriction dramatically reduces the search space without which network optimization would be computationally unfeasible.
In the first step of the algorithm, random node points are generated by simulating random packings of equally sized hard spheres. This minimal model comprises only two fit parameters (the volume fraction of spheres and their radius), and reproduces the characteristic feature of short-range repulsion in the radial distribution function g(r) for the center positions of the spheres, see Fig. 2B. The final set of node positions P sim is determined by taking a random selection of about 20% of simulated hard spheres centers, to match the measured density ρ 0 of node positions in sinusoidal networks. This random selection does not change g(r), as it scales down both denominator ρ 0 and numerator dx ρ in Eq. (1) by the same factor. This random selection reflects volume exclusion by hepatocyte cells (which comprise a volume fraction of about 80% in liver tissue 7 ) and other cell types.
In the second step, we select a subgraph of the Delaunay graph D(P sim ) of the set of surrogate node positions P sim , using multi-objective optimization for node degree and edge length distribution. We introduce the cumulative probability density functions (CDF) CDF(d 0 ) = ∑ d≤d 0 p(d) and CDF(l 0 ) = l 0 0 dl p(l) for node degree and edge length distribution of the sinusoidal network, respectively, as well as analogous CDFs CDF sim (d) and CDF sim (l) for simulated networks. Note that p(l) has dimensions of an inverse length, hence CDF(l) is dimensionless. We define two cost functions C d and C e for the node degree and edge length distribution, respectively, as the difference of the CDFs and We normalized C e , using the cut-off distance R c as a characteristic length scale. As illustration, Fig. 2C,D shows probability distribution functions (PDF) and cumulative distributive functions (CDF) for the simulated network from panel A. The difference in CDFs provides a robust distance measure, which is in particular robust against small shifts of the probability distributions p(d) and p(l). We also introduce a multi-objective cost function as weighted mean of the individual cost functions, parametrized by a weight α We determined sets of edges E * sim (α) that minimize C α for given α by simulated annealing, see methods section for details. Note that choosing either α = 0 or α = 1 would correspond to optimizing only a single objective, C e or C d , respectively. These single-objective optima set the utopia points, i.e., the best-achievable values for each individual cost function, Fig. 2E. The value of the respective other cost function define the nadir points N e and N d . The definition of the nadir points requires to take a limit of α, α 1 or α 0, as Otherwise, the values of C e and C d are not well-defined for a single-objective optimization with α = 1 or α = 0, respectively.
In the general case of multi-objective optimization with 0 ≤ α ≤ 1, the values of the individual cost functions for the optimal network are bounded by the utopia and nadir points, , defines a Pareto front that separates a region of impossible pairs of (C d ,C e ) values, from a region of possible values. Remarkably, this Pareto front deviates only little from the straight lines defined by C d = U d and C e = U e for intermediate values of α. Hence multi-objective optimization for both degree and edge length distribution can be achieved with minimal impairment on the individual cost functions. This shows that our algorithm is robust with respect to the particular choice of α.
In the following, we use the value α * = (N e − U e )/(N d − U d + N e − U e ), which corresponds to choosing equal weights for the individual cost functions after an appropriate normalization, see Methods section.

5/20
Objective 1 (node degree): Comparison of radial distribution function g(r) from sinusoidal network (black dots), and simulated node positions (red). C, C'. Probability density function (PDF) and cumulative distribution function (CDF) for the node degree distribution p(k) of a simulated network (red), and the sinusoidal network (blue). D, D'. Same as panel C, but for the edge length distribution p(l). E.
Pareto-front of multi-objective optimization. Shown are normalized costs for each of the two objectives: optimization of node-degree and optimization of edge-length distributions, for varying weight α of the multi-objective cost function.

6/20
Nematic order determines anisotropy of network permeability The sinusoidal network is not isotropic, but displays weak nematic alignment along a direction parallel to the direction of blood flow along the PV-CV axis, see Fig. 3A. Previous work has shown that the local director of sinusoidal network anisotropy follows curvilinear flow lines, where PV and CV serve as source and sink, respectively 13,25 . Inside a central region of the lobule, these flow lines are approximately straight, which simplifies our subsequent analysis. We quantify nematic alignment of the sinusoidal network inside this central region in terms of a nematic order parameter Here, e x is a unit vector pointing along the CV-PV axis (which is parallel to the x axis for the given data set), and e j is a unit vector parallel to the j-th edge. In Eq. (5), averaging is performed over all edges completely located inside the region of interest.
To compute S, we choose a cuboid region of interest of dimensions L x × L y × L z , located in the central region of the lobule, with cuboid edges aligned with the x, y, z axes of a coordinate system, see Fig. 3A. We found S = 0.12 ± 0.06 (mean±s.e., n = 3 independent data sets from different mice; S = 0.18 for the data set in Fig. 3A) 1 . Next, we modified our stochastic network generation algorithm, to generate networks that likewise show nematic alignment along a specified axis. For that aim, we compute a nematic order parameter S sim for simulated networks analogous to Eq. (5) and use an augmented cost function in the optimization with interaction parameter λ . Here, is a normalization constant introduced for numerical convenience, corresponding to the use of a normalized cost function C, see Eq. (10). For simplicity, we do not account for a possible dependence of edge lengths on their orientation in space. Fig. 3B shows an example of a simulated anisotropic network (with nematic order parameter S sim ≈ 0.16, similar to the value S ≈ 0.18 measured for the experimental sinusoidal network shown in Fig. 3A). More generally, the nematic order parameter S sim is a monotonic increasing function of the interaction parameter λ in Eq. (6), see Fig. 3C. The nematic order parameter S sim of simulated networks saturates at a maximal value S sim = 0.16 ± 0.01, suggesting that there exists a maximal value of S compatible with the given degree and edge length distribution. The geometric anisotropy of the sinusoidal network results in anisotropic transport properties. We can quantify the transport properties of spatial networks using a minimal flow model that assumes a constant resistance κ per unit length for each edge. This assumption is justified since Reynolds numbers are small 14 . Thus, we can model flow using Kirchhoff's laws (p j − p k ) = κ l j,k J j,k for all edges connecting nodes j and k (Ohm's law) , 0 = ∑ k J j,k for all nodes j (conservation of current).
Here, J j,k denotes the signed current through the directed edge connecting node j and k (with units m 3 s −1 ), l j,k the length of that edge, and κ a constant resistance per unit length (with unit Pa m −4 s −1 ), corresponding to the simplifying assumption of a constant diameter of sinusoid tubes. Fig. 3D shows the flow computed for the chosen region of interest of the sinusoidal network. Here, we impose a pressure difference ∆p at opposing boundaries of the region of interest as indicated. Specifically, we identify all nodes close to the two boundary surfaces of the cuboid region normal to the x axis as either sink and source nodes, with p i = p 0 and p i = p 0 + ∆p, respectively. Solving the linear problem defined by these boundary conditions and Eq. (7), with conservation of current at all nodes that are neither source nor sink, yields the pressure at these remaining nodes and the currents J j,k . The total inflow J x at the source nodes equals the total outflow at the sink nodes. Contrary to intuition, we find a distribution of flows that is very inhomogeneous, see Fig. 3D. As an interesting side remark, also an inhomogeneous distribution of flows can allow for a homogeneous supply of the tissue, see Fig. S3 for results from a minimal model of metabolic uptake from the sinusoidal network by surrounding tissue.
The permeability of a network in x direction is proportional to the ratio of the total current J x divided by the imposed pressure difference ∆p. We normalize the permeability by the dimension L x of the cuboid in x direction, its cross-sectional area A x = L y L z , and the resistance per unit length κ of individual edges, to obtain a normalized permeability (with units of an area density) as Analogously, we define K yy and K zz for the y and z direction, respectively. This definition of a normalized permeability defines a purely geometrical measure of the network that is independent of κ, a property known as Darcy's law. We confirmed that indeed Darcy's law holds approximately for the large samples used, i.e., the normalized permeability is indeed independent of the dimensions of the cuboid if boundary effects can be neglected, see Fig. S1. For smaller sample volumes, computed permeabilities display a high statistical variability, but still have the same mean value. To convey the geometric meaning of the normalized permeability, we consider a minimal network that consists of straight lines parallel to the x axis, running from one boundary face of the cuboid to the opposing face: there, we have K xx = K 0 , where K 0 is the area density of these lines in any cross section of the cuboid. Fig. 3E shows computed normalized permeabilities K xx and K yy of sinusoidal networks in the direction of nematic alignment and normal to it, respectively. We have K xx = 1.1 ± 0.1 10 3 mm −2 and K yy = 0.6 ± 0.1 10 3 mm −2 (mean±s.e., n = 3). The computed permeability in z direction, K zz = 0.9 ± 0.3 10 3 mm −2 , is rather unreliable due to the small dimension L z of the experimental sample in z direction. The anisotropy ratio K xx /K yy ≈ 1.83 is consistent with a previously reported value 2.2 computed for a 0.15 µm × 0.15 µm × 0.15 µm sample of the sinusoidal network imaged using µCT 14 .
For the simulated networks, we find that the permeability K xx in the direction of nematic alignment increases linearly as a function of the nematic order parameter S sim , while the permeability K yy in normal direction decreases, see Fig. 3E.

Permeability-at-risk
We used the minimal flow model to study the resilience properties of sinusoidal networks, i.e. the ability of the network to support flow even after partial damage. If edges are removed from a network, the permeability of the network decreases. 2 We quantify network resilience in terms of permeability-at-risk curves, see Fig. 4A. Specifically, we plot the permeability of perturbed networks as a function of the fraction γ of removed edges. At a critical value γ * of γ, this permeability becomes zero, indicating the percolation threshold of the networks. We consider two different scenarios for the removal of edges: (i) removal of high-current edges, i.e., those edges that carried the highest current in the unperturbed network, (ii) random removal of edges with probability proportional to their length. We find that removing high-current edges results in a much steeper decrease of the permeability (scenario i, solid line), as compared to a random removal of edges (scenario ii, dashed). This shows that the edges that carry a high current in the unperturbed network are indeed crucial for the global transport properties of the network. Transport by these high-current edges is only partially compensated by re-routing of flow through low-current edges. In the absence of additional perturbations, low-current edges are dispensable and the network permeability decreases only little if the majority of low-current edges is removed (e.g. if all edges with current below the 25% percentile among those edges that do not touch the boundary of the region of interest are removed, K xx decreases by 11.1 ± 5.1%). However, a network without these low-current edges displays permeability-at-risk curves that decay even steeper as function of γ, see Fig. S2: low-current edges contribute at least partially to network resilience.
Remarkably, the computed permeability-at-risk curves match to very good extend in both scenarios for sinusoidal networks and simulated networks. Here, the nematic order parameter S of simulated networks matches the value of the experimental sinusoidal networks.
Even for a strong reduction of network permeability, the majority of network nodes can still be connected to both the source and the sink. For example, we find for γ = 10% that only 1.3 ± 2.1% and 0.5 ± 0.1% of nodes in the simulated network are disconnected from source or sink in scenario i and ii, respectively. As expected, these fractions increase to 100% if γ approaches the respective percolation thresholds.
Next, we asked how nematic order of a network influences it resilience, using our network generation algorithm that reproduces the statistical features of sinusoidal networks, while allowing us to tune its nematic order parameter. We find that nematic simulated networks (whose nematic order parameter matches those of the sinusoidal networks) have a higher permeability along the direction of nematic alignment than isotropic networks, not only in the absence of perturbations, but also if edges are removed, see Fig. 4BC. We anticipate that the observed weak nematic alignment of sinusoidal networks represents a good compromise between increasing network permeability and decreasing network resilience.
Nonetheless, in scenario i, where high-current edges are removed first, both sinusoidal networks and simulated networks display a comparatively low resilience. To understand the origin of this vulnerability, we consider simple example networks. In Fig. 4BC, we plot resilience curves of simple regular lattices: a full cubic lattice (blue) and a layered stack of square lattices (green). Both regular lattices display a much weaker reduction of network permeability upon edge removal. We attribute this apparent resilience of regular lattices to the fact that a large subset of their edges carry the same current in the unperturbed network, and thus can be considered equally important. In contrast, disordered networks display the same vulnerability: we tested random sub-lattices of a cubic lattice and observed a similar dramatic drop in network permeability if high-current edges are removed, see Fig. 4BC (cyan). The filling fraction of 48% of this sub-lattice served as fitting parameter that allowed to change the disorder of the network in a continuous manner. In conclusion, a low resilience against removal of high-current edges, as observed in sinusoidal networks, appears to be a general property of disordered networks, which are typically characterized by an inhomogeneous distribution of currents across their edges.  . Resilience of simulated three-dimensional transport networks. A. Permeability-at-risk. Shown is the permeability of both experimental and simulated networks after perturbation as a function of the fraction γ of removed edges for two basic scenarios: removal of edges that carried the highest current in the unperturbed network (solid), and random removal of edges with removal probability proportional to their length (dashed). Results are similar for sinusoidal networks (blue) and simulated networks (red). Above a critical γ * (percolation threshold), the permeability drops to zero. Permeabilities K xx were computed along the direction of nematic order (x direction), and normalized by the permeability K 0 of the unperturbed network to facilitate comparison (K 0 = 1089.1 ± 121 mm −2 for sinusoidal networks, K 0 = 943 ± 32 mm −2 for simulated networks). Shaded areas indicate mean±s.e. (n = 3 networks for each condition; region of interest for sinusoidal networks as in Fig. 3A). B. Logarithmic plot of permeability as function of γ for simulated networks for scenario i of removal of high-current edges for isotropic networks (black; λ = 0), and nematic networks (red; λ = 0.2, resulting in S ≈ 0.129, comparable to the value of sinusoidal networks). We compared the permeability-at-risk curve of simulated networks to that of three simple networks: full cubic lattice (blue; 5 × 5 × 5, lattice spacing 31.6 µm ensuring K xx (γ = 0) = 10 3 mm −2 ), layers of square lattices (green; 5 × (5 × 5), lattice spacing 31.6 µm), random sub-lattice of cubic lattice (cyan; sub-lattice of 10 × 10 × 10, with filling fraction 48%, mean of n = 10 realizations, lattice spacing 14.3 µm).
Cartoons of the three simple networks are shown below. C. Analogous to panel B for scenario ii of random removal of edges.

Discussion and Outlook
We analyzed geometry and transport properties of microvasculatory sinusoidal networks in three-dimensional liver tissue. We drew advantage of recent advances in three-dimensional imaging and digital reconstruction of liver tissue 12,13 . These allowed us to extend previous work, which predominantly addressed two-dimensional biological transport networks 1-3, 38, 39 , to three space dimensions. To overcome the limitations of spatially restricted samples of sinusoidal networks, we developed a Monte-Carlo algorithm for multi-objective optimization to generate simulated networks that faithfully reproduce the statistical features of digitally reconstructed sinusoidal networks (degree and edge length distribution). This algorithm is robust with respect to the weighting of the individual objectives. By varying the nematic alignment of simulated networks, we quantified the relationship between the nematic order parameter and the anisotropy of network permeability: for increasing nematic order, the permeability increases along the direction of nematic alignment, yet decreases in the perpendicular directions, displaying a linear dependence on the nematic order parameter. To characterize the resilience of networks in a scale-invariant manner, we introduced the notion of permeability-at-risk, i.e., we compute the residual permeability of the network as a function of the fraction γ of edges removed. This measure is more general than the familiar percolation threshold at which the permeability becomes zero almost certainly, or the fault tolerance considered in Tero et al. 2 , which quantifies the probability of percolation upon removal of an individual edge. The permeability-at-risk curves for different networks follow a characteristic, concave shape, monotonically decaying as function of γ and becoming zero at the percolation threshold of the network. Yet, individual curves do not collapse on a master curve, but reflect properties of individual networks. Computed permeability-at-risk curves match for both sinusoid and simulated networks. This serves as additional validation of our network generation algorithm. Using simulated networks that mimic real sinusoidal networks, we were able to vary the nematic order parameter of the networks. We found that weakly nematic networks not only have a higher permeability in the absence of perturbations as compared to isotropic networks, but also retain this property even if edges are removed. Thus, weakly nematic networks optimize performance without compromising resilience. Nonetheless, according to our minimal transport model, sinusoidal networks and their simulated counterparts are rather vulnerable to the removal of high-current edges when compared to regular lattices. We attribute the lower resilience of disordered lattices to the presence of a small fraction of edges that carry substantially higher currents. Indeed, we observed the same property in random sub-lattices of regular lattices, where likewise a small subset of edges carry high currents. We note that the sub-lattice represents an example of a random resistor network, for which a rich body of theoretical results exist 27,28 . Previous studies on optimal networks emphasized building and maintenance costs of networks, which were modeled as functions of total network length 9 . In the case of sinusoidal networks, however, we anticipate that other design constraints are even more important than building costs. The multi-objective optimization used here to simulate synthetic networks adjusts both node degree and edge length distribution and thus controls the total length of the network (by fixing the total number of edges and their mean length).
To compute network permeabilities, we used a minimal transport model based on linear Kirchoff equations with constant resistance per unit length. Generally, the apparent viscosity of blood flow depends on vessel diameter by the Fahraeus-Lindquist effect 40,41 , yet the variation of sinusoid diameter in the sinusoidal network is small 12 , which validates our assumption of a constant viscosity. The assumed linear relation between currents and pressure differences is valid for Newtonian fluids at low Reynolds numbers. However, blood is a non-Newtonian, shear-thinning fluid, whose apparent viscosity decreases with flow [42][43][44] . This implies a reduced effective resistance at higher flow rates, which would favor an even more inhomogeneous distribution of currents within the network. We speculate that adaptive mechanisms, like shear-dependent regulation of sinusoid diameter [29][30][31] , or sinusoid contraction dependent on local concentration of solutes 45 , may facilitate a more homogeneous distribution instead. Moreover, the diameter of red blood cells is only slightly smaller than the diameter of the sinusoids, which can result in transient clogging of sinusoids, especially where flow is high 32,33,46 or where the diameter of sinusoids is reduced by adaptive mechanisms. We propose that transient clogging (corresponding to a time-varying network 47 ), can likewise result in a more homogeneous time-averaged flow distribution. Notwithstanding these uncertainties, our simple model serves as first approximation, which provides a suitable tool to characterize the geometry of spatial networks, even if actual rates of blood flow in liver tissue should differ.
In addition to the sinusoidal network, liver tissue comprises a second network, not considered here, the bile canaliculi network, which transports bile fluid containing digestive enzymes 48,49 (for recent digital reconstructions see 12,13,21 ). Like the sinusoidal network, the bile canaliculi network spans across the entire liver lobule and contacts every single hepatocyte. Yet, bile fluid is toxic and must not enter the blood-stream, hence the two networks should nowhere get too close to each other. In two space dimensions, these competing design requirements could only be fulfilled by 'lines' of alternating network type, connecting source and sink. Such an architecture would only possess low resilience to damage. In three space dimensions, however, parallel network layers of alternating network type are possible. Such a design confers high transport capacity and resilience to damage. Indeed, signatures of such layered order were recently identified in liver tissue 13 . It must be stressed, however, that the geometry of sinusoid and bile canaliculi networks is much more disordered than a regular design of alternating network layers, as expected from self-organized networks that form in a tissue, where cells such as hepatocytes continuously divide. Such disordered networks are characterized by an inhomogeneous distribution of currents, where, at least in our minimal transport model, a small number of edges carry substantially higher currents than the majority of edges. Low-current edges, even if they dispensable for the permeability of the network in the absence of perturbations and contribute only partially to network resilience, play an important biological role nonetheless, e.g. for the uptake of metabolites from the blood-stream, which relies on diffusion through the fenestrated surface of sinusoids. We expect that also the low-current edges contribute to the supply of hepatocytes, provided these are connected to high-current edges by short distances. Future work will address the self-organization of pairs of space-filling, mutually repulsive, intertwined networks 50 , and study their transport and resilience properties, inspired by the bile and sinusoidal network in liver tissue.

Data acquisition
As described previously 12,13,34 , fixed tissue samples of murine liver were optically cleared and treated with fluorescent antibodies for fibronectin and laminin, thus staining the extracellular matrix surrounding the sinusoids. Subsequently, samples were imaged at high-resolution using a multiphoton laser-scanning microscopy. Three-dimensional image data was segmented and network skeletons computed using MotionTracking image analysis software 34 . The data sets analyzed in this study correspond to the same used in 13 .

Hard sphere packing model
We used Monte-Carlo simulations to compute the radial distribution function g(r) = g(r; R 0 , η) of a packing of equally sized hard spheres with radius R 0 and volume fraction η (η = 0.15, 0.20, 0.30, 0.35). In these simulations, n = 10 4 spheres were initially positioned at regular grid positions, then 2000 Monte-Carlo cycles were performed, where each cycle consisted of testing a Monte-Carlo move with Gaussian displacement (standard deviation σ = R 0 /3) for each of the spheres (acceptance rate about 75%).
For efficient computation, we exploited the fact that the radial distribution functions scales as g(λ r; λ R 0 , η) = g(r; R 0 , η) if the radius R 0 of the hard spheres is changed to a new value λ R 0 . Spline interpolation was used to interpolate g(r; R 0 , η) for intermediate values of η. A fit of the radial distribution function g(r; R 0 , η) for the minimal hard sphere model and the radial distribution function of the sinusoidal network resulted in R 0 = 9.0405µm and η = 0.2135 µm −3 .
In a final step, a subset of sphere positions was selected to match the node density of the sinusoidal network. For Fig. 2, a total of 1643 nodes were selected in a region of dimensions 420 µm × 450 µm × 90 µm, corresponding to a node density of 0.9659 · 10 5 mm −3 , which equals the node density of the entire sinusoidal network data set, including the void spaces occupied by the portal and central veins. For Fig. 3, a total of 1964 nodes were selected in a region of dimensions 250 µm × 250 µm × 250 µm, corresponding to a node density of 1.2571 · 10 5 mm −3 , which equals the node density of the sinusoidal network excluding the void spaces occupied by the portal and central veins. Note that the random selection of a subset of hard sphere centers does not change the normalized radial distribution function g(r).
For the set of simulated node positions P sim , we computed the Delaunay graph D(P sim ) and the minimum spanning tree M(P sim ). As a technical note, the computation of the Delaunay graph D(P sim ) can generate artifacts at the boundaries of the simulation domain with unusually long edges. To avoid this, we first mirrored all nodes in the set P sim at the 6 planes defined by the boundaries, thereby obtaining an enlarged set of nodes and then computed the Delaunay graph for this enlarged set. Finally, we retained only edges that connected original nodes, while discarding those that connect to one or two mirrored nodes.

Multi-objective optimization
We used simulated annealing to determine an optimal set of edges E * α that minimizes the multi-objective cost function C α defined in Eq. (4) for a given value of α with 0 ≤ α ≤ 1.
The optimization was initialized by a set of edges E sim = E 0 where E 0 = M(P sim ) ∪ E rand comprises the minimum spanning tree M(P sim ) and a random selection of edges E rand ⊆ D(P sim ) \ M(P sim ) chosen from the set difference of the Delaunay graph and the minimum spanning tree of the set of node positions P sim . For the simulations shown in Fig. 2, the number of edges in E rand was chosen to match the number |E \ M(P)| of edges of the sinusoid data set that do not belong to the minimum spanning tree. For the simulations shown in Fig. 3, this number |E rand | was chosen such that volume density of edges agree for E 0 and E.
In each step of the optimization procedure, we randomly selected a Monte-Carlo move that swaps a random pair of edges e 1 ∈ E sim and e 2 ∈ D(P sim ) \ E sim . This Monte-Carlo move is accepted with probability exp(−∆C/T ) if ∆C > 0, and with probability 1 if ∆C < 0, resulting in an update of the set of edges E sim . Here, ∆C is the change in the normalized cost function that would result from this move and T an effective temperature parameter.
The temperature parameter T was initially set to 1 and kept constant for an initial melting phase of 10 4 steps. Subsequently, T was reduced by a factor of 0.998 after every 1000 steps, until a final temperature of 10 −10 was reached. As validation tests, we confirmed that (i) the cost of the final network obtained at the end the annealing procedure equals the cost of the minimal-cost network encountered during the entire MC-optimization, (ii) multiple runs gave almost identical values of the minimal cost, and (iii) the computed Pareto front is continuous and concave, as predicted by theory. To account for the difference in range of the individual cost functions, C d and C e , we used a normalized total cost function with linear rescaling C α in the simulations, defined as where the normalized individual cost functions read and C e = (C e − U e )/(N e − U e ). The utopia and nadir points, U d , U e , N d , N e were determined in preliminary simulations. The normalized cost functions satisfy 0 ≤ C d [E * sim (α)] ≤ 1 and 0 ≤ C e [E * sim (α)] ≤ 1 for the set of edges E * sim (α) that minimize C α . We have the linear relationship In Fig. 2E, we used the range 10 −4 ≤ α ≤ 1 − 10 −4 . The value α = 0.5 corresponds to the value α * = (N e − U e )/(N d − U d + N e − U e ) for the cost function defined in Eq. (4), and was used in all figures, unless stated otherwise.

Supporting Information
Supporting information for Karschau et al.: Resilience of three-dimensional sinusoidal networks in liver tissue.
Supporting information movie S1. Sinusoidal network from Fig. 1C, corresponding to part of a liver lobule spanning from central to portal vein. The centerlines of the sinusoidal network are obtained from high-resolution three-dimensional imaging of mouse liver tissue and digital reconstruction 12,13 .
Network permeability computation biased for small sample volumes. Network permeability is an intensive measure that should be independent of network size for statistically homogeneous networks. This property is the basis for effective medium theories. In reality, however, small sample volumes can introduce a systematic bias in the computation of network permeabilities, see Fig. S1. This provides additional motivation of our use of large network samples, as well as for our network generation algorithm to simulate spatially unrestricted networks with prescribed statistical properties. Permeability-at-risk if low-current edges removed. We repeated the analysis of Fig. 4C, but additionally removed lowcurrent edges. If the networks were not perturbed otherwise, we find that their permeability almost did not change. However, if together with the low-current edges also a fraction γ of edges is removed, the permeability is substantially reduced, see Fig. S2. This shows that low-current edges partially contribute to the resilience of networks, allowing for limited re-routing of flow if high-current edges are removed.  Figure S2. Low-current edges contribute to resilience. We repeated the analysis of Fig. 4C, but additionally removed low-current edges from both simulated and sinusoidal networks. A. Left: We show again the permeability-at-risk curve for simulated networks from Fig. 4C (solid red, scenario i; λ = 0.2, mean of n = 3 simulated networks), where a fraction γ of edges is removed, starting with the edges that carry the highest current in the unperturbed network. We compare this curve to a second permeability-at-risk curve, where additionally 25% of the edges that carry the lowest current in the unperturbed network are removed (dotted, marked with open circle). Alternatively, we can likewise remove 25% of all edges, choosing again those that carry the lowest current in the unperturbed network, yet restrict the selection to edges that do not touch the boundary of the simulation domain (dotted, marked with filled circle). The fraction γ is always measured relative to the total number of edges in the unperturbed network. Hence, for the dotted curves, the number of edges of the perturbed networks equals 75% − γ of the number of edges in the unperturbed networks. Right: Same as left panel, but for random removal of a fraction γ of edges (scenario ii), with removal probability proportional to edge length. B. Same as panel A, but for the sinusoidal networks. Of note, the permeability-at-risk curves for the networks with low-current edges removed differ depending on whether we restrict the selection to low-current edges that do not touch the boundary (dotted curve, marked with filled circle), or not (dotted curve, marked with open circle). This can be attributed to the small dimensions of the tissue sample, notably in z direction, which results in a number of edges at the boundary that carry almost zero current.
Minimal model of metabolic uptake. We consider a minimal model of metabolic uptake from the sinusoidal network by surrounding tissue (for refined models, see e.g. Ref. 51 ). Let c(q,t) be the concentration of a molecule in solution in the flowing blood at time t at position q in the sinusoidal network. The dynamics of the concentration profile c(q,t) in the network is governed by diffusion, convection, and uptake d dt c(x,t) = −∇ · (vc) + D∇ 2 c − β c. (S1) Here, we assumed for simplicity a constant uptake rate β with units of an inverse time. Additionally, we assume that concentration is homogeneous across the cross-section of a sinusoid with constant area A. Thus, we need to consider only concentration profiles of a one-dimensional coordinate along a sinusoid edge. Correspondingly, the flow velocity v for an edge with unit vector e can be taken as v = I/Ae, where I is the current through the edge. For efficient computation on a network, we employ an upwind Euler scheme with concentrations evaluated at the node positions q j with fixed time step ∆t c(q j ,t + ∆t) = − ∑ k max{0, I k, j } [c(q k ,t) − c(q j ,t)] ∆t V j Here, the sum ∑ k extends over all nodes connected to node j, with l j,k denoting the length of the edge connecting both nodes. For each node j, we have introduced an associated volume V j = ∑ k l j,k A/2. Note that due to Eq. (7), these update rules conserve total mass ∑ c(q j )V j if β = 0. To increase the spatial resolution, we divided each edge of the original network (see Fig. 3) into n = 4 equal parts. For nodes j with small V j , Eq. (S2) is modified by directly using the analytical solution for diffusive exchange between a pair of neighboring nodes, which speeds up computations by allowing to choose larger time steps ∆t. Fig. 3 shows steady-state concentration profiles for the network flow from Fig. 3D. Remarkably, the spatial concentration profile varies homogeneously, despite the rather inhomogeneous distribution of flows in the network. Concentration decays approximately exponential along the flow direction; the decay length δ ∼ v/β is set by the ratio of a typical flow speed v and uptake rate β . Accounting for diffusion with D > 0, homogenizes the concentration profile even further, see  Figure S3. Computed concentration profile of a prototypical metabolite using a minimal model of metabolic uptake from the sinusoidal network. A. Steady-state concentration profile computed according to Eq. (S2) for a case without diffusion. Approximately, the concentration profile decays exponentially along the pressure gradient direction (black: steady-state concentration c(q j ) normalized by maximal concentration c max as function of position x = q j · e x , red: fit exp(−x/δ )). B. Same as panel A, but accounting for diffusion. Parameters: uptake rate β = 0.1 s −1 , diffusion coefficient D = 0 in A, D = 1000 µm 2 /s (∼ diffusion coefficient of oxygen) in B; for flow computation: sinusoid radius R = 3 µm, sinusoid cross-sectional area A = πR 2 , dynamic viscosity of blood µ = 3.5 · 10 −3 Pa s 18 , resistance per unit length κ = 8µ/(πR 4 ) (assuming Poiseuille flow), pressure difference ∆p = ∆p 0 L x /L 0 with pressure difference between PV and CV ∆p 0 = 100 Pa 14,17 , typical distance between PV and CV L 0 = 500 µm, spatial dimension of region of interest L x = 210 µm.