How to Pare a Pair: Topology Control and Pruning in Intertwined Complex Networks

Recent work on self-organized remodelling of vasculature in slime-mold, leaf venation systems and vessel systems in vertebrates has put forward a plethora of potential adaptation mechanisms. All these share the underlying hypothesis of a flow-driven machinery, meant to prune primary plexi in order to optimize the system's dissipation, flow uniformity, or more, with different versions of constraints. Nevertheless, the influence of environmental factors on the long-term adaptation dynamics as well as the networks structure and function have not been fully understood. Therefore, interwoven capillary systems such as found in the liver, kidney and pancreas, present a novel challenge and key opportunity regarding the field of coupled distribution networks. We here present an advanced version of the discrete Hu-Cai model, coupling two spatial networks in 3D. We show that spatial coupling of two flow-adapting networks can control the onset of topological complexity given the system is exposed to short-term flow fluctuations. Further, our approach results in an alternative form of Murray's law, which incorporates local vessel interactions and flow interactions. This scaling law allows for the estimation of the parameters in lumped parameter models and respective experimentally acquired network skeletons.

Many recent studies on biological transportation networks have been focused on the hypothesis that vasculature is remodelled according to the flow-induced stress sensed by the cells making up the tissue (1). This self-organized process optimizes the structures for the task at hand (distributing oxygen and nutrients, getting rid of waste, carrying local secretion). The actual tissue response is dependent on the time-scales probed, as short-term changes usually concern rapid vessel diameter changes in response to pressure fluctuations or medication, while long-term effects e.g. due to metabolic changes may manifest in permanent diameter changes (2), usually leaving the vessel structure with a trade-off between efficiency and redundancy (3). Particular focus has been directed to the remodelling of the capillary plexus and other rudimentary transport systems in the early developmental stages of organisms, i.e. by studying complex signalling cascades involving growth factors like VEGF in vascular systems of mammals (4) or auxin in plants (5). Yet, the onset of refinement seems to be correlated with mechanical stresses (such as shear flow) as has been shown in a variety of model organisms like chicken embryo (1,6), zebrafish (7), leaves(8) and slime mold (9). Early theoretical approaches by Murray (10, 11) posited that diameter adaptation would minimize the overall shear and power dissipation of the system. Recent models using global optimization schemes on expanded vessel networks (where diameter adaptation may lead to link pruning) involving random damage, flow fluctuations or rescaled volume costs were able to account for the trade-off of shunting and redundancies (12)(13)(14). Further advances were made in empirical studies of local vessel dynamics, e.g. blood vessel systems (15-17) as well as derived by constructing Lyapunov functions (18) describing the networks' effective maintenance cost (19). It has further been shown that the outcomes of locally adapting networks are robust against variations of the initial topological structure (20) and that plexus growth and correlated flow fluctuations can provide elaborate hierachies (3,21). Many of these effects may also be seen in continous adaptation models in porous media (22,23). It is interesting and important to note here that these adaptation mechanisms may leave certain fingerprints in the form of scaling laws, both allometric (24) and geometric (25). Reviewing these works, it becomes apparent that only single networks were studied, involving volume or metabolic constraints. We here focus on the development and function of multicomponent flow networks, which are influencing each other based on their spatial architecture. Biologically speaking, these are systems consisting primarily of blood vessels and a secondary entangled, interacting system as found for example in liver lobule (26-29), the kidney's nephrons (30, 31), the pancreas (32-34) or the lymphatic system (35). More generally, any complex spatial network in 3D can be thought of as being entangeld with its environment via its spatial complement which itself has the character of a network. In this work we study the adaptation of two coupled spatial networks according to an advanced version of the discrete Hu-Cai model (19). Each network is subject to flow driven and volume-constrained optimization on its own. Meanwhile we introduce the networks' interaction in the form of a mutual repulsion of vessel surfaces preventing them from touching directly by their otherwise flow driven radius expansion. Though the onset of redundancy is primarily driven by the existence of fluctuations, we find mutual repulsion to greatly reduce the networks relative loop density. Time-lapse experiments (7) counting the pruning events and topology analysis on pruned structures (3,36,37) allow for some insights. We further generalize an important scaling law, which has been discussed again recently in this context (38): Murray's Law. This generalization enables us to predict the model parameters with high fidelity for Kirchhoff networks solely from a given graph topology and it's edge radii distribution. In the same manner we find a reasonable estimation for the parameters in experimentally acquired data sets of sinusoids in the mouse's liver acinus.

Results
A. Spatial coupling stabilizes spanning trees in pruned networks with flow fluctuations.. We model flow distribution networks as a composition of m rigid cylindrical tubes (referred to as edges) of given length l and radius r carrying a laminar flow, which are linked together by n branching points (vertices). Further we assume the flow to be non-pulsatile and hence a constant pressure drop ∆p over the tube defined by its entryand terminal-point hydrostatic pressures. The fluid considered has viscosity η and by including no-slip boundary conditions one can derive the volume flow rate f as the Hagen-Poiseuille law f = πr 4 8ηl ∆p (39). We assign a flow fe, a pressure-drop ∆pe and a conductivity ce = fe ∆pe to every edge. Likewise we assign a source or sink sj as well as a nodal pressure pj to every vertex. For further calculation of the flows and pressures in such lumped networks see (40). In order to model biological transport networks such as found in the liver lobule, see Figure 1a, we define a multilayer network consisting of two intertwined, yet spatially separate objects each consisting of edges e and e , with designated vessel radii re, r e . As underlying graph topology we choose from the graph skeletons of the triply-periodic minimal surfaces P ('dual' simple cubic), D ('dual' diamond cubic) or G ('dual' Laves) (41) to acquire two symmetric intertwined networks. Here, we start to model each network as a simple cubic lattice, see Fig. 1b and 1c , for further comment on the other structures see the supplement. We define a local affiliation of edges: Every quadruple q = (ei, ej, e k , e l ) of edges forming a loop in its respective network is bound to enclose an edge e n of the other network. We refer to all edges from the quadruple q as 'affiliated' with e n and vice versa. Each edge in any network is affiliated with up to four other edges in the other network (this will naturally change for different graph topologies). As all theses edges are simply tubes in our model, we define the distance between affiliated tube surfaces to be, ∆r ee = L − (re + r e ) [1] where L is the initial distance of the abstract network skeletons (equal to distance in case of simultaneously vanishing radii). In order to model a system of blood vessels which is entangled with a secondary secreting vessel network, we postulate that the respective tube surfaces must not fuse or touch directly, i.e. having ∆r ee = 0. We do not incorporate periodic boundaries. For the networks' optimization we follow the conventional ansatz of minimization of flow dissipation and network volume, while including flow fluctuations as in (19). Subsequently we begin by constructing a Lyapunov function Γ ≥ 0 for the combined system: introducing an 'energy density' E for the respective tube surfaces of the network edges [3] with F ee = 1 if edges e and e affiliated 0 else [4]  As in (19) we consider a constant length scale l for all tubes, so we may rewrite the conductivity ke = lce = π 8η r 4 e , which implies the changes of conductivity to be solely the result of radii adaptation. Additionally we set the volume penalty as ve = lr 2 e . We define a pressure gradient as ∆ϕe = k −1 e fe where ke = πr 4 e 8η are elements of the diagonal K. We do not consider cost variations via an exponent ke → k γ e as sugguested in other studies (12,19). We may minimize Γ using a gradient descent method choosing the radii dynamics to be ∂tre = −χ∇r e Γ, and vice versa for r e with auxiliary coefficients χ, χ > 0. If a given radius falls below a set threshold we prune it, i.e. we remove the respective edge from the graph. Details of the derivation are given in the supplementary material. Further we assume the adaptation of the vascular networks to depend on an averaged potential landscape instead of instantaneous configurations, which are bound to occur in real systems due to short-term metabolic changes or vessel blocking/damage. I.e. we assume a constant vessel radius between two adaptation events, while the flow rates change throughout the system due to changes in the sinks' magnitude. This incorporates the assumption of a time-scale separation between the radii adaptation (long-time changes, not to be confused with short term contraction/dilation) and changes of hydrostatic pressure. We define fluctuations in accordance to (14) and subsequently we update the ODE systems by substituting all occurencenes ∆ϕ 2 with ∆ϕ 2 . In this model all sink fluctuations are uncorrelated and follow the same probability distribution for every vertex j. Each sink has mean µ and standard deviation σ, see 'Materials and Methods'. In doing so we prevent shunting and the generation of spaning trees, which is caused by the typical 'single source -multiple sinks' D R A F T approach. Further, using this ansatz we prevent accidentally partitioning the Kirchhoff graph as opposed to realizing the sink-source configurations one by one (3). In order to perform a numerical evaluation of the resulting ODE system we define the units as follows: the radii in units of the grids' distance re = Lρe, fluctuating nodal in and outflow si = µςi, time t = h L 3 µ , the conductivity ke = η −1 L 4 κe and hence pressure ∆ϕe = µη L 4 ∆Φe and the networks' edge surface distance ∆r −ε ee = L −ε ∆ρ −ε ee . The ODEs may now be rescaled accordingly to the mean µ. The dimensionless form of the dynamical equations for each network is then of the form: [5] with the effective temporal response parameters λ0 = χ µη The coefficient matrices U , V bear the information about the sinks-source cross correlations, see 'Materials and Methods' 3. The terms λ 0 , λ 1 , λ 2 and λ 3 , describing the second network with indices e , are constructed analogously. To construct the networks in a symmetric manner we set the respective tube-lengths equally l = l . The optimization according to Eq. (5) is performed until the edge radii reach a stationary state. During this optimization we mark edges whose radius falls below a critical threshold ρc. These edges are no longer updated and are considered to have a radius of virtually zero (though for computational reasons they are here set to ρc = 10 −21 , these edges are refered to as pruned). The graphs consisting of mc and m c non-pruned edges are then analyzed for its cyclic structure. We measure the number of the independent loops (nullity) in the final network as (42), assuming we never break our graph into multiple components. We calculate the relative nullity of a graph as = mc−n+1 l 0 , where l0 is the initial number of independent loops. The coefficients λ2, λ 2 play a negligible role in the systems final topology, as was noted in (19) and is explicitly shown in Figure 2 for the case of λ1 = λ 1 = 0. Note that the parameters λ2, λ 2 influence the final vessel diameter as well as the time scales for reaching the stationary state. With that being said, we focus our numeric studies on a systematic scan of the effective network coupling and flow-fluctuation parameters. For that purpose we set the fluctuation rates identical λ3 = λ 3 and scale the repulsion relative to each other as λ1 = λ 1 . For our simulations we set the response λ0 = λ 0 = 10 −4 , and volume penalty λ2 = λ 2 = 10 6 (which will provide us reasonable computation times for reaching a stationary state and prevent the problem becoming too stiff). The repulsion's exponent is set ε = 2. We solve the equations of motion Eq. (5) numerically with the Heun scheme using a manually adjusted increment of h ≤ 10 −4 (decreasing time step for increasing λ1). The initial edge radii are chosen randomly (with all affiliated edge pairs fulfilling 0 < ρe + ρ e < 1.  All simulation results are presented here for one network only, as the results are highly symmetric. In Figure 3d we present the resulting nullity state diagram and it's respective analysis. Refer to the supplement for the results of both networks in comparison as well as the respective results on their dynamics. Figure 3 shows that the nullity may be influenced not only by the rate of fluctuation λ3, but by the repulsion of the two networks as well. The nullity's transition is continuous, starting from a tree-like state ( = 0) at small fluctuations and increasing monotonically in a logarithmic manner beyond a critical λ3 ≥ λ c . Eventually, the -trajectory saturates for large fluctuation rates λ3 ≥ λ s towards a maximal nullity max. The trajectories are altered regarding the onset of the transition as well as the saturation limit, see Figure 3e. To quantify these shifts we acquire the critical λc by identifying the logarithmic transition's root. The saturation λs is quantified by fitting the trajectories with a sigmoidal curve, providing us the possibility to extrapolate a saturation point which we define as 0.99 max, see Figure 3f. It appears that the critical point λc becomes a monotonically increasing function in dependence of the coupling parameter λ1. In Figure 3g we show the rescaled trajectories of 3e near the onset of the nullity transition. Introducing the reduced fluctuation parameter = λ 3 −λc λc we find the transition to follow a trivial logarithmic law (λ1, λ3) ≈ a (λ1) (ln − 1), with the coupling dependent scale a shown in the inset of Figure 3g. Generating these state-diagrams for the respective D, G-skeletons we found no any qualitative differences arising from altered plexus topology.

B. Generalized form Murray's law allows estimation of adap-
tation parameters of real networks. Our model of spatial coupling also points to a new form of scaling at vessel branchings. We refer hereby to the concept by Murray's law, connecting the radii ρ0 of the parent vessel splitting into at least two child branches with radii ρ1, ρ2 as: Originally a cubic scaling exponent α = 3 was predicted (10) whereas relative costs models (12, 19) sugguest α = 2 (γ + 1) [7] D R A F T We find the transition to follow the scale a is shown in the inset. while discarding flow fluctuations. Testing these scaling laws for the network-skeletons of the sinusoids and bile canaliculi in the liver-lobule, by fitting the exponent α in Eq. (6) for every branching of degree three (Y-branching), one can see that there is significant deviation from the predicted exponent α = 3, see Figure 4. The aquired values are α ≈ 3.76 for sinusoids and for bile canaliculi α ≈ 3.32. As capillary systems were already known to defy the cubic relationship (25), we suspect this to be correlated with the reticulated nature of these networks. Further, this deviation is not well described by the cost exponent γ. In accordance to Eq. (7) we would deduce γ to be smaller than one. Yet in this model's context the γ induced loop transition occurs at γ ≥ 1 (12, 19) so the network would have to be a mere spanning tree to be valid. As the sinusoids and canaliculi are highly reticulated networks we propose that these loops are actually fluctuation induced (13, 14) and subsequently altered by the networks' mutual interaction. We can deduce from our pruning model a new set of coefficients aj which are dependent on their corresponding edge's neighbourhood and the respective coupling strength, as well as the global structure of sinks and sources (which were assumed to be uncorrelated and identically distributed). This procedure greatly alters the form of Eq. (6) and we derive a new scaling law derived from the steady-states of the ODEs in Eq. (5). We recover the cubic exponent of the original model with Histograms and fits ℵ(µ, σ) presenting the estimated simulation parameters λi, rescaled axis is xi = log10λi: (a) Parameters for an ideal Kirchhoff network, dual Laves-Graph topology, single source-multi sink system, the system was initialised symmetrically with x1 = 4 , x2 = 6 and x3 = 2 (b) Parameters for the sinusoidal system near the central vein, we found the distributions with mean values x1 = 3.08 ± 1.37, x2 = 3.43 ± 0.73, x3 = 1.12 ± 0.84.
information of the respective vascular sections, see Material and Methods.

B.1. Estimating parameters from ideal Kirchhoff networks.
We tested the fesibility of ansatz Eq. (8) by simulating the pruning on a dual Laves graph topology (3-regular), with n = 206 (n = 210) vertices and m = 259 (m = 261) edges and setting the parameters symmetrically to λ1 = λ 1 = 10 4 , λ2 = λ 2 = 10 6 , λ3 = λ 3 = 10 2 . The sources were positioned in random vertices of the system. Edges of the respective networks affiliated with each other similarly as before, though not by penetrating faces but by finding the first nearest neighbours of edges inside a perimeter δ. We numerically (43) find the roots of Eq. (8) for a set of positive definite λ1, λ2, λ3 and λ 1 , λ 2 , λ 3 respectively. As we do not intend to use information on the direction of the currents at the sink-nodes (and this information is not available in the real system) we will solve Eq. (8) for the seven relevant sign permutations at each branching, see supplement. For further evaluation only the fit of highest quality (function value) is used. We use a logarithmic rescaling xi = log10λi in order to find a symmetric representation of the histogram's data. Doing so we fit a normal distribution ℵ(µ, σ) to the distributions maximum and we find strong agreement with the actual parameters for both networks, see Figure 5a. B.2. Estimating parameters for sinusoidal networks. We use the same approach to estimate the parameters λ1, λ2, λ3 for the sinusoidal system in the liver lobule of mouse, based on the available information of topology and radii distribution. We first make some simplifying assumption about the position of sinks and sources as well as cropping the network-skeletons to remove degree two vertices, see Figure 6: First we identify the vertices in the sinusoidal network which are closest to the central vein (CV) and label them as sinks. We calculate their center of mass (CMS) and use it as the center of sphere of radius R (here chosen as 390µm). Any other components,vertices or edges outside this perimeter are discarded, see Figure 6a. We next identify all branching points in the sinusoidal network and all paths p = (ei, ..., ej) consisting of edges ei which start from these points. We proceed for the canaliculi the same way, see Figure 6b, and then check for each segment of a path whether there is another segment of another network's path inside a perimeter δ (here 25µm). If so, these paths count as affiliated, see supplement for further details. We merge all edges along a path towards a single edge by using the conventional addition theorems for series of conductivities, see Figure 6b. Proceeding like this we end up with a reduced sinusoidal network, with n = 318 and m = 452. From this we calculate the repulsion and fluctuation terms in Eq. (8) and solve for the parameters again under consideration of all sign permutations for aj. The solutions' histograms are presented in Figure 5b, where we find the parameters to be roughly distributed at x1 = 3.08 ± 1.37, x2 = 3.43 ± 0.73, x3 = 1.12 ± 0.84. In context of the nullity state diagrams we could place these in the transition zone, close to the border of the nullity-transition, see supplement. Yet, the large standard deviation and asymmetric character of the error (due to the logarithmic scaling) provide a certain level of uncertainty. Further, we found the choice of perimeters R, δ to interfere with the distributions presented in 5b. The imposed distribution of sources in the experimental system is also not strictly valid as one should only consider the vertices at the periphery as such.

Discussion
We have shown that spatial coupling presents another potential mechanism of controlling the topological complexity of optimal transport systems in 3D. Though the onset of redundancy is primarily driven by the existence of flow fluctuations, we have shown that mutual repulsion of vessel surfaces reduces the networks relative loop density and provides another method by which a system may tune its ultimate architecture. It's also possible to retrieve tree-like states at high fluctuation levels, which imposes a new stabilization mechanism for spanning trees in noisy networks. We considered the special case of 'single source to multiple sinks' in combination with simple cubic lattices as plexi as the simplest case possible. No qualitative differences could be found in the phase diagrams in comparison with the dual diamond, or dual Laves graphs. The Lyapunov ansatz provides a generally applicable tool in network optimization, and should properly be tested for other boundary conditions or graph geometries which resemble realistic structures. Our model may also be seen as a toy-model for intertwined flow networks as found in the mammalian liver lobules and other related organ structures. Our approach further enabled us to derive a more general form Murray's law, directly involving flow interactions and environmental interactions. We find this technique to predict the interaction parameters with high fidelity for simulated Kirchhoff networks given their topology and respective edge radii distributions. In the same manner we find reasonable estimations for the parameters in experimentally acquired data sets of sinusoids in liver lobule of mice. Hence one could use this methodology as an effective classification of spatially adapting network structures. Apart from a refinement of the presented geometric coupling approach, an explicit involvement of hydrodynamic-chemical feedback between the distribution networks caused by the respective partner network should be taken into account: Biological systems such as the capillary networks in the liver lobule, present complex dual transport systems where the actual flow rates are not necessarily influenced by their respective partner network's flow rate (44) but rather by the concentration of bile salt components transported (27) as well as secretion rates of hepatocytes. Subsequently, future studies have to consider concentration/pump-rate dependent flow rate perturbations in the optimization model used so far, as well as a direct postulation of feedbacks in the adaptation dynamics to account for the actual network structures.

Materials and Methods
A. Sample Preparation, imaging and segmentation. Mouse livers from adult mice were fixed by trans-cardial perfusion, sectioned into 100 mm serial slices, optically cleared and immunostained, as described in (29). To visualized the different tissue components, the tissue sections were stained for nuclei (DAPI), cell borders (Phalloidin), bile canaliculi network (CD13), and the extracellular matrix (ECM, fibronectin and laminin) facing the sinusoidal network(45). High-resolution images of the liver lobule (Central vein -portal vein axis) were acquired at by using confocal microscopy with a 63x/1.3 objective (0.3µm voxel size). Finally, the resulting images were segmented and network skeletons calculated with the Motion Tracking software as described in (29) and (46).

B. Calculating currents and fluctuations in lumped networks.
Defining fluctuations in accordance to (14) we only consider s-configurations in which there exists one source-node (here j = 0) and all other nodes being sinks with the following characteristics: We assume the fluctuations are uncorrelated and follow the same probability distribution we have for the mean µj = µ , standard deviation σj = σ and correlation coefficient ρ jk = δ jk . Subsequently we may calculate the average squared pressure, writing: with coefficient matrices and with the auxiliary conductivity tensor It may be noted here that A e U describes the pressure p 2 e in case of a constant source-sink landscape in the absence of any variance σ 2 .
Hence A e V describes the pressure perturbation caused by fluctuations of strength σ 2 . For the full derivation of Eq. (11) and the respective coefficient matrices see the supplementary material.

ACKNOWLEDGMENTS.
Our thanks go to Hernan Morales-Navarrete (image reconstruction) and Fabian Segovia-Miranda (experiments) for providing us with the sinusoidal and canaliculial skeleton datasets. We'd like to thank Marino Zerial, Julien Delpierre, Quentin Vagne , Dora Tang, Andre Nadler, Mark Werner and all members of the Modes Lab and Zechner Lab for their comments, helpful discussions and resourceful feedback during the process of creating this work and draft. In this section we will give a detailed derivation of the dynamic equations for the networks' radii and an introduction to the common notifications used. Having established a fully connected graph which we will use for our simulations we choose an arbitrary direction for each edge. Every cyclic path, meaning a progression of incident edges which starts and terminates at the same node, we will refer to as circle or loop. This special class of networks are often referred to as Kirchhoff or lumped networks. We formulate the constraint, e Bjefe = sj ⇔ B · f = s [1] with s ∈ R n , f ∈ R m , B ∈ R n×m referred to as current law, which can be interpreted as a mass conversation law for every vertex. Nodes with positive values of sj are to be referred to as sources, negative ones as sinks. B denotes the graph's incidence matrix. This setup is also referred to as value boundary problem (Neumann problem (1)), with terminal nodes (sj = 0) and internal nodes (sj = 0). In linear lumped systems one can derive Ohm's law for every edge,

D R A F T
with the conductivity matrix C ∈ R m×m being diagonal C ef = ceδ ef . From the Hagen-Poiseuille law we know the conductivities to be ce =  (2) we can see that: The n × m matrix B · C is non-invertible as it is a non-square matrix of rank n − 1. A solution may however be constructed by means of the generalized inverse (denoted here as [·] † (3) ). Via the dissipation minimization argument one can deduce the potential vector to be: Proving equation Eq. (4) to be the unique solution which minimizes the system's dissipation is given in (4,5). On a further 3 note, we'd like to mention that this calculation is not just restricted to laminar viscous flows but any stationary transport 4 system, e.g. diffusion. This phenomenon is also refered to as the Thomson principle (6, 7).

6
The optimization of the flow system will require the formulation of a positive definite cost function. For convenience we will write a Lyapunov function Γ as Γ = P + P + E [5] with P = We may further use the vectorial notation for dissipation-volume terms P ,P , using Ohms law to formulate it in terms of the nodal sinks/sources (here just for P, derivation for P' is performed analogously), with q = a 8η π 1/2 and ke = πr 4 e 8η as entries of the diagonal K. We used ∆ϕ = K −1/2 B · K 1/2 † s and f = K · ∆ϕ. We calculate the (pseudo-)time derivatives of P to be

Felix Kramer, Carl Modes
The derivative of the generalized inverse B · K · B T = A being (8), Fortunately the projector terms vanish as we have, Together with the identity A † = A T † the total time-derivative of P becomes, With partial derivatives simplifying this formula as: As ∆ϕ = B T · A † s, we may write the total time-derivative as ∂ h re and re-substituting q this becomes 1 The coupling component's derivative E may be calculated with partial derivatives, We then acquire the ODE system presented in the main body of our paper by setting α = 4aη π and β = b(ε−1)η πl and α , β accordingly,

Uncorrelated and coupled flow fluctuations
In this section we give a detailed derivation of the analytic form of the mean squared pressure in case of uncorrelated, identically distributed sink fluctuations as introduced in (9). When considering the sink conditions as well as the source constraint defined in the theory section, we may write the moments as, Hence we may calculate the squared-mean pressure by using the auxillary conductivity tensor A e jk = K −1 Ordering the terms for µ and σ respectively we can acquire the coefficient matrices U and V , We further suggest to expand this ansatz by introducing additional sources which act as clones of the very first one, i.e. we will have sp = sq using the indices p, q for sources and u, v for sinks. Then conditions Eq. (24), Eq. (26),Eq. (27) will become for a sources and b sinks (with a + b = n), And hence we may calculate the mean squared pressure and its coefficient matrices respectively as, In case of more elaborate configurations as proposed in (10), i.e. spatially correlated s-configurations, the authors suggested using the arithmetic average of a set of potential sink realizations sx to calculate the mean squared flow as, which demands the calculation of each single flow realization with the Moore-Penrose inverse. This problem may be more elegantly formulated for the analogue mean squared pressure, using Eq. (28) to write, This allows us to more quickly handle averages using just one arithmetic mean over a set of sink-source realizations sq as which might improve the computation of Eq. (28) in case of small iteration numbers x.

Limit case of uncoupled flow networks
From the rescaled ODE set of two coupled networks we may easily reproduce the uncoupled case for the two networks by setting λ1 = λ 1 = 0. With focus on just one network we have the radii adaptation, which resembles the gradient flow driven optimization model proposed in the original Hu-Cai model (11) in combination with 10 the flow fluctuation notation proposed in (9). We calculate the adaptation of simple cubic networks with a single corner source 11 node (sinks otherwise) for a systematic scan of λ2,λ3, see Figure S1 and S2. In particular we can show that in agreement with

12
(11) only an increase in the flow to fluctuation ratio λ3 results in an increase of nullity and that the continuous transition 13 observed is independent of the systems effective volume penalty λ2, see Figure S1. The system's transition is logarithmic for 14 λ3 > 1 and saturates for large λ3, as was shown in the main body of this study. This transition is qualitatively conserved for 15 all other periodic lattices, i.e. fcc, bcc and diamond lattice follow the same transition type. The actual dynamics and pruned structures are represented for several sample parameters in Fig. S2. Fluctuation induced pruning into spanning trees usually takes place on shorter time-scales then their reticulated pendants and displays less avalanche-like events of a edge radii suddenly collapsing. Nevertheless, there is also a quantitative variation of the final nullity, depending on the total network size, see Figure S3a. How then is the trajectory dependent on the actual system size? An answer could be given by looking into the fluctuation generating machinery which manifests the positive feedback in Eq. (40). Adding up U and V and resorting the terms around λ3, we rewrite ∆Φ 2 e as, By introducing further abbreviations to improve readability we rewrite Eq. (41) as, h (λ3, n) = 1 + λ3 n − 1 Hence we have h (λ3, n) −→ 1, for large networks with small fluctuations λ3 n, so the solution in Eq. (43) becomes, This impact on the (positive) feedback terms suggest a size dependent change in the overall network structure. Though we can't exclude the possibility of other finite size effects, such as the ratio of surface (nodes with a degree smaller than the maximal possible for the lattice) to volume nodes(nodes with maximum degree). We can observe the transition to follow a logarithmic relation for any network size (node number n), From fitting the trajectories in S3a one can quantitatively see that the transition curves converge with increasing n at given λ3, suggesting our approximation to be a plausible cause.  Table S1. Specifics for intertwined graph topologies of the triply periodic surfaces P, D and G

26
In Figure S5 we present the nullity transition for the intertwined cubic networks, displaying yet again the retardation of the Phenomenologically, the transition resembles an asymmetrical sigmoidal for λ3 > λc, whose saturation value, steepness and 32 x-offset are functions of λ1. Note that, though qualitatively similar, the curves development for increasing λ1 is different in the 33 two networks. In particular we observe a significant change of the saturation onset in the second network, see Figure S5c, S5d.

34
Next we have listed the systematic λ1, λ3 of the final stationary state after pruning. In Figure S6a,S6b we present the nullity

Scaling laws in coupled, noisy networks
This section centers around the scaling laws resulting from the formalism. One can show that introducing fluctuations and coupling alter the classical form of Murray's law ρ 3 0 = ρ 3 1 + ρ 3 2 in the following way: Given the Kirchhoff current law we and rewriting it via Ohm's law into a reduced form, we get for all nodes i and reference flow µ, with Θie = ±1 0 distinguishing between in-and outgoing flows on the relevant edges, see Figure S13. In order to acquire the cubic form we we substitute one multiplicand ρe with the result of the stationary states equation, The structure of Eq. (58) is unchanged by introducing multiple sources (clones) into the system, with exception for the specific entries of U , V , see Eq. (36). For experimental validation of Eq. (56) it will be necessary to know the networks vessel radii as well as the sink/source pattern ( although it may be sufficient to know where the system's source is and to consider every other node as sink). One should consider simple Y-branchings of low sink/source value (points of neglectable secreting/leakage) such that |aeρ 3 e | 1, further setting the index for the largest vessel to zero ( and accordingly increasing for the other vessel pieces at the branching) may write, Doing so we can demonstrate the quality of this procedure for sink-source conditions, see Figures S17, S18, S19 showing that 57 increasing number of cloned sinks does not alter the quality of the estimation. We would like note though that the estimation 58 gets considerably worse when encountering a parameter regime which impose neglectable interactions such as λ1 < 10 3 or 59 λ3 < 10 −1 . according to the given segmentation. These vertices will be declared sinks while all other vertices in the sinusoidal section will 63 be declared sources, which will be the fluctuating elements. We also now consider multiple sinks, which are set identical in 64 their dynamics for computational simplicity, see supplement. We calculate their center of mass (CMS) and use it as the center 65 of sphere of radius R. All other vertices and edges of either the sinusoidal or the canaliculi network which lie inside this sphere 66 are the only components to be considered for future analysis. Other components outside the range of interest as well as edges 67 penetrating the sphere's surface are discarded. We next identify all branching points in the sinusoidal network and all paths 68 p = (ei, ..., ej) consisting of edges ei which start from these points. We proceed for the canaliculi the same way, and then check 69 for each segment of a path whether there is another segment of another network's path inside a perimeter δ. If so, these paths 70 count as affiliated and we define a distance between affiliated paths by calculating the minimal distance between two segments 71 of these paths. For a reference length L we use the minimal distance found between the two networks (as we did before by 72 using the lattice constants in the simulations). Then we coarse-grain the network of interest by merging all edges along a path.

73
In practice this means that we use the conventional addition theorems for series of conductivities, 1 κ ef f = i  Finally this scheme could also be considered in the case of uncoupled networks λ1 = 0, so that we may write this as,