Discontinuous transition to loop formation in optimal supply networks

The structure and design of optimal supply networks is an important topic in complex networks research. A fundamental trait of natural and man-made networks is the emergence of loops and the trade-off governing their formation: adding redundant edges to supply networks is costly, yet beneficial for resilience. Loops typically form when costs for new edges are small or inputs uncertain. Here, we shed further light on the transition to loop formation. We demonstrate that loops emerge discontinuously when decreasing the costs for new edges for both an edge-damage model and a fluctuating sink model. Mathematically, new loops are shown to form through a saddle-node bifurcation. Our analysis allows to heuristically predict the location and cost where the first loop emerges. Finally, we unveil an intimate relationship among betweenness measures and optimal tree networks. Our results can be used to understand the evolution of loop formation in real-world biological networks.

T he reliable function of supply networks is essential for biological as well as technical systems. Leaf venation networks supply plant leaves with water and nutrients 1 and vascular systems supply vertebrates with oxygen and nutrients 2 . On the other hand, society relies on man-made supply networks such as power grids 3 or hydraulic networks 4 . Finally, networks that formed over time such as drainage basins show a similar structure 5 . Understanding the design principles of such networks is a central challenge in network science 6 .
The evolution or construction of supply and transportation networks is essentially determined by the trade-off between cost and resilience [7][8][9] . Cost limits the number of connections in the network, as resources are generally scarce. Resilience requires additional connections to cope with damages or perturbations. Many actual networks contain loops to establish a certain level of topological resilience, hence they stay connected and operational even if some elements fail 10 . The interplay of topology and resilience is analysed in various disciplines including traffic networks 8 , communication networks 11 or dynamical networks 12 . Finally, a variety of results on structural resilience, that is the ability of a network to remain connected when a fraction of nodes or links fails, have been obtained in network science 13,14 .
In this article, we focus on linear flow networks modelling power grids, hydraulic networks or vascular networks 3,4,15 . Different structural patterns are observed in nature, consisting of both networks with and without loops. For instance, leaf venation networks are loopy in general, except for a few old species such as Ginkgo. In electric power systems, large-scale transmission grids are strongly meshed, while local distribution grids are topological trees (Fig. 1). Optimal network structures balancing costs and resilience have been analysed via extensive numerical simulations in the setting where a single source supplies the remaining network, such as in plant leaves [15][16][17] . The optimal structure does not contain any loops if connections are reliable and perturbations are weak, for instance in distribution grids. Loops come into being when sources or sinks fluctuate strongly or connections are subject to damages, such as in transmission grids or newer leaf species. While some work has been done in the context of networks optimising transport time 18 , the exact mechanism of loop formation in minimaldissipation networks is still not fully understood.
Here, we analyse the transition to loop formation on a theoretical basis and derive several analytical results. We consider optimal network structures in the sense that function is optimised while costs are constrained or vice versa. Two aspects of resilience are studied in detail-damage to edges and fluctuations of supply and demand. In particular, we investigate the optimal structure as a function of the severity of damage and the strength of fluctuations. In contrast to prior work, we focus on the occurrence of the very first loop, which enables an analytical approach to loop formation and yields several rigorous results. We first establish this approach for an elementary sample network and then generalise it to networks of arbitrary size and compare analytic predictions and numerical results.
In particular, we demonstrate that the transition to loop formation is generally discontinuous in the sense that optimal edgecapacities jump discontinuously when fluctuations increase or costs decrease. Loopy network structures emerge as new local minima of the dissipation function that form via a saddle-node bifurcation, and not via a bifurcation of an already existing minimum. Hence, a large number of local minima may exist simultaneously and we establish a purely topological expression based on the edge betweenness to understand their structure. As a direct application of our analysis, we derive a simple criterion to predict the location of the first loop in the transition from a tree network.

Results
Modelling supply networks. We consider a simple supply network model which was previously used to study loop formation in generic distribution networks 15,16 . Mathematically, the supply network is constructed from a graph G with node set V and edge set E. At each node n ∈ V, there is an in-or outflow with a strength P n , where P n > 0 denotes a source and P n < 0 a sink. The in-and outflows may either represent individual supply nodes or allocated demands associated with the node 19 . An edge in the network is either labelled by its index e ∈ E or by its terminal a b c d nodes e = (n, m) which we use interchangeably. For each edge, we fix an orientation which is encoded in the node-edge incidence matrix I 2 R jVj jEj with elements I n;e ¼ 1 if edge e starts at node n; À1 if edge e ends at node n; 0 otherwise: Each edge is assigned a capacity k e 2 R þ and a flow whose strength or value is denoted as F e 2 R. Fixing the orientation of an edge e = (n, m) means that F e > 0 describes a flow from node n to node m and F e < 0 describes a flow from node m to node n. The flows satisfy the continuity equation or Kirchhoff's current law (KCL) at every node of the network, X e2E I n;e F e ¼ P n ; 8 n 2 V: ð2Þ In addition to that, we assign a potential θ n to each node in the network. In terms of physical quantities, this potential θ n 2 R can represent the pressure at the nodes of a hydraulic network, the voltage in DC resistor networks or the nodal voltage phase angle in linearised AC power grids 4,15,20,21 . For these systems, the flow on a link e = (n, m) scales linearly with the potential drop θ n − θ m along the link and can be calculated as Together with the continuity equation (2), this linear set of equations determines the values of the potentials θ n up to a global constant. The resulting flows automatically satisfy Kirchhoff's voltage law (KVL) which states that the flow around any closed loop expressed in terms of the edges C ¼ fe c 1 ; e c 2 ; :::; e c max g vanishes [ref. 22 Here, the factor z e ∈ { − 1, 1} is used to keep track of the orientation of an edge e with respect to the orientation of the edge in the loop C, i.e. z e ¼ 1 if edge e ¼ ðe 1 ; e 2 Þ is oriented from e 1 to e 2 ; À1 if edge e ¼ ðe 1 ; e 2 Þ is oriented from e 2 to e 1 : & Optimising supply networks: minimum dissipation topologies. We now illustrate how to determine the optimal supply network that is described by the above the set of equations. To this end, we want to find the edge capacities that determine the network structure that is optimal for performing a given task. Throughout this manuscript, we call the network structure optimal if the edge capacities are such that the overall network dissipation is minimised, as suggested for example in refs. [15][16][17] . The network dissipation may be calculated as In addition to that, we assume that the resources to build the network are limited. This resource constraint takes the form X where the cost parameter γ > 0 depends on the type of problem under consideration. For instance, assuming Poiseuille flow through cylindrical pipes of fixed length and radius R e , k e $ R 4 e , such that γ = 1/2 fixes total fluid volume and γ = 1/4 fixes total pipe mass [15][16][17]23,24 . The parameter K corresponds to the overall available budget. Note that different definitions of optimal networks arise in other applications, e.g. in hydraulic engineering where typically the cost is minimised while the dissipation is constrained 25 . In Supplementary Fig. 5 and Supplementary Note 5, we demonstrate that the same kind of discontinuous transition is observed when extending our analysis to this setup.
In general, it is neither useful nor meaningful to allow arbitrary connections between the nodes. Geometric constraints apply to a variety of networks. For instance, leaf vascular networks or river basins are naturally planar. To take into account such constraints and keep the problem feasible one typically fixes a set of potential edges E such that E E. These edges are often taken from a square grid 16 , a triangular grid 15 , or various types of disordered tessellations 7,26 . Note that while planarity of the network described by the set of potential edges E simplifies the theoretical analysis, our results are not limited to planar networks as we demonstrate for a simple, non-planar network in Supplementary  Fig. 6.
We focus on two different models here: a model with fluctuating sources and sinks and a model of stochastic damage to the edges. Both models can be thought of as quantifying network resilience: We call a network resilient if it is able to function properly under the uncertainties induced by edge damage or fluctuating inputs. For both models, our main question will be the following: Under which conditions does the optimal network structure contain loops and how do these loops emerge?
Fluctuating sink model: First, we introduce the fluctuating sink model. In this model, we include fluctuations by treating the P n as random variables. For each random realisation, the sources and sinks are balanced, i.e. they sum to zero, X n2V P n ¼ 0: Network structures are then optimised to have a minimum average dissipation for a given set of resources. Here, the brackets 〈⋅〉 denote the expected value taken over all realisations of the random variables P n . Note that the fluctuations affect only the flows directly by virtue of Eq. (3), whereas the network topology is assumed to be fixed by the construction of the network such that the average is taken over the squared flows only. Equation (8) can be minimised analytically with respect to the k e , where the resource constraint is taken into account via the method of Lagrange multipliers. Calculating the optimal edge capacities by extremising the Lagrange function yields 16 (Supplementary Note 1) This expression depends on the second moments of the flows hF 2 e i, which in turn depend on the capacities k e . Hence, Eq. (9) can be interpreted as a self-consistency condition which has to be solved together with Eq. (3).
Edge-damage model. A second class of dissipation-optimised networks that is relevant to biology and engineering seeks to find optimal networks subject to damage. For instance, leaf vasculature might be attacked by a herbivorous insect, or a power grid might lose a power line due to an outage. In the following, we generalise the broken-bond model considered in ref. 15 by allowing partial damage to the network capacities instead of complete removal of edges.
In this edge-damage model, the sources and sinks are still balanced but do not fluctuate stochastically. Instead, we assume that all nodes but one are sinks with P j > 1 ¼ À P supplied by a single node with P 1 ¼ ðN À 1Þ P, where N is the number of nodes. To model partial damage of edge l, we modify the edge capacities according to with the damage fraction Thus, a damage parameter Δ = 1 corresponds to complete removal of the damaged edge. We now define the average over all possible damage scenarios. Specifically, if g(k e ) is some function of the capacities k e , we define where |E| is the number of edges in the network. Here and in the following, we use the notation hÁi 0 to distinguish the average over damage scenarios from the average over fluctuating sources and sinks. As before, the central objective is to minimise the average dissipation of the network, taken over all possible damaged edges under the resource constraint Eq. (6). We now proceed to study loop formation in the two models outlined above in detail.

Discontinuous transition to loop formation in small network.
As an illustrative example, let us consider an elementary network as sketched in Fig. 2a and analyse the transition to loop formation in both, the fluctuating sink model and the edge-damage model.
Disontinous transition in fluctuating sink model: The network consists of four variable sinks at nodes 2, 3, 4, 5 (circles) that are modelled as iid Gaussian random variables P 2;3;4;5 $ N ðμ; σÞ and four edges (arrows) connecting them with capacities k i and flows F i , i ∈ {1, 2, 3, 4}. A fifth, potential edge is shown as a dotted arrow. If it exists, it carries flowF 5 and has capacity k 5 = κ (Fig. 2a). The central question we will study for this setup is the following: When is the optimal network tree-like (κ = 0) and when is it loopy (κ > 0)-and how does κ behave at the transition point?
We first consider the case where the loop is not present, i.e. κ = 0. In this case, the network is a tree and we can calculate the second moments hF 2 i i; i 2 f1; 2; 3; 4g explicitly in terms of the capacities: they are determined by the statistics of the source and the sinks by virtue of the continuity equation (2). Using the optimal capacities for a tree network (6), we obtain an explicit equation for the optimal dissipation 〈D tree 〉 that only depends on the statistics of the sinks (Supplementary Note 4) How does this result change if we allow closing the loop as illustrated in Fig. 2a, i.e. if we include the corresponding edge in the set of potential edges E?
Let us assume that the loop carries a flowF 5 and has a nonzero capacity k 5 = κ > 0. In the following, we denote the flows and capacities in the loopy network with a tilde. In the presence of a loop, we can no longer determine the flows using the continuity equation (2) alone. Instead, we have an additional degree of freedom: a cycle flow f around the newly formed edge such that This approach allows us to eliminate the dependence on the cycle flow strength f, and we can evaluate the dissipation 〈D loopy 〉 of the loopy network by inserting the result into Eq. (8) (Supplementary Note 4). The new expression for the dissipation no longer contains the flows explicitly, which considerably simplifies finding the optimal topology: we no longer have to take care of the interdependence of flows and capacities, but can minimise 〈D loopy 〉 in terms of only the capacitiesk i .
We proceed to evaluate the optimal network structure fixing the mean of fluctuations to μ = −1 and the resource constraint to K = 1. To examine the effect of the two remaining parameters separately, we analyse the transition to loop formation for γ = 0.9 fixed while varying σ and for σ = 3 fixed with varying γ. We then compute the dissipation 〈D loopy 〉 as a function of the capacities κ and k 1 and compare it to the dissipation 〈D tree 〉 of the corresponding tree network. Note that the capacities in the optimum tree network are explicitly given by Eq. (9) such that 〈D tree 〉 is fixed. For the loopy network, we still need to determine the optimum structure, i.e. we compute the minima of 〈D loopy 〉 as a b c Graph set up to analyse the transition from tree networks to loopy networks. a Elementary network to study spontaneous loop formation in optimum supply networks. The network consists of five nodes (green circles) where node n = 1 has an inflow of four, P 1 = 4, and all other nodes have an outflow of unity. These in and outputs determine the flows F i , i ∈ {1, 2, 3, 4, 5} along with the links with capacities k i . The optimum topology for this set-up is a tree network. If the in-and outputs are fluctuating, an additional edge (dotted arrow) may be beneficial to reduce the average dissipation. This edge introduces a new degree of freedom expressed as a cycle flow f. b For a larger network, we generalise this setup as follows: we start from a tree network and then consider the impact of a new edge at an arbitrary position (n, m) (dotted, red arrow). We then collect the edge sets L (shaded green) and R (shaded blue) along the shortest path from the source to the newly formed edge. This edge induces a cycle flow f. c A network formed from a triangular grid with a set of potential edges E coloured in grey which we will analyse throughout the manuscript. Realised edges (black) correspond to a global minimum of the dissipation for the fluctuating sink model where a single, fluctuating source (large circle) supplies the remaining network.
a function of κ andk 1 recalling thatk 3 ¼k 1 ,k 4 ¼k 2 andk 2 is then fixed by the resource constraint Eq. (6). For both varying fluctuations σ and varying costs γ, we find that the transition to loop formation is discontinuous: the loop starts to form with a non-zero capacity κ when analysing the globally optimal network structure ( Fig. 3a, c, thick, orange line). Analogously, the capacity k 1 bifurcates (red line).
But what is the nature of this transition? In fact, we find that new minima emerge through a saddle-node bifurcation independently of the parameter we vary. Thus, new minima do not form from the existing tree minimum but instead emerge elsewhere in the energy landscape. To support this claim, we plot the capacity at the saddle in Fig. 3 (dotted, coloured lines) and analyse the dissipation landscape close to the bifurcation point ( Supplementary Fig. 3). Using these results, we can also map out the parameter region where loop formation is beneficial ( Supplementary Fig. 2). In Supplementary Fig. 7, we illustrate the nature of this transition for an even simpler network and find a closed-form solution for the region of the parameter space where loop formation is beneficial.
Discontinuous transition in edge-damage model: We now turn to the edge-damage model and analyse the optimal topology again for the graph shown in Fig. 2a. Most importantly, we find that the transition between a tree-like and a loopy optimal network is also discontinuous in the damage model in both the cost parameter γ and the damage parameter Δ, and new extrema appear again through saddle-node bifurcations (Fig. 3b, d). This demonstrates that despite the fact that in the damage model, the optima follow a different scaling law from those in the fluctuation model 15 , the mechanism and type of the transition from tree-like to loopy optimum is generic.

Discontinuous transition persists beyond the first loop.
Whereas the transition to the first loop that forms is important in many real-world supply networks, such as the Gingko leaf and the distribution grid shown in Fig. 1, other networks display several loops, such that their formation beyond the first loop becomes important. In particular, the tree has mainly theoretical importance in many applications such as hydraulic networks where spanning trees in loopy networks play an important role in modelling and optimisation 25,[27][28][29] . Remarkably, we can demonstrate numerically that the discontinuous character of loop formation persists beyond the first loop.
In Fig. 4, we analyse this transition for the fluctuating sink model with cost parameter γ = 0.5 for a larger, globally optimal tree network which was formed from a set of potential edges E corresponding to a triangular grid as shown in Fig. 2c. We map out the order in which new loops form (colour code) when decreasing the cost for new edges and slightly perturbing the previous network structure. All new loops emerge discontinuously with a non-zero capacity from an existing loopy network topology (Fig. 4c). Note that in contrast to Fig. 3, the optimal capacities are obtained here using an iterative approach for finding local minima of the dissipation that is due to ref. 16 (see "Methods" section). In a Supplementary Fig. 8, we demonstrate that an analogous transition exists for varying fluctuation strength σ and fixed cost parameter γ.
Identifying optimal trees for networks of arbitrary size. We now generalise our reasoning to larger networks with an arbitrary number of nodes N. For this analysis, we focus on the fluctuating sink model. Again we assume that all nodes j = 2, …, N act as sinks with P j being random variables and that the source j = 1 balances the sinks. We start from a tree network and analyse at what value of the cost parameter γ it becomes beneficial to add a single edge thus closing a single loop. This setup is sketched in Fig. 2b. We first demonstrate how to calculate the dissipation in such a setting and then illustrate the procedure to minimise it.
In an arbitrary tree network, the flows do not depend on the link capacities but only on the topology of the network as illustrated in the last section. This is due to the fact that for each node j = 2, 3, … there is only one path from the respective node to the root j = 1 of the tree. The flow F e on an edge e is thus directly given by the KCL Eq. (2). Here, we fix the orientation of the flows such that they point away from the source as illustrated in Fig. 2b. Therefore, flows in tree networks are always positive.
To express the flows F e in terms of the sources and sinks P j , we introduce the tree matrix T 2 R jEj jEj by T e;j ¼ þ1 if edge e is on the path from node j þ 1 to the root j ¼ 1 0 otherwise: This yields an explicit expression for the flows, We can insert this result into the network dissipation (Eq. (8)), which yields where T = E(G) is the set of all edges in the tree, i.e. before the addition of a loop.
From trees to loopy networks: optimising networks with a single loop. Remarkably, we can also find an explicit expression for the dissipation eliminating the flows in a near-tree network by exploiting the KVL to eliminate the new degrees of freedom, similar to the strategy in the previous section.
We consider a network that consists of a tree plus a single link ℓ = (m, n) with capacity κ as sketched in Fig. 2b. The edges on the paths from nodes n and m to the root node are summarised in the edge sets L and R, respectively, which we define as follows: Denote by p(m) and p(n) the set of all edges along the shortest path from the source node to the node m and n, respectively, oriented in the direction pointing away from the source. Note that these paths are unique in a tree network. Then define the following sets: such that the union of the edge set L ∪ R ∪ {(m, n)} forms a cycle. As we will see in the following, this definition turns out to be useful when studying the dissipation in the presence of a single loop. Due to the presence of the loop, we have a new degree of freedom, the cycle flow strength f. According to the KCL Eq. (2), the flows in the loopy network are given bỹ if e 2 LðnÞ F e otherwise: The value of f is fixed via the KVL: We can now evaluate the dissipation Eq. (8) in the presence of the new edge (m, n) by plugging in the relations (20) and (21), The average dissipation thus reads where we introduced the abbreviations which are functions only of the updated capacitiesk e along the sets of edges L and R. Importantly, the resulting expression no We order the loops in a colour code according to their appearance with increasing cost parameter γ: the darker the edge colour, the earlier the edge appears. For the loop that appears as the i-th loop, we denote its critical cost parameter γ ci where the loop starts to become beneficial for the dissipation-optimised network. c The transition to loop formation is discontinuous beyond the first loop: loops appearing at higher values of γ again appear with a non-zero capacity as shown in detail in the inset. Fluctuation strength is fixed to σ = 0.5 for all panels.
longer contains the updated flowsF e , but only the flows in the tree network F e , which are determined by Eq. (17). This allows us to minimise the dissipation with respect to the updated capacities k e without having to take into account the interdependence of flows and capacities. Now that we have derived an explicit equation for the dissipation in near-tree networks, we will demonstrate how to minimise the resulting expression. For tree networks, the minima of the dissipation may be calculated explicitly using the method of Lagrange multipliers (see Supplementary Note 2). In contrast to that, we have to take into account an inequality constraint κ ≥ 0 for near-tree networks as local minima may exist also at the boundaries of the domain. This can be done using the Karush-Kuhn-Tucker (KKT) conditions with the new Lagrange type functionL whereλ; μ 2 R are KKT multipliers. The minimum is then determined by the KKT conditions (see "Methods" section). This approach results in explicit equations for the optimal edge capacitiesk e , κ in near-tree networks for which we could, however, not find a closed-form solution for arbitrary networks and values of γ (Supplementary Note 2). Still, we can make use of the resulting equations to gain insight into the process of loop formation. In particular, the KKT condition for the newly added edge (m, n) with capacity κ reads B m;n ¼ ð1 þ C m;n κÞ 2 κ γÀ1 γλ _ κ ¼ 0; i.e. the capacity of the new edge either vanishes (κ = 0) or has the non-zero value given above. Importantly, we can obtain insights into the process of loop formation even without explicitly solving these equations.
How do loops emerge? We now illustrate how to make use of Eq. (26) to understand the process of loop formation. In particular, we rigorously demonstrate that loops form discontinuously as illustrated for the small tree network. Furthermore, we show that the tree remains a local minimum of the average dissipation even after the formation of a loop. We summarise these results in the following. Theorem 1 (Tree remains KKT point) Consider a linear flow network subject to the resource constraint with γ ∈ (0, 1). Then the following statements hold for the KKT points of the average dissipation 〈D loopy 〉: 1. There is always a KKT point at κ = 0, i.e. the tree is always a (potential) local minimum. 2. The KKT point at κ = 0 is isolated in the sense that we can find a real number ε > 0 such that there are no other KKT points for κ ∈ (0, ε).
The proof makes use of the fact that we can find lower and upper bounds for the functions B m,n ,λ and C m,n even without explicitly solving Eq. (26) and can be found in Supplementary Note 3. We note that the fact that the tree remains a local minimum is well-known for deterministic sources 17,24 .
We thus showed rigorously that for γ < 1, KKT points that characterise a loopy network cannot emerge through a bifurcation of the local optimum describing a tree network since the KKT point at κ = 0 is isolated. Instead, new local minima of the dissipation generally emerge elsewhere and the transition to loopy networks is discontinuous. Having understood the mechanism of loop formation, we now proceed to analyse which edges will form the first loops.
Where do loops emerge first? We now study the location of the first loop that appears in the globally optimal network as the parameters of the model are varied. We start from the regime where the global optimum is a topological tree. Consider the expression for the average loopy dissipation 〈D loopy 〉 to which a single edge (m, n) of capacity κ was added, as calculated in Eq. (23). We can find the location where loops form first by making the following approximation: assume that after the addition of the loop, the capacities of the edges e along the shortest path from the source to the loop, e ∈ R ∪ L, change only by a constant factor c (γ, e), i.e.k e ¼ cðγ; eÞ k e , whereas the other edges remain unchanged such that c(γ, e) = 1 for these edges. Looking at Fig. 3a, we can see that this is a reasonable assumption for the small network considered there. Note that the prefactor can be expected to be close to unity c(γ, e) ≈ 1 even for edges e ∈ L ∪ R if we assume that the network is very large because then the new edge will emerge with a very small capacity due to the resource Here, we defined the quantities B m,n (k e ) and C m,n (k e ) which we obtain from B m,n and C m,n (Eq. (24)) by replacing the updated capacityk e by the tree capacity k e . The last expression can then be simplified considerably by making use of Eq. (3), Thus, the emergence of loops is essentially governed by the potential drop across neighbouring vessels. Similar to how cracks in brittle materials form to relieve high elastic stresses, loop formation is determined by the relief of large pressure drops. Our explicit prediction is consistent with the idea that the reduction of pressure drops may have driven the evolution of leaf venation 30 . From a developmental perspective, it connects to work explaining plant vein formation using models where mechanical stress relief is a crucial ingredient [31][32][33] . We confirm the 'stress relief' by loop formation in terms of the potential drop by analysing the pressure drop before and after the formation of the loop in Supplementary  Fig. 1.
We now study the predictions made using Eq. (28) numerically (see "Methods" section for details). Starting from an optimal tree network, we first calculate the pressure drop (Fig. 5a). We then successively decrease the cost for new edges and monitor the order in which new loops form (Fig. 5b). Again, the transition to loop formation is discontinuous, such that loops emerge with a non-zero capacity at a critical value of the cost parameter γ c , which is highly correlated to the pressure drop (Fig. 5c). We may thus predict the location and cost parameter where loops form based on the potential drop.
Edge betweenness determines network dissipation. As we have demonstrated, all trees are-and remain-locally optimal structures and loopy networks emerge via saddle-node bifurcations. Thus, there may be a multitude of different local minima for a given set of network parameters, so a natural question that arises is the following: How can we determine which of the local minima have less dissipation than others and how can we find an order of different topologies, e.g. to find the topology that globally minimises the dissipation? Remarkably, we can trace back the answer to a purely topological property: the edge betweenness centrality.
We start by simplifying the locally optimal dissipation of the tree networks. In Eq. (17), we expressed the flows F e in a tree network using the tree matrix T . If we plug this expression into the self-consistency equation for the capacities Eq. (9), set the overall available capacity to K = 1, and plug everything into the dissipation Eq. (8) we arrive at the locally optimal dissipation in tree networks Importantly, the entries appearing in this expression only depend on the mixed moments of the sinks and their second moments. Since the sinks are i.i.d. Gaussian random variables, these moments are identical for different sinks and are given by Thus, the sum runs over identical entries and we can calculate the dissipation as Here, N p (e) is the sum over the column of the tree matrix T that corresponds to edge e. In fact, N p (e) has the following interpretation: it is the number of paths from the source s to any other node v that go through the edge e and may thus be identified as a measure of shortest-path edge betweenness 34,35 (see "Methods" section). What can we learn from this analysis for loopy networks?
To estimate the contribution of a single edge to the overall network dissipation in a loopy network, we first calculate its edge betweenness (Fig. 6a) and, based on this, the contribution it would have to the dissipation in a tree network (Fig. 6b). Adding up the resulting expressions, we arrive at the tree estimate of the dissipation in a loopy network hD Ã tree iðN p ðeÞÞ. For near-tree networks, the correlation between the estimate and the actual dissipation at local minima is almost perfect as predicted by Eq. (31) (Fig. 6c). Increasing the number of loops by tuning the noise parameter σ and the cost parameter γ, edge betweenness and dissipation remain correlated even when there is a significant number of loops present in the network (Fig. 6d). Thus, we can still understand the minimal dissipation in loopy networks based on this topological measure.
We further discuss the possibility of characterising the global tree minima of the network dissipation in Supplementary Note 6.

Discussion
In summary, we demonstrated that the transition to loop formation in optimal supply networks is discontinuous throughout different models and parameters. We explored this discontinuity in detail for a small example network and rigorously proved that the discontinuous nature of the transition persists for larger networks as well. We showed that loops emerge through a saddlenode bifurcation, explaining the discontinuous transition.
Our results shed light on recent advances in the study of optimal supply networks. While the emergence of loops through fluctuations or damage was discovered recently 15,16 , the theoretical nature of this transition was until now not well understood. Here, we closed this gap by analysing the nature of the transition to loop formation in more detail. In particular, we obtained a measure of network stress that allowed us to predict the location and parameters where loops start to form. This opens a new pathway to the understanding of loop formation in natural networks such as leaves 10 .
Our results offer a new understanding of the interplay between the structure and function of supply networks. By unveiling the relationship between the network's topological edge betweenness and its average dissipation, we established a new link between the form and function of networks. These results may aid in the understanding and design of globally optimal network structures such as biological vasculature, electrical grids, or neural networks. Our explicit prediction is consistent with the idea that the reduction of pressure drop variations may have been a factor in the evolution of leaf venation 30 . More generally, we show that globally optimal network structures may be obtained by following simple local rules for adding new links, in contrast to previous work based on pruning an existing network 26,36 .
Let us finally discuss how our results derived for linear flow models relate to other types of networks and systems. The starting point of our analysis was the fundamental trade-off between cost and resilience, which determines the optimal structure of a network, and which extends far beyond the theory of supply networks. Resilience requires additional capacity or links which can take over the load in the presence of failures or fluctuations-but these are generally costly. From a practical view of network design, the fundamental question is thus: Where and how should new connections be added that increase resilience in an optimal way? Firstly, we discuss the question where new links should be added. A large body of literature in network science approaches aspects of resilience from the viewpoint of percolation theory. The fundamental question in this purely structural treatment is: Given a network, how many nodes or links may fail before the network gets disconnected? It has been shown that a decisive quantity to assess the resilience to random failures is given by the ratio of the second and first moment of the degree distribution, 〈k 2 〉/〈k 1 〉 37,38 . These fundamental results were then used to optimise network resilience with respect to both random failures and targeted attacks 13,14 . In the case of failures, it is beneficial to add links between nodes which already have a high degree to effectively increase 〈k 2 〉. This result might appear very different from the findings of the present paper at first glance, but there are in fact common underlying principles. In supply network models, new links should be added where they will potentially attract a high flow. In percolation type models, new links should be added where they will potentially attract a high betweenness-a quantity which can also be interpreted as a flow 34,39 . As a result, one should pick nodes whose characteristic quantity, either potential θ n or degree k, differs from their surrounding.
Secondly, we consider the question of how new links should be added. The main finding of our work is that new links emerge in a discontinuous way with a finite non-zero capacity. That is, to be beneficial for the network, new links must have a certain minimum connection strength. This result has no direct equivalence in percolation approaches to network resilience since the vast majority of studies in this field considers unweighted networks only. However, there is a strong interest in network formation processes, which induce discontinuities in macroscopic connectivity of the network-including competitive percolation models 40,41 , as well as transportation network models 42 . In the context of network resilience, it has been shown that interdependencies and cascade effects can make the percolation transition discontinuous 43 .

Methods
Numerical simulation of loop formation. When analysing the transition to loop formation such as in Figs. 4 and 5, we start from an optimal tree network T for given parameters μ, σ and γ corresponding to the dissipation minimum shown in Fig. 2c. As a next step, we add all non-tree edges from the underlying triangular grid with a very small capacity that corresponds to 1% of the minimum capacity in the optimal tree, k f ¼ 0:01 Á min e2T ðk Ã e Þ, and then renormalise all capacities to make sure the resource constraint (6) holds. Finally, we then apply the iterative method described in ref. 16 to let the new topology relax to a local minimum. If this minimum contains loops despite having started very close to the (global) tree minimum, and if its dissipation is lower than the one of the tree, we conclude that the given loopy topology is favourable.
To analyse the predictive power of the pressure drop in Fig. 5, we initially consider a large optimised tree network with N = 169 nodes and cost parameter γ = 0.4 for which we calculate the pressure drop (Fig. 5a) and then increase the cost parameter from γ = 0.7 to γ = 0.99, reoptimising the network for each value of gamma. Using the procedure outlined above, we compare the predicted positions of the first loops as indicated by the initial pressure drop hðθ m À θ n Þ 2 i with the actual order in which they appear (Fig. 5b).
Evaluating edge betweenness. In Eq. (31), we derived an alternative expression for the network dissipation at local minima that is based on the edge betweenness N p (e). The edge betweenness is defined as 34,35,38 N p ðeÞ ¼ X t2V σðs; tjeÞ σðs; tÞ : Here, σ(s, t) is the number of shortest paths from node s to node t and σ(s, t|e) is the number of these shortest paths that contain the edge e. In the given setting, we consider this measure with respect to a single source s that is identified as the source node of the network. Furthermore, when analysing tree networks, there is only one path from the source to every node σ(s, t) = 1 and thus σ(s, t|e) = 1 ∨ = 0.
In the main text, the edge betweenness is calculated using a method implemented in PYTHON's NETWORKX library [44][45][46] .
Finding minima of a function with inequality constraints using KKT conditions. Consider the function D(k) of some real vector k ¼ ðk 1 ; :::; k N Þ > 2 R N that is subject to the equality constraint h(k) = 0 and the inequality constraint g(k) ≤ 0 which we assume to be described by differentiable, real-valued functions g; f : R N ! R. To identify potential maxima or minima of the function subject to the constraints, we can make use of the KKT conditions. To this end, we consider the Lagrange type functioñ whereλ; μ 2 R are called KKT multipliers. Then the following conditions, the KKT conditions, are a necessary condition for a point k * being a minimum of D(k) 47 This formulation may be used to find out whether adding a single loop to a tree network may reduce its dissipation.

Data availability
Photographs of leaf venation networks in Fig. 1 are available upon request. The topology of the Scandinavian power grid has been extracted from the open European energy system model PyPSA-Eur 49 , which is fully available online at https://doi.org/10.5281/ zenodo.3886532. Distribution grid in Fig. 1c was extracted from ref. 50 .

Code availability
Computer code will be made available at https://github.com/FNKaiser/Optimal_Supply_ Networks upon publication.