Transitions from trees to cycles in adaptive flow networks

Transport networks are crucial to the functioning of natural and technological systems. Nature features transport networks that are adaptive over a vast range of parameters, thus providing an impressive level of robustness in supply. Theoretical and experimental studies have found that real-world transport networks exhibit both tree-like motifs and cycles. When the network is subject to load fluctuations, the presence of cyclic motifs may help to reduce flow fluctuations and, thus, render supply in the network more robust. While previous studies considered network topology via optimization principles, here, we take a dynamical systems approach and study a simple model of a flow network with dynamically adapting weights (conductances). We assume a spatially non-uniform distribution of rapidly fluctuating loads in the sinks and investigate what network configurations are dynamically stable. The network converges to a spatially non-uniform stable configuration composed of both cyclic and tree-like structures. Cyclic structures emerge locally in a transcritical bifurcation as the amplitude of the load fluctuations is increased. The resulting adaptive dynamics thus partitions the network into two distinct regions with cyclic and tree-like structures. The location of the boundary between these two regions is determined by the amplitude of the fluctuations. These findings may explain why natural transport networks display cyclic structures in the micro-vascular regions near terminal nodes, but tree-like features in the regions with larger veins.


INTRODUCTION
Network structures are found in all of our everyday life, ranging from social interactions over technological infrastructure to natural systems. Networks serve vital functions on microscopic to macroscopic length scales, ranging from proteins, DNA, cells, organs and organisms [1]. An important function of networks is to transport people, goods, metabolites, and information among other [2][3][4][5][6][7]. While in technology, transport networks are relatively rigid and yield limited adaptivity, biological transport networks are capable of ensuring robust flow and operate satisfactorily over a vast range of parameters to prevent operational failure even under extreme conditions.
An example for one of the most advanced transport networks is the mammalian vascular network. Every second, the vasculature is without interruptions serving regions of the brain despite ever changing neural activity [8] or changing demands in other tissue [9]. The mammalian vasculature, composed of bifurcations (nodes) and vessels (links), is highly adaptive because vessel diameters (weights) dynamically adjust to changes in flow properties such as pressure and shear stress via a variety of vessel response mechanisms [9]. Several biophysical models have investigated the response mechanism of single vessels [10][11][12][13]. Here, we aim at understanding how an adaptive flow network may respond on a global network level. Notably, the mammalian vasculature forms a complex network [14,15] displaying both tree-like [16] and cyclic motifs [17] which also are found in the leaves of trees or in power grid networks [18].
Recent work investigated optimal topology of flow networks from the perspective of their energy efficiency, damage resilience, or cost of repair [19][20][21][22][23], features which may be argued to have evolved over long time scales. From the perspective of optimization theory, it is interesting to note that depending on the shape of the cost function, tree-like or cyclic structures may be more effective [24][25][26][27], leading to phase transitions between tree-like and cyclic structures. Moreover, cycles not only confer redundant structures and thus improve damage resilience (i.e., another cost function), but may also be more favorable in the presence of fluctuations [20,28].
However, while research on adaptive (co-evolving) networks from the perspective of theoretical physics is under active development, the understanding of adaptive flow networks in particular remains relatively un(der)explored in the field of network theory [29]. Thus, we consider the dynamic stability of particular network configurations given that vessels may adapt their diameters slowly over relatively short time scales. Specifically, we wish to address the following questions: given a flow network with dynamically adapting weights, (i) what network configurations (weights) are dynamically stable, and (ii) if non-uniform flow (load) fluctuations are present in the network, how far do these fluctuations affect adaptation into network regions where fluctuations are absent? In other words, when spatially inhomogeneous load fluctuations are present, does the network partition into clusters exhibiting tree-like motifs and cycles?
We wish to establish a fundamental understanding of the possible dynamics and gain insights from the perspective of network theory, nonlinear dynamics and physics, rather than of specifically physiological aspects. Thus, to obtain answers to the above questions, we dispose of the mathematical intricacies inherent to solving biophysically detailed models. Instead of building on physiological models of blood vessel changes which occur via acute responses in tone or the slower remodeling of the vessel [10][11][12][13][30][31][32], we defer to simple conceptual models of adaptive flow networks, based on basic physical principles, allowing for analytical tractability.
We assume that the network has one constant inlet (source) and many outlets (sinks) subject to load fluctuations. In the language of vascular physiology, we model a bifurcating arterial network where the only inlet is a feeding artery, bifurcating in a tree-structure to the terminal nodes interfacing via capillaries to the venous network. At this interface, changing supply demands constitute load fluctuations which are rapid compared to the adaptive network dynamics. A similar (though physiologically different) situation is seen in (real) trees, where the stem feeds the tree with water and nutrients at a more or less constant rate. In the leaves, so-called stoma evaporate the sap, and their periodic opening and closing correspond to fluctuating sinks [28,33]. Inspired by these natural networks, we wish to address two questions: How strong do load fluctuations need to be so that cyclic shunts (loops) 1 emerge in the network that break the topology of a spanning tree? How far do these cyclic shunts reach into the tree toward the feeding vessel, so that a nonuniform network structure emerges, divided into two subgraphs, one tree-like and the other with cycles?
The article is organized as follows. In the next section, we introduce the model and simulations that we study. Sections 3 and 4 discuss the dynamics on simple network motifs (one source and one sink, and one source and two sinks, respectively). In section 5 we investigate the emergence of cyclic structures in a larger network with one constant source and many fluctuating loads, which we conclude in the Discussion in section 6.

Network Structure
V denotes the set of nodes of the network with N = |V| < ∞ and A ⊆ N × N the set of edges. The edges are bidirectional, so (i, j) ∈ A implies (j, i) ∈ A. Each node is assigned a pressure p i . The edge flow is f ij > 0 from node i to j. Furthermore we assume that the network is resistive and linear, i.e., it is Ohmian with f ij = C ij (p i − p j ), where an edge carries the property of a conductance between nodes i and j with C ij = C ji > 0 only if (i, j) ∈ A.
Here, we study the three kinds of wirings illustrated in Figure 1: (a) one source and one sink, (b) one source and two fluctuating sinks, (c) a tree-like network with height H allowing for cross-edges on every bifurcation/branching level, l = 0, . . . , H, leading to cyclic structure. A cycle is a connected subnetwork of m nodes such that each node has exactly two neighbors.
To model sources and sinks in the network, we include nonzero nodal flows h i . Denote by S ⊂ V the set of sink nodes, |S| = n = 2 H . We define one source at the edge feeding the network, h 1 = 1, and n sinks with h i (t) < 0 at all leaves of the tree-like structure (i.e., where capillaries connect to the vein network). For all other nodes, h i = 0. Mass balance requires that k∈V h k (t) = 0 for all t ∈ R.

Mass Conservation
Mass conservation demands that edge flows, f ij from node i to adjacent nodes j, and nodal flows, h i , match the local accumulation rate at node i, dV i /dt, i.e., where ρ is the fluid density and V i is the vessel volume at node i. However, assuming that accumulation is nearly instantaneous or vessels are inelastic the nodal accumulation rate becomes negligible [15]. Mass balance then becomes Kirchhoff 's first law stating that which is re-written in vector/matrix notation by defining the nodal flow h : = (h i ) i∈V and the Kirchhoff matrix which is solved for p : = (p i ) i∈V .

Dynamically Adapting Conductances
To impose adaptive dynamics to the network, we postulate the generic ad-hoc law for the conductances: Thus, the first term on the right hand side induces growth proportional to the power dissipated along the edge, thus mitigating rising pressure differences by increasing the conductance along the edge. Thus, the network adapts itself toward minimizing power consumption. The last term prevents unlimited growth of the conductances.
where we drop the primes and omit the argument t from now on.

Fluctuating Sinks
We consider sinks with periodic and stochastic drive, compliant with k∈V h k (t) = 0. Sinks are assumed to fluctuate with a characteristic time scale T ∼ 1/ω. Periodic driving may be implemented for N = 2 with h 1,2 = ±a cos ωt and for N = 3 with h 1 = 1, h 2,3 = − 1 2 ± a cos ωt. For networks with N > 3 ( Figure 1C), we implement only stochastic driving 2 .
For stochastic driving, let (s k ) be a sequence where for each k ∈ N 0 , the random variable s k has support S and is distributed identically, uniformly and independently. Then for time t and sink i ∈ S, In plain words, this reflects the situation where at each point t in time, one of the sinks has higher load than the others; after each time interval of length T, the sink with higher load is again chosen uniformly at random. Independent of time, the source (root node) has h 1 = 1. For all other nodes j (neither source nor sink), h j = 0.
2 Note that periodic driving may be generalized to larger spanning trees with N = 2 L > 3 where L denotes the number of branching levels: h l = − 1 2 ± a l mod 2 cos ωt where l = 2 L−1 + 1, . . . 2 L are indices for the leafs of the tree. Thus, there are (2 L − 2 L−1 )/2 pairs of leaves that balance each other. To break the symmetry, we may permute the leaf indices l in h l .

Solving the Flow
Equations (5) and (6) are invariant with regards to timedependent pressure shifts, p k (t) → p k (t) + P(t) =: p ′ k (t). Thus, we may let P(t) : = −N −1 k∈V p k (t) and which we later use to obtain Equation (27). A general solution of Equation (3) is of the form p = p hm + p in ∈ R N , where p hm and p in solve the homogeneous and inhomogeneous problems, respectively. -Each row in the Kirchhoff matrix K ∈ R N×N in Equation (3) has sum zero which implies that (1, . . . , 1) ∈ ker (K). Since rank(K) = N − 1 (no isolated nodes), we have that ker(K) = span((1, . . . , 1)) and p h = c · (1, . . . , 1), c ∈ R. In particular, by letting c : = P(t) at any given time t ≥ 0 we may choose a specific instance of the homogeneous solution during the simulation. -The inhomogeneous solution p in is determined exactly by solving the reduced system K r ·p in = h r , where the reduced Kirchhoff matrix K r ∈ R N−1×N−1 with full rank is given by deleting row l and column l in K, and h r ∈ R N −1 is constructed by removing entry l in h. Finally, p in is given by p in,l = 0 and complementing all other entries fromp in ∈ R N−1 .

Simulations
We use a simple Euler scheme with time step t = 10 −3 for numerically integrating Equation (5). With a given parameter value a, we run the dynamics for a duration of τ sim = 10 4 . We take averages and standard deviations of conductances over the time interval [τ sim /2, τ sim ]. For a parameter scan, an outer loop runs over values of driving amplitude a, starting at a = 1.0 and decrementing with a = 10 −3 . For a = 1.0, the conductance of each edge (i, j) ∈ A is initialized as C ij = 1. For a < 1.0, the integration is initialized with the conductance averages obtained in the previous run at parameter value a + a.
We have checked that the results are robust under variation of t and τ sim . In the limit of fluctuations much faster than adaptation, only the distribution values h i but not their temporal order determines the conductance values obtained. This fact and the symmetry of the tree topology under swapping sink nodes are used to speed up the simulations.

One Source and One Sink (2-Node Model)
Let us consider a system with only two nodes V = {1, 2} linked by the edge with conductance C 12 , as depicted in Figure 1A.
We assume there are fluctuations driving the system but no net pumping between the two nodes, hence h 1 t = h 2 t = 0. Invoking Kirchhoff 's first law, C 12 (p 1 − p 2 ) = h 1 (t), the model reduces to a single equation, Superficially, two equilibria seem feasible: (i) the trivial solution with C 12 = 0; and (ii) a non-trivial solution, defined via the condition C 2 12 = h 2 1 > 0. To obtain more insight into these solutions, let us assume that the drive h 1 (t) has a well-defined characteristic time scale, T = T(h 1 ). Then we may consider two limiting cases: slow driving (T ≫ 1) and rapid driving (T ≪ 1). For slow driving, the conductance is slaved to the driving, i.e., C 12 → h 1 (t) as t → ∞. For rapid driving, we may average the equations and seek solutions, C kl , averaged over rapid fluctuations with characteristic time scale T, and observe that C 2 12 → C 2 12 and C 2 12 → h 1 (t) 2 as T → 0. For slow driving h 1 (t) is quasi-stationary and for fast driving h 2 1 is constant. Therefore, determining stability of the two equilibria is straightforward, as we then simply may inspect the derivative of the right hand side of Equation (11), −(h 2 1 /C 2 12 + 1). Since h 2 1 > 0 and C 2 12 > 0, the non-trivial branch is always stable; however, C 12 = 0 corresponds to a singular (and unstable) solution.
For the case of periodic driving of the form where a ≥ 0 and ω = 2π/T, we may find an explicit (positive valued) solution for Equation (11), with the transient term T : = 2Me −2t a −2 (ω 2 + 1) where M is a constant determined by the initial condition. We may study two different limiting behaviors in the asymptotic limit t → ∞. For sufficiently fast driving (ω ≫ 1), As ω −1 → 0, fluctuations in C 12 become entirely negligible, so that C 12 → a/ √ 2 and C 12 → a/ √ 2 (note also 3 ). For slow driving (ω ≪ 1), the conductances are slaved to the driving and we have

One Source and Two Sinks (3-Node
Model, H = 1) We consider a motif with one source with h 1 = 1 and two fluctuating sinks h 2 and h 3 , as depicted in Figure 1B. For T → 0 (rapid driving) the sources obey h 2 t − h 3 t → 0, i.e., there is no net pumping between nodes k = 2 and k = 3. The conductances follow the dynamics given by Eight (quasi-stationary) solutions are conceivable where either C kl > 0 or C kl = 0 for every edge (k, l); however, only solutions with C 12 > 0 and C 13 > 0 are physically meaningful, and so, only the two solutions with C 23 = 0 or C 23 > 0 are feasible. Note that unlike the case of the 2-node motif, the branch C 23 = 0 is not singular anymore. The general solution for N ≥ 3 nodes is more intricate than for N = 2 nodes. Therefore, we limit our analysis from now on to the case of rapid driving, where T ≪ 1 is very small, and we consider only the asymptotic solutions where t → ∞. In analogy to the 2-node motif, we seek solutions, C kl , averaged over rapid fluctuations with characteristic time scale T. Considering that C ij changes on a slow time scale, C ij → C ij as T → 0, and therefore we may from now on use C kl and C kl interchangeably and omit · around conductances. For symmetry reasons, rapid driving implies that C 12 = C 13 . In this limit, the dynamics of the conductances is then constrained to a two dimensional symmetry manifold and effectively reduced to the two equations given by together with Equations (19)(20)(21), in which h k and p k fluctuate rapidly but C kl may be considered quasi-stationary. Pressures may be eliminated by observing the following equalities. Subtracting Equation (20) from (21), we have (24) and substitution of this expression into Equation (20) yields We first consider the "tree-like" solution branch with a nonconducting cross edge C 23 = 0. Stationarity for C 12 > 0 implies that (p 1 − p 2 ) 2 = 1, and mass balance (Equation 25) requires that p 2 − p 1 = h 2 /C 12 . Therefore, the tree-like branch is given by Second, we consider the non-trivial branch with C 23 > 0. Again, we first look for a stationary solution C 12 > 0 implying that 1 = (p 1 − p 2 ) 2 . Using the symmetry C 12 = C 13 and Equation (19) allows us to write 1 = C 12 (−p 2 − p 3 + 2p 1 ). Equation (10) implies −p 2 − p 3 = p 1 , and we have the constant pressure p * 1 = 1/(3C 12 ) in node k = 1. Next, Equation (10) implies also p 1 − p 2 = 3/2p 1 − 1/2(p 2 − p 3 ). Thus, 1 = (p 1 − p 2 ) 2 = 9/4 p 2 1 − 3/4 p 1 (p 2 − p 3 ) − 1/4 (p 2 − p 3 ) 2 = 9/4(p * 1 ) 2 + 1/4, since by assumption (C 12 + 2C 23 ) p 2 − p 3 = h 2 − h 3 = 0. Hence, we have p * 1 = 1/ √ 3 and therefore C 12 = 1/ √ 3. Finally, to determine C 23 we use Equation (24) and (p 2 − p 3 ) 2 = 1. We obtain For periodic driving we let h 2,3 = − 1 2 ± a cos ωt and obtain the solutions From B per △ we read off the critical value for the drive amplitude, a c = 1/ √ 6 at which C 23 > 0 becomes physically viable. We now study stability for B ∧ and B △ by considering Equations (22) and (23) which restrict dynamics to the two dimensional subspace defined by C 12 = C 23 . Using Equations (24) and (25) and eliminating pressures, we may compute the Jacobian which we evaluate for the two branches to obtain the corresponding eigenvalues for periodic driving, More generally, we find λ ∧ 2 = ( h 2 3 −2 h 2 h 3 )/ h 2 2 , but λ △ 2 yields an unwieldy expression. Thus, the two branches swap stability in a transcritical bifurcation at a c = 1/ √ 6, so that B ∧ is stable for a < a c and B △ is stable for a > a c , as shown in Figures 2A,B. In Figures 2C,D, we show simulations for periodic and stochastic driving with varying time scales, ω or T, respectively. For sufficient time scale separation (i.e., T 10 −1 ), the simulation results match the analytically obtained values on the stable branch to large precision. These results demonstrate that the simulated behavior converges to our analytical predictions as the characteristic driving frequency T −1 → ∞.

One Source and Multiple Fluctuating Sinks (Augmented Tree, H > 1)
Let us now turn to larger systems with more than two sinks and nodes representing intermediate branching points between Frontiers in Physics | www.frontiersin.org source and sinks. The networks considered are complete binary trees, augmented with cross edges directly at each branching, shown as dashed edges in Figure 1C. We consider the stationary conductance values under stochastic driving with a rapid characteristic time scale T ≪ 1. For a level l ∈ {0, 1, . . . , H} in a system with height H, we consider the dynamics for conductances of tree-edges, C (l,H) ∧ , and of cross-edges, C (l,H) − , see Figure 1C. By a (l,H) c we denote the critical driving amplitude above which cross-edges with nonzero conductances emerge. Figure 3 shows the asymptotic values of the conductances of cross-edges for fluctuations with varying amplitude a. Non-zero conductances in cross-edges appear with increasing amplitude in an order that follows the hierarchy of the tree, i.e., a (k,H) c > a (l,H) c for all 0 ≤ k < l < H. To illustrate, as we increase a, the cross-edges between sinks, l = H, are the first that begin to conduct. Then as we increase a even further, the next level, l = H − 1, obtains conducting cross-edges, and so on. Accordingly, the last transition to non-zero conductance for increasing amplitude a occurs in C (1,H) − between the root's child nodes.
One may speculate if all observed transitions from zero to non-zero conductance are of the same kind at different levels, only varying in the parameter value a (l,H) c at the transition and the slope in the supercritical regime. Figure 4 shows that this is indeed the case for the cross-edge at level l = 1. Plotting the conductances of the edges of the source node as a function a/a (l,H) c , we observe a perfect collapse of all data for varying H (and l = 1 fixed ), by rescaling a → a/a (l,H) c = :ã in systems of different heights H. So the top triangle including the source node of an augmented tree of general height H behaves in the same way as the H = 1 system with driving parameter a rescaled. Thus, the analytic results Equations (28) and (29) for H = 1 allow us to express the scaled conductances explicitly, The a-dependence of conductances on the other levels is more intricate, as is shown in Figure 5 for the case of a system with H = 6 levels. The derivatives of cross-edge conductances reveal detail not apparent in the coarser plot (Figure 3) of these conductances themselves. The slope is maximal at the transition (as the onesided derivative with a approaching a c (l, 6) from above). As a is increasing, the conductance curve becomes slightly flatter. This non-linear effect hints at dependencies between the levels, l. However, one may show that using the same rescaled amplitude parameter,ã, and then C (l,H)  cross edges bridging longer distances in the tree; such a network structure is shown in Figure 1D. Two further Supplementary Videos 3, 4 (297323_Martens_Video 3.MP4, 297323_Martens_Video 4.MP4) employ damaged versions of the latter structure in which roughly half of all cross edges have been randomly selected and removed. These simulations confirm the observation that the cross edges begin to conduct in order of the tree hierarchy, starting from sinks and moving toward the source.

DISCUSSION AND OUTLOOK
We have studied how fluctuating loads affect the re-configuration of vessels in an adaptive flow network. To do this, we have introduced a minimal model, consisting of a resistive network with conductances on the edges which adapt dynamically toward minimizing fluctuations in the network. To be explicit, the conductance is up-regulated with the pressure drop squared (power dissipated), but down-regulated with a factor proportional to the conductance. Although adaptation adheres to local rules they also experience a global feedback through the coupling via the flow network. An important question, which also appears in the context of synaptic plasticity [35], is then: what are the (stable) equilibrium configurations of conductance in the network? Furthermore, we introduced load fluctuations by including fluctuating sinks in certain network nodes. For a characteristic time scale T of the fluctuations much more rapid than the time scale of adapting conductances, T ≪ 1, we can treat fluctuations in terms of their averages over time, i.e., their amplitudes. Considering both periodic and stochastic fluctuations, we analytically and numerically investigated small and large networks to determine how their equilibrium conductance configurations depends on the amplitude of the load fluctuations, a. In particular, we used this model to investigate how far into the network the fluctuations would be able to induce re-configurations.
First, we investigated two very simple network motifs regarding the space of equilibrium configurations and their stability. The dynamics for the motif consisting of one fluctuating sink and source can be solved exactly ( Figure 1A). The conductance connecting the two nodes is non-zero and stable for all fluctuation amplitudes with a > 0. The triangular motif with one constant source and two fluctuating sinks ( Figure 1B) can also be solved analytically, but exhibits a transcritical bifurcation at a critical drive amplitude a c : for sub-critical drive, a < a c , the tree-like solution branch without cross edge is stable; for super-critical drive, a > a c , the cyclic structure (triangle with non-zero tree-edge and cross-edge) is stable, see Figure 2. Numerical simulations demonstrate that the system behavior asymptotically approaches the solutions obtained for the limit of rapid driving (T → 0) as T → ∞.
Next, we ran numerical simulations of larger tree-like networks with a constant source at its root and fluctuating loads in the leaves. Using an initial configuration C ij | t=0 = 1 for all edges (i, j) ∈ A, we then investigated into which quasistationary configuration the network settles. To investigate the system behavior analytically, we studied two scenarios. First, we assumed a simplified network structure, i.e., a complete binary tree with cross-edges only at each branching, see Figure 1C and Supplementary Video 1. Conductances C (l,H) − of cross-edges undergo a transition at a critical driving amplitude. We find C (l,H) − > 0 for a > a c (l, H), and C (l,H) − = 0 for a ≤ a c (l, H). The critical driving amplitudes a c (l, H), depending on the edge level l and height H of the tree-like structure, are hierarchically ordered so that the transitions to non-zero cross-edge conductances appear successively in descending order of l as the amplitude grows. Thus, we observe that as the amplitude is increased cycles first emerge near the fluctuating sinks and then spread toward the root with the constant source. Second, we compared these observations with the behavior in a generalized network topology. In this second scenario, we ran numerical simulations of networks lacking intrinsic tree structure, i.e., where all cross-edges are allowed, see Figure 1D and Supplementary Video 2. In this case, the threshold for an edge (i, j) between sinks i and j has a dependency on the length of the tree path between i and j. The longer the tree path (distance on tree), the lower its threshold. This means that-as the amplitude is increased-minimal distance short-cuts for longer paths form before the cross-edges from the first scenario. In both scenarios, however, the network is partitioned into two regions with treelike motifs near the source at the root and cycles closer to the fluctuating sinks in the leaves of the tree. Transitions get more complicated when randomly chosen cross-edges are topologically deleted (Supplementary Videos 3, 4).
Other work has investigated the appearance of cyclic structures (loops) by discussing flow networks via optimization in terms of minimizing dissipative losses under continuous reconfiguration of the flow conditions [19,20,28]. While this point of view may be motivated by evolutionary principles and yielded many interesting insights regarding the critical emergence of cycles, we were interested in formulating a minimal dynamical model [36]. Contrasting previous studies, we did not assume spatially uniform fluctuations and investigate conditions under which cyclic structures form in general; rather, we investigated if partitions between tree-like and cyclic structures would form and where their boundaries lie. Many studies in vascular physiology focused on the dynamics of single vessels and have complex biophysical models including a large degree of complexity prohibitive to mathematical analysis, and only few have computationally investigated dynamics in networks [30,32]; here, we tried to systematically address stability of such vessel models from a mathematical perspective. The study on simple network motifs has been insightful and complemented the numerical findings that we have made for larger trees. Further research may address equilibrium configurations and their transitions (a (l,H) c ) from tree-like to cyclic structures in network structures shown in Figures 1C,D, and attempt to find complete solutions via mathematical analysis and scaling arguments. For the case of more general network topologies (see e.g., Figure 1D and for the case of heterogeneities, an open question remains whether multi-stable equilibria are possible. Furthermore, research should be conducted on reducing complex biophysical dynamics to simpler mathematical models, which could be investigated toward identifying classes of different dynamic behavior. Finally, generalizations to coevolving networks with edges being dynamically created/deleted may be considered [37].