Distribution networks achieve uniform perfusion through geometric self-organization
DDistribution networks achieve uniform perfusion through geometric self-organization
Tatyana Gavrilchenko , and Eleni Katifori Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, Pennsylvania 19104 Center for Computational Biology, Flatiron Institute, Simons Foundation, New York, NY, 10010 (Dated: September 10, 2020)A generic flow distribution network typically does not deliver its load at a uniform rate across aservice area, instead oversupplying regions near the nutrient source while leaving downstream regionsundersupplied. In this work we demonstrate how a local adaptive rule coupling tissue growth withnutrient density results in a flow network that self-organizes to deliver nutrients uniformly. Thisgeometric adaptive rule can be generalized and imported to mechanics-based adaptive models toaddress the effects spatial gradients in nutrients or growth factors in tissues.
Multicellular and macroscopic living organisms arecontinually faced with the challenge of how to uniformlydistribute nutrients throughout their entire volume tomaintain metabolic function while minimizing waste ofresources. For this, they have evolved complex flow sys-tems in the form of dense and space-filling networks ofsmall vessels, termed capillaries, which distribute fluidladen with nutrients. In such perfusable systems, nu-trients are carried with the flow through the capillar-ies and gradually diffuse across semipermeable capillarywalls where they are depleted at a metabolic absorptionrate. In the absence of fluctuations in the flow and othermitigating factors, most network architectures, heteroge-neous or uniform, will not distribute nutrients equally:in general, the tissue that is upstream will absorb morenutrients than the tissue downstream (Fig. 1).When biologically-inspired microfluidic networks havebeen considered in the past, the emphasis has been onshear stress, resistance, flow rate, and pressure distribu-tions, e.g as in [1, 2], and the functionality of nutrientdelivery has been largely ignored. Recent work has in-vestigated a network adaptation algorithm for uniformedge flow, showing that it is possible to tune edge conduc-tances while maintaining the network structure to obtainequal flow over all edges [3]. In the context of microvas-cular networks, the finite size of red blood cells may aiduniform flow [4–6]. However, uniform flow does not guar-antee uniform nutrient distribution: if nutrients are con-tinually depleted by absorbing tissue, downstream tissuewould generally still have a deficient supply.Recent work has considered the optimal architecturefor uniform nutrient perfusion in plant leaves [7]. Therethe authors have considered networks of fixed edge andnode positions and used Monte Carlo methods to discoverthe optimal edge conductance distribution that deliversnutrients uniformly. The solution for a uniform planarnetwork with a line of sources on one side and a uniformdistribution of sinks in the bulk exhibits a strong spatialgradient of edge conductance from the proximal to thedistal network sites, with an order of magnitude varia-tion in vessel radii. While adding hierarchy in the vesselradii can result in uniform perfusion, such stratificationis not always possible due to developmental or physicalconstraints. Similarly, fabricating a network design witha wide range of vessel sizes may be difficult depending on the experimental setup. For instance, the method ofcasting networks in sacrificial ink can yield a printed di-ameter of up to only 2.5 times the nozzle diameter [2].In this work, we present a model of perfusion that con-siders not only the amount of nutrients leaving an edge,but also the average supply of nutrients in the surround-ing tissue fed by multiple edges. We demonstrate that anarbitrary initial network can self-organize to achieve uni-form nutrient perfusion using a simple geometrical andbiologically-plausible adaptation rule based on local in-formation. Unlike previous models of adaptation in net-works, we choose to constrain edge radii to be constantthroughout a network, but allow freedom in the edgelengths and network connectivity. Similar vertex modelsof local cell neighbor interactions are used to model thebehavior of tissue sheets and typically incorporate me-chanical cues for tissue adaptation, such as wall tensionand cell elasticity [8]. These models are able to recoverthe geometric structure of tissue sheets [9, 10] as well aspredict complex bulk properties such as collective motionpatterns of cells [11, 12]. In this work, we introduce a dy-namical rule reminiscent of such vertex models to obtainuniformly perfusing networks in an abstract setting ofareas of tissue separated by channels where the extracel-lular fluid can flow. Starting with an arbitrary network,the algorithm tunes vertex positions under forces that,as we show, stem from differential growth of the tissue.A segment of tissue receiving more nutrients grows at afaster rate, thus effectively pushing the channels at itsboundary further apart, increasing the nutrient deliveryat underfed tissue. When the forces become locally bal-anced the system reaches equilibrium, and the final net-work state obeys uniform nutrient distribution.The first step of perfusion modeling is determining therate at which nutrients leave the network through theedges, or the edge absorption rate. The exact form ofthe absorption rate depends on the physical details of thesystem; previous work has used a rate proportional to theinitial edge concentration [13], and a more biologically-motivated absorption rate is presented in [7]. In capillarysystems modeling oxygen transport to tissue, the mostbasic Krogh model predicts a linear decay in oxygen con-centration along a capillary [14], but realistic extensionsresult in more complicated forms [15–19]. Here, for acapillary of length L ij connecting nodes i and j , we use a r X i v : . [ phy s i c s . b i o - ph ] S e p an exponentially decaying absorption rate: C ij ( L ij ) ≈ C ij (0) e − αL ij /Q ij (1)where C ij (0) and C ij ( L ij ) are the nutrient concentrationat the capillary inlet and outlet, respectively, and Q ij is the capillary flow rate. In general, the parameter α depends on the vessel radius, the nutrient diffusion rateacross the vessel membrane, and the tissue metabolism;however, as an initial simplification it is set to be constantat the value α = 0 .
1. We note that eq. 1 is consistentwith the edge absorption rate used in [7] in the high flowvelocity regime.Given the nutrient decay rate through a capillary, theedge nutrient absorption rate φ ij is calculated by φ ij = Q ij (cid:18) C ij (0) − C ij ( L ij ) (cid:19) . (2)To compute this value for every edge of a flow network,first Q ij for each edge is computed using a combinationof current conservation at nodes and Ohm’s law overedges. Current conservation implies that the flow ratethroughout a single capillary is constant. We considernetworks with a single current source and sink and set afixed initial nutrient concentration at the current source.Edge resistance is determined by the HagenPoiseuille law R ij = 8 µL ij / ( πr ), where µ is the fluid viscosity and r the vessel radius. The parameters r and µ are heldconstant in this work, although the model can easily beadapted to consider varying r as well.All capillaries leaving the node i will have the sameinitial concentration, so we simplify the notation to C i (0) = C ij (0). The drop in nutrient concentrationacross each edge is computed iteratively through the net-work, starting from current source nodes { n } , which areprescribed an initial C n . In this work all networks areset to have a single source node with an initial nutrientconcentration C = 100. At node i , conservation of nu-trient mass per time is obeyed by enforcing the relation (cid:80) k , Q ki > C ki ( L ki ) Q ki = C i (0) (cid:80) j , Q ij > Q ij , where the no-tation Q ij > i and into node j . If node n is a current sinknode and has no outgoing edges, the node concentrationis C n (0) = q n (cid:80) k , Q kn > C kn ( L kn ) Q kn , where q n is the sinkcurrent at node n .Finally, the measure of network absorption uniformitywill be captured by the nutrient absorption density Φ f ofeach face f . It is defined as the total amount of nutrientsreceived from adjacent edges divided by two (since eachedge supplies two faces) scaled by the face area A f :Φ f = 12 A f (cid:88) ( ij ) ∈ f φ ij . (3)Since this definition of Φ f implies that nutrients leavethe system via boundary edges, (cid:80) f Φ f is in general not outin (a) (b)(c) FIG. 1.
Networks with a uniform topology will havea gradient in face nutrient absorption density as thenutrient concentration decays. (a) Schematic of nutrientperfusion in a distribution network, with Φ f signifying thenutrient absorption density per face as delivered by adjacentedges. (b,c) Simulations of perfusion in networks with uni-form topology. All networks have a single flow source andsink, with Q in = Q out = 10, an initial nutrient concentration C = 100, and are constrained to lie within the 1 × f , indicated with color, for (b) a square gridwith 64 faces and (c) a Voronoi diagram with 50 faces. preserved for two arbitrary networks with the same inputnutrient concentration. The nutrient absorption densityover all faces is computed for two networks with a singlecurrent source and sink in Fig. 1(b) and (c), with colorindicating Φ f . A gradient in Φ f is clearly present, withfaces close to the source well-supplied while faces near thesink are starved. Our goal is to alter the geometry of thenetwork in a way that eliminates the nutrient densitygradient, yielding a uniformly perfusing network whichwill have all faces in the same color.This is accomplished by employing a face equaliza-tion algorithm which, similar to a vertex model, imposesforces on the network vertices, allowing their positions toshift. A vertex will experience a repulsive force from anadjacent face with a high nutrient density and an attrac-tive force from a low density face. Vertices lying on theboundary can shift as well, but the motion is restrictedto either purely in the horizontal or vertical directionso they remain on the boundary. The force on a vertexbecomes zero when all adjacent faces attain the same nu-trient absorption density. In this way, the high nutrientdensity faces grow in area and low density faces shrinkuntil perfusion is equalized across the network.We now outline the equalization procedure in detail.An important aspect of this algorithm is that the forceon a vertex is computed using only from faces contain-ing the vertex. Given a network, let { Φ f } be the set ofnutrient absorption densities over all faces f . For eachvertex i let (cid:104) Φ f (cid:105) i ≡ N f (cid:80) f,i ∈ f Φ f , the mean nutrient ab-sorption density of the N f adjacent faces to i . The forcefrom face h on i is set to have magnitude Φ h − (cid:104) Φ f (cid:105) i andis directed away from the face centroid along the angle FIG. 2.
The face equalization algorithm results in uni-formly perfusing networks.
Several sample Voronoi net-works with 50 faces are shown, with the corresponding equal-ized networks below. Equalized networks have a standard de-viation in { Φ f } equal to 0 . (cid:104) Φ f (cid:105) . The color bar and currentsource and sink locations are the same as in Fig. 1. bisector of the face, in a fashion consistent with osmoticpressure forces in vertex models [20]. Thus, a face h withΦ h > (cid:104) Φ f (cid:105) i exerts a force on vertex i pointing away fromthe face centroid (pushing the vertex away) and a face h with Φ h < (cid:104) Φ f (cid:105) i exerts a force towards the face centroid(pulling the vertex closer). The vertex coordinates areshifted by a fixed step size scaled by the net force fromall N f adjacent faces. If all faces adjacent to a vertexhave equal nutrient density, Φ h = (cid:104) Φ f (cid:105) i for all h and thenet force on the vertex will be zero. After vertex coor-dinates shift, the nutrient flow and absorption throughthe network are recalculated, and the new set { Φ f } isrecomputed. The process repeats over all vertices, stop-ping when the standard deviation of { Φ f } is less than onepercent of (cid:104) Φ f (cid:105) , averaged over all faces. The fixed stepsize (set to be 10 − ) is chosen to ensure that no edgesoverlap once all vertex coordinates are shifted. If this issatisfied, we find the algorithm to be robust and capableof achieving arbitrarily uniform networks.Ideally, this adaptation process is purely geometric,preserving the set of vertices and edges and changing onlythe vertex positions. However, problems can arise whennetwork edges overlap. To avoid spurious edge crossings,we allow two ways for the network topology to change.First, if an edge becomes shorter than a threshold length(set to 1 /
30 of the full network side length) then the edgeis removed and the two vertices connected by the edgeare merged. Second, if an angle formed by vertices abc becomes smaller than a threshold angle (set to 5 degrees),then the longer edge ab is removed and a new edge ac isadded, effectively collapsing the acute angle.To analyze the effects of equalization, we use a set of50 Voronoi networks with 50 faces, with initial vertexpositions chosen at random and adjusted using a springforce-directed positioning method. These networks areequalized to produce an ensemble of uniformly perfusingnetworks, a sample of which is shown in Fig. 2 along-side the initial networks. Figure 3(a)-(c) outlines thegeneral motion of network faces during the equalization (c) (d)(b)(a) (b) FIG. 3.
Face centroids shift towards the sink as theface absorption density is equalized.
Information from50 networks is spatially binned, showing trends across the ini-tial and equalized network ensembles. (a) The initial Voronoinetworks have the density Φ f decaying away from the source,but (b) a uniform distribution of face centroids. Arrows in-dicate the average movement of the face centroids during theequalization procedure, with arrow length proportional to themean displacement, showing a general motion of faces towardsthe sink. (c) Equalized networks have a higher density of facesby the sink. (d) Φ f serves the role of pressure in the usual ver-tex model, as shown by the trajectories of face areas over thecourse of the equalization for a single network. Faces grow orshrink faster in the early stages, but as networks are equalizedthe change in area for all faces approaches zero. algorithm, ultimately showing the trade-off between uni-formly distributed face centroid positions and uniformlydistributed nutrients. These figures show the cumulativevalue of either all networks face centroid positions or faceΦ f values over all 50 networks, binned in a 10 ×
10 gridindicating the space of the network. For the initial net-works, nutrients are concentrated by the source Fig. 3(a),showing once again a gradient in nutrient absorption den-sity. In the initial networks the face centroid distributionis approximately constant across space of the network,shown in Fig. 3(b). The motion of the face centroidsduring equalization is shown here as well, with arrows in-dicating the mean direction of motion for a single spatialbin. Thus, the equalized networks have a higher densityof face centroids by the nutrient sink, shown in Fig. 3(c),but a uniform distribution of Φ f density (not shown).This indicates that equalization under the local adaptiverule serves to exchange uniformity in face centroids foruniformity in face nutrient absorption density by shiftingthe network face centroids towards the sink.Finally, Fig. 3(d) shows how the incremental force-equalization algorithm functions similarly to a growthinduced pressure, as vertex forces arise from nutrient den-sity differences in adjacent faces. We take G f to representa quantity akin to the cumulative force a face is subjectedto by difference between its own nutrient absorption den-sity and that of its neighbors: G f = (cid:80) g (Φ f − Φ g ) ˜ L fg ,where ˜ L fg is the length of the edge separating faces f and g . We find a nearly linear relation between G f for aface f and its change in area over time, dA f /dt . Grow-ing faces start with dA f /dt > dA f /dt <
0, and as the network is equalized dA f /dt approaches 0 for all faces.We now quantify morphological network features thatcause equalized networks to be uniformly perfusing,namely the distribution of the face sizes and shapes inthe two sets of networks. We find that uniform networkshave larger faces near the source and smaller faces nearthe sink. This is demonstrated in Fig. 4 (a)-(b), whereface area is plotted as a function of Euclidean distancefrom the face centroid to the source at (0 , p = P/ √ A where P and A are the face perimeter and area. Regularpolygons have minimal shape parameter; for example, anequilateral triangle has p = 2( √ ≈ .
56 and a squarehas p = 4, whereas a non-regular triangle or quadrilat-eral will have a strictly larger p .The change in the face shape parameter distributionbetween the initial and equalized networks is shown inFig. 4(c)-(d). For the initial networks, p has no correla-tion with the location in the network, but is dictated bythe number of sides in the face. Moreover, all faces arenearly regular, since the shape parameter for each typeof polygonal face lies close to the minimal shape parame-ter for that regular polygon. For the uniformly perfusingnetworks, there is no longer a strong dependence betweenthe face p and the number of sides, but there is a corre-lation between the shape parameter and location: facesnear the source tend to be more compact and faces by thesink tend to be more elongated. This can be explained faces colored bynumber of sides FIG. 4.