Non-local transitions and ground state switching in the self-organization of vascular networks Open Access

Transport or distribution networks are crucial for the proper functioning of a wide range of natural and technological systems. An important question is what network structures are optimal in terms of energy consumption. A classical result is that optimal transport networks are minimal spanning trees, provided that loads are stationary;¹ however, many realistic settings display fluctuating loads, in which case cyclic motifs allow for short cuts (shunts), thus improving energy consumption. A number of studies demonstrated the necessity of such cycles using optimization techniques;^2–6 an alternative dynamical systems approach was put forward by Hu and Cai⁷ in a model where vessel diameters adapt dynamically in response to a power dissipation function, leading to adaptive network dynamics.⁸ The bifurcation structure for this model has been the subject of a number of studies,^7,9–11 which assumed that fluctuations occur very rapidly, i.e., time scales associated with fluctuations and the vessel conductance were perfectly separated. Here, we relax this assumption and use the energy methods to elucidate the role of slower fluctuations. In a parameter regime of bistability with both tree-like and cyclic motifs, noise renders cyclic motifs more favorable even though tree-like configurations have lower energy.

In natural and engineered systems, transport and flow networks are essential to the distribution of resources.^12,13 Transport networks range from mammalian vascular vessel networks,^14,15 leaf venation,¹⁶ and brain glymphatic fluid pathways¹⁷ to arboreal sap conduits;¹⁸ and man-made infrastructures such as roads¹⁹ and microfluidic circuits.²⁰ Transport networks are tasked with the critical function of efficiently channeling resources from sources to sinks amidst diverse demands and constraints. The efficiency of these networks is often attributed to their ability to self-organize and adapt dynamically to local flow conditions,^8,21 resulting in complex network topologies that incorporate cycles and loops^2,3,22,23 and community structures.²⁴ Such configurations are not only indicative of robustness and flexibility but also optimized in terms of energy consumption.

While the topology of technological transport networks tends to be more static and, thus, offers limited adaptability, biological networks such as the mammalian vascular system demonstrate remarkable flexibility and operate reliably over a vast parameter range,^25–27 thereby avoiding the risk of failure even in extreme scenarios. In the mammalian vasculature, vessels dynamically adjust diameters in response to changes in flow properties like pressure and shear stress.^28,29

A range of studies have been dedicated to the optimal network configuration in terms of energy consumption. While an optimal distribution network assumes the structure of a minimal spanning tree when loads are static,¹ a range of studies revealed that fluctuating loads instead result in the emergence of cyclic motifs in order to optimize energy consumption.^2–6 Intuitively, this can be explained by the realization that shunting (shortcutting) nearby nodes is cheaper than rerouting the fluctuating demands up and down a number of levels in a tree structure. An alternative line of research, introduced by Hu and Cai⁷ takes the perspective of dynamical systems theory, instead considering the time evolution of the network by invoking a dynamic rule for vessel diameters (conductances), such that they adapt in response to a power consumption function controlled by an exponent $γ$⁠. Bifurcation analysis for a minimal triangular motif and larger networks in the limit of rapid fluctuations revealed that cyclic motifs may emerge in either transcritical or saddle-node bifurcations, depending on the specific value of $γ$⁠; moreover, the saddle-node bifurcations give rise to bi-stability between tree-like and cyclic configurations.^9–11 In such a dynamical formulation, the time scales of the adaptation rule for the conductances and of fluctuations are of order $∼1$ and $∼T$⁠, respectively—thus, a time scale separation is made explicit. Crucially, biological settings do not adhere to a strict time scale separation, i.e., we cannot claim $T≪1$⁠. This raises the question how the time scale associated with fluctuations, $T$⁠, affects bifurcations and bi-stability. Specifically, previous studies used the model assumption of rapid fluctuations, $T≪1$⁠, but it is unclear to what extent this assumption accurately predicts the favored network configurations.

This article seeks to understand the response of adaptive flow networks, subject to non-rapid fluctuation. To achieve this, we first derive a framework for the rapid fluctuation limit, using an “energy function,”^7,30 in physical terms describing the energy consumption due to dissipation along the total network. A minimal energy principle³¹ is used to characterize the stability of equilibria and offers a more global view on stability. Using this methodology, we are then able to study how the timescale and level of fluctuations affect the system dynamics under various conditions. Our study reveals that fluctuations may stabilize cyclic configurations even when tree-like configurations are energetically favorable. Thus, we clarify the influence of the time scale associated with fluctuations, an aspect that was neglected in previous studies.

Let $V$ denote a set of nodes of a network with $N= |V |<∞$ and $A⊆N×N$ the set of edges. The edges are bidirectional; therefore, $(i,j)∈A$ implies $(j,i)∈A$⁠. Each node is assigned pressure $p i$⁠. The edge flow is $Q i j>0$ from node $i$ to $j$⁠. We assume that the network is resistive and linear, i.e., Ohmian with $Q i j= C i j( p i− p j)$⁠, where an edge carries the property of conductance between nodes $i$ and $j$ with $C i j= C j i>0$ only if $(i,j)∈A$⁠.

For stochastic sink fluctuations occurring on a time scale $T>0$⁠, the net flows in all sinks are constant on all time intervals $I n=[(n−1)T,nT[$ with $n∈ N$⁠. Independently for each $n∈ N$⁠, we draw the high-flow sink $s n∈S$ uniformly at random (with probabilities $1/ |S |$⁠).

In this section, we discover non-local transitions where the ranking of solutions with respect to the energy functional $E$ changes, and, thus, switching between ground (minimum $E$⁠) states occurs (see Sec. III E and III F). The impact of these transitions in a system with finite timescale fluctuations is studied in Sec III G. We provide the context for these findings by first reviewing the known saddle-node bifurcations giving rise to cyclic solutions^11,23 (Secs. III A–III C). For particular values of cost exponent $γ$⁠, we reveal additional bifurcations (Sec. III D).

We consider a triangular network motif with one source $h 1=1$ and two fluctuating sinks $h 2$ and $h 3$⁠. This motif has already been considered in the previous studies^9–11 as it serves as an elementary building block for larger networks in setups with one source and many sinks (terminal loads).

We assume that drive $h 2 , 3f(t)$ fluctuates on a rapid time scale, i.e., $T=T( h 2 , 3)≪1$⁠. The sources then obey $⟨ h 2⟩−⟨ h 3⟩→0$⁠, and no net pumping occurs between nodes $k=2$ and $k=3$⁠.

Accordingly, we seek solutions, $⟨ C k l⟩$⁠, averaged over rapid fluctuations with characteristic time scale $T$⁠. We consider that $C i j$ changes on a slow time scale, $C i j→⟨ C i j⟩$ as $T→0$⁠. From now on, we, therefore, use $⟨ C k l⟩$ and $C k l$ interchangeably and simplify notation by omitting $⟨⋅⟩$ around the conductances. Note that the symmetry of both the triangular network motif and of the (rapid) fluctuations implies that $⟨ C 12⟩=⟨ C 13⟩$⁠.

Relations between the average squared terminal fluctuations, $h 1, h 2$⁠, and $h 3$⁠, and the fluctuation amplitude $α$ relevant to evaluating the pressure difference terms in Eq. (10), are derived in Appendix A.1. in Ref. 11. They are $⟨ h 2 2⟩= 1 + α 2 4,⟨ ( h 2 − h 3 ) 2⟩= α 2,⟨ h 2 h 3⟩= 1 − α 2 4$⁠.

We summarize the results of the bifurcation behavior and stability analysis for stationary solutions occurring in the triangular network motif.^9,11 Bifurcation diagrams and phase diagrams for the selected values of $γ$ are shown in Fig. 1 and in Figs. 3(a)–3(c), respectively.

Finally, for $γ>1$⁠, the saddle-node bifurcation is annihilated in the transcritical bifurcation point TC and the tree-like and cyclic solutions stably co-exist for any value of $α$ [see Fig. 2(c)]. Note that the tree-like and cyclic solutions display nearly identical values for $C 23$ below a specific value of $α$ and only begin to substantially differ above a certain $α$-value.

The stability diagram from these bifurcations is shown in Fig. 1. We remind the reader that not all values of $α$ and $γ$ are of interest to us, due to certain limitations in their physical interpretations. First, we limit our attention to $α≤1$ to avoid sign reversal in the terminal fluctuations (see Sec. II B). Second, for the exponent value $γ>1$⁠, the maintenance power of a vessel segment scales superlinearly with its conductance. This case appears less realistic than the linear (⁠ $γ=1)$ and sublinear (⁠ $γ<1$⁠) ones. Under $γ>1$⁠, maintenance power may always be lowered by replacing a vessel segment with several parallel segments while keeping total conductance.

We now turn our attention to the model in terms of minimizing energy consumption $E$ as given by (2). Relative values for $E$ evaluated at different equilibria allow for a more global characterization of stability, which is particularly useful to discuss the system behavior in the presence of multistability (occurring for $α> α \,SN$ and $γ≤1$⁠) and of noise-driven dynamics in the regime of non-rapid slow fluctuations, which we will discuss later in Sec. II B.

We evaluate $E$ along the $C 12$-nullclines [blue curves in Figs. 3(a)–3(c)], as shown in Figs. 3(d)–3(f). In particular, maximum and minima of $E$ give the energy consumption at the saddle and stable fixed points. Fixing $α=1$ and choosing subcritical $γ=0.76< γ 0$⁠, the unique minimum of $E$ is at the tree-like solution [see Fig. 3(d)]. For values $γ=0.80> γ 0$ and $γ=0.85> γ 0$⁠, i.e., above the saddle-node bifurcation, the stable cyclic solution is a local minimum; the unstable cyclic solution (saddle) is a local maximum of $E$ as seen in Figs. 3(e) and 3(f). A change of behavior occurs between $γ=0.80$ and $γ=0.85$ at $γ= γ 1$⁠. For $γ< γ 1$⁠, the tree-like solution with $C 23=0$ is the unique global minimum of $E$⁠. For $γ> γ 1$⁠, the stable cyclic solution uniquely globally minimizes $E$⁠. Thus in addition to the saddle-node bifurcation at $γ= γ 0$⁠, the system undergoes another transition at $γ 1$⁠, where the ground state changes from the tree-like fixed point to the stable cyclic fixed point.

The phenomenon of switching ground state generalizes to parameter values $α≤1$⁠. In the stability diagram in Fig. 1, the curve $α E(γ)$ further subdivides the parameter regime of bistability. Above the curve, i.e., for $α> α E(γ)$⁠, the stable cyclic solution is the ground state configuration. Below the curve, the tree-like solution globally minimizes $E$⁠.

Let us consider a larger network with more cycles and, thus, more possibilities for the presence or absence of cross-edges. Figure 4 shows the test network with 9 nodes and 12 edges, where we denote the four vertically drawn edges as cross-edges, the remaining ones as tree-edges. There are $2 4=16$ combinations of conducting/non-conducting cross-edges; 15 of these are cyclic (with at least one cross-edge conducting); the remaining combination is the tree without any conducting cross-edges.

Out of the 15 cyclic combinations, 7 are not observed as stable equilibria for any values of $γ∈[0,1]$⁠, while $α=1.0$ remains constant. For each of the remaining eight cyclic combinations, Fig. 4 displays the energy values for the stable equilibrium on the $γ$-interval where the cyclic equilibrium exists. As the parameter $γ$ is varied, the ranking of equilibria by energy varies. In particular, which equilibrium constitutes the ground state is strongly parameter-dependent. At $γ=1.0$⁠, the ground state is the one with all cross-edges except ${2,3}$⁠. At $γ≈0.86$⁠, $γ≈0.67$⁠, the ground state configuration switches between different cyclic equilibria. Below $γ≈0.30$⁠, the ground state is the tree.

We now relax the time scale separation between the sink fluctuations and the network adaptation. Accordingly, we no longer require that $T≪1$⁠, and we analyze the tree-node motif under stochastic fluctuations with time scale $T$ as described in Sec. II B. We choose the fluctuation amplitude $α=1.0$⁠. Thus at any point in time, we either have $h 2(t)=−1$ and $h 3(t)=0$⁠, or $h 2(t)=0$ and $h 3(t)=−1$ as the net flow at sink nodes 2 and 3. We integrate the dynamics for conductances $( C 12, C 13, C 23)$ according to Eq. (6). We record the distribution of $C 23$ values and plot histograms in Fig. 5.

For fluctuation time scale $T=1.0$ (upper row of panels in Fig. 5), the same distributions are obtained regardless of the chosen initial conditions. At $γ=0.76$⁠, in the subcritical regime with respect to the saddle-node bifurcation, the configurations encountered are tree-like with $C 23$ close to zero. For the supercritical choices $γ=0.80$ and $γ=0.85$⁠, mostly cyclic configurations are obtained [see Figs. 5(b) and 5(c)]. In these panels, for comparison we demark the value of $C 23$ corresponding to the stable cyclic equilibrium obtained for the rapid fluctuation limit (⁠ $T≪1$⁠) with a vertical dashed line. The distribution of $C 23$ peaks remains (relatively) close to this value for both $γ=0.80$ and $γ=0.85$⁠. For $T=1.0$ and all values of $γ$⁠, the distribution of $C 23$ is independent of the initial condition.

At shorter fluctuation time scale $T=0.5$ and supercritical $γ∈{0.80,0.85}$⁠, the observed $C 23$ distribution does depend on the initial condition [see Figs. 5(e) and 5(f)]. Using the initial conditions close to the tree-like solution, we observe only tree-like configurations with values of $C 23$ concentrating toward values close to zero. As an initial condition sufficiently far from the tree-like configuration, we use $( C 12, C 13, C 23)=(1,1,1)$ and observe mostly tree-like configurations.

Figure 6 shows an example of the conductances evolving under slow sink fluctuations. Notice that when sink $s$ is active during a given time interval, the direct conductance $C 1 s$ from source node 1 to active sink $s$ is strengthening, while the other two conductances are exponentially decreasing. In Fig. 6, this is seen for time interval $[62.0,69.0]$ where sink $s=3$ is constantly active. After switching to sink $s=2$⁠, both $C 12$ and $C 23$ increase rapidly. The increase of $C 23$ is due to the large value of $C 13$ causing a low pressure difference between source 1 and sink 3 at this time. This makes the direct connection via $C 12$ and the indirect connection via sink 3 involving $C 23$ almost equally cost-efficient in establishing flow from source 1 to current sink 2. The mechanism explains the broad fluctuations in $C 23$ at $T≫1$⁠. The cross-conductance intermittently assumes values close to zero without being drawn into the equilibrium at $C 23=0$⁠.

Note that considerations using the energy function are only (strictly) valid for the deterministic limit (or $T≪1$⁠); in the case of finite time fluctuations, the energy function is no longer a Lyapunov function for the system.

The model by Hu and Cai describes a vascular network with adaptive conductances at two levels, namely, (i) an energy landscape established on the space of conductance values, with the energy function $E$ being the maintenance cost plus the power needed to sustain the pipe flow, summed over all vessel segments; and (ii) a system of coupled differential equations describing the time evolution of the conductances. Levels (i) and (ii) are linked through the fact that the conductances’ evolution is chosen such as to never increase the value of the energy function $E$⁠. Thus, each critical point of $E$ is a fixed point of the dynamics. In particular, each local minimum of $E$ is a stable fixed point (Fig. 3).

The model shows multistability already for the system of three nodes with one source and two sink nodes. The energy function $E$ is a means of discerning the stable solutions. This establishes another critical curve in parameter space where a switch in the ground state occurs: the two stable solutions, tree-like and cyclic, exchange their ranking in terms of energy and, thus, their thermodynamic stability. This phenomenon is illustrated in Fig. 4 for an example network with one source and eight sinks.

The presence of multistability brings up the question, which solution is approached by the dynamics. Specifically, is the energetically lower solution always preferred? We made this question precise by studying the effect of stochastic sink fluctuations, occurring on a time scale of order $T$ compared to the time scale of the conductance dynamics which is of order $1$⁠. These fluctuations were averaged over when they were assumed to occur on a faster time scale compared to the dynamics of the conductances (i.e., $T≪1$⁠), resulting in the deterministic dynamics given by Eqs. (10) and (11) (see the results in Figs. 1–3 and 7). However, when the associated time scale separation is relaxed, stochasticity is introduced in the dynamics (Fig. 5). The main result of the analysis is that the system exhibits mostly cyclic configurations even when the tree-like solution is energetically favorable in the deterministic limit for sufficiently large fluctuation strength, $α$⁠.

A question related to multistability is, thus, how large the basin of attraction is for the tree-like and cyclic solutions, respectively. The likelihood that the vascular network remains in the most energetically preferable state depends on the state’s stability against significant perturbations. While linear-stability analysis gives local information of stability, basin stability provides a more global view on this issue.³⁴ It is not straight forward to define a trapping region (such as separatrices), especially while incorporating varying values for the system parameters. Therefore, we did not attempt to devise a measure for the basin volume; however, the inspection of the phase portraits in Figs. 3(b) and 3(c) displays that the phase space is subdivided into regions of attraction by the separatrices emanating from the saddle point $S$⁠, i.e., the stable and unstable manifolds. Thus, we may estimate that the basin size³⁵ for the tree-like solution scales with the distance between the stable fixed point $Λ$ and $S$ and that the basin size for the cyclic scales like the distance between $S$ and the stable node $Δ$⁠. Note that all three fixed points more or less lie on a straight line. Thus, qualitatively we can see that the basin size for cycles $Δ$ is comparatively much larger than the one for $Λ$⁠, which at least for the shown parameters explains the higher likelihood to reside near the cyclic attractor when terminal fluctuations are sufficiently slow.

Shifting focus to larger systems, we have also analyzed the energy landscapes and dynamics on a large system based on the vascular networks empirically obtained by Blinder et al.³³ Previous analysis by the authors,¹¹ considering the rapid fluctuation limit $T≪1$⁠, had shown that a saddle-node bifurcation creates a pair of fixed point solutions with non-zero conductance for each cross-edge in a large network, analogous to the case of the triangular motif. Figure 7(a) shows an example of such a bifurcation in an empirical network³³ for $α=1$ with $T≪1$⁠. Similar to the energy function of the small system evaluated along the nullclines in Fig. 3, of Figs. 7(b)–7(d) show slices through the energy landscape for the large network. The main observation from the three-node motif is reproduced in the large network: Close to the bifurcation, the tree-like solution lies energetically lower than the cyclic solution; however, this ranking switches when tuning the parameter $γ$ farther away from the transition. We conjecture that the saddle-node bifurcations and the additional change of behavior in terms of the energy minimality of solutions are generic features of the model on all networks involving cycles.

For future research, it could be interesting to also consider other types of time-dependent networks. In particular, how are energy configurations altered in co-evolving networks where edges are dynamically created/deleted,^8,21 such as in growing supply networks?³⁶ 

K.K. acknowledges support from Project No. PID2021-122256NB-C22 funded by MCIN/AEI/10.13039/501100011033/FEDER, UE. We acknowledge the financial support from the Royal Swedish Physiographic Society of Lund.

The authors have no conflicts to disclose.

Konstantin Klemm: Conceptualization (equal); Formal analysis (equal); Writing – original draft (equal); Writing – review & editing (equal). Erik A. Martens: Conceptualization (equal); Formal analysis (equal); Writing – original draft (equal); Writing – review & editing (equal).

Data sharing is not applicable to this article as no new data were created or analyzed in this study.