1. Introduction

Our study concerns modeling and simulation of wave propagation along spatial networks, inspired by (and with intended applications to) modeling blood flow in cardiovascular networks [16]. There are also connections with the modeling of traffic flow in urban networks, supply chains and telecommunication data networks (see [10], [11], [12], [13], [19]). For a spatial network, the dynamics is typically described first at the level of individual edges, followed by a model for the junctions; the simulations can be performed on parts of the network or, if feasible, on the entire network. Current network models can use as many as $87$-edge trees, with the more common ones using $55$-edge trees, followed at each of its terminal sites (organs) by networks of small arteries, arterioles and capillaries, with as many as $23$ generations (see e.g. [25] [34]). Capillary networks typically are no longer tree-like (e.g. [8], [32]), since their role is to irrigate all tissues and organs in the body. Understanding the dynamics in such terminal networks is critical due to the role they play in prescribing the boundary conditions for the system of large arteries. The vascular system also include the venous system, whose geometry resembles (to some extent) the arterial network, but with notable differences both in network geometry and in constitutive properties. The modeling of the venous system presents different challenges due to its physiological conditions: low pressure gradients but high fluid volume, collapsible tubes, presence of venous valves [21]. The loop is completed with the inclusion of the heart chambers, heart valves and the pulmonary circulation [27].

Available computational power (nowadays and for the foreseeable future) limits the simulation of such complex network in its entirety. Still, a system-level analysis and simulation is a desirable task, due to non-local phenomena (such as hypertension or autonomic regulation) that cannot be explained or replicated by only considering portions of the network. One remedy is as follows: when only a portion of the large network is analyzed in detail (such as the large arterial network, or the vascularization of an individual organ), the rest of the network can be resolved at a coarse level, through reduced models (see e.g. [14], [22], [24], [30], [33]). This is when the issue of appropriate boundary conditions arises. Inflow and outflow conditions must realistically mimic the behavior of the network which has been removed from the model. If a full $3$D model is employed in part of the network, then one usually employs reduced $1$D or $0$D models which are coupled at the inflow/outflow. Recent studies (e.g. [20], [26]) have introduced new optimization techniques for computing optimal boundary conditions to match experimental data. The issue of modeling appropriate boundary conditions remains crucial. In this study we restrict ourselves to using the simplest boundary conditions with the sole goal of emphasizing the influence of the boundary conditions to the systemic dynamics.

Under normal conditions, the physical system being modeled is in dynamic equilibrium -with the heart beating regularly and the pressure varying from high (systolic) to low (diastolic) values throughout the network, in a quasi-periodic fashion. In presence of disturbances, such as a blockage or a release thereof, the response of the system is to return to its equilibrium, known as homeostasis, in an optimal fashion. In the network of large arteries, the process of autoregulation can be modeled using boundary controls (at the root of the network -the heart and at peripheral nodes -the ends of the large arteries), driving the system (or a part of it) back to a dynamic equilibrium in minimum time. A major factor in the controllability of the vascular network is the peripheral resistance, which is in turn a result of dynamics on much more complicated network (microcirculation at the level of individual organs). Such cardiovascular control analysis using 0D models have been described in the literature (e.g. [4]).

In this paper we focus on 1D models of such vascular networks, based on an algebraic relationship between pressure and cross-section area. This approach has already been validated in the literature (see [1], [28], [2], [29], [31]). These 1D models can also accommodate higher order terms (see e.g. [5], [6], [14]) to model visco-elastic and/or dispersive effect, for instance. Numerical discretisations of such equations (using discontinuous Galerkin techniques) can be found in [9]. The novelty here is the use of these models to study optimization tasks that are relevant in the autonomic regulation mechanism present in the physiological settings.

The paper is organized as follows. In Section 2 we introduce the governing equations of the model. Section 3 is devoted to a description of the boundary conditions used in our modeling and the connection with the Riemann problem. In Section 4 the results of several numerical studies problems are presented and discussed. The numerical discretization of the PDEs system using a discontinuous Galerkin formulation coupled with a two-steps method of Adam-Bashforth is detailed in the appendix A.

---

2. Mathematical model

In the present work we study the flow and pressure waves in a network of fluid-filled tubes with elastic walls. The focus is on various types of networks. Reduced models have been extensively employed in the literature ([24], [26], [29]). Here the starting points are the $1$D models for the pressure and flow velocity in an arterial network. Denote $A=A(x,t)$, $U=U(x,t)$, $P=P(x,t)$ and $f=f(x,t)$ the cross section area, average flow velocity, blood pressure and friction force for unit length, respectively, at location $x$ at time $t$. The model we consider here is the following:

$\begin{aligned}\label{F1} \frac{\partial A}{\partial t}+\frac{\partial (AU)}{\partial x}&=0, \\ \frac{\partial U}{\partial t}+U\frac{\partial U}{\partial x}+\frac{1}{\rho}\frac{\partial P}{\partial x}&=f, \end{aligned}$

(1)

where the first equation represents the mass conservation in the network while the second is consequence of the Navier-Stokes equations under some assumption on the flow velocity profile across a cross section. $\rho$ is the density of the blood. Throughout this paper we assume a parabolic velocity profile, that yields the friction term $f =-22\mu \pi U/A$, where $\mu$ is the fluid viscosity. We note that some papers have been devoted to improving these equations, for instance by removing the assumption on the flow velocity profile (see e.g. [5]).

The arteries (edges of the network) are modeled as fluid-filled tubes with elastic walls (with Young modulus $E$, wall thickness $h_w$, wall density $\rho_w$ and unstressed radius $r_0$), in which the hydrodynamic pressure $P$ is assumed constant across the vessel cross section and determined by the properties of the wall, measured by the wall displacement $\eta =\eta (x,t)$. Here we use a linear elastic (algebraic) model, in which the pressure $P$ is linearly dependent on $\eta$:

$P=P_{ext}+\frac{\beta }{r_{0}^{2}}\eta,$

where $P_{ext}$ is the external pressure and $\beta =\frac{Eh_w}{1-\sigma ^{2}}\sqrt{\pi}$ ($\sigma$ is the Poisson ratio, usually taken to be $\sigma =\frac{1 }{2}$). The pressure dependence on the cross section area $A$ is then

$P=P_{ext}+\frac{\beta}{A_{0}}(\sqrt{A}-\sqrt{A_{0}}), \label{F2}$

(2)

where $A_{0}=\pi r_{0}^{2}$ is the unstressed cross section area. When considering additional effects, such as visco-elasticity or wall inertia, the pressure term may include also wall velocity $\eta_t$ and wall acceleration $\eta_{tt}$, respectively, which leads to different fluid-structure interaction models. One such model is reported elsewhere [7].

---

3. Boundary conditions

In this section we describe the boundary conditions considered in our study. The model is solved using a discontinuous Galerkin numerical scheme, as reported in appendix A.

---

3.1. The Riemann problem

Due to the nature of the discontinuous Galerkin scheme used, we must specify two states (left and right) for each of the variables considered: $(A_{L},U_{L})$ and $(A_{R},U_{R})$. We first set up the Riemann Problem at each interface.

At a time $t$, each interface separates two constant states, $(A_{L},U_{L})$ and $(A_{R},U_{R})$, and we need to determine the two upwinded states, $(A_{L}^{u},U_{L}^{u})$ and $(A_{R}^{u},U_{R}^{u})$, originated on each side of interface at time $t+\Delta t$. To do this, the following equations need to be solved:

$W_{f}(A_{L},U_{L})=W_{f}(A_{L}^{u},U_{L}^{u}), \\ W_{b}(A_{R},U_{R})=W_{b}(A_{R}^{u},U_{R}^{u}), \\ A_{L}^{u}U_{L}^{u}=A_{R}^{u}U_{R}^{u}, \\ \rho \frac{(U_{L}^{u})^{2}}{2}+P(A_{L}^{u})=\rho \frac{(U_{R}^{u})^{2}}{2} +P(A_{R}^{u}).$

(3)

The first two equations come from the assumption that the flow between two initial states is inviscid, and the forward characteristic information, $W_{f}$, and the backward characteristic information, $W_{b}$, are given by

$W_{f} =U+4(c-c_{0}), \\ W_{b} =U-4(c-c_{0}),$

with

$c=\sqrt{\frac{\beta }{2\rho A_{0}}}A^{1/4},\text{ }c_{0}=\sqrt{\frac{ \beta }{2\rho }}A_{0}^{-1/4}.$

The remaining equations follow from conservation of mass and continuity of the total pressure at the interface. To obtain $(A_{L}^{u},U_{L}^{u})$ and $(A_{R}^{u},U_{R}^{u})$ an iterative Newton-Raphson method is employed.

The boundary conditions can be classified in three types, depending on the location in the domain of the arterial network: inflow, junction and terminal boundary condition.

---

3.2. Inflow boundary conditions

At the inflow, we have the option to prescribe an inflow time-dependent area, $A_{bc}(t)$, velocity, $U_{bc}(t)$, or flow rate, $Q_{bc}(t)$. Thus, at the inlet of the arterial domain, we prescribe one of the following

$A_{bc}(t):\left\{ \begin{array}{l} U_{L}=U_{R}, \\ A_{L}=\left( 2(A_{bc})^{\frac{1}{4}}-(A_{R})^{\frac{1}{4}}\right) ^{4}, \end{array} \right. \\ U_{bc}(t):\left\{ \begin{array}{l} U_{L}=2U_{bc}-U_{R}, \\ A_{L}=A_{R}, \end{array} \right. \\ Q_{bc}(t):\left\{ \begin{array}{l} U_{L}=2\frac{Q_{bc}}{A_{R}}-U_{R}, \\ A_{L}=A_{R}. \end{array} \right.$

In the real vascular system, the correct inflow conditions are dictated by the heart model that is being used when coupled with systemic network. Numerical simulations in the next section are reported for the prescribed flow rate $Q_{bc}$, with the mention that comparable results were obtained when the other inflow conditions were used.

---

3.3. Junction boundary conditions

The spatial network models allow for junctions with arbitrarily large degrees ($N$). For $N=3$ the possible types of junctions are 1$\rightarrow$2 and 2$\rightarrow$1.

Note that junctions with $N=2$, which have an incoming and an outgoing edge, are modeled using the Riemann problem in section 3.1. In this case, the upwinded states are determined through the numerical solution of (3).

If $\left( A_{1},P_{1},U_{1}\right)$, $\left( A_{2},P_{2},U_{2}\right)$ and $\left( A_{3},P_{3},U_{3}\right)$ are the initial states at the points of each elemental region adjacent to the junction, then the upwinded states $\left( A_{i}^{u},P_{i}^{u},\right.\left.U_{i}^{u}\right)$ $(i=1,2,3)$ are determined by solving non-linear systems. For example for a junction of type 1$\rightarrow$2:

$W_{f}(A_{1}^{u},U_{1}^{u})=W_{f}(A_{1},U_{1}),$

$W_{b}(A_{2}^{u},U_{2}^{u})=W_{b}(A_{2},U_{2}),$

$W_{b}(A_{3}^{u},U_{3}^{u})=W_{b}(A_{3},U_{3}),$

$A_{1}^{u}U_{1}^{u}=A_{2}^{u}U_{2}^{u}+A_{3}^{u}U_{3}^{u},$

$P(A_{1}^{u})+\frac{1}{2}\rho (U_{1}^{u})^{2}=P(A_{2}^{u})+\frac{1}{2}\rho (U_{2}^{u})^{2},$

$P(A_{1}^{u})+\frac{1}{2}\rho (U_{1}^{u})^{2}=P(A_{3}^{u})+\frac{1}{2}\rho (U_{3}^{u})^{2},$

where $W_b$ and $W_f$ are the characteristic waves obtained by the solution of the Riemann problem. Similar relations are imposed for the other types of junctions illustrated in Figure 2.

---

3.4. Terminal boundary conditions

For terminal conditions, the initial state $(A_{L},U_{L})$ is located at the end point of the arterial domain and the initial state $(A_{R},U_{R})$ is selected to produce an upwinded state $(A^{u},U^{u})$ for the next time step. We can distinguish between different terminal conditions (see [3]) considering an analogy with electrical circuits.

• Pure resistance condition: Assume a model with terminal reflection coefficient $R_{t}$, which is based on the assumption that $W_{b}$ is proportional to $W_{f}$:

$W_{b}=-R_{t} W_{f},$

where $-1\leq R_{t}\leq 1$. $R_t=1$ corresponds to a complete reflection of the characteristic (complete blockage in that terminal site, $U=0$); $R_t=0$ means there is no reflected characteristic at the terminal site, and $R_t=-1$ means that it is free terminal site ($c=c_0$). In this case we have

$A_{R} =A_{L}, \\ U_{R} =W_{f}(1-R_{t})-U_{L}.$

• RC model: we consider a lumped parameter model made of a resistance $R_{\mu }.$ We can compute $A^{u}$ by solving the nonlinear equation

$R_{\mu }\left( U_{L}+4(A_{L})^{\frac{1}{4}}\sqrt{\frac{\beta _{L}}{2\rho A_{0_{L}}}}\right) A^{u}-4R_{\mu }(A^{u})^{\frac{5}{4}}\sqrt{\frac{\beta _{L} }{2\rho A_{0_{L}}}}+ \\ -P_{0}-\frac{\beta _{L}}{A_{0_{L}}}(\sqrt{A^{u}}-\sqrt{A_{0_{L}}})+P_{out}=0,$

by a Newton's Method, starting from $A^{u}=A_{L},$ and

$U^{u}=\frac{P(A^{u})-P_{out}}{R_{\mu }A^{u}},$

where $P_{out}$ is the pressure at the end of the arterial domain and $R_{\mu }$ is the resistance. In this case, the terminal conditions become:

$A_{R} =A_{L}, \\ U_{R} =2U^{u}-U_{L}.$

• Compliance model: in this model we add a compliance C to the previous one. We need to compute $A^{u}=A_{in}^{n}$ by solving the recursive equation

$A_{in}^{n}=\left( \sqrt{A_{in}^{n-1}}+\frac{\Delta tA_{0_{L}}}{C\beta _{L}} \left( Q_{in}^{n-1}+\right. \right. \\ \left. \left. +\frac{1}{R_{\mu }} \left( P_{out}-P_{0}-\frac{\beta _{L}}{A_{0_{L}}}\left( \sqrt{A_{in}^{n-1}}-\sqrt{A_{0_{L}}}\right) \right) \right) \right) ^{2},$

where $C$ is the compliance, so the terminal conditions become

$A_{R} =\left( 2((A^{u})^{n})^{\frac{1}{4}}-(A_{L})^{\frac{1}{4}}\right) ^{4}, \\ U_{R} =U_{L}.$

The results in the next Section use the pure resistance model presented above. Simulations have also been performed using the RC model and the Compliance model, with similar outcomes. Since an additional iteration (Newton or recursive) is required for each time step of the simulation, thus increasing additional computational time, we chose not to pursue them in this study.

---

4. Results

In this section we present our simulation results addressing several scenarios: (1) the effect that truncation in a fractal network has on the flow in the root edge; (2) the effect that adding or subtracting an edge has to the network dynamics; and (3) optimization of the heart rate in the event of a blockage/unblockage of an edge or of an entire subtree.

For the simulations, we prescribe a periodic inflow (at the root of the network), with period $T=60/HR$, as $Q=Q_{bc}(t)$ in liters/sec, where

$Q_{bc}(t)=\begin{cases} \begin{aligned} &6.75\times 10^{-4}\sin \left(\frac{\pi t}{\tau }\right), &\text{ for } t\in[0,\tau],\\ &0,&\text{ for } t\in [\tau,T], \end{aligned} \end{cases}$

with $\tau =T/4$ a quarter of the heart beat period. The following parameter values are used throughout the sequel: $\mu =4\times 10^{-3}\text{ }Pa\cdot s$ (viscosity of the blood), $\rho =1050\text{ }Kg/m^{3}$ (blood density), $P_{ext} =1862\text{ }Pa$ (external pressure), and $\beta=1418\text{ } kg/s^2$.

Multiple runs were performed in an asymmetric tree with $2$ generations for various values of the resistance($R_{t}=0.3,R_{t}=0.8,R_{t}=1$). In the case of low resistance values the steady state is reached very quickly and no oscillations are noticed. When moderately high resistance values ($R_t= 0.8$) are used, slow oscillations can be observed (Figure 4) during the steady state phase ($\sim$ 0.4 Hz), while for maximum resistance ($R_t=1$) slow oscillations ($\sim$ 0.1 Hz) become present even during the pressure build-up stage (Figure 5).

The presence of slower oscillations resembles the physiological phenomenon of the Mayer waves present in the vascular system [23]. Mayer waves are oscillations at a much slower frequency (0.1 Hz) than the heart beat or even respiration (0.25-0.3 Hz), and are partly responsible for the variability of the heart rate, manifested in the autoregulation mechanism. Mayer's waves are sometimes attributed to the action of the nervous system. The exact frequencies (0.4 Hz and 0.1 Hz for the two resistance values used in simulations here $R_t=0.8$ and $R_t=1$, respectively) where these oscillations manifest is not meant to match the frequency of the Mayer waves observed in the real system. The latter simulation does not reach a steady state because of the complete blockage of the terminal edges, so this phenomenon is transitory in time. In reality, a controlling mechanisms dictates what resistance values are applied in order to regulate the pressure, and therefore the slow oscillations are also transitory in nature.

Our simulations suggest that these oscillations may be in part due to the spatial network itself, by triggering network-induced oscillations. In fact we go one step further and conjecture that the presence of these network induced oscillations at certain frequencies has determined the nervous system to tune its autonomic regulation around this 0.1 Hz frequency, and not the other way around. Such conjecture requires a much thorough experimental validations, which is outside the scope of this numerical study.

---

4.1. Effect of truncation of the network

We study the effect of truncation (by considering fewer generations of the same network) to the flow in the root edge. This may suggest the effectiveness of modifications in the terminal conditions for the simplified tree in order to mimic the behavior in the larger tree.

We construct a $2$-generation self-similar tree based on a simple junction (3-edge tree) with the physical sizes given in Table 1 and run the simulation using outflow conditions with $R_{t}=0.8$.

The length-to-radius ratio is not representative of what is reported in real systems and thus we are not concerned here with the variations in this ratio or in the sensitivity of the model to this ratio. This will be reported elsewhere.

For the truncation we use a $3$-edge tree with the following lengths and radii:

Note that the truncated tree has the same length and radii for the root edge, while the lengths of the children edges are chosen to match the lengths of the longest path and shortest path, respectively. The outflow area for the terminal edges is the same, so overall, the truncated tree has a smaller outflow area than the fractal tree.

We record first the flow in the middle of the root edge when the number of generations of the self-similar network change (from 0 to 2), using the same outflow conditions ($R_{t}$ value for the terminal resistance). The reference flow profile in the root edge was obtained for terminal resistance $R_{t}=0.3$. In Figure 6 temporal recordings for 0-generations and 2-generations are included respectively.

Notable differences are due to the change in the reflected waves. Even though there is some pressure increase from beat to beat, the flow is close to being periodic and therefore the comparisons can be made even at these stages of the simulation. We note that performing the simulation past the first 2 seconds does not provide any additional features relevant to this comparison, since a quasi-steady-state is achieved within the 1st second, so only the first couple of beats are reported here.

Next we perform numerical optimization on the value of the terminal resistance $R_{t}$ to be used for the outflow conditions in the simplified tree in order to best match (in least square sense) the flow in the middle of the root edge of the simplified tree with the flow in the middle of the root edge of the larger tree.

The result of the optimization, performed using MATLAB's fminsearch implementation of the derivative-free Nelder-Mead algorithm, shows that the optimal value for the resistance in the simplified tree is $R^*_{t}=-0.25$. Negative values of $R_t^*$ indicates that in the truncated tree the optimal terminal conditions are in such a way that there is no reflected backward characteristics, hence the resistance to flow out is much smaller than in the referenced 2-generation tree.

---

4.2. Effect of adding/removing edges in a network

We investigate the effect of adding or subtracting an edge to a given network. In physiology this phenomenon has been observed (see e.g. [8] where the brain vasculature is considered). Here we study the "efficiency" of the resulting network by measuring the total outflow of the network during a given period of time. For the simulations, we use a $7$-edge network with cycles, that is a network in which there are some nodes that are connected in a closed loop and then we remove the middle edge.

The inflow and outflow recordings for both networks are displayed below. Note that the at the inflow both networks exhibit the same behavior up until reflected waves return to the root of the network. Likewise, at the outflow the flow and pressure stay constant for the first 0.5 sec, indicating when the first waves arrive at the outflow.

For this very simple network, we compute a total outflow of $0.9552$ $cm^{3}$ in the presence of the central edge, while in its absence the total outflow is increased to $1.0435$ $cm^{3}$. The inflow during the $3$ sec simulation was $3.0706$ $cm^{3}$. These numerical values are obtained by integrating in time (over the 3 sec time interval) the flow $Q=AU$. Note that the outflow in both cases is roughly one third of the inflow, the difference being accounted in the volume of fluid that has been "stored" in the network itself, which caused the pressurization of the network at equilibrium.

The results of simulation show that the flow through the network is enhanced when fewer edges (cycles) are present, similar to Braess's paradox in traffic flow, which states that adding extra capacity to a network can in some cases reduce overall performance.

---

4.3. Effect of blockages in a network on the flow

Here we study the time-optimal problem of returning to basal flow and pressure conditions on a network after a temporary blockage of a subnetwork has been removed. As initial network we consider again the fractal tree with $15$ edges ($2$-generation tree) introduced in section 4.1. The length and radius ratios between parent and children edges are 0.9 and 0.8 respectively. We run a 20-sec simulation for the flow in the entire network before any blockage is applied (see Figure 3), with a terminal resistance model and $R_t=0.8$ for each of the terminal edges. The transition time for the flow to settle varies slightly with the heart rate, but it is significantly less than the 20 seconds; e.g. for constant heart rate of $HR=75$ beats/min, the periodic state settles after roughly the first 5 seconds.

Exactly at 20 seconds, we introduce an instantaneous blockage at the end of edge 3 (hence the subnetwork having edge 3 as root is also blocked off). This makes edge 3 a terminal edge (with $R_t=1$), while the other terminal edges remain at $R_t=0.8$. We run a 5-sec simulation with this modified network, keeping the heart rate the same as before blockage. After 5 seconds of complete blockage at end of edge 3 (so at second 25 in the total simulation), we instantaneously remove the blockage and run another 15-sec simulation on the original network. We observe a return to the previous equilibrium in a certain amount of time, which we call the recovery time. This is defined as the time that the system take to return at its initial state when the blockage is removed. In these simulations the mean pressure $AMP$ was computed using a physiologically-inspired formula: $AMP=(P_{h}+2P_{l})/3$, where $P_{h}$ and $P_{l}$ are the highest and the lowest value of the pressure, although the time averaged mean blood pressure can also be used, with similar results. For $HR=75$, note that the recovery time is slightly higher than 5 seconds.

The entire 40-sec sequence of the pressure and flow dynamics (for $HR=75$) in the middle of edges 1, 3 and 4 is presented in Figure 8. Note that during the blockage/unblockage sequence, pressure and flow behave differently in various parts of the network. Most notably, we see that the blockage happens during the diastolic period, and since it is instantaneous (like a clamping), all the fluid mass past edge 3 is lost and no backflow is generated, hence the upstream pressure (of edge 1) is elevated only for the next systolic period, after which the pressure settles around its new equilibrium value (in this case lower for edge 1). This decrease in steady state pressure values may be different than what is observed in the real vascular system. In that case, an amputation of the distal sub-network causes a significant increase in pressure in the upstream vessels [15].

In edge No. $3$ the systolic pressure increases during the entire blockage segment, while in other edges (e.g. No. $4$) there is a sustained decrease in systolic pressure. This is consistent with the fact that edge 3 is a terminal edge during the 5-sec blockage. Also worth noting is the fact that, for short period of times during the blockage, the pressure becomes negative. This simply means that the wall displacement becomes sufficiently large (in the negative -or inward -direction) for the force per unit area to change its direction from being outward (positive pressure) to inward (negative pressure). From a modeling perspective, this is simply a reflection of the non-collapsible nature of the tube walls, and the fact that longitudinal restoring forces can balance any inward (or outward) force. The original algebraic assumption of the pressure difference between hydrodynamic pressure $P$ and external pressure $P_{ext}$ being proportional to wall displacement can be made responsible for this behavior. Improved models that take into account the mechanical wall properties can better explain these negative values, but we defer this discussion elsewhere. The fact that the bulk of the network is now past edge 4 explains why the new equilibrium pressures in the new network are lower while the flow is higher. This suggests that the presence of significant bifurcations more distal in the network is creating less pressure buildup in the upstream vessels.

These 40-sec simulations were then employed to assess the recovery time as a function of the heart rate. An optimization using MATLAB's fminsearch yielded the optimal value of the heart rate parameter $HR=72.2$ beats/min, corresponding to a minimum recovery time of $3.5$ sec. Values of $HR$ below $70$ or above $75$ gave much longer recovery times. This result indicates that increasing $HR$ (hence cardiac output) does not necessarily translate to a quicker recovery time, which can be attributed to the reflections at junctions and at the terminal sites. This is also suggested by the the flows in edge $1$ & $3$, which show an increase in magnitude of the reflected waves.

---

5. Discussion

The results presented here are based on very restrictive assumptions on the wall properties, boundary conditions etc, and hence have limited range of applicability as far as quantitative analysis or patient specific applications are concerned. The intent here has been to lay the groundwork for the inclusion of the spatial features in the optimization of the arterial network and to apply such analyses when more realistic assumptions are made. In particular, the issue of imposing appropriate outflow boundary conditions remains crucial as described earlier. While inflow conditions can be harvested from MRI data, outflow conditions are difficult/impractical to obtain, due to the limited access to terminal sites of the spatial network. Hence the importance of using appropriate models for the peripheral circulation or micro-vasculature network. Such further studies will hopefully lead to novel nonlinear tools for assessing the autonomic control mechanism and ultimately to patient specific assessment tools which can be used in clinical setting.

Simulations of the mathematical models involving PDEs on networks such as those presented in this paper reveal macroscopic phenomena that cannot be anticipated by looking at individual edge dynamics. The nature of phenomena such as appearance of low frequency oscillations in the network, time of recovery after blockage in a network is removed, truncation of a fractal tree network are revealed through these simulations and can help further study the real phenomena observed in physiological conditions. Interpreting physically correct boundary conditions remains crucial for modeling for long term behavior of the network dynamics and this is where improvements of the models can and will be done in future studies. The pulsatile nature of the dynamics on networks and the particle-like behavior of the pulses can also lead to higher order models, such as nonlinear dispersive models posed on a network (tree).

---

Appendix A. Numerical discretization

The system (1) can be written in conservation form:

$\frac{\partial \mathbf{U}}{\partial t}+\frac{\partial \mathbf{F}}{\partial x} =\mathbf{S},$

(4)

with

$\mathbf{U}=\left[ \begin{array}{c} A \\ U \end{array} \right] ,\text{ }\mathbf{F(U)}=\left[ \begin{array}{c} AU \\ \frac{U^{2}}{2}+\frac{P}{\rho } \end{array} \right] \text{ and }\mathbf{S(U)}=\left[ \begin{array}{c} 0 \\ \frac{1}{\rho }\left( \frac{f}{A}-\frac{\partial P}{\partial \beta }\frac{ d\beta }{dx}-\frac{\partial P}{\partial A_{0}}\frac{dA_{0}}{dx}\right) \end{array} \right] .$

We use a discontinuous Galerkin scheme, which has been successfully applied in this context (see [3], [29]) to solve the system. This method has advantages over finite difference schemes, especially since the geometry of our spatial domain inherently leads to discontinuities at junctions. Moreover, discontinuous Galerkin methods are better suited to accommodate higher order terms, such as those present in the visco-elastic and inertial models.

We now describe the discretization scheme, which follows closely [29]. First we discretize the domain $\Omega =[0,l]$ into a mesh of $N_{el}$ elemental non-overlapping regions $\Omega _{e}=(x_{e}^{L},x_{e}^{R})$, such that $x_{e}^{R}=x_{e+1}^{L}$ for $e=1,...,N_{el}$ and

$\bigcup\limits_{e=1}^{N_{el}}\Omega _{e}=\Omega .$

The weak form of the system is obtained by multiplying the equation (4) by a vector of test functions $\mathbf{\Phi }$ and integrating over $\Omega$:

$\left( \frac{\partial \mathbf{U}}{\partial t},\mathbf{\Phi }\right) _{\Omega }+\left( \frac{\partial \mathbf{F}}{\partial x},\mathbf{\Phi }\right) _{\Omega }=(\mathbf{S,\Phi })_{\Omega },$

where

$(\mathbf{w,v})_{\Omega }=\int_{\Omega }\mathbf{wv}dx. \label{F3}$

(5)

The integrals are decomposed into elemental regions as follows:

$\sum\limits_{e=1}^{N_{el}}\left( \left( \frac{\partial \mathbf{U}}{\partial t}, \mathbf{\Phi }\right) _{\Omega _{e}}+\left( \frac{\partial \mathbf{F}}{ \partial x},\mathbf{\Phi }\right) _{\Omega _{e}}\right) =\sum\limits_{e=1}^{N_{el}}( \mathbf{S,\Phi })_{\Omega _{e}}, \label{F10}$

(6)

and the second member of (6) is integrated by parts:

$\sum\limits_{e=1}^{N_{el}}\left( \left( \frac{\partial \mathbf{U}}{\partial t}, \mathbf{\Phi }\right) _{\Omega _{e}}-\left( \mathbf{F},\frac{d\mathbf{\Phi } }{dx}\right) _{\Omega _{e}}+[\mathbf{F\cdot \Phi ]}_{x_{e}^{L}}^{x_{e}^{R}} \right) =\sum\limits_{e=1}^{N_{el}}(\mathbf{S,\Phi })_{\Omega _{e}}.$

The solution $\mathbf{U}(x,t)$ is approximated by a discretised expansion denoted by $\mathbf{U}^{\delta }(x,t)$ and, in the same way, $\mathbf{\Phi } (x)$ is approximated by $\mathbf{\Phi }^{\delta }(x).$ As basis for the expansion we have chosen polynomials of degree $K$ on each elemental region $\Omega _{e}$. In addition, to obtain a global solution in the domain $\Omega$, information must propagate between elemental regions $\Omega _{e}$ and this is achieved by upwinding the boundary flux, which is denoted as $\mathbf{F}^{u}$.

The upwinded fluxes on each side of the interface, $\mathbf{F}_{L}^{u}$ and $\mathbf{F}_{R}^{u}$, are calculated by solving the Riemann problem at each interface (see Section 3) and setting as $\mathbf{F}_{L}^{u}=\mathbf{F}(A_{L}^{u},U_{L}^{u})$ and $\mathbf{F}_{R}^{u}= \mathbf{F}(A_{R}^{u},U_{R}^{u})$.

In this way we obtain:

$\sum\limits_{e=1}^{N_{el}}\left( \left( \frac{\partial \mathbf{U}_{e}^{\delta }}{ \partial t},\mathbf{\Phi }_{e}^{\delta }\right) _{\Omega _{e}}-\left( \mathbf{F(\mathbf{U}_{e}^{\delta })},\frac{d\mathbf{\Phi }_{e}^{\delta }}{dx} \right) _{\Omega _{e}}+[\mathbf{F}^{u}\mathbf{\cdot \Phi }_{e}^{\delta } \mathbf{]}_{x_{e}^{L}}^{x_{e}^{R}}\right) =\sum\limits_{e=1}^{N_{el}}(\mathbf{S( \mathbf{U}_{e}^{\delta }),\Phi }_{e}^{\delta })_{\Omega _{e}}.$

Integrating again the second term by parts we get:

$\sum\limits_{e=1}^{N_{el}}\left( \left( \frac{\partial \mathbf{U}_{e}^{\delta }}{ \partial t},\mathbf{\Phi }_{e}^{\delta }\right) _{\Omega _{e}}+\left( \frac{ \partial \mathbf{F(\mathbf{U}_{e}^{\delta })}}{\partial x},\mathbf{\Phi } _{e}^{\delta }\right) _{\Omega _{e}}+[(\mathbf{F}^{u}-\mathbf{F(\mathbf{U} _{e}^{\delta }))\cdot \Phi }_{e}^{\delta }\mathbf{]}_{x_{e}^{L}}^{x_{e}^{R}} \right) = \notag \\ =\sum\limits_{e=1}^{N_{el}}(\mathbf{S(\mathbf{U}_{e}^{\delta }),\Phi }_{e}^{\delta })_{\Omega _{e}}.$

(7)

To simplify the method, we have mapped each elemental region onto the standard element $\Omega _{st}=\{\xi \in\mathbb{R} :-1\leq \xi \leq 1\}$. This mapping is defined as

$\chi _{e}(\xi )=x_{e}^{L}\frac{1-\xi }{2}+x_{e}^{R}\frac{1+\xi }{2}, \xi \in \Omega _{st},$

and its inverse is given by

$\xi =\chi _{e}^{-1}(x)=2\frac{x_{e}-x_{e}^{L}}{x_{e}^{R}-x_{e}^{L}}-1, x_{e}\in \Omega _{e}.$

We selected as expansion basis the Legendre polynomials $L_{k}(\xi )$, with $k$ the polynomial order, because they are orthogonal with respect to the product (5). In this way, the solution is expanded on each elemental region $\Omega _{e}$ as

$\mathbf{U}_{e}^{\delta }(\chi _{e}(\xi ),t)=\sum\limits_{k=0}^{K}L_{k}(\xi ) \widehat{\mathbf{U}_{e}^{k}}(t), \label{F5}$

(8)

with $\widehat{\mathbf{U}_{e}^{k}}(t)$ the time-varying coefficients of the expansion.

Replacing (8) in (7) and letting $\mathbf{\Phi }_{e}^{\delta }= \mathbf{U}_{e}^{\delta }$, we obtain $2K$ differential equations to be solved for each $\Omega _{e},e=1,...,N_{el}:$

$\frac{d\widehat{U_{i,e}^{k}}}{dt}=\mathcal{F}(\mathbf{U}_{e}^{\delta }), k=0,...,K,\text{ }i=1,2,$

where $\widehat{U_{i,e}^{k}},i=1,2,$ are each of the two components of $\widehat{\mathbf{U}_{e}^{k}}(t)$ and

$\mathcal{F}(\mathbf{U}_{e}^{\delta })=-\left( \frac{\partial F_{i}}{\partial x},L_{k}\right) _{\Omega _{e}}-\frac{2}{x_{e}^{R}-x_{e}^{L}}[L_{k}\ast (F_{i}^{u}-F_{i}(\mathbf{U}_{e}^{\delta }))]_{x_{e}^{L}}^{x_{e}^{R}}+(S_{i}( \mathbf{U}_{e}^{\delta }),L_{k})_{\Omega _{e}}.$

The method is completed with a second-order Adams-Bashforth time-integration scheme:

$\left( \widehat{U_{i,e}^{k}}\right) ^{n+1}=\left( \widehat{U_{i,e}^{k}} \right) ^{n}+\frac{3\Delta t}{2}\mathcal{F}\left( (\mathbf{U}_{e}^{\delta })^{n}\right) -\frac{\Delta t}{2}\mathcal{F}\left( (\mathbf{U}_{e}^{\delta })^{n-1}\right) , \\ k=0,...,K,\text{ }i=1,2,\text{ }e=1,...,N_{el},$

in which $\Delta t$ is the time step and $n$ the number of every time step. To calculate the integrals we use a Gauss quadrature formula of order $q\geq K+1$.

---

Acknowledgments

We thank the anonymous reviewers for very substantive recommendation which helped improved the presentation. R.C. was supported by a grant from the UCCS BioFrontiers Center of the University of Colorado. R.C. is also grateful for the hospitality of the Universitá degli Studi di Salerno during his stay there.

---