Transport networks are crucial for the functioning of natural and technological systems. We study a mathematical model of vascular network adaptation, where the network structure dynamically adjusts to changes in blood flow and pressure. The model is based on local feedback mechanisms that occur on different time scales in the mammalian vasculature. The cost exponent $\gamma $ tunes the vessel growth in the adaptation rule, and we test the hypothesis that the cost exponent is $\gamma =1/2$ for vascular systems [D. Hu and D. Cai, Phys. Rev. Lett. **111**, 138701 (2013)]. We first perform bifurcation analysis for a simple triangular network motif with a fluctuating demand and then conduct numerical simulations on network topologies extracted from perivascular networks of rodent brains. We compare the model predictions with experimental data and find that $\gamma $ is closer to 1 than to 1/2 for the model to be consistent with the data. Our study, thus, aims at addressing two questions: (i) Is a specific measured flow network consistent in terms of physical reality? (ii) Is the adaptive dynamic model consistent with measured network data? We conclude that the model can capture some aspects of vascular network formation and adaptation, but also suggest some limitations and directions for future research. Our findings contribute to a general understanding of the dynamics in adaptive transport networks, which is essential for studying mammalian vasculature and developing self-organizing piping systems.

Most modeling studies of vascular transport compare model outcomes and real networks in terms of summary statistics, such as the existence of cycles, their length distribution, etc. Here, we take a step toward alignment of network data and models. We explore the parameter range in which the model sustains a given empirical network as a stable solution. The present work is, thus, taking a step toward model calibration for vascular networks. The calibration is facilitated by detailed and mostly rigorous analysis we perform for the stationary solutions of the Hu–Cai model on a small test network.

## I. INTRODUCTION

Transport and flow networks are ubiquitous in nature and engineering, ranging from blood vessels in mammals,^{1,2} glymphatic fluid networks in the brain,^{3} sap flow in trees, plant roots, and leaf venation^{4} to roads, pipelines, and fluid microcircuits.^{5} 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.^{12} 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.^{14} 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 $\gamma $ in Eq. (2), which tunes between a concave vs convex vessel growth in the adaptation rule. Hu and Cai argued^{14} that this exponent should be $\gamma =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 $\gamma $. We then conducted numerical simulations on the model on network topologies extracted from perivascular networks of rodent brains^{17} to study the emergent network structure as vessels appeared or vanished while varying the cost exponent $\gamma $. 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?

^{17}However, our comparison of model predictions with experimental data suggests that—in contrast to Hu and Cai’s study

^{14}— $\gamma $ is closer to 1 than to $1/2$ when chosen such that the model is able to sustain the topological structure found in real world networks.

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.

## II. MODEL, NETWORK DATA, AND METHODS

### A. The model by Hu and Cai

Let $V$ denote a set of nodes of a network with $N= |V |<\u221e$ and $A\u2286N\xd7N$ the set of edges. The edges are bidirectional; therefore, $(i,j)\u2208A$ implies $(j,i)\u2208A$. 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\u2212 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)\u2208A$. Furthermore, we assign length $ L i j$ to each vessel segment between node $i$ and node $j$.

^{14}proposed an adaptation model for biological transport networks,

### B. Sink fluctuations

### C. Pial network of neocortex in rodents

As an empirical testbed for the model, we employ a vascular network from the neocortex in a rodent.^{17} 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\u2208S$, 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 \u2212 x j ) 2 + ( y i \u2212 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.

## III. RESULTS

### A. A system with one source and two sinks

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)\u226a1$. The sources then obey $\u27e8 h 2\u27e9\u2212\u27e8 h 3\u27e9\u21920$, and no net pumping occurs between nodes $k=2$ and $k=3$.

#### 1. Stationary solutions

We denote solution branches by $ B=( C 12, C 23)$ where we generally assume $ C 12, C 23\u22650$. 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 \u2227=( C 12,0)$, or “triangular” if $ C 23>0$, i.e., $ B \Delta =( C 12, C 23)$. We consider only $\gamma >0$, and Eqs. (6a) and (6b) inform us that tree-like solutions exist only if $0<\gamma \u22642$ (note that $ x 0\u21921$ as $x\u21920$); if, on the other hand, $\gamma >2$, the right-hand sides of (6b) are indeterminate $ C 23=0$ (or infinite if $ C 23\u2192 0 +$). Thus, we consider tree-like solutions for $0<\gamma \u22642$.

#### 2. Tree-like solution

Applying the tree-like solution to (8a), we have $ p 2\u2212 p 1= h 2/ C 12$. Imposing further a stationary nontrivial solution for (6a) yields $1= C 12 1 \u2212 \gamma \u27e8 ( p 1 \u2212 p 2 ) 2 \u27e9= C 12 \u2212 1 \u2212 \gamma \u27e8 h 2 2 \u27e9$ so that $ \u27e8 h 2 2 \u27e9= x 1 + \gamma = C 12 1 + \gamma $ or $x= C 12= \u27e8 h 2 2 \u27e9 1 1 + \gamma $. The tree-like branch is, therefore, given by $ B \u2227=( C 12, C 23)= ( \u27e8 h 2 2 \u27e9 1 1 + \gamma , 0 )$.

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

Similarly, for $ f y$, we only retain the term proportional to $ \xi y=1/x$; we use that $ \xi y\u22c5 x \u2212 \gamma = ( \u27e8 h 2 2 \u27e9 ) \u2212 1$ and obtain $ f y B \u2227=\u22122a\u2212a\u27e8 h 2 h 3\u27e9/ \u27e8 h 2 2 \u27e9$.

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 $\xi , \xi x, \xi y\u221dy$, thus resulting in terms of the form $ y p$ with some exponent $p$. If $p>0$, associated factors become zero as $y\u21920$; but if $p<0$, the associated term is indeterminate. Consider first $ g x\u221d y 2 \u2212 \gamma $. Since $y=0$ on the tree-like branch, we have $ g x B \u2227=0$ for $\gamma \u22642$. It follows that the associated Jacobian is an upper triangular matrix.

Next, consider $ g y$, for which the first term is $ y \u2212 1 \u2212 \gamma \xi 2\u221d y 1 \u2212 \gamma $ and the second is $ y \u2212 \gamma \xi \xi y\u221d y 1 \u2212 \gamma + y 2 \u2212 \gamma $. Provided that $\gamma <1$, $ g y$ is determinate for $y=0$, in particular, all terms go to 0 as $y\u21920$, and we obtain $ g y B \u2227=\u22121$. For $\gamma =1$, we use that the term with $ y 1 \u2212 \gamma = y 0\u21921$ as $y\u2192 0 +$ (all other terms containing $y$ are 0), and we get $ g y B \u2227=\u22121+ \u27e8 h 2 \u2212 h 3 ) 2 \u27e9/ \u27e8 h 2 2 \u27e9 2 / ( 1 + \gamma )$. Finally, considering $1\u2264\gamma <2$ the expression $ y 1 \u2212 \gamma \u2192+\u221e$ (as long as $\gamma <2$) as $y\u2192 0 +$ so that we have $ g y\u2192+\u221e$.

Note that letting $\gamma =1$ and using (A2a), we have $ \lambda 2 B \u2227= 3 \alpha 2 \u2212 1 \alpha 2 + 1$. This determines critical strength of fluctuations for $ \alpha 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.^{15}

To sum up, the tree-like branch is always stable for $0<\gamma <1$; for $\gamma =1$ only if $\alpha < \alpha c$; and for $1<\gamma \u22642$, it is “hyper-repellent.”

#### 3. Triangular solutions

##### 1. Case *γ* = 1

^{15,16}

##### 2. Case *γ* ≠ 1

##### 3. Bifurcations varying *γ* for different length *L*_{23}

We complete the consideration of the system with three nodes by observing response to varying $\gamma $ is varied while keeping $\alpha $ 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 $\gamma $ decreases through a critical value. Solutions are obtained numerically from (28) and (29) by tuning $ C 12$ such that $\alpha =1$ results for $\gamma $ and $ L 23$.

### B. Parameter scan on an empirical vascular network

We analyze the stationary solutions of the model on the empirical vascular network described in Sec. II C. While keeping fluctuations at maximum ( $\alpha =1$, single moving sink), we vary the cost exponent $\gamma $ and observe the conductances after settling to stable stationary solutions of (4). We observe that for $\gamma =1$, all vessel segments have positive conductivity, $ C i j>0$. As $\gamma $ 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 $\gamma $-dependence of five conductances involved in bifurcations earliest, i.e., at relatively large $\gamma $.

As $\gamma $ 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 $\gamma \u21920$ 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 $\gamma =0.5$^{14} does not reproduce the empirical network with all positive conductances. However, choices of $\gamma $ 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, $\alpha =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 $\gamma $ 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.

## IV. DISCUSSION

The model for vascular network adaptation by Hu and Cai^{14} 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<\gamma <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 $\gamma =1$ than the previously suggested^{14} value $\gamma =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 $\gamma =1/2$, a value of zero conductance is assigned to several vessel segments, incompatible with the existence of these connections. Increasing $\gamma $, 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^{18} 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.^{14} In the latter, each sink $i$ independently opens with probability $p$; otherwise, it fully closes. The case of maximal fluctuations ( $p\u21920$), however, is that of a single moving (open) sink and, thus, identical to a maximum fluctuation amplitude ( $\alpha =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^{17} and elsewhere.^{21} 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.^{22} Large scale whole brain vascular networks measured in rodents have recently become available.^{11} 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.^{5}

## ACKNOWLEDGMENTS

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

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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 AVAILABILITY

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

### APPENDIX: “POWER” TERMS

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

#### 1. Periodic drive

#### 2. Stochastic drive

## REFERENCES

*et al.*, “

*Progress in Industrial Mathematics at ECMI 2018*(Springer, 2019), pp. 147–155.

*et al.*, “