{'fraction_non_alphanumeric': 0.023528728286760007, 'fraction_numerical': 0.0050252715970487425, 'mean_word_length': 4.601529450048813, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Distribution efficiency and structure of complex networks

Flow networks efficiently transport nutrients and other solutes in many physical systems, such as plant and animal vasculature. In the case of the animal circulatory system, an adequate oxygen and nutrient supply is not guaranteed everywhere: as nutrients travel through the microcirculation and get absorbed, they become less available at the venous side of the vascular network. Ensuring that the nutrient distribution is homogeneous provides a fitness advantage, as all tissue gets enough supply to survive while waste is minimized. How do animals build such a uniform perfusing flow system? We propose a local adaptation rule for the vessel radii that is able to equalize perfusion, while minimizing energy dissipation to circulate the flow and the material cost. The rule is a combination of different objective cost functions that compete to produce complex network morphologies ranging from hierarchical architectures to uniform mesh grids, depending on how each cost is weighted. We find that our local adaptation rules are consistent with experimental data of the rat mesentery vasculature.

Flow networks efficiently transport nutrients and other solutes in many physical systems, such as plant and animal vasculature. In the case of the animal circulatory system, an adequate oxygen and nutrient supply is not guaranteed everywhere: as nutrients travel through the microcirculation and get absorbed, they become less available at the venous side of the vascular network. Ensuring that the nutrient distribution is homogeneous provides a fitness advantage, as all tissue gets enough supply to survive while waste is minimized. How do animals build such a uniform perfusing flow system? We propose a local adaptation rule for the vessel radii that is able to equalize perfusion, while minimizing energy dissipation to circulate the flow and the material cost. The rule is a combination of different objective cost functions that compete to produce complex network morphologies ranging from hierarchical architectures to uniform mesh grids, depending on how each cost is weighted. We find that our local adaptation rules are consistent with experimental data of the rat mesentery vasculature.
Large scale organisms need a distribution network to deliver oxygen and other nutrients that are crucial for the function of their tissues. The supply efficiency of the vasculature will dictate the survival or death of the cells surrounding the vessels of the network. Thus, biological transport networks are thought to have undergone a process of gradual optimization through evolution, culminating in organizational principles, such as hierarchy or Murray's law [1,2]. The broad strokes of biological distribution networks are to a large extent genetically determined. However, the vast number of blood vessels in the circulatory system (more than 10 9 ) suggest that the position and diameter of each individual link cannot be encoded in the DNA. Instead, development relies on local feedback mechanisms, in which each vessel is measuring quantities such as the flow and the nutrient density and responds by adapting the diameter. For example, increased flow through a vascular segment will result in improved conductivity of the vessel [3]. This adaptive mechanism allows the network to modify its structure and respond to environmental changes.
Murray was the first one to mathematically cast vascular adaptation as a minimization problem of the energy dissipated by the flow. Most of the recent studies still focus on the optimization of transport cost subject to a material cost constraint [4]. In the absence of fluctuations, it was shown that minimum dissipation networks have a treelike topology and cannot maintain loops [3,5]. Other attempts managed to preserve loops by imposing uniformity in flow field, however this leads in a uniform conductance network that does not have any hierarchy [6] drastically increasing the transport costs.
Experimental evidence show that biological flow networks exhibit both hierarchy and loopy structure [7], therefore suggesting the need for new optimization principles that go beyond dissipation and material cost [8]. Indeed, the main purpose of capillary networks is to efficiently supply tissues with nutrients such as oxygen, transported by the fluid. Details matter (where the tissue is, what the tissue function is), and uniformly distribut- ing nutrients (equidistribution) is important as underfed areas will become hypoxic and eventually die. The importance of modeling certain aspects of vasculature design and its effects on flow and nutrient transport (e.g in the case of fetoplacental networks [9], the liver [10], the kidneys [11], the torso and head [12] or brain microvasculature [13,14]) has been long recognized by the community.
Equidistribution of load has been examined in vertex models where the radii of the supply vessels are constrained to be fixed but the vertices are allowed to move [15]. However, solute perfusion networks have not been studied systematically from the vascular optimization perspective, in which each vessel adapts its diameter in order to regulate the nutrient supply, leading to hierarchical networks. This is crucial because often efficient nutrient perfusion is incompatible with network structures that minimize the cost of flow transport. This comes from the fact that in some systems the nutrients are dis-tributed in a different manner compared to the flow. For instance, anisotropy in the distribution of nutrients occurs due to non newtonian rheology effects like plasma skimming. In other cases, distant vessels from the solute source (towards the venule side) are often under-supplied with oxygen. As a result, the network must readjust, often increasing the diameters of the vessels in the hypoxic areas in order to meet the metabolic demand [16].
In this work, we derive from first principles a local adaptation rule that optimizes nutrient transport. At first, we describe a model for nutrient absorption through the vessels of a flow network. Then, we propose an energy cost function that consists of a perfusion term, the energy dissipation to circulate the flow and the material cost. We derive a local information adaptation rule for the diameters of the network links that minimizes the cost function. The different terms of the cost function compete to produce a variety of complex network morphologies. We explore the phenotypic transitions between the undersupplied and over-supplied regime. Finally, we validate the optimization rule in the rat mesentery network and we provide evidence that uniform perfusion optimization explains the angio-adaptation in realistic networks.

NETWORKS OF UNIFORM PERFUSION
We model the microcirculation as a network of N e linear fluidic resistors, connected at N v nodes. The link connecting nodes i and j corresponds to the vessel e ij of a specific conductance k ij . In the vascular networks that we consider the Reynolds number is small, and the vessels are long compared to their diameter. In this regime the conductances of individual vessels can be approximated by Hagen-Poiseuille's law k ij = πR 4 ij 8µlij , where R ij is the radius of the vessel e ij , l ij is its length and µ the blood viscosity. The viscosity does depend on the blood hematocrit and vessel diameter, but for the theoretical part of this work for simplicity we used the same average value for the entire network (see also Materials and Methods). Each vessel carries a current Q ij . When Q ij > 0 the liquid flows from node i to node j. The flow is conserved at every node i in the network: where the summation is over all nodes j that are connected to node i. Q b i is the external flow entering the node i (or leaving it if Q b i < 0). We define the weighted Laplacian matrix of the network L = D − A, D is the weighted degree matrix and A the weighted adjacency matrix such that A ij = k ij . Calculating the pseudo-inverse of the Laplacian and using the boundary condition vector of the flow allow us to solve for the pressures p i for all the the nodes: P = L −1 Q B , where P and Q B are N v ×1 vectors whose i th element is the pressure p i and net current Q B i respectively (Q B i is non zero only if node i is boundary flow input or output). Using the pressures and the conductances we can solve for the flow of each segment Q ij = k ij (p i − p j ). To include the nutrients, we need to expand the model by introducing a nutrient density field ρ i , defined at each node i in the network. The nutrient density is transported by the flow so that the product of the density with the current Q ij gives the nutrient current J in ij that enters from node i and propagates downstream to node j through the segment e ij , as seen in Fig. 1: where Θ(p i − p j ) is the Heaviside function, necessary to ensure that the nutrient current has the same direction with the flow entering the segment e ij . For simplicity in the results in the main text we ignore phase separation of red blood cells at the vessel bifurcations, however, the model is generalized to account for plasma skimming (see Materials and Methods 6). In perfusion networks, often nutrients diffuse through the vessels to the perivascular spaces. Therefore, the current J out ij at the exit of the e ij vessel is reduced by ∆J ij = J in ij −J out ij . Here the "in" and "out" superscript is distinguishing the nutrient current entering and leaving the vessel respectively. In the case of nutrient diffusion across a long pipe, the absorption depends on the pipe characteristics such as the length, the diameter and the vessel's absorption rate. Furthermore, the absorption depends on the fluid velocity and the density of solute. In this article we consider diffusion in low Reynolds number and we use the expression for nutrient absorption current as estimated analytically in [17]: where ξ is the absorption rate of the vessel e ij and λ = (2 3 µπ 3 ) 1/4 is a constant characteristic of Poiseuille conductance. The dimensionless ratio (πRij ) −1 Qij ξlij = t d ta is the ratio of the time t d needed to diffuse radially through the cross section to the time t a needed to be advected out of the tube. Our method for direct calculation of the nutrient distribution can be applied to other absorption functions and as a result it is possible to be generalized to a broad array of systems. We proceed to calculate the nutrient distribution and absorption through the whole network. Nutrients enter the network through the N B source nodes. The density ρ B m of the nutrients at every vertex m that belongs to those boundary nodes is known. The particles are transported by the flow and propagate to the rest of the network. Despite the nutrient absorption at the edges, at each node there is conservation of nutrient mass so that the nutrient current entering the node j equals the nutrient current exiting the node j, i J out ij = n J in jn . Using the conservation of mass on the nodes, we can define an effective Laplacian matrix L eff which contains the nutrient absorption as a weight on the edges for a general absorption function ∆J (see Materials and Methods). Next, we obtain the nutrient field in the whole network by inverting L eff and multiplying with the boundary condition vector J B that contains information of the nutrient external input on the nodes. The elements of J B i are zero for internal nodes, and equal to Obtaining the nutrient field ρ, we can calculate the nutrient absorption in each segment. The average absorption of the network is defined as . M is considered as the characteristic metabolism of the network and depends on the graph structure and the absorption rate ξ that is the same for all the vessels.

THE ADAPTATION RULE FOR UNIFORM PERFUSION NETWORKS
Our objective is to create networks that optimize a combination of cost functions: solute uniformity, dissipation, and material cost. This is achieved through a gradient descent algorithm that at each step modifies the hydraulic conductances of the vessels in such fashion so that the system moves closer to a minimum of an energy functional equal to a weighted sum of all three cost functions. The energy functional reads (see Supplement): The energy explicitly depends on the pressures, conductances and the nutrient current densities written as p i , k ij , and ρ j respectively. However, the pressures and nutrient densities are an implicit function of the conductances, as they can be calculated using the constitutive relations and Ohm's law. The first sum is the variance of nutrient absorption in the system.
The second sum treats the mechanical energy of the system and consists of the dissipation power required to pump the fluid from the input to the output i>j, i,j =1 k ij (p i − p j ) 2 and the material cost to maintain the vessels i>j, i,j =1 k γ ij to the uniform perfusion term is quantified by the coefficient ω. The coefficient α compares the dissipation to the material cost. In this work we consider γ = 1 2 so that the material cost is proportional to the total volume of the system [6].
The optimized network is the one that corresponds to a local minimum of the energy 4. This minimum should satisfy ∇H = 0 with respect to all the independent variables. We propose that vascular systems adapt and remodel based on a local information adaptation rule that acts on each edge and that eventually leads the system to a steady state close to a local minimum of the energy 4. A similar idea in the absence of a perfusion term was presented in [18]. To find an adaptive rule that results to an optimized steady state, i.e. a minimum of H, we perform a projected gradient descent. The {k ij } are treated as adiabatic variables [6], so at each step n of the iterative process we update the conductances as k n+1 ij = k n ij +δk ij , where δk ij = −ν ∂H ∂kij and ν is the characteristic step size of the gradient descent algorithm [6]. The rest of the variables of the system (p i , ρ i ) are assumed to vary more rapidly to remain at local equilibrium so ∂H ∂p k = ∂H ∂ρ k = 0. At the end of each step, we update the rest of the variables p n+1 i , ρ n+1 i according to the constitutive relations and the new {k n+1 ij }. Using the above formulation we derive the following local information adaptation rule (see Supplement): The growth or decrease in vessel ij conductance depends only on the vessel flow Q ij , local nutrient density ρ i , and geometric parameters of the vessel. Implicitly, this adaptation function also depends on the boundary conditions of the flow and the nutrient density.

RESULTS
To get an insight on the structure of the optimally perfused networks we use the proposed local adaptation relation 5 to numerically search for the minimum of the energy H. The simulations are performed on a triangular lattice of randomized edge lengths. We initialize the vessel conductances by assigning a value to each vessel drawn from a uniform distribution [0.9, 1.1]. Then we initiate the adaptation process and at each step we update the conductances according to the energy gradient. For the simulations we use the non-dimensional energy 4 (see Supplement). . C: Heatmap of the relative error between the experiment and the simulation results. Networks that operate in the advective regime (low absorption rate ξ) show good agreement with the experimentally measured network: there is a range of parameters α, ω for which the relative error is less than 2.5 %. D: Optimization result only for the minimization of dissipation without the uniform perfusion term. This network as expected does not preserve the topology of the experiment since loops vanish. E: Optimization only of the uniform perfusion term without dissipation and material cost. This network maintains the experimental topology however it does not preserve hierarchy, e.g. distant vessels with respect to the nutrient source have dilated diameters to increase their supply.
The adaptation process eventually leads to a minimization of the total energy H. During this process the uni-form perfusion term and the dissipation compete to impose space filling and hierarchy respectively. In the absence of flow fluctuations for γ < 1 loops are not stable [3] when dissipation energy is minimized. However, here we demonstrate that it is possible for a network to maintain loops for stationary, non fluctuating flow when the target function to be optimized includes a term that minimizes uniform perfusion. The absorption rate ξ affects the competition between the perfusion term and the mechanical costs changes, and as a result the network morphology. For high ξ the network is organised in a space-filling, non-hierarchical mesh while for low ξ the network becomes hierarchical. This behaviour is shown in Fig. 2. In the first set of panels A l , B l , C l , D l the absorption rate is low (ξ = 0.2). For small absorption rate the network organises in long channels and the number of loops drastically decreases, optimizing the transport cost. For higher absorption rate ξ = 0.3 in Fig. 2(C h ) the network transitions from the advective regime to the diffusive regime and loops coexist with long channels. In this regime, uniform distribution of nutrients is prioritized while the transport cost is maintained relatively low. For very high absorption rate the average metabolism increases and the uniform perfusion term dominates. As a result the network is organized in a uniform grid and the vessels that are under-supplied tend to increase their cross sectional area.
In order to classify the spectrum of micro-vascular networks we use the entropy of the flow as an order parameter. Informational entropy measures the amount of uncertainty in a situation or system [19]. For a water distribution network, the entropy S represents the uncertainty faced by the water molecule as to which path it takes to go from the input source to the exit sink. The total entropy of the flow system is calculated as S = Nv n=1 w n S n , where S n is the entropy of node n S n = − n,j P nj ln(P nj ) for Q nj > 0 and w n = Tn T0 is the the relative weight of node n, calculated as the ratio of the total outflow T n of the node n to the the total supply T 0 . Finally, here P nj = Qnj Tn , is the probability that a molecule that arrives at n will flow in the vessel e nj .
The entropy of the optimized networks for a range of values of α and ω, for low and high ξ as shown in Fig. 3. The eight networks for low and high ξ are displayed in Fig. 2 are labeled in the entropy heatmaps. Hierarchical networks, as for example 'C l ' exhibit lower entropy, as most of the loops disappear reducing the uncertainty in the possible flow paths. In this case, high ω means dissipation cost is prioritized to uniform perfusion. On the contrary networks in which uniform perfusion prevails, e.g. 'A h ' , have higher connectivity and as a result the multiplicity of paths is higher leading to higher entropy.

UNIFORM PERFUSION IN THE RAT MESENTERY MICROVASCULAR NETWORK
Structural adaptation of vascular beds involves responses not only to mechanical stimuli but also to the metabolic state of the tissue. In this section we show that the local adaptation rule of Eq. 5 that optimizes a linear combination of perfusion and mechanical cost is consistent with real vascular networks. To properly simulate perfusion and test the adaptive rule, data about the geometry (diameters and lengths of each segment), topology (connection matrix of vessel segments), and boundary conditions (volume flow rates, the discharge hematocrit in all vessel segments feeding the network and the volume flow rates for those segments leaving the network) are necessary. These data was provided for a 2-D section of the rat mesentery vascular network in [16]. This network was supplied by arterioles with inner diameters of 30 µm and drained by venules of 45 µm.
We first test the stability of the experimental network to an adaptation relation of the form 5. We start from the experimentally measured configuration of conductances and flows and we apply the adaptation algorithms for different sets of parameters {α, ω, ξ}.
The adaptation is terminated once the network stops evolving and (δκ ij drops below a certain threshold for every edge e ij ). Then we compare the final adapted radii R sim n of each vessel n with their initial, experimental values R ex n and calculate the relative error: δR = We find that there is set of coefficients for which the simulation result remains close to the experimental solution and that is ω = 12 * 10 −4 , α = 2.3 and ξ = 1.3 * 10 −3 , with error δR = 1.6%. This implies that the experimental network is an equilibrium point of our optimization scheme. In Fig. 4A we present the experimental network and in Fig. 4B we plot the best fit result. The agreement between the two images suggests that the experimental vascular network is a local minimum for our energy functional, and stable under the adaptation rule. This adaptation result combines hierarchy and mesh structure. The network is operating at the advective regime because ξ 1, and since the dissipation coefficient is large the transport of nutrients is more prioritised over the absorption [13]. Nonetheless, the inclusion of the perfusion variance term in the cost function is crucial to maintain the loops and prevent distant vessels from vanishing Fig. 4D.That maintains the distribution of nutrients homogeneous.
This result is not due to fine tuning, and carefully selected values for α, ω and ξ. As shown in Fig. 4C, there is a range of parameters α, ω for which the steady state network is close to the experimental state with an error less than 3 %.
To further verify that the experimental state is indeed a minimum energy configuration of our model, we go beyond checking for stability and instead investigate if the adaptation algorithm would reproduce the network starting from a different, but nearby initial condition. For this, we perturb the initial experimental radii of order (80%) so that the disturbed diameter is:R i = R i + η i R i , whereR i is the initial distribution and η i is a random variable sampled from [−0.8, 0.8], as seen in Fig. 5A, and use these as the initial diameters of the adaptation algorithm. For the adaptation process we choose in Fig. 5B we show the steady state network obtained by an adaptation process using the parameters α, ω, ξ that best agree with the experimental data in the stability test.The steady state network recovers the experimental topology and conductance distribution. Under this adaptation rule the average difference from the experimental diameters decreased from (80%) to (20%). The rugged energy landscape of this energy functional results in the adaptation algorithm getting trapped in local minima. In the absence of any annealing type process, exactly reproducing the experimental state is unlikely.
FIG. 5. Mesentery vascular network with randomized initial vessel diameters. A: Randomized initial conditions. The added noise on the experimental diameters is 80% of their initial value. B: The optimized network ends in a minimum energy state that recovers the experimental distribution of experimental diameters within an error is δR = 20%. For the adaptation process we chose the set of coefficients ω, α, ξ that produced the optimum network calculated at Fig. 4B.

DISCUSSION
In this work we investigate how vascular networks remodel under a local information adaptation rule directly derived from optimization principles. Former studies that derive the local angioadaptation rules from optimization principles consider only the fluid flow in the adaptation process, and ignore the nutrient concentration. Either they minimize the flow dissipation cost [3], [1] or they impose uniform flow distribution [6], [20]. These optimization schemes produce either tree-like networks in the case of dissipation, with loops only appearing in the presence of fluctuations [1], or a uniform grid in the case of uniform flow. However, often just the flow itself cannot capture the behaviour of the nutrients. In the microvasculature, nutrients typically diffuse to the tissue along the length of the vessels, and as a result one would typically expect some extent of nutrient deficit downstream, distally to network's input. Moreover, when plasma skimming is present, the red blood cell concentration is no longer directly proportional to the plasma current [21]. Unlike prior studies, our optimization principle results in an adaptation rule that includes the effect of hypoxia-induced angiogenic factors. [22]. To apply our system of adaptation equations iteratively required the solution of a system of flow equations that aimed to reproduce how nutrients are carried by the flow and how they diffuse through the vessels. We solve this system and obtain the nutrient field by inverting a sparse matrix instead of performing multiple iterations to propagate the solution node-wise starting from the input node to the exit (Materials and Methods). The computational efficiency of our method would permit generalization of our method to larger network systems.
Minimizing the variance of perfusion , alongside the dissipation and material cost, can produce a spectrum of transport networks depending on the value of the constants α, ω and ξ. This supports the alluring hypothesis that all the diversity of physiological vasculature can be produced by the same adaptation equation, and a handful of fitting constants. This diversity is produced by competition of the mechanical costs, which are optimized by a tree like architecture, and the uniform perfusion term , which nudges the network to remain tightly meshed (Materials and Methods).
In general, perfusion networks can be roughly classified in two categories: advective and diffusive [13]. The model in this work predicts a transition as we vary the absorption demand ξ. In the advective regime that appears for low ξ, the overall nutrient absorption is low, consistent with a low metabolism. In this scenario favors formation of hierarchical and long channels, as for the same parameters α and ω dissipation is relatively more costly than the perfusion variance term, which is smaller. For intermediate ξ, the metabolic demand is increased and therefore the network transitions to a regime in which perfusion becomes more important. In this case a mesh grid of small vessels is maintained, but can coexist with the formation of larger channels. If we further increase the absorption so that ξ becomes large, uniformity is favored over hierarchy and the network transitions to a well-connected capillary bed. The system's ability to maintain loops comes from the perfusion variance cost term competing with the material cost (Materials and Methods).
The crossover between advection to diffusion can be successfully captured from the entropy of the flow, an order parameter that we use to quantify the structure of the various micro-circulation networks. Near the critical metabolism there is a discontinuous jump in the entropy while moving from lower values of absorption to higher values, as shown in Supplementary Fig 3. The flow entropy is an efficient parameter to capture the transition from a hierarchical to a uniform network while other parameters such as graph connectivity fail. For instance in Fig. 2 networks B h and C h have the same number of links however in B h the vessel conductances are of the same order of magnitude, while in C h there is coexistence of two vessel groups those that are forming long channels high in flow and those which form a uniform grid, the difference in their conductances is almost two orders of magnitude.
To test whether the optimization rule is consistent with experimentally obtained data, we apply the hemodynamic adaptation rule to the rat mesenteric vasculature. We show that minimizing the energy cost function 4 that accounts for uniform perfusion at the advective regime reproduces the mesenteric vasculature structure to a relatively high degree. Remarkably, the experimental architecture was captured only with the use of three fitting parameters {ξ, ω, α}, and for a broad range of values. Pries et al in [23] manage to reproduce the mesentery structure with very high accuracy. However, in that work, the adaptation rule was ad-hoc, not directly connected or derived from optimization principles, and included 9 fitting parameters.
The adaptation model presented in our work is not necessarily complete: other mechanisms not included in the analysis may be relevant or even essential in the control of micro-vascular network structure. Specifically, highly metabolically active organs such as the brain, where fluctuations play an important role, may require additional elements in the model [24], or in the case of capillary beds, in which rheological effects of blood are more significant [6]. Nonetheless, the theory developed here explains to a good degree the vascular structure of the mesentery and paves the way to include uniform perfusion of nutrients as dominant mechanism for angio-vascular adaptation. The model merits further experimental investigation, in other animal tissues and organs, and even plants, which have nutrient distributing flow networks. Furthermore, understanding how biological networks optimize transport can be educative in other fields of research and technology. For instance in the field of energy it could suggest new optimized designs of the channels of redox flow batteries or give new insights in the optimisation of supply distribution networks.
In the course of preparing this manuscript we became aware of related unpublished work by Kramer et al, concurrently posted on the Arxiv. In the results presented in this work, it is assumed that the nutrients at each junction are distributed proportionally to the flow. This implies that if a flow density ρ o enters a junction, the outgoing currents from that junction will carry the same nutrient density, and that the total amount of nutrients in each vessel will be proportional to the flow. However, often in micro-circulation the nutrients are not distributed proportionally to the flow: the daughter vessel with the larger diameter receives disproportionately more red blood cells. If at node 'i' the nutrients are ρ i , the nutrient density at the daughter vessels {ij} that are supplied by i is determined by ρ j = Φ ij ρ i , where Φ ij is coefficient equal to: The sum in the denominator represents the total incoming current to node i and η is the skimming exponent, a measure of the non linear splitting of the RBC's at the bifurcations. For η = 0 we have no skimming effect. For the human vasculature a reference value is η = 0.26. For simplicity in our simulations we ignored the plasma skimming effect. Nevertheless the model is designed to account for it. It is worth mentioning that according to [25] for diameters over 30 µm the plasma phase separation effect is negligible for the rat mesentery network however in the work of Chang [? ]. In the equations that follow, we present the most general model with η = 0.
The nutrient current, flowing from node i to node j reads: Nutrients like oxygen get absorbed as they are transported through the vessels so there is a drop between the incoming and exiting currents is the amount of nutrients that is being absorbed along the length of the vessel. Summing 8 over the index i we obtain: Furthermore, we imply conservation of nutrient mass at each node so that the nutrient current entering the node j equals the current exiting the node j i J out ij = n J in jn . Then, Eq. 9 can be written as: where A ij = . For simplicity, we define B ij ≡ (1 − A ij )Φ ij Q ij Θ(p i − p j ) and F ij ≡ Φ ij Q ij Θ(p i − p j ), which simplifies Eq. 9 to i ρ i B ij − ρ j k F jk =0. This final expression allows us to solve the system in matrix form. Q ij is a square matrix with dimension N v × N v containing the information of the currents, A ij is of same size but contains information about the loss among the different edges. Φ ij is the matrix element that refers to the splitting law of nutrients at each bifurcation. Finally, ρ i is a vector element which corresponds to nutrient number at each node i of the network. This is the variable that we want to solve for in the system.
The matrix form of the system becomes: HereD is a diagonal matrix that has the effect of a diode. It contains information about the direction of the nutrient current ensuring that we keep only the positive components, also it contains the weights for the case of nonlinear splitting of RBCs at the bifurcations.
We can solve for ρ by imposing the boundary conditions on the density ρ B . For each boundary node belonging to the set of n b boundary nodes where current enters or leaves the system we introduce a pseudo node and we connect them with an edge. So the total number of nodes is nowÑ = (N v + n b ). We define the incoming nutrient current on the input node as J B m = −ρ B m Q b m . Then we write the transport equations in a compact form using the effective Laplacian matrix as defined below: where I is the unitaryÑ ×Ñ matrix. In order to invert the matrix expression in the brackets on the left side of the equation and to solve for all the values of ρ i we delete the columns and rows of the Laplacian that correspond to the boundary nodes. ment [25]. Our study required the network structure, the vessel diameters and only the flows boundary conditions and the discharge hematocrit at the feeding nodes from which we obtained the oxygen supply. The actual flows and solute distribution were calculated using the constitutive relations for the flow and nutrients, and found in agreement with the experimentally provided ones. In order to solve for the flows and the solute distribution, in principle one has to take into account the dependence of the relative viscosity on the discharge hematocrit and the diameters of the vessels, as described in [26]. However, we found that for the rat mesentery the corrections to the relative viscosity due to the discharge hematocrit are negligible. Therefore, in order to avoid the computational cost that is required to solve a nonlinear set of equations, we use the average viscosity over all the vessels of the network. For more details see Supplement.