Faisal Amlani

Niema M Pahlevan

A stable high-order FC-based methodology for hemodynamic wave propagation

Keywords: Fourier continuation, hemodynamic wave simulation, mathematical physiology, high-order methods

A high-order numerical algorithm is proposed for the solution of one-dimensional arterial pulse wave propagation problems based on use of an accelerated "Fourier continuation" (FC) methodology for accurate Fourier expansion of non-periodic functions. The solver provides high-order accuracy, mild Courant-Friedrichs-Lewy (CFL) constraints on the time discretization and, importantly, results that are essentially free of spatial dispersion errors-enabling fast and accurate resolution of clinically-relevant problems requiring simulation of many cardiac cycles or vascular segments. The left ventricle-arterial model that is employed presents a particularly challenging case of ordinary differential equation (ODE)-governed boundary conditions that include a hybrid ODE-Dirichlet model for the left ventricle and an ODE-based Windkessel model for truncated vasculature. Results from FC-based simulations are shown to capture the important physiological features of pressure and flow waveforms in the systemic circulation. The robustness of the proposed solver is demonstrated through a number of numerical examples that include performance studies and a physiologically-accurate case study of the coupled left ventricle-arterial system. The results of this paper imply that the FC-based methodology is straightforwardly applicable to other biological and physical phenomena that are governed by similar hyperbolic partial differential equations (PDEs) and ODE-based time-dependent boundaries.

Introduction

Fourier continuation (FC) methods [START_REF] Bruno | High-order unconditionally stable FC-AD solvers for general smooth domains I. Basic elements[END_REF][START_REF] Lyon | High-order unconditionally stable FC-AD solvers for general smooth domains II. Elliptic, parabolic and hyperbolic PDEs; theoretical considerations[END_REF][START_REF] Albin | A spectral FC solver for the compressible Navier-Stokes equations in general domains I: Explicit time-stepping[END_REF][START_REF] Bruno | Spatially dispersionless, unconditionally stable FC-AD solvers for variablecoefficient PDEs[END_REF][START_REF] Amlani | An FC-based spectral solver for elastodynamic problems in general threedimensional domains[END_REF][START_REF] Bruno | Higher-order implicit-explicit multi-domain compressible Navier-Stokes solvers[END_REF] broaden the applicability of Fourier-based partial differential equation (PDE) solvers by resolving the well-known Gibbs "ringing" phenomenon to enable high-order convergence of Fourier series approximations of non-periodic functions. This is accomplished by constructing an interpolating Fourier series representation via a discrete periodic extension that approximates a given function to high-order and renders it periodic on a slightly larger domain. FC-based PDE solvers have already been introduced for a number of problems including the classical wave and diffusion equations with constant and variable coefficients [START_REF] Bruno | High-order unconditionally stable FC-AD solvers for general smooth domains I. Basic elements[END_REF][START_REF] Lyon | High-order unconditionally stable FC-AD solvers for general smooth domains II. Elliptic, parabolic and hyperbolic PDEs; theoretical considerations[END_REF][START_REF] Bruno | Spatially dispersionless, unconditionally stable FC-AD solvers for variablecoefficient PDEs[END_REF][START_REF] Elling | GPU-accelerated Fourier-continuation solvers and physically exact computational boundary conditions for wave scattering problems[END_REF]; the nonlinear Burgers system [START_REF] Bruno | Higher-order linear-time unconditionally stable alternating direction implicit methods for nonlinear convection-diffusion partial differential equation systems[END_REF]; the Euler equations [START_REF] Shahbazi | Multi-domain Fourier-continuation/WENO hybrid solver for conservation laws[END_REF][START_REF] Shahbazi | Multi-dimensional hybrid Fourier continuation -WENO solvers for conservation laws[END_REF]; the 2D compressible Navier-Stokes equations [START_REF] Albin | A spectral FC solver for the compressible Navier-Stokes equations in general domains I: Explicit time-stepping[END_REF][START_REF] Albin | Fourier continuation methods for high-fidelity simulation of nonlinear acoustic beams[END_REF][START_REF] Bruno | Higher-order implicit-explicit multi-domain compressible Navier-Stokes solvers[END_REF]; and the 3D elasticity equations with variable coefficients [START_REF] Amlani | An FC-based spectral solver for elastodynamic problems in general threedimensional domains[END_REF]. These solvers maintain the well-established qualities of other spectral solvers for time-domain equations including accuracy by means of coarse discretizations and a faithful preservation of dispersion characteristics of the underlying continuous problems. Additionally, although classical spectral solvers can include restrictive Courant-Friedrichs-Lewy (CFL) conditions on time discretizations as well as restrictions to periodic geometries, explicit FC-based solvers-like the one presented in this paper-rely on equispaced Cartesian grids (which can be extended to curvilinear formulations [START_REF] Albin | A spectral FC solver for the compressible Navier-Stokes equations in general domains I: Explicit time-stepping[END_REF][START_REF] Amlani | An FC-based spectral solver for elastodynamic problems in general threedimensional domains[END_REF]) that lead to mild, linear CFL constraints on the timestep while treating physically-realistic non-periodic domains.

The objective of this work is to introduce a new stable high-order algorithm for the numerical solution of one-dimensional (1D) arterial pulse wave problems based on this Fast Fourier Transform (FFT)-speed FC methodology. In particular, the physical model employed represents the most complicated case found in wave propagation models of the cardiovascular system: an arterial segment governed by a hyperbolic PDE whose left boundary is governed by a complex hybrid ordinary differential equation (ODE)-Dirichlet condition for the left ventricle (LV) and whose right boundary is governed by a Windkessel ODE. The former boundary condition, in particular, is "non-stationary" in the sense that the time at which the ODE-governed condition transitions to a Dirichlet condition is itself determined by the solution of the PDE. This complex nonlinear system is representative of the myriad formulations of governing equations that have been constructed to reasonably capture the complete physics of pressure and flow propagation in the cardiovascular system; hence the numerical methodology herein can be easily extended to other cases including those utilizing simpler boundary conditions (e.g., a given inflow function, junction-matching conditions or one-sided lumped terminal models). The method described in this paper enables fast, accurate and essentially dispersionless resolution of arterial wave propagation systems, ultimately yielding errors within a prescribed tolerance using a fixed number of points per wavelength for increasingly large problems and yielding stable simulation for metrics that require information over a large number of cardiac cycles. Such a fast and dispersionless method may also be suitable for feedback regulatory models necessitating long-time simulations as well as for iteratively solving inverse problems related to arterial wave reflections.

One-dimensional arterial models are noted for their relatively inexpensive computational complexity when compared to 3D models. Several numerical methodologies have been constructed for the various formulations (e.g., area-velocity, area-flow, pressure-flow). These include characteristic methods [START_REF] Wang | Wave propagation in a model of the arterial circulation[END_REF] as well as finite element (FE) [START_REF] Formaggia | One-dimensional models for blood flow in arteries[END_REF][START_REF] Matthys | Pulse wave propagation in a model human arterial network: assessment of 1-D numerical simulations against in vitro measurements[END_REF][START_REF] Steele | In vivo validation of a one-dimensional finiteelement method for predicting blood flow in cardiovascular bypass grafts[END_REF][START_REF] Mynard | A 1D arterial blood flow model incorporating ventricular pressure, aortic valve and regional coronary flow using the locally conservative galerkin (LCG) method[END_REF][START_REF] Sherwin | One-dimensional modelling of a vascular network in spacetime variables[END_REF], finite volume (FV) [START_REF] Müller | Well-balanced high-order solver for blood flow in networks of vessels with variable properties[END_REF] and finite difference (FD) [START_REF] Berger | Differential effects of wave reflections and peripheral resistance on aortic blood pressure: a model-based study[END_REF][START_REF] Olufsen | Numerical simulation and experimental validation of blood flow in arteries with structured-tree outflow conditions[END_REF][START_REF] Azer | A one-dimensional model of blood flow in arteries with friction and convection based on the womersley velocity profile[END_REF][START_REF] Saito | One-dimensional model for propagation of a pressure wave in a model of the human arterial network: Comparison of theoretical and experimental results[END_REF][START_REF] Du | A fast algorithm for the simulation of arterial pulse waves[END_REF] methods. A number of these tools have been successfully utilized for blood flow analysis and are wellestablished. Some of these methods, however, can have restrictive spectral radii resulting in very fine timesteps relative to spatial mesh sizes. Additionally, FE/FV/FD-based methods have been well-known for a long time to suffer from high numerical dispersion [START_REF] Mullen | Dispersion analysis of finite element semidiscretizations of the two-dimensional wave equation[END_REF]: phase errors of a solution can accumulate over subsequent periods of a propagating wavetrain, requiring an ever-increasing numbers of points per wavelength to resolve larger problems within a given accuracy. For excellent reviews of the derivations and implementation of both analytical and numerical models used in 1D arterial propagation, the reader is referred to [START_REF] Boileau | A benchmark study of numerical schemes for one-dimensional arterial blood flow modelling[END_REF][START_REF] Shi | Review of zero-D and 1-D models of blood flow in the cardiovascular system[END_REF][START_REF] Wang | Verification and comparison of four numerical schemes for a 1D viscoelastic blood flow model[END_REF].

For hemodynamics applications, fast and accurate 1D solvers are necessary to explore all possible dynamics that result from changes brought about by diseases. Stable and high-order resolution can extend the applicability of reduced-order systems towards more sophisticated models (e.g., that of this paper) and towards multi-dimensional coupling (e.g., 1D-3D coupling) in an ultimate effort to improve physiological relevancy and utility. Robust methods are also essential for maintaining cost reductions promised by 1D formulations; large numbers of simulations may be required to study, for example, alterations in geometric and elastic properties of the arterial tree, or PDE-constrained optimization problems [START_REF] Pahlevan | Intrinsic frequency for a systems approach to haemodynamic waveform analysis with clinical applications[END_REF]. To this end, the FC-based methodology proposed in what follows is a competitive approach for forward problems. In Section 2, the full governing equations are presented for the LV-arterial system. Section 3 describes the full solver including a detailed presentation of the construction of Fourier continuations as well as a full summary of the overall algorithm and implementation of the complex boundary ODEs. Section 4 provides a variety of illustrative numerical experiments attesting to high-order convergence, limited numerical dispersion and long-time stability. Section 5 presents a case study that simulates waveforms from physiologically-accurate data by an LV-arterial system model; comparisons with a finite difference solver are also provided. Concluding remarks are ultimately discussed in Section 6.

The theoretical model

One-dimensional arterial pulse wave models have been demonstrated to capture important features of pressure, flow and area wave propagation in large arteries, and have been demonstrated to illustrate many features of circulation under both normal and diseased conditions. The problem of propagation of pressure and flow waves in an arterial system is often modeled by decomposing the domain into segments of elastic tubes of given lengths whose properties can be described by single axial coordinates [START_REF] Van De Vosse | Pulse wave propagation in the arterial tree[END_REF][START_REF] Berger | Differential effects of wave reflections and peripheral resistance on aortic blood pressure: a model-based study[END_REF][START_REF] Azer | A one-dimensional model of blood flow in arteries with friction and convection based on the womersley velocity profile[END_REF][START_REF] Hughes | On the one-dimensional theory of blood flow in the larger vessels[END_REF][START_REF] Berger | Wave propagation in coupled left ventricle-arterial system[END_REF]. A number of equivalent formulations exist [START_REF] Boileau | A benchmark study of numerical schemes for one-dimensional arterial blood flow modelling[END_REF][START_REF] Shi | Review of zero-D and 1-D models of blood flow in the cardiovascular system[END_REF][START_REF] Wang | Verification and comparison of four numerical schemes for a 1D viscoelastic blood flow model[END_REF]: those in terms of cross-sectional area of the aorta A(x, t), pressure P (x, t) and axial velocities U (x, t) (the so-called AU formulations); those in terms of the area and the flow rate Q(x, t) = AU (the so-called AQ formulations); as well as formulations in terms of Z 0 (aorta)

Z 0 C R s load impedance (arterial system)
left ventricle

(with aortic valve)

Figure 1: The coupled LV-arterial model consisting of a linear transmission tube of length (aorta); the input boundary condition (left ventricle) modeled as a time-varying compliance; and the outflow boundary modeled as a terminal load impedance represented by a three-element Windkessel [START_REF] Berger | Differential effects of wave reflections and peripheral resistance on aortic blood pressure: a model-based study[END_REF]. The parameters Z 0 , Cc and Rc define the characteristic impedance, the compliance and the resistance of the eliminated vascular network, respectively.

pressure and flow rates (the so-called P Q formulations). All of these versions are constructed to account for the fluid-structure interaction of the problem and can be solved with similar numerical strategies. This paper is concerned with a coupled P Q formulation for purposes of modeling the LV-arterial system illustrated by Figure 1. This system represents a most complicated configuration found in modeling segments in an arterial system (a hyperbolic PDE coupled to a complex hybrid ODE-Dirichlet boundary condition for the left ventricle and a Windkessel ODE for truncated vasculature). This simplified LV-arterial model has been shown to capture physiologically-relevant characteristics of the systemic circulation [START_REF] Berger | Differential effects of wave reflections and peripheral resistance on aortic blood pressure: a model-based study[END_REF]. The numerical methodology subsequently described in Section 3 can be also applied to other formulations.

Governing equations

In the P Q formulation [START_REF] Van De Vosse | Pulse wave propagation in the arterial tree[END_REF][START_REF] Berger | Differential effects of wave reflections and peripheral resistance on aortic blood pressure: a model-based study[END_REF][START_REF] Berger | Wave propagation in coupled left ventricle-arterial system[END_REF][START_REF] Campbell | Pulse reflection sites and effective length of the arterial system[END_REF], the pressure P (x, t) and flow Q(x, t) in a single uniform elastic tube of a given length (e.g., the aorta in Figure 1) can be governed by

  ∂P ∂t (x, t) ∂Q ∂t (x, t)   = -   1 C ∂Q ∂x (x, t) 1 L ∂P ∂x (x, t) + R L Q(x, t)   , 0 ≤ x ≤ , t ≥ 0, (1)
with initial conditions P (x, t = 0) and Q(x, t = 0) prescribed at some initial time t = 0 (and evolved for a given number of cardiac cycles of period T per ). The system (1) derives from the assumptions that 1) blood fluid is Newtonian (described by linearized Navier-Stokes equations); 2) that velocity fields are axisymmetric; and 3) that the tube walls displace elastically. Material parameters are specified by the fluid inductance L ∈ R + , the hydraulic resistance R ∈ R + and the vessel compliance C ∈ R + . The pulse wave velocity (needed for the CFL condition of the explicit solver presented in Section 3) can be approximated as

c PW = 1 √ LC . (2)

Inlet boundary conditions

The inlet (left) boundary condition at x = 0 represents the dynamics of the left ventricle in Figure 1 and can be accurately modeled as a time varying compliance C v (t) [START_REF] Berger | Differential effects of wave reflections and peripheral resistance on aortic blood pressure: a model-based study[END_REF][START_REF] Senzaki | Single-beat estimation of end-systolic pressure-volume relation in humans: A new method with the potential for noninvasive application[END_REF][START_REF] Kass | Ventriculo-arterial coupling: Concepts, assumptions, and applications[END_REF] defining a relationship between the pressure P v (t) inside the left ventricle and the corresponding LV-volume

V v ∈ R + as V v (t) -V dead = C v (t)P v (t). (3)
The constant V dead , known as the dead volume, defines the limit for pressure generation: a ventricle whose volume falls below this value can no longer generate any pressure [START_REF] Suga | Instantaneous pressure-volume relationships and their ratio in the excised, supported canine left ventricle[END_REF]. Recalling that the flow into the aorta is defined in terms of this volume as Q(x = 0, t) = -∂V v /∂t, differentiating equation (3) with respect to t yields the corresponding differential equation for the pressure in the ventricle as Hence once the pressure P v (t) in the ventricle is greater than that of the tube boundary P (x = 0, t), the valve opens and P (x = 0, t) = P v (t) with the corresponding flow condition given naturally by the governing arterial wave equation [START_REF] Bruno | High-order unconditionally stable FC-AD solvers for general smooth domains I. Basic elements[END_REF]. Once the inflow reaches zero (or, numerically, the time at which Q(x = 0, t) ≤ 0 or P v (t) ≤ P (x = 0, t)), the valve closes and the left boundary condition remains Q(x = 0, t) = 0 (a Dirichlettype condition). Generally, C v (t) is given by clinical parameters either through a closed-form approximation to data or a look-up table [START_REF] Berger | Differential effects of wave reflections and peripheral resistance on aortic blood pressure: a model-based study[END_REF][START_REF] Senzaki | Single-beat estimation of end-systolic pressure-volume relation in humans: A new method with the potential for noninvasive application[END_REF][START_REF] Kass | Ventriculo-arterial coupling: Concepts, assumptions, and applications[END_REF].

∂P v ∂t (t) = - 1 C v (t) ∂C v ∂t (t)P v (t) + Q(x = 0, t) , t ≥ 0. (4)

Outlet boundary conditions

At the "terminal boundary" x = (right boundary in Figure 1), the physical structure approximates the effect of the eliminated peripheral vessels. This effect can be reasonably modeled using a zero-dimensional lumped parameter model, known as a Windkessel, that couples to the vessel governed by equation ( 1) through a matched impedance that is related to the vessel parameters by [START_REF] Berger | Differential effects of wave reflections and peripheral resistance on aortic blood pressure: a model-based study[END_REF]

Z 0 = L C , (5)
for fluid inductance L and vessel compliance C. This circuit-like Windkessel model relates the pressure and the flow at the outflow boundary. The pressure P c in the terminal compliance chamber is related to the aortic pressure P (x = , t) at the outlet of the vessel via the ODE given by

∂P c ∂t (t) = 1 C c Z 0 P (x = , t) - R c + Z 0 C c R c Z 0 P c (t) (6)
for effective compliance C c ∈ R + , characteristic impedance Z 0 ∈ R + (approximated by equation ( 5)) and a total peripheral resistance R c ∈ R + -all chosen so that the load is matched to the domain (tube) at high frequencies [START_REF] Berger | Differential effects of wave reflections and peripheral resistance on aortic blood pressure: a model-based study[END_REF]. The corresponding outflow Q(x = , t) at the terminal boundary is given by

Q(x = , t) = 1 Z 0 (P (x = , t) -P c (t)) . (7)
The solver implementation of this complex boundary condition-and that of the inlet dictacted by equation (4)-is described in detail in Section 3.5.

A summary of the physical parameters

The physical parameters (often given by clinical data) that are necessary for the complete modelincluding the boundary ODEs-are summarized in Table 1.

The numerical model

This section describes the proposed FC-based arterial pulse wave solver, including the construction of continuations and their subsequent spatial differentiation; a corresponding explicit treatment for temporal derivatives; and a spectral filter design that eleminates growth of high-frequency errors without unduly deteriorating accuracy. The general principles of Fourier continuation are presented in Section 3.1 and a relatively detailed summary for constructing (accelerated) Fourier continuations is presented in Section 3.2. The implementation of the complex boundary conditions for the inlet and outlet is also presented (Section 3.5) along with a pseudo-code of the full algorithm (Section 3.6).

The Fourier continuation approach

Considering a uniform discretization

x i = i∆x, i = 0, . . . , N -1, ∆x = 1/(N -1) and corresponding point values f (x i ) of a given smooth function f (x) : [0, 1] → R, the FC method constructs a rapidly-convergent interpolating Fourier series representation f cont : [0, b] → R on a region [0, b]
that is slightly larger than the given physical domain [0, 1]. This representation is given by

f cont = M k=-M a k e 2πikx b s.t. f cont (x i ) = f (x i ), i = 0, ..., N -1, (8)
where M = (N + N cont )/2 is the bandwidth for an N cont (defined as the number of points added to the original domain such that b = (N + N cont )∆x). The Fourier continuation function f cont approximately renders the original function f periodic, i.e., f cont approximates f to high-order in the original domain [0, 1] and is periodic on the slightly larger domain [0, b], b > 1. Spatial derivatives in the governing equations of (1) can then be produced by exact termwise differentiation of the series (8) as

∂f cont ∂x (x) = M k=-M 2πik b a k e 2πikx b . (9)
This ultimately provides the numerical derivatives of f to high-order by restricting the domain of ∂f cont /∂x to the original unit interval. Hence the approximation rests in the construction of (8) from which the computation of the derivative (9) can be facilitated by the Fast Fourier Transform (FFT). Note here that f has been defined on [0, 1] without loss of generality: it is easily extended to any interval by affine transformations [START_REF] Bruno | High-order unconditionally stable FC-AD solvers for general smooth domains I. Basic elements[END_REF].

Accelerated Fourier continuation: FC(Gram)

The coefficients a k of (8) are found in the simplest treatment [START_REF] Bruno | Accurate, high-order representation of complex three-dimensional surfaces via Fourier continuation analysis[END_REF][START_REF] Boyd | Exponentially-convergent strategies for defeating the Runge phenomenon for the approximation of non-periodic functions, part I: single-interval schemes[END_REF] through solving the least squares system min

a k N -1 i=0 |f cont (x i ) -f (x i )| 2 (10)
by the Singular Value Decomposition (SVD). This can become rather costly for time-dependent solutions of complex boundary-valued PDEs such as equation (1) (each timestep requires application of the corresponding SVDs). In order to reduce this computational expense, accelerated techniques have been introduced [START_REF] Bruno | High-order unconditionally stable FC-AD solvers for general smooth domains I. Basic elements[END_REF][START_REF] Lyon | High-order unconditionally stable FC-AD solvers for general smooth domains II. Elliptic, parabolic and hyperbolic PDEs; theoretical considerations[END_REF][START_REF] Albin | A spectral FC solver for the compressible Navier-Stokes equations in general domains I: Explicit time-stepping[END_REF]. Known as FC(Gram), they rely on use of small vectors of only a handful of function values near the left and right endpoints at x = 0 and x = 1 which are projected onto a Gram polynomial basis whose continuations are precomputed through solving the corresponding least squares problem [START_REF] Shahbazi | Multi-dimensional hybrid Fourier continuation -WENO solvers for conservation laws[END_REF] by means of a high-precision SVD. More precisely, the accelerated FC(Gram) method uses a subset of the given function values on small numbers d and d r of matching points {x 0 , .., x d -1 } and {x N -dr , ..., x N -1 } contained in small subintervals on the left and right ends of the interval [0, 1] to produce a discrete periodic extension of size N cont . This is accomplished by projecting these end values onto a Gram basis up to degree d -1 (resp. d r -1) of polynomials (producing a polynomial interpolant) whose FC extensions are precomputed and whose orthogonality is enforced by the natural discrete scalar product defined by the discretization points. This effectively forms a "basis" of continuation functions with which to rapidly extend the given function f and provide a smooth transition from f (x = 0) to f (x = 1) over the interval [0, b].

In particular, this work employs a "blend-to-zero" version of FC(Gram) [START_REF] Albin | A spectral FC solver for the compressible Navier-Stokes equations in general domains I: Explicit time-stepping[END_REF][START_REF] Amlani | An FC-based spectral solver for elastodynamic problems in general threedimensional domains[END_REF] in order to obtain the additional N cont values. The strategy takes the function values defined at the matching points on the left (resp. the right) and provides a smooth transition to the left (resp. the right) that terminates in zero values at the the discrete points {x -Ncont-dr , ..., x -Ncont-1 } (resp. {x N +Ncont , ..., x N +Ncont+d -1 }). The sum of the leftward and rightward extensions (to zero) provides the needed complete Fourier continuation. Each blendto-zero transition is precomputed on a Gram polynomial basis (to interpolate the left d and right d r values of f ) by solving the corresponding minimization problems via SVD. This basis can be constructed for the left matching problem, for example, by orthogonalizing the Vandermonde matrix

M =      1 x 0 (x 0 ) 2 ... (x 0 ) d -1 1 x 1 (x 1 ) 2 ... (x 1 ) d -1 . . . . . . . . . . . . . . . 1 x d-1 (x d-1 ) 2 ... (x d -1 ) d -1      (11)
of point values at the discrete points {x 0 , ..., x d -1 } of the monomials {x 0 , x 1 , .., x d -1 } (the right-matching problem uses a substitution of 1 -x). Orthogonalization is performed via a stabilized Gram-Schmidt orthogonalization (M = QR) in high-precision arithmetic and with the monomial basis oversampled (see [START_REF] Albin | A spectral FC solver for the compressible Navier-Stokes equations in general domains I: Explicit time-stepping[END_REF][START_REF] Amlani | An FC-based spectral solver for elastodynamic problems in general threedimensional domains[END_REF] for details). This ultimately produces the left and right projection operators Q and Q r to determine the coefficients in the polynomial basis used in interpolating the left and right matching values of the function f . Each basis has its corresponding blend-to-zero continuation precomputed. Defining the vector of matching points for the left and right as

f = (f (x 0 ), f (x 1 ), ..., f (x d -1 )) T , f r = (f (x N -dr ), f (x N -dr+1 ), ..., f (x N -1 )) T , (12)
the continuation operation can hence be expressed in a block matrix form as

f cont = f A Q T f + A r Q T r f r , (13)
where f = (f (x 0 ), . . . , f (x N -1 )) T is a column vector containing the discrete point values of f ; f cont is a vector of the N + N cont continued function values; I is the N × N identity matrix; and A , A r contain the corresponding N cont values that blend the left and the right bases to zero. Again, the columns of Q , Q r contain the d , d r point values of each element of the corresponding Gram polynomial basis and, for d = d r , differ only by column-ordering (and A , A r by row-ordering). As in other FC-based solvers, a number of d , d r = 5 matching points are used for all simulations in this paper with a periodic extension comprised of N cont = 25 points. The matrices A , A r , Q and Q r need to be computed only once and stored in file for use each time the Fourier continuation of a given function is constructed. Further technical details on the construction of the blend-to-zero continuations employed in this paper can be found in [START_REF] Bruno | High-order unconditionally stable FC-AD solvers for general smooth domains I. Basic elements[END_REF][START_REF] Albin | A spectral FC solver for the compressible Navier-Stokes equations in general domains I: Explicit time-stepping[END_REF][START_REF] Amlani | An FC-based spectral solver for elastodynamic problems in general threedimensional domains[END_REF]. In summary, the FC(Gram) algorithm appends N cont values to a given discretized function in order to form a periodic extension in [1, b] that transitions smoothly from f (x N -1 ) back to f (x 0 ). Figure 2 illustrates an example Fourier continuation of a non-periodic function: the original discretized function on [0, 1] is translated by a distance N cont ∆x with the subsequent interval filled-in by the sum of leftward and rightward blend-to-zero continuations that ultimately renders the function periodic. Hence the resulting continued vector f cont can be seen as a set of discrete values of a smooth and periodic function that can be approximated to high-order via FFT on an interval of size (N + N cont )∆x.

Remark 3.1. Since the computation of derivatives for continued functions are spectrally accurate, the error of a Fourier continuation-based solver is dominated by the polynomial approximation used to project the end function values onto the blend-to-zero basis. Hence for the choices of d , d r = 5 matching points that are used throughout this paper, the corresponding fourth-order interpolating Gram polynomial basis results in essentially fifth-order convergence in space.

Explicit integration in time

For marching forward-in-time, the explicit Adams-Bashforth scheme of order four (AB4) is employed similarly to other FC-based solvers [START_REF] Albin | A spectral FC solver for the compressible Navier-Stokes equations in general domains I: Explicit time-stepping[END_REF][START_REF] Amlani | An FC-based spectral solver for elastodynamic problems in general threedimensional domains[END_REF]. Other explicit timesteppers can be used, including the fourth-order Runge-Kutta (RK4) method. Both methods provide adequate regions of absolute stability [START_REF] Bashforth | An attempt to test the theories of capillary action: by comparing the theoretical and measured forms of drops of fluid[END_REF][START_REF] Saidu | A simplified derivation and analysis of fourth order Runge-Kutta method[END_REF], but each timestep for RK4 entails four evaluations of the right-hand-side. For 1D problems, this is not particularly burdensome; however, enforcement of boundary conditions at intermediate RK steps may be problematic, especially for time-dependent boundary conditions [START_REF] Abarbanel | On the removal of boundary errors caused by Runge -Kutta integration of nonlinear partial differential equations[END_REF][START_REF] Carpenter | The theoretical accuracy of Runge-Kutta time discretizations for the initial boundary value problem: a study of the boundary error[END_REF] which are often utilized for hemodynamics models.

The system defined by the governing PDE in equation ( 1) and the governing boundary ODEs of equations ( 4) and ( 6) can be written as

∂ ∂t         P v (t) P (x, t) Q(x, t) P c (t)         = -         1 Cv(t) ∂Cv ∂t (t)P v (t) + Q(x = 0, t) 1 C ∂Q ∂x (x, t) 1 L ∂P ∂x (x, t) + R L Q(x, t) -1 CcZ0 P (x = , t) + Rc+Z0 CcRcZ0 P c (t)         =         F 1 (P v , Q, t) F 2 (P, Q, x, t) F 3 (P, Q, x, t) F 4 (P c , P, t)         . ( 14
)
Integration in time is performed on the basis of the explicit AB4 discretization with uniform timestep ∆t > 0.

The resulting discrete equations for a time t + ∆t are given by

        Pv(t + ∆t) P (x, t + ∆t) Q(x, t + ∆t) Pc(t + ∆t)         =         Pv(t) P (x, t) Q(x, t) Pc(t)         + ∆t 24         55         F1(t) F2(t) F3(t) F4(t)         -59         F1(t -∆t) F2(t -∆t) F3(t -∆t) F4(t -∆t)         + 37         F1(t -2∆t) F2(t -2∆t) F3(t -2∆t) F4(t -2∆t)         -9         F1(t -3∆t) F2(t -3∆t) F3(t -3∆t) F4(t -3∆t)                 , (15)
with the notational license F 1 (P v , Q, t) = F 1 (t), etc., for clarity. The boundary conditions calculated on the basis of the ODEs for P v , P c are then injected after each timestep t → t + ∆t; the resulting boundary values at time t + ∆t are then used to produce the FC approximations needed to evaluate the right hand sides at the next iteration for time t + ∆t. For all simulations in this paper, a timestep of

∆t ≤ √ 2 2 1 c PW ∆x ( 16
)
has provided absolute stability in all cases, including runs of hundreds of thousands of timesteps.

Filtering in frequency-space

Similarly to other spectral solvers, a spatial filter [START_REF] Albin | A spectral FC solver for the compressible Navier-Stokes equations in general domains I: Explicit time-stepping[END_REF][START_REF] Amlani | An FC-based spectral solver for elastodynamic problems in general threedimensional domains[END_REF] is employed to control the growth of the error in unresolved modes. This is applied in Fourier space on a function f (x) with Fourier coefficients fk by

N 2 k=-N 2 fk exp(ikx) -→ N 2 k=-N 2 exp (-α (2k/N ) p ) fk exp(ikx), ( 17
)
where p is a positive integer that determines the rate of decay of the coefficients, and where the real valued α > 0 determines the level of suppression such that highest-frequency modes are multiplied by e -α . For all numerical examples that follow, a mild choice of parameters inspired by [START_REF] Amlani | An FC-based spectral solver for elastodynamic problems in general threedimensional domains[END_REF][START_REF] Elling | GPU-accelerated Fourier-continuation solvers and physically exact computational boundary conditions for wave scattering problems[END_REF] is used and is given by (p, α) = (8, -Dc PW ∆t/∆x ln(10 -2 )), [START_REF] Müller | Well-balanced high-order solver for blood flow in networks of vessels with variable properties[END_REF] where c PW = 1/( √ LC) is the approximate pulse wave velocity and where ∆x is the spatial step size used in the computational domain. In general, a choice of D ≈ 1 provides absolute stability, including for those boundary conditions employed in Section 5 that result in possibly discontinuous-in-time solutions. Larger (stronger) values of α in equation [START_REF] Müller | Well-balanced high-order solver for blood flow in networks of vessels with variable properties[END_REF] have been demonstrated to still maintain the convergence and dispersion characteristics promised by FC-based solvers (e.g., the choice of α = 16 log(10) in [START_REF] Albin | A spectral FC solver for the compressible Navier-Stokes equations in general domains I: Explicit time-stepping[END_REF] eliminates the highest modes completely). However, the milder, CFL-based value chosen in [START_REF] Müller | Well-balanced high-order solver for blood flow in networks of vessels with variable properties[END_REF] ensures that the highest frequency terms aren't entirely eliminated and that the filter approaches unity as ∆t → 0 for a fixed ∆x (see [START_REF] Amlani | An FC-based spectral solver for elastodynamic problems in general threedimensional domains[END_REF][START_REF] Elling | GPU-accelerated Fourier-continuation solvers and physically exact computational boundary conditions for wave scattering problems[END_REF] for comparisons and more detailed discussions).

Implementation of the complex boundary conditions

The complex boundary conditions described in Sections 2.2 and 2.3 are governed by an ODE for the left ventricle (and at times a Dirichlet condition) and an ODE for the eliminated vasculature, respectively. At the inlet in particular, the ODE solution is only used when the valve is open, i.e., the pressure inside the ventricle is larger than the pressure at the inlet of the aorta. When the valve is closed, the boundary condition is simply given as a null inflow Q(x = 0, t) = 0. The corresponding valve condition must be tracked by a binary variable indicating whether the valve is open or closed. A summary of the implementation for the boundary updates governed by the inlet and outlet ODEs is given by the diagram of Figure 3.

Solver summary

The complete solver for the LV-arterial system is summarized in the pseudo-code of Algorithm 1.

Numerical performance: convergence, accuracy, dispersion and stability

Convergence and accuracy

Similarly to other verification procedures that have been used in literature (e.g., [START_REF] Raghu | Verification of a one-dimensional finite element method for modeling blood flow in the cardiovascular system incorporating a viscoelastic wall model[END_REF][START_REF] Raghu | Comparative study of viscoelastic arterial wall models in nonlinear one-dimensional finite element simulations of blood flow[END_REF] for hemodynamics and [START_REF] Amlani | An FC-based spectral solver for elastodynamic problems in general threedimensional domains[END_REF][START_REF] Roache | Code verification by the method of manufactured solutions[END_REF][START_REF] Roy | Verification of a compressible CFD code using the method of manufactured solutions[END_REF][START_REF] Vedovoto | Application of the method of manufactured solutions to the verification of a pressure-based finite-volume numerical scheme[END_REF] for higher-dimensional PDE systems), the method of manufactured solutions (MMS) can be employed to verify the implementation and numerical accuracy of the proposed FC-based methodology. In MMS, one postulates a smooth solution and algebraically derives the corresponding right-hand forcing terms and boundary conditions to render the proposed function an exact solution of the PDE. For the following experiment, solutions to the pressure P (x, t) and flow Q(x, t) are postulated as

P Q = P 0 sin(k 1 x -ω 1 t) cos(k 2 x -ω 2 t) Q 0 sin(k 2 x -ω 2 t) cos(k 1 x -ω 1 t) , (19)
for P 0 = 1/2, Q 0 = 3/2, k 1 = 2.5, k 2 = 6.5 and temporal frequencies ω 1 = 420Hz and ω 2 = -ω 1 . The domain considered is an elastic tube of length = 21cm with fluid parameters C = 0.0221ml/(mmHg•cm), L = 0.000256mmHgs 2 /cm 4 and R = 0.000024mmHg s/cm 2 . Figure 4 presents a snapshot of the solutions for both pressure P (x, t) and flow Q(x, t) at an arbitrarily chosen time. Using discretizations sizes that are integer multiples of the coarsest one used (N = 128 points), the simulation is advanced for 50 000 timesteps in all cases at a fine timestep size of ∆t = 2 • 10 -4 s so that errors are dominated by the spatial discretization. The maximum absolute errors over all space and for all timesteps are displayed in Figure 5 for both P and Q. The overlaid slopes in the plots illustrate the expected fifth-order accuracy of the Fourier continuation parameters employed in this paper (see Remark 3.1).

Algorithm 1 Summary of the FC-based solver for the LV-arterial system model. Construct the continuations P cont (x, t), Q cont (x, t) // via eq (13)

Input
4:

Compute Discrete Fourier Transforms of P cont , Q cont 5:

Apply the filter and compute termwise spatial derivatives // via eq (9)

6:

Invert the transforms to obtain ∂Pcont ∂x (x, t), ∂Qcont ∂x (x, t)

7:
Restrict the functions to the original interval // ∂P ∂x := ∂P cont ∂x {x i } , etc.

8:

Advance the solution including the boundary to t + ∆t // via eqs (15)

9:

Advance the boundary ODEs/PDEs (if they exist) to t + ∆t 10:

Update the boundary values on P (x, t), Q(x, t) to t + ∆t 11: end while Output the numerical solutions of P (x, t), Q(x, t) on 0 ≤ t ≤ t f , x = 0 ≤ x ≤ Figure 4: The numerical values at the indicated snapshot in time of the the solutions P (x, t) and Q(x, t) produced by the right hand side derived from solution [START_REF] Berger | Differential effects of wave reflections and peripheral resistance on aortic blood pressure: a model-based study[END_REF].

Figure 5: The maximum (L ∞ ) errors over all time and space after 50 000 timesteps for both P (x, t) and Q(x, t). Errors are calculated by comparison with the exact solution (equation ( 19)) and with the timestep fixed finely so that the errors are dominated by spatial convergence. The overlaid slopes illustrate the fifth-order convergence expected from employing d , dr = 5 matching points to construct Fourier continuations.

Dispersion and stability

Dispersion and stability characteristics of the proposed method can be studied by considering test problems for which waves propagate over long distances or long times. This is accomplished by considering a domain of size = 1cm and a number N λ of wavelengths across the plate for the MMS wave solution of period 1s given by

P Q = sin(2πN λ (x -ct)) sin(2πN λ (x -ct)) . (20)
The parameter c = 1/ √ LC is taken to be the pulse wave velocity and the physical domain is defined by C = 0.0221ml/(mmHg•cm), L = 0.000256mmHgs 2 /cm 4 and R = 0.000024mmHg s/cm 2 . Even though the given solution in equation ( 20) is periodic, the FC solver still extends the domain and creates a new periodic extension. Examples of the corresponding solutions for N λ = 5 and N λ = 25 are presented in Figure 6.

Simulation of propagation over long distances can be achieved by increasing the number N λ of wavelengths across the complete domain. With a timestep chosen small enough so that numerical errors are dominated by those arising from the spatial discretization, Figure 7 presents maximum numerical errors in the FCbased simulations of P and Q (labeled "FC") over all time and space as functions of the number N λ of Figure 7: The maximum numerical errors (FC and second-order FD) over all space and over one full temporal cycle (defined as the time required for one peak to travel the length of the aorta) for a plane wave solution with increasing number of wavelengths N λ .

wavelengths. Each wavelength is discretized by a fixed density of d λ = 10 and d λ = 15 points per wavelength. Overlaid in the figures are the corresponding maximum errors from the same solutions computed by a stable finite-difference solver (second-order central differences in space, explicit AB4 integration in time) using fixed densities d λ = 75 and d λ = 200 (labeled "FD"). Clearly, a (linear) degradation in accuracy of the FD solution is seen as N λ grows: even for low numbers N λ of wavelengths, larger and larger densities d λ would be necessary in order to produce reasonable experimentally-relevant accuracies. By contrast, the accuracy resulting from the FC algorithms remains essentially constant as N λ grows: the errors are effectively independent of frequency when the density d λ of points per wavelength is kept constant in the spatial discretization (that is, virtually no numerical dispersion). Similar results have been presented for other FC-based solvers in [START_REF] Albin | A spectral FC solver for the compressible Navier-Stokes equations in general domains I: Explicit time-stepping[END_REF][START_REF] Elling | GPU-accelerated Fourier-continuation solvers and physically exact computational boundary conditions for wave scattering problems[END_REF][START_REF] Lyon | High-order unconditionally stable FC-AD solvers for general smooth domains II. Elliptic, parabolic and hyperbolic PDEs; theoretical considerations[END_REF][START_REF] Albin | Fourier continuation methods for high-fidelity simulation of nonlinear acoustic beams[END_REF][START_REF] Amlani | An FC-based spectral solver for elastodynamic problems in general threedimensional domains[END_REF] when compared to high-order Padé schemes, (up to eighth-order) FD schemes or various orders of finite-volume methods (including hybrids with Discontinuous Galerkin methods). The advantages of FC demonstrated in all these cases result from their mild CFL restrictions or their ability to accurately approximate the dispersion characteristics of the underlying continuous operators at FFT speeds-a distinct advantage in modeling arterial pulse waves for hemodynamics research.

Stability can be studied by additionally considering solutions over several temporal cycles (long times), where the length of a cycle is the time required for any one crest to travel through the entire physical domain. Figure 8 presents the maximum numerical errors over all time and space resulting from an N λ = 25 wavelength solution of equation ( 20) with a cycle discretized by 5 300 timesteps (for a total of 80 000 timesteps in the simulation). As can be observed, the maximum error remains constant for arbitrarily long times and hence demonstrates long-term stability of the FC-based arterial wave solver.

FC-based simulations of the LV-arterial model and physiological relevancy

A number of factors influence the pumping ability of the left ventricle, including those related to the direct coupling between the left ventricle and the arterial system [START_REF] Berger | Differential effects of wave reflections and peripheral resistance on aortic blood pressure: a model-based study[END_REF][START_REF] Berger | Wave propagation in coupled left ventricle-arterial system[END_REF][START_REF] Campbell | Pulse reflection sites and effective length of the arterial system[END_REF]. The complexity of the dynamics that result from the interactions between these two elements can affect blood pressure, blood flow and pumping load. The theoretical model addressed in this paper (Section 2) has been selected not only for its physiological relevancy, but for its particular ability to capture these non-stationary and nonlinear dynamics of the complex coupling (which manifests as, alternatingly, boundary conditions determined by an ODE given by (4) during systole, and a zero-valued Dirichlet-type boundary condition during diastole when the aortic valve is closed). The resulting wave reflections within the coupled system play an important role in pressure and flow dynamics as well as certain diagnostic metrics.

For the simulations that follow, Table 2 presents the input values and corresponding physical dimensions for the LV-arterial system governed by equations ( 1), ( 4) and [START_REF] Bruno | Higher-order implicit-explicit multi-domain compressible Navier-Stokes solvers[END_REF]. These parameters are determined for their physiological relevance. The empirically given time-varying compliance C v (t) of the left ventricle (whose reciprocal is the corresponding elastance) is shown in Figure 9, and it can be determined by an analytical approximation to clinical data [START_REF] Senzaki | Single-beat estimation of end-systolic pressure-volume relation in humans: A new method with the potential for noninvasive application[END_REF][START_REF] Kass | Ventriculo-arterial coupling: Concepts, assumptions, and applications[END_REF]. Figures 10 and 11 demonstrate the physiological accuracy of the simulated results produced by the FC method applied to the LV-arterial model governed by the hyperbolic PDE in equation ( 1) and the boundary conditions defined by the ODEs in equations ( 4) and [START_REF] Bruno | Higher-order implicit-explicit multi-domain compressible Navier-Stokes solvers[END_REF]. These simulations are computed on a discretization composed of N = 64 points. The figures illustrate expected physiological features including the equality between the aortic pressure P (x = 0, t) and the ventricular pressure P v (t) during the systolic phase in the absence of a diseased valve (Figure 10). Additionally, the simulations capture the expected pulse pressure amplification as the pressure wave propagate downstream (Figure 11(a)) and the expected decrease of the flow amplitude as the flow wave propagates downstream through the vessel (Figure 11(b)). Figure 12 compares numerical values of pressure and flow at the aorta (inlet) for the FC-based arterial wave solver (labeled "FC") with a second-order central difference scheme in space ("FD") for N = 30 discretization points in the spatial domain. The numerical noise from the finite difference scheme is apparent in the provided inlet figures. This noise is particulary evident in the spurious non-physical numerical reflections that begin from the indicated time t = T d of the "dichrotic notch". The location of this notch represents the point of a small and brief increase in arterial pressure caused by the closing of the aortic valve. Accurate identification of the notch is important for advanced hemodynamics analysis such as those related to the buckberg index [START_REF] Hoffman | The myocardial supply: Demand ratio-a critical review[END_REF] as well as the recently-introduced cardiovascular Intrinsic Frequency (IF) method [START_REF] Pahlevan | Intrinsic frequency for a systems approach to haemodynamic waveform analysis with clinical applications[END_REF][START_REF] Pahlevan | Noninvasive iPhone measurement of left ventricular ejection fraction using intrinsic frequency methodology[END_REF]. Figure 10: Aortic input pressure P (x = 0, t) (blue circles) and left ventricular pressure Pv(t) (yellow) produced by the FC-based method using N = 64 discretization points. As expected, for a non-diseased aortic valve, the aortic pressure is equal to the ventricular pressure during the systolic phase (the time between the opening and the closing of the aortic valve).

The information contained in one cycle (of approximate length T per ) of the the aortic pressure P (x = 0, t) and flow Q(x = 0, t) at the inlet (e.g., Figure 12) is of particular interest in hemodynamics research and clinical practice since global cardiovascular indices such as the left ventricular pulsatile power and stroke work can be computed from these waveforms. Pulsatile power can be computed for a cardiac cycle of period T per from the numerical data by the expression pulsatile power := 1 T per Tper P (x = 0, t)Q(x = 0, t)dt -P (x = 0, t)Q(x = 0, t),

where P and Q represent the mean pressure and flow in the cycle. The left ventricular stroke work is defined as

stroke work := - V (Tper) V (0) P v (t)dV v = Tper P v (t)Q(x = 0, t). (22)
Another clinically-important quantity related to global LV function is the ejection fraction (EF). EF is the ratio of the volume pumped out of the left ventricle during each cardiac cycle to the volume of LV at the end diastole. Although EF is a geometric ratio, it is a surrogate for the LV contractility and is the most common index used in measuring LV systolic function in clinical practice (e.g., a low value of EF that is < 0.4 is indicative of heart failure). The ejection fraction is defined as

ejection fraction := V v (t) LV EDV = Tper Q(x = 0, t)dt, (23)
for time-varying ventricle volume V v (t) and preload LV end diastolic volume LV EDV . Figure 13 presents the numerical values of pulsatile power (left) and stroke work (right) comparing FC solutions (of sizes N = 48 and N = 64) with second-order (central) finite differences (of sizes N = 128 and N = 144). As is well-known, these values must be derived from the steady-state oscillatory dynamics of the underlying PDE (reached after approximately ten cycles in these simulations). Figure 13 additionally demonstrates that the measurements from finite differences are not reliable for many more preliminary cycles when compared to FC: the steady-state is numerically achieved several cycles after that of the FC-based solutions. This renders the FC method competitive for robustness of values for an arbitrary number of simulated cycles; these values stay relatively even for both N = 48 and N = 64 discretization points. In contrast, finite difference based solutions do not achieve reliable values (if at all) for many more cycles. This is important for feedback regulatory models requiring long-time simulations as well as for solutions to inverse problems related to arterial wave reflections. Figure 14 additionally presents the ejection fractions computed during each cycle for FC and FD. As is the case with the stroke work and pulsatile power, the steady state in values is achieved far earlier and remains stable for both discretizations of the FC-based simulations. This indicates that after the 11 th or 12 th cycle, EF is robustly calculated. In contrast, the finite difference method oscillates around the mean, leading to uncertainty in the computed EF as a function of cycle.

Conclusions

This work introduces a new high-order Fourier continuation-based method for resolving 1D arterial wave propagation problems in hemodynamics applications. The methodology is based on the Fourier continuation approach for resolving the Gibbs phenomenon, providing fast and accurate solutions with mild linear CFL constraints for stability and with effectively no numerical dispersion. This enables robust and precise solutions for long-time or many wavelength simulations that are required in the study of pressure and flow waveforms in the mammalian circulatory system. These unique advantages are particularly important when considering feedback regulatory models that require long-time simulations and when considering inverse problems (by iterative techniques) that are related to arterial wave reflections. A variety of numerical tests, including a demonstration of the physiological accuracy of the LV-arterial model (with complex hybrid ODE-Dirichlet boundary conditions), demonstrates the efficacy and versatility of the proposed solver. The methodology described herein can be easily extended to other 1D hemodynamics formulations including those with variable coefficients (e.g., localized changes to arterial wall compliance) and those exhibiting the nonlinear phenomena inherent to flow dynamics in large arteries. These features can be captured in a manner similar to other FC-based methods that have treated them successfully, including previous work by co-authors of this paper on spatially-varying material parameters for elastodynamic equations [START_REF] Amlani | An FC-based spectral solver for elastodynamic problems in general threedimensional domains[END_REF] as well as by others on nonlinear Burgers systems [START_REF] Bruno | Higher-order linear-time unconditionally stable alternating direction implicit methods for nonlinear convection-diffusion partial differential equation systems[END_REF] and nonlinear Euler equations [START_REF] Shahbazi | Multi-domain Fourier-continuation/WENO hybrid solver for conservation laws[END_REF][START_REF] Shahbazi | Multi-dimensional hybrid Fourier continuation -WENO solvers for conservation laws[END_REF]. This manuscript also represents a first step towards future work entailing coupling with 3D and multiphysics PDEs (e.g., mass transport equations) as well as building towards the 1D hemodynamic modelling of the entire circulation [START_REF] Boileau | A benchmark study of numerical schemes for one-dimensional arterial blood flow modelling[END_REF][START_REF] Mynard | One-dimensional haemodynamic modeling and wave dynamics in the entire adult circulation[END_REF]. The corresponding challenges, which include treatment of multiple segments (branching) and the subsequent numerical representation of their boundary conditions, may be possibly addressed by the domain decomposition strategies that have been developed for FC-based solvers and that have been shown to have negligible artificial reflections at interfaces (e.g., junctions) [START_REF] Albin | A spectral FC solver for the compressible Navier-Stokes equations in general domains I: Explicit time-stepping[END_REF][START_REF] Amlani | An FC-based spectral solver for elastodynamic problems in general threedimensional domains[END_REF][START_REF] Bruno | Higher-order implicit-explicit multi-domain compressible Navier-Stokes solvers[END_REF]. Additionally, in certain diagnostic metrics, it has been suggested that proper consideration of the geometry of the vasculature may be very important [START_REF] Schrauwen | Geometry-based pressure drop prediction in mildly diseased human coronary arteries[END_REF]; future work also entails the extension of the proposed methodology to curved elastic tube domains with curvilinear formulations of FC-based PDE solvers that are similar to [START_REF] Albin | A spectral FC solver for the compressible Navier-Stokes equations in general domains I: Explicit time-stepping[END_REF][START_REF] Amlani | An FC-based spectral solver for elastodynamic problems in general threedimensional domains[END_REF][START_REF] Bruno | Higher-order implicit-explicit multi-domain compressible Navier-Stokes solvers[END_REF]. Furthermore, the results of this paper imply that the FC-based methodology is straightforwardly applicable to other biological and physical phenomena that are governed by hyperbolic PDEs with nonlinear, highly-complex time-dependent boundary ODEs.

Figure 2 :

 2 Figure 2: An example Fourier continuation of a non-periodic function. The original function on [0, 1] is translated by a distance of length Ncont∆x whose values are filled-in by the sum of the blend-to-zero continuations (dashed lines) in order to render the function periodic. Green triangles and circles represent the discrete d , dr = 5 matching points, and the blue squares represent the discrete Ncont = 25 continuation points that comprise the extension.

Figure 6 :

 6 Figure 6: The numerical values at the indicated snapshop in time of the solution (20) for N λ = 5 (left) and N λ = 25 (right) wavelengths.

Figure 8 :

 8 Figure 8: The numerical errors for FC solutions considered over many temporal cycles (where one cycle spans 5 300 timesteps) for a tube N λ = 25 wavelengths in size, demonstrating long-time stability for both P (x, t) (left) and Q(x, t) (right).

Figure 9 :

 9 Figure 9: The time-varying left ventricular compliance Cv(t) employed in this section.

Figure 11 :

 11 Figure 11: (a) Pressure at the left and right boundaries (x = 0 and x = , respectively) capturing the expected amplification as the waves propagate downstream through the vessel. (b) The corresponding flow at the left and right boundaries that demonstrate the expected decrease in amplitude as the flow wave propagates downstream.

Figure 12 :

 12 Figure 12: Numerical values of pressure (left) and flow (right) at the inlet from a simulation computed by Fourier continuation and a (central) finite difference scheme. The numerical noise from the finite difference scheme is evident in the inset.

Figure 13 :

 13 Figure 13: The numerical values of the stroke work pulsatile power (right) computed for period of a thirty cycle simulation.

Figure 14 :

 14 Figure 14: The numerical values of the ejection fraction computed for each period of a thirty cycle simulation.

Table 1 :

 1 A summary of the physical parameters for the LV-arterial model that includes the left ventricle and the terminal compliance chamber. The end diastolic volume LV EDV and period of the cardiac cycle Tper are additional physiological input parameters that are used to fully describe the physical model.

	Name	Variable	Units
	Tube (domain) length		cm
	Tube (domain) compliance	C	ml/(mmHg cm)
	Tube (domain) inductance	L mmHg sec 2 /cm 4
	Tube (domain) resistance	R	mmHg sec/cm 4
	Terminal load compliance	C c	ml/mmHg
	Terminal load resistance	R c	mmHg sec/ml
	LV end diastolic volume (preload) LV EDV	ml
	LV dead volume	V dead	ml
	LV compliance	C v (t)	ml/mmHg
	Period of a cardiac cycle	T per	sec

  aorta parameters L, C, R, Input left boundary parameters C v (t), LV EDV, V dead Input right boundary parameters C c , R c Input size of discretization N Input number of cardiac cycles to simulate ( =⇒ final time t f )

	1: Initialize the pressure P (x, 0) and flow Q(x, 0)	// initial time t = 0
	2: while t < t f do
	3:

Table 2 :

 2 Input parameters for the simulations constructed by the LV-arterial system solver in Section 5.

	Name	Variable	Value	Units
	Tube (domain) length		21	cm
	Tube (domain) compliance	C	0.0221 ml/(mmHg cm)
	Tube (domain) inductance	L 0.000256 mmHg sec 2 /cm 4
	Tube (domain) resistance	R 0.000024	mmHg sec/cm 4
	Terminal load compliance	C c	1.05	ml/mmHg
	Terminal load resistance	R c	1.47	mmHg sec/ml
	LV end diastolic volume (preload) LV EDV	112	ml
	LV dead volume	V dead	4	ml
	LV compliance	C v (t) Figure 9	ml/mmHg
	Period of a cardiac cycle	T per	0.8	sec

Input: the valve condition and solutions P, Q at t, x = 0 ≤ x ≤ Advance the PDE by one timestep (lines 4-9 in Algorithm 1) Output: P (x, t + ∆t) and Q(x, t + ∆t), x = 0 ≤ x ≤ Advance boundary ODEs ( 4) and ( 6) by one timestep Output: P v (t + ∆t) and P c (t + ∆t) Solve for outflow via equation ( 7)

Output: the valve condition and solutions P, Q at t + ∆t, x = 0 ≤ x ≤