Bifurcations in adaptive vascular networks: Toward model calibration

Transport and flow networks are ubiquitous in nature and engineering, ranging from blood vessels in mammals,^1,2 glymphatic fluid networks in the brain,³ sap flow in trees, plant roots, and leaf venation⁴ to roads, pipelines, and fluid microcircuits.⁵ These networks have to efficiently distribute resources from sources to sinks under varying demands and constraints. One way to achieve this efficiency is by self-organizing adaptive networks, where the network structure dynamically adjusts to the local flow conditions. Such networks can exhibit complex topologies that include cycles and loops, which may enhance robustness and flexibility, as shown by previous studies using optimization methods.^6–9 

However, there are two main challenges in understanding and modeling self-organizing adaptive networks. The first one is to obtain reliable measurements of the network structure and dynamics in real systems, such as the vascular system of animals or plants.^1,10,11 The second one is to develop and validate mathematical models that can reproduce the observed network features and predict their behavior under different scenarios.

In this article, we address these challenges for a specific class of self-organizing adaptive networks: vascular network systems. We focus on the calibration of a model that describes how the network adapts to changes in blood flow and pressure based on local feedback mechanisms.¹² Such mechanisms are known to occur—in general—on a range of time scales in the mammalian vasculature, ranging from seconds (myogenic responses induced by smooth muscle cells encompassing arterioles) to months (remodeling of blood vessels).^12,13 In particular, we build on a mathematical model of vessel adaptation originally proposed by Hu and Cai.¹⁴ This model produces results from our previous studies as a limiting case.^15,16 The network structure emerging from this adaptive model depends on the spatial distribution of sources and sinks and on the amplitude and frequency of flow fluctuations. The model can capture certain aspects of vascular network formation and adaptation, such as the emergence of cycles, the scaling of vessel diameters, and the balance of pressure fluctuations.^14–16 An important model parameter is the exponent $γ$ in Eq. (2), which tunes between a concave vs convex vessel growth in the adaptation rule. Hu and Cai argued¹⁴ that this exponent should be $γ=1/2$ for vascular systems. To test this hypothesis, we first carried out bifurcation analysis for a simple triangular network motif, composed of one constant source and two sinks subject to a fluctuating demand, in dependence on fluctuation amplitudes and the cost exponent $γ$⁠. We then conducted numerical simulations on the model on network topologies extracted from perivascular networks of rodent brains¹⁷ to study the emergent network structure as vessels appeared or vanished while varying the cost exponent $γ$⁠. Therefore, we aim at addressing two interlinked questions:

• Is a specific measured flow network consistent in terms of physical reality?

• Is the adaptive dynamic model consistent with measured network data?

This article is structured as follows. Section II explains the mathematical model of adaptation and briefly discusses the nature of the vascular network data used for comparing model predictions. Section III presents detailed bifurcation analysis for the triangular motif and discusses the numerical simulation on the perivascular rodent brain networks, and Sec. IV concludes our study with a discussion.

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$⁠. Furthermore, we assign length $L i j$ to each vessel segment between node $i$ and node $j$⁠.

As an empirical testbed for the model, we employ a vascular network from the neocortex in a rodent.¹⁷ The network with identifier 012208 has 855 connections (vessel segments) connecting 826 nodes, 416 of which are sink nodes. A node $i$ is a sink node, $i∈S$⁠, if and only if $i$ has degree 1 or degree 2 meaning that fewer than three vessel segments meet at $i$⁠. In addition to the list of node pairs forming connections, the data contain the coordinates $( x i, y i)$ in two-dimensional Euclidian space for each node $i$⁠. The length of a vessel segment between nodes $i$ and $j$ is given by $L 23= ( x i − x j ) 2 + ( y i − y j ) 2$ being the Euclidian distance between the nodes. Here, we assume that vessel segments have negligible curvature, thus extending along a straight line.

We consider a triangular motif with one source $h 1=1$ and two fluctuating sinks $h 2$ and $h 3$⁠. We assume that the drive $h 2 , 3(t)$ is rapid with a well defined 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$⁠.

We denote solution branches by $B=( C 12, C 23)$ where we generally assume $C 12, C 23≥0$⁠. The following stationary solutions are possible. We always require $C 12>0$⁠; otherwise, no net transfer of mass is possible along the network, contradicting mass balance. Solutions can either be “tree-like” if $C 23=0$⁠, i.e., $B ∧=( C 12,0)$⁠, or “triangular” if $C 23>0$⁠, i.e., $B Δ=( C 12, C 23)$⁠. We consider only $γ>0$⁠, and Eqs. (6a) and (6b) inform us that tree-like solutions exist only if $0<γ≤2$ (note that $x 0→1$ as $x→0$⁠); if, on the other hand, $γ>2$⁠, the right-hand sides of (6b) are indeterminate $C 23=0$ (or infinite if $C 23→ 0 +$⁠). Thus, we consider tree-like solutions for $0<γ≤2$⁠.

Applying the tree-like solution to (8a), we have $p 2− p 1= h 2/ C 12$⁠. Imposing further a stationary nontrivial solution for (6a) yields $1= C 12 1 − γ ⟨ ( p 1 − p 2 ) 2 ⟩= C 12 − 1 − γ ⟨ h 2 2 ⟩$ so that $⟨ h 2 2 ⟩= x 1 + γ= C 12 1 + γ$ or $x= C 12= ⟨ h 2 2 ⟩ 1 1 + γ$⁠. The tree-like branch is, therefore, given by $B ∧=( C 12, C 23)= ( ⟨ h 2 2 ⟩ 1 1 + γ , 0 )$⁠.

We evaluate the Jacobian (11a)–(11d) for the tree-like branch. Using $y= C 23=0$⁠, we have $ξ=0$⁠, $ξ x=0$⁠, and $ξ y=1/x$⁠. For $f x$⁠, all terms cancel that are proportional to both powers of $x$ and $ξ$⁠, $ξ 2$⁠; thus, $f x B ∧=−1−γ$⁠.

Similarly, for $f y$⁠, we only retain the term proportional to $ξ y=1/x$⁠; we use that $ξ y⋅ x − γ= ( ⟨ h 2 2 ⟩ ) − 1$ and obtain $f y B ∧=−2a−a⟨ h 2 h 3⟩/ ⟨ h 2 2 ⟩$⁠.

More care is needed to evaluate $g x$ and $g y$⁠, where products and powers of $y$ appear together with terms proportional to $y$⁠, such as $ξ, ξ x, ξ y∝y$⁠, thus resulting in terms of the form $y p$ with some exponent $p$⁠. If $p>0$⁠, associated factors become zero as $y→0$⁠; but if $p<0$⁠, the associated term is indeterminate. Consider first $g x∝ y 2 − γ$⁠. Since $y=0$ on the tree-like branch, we have $g x B ∧=0$ for $γ≤2$⁠. It follows that the associated Jacobian is an upper triangular matrix.

Next, consider $g y$⁠, for which the first term is $y − 1 − γ ξ 2∝ y 1 − γ$ and the second is $y − γξ ξ y∝ y 1 − γ+ y 2 − γ$⁠. Provided that $γ<1$⁠, $g y$ is determinate for $y=0$⁠, in particular, all terms go to 0 as $y→0$⁠, and we obtain $g y B ∧=−1$⁠. For $γ=1$⁠, we use that the term with $y 1 − γ= y 0→1$ as $y→ 0 +$ (all other terms containing $y$ are 0), and we get $g y B ∧=−1+ ⟨ h 2 − h 3 ) 2 ⟩/ ⟨ h 2 2 ⟩ 2 / ( 1 + γ )$⁠. Finally, considering $1≤γ<2$ the expression $y 1 − γ→+∞$ (as long as $γ<2$⁠) as $y→ 0 +$ so that we have $g y→+∞$⁠.

Note that letting $γ=1$ and using (A2a), we have $λ 2 B ∧= 3 α 2 − 1 α 2 + 1$⁠. This determines critical strength of fluctuations for $α c=1/ 3$⁠, where the tree-like branch loses stability in a transcritical bifurcation. Indeed, this is the same result we obtained in a previous study.¹⁵ 

To sum up, the tree-like branch is always stable for $0<γ<1$⁠; for $γ=1$ only if $α< α c$⁠; and for $1<γ≤2$⁠, it is “hyper-repellent.”

We complete the consideration of the system with three nodes by observing response to varying $γ$ is varied while keeping $α$ constant. This scenario, for different lengths of the cross-connection $L 23$⁠, is the most relevant one to the parameter scan with the empirical network in Sec. III B. Figure 3 shows that a pair of a stable and an unstable triangular solution disappear in a saddle-node bifurcation as $γ$ decreases through a critical value. Solutions are obtained numerically from (28) and (29) by tuning $C 12$ such that $α=1$ results for $γ$ and $L 23$⁠.

We analyze the stationary solutions of the model on the empirical vascular network described in Sec. II C. While keeping fluctuations at maximum (⁠ $α=1$⁠, single moving sink), we vary the cost exponent $γ$ and observe the conductances after settling to stable stationary solutions of (4). We observe that for $γ=1$⁠, all vessel segments have positive conductivity, $C i j>0$⁠. As $γ$ decreases, there is a sequence of bifurcations in each of which a single connection turns from positive to zero conductance. Figure 4(a) shows the $γ$-dependence of five conductances involved in bifurcations earliest, i.e., at relatively large $γ$⁠.

As $γ$ decreases toward zero, a total of 30 bifurcations occur. In each of these, a single connection loses its non-zero conductance. The 30 connections involved in bifurcations are drawn as thick lines in the network illustration in Fig. 4(b). The 825 connections retaining non-zero conductance for $γ→0$ form a tree on the 826 nodes: by following these remaining connections, all nodes are mutually reachable by a unique path so that there are not any cycles.

Empirically detected vessel segments naturally have a non-zero conductance. Assigning positive conductances to all existing vessel segments is required for a realistic parameter choice of the model. The result in Fig. 4 suggests that a choice of $γ=0.5$¹⁴ does not reproduce the empirical network with all positive conductances. However, choices of $γ$ closer to a value of $1$ do provide stable stationary solutions with all conductances positive.

Here, we have assumed that the case of maximal fluctuation amplitude, $α=1$⁠, maximizes conductances in the stable stationary solutions. While this assumption is plausible in the light of the results obtained for the model so far, a demonstration of the assertion is left for future work.

Since empirical data—in general—cannot be guaranteed to be free of error, one may ask if the conductances vanishing at large values of $γ$ belong to erroneously registered vessel segments. This appears unlikely when considering the length of the connections. The lengths of the five vessel segments featured in Fig. 4(a) lie in the range $[28.6,79.6]$ to be compared to median $15.3$ and mean $20.7$ for the distribution of all 855 segment lengths in the network.

The model for vascular network adaptation by Hu and Cai¹⁴ exhibits rich phenomenology in terms of bifurcations even in the simplest non-trivial system with one source and two fluctuating sinks. In the biologically most relevant regime of cost exponent $0<γ<1$⁠, tree-like solutions are stable irrespective of other parameters. Cyclic (triangular) behavior appears as a pair of a stable and an unstable solution in a saddle-node bifurcation. In the supercritical regime, there is multi-stability. Tree-like or cyclic behavior of the system depends on initial conditions and basins of attraction in the space of conductances.

Testing the model on an empirical network with hundreds of nodes reveals cyclic solutions appearing by saddle-node bifurcations, just as in the three-node system. In contrast to the rigorous analysis of small systems, we were able to find stable stationary solutions only by integrating the equations of motion from a random sample of initial conditions. Asymptotics of the trajectories being consistent across initial conditions makes it plausible that we are identifying the stable fixed point whose basin of attraction is dominating the phase space.

In terms of model calibration, our test of the model on an empirical network suggests that the cost exponent is closer to a value of $γ=1$ than the previously suggested¹⁴ value $γ=1/2$⁠. At the latter value, the model’s prediction is incompatible with the simultaneous presence of all empirically observed vessel segments. In the stable stationary solutions found at $γ=1/2$⁠, a value of zero conductance is assigned to several vessel segments, incompatible with the existence of these connections. Increasing $γ$⁠, however, the model assigns realistic non-zero conductances to all connections in agreement with them being observed as vessel segments in the empirical network.

Future work is to improve methods and aim at exhaustive identification of all stationary solutions also in large networks. Structural features of vascular networks are to be exploited by a suitable method. In particular, network sparseness may facilitate efficient search of solutions. Typical vascular networks have a small cyclomatic number, meaning all cycles can be removed by breaking a relatively small number of connections. The structural similarity to trees is also seen in low tree-width¹⁸ of vascular networks that we have identified in preliminary work.^19,20 Ideally, the tree-like structure can be exploited to recursively compute the set of stationary solutions on a network.

The robustness of results is to be tested by varying the type of sink fluctuations. While the implementation of fluctuations in the present paper is in line with our own previous work,^15,16 it differs from the original definition by Hu and Cai.¹⁴ In the latter, each sink $i$ independently opens with probability $p$⁠; otherwise, it fully closes. The case of maximal fluctuations (⁠ $p→0$⁠), however, is that of a single moving (open) sink and, thus, identical to a maximum fluctuation amplitude (⁠ $α=1$⁠) in the present work. Going beyond these two types of fluctuations, spatiotemporally correlated sinks and their effects on a tree-like vs cyclic network structure are worth exploring.

Several situations are conceivable that prompt the need to use an model that adapts vessel diameters (2), toward finding stable configurations (vessel diameters) that are consistent with flow parameters. First, some experimental network data only provide information about topology (connectivity), but not about vessel diameters, including the data used in this study¹⁷ and elsewhere.²¹ However, in some circumstances, such data may be relevant to add to simulations—in such cases, the adaptive model (2) can be used to generate data with vessel diameters that are consistent with prevalent flow parameters.²² Large scale whole brain vascular networks measured in rodents have recently become available.¹¹ While these data contain information about vessel widths, the data are extracted from animals post mortem. This raises two questions: first, does the data contain erroneous vessels (presumably small or short)—this question relates to the consistency of the network data. Second, what are the actual in vivo diameters that appear when blood pressure corresponds to the alive state of the investigated organism. These questions are also relevant when constructing realistic models of vasculature in whole organs,^2,23 and adaptive models, such as the one studied here, may be used to check for consistent edge weights (conductivities) along the network.

In future work, one may raise the question of relevance of certain vessels’ large networks where the size of the network is prohibitive for the computation/simulation of the entire Kirchhoff network; i.e., is it possible to simplify (i.e., coarse grain) the network toward computational accessibility? Moreover, myogenic responses are known to lead to oscillatory behavior in vessel diameters. Are such oscillations a feature of the self-organization intrinsic to the adaptive dynamic nature of the network? Or are these oscillations due to external forcings due to neuronal activity (astrocytes and pericytes)? Finally, it would also be interesting to conduct experiments with flow networks. An example of a potential candidate system is an experimental system with adaptive memristor flow circuits proposed very recently.⁵ 

K.K. acknowledges financial support from MCIN/AEI/10.13039/501100011033/FEDER, UE (Project No. PID2021-122256NB-C22).

The authors have no conflicts to disclose.

Konstantin Klemm: Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). Erik A. Martens: Formal analysis (equal); Investigation (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.

It is useful to determine expressions for the forcing terms appearing in the triangular network system (6b), $⟨ h 2 2 ⟩$⁠, $⟨ ( h 2 − h 3 ) 2 ⟩$⁠, and $⟨ h 2 h 3⟩$⁠, as functions of the fluctuation amplitude, $a$⁠. To this end, we consider two types of forcing.