Optimal Transport Flows for Distributed Production Networks

Network flows often exhibit a hierarchical tree-like structure that can be attributed to the minimisation of dissipation. The common feature of such systems is a single source and multiple sinks (or vice versa). In contrast, here we study networks with only a single source and sink. These systems can arise from secondary purposes of the networks, such as blood sugar regulation through insulin production. Minimisation of dissipation in these systems lead to trivial behaviour. We show instead how optimising the transport time yields network topologies that match those observed in the insulin-producing pancreatic islets. These are patterns of periphery-to-center and center-to-periphery flows. The obtained flow networks are broadly independent of how the flow velocity depends on the flow flux, but continuous and discontinuous phase transitions appear at extreme flux dependencies. Lastly, we show how constraints on flows can lead to buckling of the branches of the network, a feature that is also observed in pancreatic islets.

Network flows often exhibit a hierarchical tree-like structure that can be attributed to the minimisation of dissipation.The common feature of such systems is a single source and multiple sinks (or vice versa).In contrast, here we study networks with only a single source and sink.These systems can arise from secondary purposes of the networks, such as blood sugar regulation through insulin production.Minimisation of dissipation in these systems lead to trivial behaviour.We show instead how optimising the transport time yields network topologies that match those observed in the insulin-producing pancreatic islets.These are patterns of periphery-to-center and centerto-periphery flows.The obtained flow networks are broadly independent of how the flow velocity depends on the flow flux, but continuous and discontinuous phase transitions appear at extreme flux dependencies.Lastly, we show how constraints on flows can lead to buckling of the branches of the network, a feature that is also observed in pancreatic islets.
Transport networks are essential for life to function on large multicellular scales.In vertebrates, blood flow delivers energy and nutrients, and removes waste through the branched network of the vascular system.The separation of vessels into arteries and veins ensures that oxygen-rich and oxygen-depleted parts of the network are kept separate.In plants the separation into xylem and phloem provides a similar function.In nature, these systems typically exhibit tree-like hierarchical structures, as e.g. in the aorta which splits into increasingly smaller arteries all the way to capillaries, the smallest vessels of the circulatory system.The tree structure can be understood as minimising the dissipation of the blood flow through the system [1][2][3].Furthermore, loop-redundant tree-like structures, as is evident from observing e.g. the veins of a tree leaf, is explained by minimising dissipation under damage or under fluctuating needs [4,5].Perhaps because of its origin as a power-minimising network, treelike structures are observed not just in vascular networks, but also for instance in both natural and artificial river networks [4,6].The shared feature of these systems is that they consists of a single fluid source with a lot of sinks, or a single sink with a lot of sources (these are equivalent, dual formulations).For instance, in the arterial system the heart is the source and the body cells the sinks, whereas the roles are reversed in venous systems.
The study of flow networks that have a single source and a single sink have received much less attention, perhaps because it is not clear what properties can be optimised over such networks -indeed dissipation minimisation lead to trivial behaviour.However, such networks are critical in flow systems that have secondary products (e.g.other than oxygen in blood) that are being produced or delivered along the flow path.Here, we exemplify such a system by the Islets of Langerhans in the pancreas.In these islets, beta cells release insulin and alpha cells glucagon into the blood stream based on blood glucose levels [7,8].In such systems there is no need for an arterial-venous separation, and the transport from any hormone-producing cell directly couples to both the downstream and upstream transport systems.
The vasculature of pancreatic islets differs from species to species.In particular, various topologies have been observed: periphery-to-center flow, straight through, and center-to-periphery flow, as illustrated in Fig. 1(a-c) [10][11][12].For instance in rodents, the center-to-periphery topology is the most common [11].furthermore, the vasculature of these islets is often very tortuous compared to the vasculature of other organs [9, 13] as shown in Fig. 1d.In this Letter we propose an optimisation problem over single-source, single-sink flows and use this to provide possible explanations for the type of features observed in pancreatic islets.
Model -To study systems of blood flow optimisation we consider, as in previous studies [1][2][3][4][5], flows on networks.While pancreatic islets indeed can have more than one inlet and outlet, we simplify the system and consider the network shown in Fig. 1e.Nevertheless, our approach works for any number of sources and sinks.In Fig. 1e cells are represented by hexagons, and edges betweens between these cells indicate where fluid can potentially flow.Inlet and outlet are indicated by arrows.This specific graph has n N = 130 nodes and n E = 357 edges, and each edge e of the graph has a conductivity C e and a length L e .An n E × n N oriented incident matrix ∆ of the graph gives each edge a unique direction, and we can thus attach to each edge a (signed) flow F e .We furthermore define the source vector S with S source = 1, S sink = −1 and S i = 0 elsewhere and require that the flow obeys Kirchoff's current law Since n N < n E , the flow is far from determined by this condition alone.
To make the flow well-defined, we furthermore require that it is derivable from a potential based on the effective conductivities, where p i is the potential (pressure) defined at the nodes, and C eff is a diagonal matrix with entries C eff ee = C e /L e .Combining Eq. ( 1) and ( 2) we can solve for the potential as where † denotes the pseudo-inverse, which is needed since the system of equations is singular, but can be solved by the pseudo-inverse if i S i = 0, which is indeed the case here.From the potentials p, the flows are immediately obtained from Eq. ( 2).The total power dissipation of the system is P = e F 2 e /C eff ee [1][2][3].As mentioned, minimising this term leads to tree topologies for a single source and sinks everywhere (S source = 1, S i = −1/(n N − 1) elsewhere).This optimum on our network topology is shown in Fig. 2a.The minimisation is done under constant "material cost" e L e C γ e , and the tree-like structures are obtained for γ < 1.In the optimum, the conductivities scale with the flow as [14] i.e. a larger conductivity is needed where there is a lot of flow.
Trivially, minimising power dissipation in the network with just a single source and single sink leads to flow only going along the shortest path between the source and sink as shown in Fig. 2b.This is also, naturally, the time minimising network for flow between the source and the sink.We are interested in network structures that visit all nodes in an "optimal" way.Indeed proper distribution of vessels is far more important in pancreatic islets and other systems than the (potentially miniscule) power being dissipated.
We consider instead graphs that minimise the average time for the product (e.g.insulin) being produced at the nodes to reach the outlet (the opposite inlet-centric definition will be discussed later).This time-optimised graph we define as follows: Denote for each node T i the average time taken from that node for product to reach the outlet.The average time that we intend to minimise is thus where T i is found by the recursive relation which follows by letting product flow in proportion to the fluid flow, with the special case T sink = 0.This is a linear equation for T i , where O i is the set of nodes that are outgoing from node i.Whether an edge e is outgoing from node i can be identified by the criteria F e ∆ ei > 0. T ij is the time taken for the product to flow from node i to neighbouring node j.The physics of the system is determined through the relation between T ij and F , C and L.
For blood vessels each edge corresponds to a tube of a given radius r e .We can thus relate F e ∼ v e r 2 e , where v e is the fluid velocity along the tube.Furthermore, in Poiseuille flow C e ∼ r 4 e , and thus √ C e /F e .This shows that the time travelling across an edge is minimised by having a small conductivity, since in small tubes the liquid will have to move faster for the same flux F .So while we are interested in solutions that minimise time, this formulation yields solutions that are severely dissipation inefficient in which case power optimisation becomes more relevant than time minimisation.Instead, we reformulate T ij in terms of only the flows F such that we stay in reasonably power-efficient regimes, as established by Eq. ( 4).Taking γ = 1, we have Periphery-center optima -The result of our minimisation scheme is shown in Fig. 2c.First we note that the single-source-sink system prevents self-similar branching solutions that are known from single-source, multiplesink systems.Instead the solution only has a few main branches, in particular one at the periphery and one at the center.Interestingly, this exact pattern is one of the three topologies [Fig.1(a-c)] of pancreatic islet blood flow observed in nature [11,12].
In some species, such as rodents, the glucagonproducing alpha cells and the insulin-producing beta cells comprising the islets are heterogeneously distributed with the beta cells in the center and alpha cells at the periphery.It has thus been suggested that the order of the flow suggests intercellular communication and regulation [11], i.e. beta cells regulating alpha cells or vice versa, depending on the flow topology.Our results show a separate, but non-exclusive explanation that the patterns can appear due to an optimisation of the flow itself, independent of any heterogeneous distribution of cells.Interpreted in an evolutionary sense, one could imagine that optima such as the one derived here led to the present flow topology before the advent of alpha and beta cells and thus drove the heterogeneous distribution thereof.The equal distribution of flow to each node and the sparsity of the network in the optimum [Fig.2c] is possible due to the scaling of Eq. ( 7).In fact, considering instead T ij ∼ L ij /|F ij | δ our model works for a broad range of δ.To illustrate this, consider the simple network of three nodes shown in Fig. 3a.
1 , which is shown in Fig. 3b for various values of δ.Minimising T , the optimal F 1 is shown in Fig. 3c, which demonstrates that this system has a discontinuous phase transition at δ ≈ 0.275 and a continuous phase transition at δ = 1 in F 1 .Between these values the optimal F 1 is independent of δ, which is the regime we study.
The discontinuous phase transition occurs because the left and top node of Fig. 3a have conflicting optima: the left node minimises its time by having F 2 large, whereas the top node needs F 1 large.Large δ favour large F 1 because then a larger flow velocity compensates for the larger length of the upper branch.In contrast, a small δ favours the shorter path of the lower branch, while only leaving a small flux through the upper branch required to transport the product from the top node.
For more complex graphs such as the one we consider [Fig.1e], the phase transition behaviour is much more rich and complex and the discontinuous transition is pushed to lower values, and the continuous transition happens at a δ smaller than (but close to) 1.Our choice of δ = 1 /2 lies safely within the regime, where the results are independent of the precise value of δ and where the resulting networks are sparse graphs.Thus our results are to a large degree independent of the scaling in Eq. (7).
The left-right asymmetry and thus the periphery-tocenter flow of Fig. 2c comes from the fact that we are minimising the time for the product to reach the outlet from the nodes.The large collection branch emerging from the center thus minimises the time for many nodes by providing a fast route.Had we instead minimised the time from the source to the nodes T r the solution would be reversed, since it would be important to reach the nodes fast.Indeed this opposite situation, with flow from center to periphery, is the most commonly observed topology in rodents [11], and could perhaps hint that the time for "information" (e.g.glucose levels) to reach the cells is more important than the product to exit the islet.
Naturally, the combination of time from inlet to nodes and from nodes to outlet T c = (1 − α) T + α T r (0 ≤ α ≤ 1) can also be optimised over.Fig. 4 shows examples of minimising T c .α = 0 corresponds to situation we have already studied, and α = 1 simply left-right mirrors this solution.In-between a compromise of the two optima is reached, e.g. as shown for α = 0.25 in Fig. 4a.The value of T c is minimised for α = 0 and α = 1, since for any other value the solution will be the average of two solutions that both compromise.At precisely α = 0. an islet topology similar to that of Fig. 1b.
Minimal flow constraints -While the solution of Fig. 2c does indeed distribute the flow to all nodes, some nodes see much more fluid flowing through them than others.As the pancreatic islets grow and if the blood flow source cannot keep up with this growth, each cell/node will see a decrease in flow permeating them.In this way, growth in networks flows can be emulated by a decreasing S source [3].At the same time a minimum amount of flow might be required at each cell e.g. for accurate estimation of glucose levels.This leads us to consider adding minimal flow constraints to the optimisation problem.
The flow through a node i with S i = 0 is and we now require F ≥ F m for all nodes for a given value of F m .Obtained minima are shown in Fig. 5 under variations of F m .The optimal configuration of Fig. 2c has min i F = 0.041 and thus for any F m smaller than this the same solution is obtained (blue section in Fig. 5f).As F m is increased above this, the network adapts to more equally divide the flow.First the two middle branches are lost [Fig.5a] and then the collection point starts moving towards the outlet [Fig.5b].In the end of this process (green section in Fig. 5f), the collection branch has moved all the way to the right [Fig.5c].As can be seen in Fig. 5f, naturally, the average time T increases as F m increases.Fig. 5f also shows that F i decreases rapidly.So as this reordering occurs, the minimum flow rate increases at the expense of the average flow rate.
As F m is increased further than reordering can accomodate for, buckling occurs [Fig.5d], the degree of which increases as F m is increased [Fig.5e].As soon as buckling occurs, F starts increasing as seen in Fig. 5f.This increase scales linearly with F m , i.e. the average flow rate, after buckling, stays at a fixed level above F m .We note that the noise in F in Fig. 5f most likely indicates that our optimisation scheme in some cases fails to find the true global optimum but ends in neighbouring local minima.
As shown in Fig. 1d, pancreatic islets do indeed have a severely tortuous vasculature [9] similar to that obtained in Fig. 5.The mechanical reason for buckling in real islets could be due to growth-induced buckling [13].Our analysis shows how such buckling could actually be of benefit to the system under growth.In comparing our result with blood vessels in pancreatic islets, we implicitly assume that transport time is of main concern.In fact the average blood flow velocity in Langerhans islets is ∼ 1.4 mm/s [9], implying that transport across an islet takes about ∼ 0.5 second.Activity of both alpha and beta cells are pulsatile, and in-vitro experiments show coherent oscillations of whole islets with periods down to three seconds [17].In order for the organism to utilize coherent release of insulin downstream of islets, it is therefore plausible that time optimization on the sub-second scale is functional.
In conclusion, by simply varying how much time to nodes and time from nodes is important through the parameter α, our proposed optimisation problem on network flows directly leads to flow topologies similar to those observed in real pancreatic islet and illustrated in Fig. 1(a-c).We have furthermore shown how minimum flow constraints leads to buckling and thus tortuous vessels as found in real islets.Our approach applies generally to all laminar network flows that have single-source and single-sink characteristics.
This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation Programme, = 0 from which the scaling follows.See e.g.Ref. [1] for details.

FIG. 1 .
FIG. 1. Vasculature of pancreatic islets.(a) Center-toperiphery flow.(b) Left-right symmetric flow.(c) Peripheryto-center flow.(d) Microscopy of vasculature, from Ref. [9], showing the high degree of tortuousness in the vasculature of pancreatic islets.(e) The network structure considered.This is constructed by taking a circular subsection of a hexagonal grid of nodes.Neighbours are then induced from a Delaunay triangulation.Each edge has a conductivity Ce associated to it, which are the parameters to be optimised over.Inlet and outlet are shows by arrows.

FIG. 2 .
FIG. 2. Optimised network structures.(a) Minimal dissipation network with sinks at all nodes.(b) Minimal dissipation network with a single sink at edge.(c) Per-node time minising network with a single sink at edge.Cell colours indicate product concentration flowing through that node as calculated by solving ρi = 1 + j∈I i |Fij| ρj/ k∈O j |F jk | (Ii denoting incoming and Oi outgoing nodes of node i).Flow lines are coloured by pressure, their thickness indicating (square root of) flow magnitude.
which is the definition we will use, and critically what enables the optimisation problem to yield non-trivial results.Variations induced by different values of γ also work, as will be shown.We thus minimise the average time T in the regime where large fluid flows are associated with large conductivities.Equations (5-7) define the optimisation problem, which we solve by a momentum-based version of gradient descent [15].

FIG. 4 .
FIG. 4. Varying the ratio α between time from nodes to outlet and inlet to nodes.(a) Obtained optimum at α = 0.25.(b) The left-right symmetric solution obtained at α = 0.5.

FIG. 5 .
FIG. 5. Optimal graphs under minimal flow constraints.(a) At Fm = 0.042, the network of Fig. 2c loses its front collection channels.(b) The collection point starts moving to the right as Fm is increased to 0.062.(c) At Fm = 0.073 the collection point has moved to the far right.(d) One horizontal vessel has been removed by introducing kinks at F = 0.078 which increases the flow in the remaining vessels.(e) Buckling increases; Fm = 0.125.As Fm → 1 the global optimum becomes a Hamiltonian path through the nodes.(f) In dark blue (left axis) the average time T is shown as function of Fm.The average flow through the nodes F is shown in light brown (right axis).Dashed line shows f (x) = x, which matches the slope of F in the buckling section.Background colours denote section: global optimum in blue, unbuckled solution in green, and buckled in red.