Resonances in the response of fluidic networks inherent to the cooperation between elasticity and bifurcations

A global response function (GRF) of an elastic network is introduced as a generalization of the response function (RF) of a rigid network, relating the average flow along the network with the pressure difference at its extremes. The GRF can be used to explore the frequency behaviour of a fluid confined in a tree-like symmetric elastic network in which vessels bifurcate into identical vessels. We study such dynamic response for elastic vessel networks containing viscous fluids. We find that the bifurcation structure, inherent to tree-like networks, qualitatively changes the dynamic response of a single elastic vessel, and gives resonances at certain frequencies. This implies that the average flow throughout the network could be enhanced if the pulsatile forcing at the network’s inlet were imposed at the resonant frequencies. The resonant behaviour comes from the cooperation between the bifurcation structure and the elasticity of the network, since the GRF has no resonances either for a single elastic vessel or for a rigid network. We have found that resonances shift to high frequencies as the system becomes more rigid. We have studied two different symmetric tree-like network morphologies and found that, while many features are independent of network morphology, particular details of the response are morphology dependent. Our results could have applications to some biophysical networks, for which the morphology could be approximated to a tree-like symmetric structure and a constant pressure at the outlet. The GRF for these networks is a characteristic of the system fluid-network, being independent of the dynamic flow (or pressure) at the network’s inlet. It might therefore represent a good quantity to differentiate healthy vasculatures from those with a medical condition. Our results could also be experimentally relevant in the design of networks engraved in microdevices, since the limit of the rigid case is almost impossible to attain with the materials used in microfluidics and the condition of constant pressure at the outlet is often given by the atmospheric pressure.

RDMT, 0000-0001-6078-0721; ECP, 0000-0002-4688-6922 A global response function (GRF) of an elastic network is introduced as a generalization of the response function (RF) of a rigid network, relating the average flow along the network with the pressure difference at its extremes. The GRF can be used to explore the frequency behaviour of a fluid confined in a tree-like symmetric elastic network in which vessels bifurcate into identical vessels. We study such dynamic response for elastic vessel networks containing viscous fluids. We find that the bifurcation structure, inherent to tree-like networks, qualitatively changes the dynamic response of a single elastic vessel, and gives resonances at certain frequencies. This implies that the average flow throughout the network could be enhanced if the pulsatile forcing at the network's inlet were imposed at the resonant frequencies.
The resonant behaviour comes from the cooperation between the bifurcation structure and the elasticity of the network, since the GRF has no resonances either for a single elastic vessel or for a rigid network. We have found that resonances shift to high frequencies as the system becomes more rigid. We have studied two different symmetric tree-like network morphologies and found that, while many features are independent of network morphology, particular details of the response are morphology dependent. Our results could have applications to some biophysical networks, for which the morphology could be approximated to a tree-like symmetric structure and a constant pressure at the outlet. The GRF for these networks is a characteristic of the system fluid-network, being independent of the dynamic flow (or pressure) at the considerations to apply the model to elastic vessel networks. In §4, we introduce the global RF for treelike symmetric elastic networks. In §5, we describe the two tree-like network morphologies that are used in this work. In §6, we find that the GRF is independent of the dynamics of the inflow, for networks that have constant pressure at the outlets, making the GRF a good quantity to study the network's dynamics. In §7, we show that the bifurcation structure of tree-like elastic networks causes the GRF to have resonances, which do not exist for rigid networks, nor for single elastic vessels. This implies that the flow magnitude across the network could be enhanced at certain frequencies due to the cooperation between the bifurcation structure and the elasticity of the network, via pulsatile forcing. We do a systematic study varying the networks' elasticity, and find features that are common to different network morphologies, and features that are morphology dependent. In §8, we present an analytical study of a single elastic bifurcation that demonstrates the emergence of the resonant behaviour. We present our conclusion in §9.

Flow through a single elastic vessel
To study the dynamics of each elastic vessel forming the network, we use the generalized Darcy's elastic model (GDEM) [41], which assumes the following: First, a generalized Darcy's model for rigid cylindrical vessels, which gives a linear relation between flow,q, and pressure gradient, @p=@x, in frequency domain, is valid locally, that is, at any point x along the flow direction,q where K(ω) = −(η/iωρ)[1 − 2 J 1 (βr)/βrJ 0 (βr)] stands for the dynamic permeability of a Newtonian fluid, which is a measurement of the resistance to flow; J 0 and J 1 are Bessel functions of the first kind of orders 0 and 1, respectively; β 2 = iρω/η, r is the average radius of the elastic vessel; η is the fluid viscosity; and i ¼ ffiffiffiffiffiffi ffi À1 p . Hats over quantities denote their Fourier transforms. Flow has been approximated as a velocity averaged over the average cross-sectional area, A, times this area.
Second, a linear relationship between pressure and flow gradient, coming from mass conservation and a Hooke-like linear model for vessel elasticity, is valid locally, namelŷ where C = 3πr 3 /2Eh is called the vessel compliance, h is the vessel wall thickness, E is the Young's modulus of the vessel given by where c is the pulse wave velocity. The limit of a rigid vessel is obtained when the vessel compliance vanishes, that is, when C → 0. Combining equations (2.1) and (2.2), one obtains a harmonic oscillator equation for the pressure in frequency domain, namely, where κ 2 = iωCη/AK(ω). Solving equation (2.4) for a single vessel of length l, with knownp in andp out for the inlet and outlet pressures on the vessel extremes, respectively, yields an analytical expression for the pressure along the vessel in frequency domain given bŷ p(x) ¼p in cos (kx) þp out Àp in cos (kl) sin (kl) sin (kx), (2:5) which, when substituted in equation (2.1), provides the flow along the vessel in frequency domain

Flow through an elastic vessel network
GDEM [41] considers elastic networks of tubes in which a viscous fluid flows subject to a time-dependent inlet flow (or pressure). In this paper, we focus on tree-like symmetric networks, in which cylindrical vessels bifurcate into two identical vessels, as illustrated in figure 1 for a 4-level network. The model deliberately excludes trifurcations or higher order branching, since bifurcations are the most probable and common forms of branching structures in nature. Articles in the literature that classify normal vascular networks consider only bifurcations [42][43][44]. More complex branching morphologies, such as trifurcations, are present in non-hierarchical pathological vasculatures, such as the ones irrigating tumours [45]. Tree-like networks will have two types of vessels: the ones in which boundary conditions are imposed, namely vessels that take inlet flows or pressures as boundary condition, or vessels in which outlet pressures are imposed; and vessels that are internal to the network. For these ones, the pressures at the vessels' extremes have to be determined as part of the solution.
Flow conservation is imposed at the nodes ( points of bifurcation) as well as equal pressures on the extremes of the vessels connected to a node. For instance, for a vessel at level j, with length l j , that bifurcates into vessels at level j + 1, with lengths l j+1 , flow conservation at the node joining these three vessels, requires the output flow at level j,q j (x ¼ l j ) to be equal to twice the input flow at level j + 1, q jþ1 (x ¼ 0), that is,q j (x ¼ l j ) ¼ 2q jþ1 (x ¼ 0). This allows one to write an equation, by means of equation (2.6) (or equation (2.8)), that involves: pressure (or flow) at the entrance of vessel at level j, pressure at the node, and pressures at the outlet (x = l j+1 ) of vessels at level j + 1. For each node in the network, there will be an equation coming from flow conservation at that node. The result is a system of equations for the pressure at the nodes, that can be solved, analytically or numerically, using the methodology described in the literature [41].
Once the pressure at the nodes is known, pressure and flow as a function of the position, x, along any vessel of the network, can be obtained by making use of equations (2.5)-(2.8). For a rigid network, the fact that flow is constant along the network allows one to describe the network behaviour in terms of a dynamic RF, χ(ω), which is a measure of how much flow, q 1 (t), goes through the network, given a pressure difference, Δp(t) = p out − p in , between the network inlet and the network outlet [25]. In frequency domain, this can be written aŝ where L ¼ P i l i , is the sum of the vessels' length, l i , at each level i. χ has been computed for a tree-like symmetric rigid network by means of an electrical analogy, adding resistances l i /A i K i in series or parallel, according to the network morphology [25], giving the following result: In equation (4.1),q 1 is the flow at the inlet of the network, but given that flow is constant along a rigid network, it is also the total flow at any network level, i, containing 2 i−1 vessels, that is,q 1 ¼ 2 iÀ1q i , whereq i is the flow in each of the vessels at level i.
For an elastic network, on the other hand, flow is not uniform along the network. Hence, we define a global response function (GRF), that would give information of the flow, integrated throughout the network, that is, the flow averaged along the flow direction, 〈Q〉 x (t). A GRF for a tree-like symmetric elastic network is defined, in frequency domain, as [38,39] x global ; À hLhQi x Dp , ( 4 :3) where Hereq i (x) is the flow at position x in one of the vessels at level i; 2 iÀ1q i (x) is the total flow at position x in level i; and is the average flow along a single vessel at level i. The second equality in this expression is obtained using equation (2.1). Dp ¼ P N i¼1 Dp i . Note that for a rigid network, hQi x reduces toq 1 , which is the flow at point x, of the network's cross-sectional area. In the limit of a rigid network, equation (4.3) reduces to equation (4.1), with χ global = χ, given by equation (4.2).

Network morphologies
In order to consider the impact that network morphology has on the dynamics, we study two different tree-like symmetric networks.
One of them, already used in previous works [25,34], approximates the morphological properties of the vasculature of a dog. The vessel morphology of the dog's circulatory system is described in detail in the literature [46], and we have approximated the number of vessels, lengths and radii to the closest numbers that satisfy branching into identical vessels (table 1). For convenience in the discussion, we will refer to this network as the 'dog's network'.
We also analyse networks whose vessels follow Murray's Law [49], since physiological studies have validated the agreement of Murray's Law in vascular systems of mammals [50] and even in networks that transport water in plants [51]. According to Murray's Law, when a vessel of radius r p bifurcates into two identical vessels, the next level vessels are narrower, with radii r d , that obey r 3 p ¼ 2r 3 d . Therefore, the radii of vessels at level i, are related to the inlet vessel radius, r 1 , by

Independence on boundary conditions
For a network response like the one in equation (4.3) to be useful, it must prove to be independent on boundary conditions. We limit the study to the case of networks whose pressure is constant at the outlet, and consider three different inflows as incident boundary condition. Namely, an aortic physiological flow measured in vivo [52], a simple time-dependent one-mode cosine function, and a simple pulse consisting of a delta function at time zero. We have chosen a constant pressure at the network outlet of 4.2 kPa for the three calculations, a value that is consistent with physiological situations. Following the methodology sketched out in §3, and explained in detailed in reference [41], we have computed flows and pressures all along the network. We have used them to compute the RF defined in equation (4.3). We have verified that the RFs of the network, using these three different inflows, and constant pressure at the network outlet, are identical. Since, by construction, the inflow consisting of a simple pulse in time domain, is constant in frequency domain, it has less noise than the other signals, and it is the one reported in all figures. We have also verified that halving or doubling the inflows or the outlet pressures do not alter the RF. Figure 2 shows the absolute value of the RF for two different networks. The first thing to notice is that both responses are non-monotonic functions of frequency and therefore present resonances (maxima at  [25], and based on the anatomical measurements collected in [46]. Typical dimensions of vessel 1 are those of the aorta, typical dimensions of vessels 2-5 are those of large arteries; of vessels 6-9 are those of main arterial branches, of vessels 10-11 are those of terminal branches, of vessels 12-25 are those of arterioles and of vessels 26- certain frequencies). This implies that the flow magnitude would be enhanced at certain frequencies of the pulsatile forcing. This enhancement comes from the cooperation between the bifurcation structure and the elasticity of the network via pulsatile forcing. Elasticity by itself does not cause resonances in the flow of Newtonian fluids subject to pulsatile forcing. That is, for a single vessel, the RF, χ global , given by equation (4.3), decreases monotonically with frequency (this will be explicitly shown in §8). The structure of bifurcations is necessary to have resonances. There is an interplay between network geometry and network elasticity to give the exact location of resonances. In figure 2, we can observe that the main resonance (the one with the highest value of the RF) for Murray's network occurs around ω = 2 rad s −1 , while for the dog's network occurs around ω = 4 rad s −1 . Values of the RF at resonance are 4 × 10 −11 m 4 and 5 × 10 −11 m 4 , respectively. These values are surprisingly close to each other, given that the steady-state network responses (at low frequencies, out of the range of figure 2) are five orders of magnitude different. This is because the outer, wider vessels, which are of a similar calibre in both networks, determine the high-frequency behaviour of the response. The elastic properties of wider vessels allow for better accumulation and release of fluid during oscillatory flow, causing the non-monotonic high-frequency behaviour of the RF.

Results
The presence of vessel elasticity changes qualitatively the dynamic response of rigid networks, for which there are no resonances. Figure 3 shows how the RF changes as the Young's moduli of the networks, E, are varied as multiples of E 0 . For both networks, Murray's in figure 3a and the dog's one in figure 3b, resonances shift to high frequencies as the system becomes more rigid. For both networks, resonance frequencies are proportional to the square root of the network's Young's moduli, that is, ω res ∼ (E/E 0 ) 1/2 (see §8). Resonance frequencies as a function of network's elasticity are shown in red in figure 4a for Murray's network. The continuous line, shown for reference, has a slope of 1/2. For the dog's network, frequencies for the first two maxima are shown in figure 4b. In between the points indicating the frequencies of the first and second maxima, there is a continuous line, shown for reference, with a slope of 1/2. As the network becomes more rigid, the first maximum takes over and becomes the main resonance of the system. The main resonance is plotted in red.
We have also found that the magnitude of the responses at resonance decreases for increasing network rigidity. This can be appreciated in figure 3a. Resonances disappear, for large network rigidities, when the values of the response function, at all frequencies, are below the value of the response at zero frequency (not shown). These are features of the RF that are independent of the network morphologies studied. Particular details of the RF are morphology-dependent. For Murray's network (figure 3a), there is a frequency region around the first resonance where the value of the response is larger than for a rigid network, while the high-frequency behaviour presents responses below the one of the rigid network and the low-frequency behaviour is quite similar to the one of the rigid network. By contrast, for the dog's network (figure 3b), vessel elasticity changes dramatically the low-frequency behaviour of the response, causing the network response to increase as a function of frequency. This is because inner, thinner vessels (with high resistance and small permeability), which are much smaller for the dog's network than for the Murray's network, determine the low-frequency behaviour of the RF. To accommodate the low-frequency behaviour of the response (given by the small vessels, with large resistance and low permeability) and the high-frequency behaviour of the response (given by the outer wider vessels, with low resistance and high permeability), the network's response increases as a function of frequency. This behaviour is qualitatively different from the monotonic decrease of the RF for rigid networks. Also, for the dog's network, the value of the response increases by several orders of magnitude, for most of the frequency spectrum, relative to the one of a rigid network. This implies that the maximum amplitudes of the overall flow, regardless of the frequencies involved in the pressure drop across the network, would increase significantly with respect to the rigid case.

Origin of resonances and scaling behaviour
As stated in §7, we have found resonance frequencies for the GRF of a Newtonian fluid in an elastic network, using a model that for a Newtonian fluid in a single elastic vessel does not exhibit resonances. This implies that the resonant behaviour is due to the structure of bifurcations, inherent to tree-like networks. In order to better see this, we compute analytically the GRF for a network consisting of a single bifurcation, and show the emergence of the non-monotonic behaviour as a function of frequency. The response for a Newtonian fluid in a single elastic vessel, with average cross-sectional area, A, can be obtained analytically from the equation for flow along the vessel (equation (2.6)) that, when integrated along the flow direction, gives   which gives an RF, χ global , defined in equation (4.3), that is equal to the area, A, times the dynamic permeability, K(ω), i.e. χ global = A K(ω), with K(ω) given by the analytical expression after equation (2.1) and, for a Newtonian fluid, has a monotonic decay as a function of frequency [26,31]. We now consider a network with a single bifurcation and compute its RF. We consider a vessel at level 1, that bifurcates into two identical vessels at level 2. Flow conservation at the bifurcation implies that outflow in vessel at level 1, is equal to the sum of inflows of vessels after the bifurcation, namelyq 1 (x ¼ l 1 ) ¼ 2q 2 (x ¼ 0). For vessel at level 1, the pressure at the vessel's entrance, is the pressure at the entrance to the system,p in , and the pressure at the exit is an unknown pressure at the node,p N . For vessels at level 2, the pressure at the entrance is the pressure at the node,p N , and the pressure at the vessels' exit is the output pressure of the system,p out . Using equation (2.6), in the equation for flow conservation at the node, we can write an equation for the pressure at the node,p N . This one is given bŷ : On the other hand (using equation (4.4)), the average flow along the network with a single bifurcation is which, using equation (8.1), can be written as ( 8 :4) withp N given by equation (8.2). Accordingly, the RF, defined in equation (4.3), is given by From this expression, it becomes clear why there is a non-monotonic behaviour as a function of frequency coming from the bifurcation. That is, the pressure at the node,p N , needed to compute χ global , and given by equation (8.2), contains non-monotonic sinusoidal terms in κ 1 l 1 and κ 2 l 2 , that are functions of frequency and of the mechanical properties of the vessels (see expression for κ after equation (2.4)). These terms were not averaged out when the integration along the flow was performed, as happened for a single vessel. Figure 5 illustrates the RF for a network, that consists of a single bifurcation, with the characteristics of the first three vessels for both the dog's and Murray's networks. We have chosen a zero pressure at the outlet. The results clearly show the non-monotonic behaviour coming from the bifurcation. For this particular example, : For this example, the RF, χ global , is explicitly independent from the pressure at the inlet. It is straightforward to prove that, for constant p out (t), χ global is independent of the pressure drop.
In order to understand the scaling behaviour observed in figure 4, we notice that the GRF has terms in K i κ i , which represent slowly varying modes of frequency, and terms in cos(κ i l i ) and sin(κ i l i ), which are rapid modes, that will determine the GRF extremes, and therefore the resonances. From the expressions for K, C and κ (after equations (2.1), (2.2) and (2.4)), and since κl is a non-dimensional quantity, we can obtain a characteristic frequency of each vessel in the system, given by ω i = (1/l i )(E i h i /ρr i ) 1/2 . As in many elastic systems, for instance, a forced harmonic oscillator, resonances appear when the forcing is made at the smallest characteristic frequency of the system. In this case, the frequency characteristic of the largest vessel in the network. This also explains why resonances shift to high frequencies as the system becomes more rigid (with higher values of Young's moduli).

Conclusion
A GRF of a tree-like symmetric elastic network is introduced as a generalization of the RF of a rigid network [38,39]. The GRF relates the network's flow, averaged along the flow direction, with the pressure difference at the network's extremes. It can be used to explore the frequency behaviour of a fluid confined in an elastic network. The GRF indicates which frequencies, involved in the dynamic pressure drop, maximize the magnitude of flow averaged along the flow direction. We have found resonance frequencies of the GRF for Newtonian fluids in elastic networks using a model that for a single elastic vessel, and for rigid networks, does not give resonances, and proved that this resonant behaviour is due to the cooperation between elasticity and bifurcations.
Some of the features of the GRF are common to networks of different morphologies, for instance, for all networks, resonance frequencies shift to high frequencies as the system becomes more rigid. For all of them, responses at resonance decrease for increasing network rigidity. Particular details of the RF are morphology-dependent. For example, in the dog's network studied here, vessel elasticity changes dramatically the low-frequency behaviour of the GRF, causing this one to increase as a function of frequency. This behaviour could be experimentally important for certain networks engraved in microdevices, since the limit of the rigid case is almost impossible to attain with the materials used in microfluidics.
For networks in which pressure is constant at the outlets, the GRF is characteristic of the system fluidnetwork, and independent of the dynamics of the inflow and of the value of pressure at the network's outlet. It might therefore represent a good quantity to differentiate healthy vasculatures from those with a medical condition. Abnormalities in large vessels could possibly be observed in the high-frequency behaviour of the GRF, while abnormalities in small vessels would in principle be observable in the lowfrequency behaviour of the GRF. Whether or not this quantity might be clinically relevant to discriminate vasculatures with a medical condition, from those of a control group, is yet to be explored.
Our methodology could also be applicable to the domain of microfluidics, where branched symmetric structures are often engraved in microchips whose materials range from elastomeric to rigid. For a possible experimental verification of our results, it would be worth recalling that, for given values of the elastic and fluid parameters of a microfluidic device, one can always attain the linear flow regime by decreasing the amplitude of the dynamic pressure drop.