Tradeoffs between energy efficiency and mechanical response in fluid flow networks
TTradeoffs between energy efficiency and mechanical response in fluid flow networks
Sean Fancher ∗ and Eleni Katifori Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA
Transport networks are typically optimized, either by evolutionary pressures in biological systems,or by human design, in engineered structures. In the case of systems such as the animal vasculature,the transport of fluids is not only hindered by the inherent resistance to flow but also kept in adynamic state by the pulsatile nature of the heart and elastic properties of the vessel walls. Whilethis imparted pulsatility necessarily increases the dissipation of energy caused by the resistance,the vessel elasticity helps to reduce overall dissipation by attenuating the amplitude of the pulsatilecomponents of the flow. However, we find that this reduction in energy loss comes at the priceof increasing the time required to respond to changes in the flow boundary conditions for vesselslonger than a critical size. In this regime, dissipation and response time are found to follow a simplepower law scaling relation in single vessels as well as hierarchically structured networks. Applyingbiologically relevant parameters to the model suggests that animal vascular networks may haveevolved to maintain a minimal response time over lesser dissipation.
I. INTRODUCTION
The material transport network is a ubiquitous modelthat can be used to describe a large array of naturallyand artificially occurring systems. The various proper-ties of such networks are frequently either thought to beadapted through natural selection [1, 2] or specifically de-signed [3–5] to optimize certain aspects of the local flowof material or the global performance of the entire net-work. One particularly common optimized quantity isthe rate at which the energy to maintain the flow is dis-sipated throughout the network. Minimizing dissipationfor a certain delivered flow can be constrained, however,by other costs to construct or adapt the network such asthe metabolic cost to maintain the vessels, the require-ment for resilience to damage, and more [6–11].One easily identifiable factor that complicates mini-mizing dissipation is pulsatility in the driving force thatmaintains the flow of material. For systems such as fluidundergoing laminar flow or electrical current traversinga linear resistor, the rate of energy dissipation per unitlength along a vessel is proportional to the square of thelocal current. In a pulsatile system with a well defined pe-riod, the time-averaged value of the square of the currentwill necessarily be greater than the square of the average,meaning the inclusion of any pulsatile components to theflow can only cause the time-averaged dissipation to in-crease. Thus, such linear systems utilizing a pulsatiledriving force will necessarily have greater total dissipa-tion over the course of a full period than an equivalentsystem utilizing a continuous, steady flow pump with thesame average flow rate.Despite this drawback, many transport networks relyon pulsatile drivers to maintain the flow of material. Onecommon such example is the animal vascular network.In humans, blood pressure and flow have pulsatile com-ponents imparted by the periodic beating of the heart. ∗ [email protected] However, the amplitude of this pulsatility has been ob-served to decrease as the blood is transported furtheraway from the heart [12, 13]. The damping of pulsatilitycan be attributed to the inherent resistance in the flow ofblood dissipating the energy contained within the kineticmomentum of the blood itself and the elastic energy dueto the nonzero compliance of arterial walls, resulting intraveling waves in flow and pressure attenuating [14–18].The culmination of these effects allows the flow of bloodto dissipate less energy than it would without dampingof pulsation, which is evolutionarily advantageous.Pulsatility is not absolutely damped, however, as itremains to some degree even at the capillary level [19].We might then wonder why evolution did not do awaywith pulsatility entirely, by developing ways to removeit altogether. One such obvious reason is that removingpulsatility entirely will require a vessel architecture thatis fundamentally incompatible with what can be built bi-ologically. In this work, we posit the existence of anothereffect, an increase in the network’s response time, whichmight prevent evolution from further decreasing pulsatil-ity and the energy necessary to maintain a certain meanblood flow.Vessel compliance in particular has been shown to af-fect the strength of pulsatile damping with stiffer, lesscompliant vessels leading to weaker damping [14, 20].While several theoretical and modelling studies have in-vestigated the effects of vessel compliance on pressurewaveforms and energetic costs in subsections of biologicalfluid flow networks such as the major arteries [14–18, 20–23], exactly how the effects of compliance might alter theglobal properties of the entire vascular network, particu-larly when the pulsatile driver relinquishes its periodic-ity to transition between different steady states, remainspoorly understood. Such an understanding could pro-vide deeper insight into the evolutionary pressures thatformed systems such as the animal vasculature as well asmore reliable models for generating artificial networks orprosthetic devices.Here, we investigate the effects of vessel complianceon the material flow in transport networks with time- a r X i v : . [ phy s i c s . b i o - ph ] F e b dependent boundary conditions. In Sec. II A we obtain aset of dynamic equations for the pressure and flow withina vessel by linearizing the Navier-Stokes equations forflow within an elastic, cylindrical tube [15, 16, 24, 25].The resulting equations are a special form of the tele-grapher’s equations [26], from which we construct net-works with well defined boundary conditions. We findthat while increasing the compliance of the vessels withina network does decrease the rate of energy dissipation,it also has the consequence of increasing the amount oftime required for the system to adapt to changes in theboundary conditions, such as a change in the heart rateof a vascular system. These effects follow a simple powerlaw relation for vessel sizes above a critical, compliancedependent value. We show that this result holds for sin-gle vessels, for which the equations can be solved analyt-ically, as well as large networks, for which we use numer-ical integration in Sec. II B. Finally, we obtain a general-ized method for approximating the time scale over whichthe flow and pressure will adapt to a given set of changesin the boundary conditions of a network. Our work high-lights the importance of the response time, the time forthe network to adapt to the new pump flow conditions,as an important design consideration for networks com-posed of elastic vessels. II. RESULTSA. Single Vessel Mechanics
We begin by considering an incompressible fluid withdensity ρ and viscosity µ flowing through a cylindricalvessel (Fig. 1A). We assume the system is rotationallysymmetric so that the flow rate and fluid pressure dependonly on the axial and radial positions z and r . Denot-ing the axial and radial fluid velocities as u z ( z, r, t ) and u r ( z, r, t ) respectively and the fluid pressure as p ( z, r, t ),the incompressibility condition and Navier-Stoke equa-tion are [15, 16, 24, 25] ~ ∇ · ~u = ∂u z ∂z + 1 r ∂∂r ( ru r ) = 0 , (1) ~ ∇ p + ρ ∂~u∂t + ρ (cid:16) ~u · ~ ∇ (cid:17) ~u − µ ∇ ~u = 0 . (2)We now reexpress Eqs. 1 and 2 in terms of the totalvolumetric flow rate I ( z, t ) = R dA u z ( z, r, t ) and area av-eraged pressure P ( z, t ) = A ( z, t ) − R dA p ( z, r, t ), wherethe integration is performed over the cross sectional area, A ( z, t ), of the vessel at axial position z and time t (Fig.1B). We must first make the approximation that the term ρ ( ~u · ~ ∇ ) ~u is negligible as it is nonlinear in the flow veloc-ity. Eq. 2 now represents two separate linear equations,one for the axial component and one for the radial com-ponent. The angular component is trivial due to ourassumption of rotational symmetry. u r ( r,z,t ) u z ( r,z,t ) u φ = 0 I ( z,t ) P ( z,t ) A FIG. 1. Diagrams of 3D and 1D flow and pressure mechanics.A) Fluid flowing through a compliant cylindrical vessel withaxial and radial velocities, u z and u r , cause changes in theradius and cross sectional area of the vessel in both space andtime. Rotational symmetry enforces u φ = 0. B) The radialcomponent of the flow is integrated out so as to write thedynamics purely in terms of the volumetric flow, I , and fluidpressure, P , dependent only on time and the axial dimension.Current and pressure pulses can travel through the vessel withdynamics dictated by Eqs. 5 and 6, resulting in exponentiallydecaying pulse amplitudes. C) The flow dynamics resultingfrom Eqs. 5 and 6 can also be interpreted as the continuumlimit of a series of inductors and resistors connected in parallelto a ground through capacitors. We now integrate the incompressibility condition (Eq.1) and average the axial component of the Navier-Stokeequation (Eq. 2) over the vessel cross section to produceequations for ∂I/∂z and ∂P/∂z . By equating the radialfluid velocity at the vessel wall to the expansion rate ofthe wall, the second term in Eq. 1 can be shown tosimply become ∂A/∂t once this integration is performed.Similarly, once again using the laminar flow solution forthe fluid velocity allows the radial term in −∇ u z to bereexpressed as 8 πI/A once it is area averaged. The axialterm is negligible whenever the wavelength of any pulsestravelling through the fluid is significantly larger than thevessel radius, which is the regime we restrict ourselvesto here. The culmination of these manipulations andapproximations is to transform Eqs. 1 and 2 into ∂I∂z + ∂A∂P ∂P∂t = 0 , (3) ∂P∂z + ρA ∂I∂t + 8 πµA I = 0 . (4)Eq. 3 can be simplified by assuming that the vesselcross sectional area scales linearly with the fluid pressureas A ( z, t ) = A + cP ( z, t ), where c is the compliance of thevessel. Additionally, we make the assumption that thevessel cross section, A ( z, t ), does not significantly change( A (cid:29) cP ( z, t )) so as to allow the factors of A in Eq. 4to be sufficiently approximated by the constant A . Wecan now define the fluid inertia and the flow resistanceper unit length as ‘ = ρ/A and r = 8 πµ/A respectively.These three parameters, c , ‘ , and r , thus characterize thesystem and allow us to define three distinct derived pa-rameters: the characteristic length scale λ = 2( p ‘/c ) /r ,the characteristic time scale τ = 2 ‘/r , and the character-istic admittance α = p c/‘ . Reformulating Eqs. 3 and 4to be expressed in terms of these characteristic parame-ters produces the more symmetric looking versions λ ∂I∂z + ατ ∂P∂t = 0 , (5) αλ ∂P∂z + τ ∂I∂t + 2 I = 0 . (6)Eqs. 5 and 6 represent a simple form of the telegra-pher’s equations with spatially and temporally indepen-dent parameters [26]. The derivative terms create travel-ing waves of current and pressure while the existence ofthe resistive term 2 I in Eq. 6 causes the waves to expo-nentially decay, as depicted in Fig. 1B. These equationscan also be derived via an analogous transmission linecircuit with no shunt resistor, as depicted in Fig. 1C.With the dynamic equations for the current, I ( z, t ),and pressure, P ( z, t ), given by Eqs. 5 and 6, we nowinvestigate the rate of energy dissipation over a vesselof total length L . In particular, we consider a pulsatilesystem with a well defined period T and correspond-ing angular frequency ω = 2 π/T . This allows the cur-rent to be expanded into a Fourier series of the form I ( z, t ) = P n ˜ I ( n ) ( z )exp( inωt ), where the summation isover all integers. The same expansion can also be per-formed for the pressure. Assuming the current bound-ary conditions, ˜ I ( n ) (0) and ˜ I ( n ) ( L ), are known for all n , Eqs. 5 and 6 can be solved exactly with the excep-tion of ˜ P (0) , which must be determined by the choiceof gauge (see Appendix B). When written in terms ofthese Fourier components, the solution can be easily sep-arated into forward and backward moving waves that ex-ponentially decay over a distance of λ/ Re( k ( nωτ )) andtravel with a phase velocity of ωλ/ Im( k ( nωτ )), where k ( y ) = p iy (2 + iy ) with the principle root being taken.If we consider the symmetric case in which ˜ I ( n ) (0) =˜ I ( n ) ( L ) for all n , then the total time averaged dissipation n D i ss i p a t i o n S c a li n g F un c t i o n , f ( L / , n ) L /10 f V e ss e l S i z e , L / FIG. 2. Plot of dissipation scaling function (Eq. 7) as a func-tion of frequency for a variety of vessel lengths. The meandissipation is given by the dashed lines in the main plot andblack curve in the inset. Triangle markers in both plots alsodesignate the value of the mean dissipation for each vessellength shown. While short vessels experience large fluctua-tions in dissipation values, the mean dissipation monotoni-cally decreases with increasing vessel size. rate, h D i = T − R T dt R L dz rI ( z, t ), can be expressedas rL P n | ˜ I ( n ) (0) | f ( L/λ, nωτ ) where f ( x, y ) = k I ( y ) sinh (cid:0) xk R ( y ) (cid:1) + k R ( y ) sin (cid:0) xk I ( y ) (cid:1) xy (cid:16) cosh (cid:0) xk R ( y ) (cid:1) + cos (cid:0) xk I ( y ) (cid:1)(cid:17) . (7)and k R ( y ) and k I ( y ) are the real and imaginary partsof k ( y ) respectively (see Appendix C). The function f ( L/λ, nωτ ) is identified as the dissipation scaling func-tion for the chosen system of a single vessel with symmet-ric current boundary conditions. Specifically, it dictatesthe factor by which individual Fourier components of thecurrent deviate from the baseline dissipation of rL | ˜ I ( n ) | .We can examine the effects the various parametershave on the overall dissipation rate by observing theproperties of the dissipation scaling function, as shownin Fig. 2. Specifically, we are most interested in the ef-fects of varying the compliance, c , while holding fixed theresistance, r , to maintain a constant baseline dissipationand inertia, ‘ , to maintain a constant fluid density. Interms of the derived parameters λ , τ , and α , this is equiv-alent to holding τ and the product αλ constant. Since α is absent from the expression f ( L/λ, nωτ ) though, wecan simply focus on the effects of varying λ .While the dissipation scaling function as given in Eq.7 limits to a value of 1 whenever either of its argumentsare taken to 0, for nonzero frequencies it is seen to gen-erally take on smaller values when L/λ is increased, orequivalently when the compliance is increased. This isnot always true due to the existence of resonant and an-tiresonant frequencies, causing the dissipation to increaseor decrease for certain values of ωτ relative to the meandissipation. However, these variations around the meanbecome negligibly small for long vessels, L (cid:29) λ . Onelimit of particular interest is ωτ (cid:29)
1. In this regime, theoscillations in the fluid travelling through the vessel areof a high enough frequency that inertial forces dominateover resistive forces. The dissipation scaling function alsotakes on a relatively simple form due to the fact that k ( nωτ ) ≈ inωτ in this regime. Using this approxi-mate form in Eq. 7 and ignoring the trigonometric termsso as to focus on the mean value, the dissipation scalingfunction is seen to oscillate around tanh( L/λ ) / ( L/λ ) as nωτ is varied, as seen in the dashed lines of Fig. 2. Thiscan further be approximated as (
L/λ ) − if the vessel islong enough that tanh( L/λ ) ≈ rL to show that the total dissipation scaleslinearly with vessel length and is independent of λ forshort vessels as rLf ( L/λ, nωτ ) ≈ rL when L (cid:28) λ . Thisreflects the simple fact that the current does not ap-preciably change over the span of a short vessel whilethe total resistance of the vessel scales linearly with L .Conversely, the total dissipation becomes independentof L and scales linearly with λ for long vessels when rLf ( L/λ, nωτ ) ≈ rλ due to the fact that λ dictates thedistance over which the resulting waves are damped be-fore becoming negligible. These results are valid in thelimit of high frequences. For ω →
0, the dissipation isdominated not by the pulse wave but by the usual Pois-seulle flow viscous losses.With an understanding of how the compliance affectsdissipation in a single vessel, we now focus on the dy-namic properties. In particular, it has been previouslynoted in computational models that an increase in ves-sel compliance correspondingly causes an increase in thetime required to establish a new steady state when thevessel boundary conditions change [20]. To investigatethis, we consider the case in which the current at thestart of a vessel obeys I (0 , t ) = ˆ I Θ( t ), where Θ( t ) is theHeaviside step function. Additionally, as was done forthe determination of the dissipation scaling function wemaintain symmetry in the boundary conditions by alsoenforcing I ( L, t ) = ˆ I Θ( t ). With these boundary condi-tions set, we can derive an expression for the current atany point within the vessel and at any time (see supple-mental material). While the exact solution results in anintegral term that is not analytically solvable, it can beexpressed in the form I ( z, t ) = ˆ I (cid:0) − ζ ( z, t ) (cid:1) , (8)where ζ is a residual normalized current function whoseexplicit form is given by Eq. S54.This residual current function is derived from the sum-mation of a pair of current wavefronts travelling throughthe vessel from both sides with a velocity of λ/τ . Thesetwo wavefronts reflect off the ends of the vessel to form aseries of ever increasing number of images that must be summed together to calculate the residual current func-tion. While ζ can be expressed exactly by such a summa-tion of images, we can better understand how it behavesin time by taking a long time limit and neglecting allsubdominant terms. When this is done, we can obtainthe scaling law ζ ∝ exp( − t/ ( τ β ( L/λ ))), where β (cid:18) Lλ (cid:19) = Lλ ≤ π − r − (cid:16) πλL (cid:17) ! − Lλ > π . (9)The function β ( L/λ ) is identified as the response scal-ing function and dictates the time scale in units of τ overwhich the flow within the vessel responds to a changein the boundary conditions. This can be seen in Fig.3A, which shows the current at the midpoint of the ves-sel normalized by the size of the input shift, I ( L/ , t ) / ˆ I ,for a variety of different vessel lengths and with timerescaled by the response time τ β ( L/λ ). For short ves-sels the reflected wavefronts are highly nonnegligible andcause the current to undergo decaying oscillations aroundthe long-time steady state value of ˆ I . These reflectionstake significantly more time to traverse long vessels andas such decay too much to cause oscillations. The insetof Fig. 3A shows the magnitude of the residual current,1 − I ( L/ , t ) / ˆ I , decaying with time. The residual currentin long vessels is seen to decay almost perfectly exponen-tially while that in short vessels fluctuates around anexponentially decaying value.A similar analysis can be performed in the pulsatilecase in which I (0 , t ) = I ( L, t ) = ˆ I cos( ωt )Θ( t ). Since theresulting steady state created by these boundary condi-tions oscillates in time, we measure the response of thesystem by calculating the amplitude of the current ratherthan the current itself. At any given point within thevessel, this amplitude, Q ( z, t ), will converge to a singlevalue, Q ∞ ( z ) (see Appendix B), which is independent oftime and given by Q ∞ ( z ) = lim t →∞ Q ( z, t )= ˆ I (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) sinh (cid:16) L − zλ k ( ωτ ) (cid:17) + sinh (cid:0) zλ k ( ωτ ) (cid:1) sinh (cid:16) Lλ k ( ωτ ) (cid:17) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (10)We numerically calculate the amplitude at any finitetime by considering not only the current at that mo-ment but also its derivative. Specifically, we define Q ( z, t ) = [ I ( z, t ) + ( ω − ˙ I ( z, t )) ] / , where ˙ I is the par-tial time derivative of I . This definition of the ampli-tude is exact in the case of I ( z, t ) oscillating sinusoidallywith frequency ω . Fig. 3B shows the amplitude at themidpoint of the vessel normalized by the long time limitgiven in Eq. 10 for a variety of different vessel lengthsand a frequency of ωτ = π/
4. The inset then depictsthe magnitude of the difference between the amplitude t / ( L / )0.000.250.500.751.001.251.501.75 N o r m a li z e d M i dp o i n t C u rr e n t , I ( L / , t ) / I A = 0 t / ( L / )10 | I / I | L / = /16 L / = /4 L / = L / = 4 0 1 2 3 4 5 6 7 8Normalized Time, t / ( L / )0.000.250.500.751.001.251.501.75 N o r m a li z e d M i dp o i n t C u rr e n t A m p li t u d e , Q ( L / , t ) / Q ( L / ) B = /4 t / ( L / )10 | Q / Q | L / = /16 L / = /4 L / = L / = 410 Vessel Size, L /10 R e s p o n s e S c a li n g F un c t i o n , ( L / ) C Eq. 9= 0= /16= /4== 4 10 Response Scaling Function, ( L / )10 M e a n D i ss i p a t i o n S c a li n g F un c t i o n , f ( L / , ) D = 0= /16= /4== 4 10 V e ss e l S i z e , L / FIG. 3. Calculation of response time for a single vessel. A) When symmetric, step function boundary conditions are appliedto the vessel, the current at the midpoint is seen to undergo decaying oscillations if the vessel is short or simply slowly climbtowards the final steady state value if the vessel is long. Inset plot shows the magnitude of the difference between the currentand its final steady state value. See supplement for plot with the time axis not normalized by β ( L/λ ). B) For step functionpulsatile boundary conditions, the current amplitude is seen to similarly undergo decaying oscillations that become smaller yetlast longer as the vessel becomes longer. Inset plot also shows the magnitude of the difference from the final steady state value.C) By performing a linear fit to the data presented in the insets of A and B, the response time can be numerically extracted.Eq. 9 is seen to provide a very accurate prediction for long vessels. D) By comparing the response time calculated in C to themean dissipation, a simple power law relation is seen to exist between dissipation and response time for long vessels. and its long time limit value normalized in the same wayand shows that it approaches zero exponentially with theexact same response scaling function as the nonpulsatilesystem ( ωτ = 0), Eq. 9. This is due to the fact thatthe inclusion of pulsatility in the boundary conditionsonly affects the overall magnitude of the residual currentfunction, not the temporal scaling (see supplemental ma-terial).We can obtain a numerically derived estimate of β forboth the nonpulsatile and pulsatile boundary conditionsby performing a linear fit to the insets seen in Fig. 3Aand 3B. Fig. 3C compares Eq. 9 to this fitted valuefor a wide range of L/λ and ωτ values with good agree- ment. For long vessels in particular, L (cid:29) πλ , the re-sponse scaling function can be well approximated simplyby β ( L/λ ) ≈ L/πλ ) . This quadratic scaling is to beexpected, as for long vessels the current changes slowlyand the inertial term in Eq. 6, τ ∂I/∂t , becomes negligi-ble. When this happens, Eq. 5 can be differentiated withrespect to z and combined with Eq. 6 to produce the dif-fusion equation ( λ / τ ) ∂ I/∂z = ∂I/∂t , which in turncauses any steep current waveform to expand diffusivelyand time scales to increase quadratically with distance.For short vessels, the invariance of the response time tochanges in vessel length is simply a consequence of τ dic-tating the time scale over which the reflecting wavefrontsdecay and thus also being a minimum possible responsetime.The dissipation scaling function can be directly com-pared to the response scaling function in how they eachdepend on L/λ . Both scaling functions simply limit to1 for very short vessels, but for long vessels the dissi-pation and response scaling functions scale as (
L/λ ) − and ( L/λ ) respectively. This provides a powerful pre-diction: for a single vessel of sufficiently large length andfixed r and ‘ , the energy dissipation rate due to any pul-satile components of the current should scale as the in-verse square root of the response time of the current tochanges in its boundary conditions, as shown in Fig. 3D.Alternatively, if the compliance of the vessel is increasedto cause a decrease in the characteristic length scale λ ,greater wave damping, and less total dissipation, the cur-rent will necessarily require greater time to respond toany changes. Thus, a tradeoff is created wherein mini-mizing dissipation comes at the cost of slower response. B. Network Mechanics
We now consider a multitude of compliant, fluid carry-ing vessels interconnected to form a fluid transport net-work of nodes and edges. For bookkeeping purposes andwithout loss of generality we can assign a directionalityto each edge, e.g. edge e = ( µ, ν ) is traversed from µ to ν . We identify the location along each vessel with thevariable z , which is z = 0 at the node µ where the edgeis outgoing from, and z = L µν where the edge is incom-ing to. Each individual vessel is still considered to obeyEqs. 5 and 6, but each vessel may have its own uniquevalues for λ , τ , and α , which are z independent. Specifi-cally, we identify λ µν as the value of λ within the vesselthat begins at network node µ and ends at node ν withthe same index notation also being applied to all otherparameters and variables. Scalar quantities such as thecharacteristic length scale λ µν are independent of the di-rection traversed between the nodes and thus symmetricwith respect to an interchange of indices. For spatiallydependent quantities such as the pressure P µν ( z, t ), theorder of the indices implies the directionality of the edgeand P µν ( z, t ) = P νµ ( L µν − z, t ), so that P µν (0 , t ) canalways be identified with the pressure P µ ( t ) at node µ .Quantities which depend on the direction the edge is tra-versed such as I µν ( z, t ) change sign when the beginningand ending nodes are switched and are thus antisym-metric with respect to an interchange of indices, so that I µν ( z, t ) = − I νµ ( L µν − z, t ).The connectivity laws of the network are taken to betwo-fold: 1) pressure is continuous across networks nodesand 2) the total current being inputted into a node mustequal the total current flowing away from it through thenetwork. These are manifested mathematically by en-forcing that P µν (0 , t ) = P µ ( t ) ∀ ν ∈ N µ , (11) ~ H ( t ) FIG. 4. Diagram of toy network structure. Vessels repre-sented by red lines connect at nodes represented by black dots.Solid lines show base vessels that connect the inlet (far left)and outlet (far right) nodes to the central nodes, and dashedlines show looping vessels that allow fluid to flow betweenbranching generations. A current driver (blue) provides anarbitrary externally imposed current H ( t ) into the inlet nodeand out of the outlet node via external vessels (dotted line). X ν ∈N µ I µν (0 , t ) = H µ ( t ) , (12)where H µ ( t ) is the current being inputted into node µ by an external source and N µ is the set of all nodes con-nected to node µ by a single vessel. For example, thenetwork depicted in Fig. 4 has H µ ( t ) = 0 for all internalnodes while H inlet = H ( t ) and H outlet = − H ( t ). Thus,for each node there are two possible boundary conditionsto specify, P µ ( t ) and H µ ( t ), creating a total of 2 N possi-ble boundary conditions to specify for the whole network,where N is the number of nodes.However, there exists and interdependence between theset of P µ ( t ) and H µ ( t ) that reduces the necessary num-ber of specified boundary conditions to a subset of size N . This can be shown by considering the relationshipbetween P µ ( t ) and H µ ( t ). By expanding P µν ( z, t ) and I µν ( z, t ) into their respective Fourier series in time, Eqs.5 and 6 can be solved to express ˜ P ( n ) µν ( z ) and ˜ I ( n ) µν ( z ),the Fourier transformed vessel pressures and currents, aslinear combinations of ˜ P ( n ) µ and ˜ P ( n ) ν , the Fourier trans-formed node pressures (Eqs. B15 and B16). This inturn allows the second connectivity law to be rewrittenas ˜ H ( n ) µ = P ν L ( n ) µν ˜ P ( n ) ν , where L ( n ) is the network Lapla-cian matrix in Fourier space and is given by Eq. D2. Thiscreates a set of N linearly independent equations for eachnonzero frequency. For the zero frequency mode the ma-trix L (0) only has rank N − P (0) µ . Thus, once a gaugeis defined the space of undetermined variables per fre-quency is reduced from the full 2 N values of ˜ P ( n ) µ and˜ H ( n ) µ down to a subset of only N values.With the mechanics and necessary boundary condi-tions defined, we now construct a simple toy network totest how its properties compare to those derived for thesingle vessel. Fig. 4 depicts a hierarchical network inwhich the input and output vessels branch inward overseveral generations to meet at a central set of nodes.Solid lines represent the base vessels that allow for fluidto reach each of the central nodes while dashed lines rep-resent looping vessels that allow for greater mixing of thefluid flow. For now we consider a relatively simple casewith no looping vessels in which all remaining vessels areof equal length and each parent vessel branches into twodaughter vessels each with their own values of λ , τ , and α . The size of the daughter vessels are taken to obey a γ parent = ba γ daughter 1 + (1 − b ) a γ daughter 2 , (13)where each a in Eq. 13 is the radius of the cross sectionof the respective vessel, γ is the branching exponent, and b is the branching ratio. Thus, γ = 2 corresponds tobranching in which the total cross sectional area is pre-served. How the cross sectional area of the daughterscompares to that of the parent is relevant due to λ , τ ,and α each scaling linearly with this area when the prop-erties of the fluid and vessel wall material are held fixed(see Appendix A).We now investigate the properties of such networkscompared to those of the single vessel derived previ-ously. To begin, we note that in a single vessel anywavefront that isn’t part of a perfectly periodic signalwill travel a distance z in a time zτ /λ and decay asexp( − z/λ ). In a network, such wavefronts will neces-sarily split when they encounter branching nodes. Sincethese wavefronts travel with velocity λ/τ within each in-dividual vessel, the time to travel from one position toanother along a specific path S through the network canbe expressed as t S = R S dz τ ( z ) /λ ( z ), where the integralis over the path S . Similarly, the wavefront will decay asexp( − R S dz /λ ( z )). For a network such as the one shownin Fig. 4, we can use this decay function to generalize theexpression L/λ , the nondimensionalized size of the singlevessel, to the path dependent σ S = R S dz /λ ( z ), where S is any path from the far left node to the far right nodethat does not backtrack. This gives us a spatial scaleof a single path to use in the same way L/λ was usedfor the results presented in Fig. 3. To obtain a similarscale for the entire network, we perform a weighted av-erage of σ S over all possible paths in which the weightof each path is the current that runs through that par-ticular path when a steady, nonpulsatile flow in inputted into the far left node and outputted out of the far rightnode. We denote this current averaged value as ¯ σ . Wecan further obtain a temporal scale, ν S , for a single pathby considering the time required for a wavefront to tra-verse the path normalized by the σ S value of that path, ν S = ( R S dz τ ( z ) /λ ( z )) /σ S . This is equivalent to defining τ by the relation τ = ( zτ /λ ) / ( z/λ ) in the single vesselcase. The current averaged ¯ ν can then be calculated us-ing the same weighting scheme as the current averaged¯ σ . As an example, we consider the network shown inFig. 4. We numerically calculate the dissipation scal-ing function by first considering a boundary current ofthe form H ( t ) = √ ωt ), as this expression has twoFourier components of opposite frequency that each havea squared magnitude of 1 /
2. This boundary current isinputted into the far left node and simultaneously out-putted out of the far right node. All other nodes aretaken to have no inputted current, which is sufficient tocompletely characterize the current throughout the net-work and in turn calculate the time averaged dissipationtotalled over the entire network. This frequency depen-dent dissipation is then normalized by the dissipation inthe system when a constant, unit magnitude current isrun through the network with the same input and outputnodes. The normalized network dissipation, denoted as f net ( ω ), generalizes the concept of the dissipation scal-ing function to a network of vessels, as performing thisprocedure on a single vessel yields Eq. 7. To obtain afrequency independent value, we define h f net i for a par-ticular network as h f net i − = lim Ω →∞ Z Ω0 dω f − ( ω ) . (14)When applied to the single vessel case, this definitionreproduces the tanh( L/λ ) / ( L/λ ) limiting form of Eq. 7.With the dissipation properties determined, we nextlook to calculate the response time of the network. Weagain consider step function boundary conditions for thecurrent in the far left and right nodes, similar to thoseused to produce Eq. 8. Specifically, a current of the formΘ( t ) is inputted into the far left node and outputted outof the far right node. We neglect the case of pulsatileboundary conditions as it was shown in the single vesselcase to produce identical response times. The system issimulated by numerically evolving it through discretizedversions of Eqs. 5 and 6. We monitor the total current inall the central nodes and denote the residual current asthis total central current subtracted from the unit longtime limit. We then perform a linear fit to the logarithmof the magnitude of the residual current, as was done forthe data presented in Fig. 3, to obtain a measure of theresponse scaling function.These numerically derived values of the dissipationscaling function and response time are then compared inthe same manner as was done for the single vessel (Fig.3D). Specifically, we continue to examine the effects ofaltering the compliance of each vessel by fixing τ andvarying λ and α while holding their product constant.These parameters are changed throughout the networkin order to maintain the linear scaling of λ , τ , and α withcross sectional area. Additionally, the methods describedpreviously to numerically calculate the dissipation andresponse scaling functions can be applied to any net-work. Along with the simple symmetrically branching,area preserving network ( b = 1 / γ = 2), we also con-sider an asymmetric network in which one of the daugh-ter vessels has twice the area of the other each time abranching occurs ( b = 2 / γ = 2) as well as a network inwhich total area is distributed randomly between the twodaughters via a triangular distribution at each branchingnode ( b ∈ [0 . , . γ = 2). The symmetric, asym-metric, and random network construction methods areapplied to both the area preserving case ( γ = 2) and theMurray’s law case ( γ = 3). Finally, we also reconstructthe symmetric and asymmetric networks for γ = 2 and γ = 3 but with the inclusion of the looping vessels de-picted as dashed lines in Fig. 4. The area of each loopingvessel is taken to be the average of the two parent vesselsleading to the nodes the looping vessel connects. Therelation between the mean dissipation scaling functionand response scaling function for each of these networksis shown in Fig. 5A and in each case the inverse squareroot power law relation is seen to hold whenever the net-work size ¯ σ is sufficiently large. Additionally, we observethat the γ = 2 and γ = 3 networks naturally separateinto two distinct bands at large values of ¯ σ with larger γ leading to larger dissipation.This particular power law is seen to hold for each of thenetworks represented in 5A, but only when the value of ¯ σ is sufficiently large. To further explore this effect, we nextconsider a homogenous network constructed in the sameway as that depicted in Fig. 4 but with all vessels beingmade completely identical ( a =const. in Fig. 5B). Fromthis we also construct a randomized network in whichthe cross sectional area of each vessel in the homoge-nous network is chosen from a log-normal distribution( a =random in Fig. 5B). Loopy versions of these net-works are then constructed by again including the loop-ing vessels depicted in Fig. 4. Finally, the area of thelooping vessels are each increased by a factor of 100 tocreate a network in which the loops hold a significantamount of the volume (“dominant loops” in Fig. 5B).The dissipation-response scaling curves for these six net-works are shown in 5B, in which it is seen that the in-verse square root law yet again holds for each networkonce the value of ¯ σ passes a certain threshold. However,the networks with large loops in particular also show a“shoulder” regime in which dissipation is roughly con-stant while response time increases drastically. Despitethis, the trend of decreasing dissipation correlating withincreasing response time is still seen for all six networksconsidered in 5B when the networks are sufficiently large,and as such the tradeoff between dissipation and responsetime still exists. Response Scaling Function, 10 M e a n D i ss i p a t i o n S c a li n g F un c t i o n , f A f b =1/2, =2 b =2/3, =2 b =random, =2 b =1/2, =3 b =2/3, =3 b =random, =3loopy, b =1/2, =2loopy, b =2/3, =2loopy, b =1/2, =3loopy, b =2/3, =3 C u rr e n t a v e r a g e d Response Scaling Function, 10 M e a n D i ss i p a t i o n S c a li n g F un c t i o n , f B f a =const., no loops a =random, no loops a =const., with loops a =random, with loops a =const., dominant loops a =random, dominant loops C u rr e n t a v e r a g e d FIG. 5. Dissipation vs. response scaling functions for a vari-ety of different networks. A) Networks in which vessel branch-ing obeys Eq. 13 for γ = 2 or 3 show an overall relationshipbetween the dissipation and response scaling functions verysimilar to that seen in the single vessel (Fig. 3D), includingthe inverse square root power law relation for large networks.B) Networks with more homogeneous vessels also display theinverse square root power law relation. However, those withlarge looping vessels in particular appear to also contain aregime of nearly constant dissipation and rapidly increasingresponse time. III. DISCUSSION
We have shown that fluid flow through a single cylin-drical vessel comprised of compliant walls obeys Eqs. 3and 4 in the linear regime. By analyzing these equationsin Fourier space and transforming back to real space, wewere able to derive the dissipation and response scalingfunctions for such a vessel under symmetric boundaryconditions. These functions showed that so long as thecharacteristic time scale, τ , is held fixed, the dissipationand response scaling functions are related via an inversesquare root power law. Generalizing our theory to a net-work of such vessels and numerically calculating the dis-sipation and response scaling functions, we see that thispower law relation remains valid for a wide variety ofdifferent branching networks.The ubiquity of this scaling law can be understood byconsidering that in order for a network to be large enoughto fall into the regime in which this law applies, the pul-satile components of the current and pressure must sig-nificantly decay over a length scale shorter than the sizeof the network. When this occurs, it is the decay length,not the physical size of the network, that determines thetotal rate of energy dissipation from the pulsatile compo-nents. For the single vessel the dissipation scaling func-tion was shown to depend linearly on the characteristiclength scale λ in the long vessel limit, L/λ (cid:29)
1, where λ = 2 p ‘/c/r . These relations in turn cause the dissipa-tion to scale as c − / when r and ‘ are held fixed.The scaling of the response time of the system can beunderstood by considering the energy contained withinthe system. By multiplying Eq. 5 by αP and Eq. 6 by I then adding them together we can derive λ ∂∂z ( αP I ) + 2 I + τ ∂∂t (cid:16) I + α P (cid:17) = 0 . (15)Eq. 15 can be integrated over the length of a vessel toshow that the difference in power being delivered intoand out of each end of the vessel (the first term) can beaccounted for by the total rate of energy dissipation (thesecond term) and the rate of change in the energy storedin the kinetic energy of the fluid and elastic energy of thevessel walls (the third term). In the long vessel limit, theelastic term will dominate over the kinetic term and thetotal stored energy will scale as α = c/‘ . When a shift inboundary conditions occurs, this stored energy will nec-essarily need to shift to a new value. However, a bound-ary condition shift of a set magnitude can only deliverenergy into the system at a limited rate that scales inde-pendently of c as λα = 2 /r . Since the change in storedenergy will necessarily scale as α , the time needed toaccumulate a sufficient amount of energy will also scaleas α , or equivalently c if ‘ is held fixed. Thus, by com-paring this scaling relative to the vessel compliance, c , tothat of the dissipation rate, we obtain the inverse squareroot power law observed in the various networks studiedhere.The response scaling function has also been shown todictate the response time not just for the mean of the cur-rent but the pulsatile components as well. This has animportant consequence in that for any given network ofcompliant vessels there is a set timescale over which anysection of the network will be able to respond to changesin any other section. Specifically, the time for any gen-eral wavefront to traverse a path S is simply given by σ S ν S = R S dz τ ( z ) /λ ( z ), but if the path is too long rel-ative to the values of λ found along that path then theresulting wavefront will be substantially decayed and thetime required to see a significant change in the currentor pressure will be increased above the value σ S ν S . Forhighly hierarchical networks, this increase is well charac- terized by β (¯ σ ), as defined in Eq. 9, but as more and/orstronger loops are added β (¯ σ ) is seen to become less ac-curate, as exemplified by the creation of the“shoulder”regime seen in the large loop networks of Fig. 5. Forindividual paths within a network, the timescale σ S ν S sets a limit on how quickly mechanical information inthe form of fluid wavefronts can be transmitted alongthat path. In order for fluid pressure and flow to respondmore quickly, external information transmission is neces-sary. In biological contexts, this can be achieved throughelectrical signals in the nervous system, and the valuesof σ S and ν S along certain paths may dictate where inthe body such electrical information transmission is mostnecessary to maintain proper blood flow in response tosudden changes such as shifts in gravity from differentbody positions. The transmission of mechanical infor-mation in biological networks has some conceptual sim-ilarity to the propagation of the effects of a link failurein power grids [27, 28] - a concrete investigation of theanalogies could help understand more subtle aspects ofthe effects of topology on the response of vascular net-works to mechanical perturbations.The theory presented here captures the qualitative be-haviors of a network of compliant vessels transmittingpulses in flow and pressure, but it has limitations re-lated to the various approximations made to derive Eqs.5 and 6. Therefore, we do not expect a strict quanti-tative agreement in matters regarding the exact propa-gated waveform shape or the transmission and reflectioncoefficients at the nodes. Two distinct instances of suchlimitations come from assumptions made during the han-dling of −∇ u z . The axial term was neglected completelyas we assumed the wavelength of the waves traversing thesystem would be significantly larger than the vessel ra-dius. While this is typically valid in biological contexts,where λ can be tens of meters when the vessels are onlycentimeters across, any system in which this is not validwould require a fourth parameter beyond r , ‘ , and c or λ , τ , and α as well as a more complicated form of the func-tion k ( ωτ ). Additionally, the radial term was simplifiedunder the assumption that laminar flow with a quadraticvelocity profile is perpetually established within the ves-sel. This is not true in general, especially in a biologicalcontext, where the Womersley number can range fromorder ∼
10 in the major blood vessels, implying theflow oscillates too rapidly to maintain a quadratic ve-locity profile, to order ∼ − in the minor blood vessels,implying a quadratic velocity profile can be maintained[29, 30]. For nonquadratic profiles, the resistive term ofEq. 4 would have to be reevaluated based on the alter-nate profile used. However, for the human vasculaturespecifically, the quadratic approximation is only violatedin the largest vessels and is thus valid for the majority ofthe network.Possibly the most significant limitation of this theoryare the linear approximations. In the single vessel case,the nonlinear term of Eq. 2 is neglected. Previous com-putational studies have shown that nonlinear models of0fluid flow through compliant vessels are better able to re-produce experimentally measured pressure waveforms inthe major arteries [20, 31]. The model presented here canbe made to reflect these nonlinearities by reincorporat-ing the neglected term from Eq. 2. In the network case,another possible source of nonlinearity is in the pressureconnectivity law, Eq. 11. This can be made to moreaccurately incorporate Bernoulli’s principle by enforcingthat the sum of the pressure and kinetic energy densityis continuous across a network node rather than just thepressure. It will be interesting to explore the relationbetween dissipation and response time in these nonlinearregimes in future works.Despite these limitations, present in the majority ofvascular models that also linearize the flow equations [15–18, 25], the theory presented here provides many valuableinsights. In particular, the values of ¯ σ and ¯ ν can helpdetermine the behavior of complex networks for whichdata is available. As an example, we consider the humanvascular system. We can estimate the aorta to have adistensibility of 8 . × − mm Hg − and cross sectionalarea of 515 mm [32], while the blood in the aorta hasa dynamic viscosity of 3 . × − Pa · s and density of1050 kg · m − [33]. Using these values of calculate theresistance, inertia, and compliance per unit length thenallows us to derive the values λ ≈ . τ ≈ . λ and τ will scalelinearly with vessel area and that the vessel area itselfdecreases linearly from 515 mm to 50 µ m as the bloodflows from the aorta to a capillary that is 1 m away, wecan calculate that σ S ≈ . ν S ≈ . σ ≈ σ S and ¯ ν ≈ ν S .Specifically, by comparing this value of ¯ σ to those seenwithin Fig. 5, we can extrapolate that the human vascu-lar network is within the constant response regime ratherthan the inverse square root regime, but is very near thetransition point. Additionally, the value of ¯ ν dictatesthat the response time in this region is on the order of1 s, which is in agreement with timescales found in ratsfor changes in local oxygen concentration after a shift inheart rate occurs [34] (while rats are much smaller thanthe 1 m distance used to obtain the estimated value of¯ ν , ν S along any single path is only weakly dependenton distance due to it being normalized by σ S and dis-tance from any cell to its nearest capillary is relativelyconstant across organismal size, thus allowing the com-parison to be made). From these findings, we can makethe prediction that the vascular network may be adaptedto minimize dissipation while restricting itself to the re- gion of minimal possible response time. More detailedmeasurements of the value of ¯ σ and ¯ ν within the vas-culature of humans as well as other animals may revealthat existing at or near this transition point is a universaltrend. Moreover, they can shed light in situations wherethe body appears to have evolved to harness pulsatilityto perform specific functions, such as the movement ofcerebrospinal fluid [35].By examining the single vessel case, we can understandseveral observable effects of existing at this transitionpoint. In single vessels at the critical value of L/λ = π ,the wavefronts traveling through the vessel decay enoughthat the reflections do not cause sudden spikes in flowand dissipation, such as those seen in the blue and greencurves of Fig. 3A and B, but are not so decayed as toenter the regime in which flow and pressure expand diffu-sively. We can extrapolate these findings to networks topredict that when ¯ σ is too small wave reflections becomehighly significant and large spikes in flow and pressureshould be visible. This phenomenon is seen in the effectsof arterial stiffening. As blood vessels become stiffer,or equivalently their compliance lowers, the value of λ within each vessel must increase, in turn causing an in-crease in pulse wave velocity ( λ/τ ) and decrease in ¯ σ .This increase in wave velocity as well as increased ampli-tude of reflected waves can be directly observed in hyper-tensive patients with increased arterial stiffness [36, 37].Conversely, our theory predicts that should arterialstiffness be lowered not only would wave velocity and re-flected wave amplitude also decrease, but if the increasedcompliance causes the ¯ σ value of the network to increasepast the transition point then a marked increase wouldalso occur in the time required to establish an increasedblood flow in the capillary bed. Moving past this transi-tion point could be particularly detrimental for networkswith significant loops due to the sudden increase in re-sponse time created by the “shoulder” regime seen inFig. 5B. However, since arterial stiffness tends to increaserather than decrease as a consequence of age and/or dis-ease, this prediction is far more difficult to verify withexisting measurements. A more complete understandingof how these effects might apply to the human vascu-lar network could represent a step towards being ableto better diagnose disease and construct prosthetics andcardiac aids that work natively with the existing bloodvessels. ACKNOWLEDGMENTS
This research was supported by the NSF Award PHY-1554887 and the Simons Foundation through Award568888. [1] W. J. Lucas, A. Groover, R. Lichtenberger, K. Furuta, S.-R. Yadav, Y. Helariutta, X.-Q. He, H. Fukuda, J. Kang, S. M. Brady, et al. , “The plant vascular system: evolu- tion, development and functions f,” Journal of integrativeplant biology , vol. 55, no. 4, pp. 294–388, 2013.[2] R. Monahan-Earley, A. M. Dvorak, and W. C. Aird,“Evolutionary origins of the blood vascular system andendothelium,”
Journal of thrombosis and haemostasis ,vol. 11, pp. 46–66, 2013.[3] R. A. Gallego, A. J. Monticelli, and R. Romero, “Opti-mal capacitor placement in radial distribution networks,”
IEEE Transactions on Power Systems , vol. 16, no. 4,pp. 630–637, 2001.[4] G. A. Jim´enez-Est´evez, L. S. Vargas, and V. Marianov,“Determination of feeder areas for the design of large dis-tribution networks,”
IEEE Transactions on Power Deliv-ery , vol. 25, no. 3, pp. 1912–1922, 2010.[5] I. Ziari, G. Ledwich, A. Ghosh, and G. Platt, “Opti-mal distribution network reinforcement considering loadgrowth, line loss, and reliability,”
IEEE Transactions onPower Systems , vol. 28, no. 2, pp. 587–597, 2012.[6] C. D. Murray, “The physiological principle of minimumwork: I. the vascular system and the cost of blood vol-ume,”
Proceedings of the National Academy of Sciencesof the United States of America , vol. 12, no. 3, p. 207,1926.[7] T. F. Sherman, “On connecting large vessels to small.the meaning of murray’s law.,”
The Journal of generalphysiology , vol. 78, no. 4, pp. 431–453, 1981.[8] E. Katifori, G. J. Sz¨ollosi, and M. O. Magnasco, “Dam-age and fluctuations induce loops in optimal transportnetworks.,”
Physical Review Letters , vol. 104, no. 4,p. 048704, 2010.[9] H. Ronellenfitsch and E. Katifori, “Global Optimization, Local Adaptation , and the Role of Growth in Distri-bution Networks,”
Phys Rev Lett , vol. 117, p. 138301,2016.[10] H. Ronellenfitsch and E. Katifori, “Phenotypes of Vas-cular Flow Networks,”
Physical Review Letters , vol. 123,no. 24, p. 248101, 2019.[11] J. B. Kirkegaard and K. Sneppen, “Optimal TransportFlows for Distributed Production Networks,”
PhysicalReview Letters , vol. 124, no. 20, p. 208101, 2020.[12] E. Jones, M. Anliker, and I.-D. Chang, “Effects of vis-cosity and constraints on the dispersion and dissipationof waves in large blood vessels: Ii. comparison of analysiswith experiments,”
Biophysical journal , vol. 11, no. 12,pp. 1121–1134, 1971.[13] J. Hashimoto, “Central hemodynamics and target organdamage in hypertension,”
The Tohoku Journal of Exper-imental Medicine , vol. 233, no. 1, pp. 1–8, 2014.[14] R. Holenstein and D. N. Ku, “Reverse flow in the majorinfrarenal vessels–a capacitive phenomenon,”
Biorheol-ogy , vol. 25, no. 6, pp. 835–842, 1988.[15] S. Sherwin, V. Franke, J. Peir´o, and K. Parker, “One-dimensional modelling of a vascular network in space-time variables,”
Journal of engineering mathematics ,vol. 47, no. 3-4, pp. 217–250, 2003.[16] J. Alastruey, T. Passerini, L. Formaggia, and J. Peir´o,“Physical determining factors of the arterial pulse wave-form: theoretical analysis and calculation using the 1-d formulation,”
Journal of Engineering Mathematics ,vol. 77, no. 1, pp. 19–37, 2012.[17] J. Flores, J. Alastruey, and E. C. Poir´e, “A novel analyt-ical approach to pulsatile blood flow in the arterial net-work,”
Annals of biomedical engineering , vol. 44, no. 10,pp. 3047–3068, 2016. [18] B. Yigit and K. Pekkan, “Non-dimensional physics ofpulsatile cardiovascular networks and energy efficiency,”
Journal of The Royal Society Interface , vol. 13, no. 114,p. 20151019, 2016.[19] J. M. Harazny, C. Ott, U. Raff, J. Welzenbach, N. Kwella,G. Michelson, and R. E. Schmieder, “First experience inanalysing pulsatile retinal capillary flow and arteriolarstructural parameters measured noninvasively in hyper-tensive patients,”
Journal of hypertension , vol. 32, no. 11,pp. 2246–2252, 2014.[20] A. Bui, I. D. ˇSutalo, R. Manasseh, and K. Liffman, “Dy-namics of pulsatile flow in fractal models of vascularbranching networks,”
Medical & biological engineering &computing , vol. 47, no. 7, pp. 763–772, 2009.[21] Q. Pan, R. Wang, B. Reglin, G. Cai, J. Yan, A. R. Pries,and G. Ning, “A one-dimensional mathematical modelfor studying the pulsatile flow in microvascular net-works,”
Journal of biomechanical engineering , vol. 136,no. 1, 2014.[22] P. Perdikaris, L. Grinberg, and G. E. Karniadakis, “Aneffective fractal-tree closure model for simulating bloodflow in large arterial networks,”
Annals of biomedical en-gineering , vol. 43, no. 6, pp. 1432–1442, 2015.[23] F. K. B¨auerle, S. Karpitschka, and K. Alim, “Living sys-tem adapts harmonics of peristaltic wave for cost-efficientoptimization of pumping performance,”
Physical reviewletters , vol. 124, no. 9, p. 098102, 2020.[24] A. Barnard, W. Hunt, W. Timlake, and E. Varley, “Atheory of fluid flow in compliant tubes,”
Biophysical jour-nal , vol. 6, no. 6, pp. 717–724, 1966.[25] W. Cousins, P. A. Gremaud, and D. M. Tartakovsky,“A new physiological boundary condition for hemody-namics,”
SIAM Journal on Applied Mathematics , vol. 73,no. 3, pp. 1203–1223, 2013.[26] J. Masoliver and G. H. Weiss, “Telegrapher’s equationswith variable propagation speeds,”
Physical review E ,vol. 49, no. 5, p. 3852, 1994.[27] B. Sch¨afer, D. Witthaut, M. Timme, and V. Latora, “Dy-namically induced cascading failures in power grids,”
Na-ture Communications , vol. 9, no. 1975, 2018.[28] X. Zhang, D. Witthaut, and M. Timme, “TopologicalDeterminants of Perturbation Spreading in Networks,”
Physical Review Letters , vol. 125, no. 21, p. 218301, 2020.[29] J. R. Womersley, “Method for the calculation of velocity,rate of flow and viscous drag in arteries when the pressuregradient is known,”
The Journal of physiology , vol. 127,no. 3, p. 553, 1955.[30] G. B. West, J. H. Brown, and B. J. Enquist, “A generalmodel for the origin of allometric scaling laws in biology,”
Science , vol. 276, no. 5309, pp. 122–126, 1997.[31] J. Alastruey, A. W. Khir, K. S. Matthys, P. Segers, S. J.Sherwin, P. R. Verdonck, K. H. Parker, and J. Peir´o,“Pulse wave propagation in a model human arterial net-work: assessment of 1-d visco-elastic simulations againstin vitro measurements,”
Journal of biomechanics , vol. 44,no. 12, pp. 2250–2258, 2011.[32] I. Voges, M. Jerosch-Herold, J. Hedderich, E. Pardun,C. Hart, D. D. Gabbert, J. H. Hansen, C. Petko, H.-H.Kramer, and C. Rickers, “Normal values of aortic dimen-sions, distensibility, and pulse wave velocity in childrenand young adults: a cross-sectional study,”
Journal ofCardiovascular Magnetic Resonance , vol. 14, no. 1, p. 77,2012. [33] D. Kumar, R. Vinoth, and V. S. Raviraj Adhikari, “Non-newtonian and newtonian blood flow in human aorta: Atransient analysis,” Biomedical Research , 2017.[34] K. Masamoto, J. Kershaw, M. Ureshi, N. Takizawa,H. Kobayashi, K. Tanishita, and I. Kanno, “Apparentdiffusion time of oxygen from blood to tissue in ratcerebral cortex: implication for tissue oxygen dynamicsduring brain functions,”
Journal of Applied Physiology ,vol. 103, no. 4, pp. 1352–1358, 2007.[35] H. Mestre, J. Tithof, T. Du, W. Song, W. Peng, A. M.Sweeney, G. Olveda, J. H. Thomas, M. Nedergaard, andD. H. Kelley, “Flow of cerebrospinal fluid is driven by ar-terial pulsations and is reduced in hypertension,”
NatureCommunications , vol. 9, no. 1, 2018.[36] T. Weber, J. Auer, M. F. O’Rourke, E. Kvas, E. Lass-nig, R. Berent, and B. Eber, “Arterial stiffness, wavereflections, and the risk of coronary artery disease,”
Cir-culation , vol. 109, no. 2, pp. 184–189, 2004.[37] W. W. Nichols, S. J. Denardo, I. B. Wilkinson, C. M.McEniery, J. Cockcroft, and M. F. O’Rourke, “Effects ofarterial stiffness, pulse wave velocity, and wave reflectionson the central aortic pressure waveform,”
The journal ofclinical hypertension , vol. 10, no. 4, pp. 295–303, 2008.
Appendix A: Linearizing the Navier-Stokes Equation
Here we derive Eqs. 3 and 4 from Eqs. 1 and 2. Webegin by making two important assumptions. The firstis rotational symmetry, meaning that all dynamic vari-ables must be independent of angular position within thecylindrical vessel and the angular flow velocity must be 0.The second is that the fluid is incompressible, meaningthe flow velocity must obey Eq 1. From here we definethe volumetric flow rate I ( z, t ) = Z dA u z ( z, r, t ) = 2 π Z a dr ru z ( z, r, t ) , (A1)where R dA represents integration over the circular crosssection and a is the radius of the cross section at axialposition z . Integrating Eq. 1 over the cross sectionalarea thus yields0 = Z dA (cid:18) ∂u z ∂z + 1 r ∂∂r ( ru r ) (cid:19) = ∂I∂z + 2 π Z a dr ∂∂r ( ru r ) = ∂I∂z + 2 πau r ( z, a, t ) . (A2)Finally, we note that if the radial velocity at r = a isnonzero, then a itself must be changing at the same ratein order to accommodate the expanding or contractingfluid. This can be expressed as ∂A∂t = ∂∂t (cid:16) πa (cid:17) = 2 πa ∂a∂t = 2 πau r ( z, a, t ) . (A3)Inserting Eq. A3 into Eq. A2 produces ∂I∂z + ∂A∂t = ∂I∂z + ∂A∂P ∂P∂t = 0 . (A4)Thus, we can see that Eq. 3 can be produced solely viathe assumptions that the fluid is rotationally symmetricand incompressible.Next we turn to the Navier-Stokes equation itself,which for an incompressible fluid and no external forces,can be written as Eq. 2. The term ( ~u · ~ ∇ ) ~u is nonlinearand thus we will make the approximation that it is neg-ligible. This allows the axial component of Eq. 2 withthe assumed rotational symmetry to be extracted in theform ∂p∂z + ρ ∂u z ∂t − µ ∂ u z ∂z + 1 r ∂∂r (cid:18) r ∂u z ∂r (cid:19)! = 0 . (A5)We then average Eq. A5 over the cross section area toproduce1 A Z dA ∂p∂z + ρ ∂u z ∂t − µ ∂ u z ∂z − µ r ∂∂r (cid:18) r ∂u z ∂r (cid:19)! = ∂P∂z + ρA ∂I∂t − µA ∂ I∂z − πµA a ∂u z ∂r (cid:12)(cid:12)(cid:12)(cid:12) r = a = 0 , (A6)where the area averaged pressure, P , is defined as P ( z, t ) = 1 A Z dA p ( z, r, t ) . (A7)From here we assume that the fluid is Newtonian andflows in a laminar way. This means u z must follow theHagen-Poiseuille equation and be expressable as u z ( z, r, t ) = U ( z, t ) − r a ! . (A8)Inserting Eq. A8 into Eq. A1 then yields I ( z, t ) = 2 π Z a dr rU ( z, t ) − r a ! = 12 πa U ( z, t ) = 12 AU ( z, t ) . (A9)We can also differentiate Eq. A8 and combine it with Eq.A9 to produce − a ∂u z ∂r (cid:12)(cid:12)(cid:12)(cid:12) r = a = 2 U ( z, t ) = 4 A I. (A10)Inserting Eq. A10 into Eq. A6 yields3 ∂P∂z + ρA ∂I∂t + 8 πµA I − A π ∂ I∂z ! = 0 . (A11)We now make the final assumption that the term( A/ π )( ∂ I/∂z ) is negligibly small compared to I . Thisis valid when I has a roughly sinusoidal or exponentialform in space, thereby causing the second derivative tointroduce a factor that scales inversely with the squareof the wavelength or exponential constant ( ∼ λ − ), andwhen the wavelength or exponential length scale is sig-nificantly larger than the vessel radius, thereby causing A/λ to be a negligibly small parameter. By comparingEqs. A4 and A11 with this term excluded to Eqs. 5 and6 we see that they are equivalent assuming r = 8 πµA , (A12a) ‘ = ρA , (A12b) c = ∂A∂P , (A12c) λ = 2 r r ‘c , (A13a) τ = 2 ‘r (A13b) α = r c‘ . (A13c)Of important note is that while A is considered a dy-namic variable, r , ‘ , and c can all be made to be ap-proximately static when A ( z, t ) ≈ A + cP ( z, t ) and A (cid:29) cP ( z, t ). This approximation also allows us toexpress c as c ≈ A E ah , (A14)where E is the Young’s modulus of the vessel walls and h is the wall thickness [24]. Further assuming the a/h is approximately constant for varying A , Eq. A14 showsthat c scales linearly with A while Eqs. A12a and A12bshow that r and ‘ scale as A − and A − respectively.Inserting these scalings into Eq. A13 in turn shows that λ , τ , and α all scale as A . Appendix B: Fourier Space Solutions
We now denote ˜ I and ˜ P as the inverse Fourier trans-form of I and P with respect to time. This provides therelations ˜ I ( z, ω ) = Z dt π e − iωt I ( z, t ) , (B1a) I ( z, t ) = Z dω e iωt ˜ I ( z, ω ) , (B1b)˜ P ( z, ω ) = Z dt π e − iωt P ( z, t ) , (B2a) P ( z, t ) = Z dω e iωt ˜ P ( z, ω ) . (B2b)Substituting Eqs. B1b and B2b into Eqs. 5 and 6 thenperforming the inverse Fourier transform operation oneach equation thus yields λ ∂ ˜ I∂z + iωτ α ˜ P = 0 , (B3) λ ∂∂z (cid:16) α ˜ P (cid:17) + iωτ ˜ I + 2 ˜ I = 0 . (B4)To obtain solutions for ˜ I and ˜ P , we first solve Eq. B3for α ˜ P and substitute that into Eq. B4 to produce − λ iωτ ∂ ˜ I∂z + (2 + iωτ ) ˜ I = 0 . (B5)For a vessel of length L , Eq. B5 has the general solution˜ I ( z, ω ) = ˜ I F ( ω ) e − zλ k ( ωτ ) − ˜ I B ( ω ) e − L − zλ k ( ωτ ) , (B6)where k ( ωτ ) = p iωτ (2 + iωτ ) , (B7)and ˜ I F ( ω ) = ˜ I (0 , ω ) e Lλ k ( ωτ ) − ˜ I ( L, ω ) e Lλ k ( ωτ ) − e − Lλ k ( ωτ ) , (B8a)˜ I B ( ω ) = − ˜ I ( L, ω ) e Lλ k ( ωτ ) − ˜ I (0 , ω ) e Lλ k ( ωτ ) − e − Lλ k ( ωτ ) , (B8b)4are the forward and backward propagating current waveamplitudes. Substituting Eq. B6 back into Eq. B3 andsolving for ˜ P then yields˜ P ( z, ω ) = k ( ωτ ) iωτ α (cid:16) ˜ I F ( ω ) e − zλ k ( ωτ ) + ˜ I B ( ω ) e − L − zλ k ( ωτ ) (cid:17) . (B9)Eqs. B6 and B9 can be fully solved once sufficientboundary conditions are given. In the case where the cur-rent boundary conditions are known, ˜ I F ( ω ) and ˜ I B ( ω )can be calculated directly from Eq. B8, thus allowing˜ I ( z, ω ) and ˜ P ( z, ω ) to be calculated from Eqs. B6 andB9. Alternatively, when the pressure boundary condi-tions are known, a similar process yields the relations˜ I ( z, ω ) = iωτ αk ( ωτ ) (cid:16) ˜ P F ( ω ) e − zλ k ( ωτ ) − ˜ P B ( ω ) e − L − zλ k ( ωτ ) (cid:17) , (B10)˜ P ( z, ω ) = ˜ P F ( ω ) e − zλ k ( ωτ ) + ˜ P B ( ω ) e − L − zλ k ( ωτ ) , (B11)where ˜ P F ( ω ) = ˜ P (0 , ω ) e Lλ k ( ωτ ) − ˜ P ( L, ω ) e Lλ k ( ωτ ) − e − Lλ k ( ωτ ) , (B12a)˜ P B ( ω ) = ˜ P ( L, ω ) e Lλ k ( ωτ ) − ˜ P (0 , ω ) e Lλ k ( ωτ ) − e − Lλ k ( ωτ ) . (B12b)Of important note is that Eqs. B6-B12 are definedassuming the forward and backward waves move and de-cay in the forward and backward z direction respectively.This forces the choice of which root to use for evaluating k ( ωτ ) in Eq. B7 to be the principle root for all real ω .Eqs. B6 and B9 can be equivalently expressed in a waythat is even in k ( ωτ ) and thus independent of which rootis taken. These take the forms˜ I ( z, ω ) =˜ I (0 , ω ) sinh (cid:16) L − zλ k ( ωτ ) (cid:17) + ˜ I ( L, ω ) sinh (cid:0) zλ k ( ωτ ) (cid:1) sinh (cid:16) Lλ k ( ωτ ) (cid:17) , (B13)˜ P ( z, ω ) = k ( ωτ ) iωτ α · ˜ I (0 , ω ) cosh (cid:16) L − zλ k ( ωτ ) (cid:17) − ˜ I ( L, ω ) cosh (cid:0) zλ k ( ωτ ) (cid:1) sinh (cid:16) Lλ k ( ωτ ) (cid:17) . (B14) When the current boundary conditions are symmetric( ˜ I (0 , ω ) = ˜ I ( L, ω )), Eq. B13 also shows that the ampli-tude of current oscillations at position z and frequency ω obeys Eq. 10.In the limit ω →
0, Eq. B14 diverges unless ˜ I (0 ,
0) =˜ I ( L, ω → I ( z,
0) to be invariant with respect to changesin z . However, in cases where the limit of ˜ I ( z, ω ) as ω → I ( z, t ) is a pul-satile function with discrete frequencies, Eq. B14 be-comes equally ill-defined. This is a consequence of thefact that pressure is a gauge variable and globally chang-ing the pressure across the whole vessel and at all timesonly affects ˜ P ( z,
0) has no impact on the current. Thus,while Eq. B14 retains a linear dependence on z fromthe hyperbolic cosine terms, the z -independent constantterm must be determined by choice of gauge. This issuecan be circumvented by defining the pressure boundaryconditions instead, as the gauge choice would be includedinto the boundary conditions. This allows for the currentand pressure Fourier transforms to be expressed as˜ I ( z, ω ) = iωτ αk ( ωτ ) · ˜ P (0 , ω ) cosh (cid:16) L − zλ k ( ωτ ) (cid:17) − ˜ P ( L, ω ) cosh (cid:0) zλ k ( ωτ ) (cid:1) sinh (cid:16) Lλ k ( ωτ ) (cid:17) , (B15)˜ P ( z, ω ) =˜ P (0 , ω ) sinh (cid:16) L − zλ k ( ωτ ) (cid:17) + ˜ P ( L, ω ) sinh (cid:0) zλ k ( ωτ ) (cid:1) sinh (cid:16) Lλ k ( ωτ ) (cid:17) . (B16)Unlike Eq. B14, Eq. B15 is well defined in the limit ω → P (0 , ω ) and ˜ P ( L, ω ). Appendix C: Single Vessel Dissipation
We now consider the time averaged dissipation ratewithin a vessel with symmetric current boundary condi-tions in which I (0 , t ) = I ( L, t ). We further fix I (0 , t ) and I ( L, t ) to be pulsatile with a period of T and frequency ω = 2 π/T . This allows I ( z, t ) to expressed as the Fouriersum I ( z, t ) = ∞ X n = −∞ ˜ I ( n ) ( z ) e inωt , (C1)where ˜ I ( n ) ( z ) = 1 T Z T dt I ( z, t ) e − inωt . (C2)5We can then define the time averaged dissipation rate, h D i as h D i = 1 T Z T dt Z L dz rI ( z, t )= 1 T Z T dt Z L dz r X n,n ˜ I ( n ) ( z ) e inωt ˜ I ( n ) ( z ) e in ωt = r Z L dz X n ˜ I ( n ) ( z ) ˜ I ( − n ) ( z ) (C3)The same analysis used in Appendix B can be ap-plied to the case of discrete frequencies considered hereto show that Eqs. B13-B16 still hold when ωτ , ˜ I ( z, ω ),and ˜ P ( z, ω ) are replaced by nωτ , ˜ I ( n ) ( z ), and ˜ P ( n ) ( z )respectively. Additionally, since I ( z, t ) is a real variable,˜ I ( n ) ∗ ( z ) = ˜ I ( − n ) ( z ). Thus, in the case of the aforemen-tioned symmetric boundary conditions, Eq. B13 can becombined with Eq. C3 to produce h D i = r Z L dz X n (cid:12)(cid:12)(cid:12) ˜ I ( n ) ( z ) (cid:12)(cid:12)(cid:12) = rL X n (cid:12)(cid:12)(cid:12) ˜ I ( n ) (0) (cid:12)(cid:12)(cid:12) · Z L dzL (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) sinh (cid:16) L − zλ k ( nωτ ) (cid:17) + sinh (cid:0) zλ k ( nωτ ) (cid:1) sinh (cid:16) Lλ k ( nωτ ) (cid:17) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (C4)The integral in the last line of Eq. C4 is the dissipa-tion scaling function and can be solved by separating thehyperbolic sine terms into exponentials and evaluatingthe magnitude. The resulting summation of exponentialterms can then be easily integrated and recombined toform Eq. 7. Appendix D: Network Laplacian
For a network of vessels that obey the connectivitylaws given by Eqs. 11 and 12, we can define the pressureat node µ as P µ ( t ) = P µν (0 , t ) = P νµ ( L µν , t ) for ν ∈ N µ .This notation allows Eq. 12 to be Fourier transformedand Eq. B15 to be substituted in to produce˜ H µ ( ω ) = X ν ∈N µ ˜ I µν (0 , ω )= X ν ∈N µ iωτ µν α µν k (cid:0) ωτ µν (cid:1) ˜ P µ ( ω ) cosh (cid:16) L µν λ µν k (cid:0) ωτ µν (cid:1)(cid:17) − ˜ P ν ( ω )sinh (cid:16) L µν λ µν k (cid:0) ωτ µν (cid:1)(cid:17) = X ν ∈N µ L µν ( ω ) ˜ P ν ( ω ) , (D1) where L µν ( ω ) = δ µ,ν X ξ ∈N µ iωτ µξ α µξ cosh (cid:16) L µξ λ µξ k (cid:0) ωτ µξ (cid:1)(cid:17) k (cid:0) ωτ µξ (cid:1) sinh (cid:16) L µξ λ µξ k (cid:0) ωτ µξ (cid:1)(cid:17) − iωτ µν α µν k (cid:0) ωτ µν (cid:1) sinh (cid:16) L µν λ µν k (cid:0) ωτ µν (cid:1)(cid:17) . (D2)For nonzero ω , L µν defines an invertible matrix, mean-ing that for any set of input ˜ H µ there exists a unique setof potentials ˜ P µ that solves Eq. D1 and vice versa. Thus,there is no restriction on the oscillatory components of H µ ( t ). However, in the ω → ω → L µν ( ω )= δ µ,ν X ξ ∈N µ α µξ λ µξ L µξ · lim ω → iωτ µξ (cid:16) k (cid:0) ωτ µξ (cid:1)(cid:17) − α µν λ µν L µν · lim ω → iωτ µν (cid:16) k (cid:0) ωτ µν (cid:1)(cid:17) = δ µ,ν X ξ ∈N µ α µξ λ µξ L µξ − α µν λ µν L µν = δ µ,ν X ξ ∈N µ r µξ L µξ − r µν L µν . (D3) L µν as defined in Eq. D3 gives a noninvertible matrixwhich describes a transformation that takes the vectorspace R N , where N is the number of nodes in the net-work, to the subspace S ⊂ R N defined such that S is thespace of all vectors v ∈ R N whose components sum to0. Thus, the ω → H µ (0) sum to 0. This is equivalentto the restriction that the constant components of H µ ( t )must all sum to 0 so that on average the current goinginto the network equals the current coming out. upporting information for “Tradeoffs between energy efficiency and mechanicalresponse in fluid flow networks” GENERAL SOLUTIONS
From the main text we have that the volumetric current, I ( z, t ), and area-averaged pressure, P ( z, t ), in a singlevessel must obey λ ∂I∂z + ατ ∂P∂t = 0 , (S1) αλ ∂P∂z + τ ∂I∂t + 2 I = 0 . (S2)These can be rewritten in a more simplified version by defining the volume function, V ( z, t ), such that I ( z, t ) = τ ∂V∂t , (S3a) αP ( z, t ) = − λ ∂V∂z . (S3b)Under these definitions, Eq. S1 is automatically satisfied and Eq. S2 takes the form − λ ∂ V∂z + τ ∂ V∂t + 2 τ ∂V∂t = 0 . (S4)Before proceeding, we will simplify the notation by enforcing units in which λ = τ = α = 1. Proper units can bereinstated at any point by making the substitutions z → z/λ , t → t/τ , and P → αP . These simplified units allowEq. S4 to be expressed as − ∂ V∂z + ∂ V∂t + 2 ∂V∂t = 0 . (S5)Eq. S5 has a set of very simple solutions that come from assuming that V ( z, t ) can be separated into a sum ofterms that are each either spatially or temporally dependent, but not both: V ( z, t ) = 1 , V ( z, t ) = z, V ( z, t ) = e − t , V ( z, t ) = t + z . (S6)Each of these solutions represents a specific state of the flow and pressure within the vessel. The first solution simplyreflects the gauge symmetry of V itself by allowing any constant to be added to V without changing the physics ofthe system. The second solution represents the gauge symmetry of P by allowing any constant to be added to P without affecting I . The third solution represents the loss of flow energy to friction and shows that the current willexponentially decay away if there is no pressure gradient to drive it. The fourth solution simply dictates Ohm’s lawin the case of steady, time-independent flow.To derive the more complicated set of solutions to Eq. S5, we first rewrite V ( z, t ) in the form V ( z, t ) = W ( z, t )exp( − t ). This allows Eq. S5 to be reexpressed as " − ∂ ∂z + ∂ ∂t + 2 ∂∂t W ( z, t ) e − t = e − t − ∂ W∂z + ∂ W∂t − W ! = 0 . (S7) a r X i v : . [ phy s i c s . b i o - ph ] F e b Assuming W can be separated into a product of one term dependent on z and another dependent on t yields thesolution V ( z, t ) = e − t e at e bz ∀ a − b = 1 . (S8)Of note is that the case of a = 1 and b = 0 is simply the first solution in Eq. S6 while that of a = − b = 0 isthe third. Both a and b can take on complex values with the real parts determining whether the current and pressurewith exponentially grow or decay in space and time and the imaginary parts determining how quickly they oscillatein space and time.Another set of solutions can be obtained by defining the quantities q ( z, t ) = p t − z , (S9) s ( z, t ) = r t − zt + z , (S10)and their derivatives ∂q∂z = − z √ t − z = s − s − , (S11a) ∂q∂t = t √ t − z = s + s − , (S11b) ∂s∂z = − √ t − z − s t − z ( t + z ) = − sq s + s − , (S12a) ∂s∂t = 12 √ t − z − s t − z ( t + z ) = − sq s − s − . (S12b)Neglecting the global factor of exp( − t ) in Eq. S7 and allowing q and s to be the independent variables then produces " − ∂ ∂z + ∂ ∂t − W ( z, t ) = " − (cid:18) ∂q∂z ∂∂q + ∂s∂z ∂∂s (cid:19) + (cid:18) ∂q∂t ∂∂q + ∂s∂t ∂∂s (cid:19) − W ( q, s )= (cid:18) ∂q∂t (cid:19) − (cid:18) ∂q∂z (cid:19) ! ∂ ∂q + 2 (cid:18) ∂q∂t ∂s∂t − ∂q∂z ∂s∂z (cid:19) ∂ ∂q∂s + (cid:18) ∂s∂t (cid:19) − (cid:18) ∂s∂z (cid:19) ! ∂ ∂s + ∂q∂t ∂∂q (cid:18) ∂q∂t (cid:19) + ∂s∂t ∂∂s (cid:18) ∂q∂t (cid:19) − ∂q∂z ∂∂q (cid:18) ∂q∂z (cid:19) − ∂s∂z ∂∂s (cid:18) ∂q∂z (cid:19)! ∂∂q + ∂q∂t ∂∂q (cid:18) ∂s∂t (cid:19) + ∂s∂t ∂∂s (cid:18) ∂s∂t (cid:19) − ∂q∂z ∂∂q (cid:18) ∂s∂z (cid:19) − ∂s∂z ∂∂s (cid:18) ∂s∂z (cid:19)! ∂∂s − W ( q, s )= 1 q " q ∂ ∂q − s ∂ ∂s + q ∂∂q − s ∂∂s − q W ( q, s ) = 0 . (S13)We can again use separation of variables and assume W is the product of a q dependent factor and an s dependentfactor. The s dependent factor must be an eigenfunction of the s -derivative terms in Eq. S13, s ( ∂ /∂s ) + s ( ∂/∂s ),which means it must be of the form s n with n being the eigenvalue. Letting W ( q, s ) = s n X ( q ) then transforms Eq.S13 into 1 q " q ∂ ∂q − s ∂ ∂s + q ∂∂q − s ∂∂s − q s n X ( q ) = s n q q ∂ X∂q + q ∂X∂q − (cid:16) q + n (cid:17) X ! = 0 . (S14)The differential equation for X in Eq. S14 is exactly the differential equation for the n th modified Bessel function.Thus, the volume function can take the form V ( z, t ) = e − t W ( z, t ) = e − t s n ( z, t ) I n (cid:0) q ( z, t ) (cid:1) ≡ V n ( z, t ) , (S15)where I n ( q ) is the n th modified Bessel function of the first kind. Replacing I n ( q ) with K n ( q ), the n th modified Besselfunction of the second kind, is also a valid solution, but here we will work exclusively with solutions involving I n ( q ).In the case where q and s take on imaginary values ( | z | > | t | ), the modified Bessel functions become regular Besselfunctions and V ( z, t ) retains a real value.The solutions presented in Eq. S15 are particularly useful when considering the relation (cid:20) ± ∂∂z + ∂∂t + 1 (cid:21) e − t s n ( z, t ) I n (cid:0) q ( z, t ) (cid:1) = e − t (cid:20) ± ∂∂z + ∂∂t (cid:21) s n ( z, t ) I n (cid:0) q ( z, t ) (cid:1) = e − t "(cid:18) ± ∂q∂z + ∂q∂t (cid:19) ∂∂q + (cid:18) ± ∂s∂z + ∂s∂t (cid:19) ∂∂s s n I n ( q ) = e − t s n ± (cid:18) ∂I n ( q ) ∂q ∓ nI n ( q ) q (cid:19) = e − t s n ± ( z, t ) I n ± (cid:0) q ( z, t ) (cid:1) , (S16)where the recursion properties of the modified Bessel functions of the first kind are used to form the final equality.Thus, the differential operator presented in the first form on Eq. S16 can be taken to be the raising (+) or lowering( − ) operator for V n ( z, t ). Denoting these operators as N ± = ± ∂/∂z + ∂/∂t + 1 also allows for the relations ∂∂z = 12 (cid:0) N + − N − (cid:1) , (S17a) ∂∂t = 12 (cid:0) N + + N − (cid:1) − , (S17b) − ∂ ∂z + ∂ ∂t + 2 ∂∂t = N + N − − . (S17c) SPECIFIC SINGLE VESSEL SOLUTIONS
We can now consider the case in which V ( z = 0 , t ) = Θ( t ), where Θ( t ) is the Heaviside step function. Thiscorresponds to the current at position z = 0 taking the form of a δ -pulse, I ( z = 0 , t ) = ∂V /∂t = δ ( t ). This particularform of V ( z, t ) can be expressed in terms of V n ( z, t ) by noting that0 = " − ∂ ∂z + ∂ ∂t + 2 ∂∂t Θ ( t − z ) V n ( z, t ) = Θ ( t − z ) (cid:0) N + N − − (cid:1) V n ( z, t ) + 2 δ ( t − z ) N + V n ( z, t )= 2 δ ( t − z ) V n +1 ( z, t ) . (S18)In the limit z → t imposed by the δ -function, q ( z, t ) goes to 0 and s n I n ( q ) for n ≥ sq/ n / Γ( n + 1), where Γ( x ) is the gamma function. Since sq = t − z , Eq. S18 must be satisfied for any n > − z = 0, q = | t | and s = 1. This allows us to use the generating function for modifiedBessel functions of the first kind, e t = e q ( z,t )2 ( s ( z,t )+ s − ( z,t ) ) = ∞ X n = −∞ s n ( z, t ) I n (cid:0) q ( z, t ) (cid:1) , (S19)to express Θ( t ) asΘ ( t ) = Θ ( t ) e − t e t = Θ ( t ) e − t I ( t ) + 2 ∞ X n =1 I n ( t ) = Θ ( t − z ) V ( z, t ) + 2 ∞ X n =1 V n ( z, t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) z =0 . (S20)Since each term in Eq. S20 satisfies the form and n > − z = 0 represents a solution to the boundary condition problem V ( z = 0 , t ) = Θ( t ).Another solution is possible, but is represented by a leftward moving wave, which is nonphysical in the case of a vesselthat starts at z = 0 and expands to the right. This in turn produces current and pressure functions of the form I ( z, t ) = ∂V∂t = δ ( t − z ) V ( z, t ) + 2 ∞ X n =1 V n ( z, t ) + Θ ( t − z ) (cid:18) (cid:0) N + + N − (cid:1) − (cid:19) V ( z, t ) + 2 ∞ X n =1 V n ( z, t ) = δ ( t − z ) lim z → t V ( z, t ) + 2 ∞ X n =1 V n ( z, t ) + Θ ( t − z ) (cid:0) V ( z, t ) + V − ( z, t ) (cid:1) − V ( z, t ) + ∞ X n =1 (cid:0) V n +1 ( z, t ) + V n − ( z, t ) − V n ( z, t ) (cid:1) = δ ( t − z ) e − t + Θ ( t − z ) e − t s − ( z, t ) − s ( z, t )2 I (cid:0) q ( z, t ) (cid:1) = δ ( t − z ) e − t + Θ ( t − z ) ze − t I (cid:0) q ( z, t ) (cid:1) q ( z, t ) , (S21) P ( z, t ) = − ∂V∂z = δ ( t − z ) V ( z, t ) + 2 ∞ X n =1 V n ( z, t ) − Θ ( t − z ) 12 (cid:0) N + − N − (cid:1) V ( z, t ) + 2 ∞ X n =1 V n ( z, t ) = δ ( t − z ) lim z → t V ( z, t ) + 2 ∞ X n =1 V n ( z, t ) − Θ ( t − z ) (cid:0) V ( z, t ) − V − ( z, t ) (cid:1) + ∞ X n =1 (cid:0) V n +1 ( z, t ) − V n − ( z, t ) (cid:1) = δ ( t − z ) e − t + Θ ( t − z ) e − t I (cid:0) q ( z, t ) (cid:1) + s − ( z, t ) + s ( z, t )2 I (cid:0) q ( z, t ) (cid:1)! = δ ( t − z ) e − t + Θ ( t − z ) e − t I (cid:0) q ( z, t ) (cid:1) + t I (cid:0) q ( z, t ) (cid:1) q ( z, t ) ! . (S22)Eqs. S21 and S22 show how the δ -pulse of current travels through and decays along the length of the vessel as wellas how the current and pressure relax back to 0 after the pulse passes.With the current and pressure from a δ -pulse of current solved for, the solution for the current given any arbitraryboundary condition I ( z = 0 , t ) = H ( t ) can be expressed as I ( z, t ) = Z ∞−∞ dt H (cid:0) t − t (cid:1) δ (cid:0) t − z (cid:1) e − t + Θ (cid:0) t − z (cid:1) ze − t I (cid:16) q (cid:0) z, t (cid:1)(cid:17) q ( z, t ) = H ( t − z ) e − z + z Z ∞ z dt H (cid:0) t − t (cid:1) e − t I (cid:16) q (cid:0) z, t (cid:1)(cid:17) q ( z, t ) . (S23)Eq. S23 shows that the current at an arbitrary time can be separated into one term that represents the decayed andtime delayed input current value (the first term) and one term that represents the residual current from all previousinput current values (the second term).We will now consider the specific case of a Heaviside-pulsatile boundary condition, H ( t ) = Θ( t )exp( iωt ). Addi-tionally, by rewriting the factor of z/q ( z, t ) as ( s − ( z, t ) − s ( z, t )) /
2, we can take advantage of Eq. S19 and theresidue theorem to express ( s − − s ) I ( q ) = s − I − ( q ) − sI ( q ) as a contour integral over the variable x by makingthe substitution s → s/x . This allows Eq. S23 to take the form I ( z, t ) = Θ ( t − z ) e iω ( t − z ) e − z + Z ∞ z dt Θ (cid:0) t − t (cid:1) e iω ( t − t ) e − t s − (cid:0) z, t (cid:1) − s (cid:0) z, t (cid:1) I (cid:16) q (cid:0) z, t (cid:1)(cid:17) = Θ ( t − z ) e iωt e − z (1+ iω ) + Z tz dt e − iωt e − t πi I dx (cid:16) x − − (cid:17) e q ( z,t ) s ( z,t ) x + xs ( z,t ) = Θ ( t − z ) e iωt e − z (1+ iω ) + 14 πi Z tz dt I dx e − t (1+ iω ) (cid:16) x − − (cid:17) e x ( t − z ) + x ( t + z ) ! = Θ ( t − z ) e iωt e − z (1+ iω ) + 14 πi I dx x − − x + x − − iω e z ( x − x − ) (cid:16) e t ( x + x − − iω ) − e z ( x + x − − iω ) (cid:17)! = Θ ( t − z ) e iωt e − z (1+ iω ) + 12 πi I dx x − x (1 − x ) − iωx (cid:16) e x ( t + z )+ x ( t − z ) − t (1+ iω ) − e z ( x − − iω ) (cid:17)! . (S24)The contour of integration around x can be any path that fully encloses the essential singularity at x = 0. To aidin later calculations, we will choose the path to be the unit circle. This path also contains a pole at x = ‘ ( ω ) =1 + iω − k ( ω ), where k ( ω ) = p iω (2 + iω ) with the principle root being taken as in the main text, as this value of x satisfies (1 − x ) − iωx = ( x − ‘ ( ω ))( x − ‘ − ( ω )) = 0. Given this pole as well as the necessary pole at x = 0, theresidue theorem can be easily applied to the second term within the integral as exp( z ( x − − iω )) has no poles atfinite x . Evaluating this term separately yields I ( z, t ) = Θ ( t − z ) e iωt e − z (1+ iω ) − lim x → − x (1 − x ) − iωx e z ( x − − iω ) ! − lim x → ‘ ( ω ) x − x x − ‘ − ( ω ) e z ( x − − iω ) ! + e − t (1+ iω ) πi I dx x − x (1 − x ) − iωx e t ( x + x − ) e z ( x − x − ) ! = Θ ( t − z ) e iωt e − zk ( ω ) + e − t πi I dx x − x (1 − x ) − iωx e t ( x + x − ) e z ( x − x − ) ! . (S25)Eq. S25 can be further simplified by reexpressing exp( t ( x + x − ) / z ( x − x − ) /
2) yet again asexp( q ( z, t )( s ( z, t ) /x + x/s ( z, t )) /
2) and using Eq. S19 to further expand it into powers of x . Similarly, we canutilize the generating function 1 − x ( x − ‘ ) ( x − ‘ − ) = 1 + ∞ X n =1 (cid:0) ‘ n + ‘ − n (cid:1) x n (S26)to expand the factor of (1 − x ) / ((1 − x ) − iωx ). The poles of this factor exist at x = 1 + iω ± k ( ω ), where k ( ω ) = p iω (2 + iω ) with the principle root being taken as in the main text. Since (1 + iω + k ( ω ))(1 + iω − k ( ω )) = 1,the expansion presented in Eq. S26 is valid for ‘ ( ω ) = 1 + iω − k ( ω ). These expansions allow the residue theorem to beapplied around the poles at x = 0, which will cause all terms with powers of x other than − iωt )exp( − zk ( ω )) was introduced by including the artificially created pole at x = ‘ ( ω ), it willnecessarily be canceled by the same pole in the remaining integral. This term can thus be neglected by also neglectingthe same pole in the remaining integral, thus simplifying Eq. S25 to I ( z, t ) = Θ ( t − z ) e − t πi I dx x − x (cid:0) x − ‘ ( ω ) (cid:1) (cid:0) x − ‘ − ( ω ) (cid:1) e q ( z,t )2 (cid:16) s ( z,t ) x + xs ( z,t ) (cid:17) = Θ ( t − z ) e − t πi I dx x ∞ X m =1 (cid:0) ‘ m ( ω ) + ‘ − m ( ω ) (cid:1) x m ∞ X n = −∞ (cid:18) s ( z, t ) x (cid:19) n I n (cid:0) q ( z, t ) (cid:1) = Θ ( t − z ) e − t I (cid:0) q ( z, t ) (cid:1) + ∞ X n =1 (cid:0) ‘ n ( ω ) + ‘ − n ( ω ) (cid:1) s n ( z, t ) I n (cid:0) q ( z, t ) (cid:1) = Θ ( t − z ) V ( z, t ) + ∞ X n =1 (cid:0) ‘ n ( ω ) + ‘ − n ( ω ) (cid:1) V n ( z, t ) . (S27) FINITE VESSEL STEP-FUNCTION PROPAGATION
We can apply this understanding of a step function boundary condition to a finite sized single vessel of length L aswell. To do so we first consider a more general case of a single vessel with arbitrary and symmetric current boundaryconditions denoted as H ( t ) = I ( z = 0 , t ) = I ( z = L, t ). This has the benefit of automatically ensuring that the0-frequency components of the currents input into both sides of the vessel sum to 0, thus satisfying the conditionimposed by Eq. D3 of the main text. Since a single vessel is also a two node system, Eqs. B15 and D1 (in the units λ = τ = α = 1) can be used to express the Fourier current in the matrix form˜ I ( z, ω ) = (cid:20) iω cosh ( ( L − z ) k ( ω ) ) k ( ω ) sinh ( Lk ( ω ) ) − iω cosh ( zk ( ω ) ) k ( ω ) sinh ( Lk ( ω ) ) (cid:21) " ˜ P ( z = 0 , ω )˜ P ( z = L, ω ) = iωk ( ω ) (cid:20) cosh ( ( L − z ) k ( ω ) ) sinh ( Lk ( ω ) ) − cosh ( zk ( ω ) ) sinh ( Lk ( ω ) ) (cid:21) iω cosh ( Lk ( ω ) ) k ( ω ) sinh ( Lk ( ω ) ) − iωk ( ω ) sinh ( Lk ( ω ) ) − iωk ( ω ) sinh ( Lk ( ω ) ) iω cosh ( Lk ( ω ) ) k ( ω ) sinh ( Lk ( ω ) ) − " ˜ H ( ω ) − ˜ H ( ω ) = ˜ H ( ω ) (cid:16) cosh (cid:0) ( L − z ) k ( ω ) (cid:1) + cosh (cid:0) zk ( ω ) (cid:1)(cid:17) (cid:16) cosh (cid:0) Lk ( ω ) (cid:1) − (cid:17) sinh (cid:0) Lk ( ω ) (cid:1) . (S28)From here, Eq. S28 can be simplified by noting that Lk ( ω ) = ( L − z ) k ( ω ) + zk ( ω ) and using the minus mode ofthe relation cosh ( a + b ) ± a + b ) = (cid:18) cosh ( a ) + cosh ( b )sinh ( a ) + sinh ( b ) (cid:19) ± . (S29)The hyperbolic functions can then be rewritten as exponentials and Taylor expanded into powers of exp( − Lk ( ω ))to produce˜ I ( z, ω ) = ˜ H ( ω ) sinh (cid:0) ( L − z ) k ( ω ) (cid:1) + sinh (cid:0) zk ( ω ) (cid:1) sinh (cid:0) Lk ( ω ) (cid:1) = ˜ H ( ω ) e − zk ( ω ) − e − (2 L − z ) k ( ω ) + e − ( L − z ) k ( ω ) − e − ( L + z ) k ( ω ) − e − Lk ( ω ) = ˜ H ( ω ) (cid:16) e − zk ( ω ) − e − (2 L − z ) k ( ω ) + e − ( L − z ) k ( ω ) − e − ( L + z ) k ( ω ) (cid:17) ∞ X n =0 e − nLk ( ω ) . (S30)Each term in Eq. S30 has a clear interpretation as a normal or reflected wave. The first term, exp( − zk ( ω )), is apositive current wave that has traveled a distance z and thus represents the primary wave originating at the z = 0 endof the vessel. The second term is has an additional factor of − L − z , which represents anegative current wave that originated at z = 0, traveled the entire length of the vessel, was reflected at the z = L end,and traveled the distance L − z back to the position z . Similarly, the third term is evaluated at L − z and representsthe positive current wave that originated at z = L and traveled backwards to the position z while the fourth term isthis wave reflected off the z = 0 end. Each additional factor of exp( − Lnk ( ω )) from the summation then representseach of these waves reflecting off both ends of the vessel to return to position z from the same direction from whichthey each originally approached.We can apply this paradigm of reflecting waves to the final form of Eq. S25, the solution for the single boundarycondition I ( z = 0 , t ) = Θ( t ) with no reflections, to produce the solution for the finite vessel case in which H ( t ) = I ( z = 0 , t ) = I ( z = L, t ) = Θ( t ): I ( z, t ) = ∞ X n =0 (cid:16) Θ (cid:0) t − ( z + 2 nL ) (cid:1) e iωt e − ( z +2 nL ) k ( ω ) − Θ (cid:0) t − (2 L − z + 2 nL ) (cid:1) e iωt e − (2 L − z +2 nL ) k ( ω ) +Θ (cid:0) t − ( L − z + 2 nL ) (cid:1) e iωt e − ( L − z +2 nL ) k ( ω ) − Θ (cid:0) t − ( L + z + 2 nL ) (cid:1) e iωt e − ( L + z +2 nL ) k ( ω ) (cid:17) + ∞ X n =0 Θ (cid:0) t − ( z + 2 nL ) (cid:1) e − t πi I dx x − x (1 − x ) − iωx e t ( x + x − ) e z +2 nL ( x − x − ) − Θ (cid:0) t − (2 L − z + 2 nL ) (cid:1) e − t πi I dx x − x (1 − x ) − iωx e t ( x + x − ) e L − z +2 nL ( x − x − )+ Θ (cid:0) t − ( L − z + 2 nL ) (cid:1) e − t πi I dx x − x (1 − x ) − iωx e t ( x + x − ) e L − z +2 nL ( x − x − ) − Θ (cid:0) t − ( L + z + 2 nL ) (cid:1) e − t πi I dx x − x (1 − x ) − iωx e t ( x + x − ) e L + z +2 nL ( x − x − ) ! = ∞ X n =0 (cid:18) Θ (cid:0) t − ( z + nL ) (cid:1) e iωt e − zk ( ω ) (cid:16) − e − Lk ( ω ) (cid:17) n + Θ (cid:0) t − ( L − z + nL ) (cid:1) e iωt e − ( L − z ) k ( ω ) (cid:16) − e − Lk ( ω ) (cid:17) n (cid:19) + ∞ X n =0 Θ (cid:0) t − ( z + nL ) (cid:1) e − t πi I dx x − x (1 − x ) − iωx e t ( x + x − ) e z ( x − x − ) (cid:16) − e L ( x − x − ) (cid:17) n +Θ (cid:0) t − ( L − z + nL ) (cid:1) e − t πi I dx x − x (1 − x ) − iωx e t ( x + x − ) e L − z ( x − x − ) (cid:16) − e L ( x − x − ) (cid:17) n ! . (S31)The second equality in Eq. S31 is obtained by combining the first and fourth terms of each sum as well as the secondand third terms. This is done by noting that the first and fourth (second and third) terms differ only by a factorof − L in the spatial arguments, thus allowing them to be combined by making the substitution2 n → n and introducing a factor of ( − n .From here, we define the function N ( z, t ) to be the number of waves that have passed through the point z at time t coming from the z = 0 end of the vessel. Specifically, N ( z, t ) = − t < z while for all t ≥ z it is the largestpossible integer that satisfies t ≥ z + N ( z, t ) L . This allows the Θ-functions in Eq. S31 to be removed by simplyrestricting n to sum from 0 to N ( z, t ) and N ( L − z, t ) respectively. The sums over n are then simple geometric sumsthat can be easily evaluated to express Eq. S31 as I ( z, t ) = N ( z,t ) X n =0 e iωt e − zk ( ω ) (cid:16) − e − Lk ( ω ) (cid:17) n + N ( L − z,t ) X n =0 e iωt e − ( L − z ) k ( ω ) (cid:16) − e − Lk ( ω ) (cid:17) n + N ( z,t ) X n =0 e − t πi I dx x − x (1 − x ) − iωx e t ( x + x − ) e z ( x − x − ) (cid:16) − e L ( x − x − ) (cid:17) n + N ( L − z,t ) X n =0 e − t πi I dx x − x (1 − x ) − iωx e t ( x + x − ) e L − z ( x − x − ) (cid:16) − e L ( x − x − ) (cid:17) n = e iωt e − zk ( ω ) − (cid:16) − e − Lk ( ω ) (cid:17) N ( z,t )+1 e − Lk ( ω ) + e iωt e − ( L − z ) k ( ω ) − (cid:16) − e − Lk ( ω ) (cid:17) N ( L − z,t )+1 e − Lk ( ω ) + e − t πi I dx x − x (1 − x ) − iωx e t ( x + x − ) e z ( x − x − ) 1 − (cid:16) − e L ( x − x − ) (cid:17) N ( z,t )+1 e L ( x − x − )+ e − t πi I dx x − x (1 − x ) − iωx e t ( x + x − ) e L − z ( x − x − ) 1 − (cid:16) − e L ( x − x − ) (cid:17) N ( L − z,t )+1 e L ( x − x − ) . (S32)From here, Eq. S29 can be used to simplify part of Eq. S32 via the relation e − zk ( ω ) e − Lk ( ω ) + e − ( L − z ) k ( ω ) e − Lk ( ω ) = sinh (cid:0) zk ( ω ) (cid:1) + sinh (cid:0) ( L − z ) k ( ω ) (cid:1) sinh (cid:0) Lk ( ω ) (cid:1) . (S33)To simplify the integral terms of Eq. S32, we will first separate each integral by splitting the factors of the form 1 − ( − exp( L ( x − x − ) / N ( z,t )+1 into the “constant” term, 1, and the “ N -dependent” term, ( − exp( L ( x − x − ) / N ( z,t )+1 .Focusing on the constant terms and ignoring the N -dependent terms for the time being, we will first add a copy ofeach integral while introducing a factor of 1 / x → x − in each of the copies. Since x is integrated counter-clockwise around theunit circle before the substitution, x will necessarily be integrated clockwise around the unit circle after. The directionof integration can be changed back to counter-clockwise though by absorbing the negative from dx/x → − dx/x , thusleaving all integration paths to be counter-clockwise around the unit circle. This provides the relation e − t πi I dx x − x (1 − x ) − iωx e t ( x + x − ) e z ( x − x − )1 + e L ( x − x − ) + e − t πi I dx x − x (1 − x ) − iωx e t ( x + x − ) e L − z ( x − x − )1 + e L ( x − x − )= e − t πi I dx x − x (1 − x ) − iωx e t ( x + x − ) e z ( x − x − )1 + e L ( x − x − ) + e − t πi I dx x − x − (1 − x − ) − iωx − e t ( x − + x ) e z ( x − − x )1 + e L ( x − − x )+ e − t πi I dx x − x (1 − x ) − iωx e t ( x + x − ) e L − z ( x − x − )1 + e L ( x − x − ) + e − t πi I dx x − x − (1 − x − ) − iωx − e t ( x − + x ) e L − z ( x − − x )1 + e L ( x − − x ) . (S34)In this form, the first and fourth integrals can be seen to perfectly cancel each other due to the relation1 x − x − (1 − x − ) − iωx − e t ( x − + x ) e L − z ( x − − x )1 + e L ( x − − x ) = − x − x (1 − x ) − iωx e t ( x + x − ) e z ( x − x − )1 + e L ( x − x − ) . (S35)Since the integrands of the first and fourth integrals in Eq. S34 are exact negatives of each other, their sum mustvanish. The second and third integrals will also cancel each other in a similar manner.After cancelling the constant terms and applying Eq. S33, Eq. S32 can be reexpressed as I ( z, t ) = e iωt sinh (cid:0) zk ( ω ) (cid:1) + sinh (cid:0) ( L − z ) k ( ω ) (cid:1) sinh (cid:0) Lk ( ω ) (cid:1) + ( − N ( z,t ) e iωt e − (cid:16) z + L ( N ( z,t )+1 ) (cid:17) k ( ω ) e − Lk ( ω ) + e − t ( − N ( z,t ) πi I dx x − x (1 − x ) − iωx e t ( x + x − ) e (cid:16) z + L ( N ( z,t )+1 ) (cid:17) ( x − x − )1 + e L ( x − x − )+ ( − N ( L − z,t ) e iωt e − (cid:16) L − z + L ( N ( L − z,t )+1 ) (cid:17) k ( ω ) e − Lk ( ω ) + e − t ( − N ( L − z,t ) πi I dx x − x (1 − x ) − iωx e t ( x + x − ) e (cid:16) L − z + L ( N ( z,t )+1 ) (cid:17) ( x − x − )1 + e L ( x − x − ) . (S36)The remaining integral terms found in Eq. S36 have poles at x = 0, x = ‘ ( ω ), and 1 + exp( L ( x − x − ) /
2) = 0. Theessential singularity at x = 0 will be ignored as (1 + exp( L ( x − x − ) / − has no well defined Laurent series. Thepole at x = ‘ ( ω ) can be easily accounted for via e iωt e − (cid:16) z + L ( N ( z,t )+1 ) (cid:17) k ( ω ) e − Lk ( ω ) + lim x → ‘ ( ω ) e − t x − x x − ‘ − ( ω ) e t ( x + x − ) e (cid:16) z + L ( N ( z,t )+1 ) (cid:17) ( x − x − )1 + e L ( x − x − ) = 0 . (S37)A similar cancellation exists for the z → L − z terms. This cancellation is to be expected given that the pole at x = ‘ ( ω ) was unnecessarily included by our choice to make the contour of integration be the unit circle.This leaves only the poles at 1 + exp( L ( x − x − ) /
2) = 0, which occur whenever L ( x − x − ) / iπ (2 m + 1) for anyinteger m . This condition is in turn met whenever x = x ± m , where x ± m = iπL (2 m + 1) ± s − (cid:18) πL (2 m + 1) (cid:19) . (S38)The various possible values of x ± m can be separated into three categories. First, when (cid:12)(cid:12) π (2 m + 1) /L (cid:12)(cid:12) < x ± m exists onthe unit circle itself along with x ±− ( m +1) . Second, when (cid:12)(cid:12) π (2 m + 1) /L (cid:12)(cid:12) > x ± m is purely imaginary with x − m existinginside (outside) the unit circle and x + m existing outside (inside) for m ≥ m < (cid:12)(cid:12) π (2 m + 1) /L (cid:12)(cid:12) = 1, x ± m is exactly i for m ≥ − i for m < / (1 + exp( L ( x − x − ) / x is near x ± m by11 + e L ( x − x − ) ≈ x − x ± m lim x → x ± m x − x ± m e L ( x − x − ) = 1 x − x ± m lim x → x ± m L x ( x + x − ) e L ( x − x − )= − x − x ± m x ± mL (cid:18) x ± m + (cid:16) x ± m (cid:17) − (cid:19) , (S39)where L’Hopital’s rule has been used to determine the limit. In the third category in which x ± m = ± i , the pole is ofsecond order and 1 / (1 + exp( L ( x − x − ) / e L ( x − x − ) ≈ (cid:16) x − x ± m (cid:17) lim x → x ± m (cid:0) x − x ± m (cid:1) e L ( x − x − ) = 1 (cid:16) x − x ± m (cid:17) lim x → x ± m (cid:18)(cid:16) L x ( x + x − ) (cid:17) − Lx − (cid:19) e L ( x − x − )= 2 (cid:0) x ± m (cid:1) L (cid:16) x − x ± m (cid:17) . (S40)0To determine if there is still a residual first order pole we expand Eq. S40 further to11 + e L ( x − x − ) ≈ (cid:0) x ± m (cid:1) L (cid:16) x − x ± m (cid:17) + 1 x − x ± m lim x → x ± m (cid:0) x − x ± m (cid:1)
11 + e L ( x − x − ) − (cid:0) x ± m (cid:1) L (cid:16) x − x ± m (cid:17) = 2 (cid:0) x ± m (cid:1) L (cid:16) x − x ± m (cid:17) + 1 x − x ± m lim x → x ± m L (cid:0) x − x ± m (cid:1) − (cid:0) x ± m (cid:1) (cid:16) e L ( x − x − ) (cid:17) L (cid:16) x − x ± m (cid:17) (cid:16) e L ( x − x − ) (cid:17) = 2 (cid:0) x ± m (cid:1) L (cid:16) x − x ± m (cid:17) + 1 x − x ± m lim x → x ± m L (cid:0) x − x ± m (cid:1) − L (cid:0) x ± m (cid:1) (cid:0) x − (cid:1) e L ( x − x − ) L (cid:16) e L ( x − x − ) (cid:17) + L (cid:16) x − x ± m (cid:17) (1 + x − ) e L ( x − x − )= 2 (cid:0) x ± m (cid:1) L (cid:16) x − x ± m (cid:17) + 1 x − x ± m lim x → x ± m L − L (cid:0) x ± m (cid:1) (cid:16) L (cid:0) x − (cid:1) − x − (cid:17) e L ( x − x − ) L (1 + x − ) e L ( x − x − ) + L (cid:16) x − x ± m (cid:17) (cid:16) L (1 + x − ) − x − (cid:17) e L ( x − x − )= 2 (cid:0) x ± m (cid:1) L (cid:16) x − x ± m (cid:17) + 1 x − x ± m · lim x → x ± m − L (cid:0) x ± m (cid:1) (cid:16) L (cid:0) x − (cid:1) − Lx − (cid:0) x − (cid:1) + 6 x − (cid:17) e L ( x − x − ) L (cid:18) L (1 + x − ) − x − + (cid:16) x − x ± m (cid:17) (cid:16) L (1 + x − ) − Lx − (1 + x − ) + 6 x − (cid:17)(cid:19) e L ( x − x − )= 2 (cid:0) x ± m (cid:1) L (cid:16) x − x ± m (cid:17) + 2 (cid:0) x ± m (cid:1) L (cid:16) x − x ± m (cid:17) . (S41)Thus, we see that 1 / (1 + exp( L ( x − x − ) / x ± m = ± i .We now explore the on contour case of (cid:12)(cid:12) π (2 m + 1) /L (cid:12)(cid:12) <
1. For the time being we will assume m ≥
0. In this regime( x ± m ) − = ( x ± m ) ∗ . We can thus use Eq. S39 to determine the residual value y ± m = lim x → x ± m e − t x − x (1 − x ) − iωx e t ( x + x − ) e (cid:16) z + L ( N ( z,t )+1 ) (cid:17) ( x − x − ) − x ± mL (cid:18) x ± m + (cid:16) x ± m (cid:17) − (cid:19) = e − t i Im (cid:0) x ± m (cid:1) Re (cid:16) x ± m (cid:17) − − iω e t Re ( x ± m ) e i (cid:16) z + L ( N ( z,t )+1 ) (cid:17) Im ( x ± m ) 1 L Re (cid:16) x ± m (cid:17) . (S42)By noting that Re( x ± m ) = Re( x ±− ( m +1) ) and Im( x ± m ) = − Im( x ±− ( m +1) ) we can next evaluate the sum y ± m + y ±− ( m +1) = − e − t e t Re ( x ± m )Im (cid:0) x ± m (cid:1) L Re (cid:16) x ± m (cid:17) (cid:18) Re (cid:16) x ± m (cid:17) − − iω (cid:19) sin (cid:18)(cid:16) z + L (cid:0) N ( z, t ) + 1 (cid:1)(cid:17) Im (cid:0) x ± m (cid:1)(cid:19) . (S43)Similarly, by noting that Re( x + m ) = − Re( x − m ) and Im( x + m ) = Im( x − m ) we can evaluate the sum1 y + m + y + − ( m +1) + y − m + y −− ( m +1) = − e − t e t Re ( x + m )Im (cid:0) x + m (cid:1) (cid:16) Re (cid:0) x + m (cid:1) + 1 + iω (cid:17) L Re (cid:16) x + m (cid:17) (cid:18) Re (cid:16) x + m (cid:17)(cid:19) − (1 + iω ) ! sin (cid:18)(cid:16) z + L (cid:0) N ( z, t ) + 1 (cid:1)(cid:17) Im (cid:0) x + m (cid:1)(cid:19) + e − t e − t Re ( x + m )Im (cid:0) x + m (cid:1) (cid:16) − Re (cid:0) x + m (cid:1) + 1 + iω (cid:17) L Re (cid:16) x + m (cid:17) (cid:18) Re (cid:16) x + m (cid:17)(cid:19) − (1 + iω ) ! sin (cid:18)(cid:16) z + L (cid:0) N ( z, t ) + 1 (cid:1)(cid:17) Im (cid:0) x + m (cid:1)(cid:19) = 2 e − t Im (cid:0) x + m (cid:1) sin (cid:18)(cid:16) z + L (cid:0) N ( z, t ) + 1 (cid:1)(cid:17) Im (cid:0) x + m (cid:1)(cid:19) L (1 + iω ) − (cid:18) Re (cid:16) x + m (cid:17)(cid:19) ! cosh (cid:16) t Re (cid:0) x + m (cid:1)(cid:17) + 1 + iω Re (cid:16) x + m (cid:17) sinh (cid:16) t Re (cid:0) x + m (cid:1)(cid:17) = − ( − (2 m +1) N ( z,t ) π (2 m + 1) e − t sin (cid:0) πzL (2 m + 1) (cid:1) π (2 m + 1) + (cid:0) Lk ( ω ) (cid:1) · cosh t s − (cid:18) πL (2 m + 1) (cid:19) + 1 + iω q − (cid:0) πL (2 m + 1) (cid:1) sinh t s − (cid:18) πL (2 m + 1) (cid:19) , (S44)where Eq. S38 and the fact that sin( x + nπ ) = ( − n sin( x ) for any integer n were used to simplify the final form.We next explore the case of (cid:12)(cid:12) π (2 m + 1) /L (cid:12)(cid:12) >
1, again assuming m ≥
0. In this regime, x − m is contained within theunit circle while x + m is not, and thus only the pole at x − m need be considered between the two. To aid in calculatingthis pole, we will use the relations x − m + (cid:0) x − m (cid:1) − = 1 − (cid:18) πL (2 m + 1) − q(cid:0) πL (2 m + 1) (cid:1) − (cid:19) i (cid:18) πL (2 m + 1) − q(cid:0) πL (2 m + 1) (cid:1) − (cid:19) = − i s(cid:18) πL (2 m + 1) (cid:19) − , (S45a) x − m − (cid:0) x − m (cid:1) − = − (cid:18) πL (2 m + 1) − q(cid:0) πL (2 m + 1) (cid:1) − (cid:19) i (cid:18) πL (2 m + 1) − q(cid:0) πL (2 m + 1) (cid:1) − (cid:19) = 2 i πL (2 m + 1) . (S45b)With these, we can more easily calculate the residual value y − m = lim x → x − m e − t x − x (1 − x ) − iωx e t ( x + x − ) e (cid:16) z + L ( N ( z,t )+1 ) (cid:17) ( x − x − ) − x − mL (cid:18) x − m + (cid:16) x − m (cid:17) − (cid:19) = e − t iπL (2 m + 1) Li q(cid:0) πL (2 m + 1) (cid:1) − (cid:18) i q(cid:0) πL (2 m + 1) (cid:1) − iω (cid:19) e − it q ( πL (2 m +1) ) − e (cid:16) z + L ( N ( z,t )+1 ) (cid:17) iπL (2 m +1) = − ( − (2 m +1) N ( z,t ) π (2 m + 1) e − t (cid:18) iω − i q(cid:0) πL (2 m + 1) (cid:1) − (cid:19)(cid:16) π (2 m + 1) + (cid:0) Lk ( ω ) (cid:1) (cid:17) q(cid:0) πL (2 m + 1) (cid:1) − e − it q ( πL (2 m +1) ) − e iπzL (2 m +1) . (S46)2From here we note that in this regime ( x − m ) ∗ = x + − ( m +1) . Thus, the residual at x + − ( m +1) must simply be the conjugateof Eq. S46 with the exception that the term iω not be conjugated since it originated from the boundary conditionsas opposed to the location of the pole. The total of these two residuals thus takes the form y − m + y + − ( m +1) = − ( − (2 m +1) N ( z,t ) π (2 m + 1) e − t (cid:18) iω − i q(cid:0) πL (2 m + 1) (cid:1) − (cid:19)(cid:16) π (2 m + 1) + (cid:0) Lk ( ω ) (cid:1) (cid:17) q(cid:0) πL (2 m + 1) (cid:1) − e − it q ( πL (2 m +1) ) − e iπzL (2 m +1) − ( − (2 m +1) N ( z,t ) π (2 m + 1) e − t (cid:18) iω + i q(cid:0) πL (2 m + 1) (cid:1) − (cid:19)(cid:16) π (2 m + 1) + (cid:0) Lk ( ω ) (cid:1) (cid:17) q(cid:0) πL (2 m + 1) (cid:1) − e it q ( πL (2 m +1) ) − e − iπzL (2 m +1) = − ( − (2 m +1) N ( z,t ) π (2 m + 1) e − t π (2 m + 1) + (cid:0) Lk ( ω ) (cid:1) sin πzL (2 m + 1) − t s(cid:18) πL (2 m + 1) (cid:19) − + 1 + iω q(cid:0) πL (2 m + 1) (cid:1) − πzL (2 m + 1) − t s(cid:18) πL (2 m + 1) (cid:19) − . (S47)Finally, we consider the case of (cid:12)(cid:12) π (2 m + 1) /L (cid:12)(cid:12) = 1. We will continue to assume m ≥ L can be freely replacedwith π (2 m + 1). Since the pole in this case is on the contour and a combination of a second and first order pole, wecan use Eq. S41 to calculate the residual y m = lim x → i i π (2 m + 1) + 2 i π (2 m + 1) ∂∂x ! e − t x − x (1 − x ) − iωx e t ( x + x − ) e (cid:16) z + π (2 m +1) ( N ( z,t )+1 ) (cid:17) ( x − x − )= ( − (2 m +1) N ( z,t ) e − t e iz π (2 m + 1) (1 + iω )+ i e − t π (2 m + 1) lim x → i t (cid:0) − x − (cid:1) (cid:0) − x (cid:1) x (cid:16) (1 − x ) − iωx (cid:17) e t ( x + x − ) e (cid:16) z + π (2 m +1) ( N ( z,t )+1 ) (cid:17) ( x − x − )+ i e − t π (2 m + 1) lim x → i (cid:16) z + π (2 m + 1) (cid:0) N ( z, t ) + 1 (cid:1)(cid:17) (cid:0) x − (cid:1) (cid:0) − x (cid:1) x (cid:16) (1 − x ) − iωx (cid:17) e t ( x + x − ) e (cid:16) z + π (2 m +1) ( N ( z,t )+1 ) (cid:17) ( x − x − )+ i e − t π (2 m + 1) lim x → i − xx (cid:16) (1 − x ) − iωx (cid:17) e t ( x + x − ) e (cid:16) z + π (2 m +1) ( N ( z,t )+1 ) (cid:17) ( x − x − ) − i e − t π (2 m + 1) lim x → i (cid:0) − x (cid:1) (cid:0) x − x (1 + iω ) (cid:1) x (cid:16) (1 − x ) − iωx (cid:17) e t ( x + x − ) e (cid:16) z + π (2 m +1) ( N ( z,t )+1 ) (cid:17) ( x − x − )= − ( − (2 m +1) N ( z,t ) e − t e iz iπ (2 m + 1) (1 + iω ) (cid:0) t (1 + iω ) (cid:1) . (S48)We then calculate the residual from the pole at x = − i by again conjugating everything except for the iω terms inEq. S48, resulting in the sum y m + y − ( m +1) = − ( − (2 m +1) N ( z,t ) e − t sin ( z ) π (2 m + 1) (1 + iω ) (cid:0) t (1 + iω ) (cid:1) . (S49)Of note is the fact that when L = π (2 m + 1), the expression π (2 m + 1) (1 + iω ) can also be written as π (2 m +1) + ( Lk ( ω )) so as to make the prefactor of Eq. S49 appear more similar to those of Eqs. S44 and S47.3We can now fully express the integral terms found in Eq. S36, less the incalculable singularity at x = 0. Bycombining Eqs. S37, S44, S47, and S49 we can show( − N ( z,t ) e iωt e − (cid:16) z + L ( N ( z,t )+1 ) (cid:17) k ( ω ) e − Lk ( ω ) + e − t ( − N ( z,t ) πi I dx x − x (1 − x ) − iωx e t ( x + x − ) e (cid:16) z + L ( N ( z,t )+1 ) (cid:17) ( x − x − )1 + e L ( x − x − )= − e − t ∞ X m =0 π (2 m + 1) π (2 m + 1) + (cid:0) Lk ( ω ) (cid:1) sin (cid:18) πzL (2 m + 1) (cid:19) Ω m ( t, L, ω ) − cos (cid:18) πzL (2 m + 1) (cid:19) Υ m ( t, L, ω ) ! , (S50)whereΩ m ( t, L, ω ) = cosh (cid:18) t q − (cid:0) πL (2 m + 1) (cid:1) (cid:19) + iω q − ( πL (2 m +1) ) sinh (cid:18) t q − (cid:0) πL (2 m + 1) (cid:1) (cid:19) πL (2 m + 1) <