On the impact of modelling assumptions in multi-scale, subject-specific models of aortic haemodynamics

Simulation of haemodynamics has become increasingly popular within the research community. Irrespective of the modelling approach (zero-dimensional (0D), one-dimensional (1D) or three-dimensional (3D)), in vivo measurements are required to personalize the arterial geometry, material properties and boundary conditions of the computational model. Limitations in in vivo data acquisition often result in insufficient information to determine all model parameters and, hence, arbitrary modelling assumptions. Our goal was to minimize and understand the impact of modelling assumptions on the simulated blood pressure, flow and luminal area waveforms by studying a small region of the systemic vasculature—the upper aorta—and acquiring a rich array of non-invasive magnetic resonance imaging and tonometry data from a young healthy volunteer. We first investigated the effect of different modelling assumptions for boundary conditions and material parameters in a 1D/0D simulation framework. Strategies were implemented to mitigate the impact of inconsistencies in the in vivo data. Average relative errors smaller than 7% were achieved between simulated and in vivo waveforms. Similar results were obtained in a 3D/0D simulation framework using the same inflow and outflow boundary conditions and consistent geometrical and mechanical properties. We demonstrated that accurate subject-specific 1D/0D and 3D/0D models of aortic haemodynamics can be obtained using non-invasive clinical data while minimizing the number of arbitrary modelling decisions.


Introduction
Computational modelling of cardiovascular dynamics has received notable attention over the last two decades. The modelling approaches range from lumped parameter, zero-dimensional (0D) models of the circulation [1,2], to one-dimensional (1D) models of blood pressure and flow propagation [3][4][5][6][7], to three-dimensional (3D) fluid-structure interaction techniques [8][9][10][11]. Each approach has its own merits and limitations. Lumped parameter methods provide a computationally inexpensive, mathematically accessible and intuitive framework to study whole-system dynamics. However, they are not suitable for studying pulse propagation phenomena or complex flows. Nonlinear 1D methods can accurately describe pulse wave propagation phenomena in extensive vascular networks while keeping the computational cost down. However, these methods are not appropriate to describe complex 3D flow features, like those observed in stenosis and aneurysms. Lastly, 3D methods are capable of representing complex flows, wave propagation and blood flow-vessel wall interactions. On the downside, a 3D approach is computationally expensive.
Irrespective of the modelling approach, clinical measurements are required to personalize the geometry, material properties (e.g. vessel wall stiffness) and inflow and outflow boundary conditions of the computational framework. In vitro studies have shown the ability of 1D and 3D modelling to reproduce the main features of pressure and flow waveforms in the aorta and larger systemic arteries, if accurate measurements of in vitro data are available (e.g. [12,13]). There are, however, numerous challenges in acquiring haemodynamic data in vivo, particularly when using noninvasive techniques. For instance, magnetic resonance imaging (MRI) is limited by the acquisition time of certain sequences, their spatial and temporal resolution, and the variability in the physiology of the subject over the imaging study [14]. Ultrasonic techniques are affected by tissue signal attenuation and thus cannot reliably be used to measure flow in deep vessels [15]. Lastly, sphygmomanometry or tonometry are operator dependent, usually restricted to superficial vessels, and arguably introduce perturbations in the measured waveform, as they involve applying pressure in the region of interest [16]. Additional challenges arise from using measurements derived from different sources, which is common clinical practice. Indeed, these measurements are often not concurrent in time and, hence, subject to potential variability in heart rate (HR) and haemodynamic state, and may violate physical principles such as mass conservation.
The above limitations in data acquisition often result in insufficient information to determine the material and boundary condition parameters required for the simulation. Therefore, numerous modelling assumptions are usually made that are not directly informed by data. These modelling assumptions rely on morphometric considerations, such as more flow goes to larger vessels, specification of nonreflective outflow boundaries, and population-based rules for assigning tissue stiffness or pulse wave velocity.
In this article, we present a computational study of aortic pulse wave haemodynamics in a young healthy volunteer using multi-scale 1D/0D and 3D/0D frameworks and a rich array of non-invasive haemodynamic data, including magnetic resonance angiography, phase-contrast MRI, dynamic area and tonometry measurements. The clinical data and the set-up of the 1D and 3D frameworks were chosen to minimize the number of arbitrary modelling assumptions on boundary conditions and material parameters. This article has three main objectives: (i) to investigate the impact of different modelling assumptions on the simulated waveforms, (ii) to verify haemodynamic predictions obtained with 1D and 3D computational frameworks against clinical data, and (iii) to identify inconsistencies in the clinical data and discuss strategies to mitigate their impact on the simulation workflow.

In vivo data: measurement and post-processing
We acquired a series of non-invasive haemodynamic measurements on a young healthy male volunteer (age 27 years, weight 80 kg, height 188 cm) using a 1.5 T MRI scanner (Philips Achieva, Philips Healthcare) and applanation tonometry (SphygmoCor system, AtCor, Australia) at Guy's Hospital, London, UK. All measurements were taken at rest with the volunteer in the supine position. We followed a protocol approved by the local ethics committee and the volunteer provided written informed consent. A brief description of the data acquisition technology is given in appendix A.

Arterial geometry
A non-contrast, free-breathing, two-phase (systole and diastole) 3D SSFP scan (TE: 3 ms, TR: 5 ms, flip angle: 908, field of view: 320 mm, slice thickness: 2 mm) was carried out on the MRI scanner. The diastolic phase was used to reconstruct the 3D aortic geometry and generate finite-element meshes for the 1D and 3D simulations (figure 1), as detailed in §2.2.1.  measurement planes were approximately perpendicular to the aortic axis. Dynamic area images were acquired using gated twodimensional (2D) cine SSFP MRI on the same planes (TE: 1 ms, TR: 4 ms, flip angle: 508, field of view: 370 mm, slice thickness: 8 mm) to obtain aortic luminal area waveforms (figure 2c). Flow measurements to the upper branches were obtained via PC-MRI in a separate study in the same volunteer. This prior information was used to determine the flow distributions in the supra-aortic vessels (figure 2a). All flow and area measurements were acquired with 40 phases. The volunteer held his breath during each measurement to minimize movement artefacts. The automatic image segmentation and analysis tools of ViewForum (Philips Healthcare) were used to segment luminal cross-sectional areas from the dynamic area images. The software GTFlow (GyroTools LLC) was used to analyse the PC-MRI data and reconstruct flow waveforms.

Pressure pulse waves
Blood pressure waveforms were recorded in the left common carotid, left brachial and right radial arteries using applanation tonometry performed by an experienced operator using a SphygmoCor system (figure 3). Measurements were taken in a quiet environment just before and after the MRI study. A sliding stretcher was used to move the volunteer to maintain the HR as constant as possible.
For each measurement site, at least 10 cardiac cycles were obtained at a sampling frequency of 128 Hz and ensemble averaged by the SphygmoCor system. Waveforms that did not meet the in-built quality control criteria in the SphygmoCor system were rejected. All pressure waveforms were calibrated with brachial blood pressure measured in triplicate by a validated oscillometric method (Omron 705CP, Omron Health Care, Japan). Figure 3 shows the resulting left carotid, left brachial and right radial pressure waveforms calculated by averaging the corresponding measurements before and after placing the volunteer in the MRI scanner.

Numerical formulations
Detailed descriptions of the 1D/0D and 3D/0D formulations used in this study are given in [7] and [11], respectively. Here we summarize their main assumptions and describe their input physical parameters and simulated haemodynamic quantities that are relevant to this study.  rsif.royalsocietypublishing.org J. R. Soc. Interface 13: 20160073 In the 1D/0D framework, the aortic geometry was represented by interconnected segments, each one modelled as a deformable vessel with properties depending on a single axial coordinate x. The following haemodynamic quantities were simulated as a function of x and time t: the cross-sectional luminal area, Aðx, tÞ, and the cross-sectional averages of the axial blood flow velocity, Uðx, tÞ, and pressure, Pðx, tÞ. The axial velocity was assumed to have an axisymmetric profile.
In the 3D/0D formulation, vessel wall displacement and blood velocity were simulated as 3D vectors varying with time and space. Blood pressure was modelled as a function of time and the three spacial coordinates. Vessel wall displacements were described as a function of the blood velocities and pressures at the fluid -wall interface, using an 'enhanced' membrane formulation with a fixed fluid domain and linearized kinematics of the vessel wall.
In both the 1D/0D and 3D/0D models, blood was assumed to be incompressible and Newtonian, with a constant density r f and viscosity m. The arterial wall was modelled as a thin, incompressible, homogeneous, isotropic, linear elastic, impermeable membrane characterized by an elastic modulus E, Poisson's ratio n ¼ 0:5 and thickness h. In the 1D model, the membrane was also assumed to deform axisymmetrically, each cross section independently of the others. The following explicit algebraic relationship between P and A (or tube law) was used [4]: where A d ðxÞ is the luminal area at diastolic pressure (P d ) and bðxÞ accounts for the elastic material properties of the arterial wall. This tube law leads to the following expression for the pulse wave velocity (c) [7]: which relates the elastic parameter b(x) required to run the 1D model with the pulse wave velocity at diastolic pressure, c d , calculated from the in vivo data.
In the 3D formulation, no assumptions regarding axisymmetry were made. The arterial wall was modelled as a thin pre-stressed membrane with the following parameters, in addition to the elastic modulus and Poisson's ratio: density r s and transverse shear factor k [11]. The mechanical behaviour of the external tissue was modelled using stiffness (k s ) and damping (c s ) coefficients. This was important to account for the tethering exerted by the external tissue on the arterial walls and eliminate spurious and non-physiological oscillations [17].

Three-dimensional and one-dimensional mesh generation
The 3D aortic geometry was reconstructed from the magnetic resonance angiography data (figure 1a). Centrelines were first defined for the aorta and three supra-aortic vessels (figure 1b). Then, perpendicular to the centrelines, a series of 2D contours were created by segmenting the lumen of the vessel (figure 1b). These contours were used to generate a 3D parametric surface (figure 1c). A finite-element mesh (159 392 linear tetrahedral elements) for the 3D simulations was created using the MeshSim library (Simmetrix Inc., NY, USA).
To obtain the 1D network topology, we generated four centrelines (figure 1d ): one through the aorta and three through the supra-aortic vessels. Subsequently, cross sections were computed perpendicularly to the centrelines and spaced 20 mm apart, defining the 26 segments that form the 1D model mesh (figure 1e). Table 1 shows the length and radii at the inlet and outlet of each 1D model segment. Finite-element meshes for each segment were created with quadrature rules and polynomial functions of order 3.

Boundary conditions
Inlet and outlet boundary conditions were chosen to be consistent between the 1D and 3D schemes. At the inlet, we prescribed the flow waveform, Q in ðtÞ, measured by phase-contrast MRI, with the time period corrected as described in §2.4 (figure 2b, colour line). This waveform was repeated for at least four cardiac cycles in the 1D simulations and 10 cardiac cycles in the 3D simulations to achieve periodic solutions. The velocity profile at the inlet of the 3D domain was assumed to be axisymmetric in agreement with the 1D formulation profile.
The outlets of the descending aorta and three supra-aortic arteries were coupled to three-element Windkessel models. These 0D electrical circuit analogues of the downstream vasculature consist of proximal resistance, R 1 , connected in series with a parallel combination of a distal resistance, R 2 , and a compliance, C (figure 4). A non-zero outflow pressure, P out , was considered, which can be interpreted as the pressure at which flow to the microcirculation ceases.  In addition to the aortic geometry and the inlet flow boundary condition described above, material and outflow boundary parameters are required to run 1D and 3D simulations. This section describes the strategy for specifying the latter considering the available non-invasive data and a number of different modelling approaches. In particular, the following three independent parameters must be characterized: (i) pulse wave velocity at diastolic pressure (c d ), (ii) outflow pressure (P out ), and (iii) proximal resistance (R 1 ) for each Windkessel. We explored two different strategies to determine the value of each of these three parameters using the available data.
The following parameters of the 1D model were not directly calculated or estimated from the in vivo data: the blood density r f ¼ 1060 kg m À3 and viscosity m ¼ 3:5 m Pa s [18], the capillary pressure P cap ¼ 4:4 kPa [19] and the polynomial order of the 1D model velocity profile (z ¼ 9) [20]. For the 3D model, in addition to these four parameters, one must further specify values of wall thickness h, which was assumed to be 10% of the luminal radius, density of the wall r s ¼ 1000 kg m À3 , transverse shear factor of the arterial wall k ¼ 5=6 [11] and coefficients related to the stiffness, k s ¼ 0:01 Pa mm À1 and damping, c s ¼ 0:15 Pa s mm À1 , of the external tissue support [11,17].

Elastic mechanical properties
The elastic mechanical properties of the arterial wall were calculated from the pulse wave velocity, c. This was estimated from the data using the following two methods. Table 1. Geometric properties and pulse wave velocities of the upper aorta model measured non-invasively in a young healthy volunteer. r in ! r out : diastolic cross-sectional radii at the inlet and outlet of the arterial segment. c in ! c out : pulse wave speed (at diastolic pressure) at the inlet and outlet of the arterial segment calculated using either the foot-to-foot or QA -loop method.
foot to foot QA -loop

. Foot-to-foot method
A uniform c d ¼ 4:56 m s À1 was calculated and mapped to all 1D model segments. This value was obtained as the ratio of the centreline distance between the PC-MRI measurement sites at the ascending aorta and the most distal descending aorta (figure 1d) to the transit time of the foot of the flow waveforms at these two sites (figure 2b). The foot of the wave was defined using Gaddum et al.'s algorithm [21], as shown in figure 5a.

QA -loop method
We applied the algorithm described in [22] to calculate a local c d at the five aortic sites where blood flow and luminal area waveforms were measured. Following the approach described in [23], we related c d to the local luminal diameter at diastolic pressure, D d , calculated from the in vivo area waveforms by assuming a circular cross section (figure 5b). We then obtained a least-squares fitting of the form ð2:4Þ with j ¼ 0:72 m 3=2 s À1 . Table 1 shows the values of c d for each of the 1D model arterial segments. Equation (2.4) is in agreement with the form of equation ( It enables calculation of c d ðxÞ at any point in the 1D domain from D d ðxÞ under the assumption of a uniform b. Thus, with this approach, we obtain a spatially varying distribution of pulse wave velocities from multiple PC-MRI and area measurements down the aorta. Once c d is known, the elastic parameter bðxÞ required for the 1D models follows from equation (2.3). For the 3D models, the elastic modulus, EðxÞ, is required instead and was calculated using

Parameters of the Windkessel outlet models
We calculated the parameters of the four 0D Windkessel models at the outlets of the 1D domain using a linear analysis of the 1D/0D formulation previously presented in [24,25]. This analysis provides analytical equations relating the aortic root inflow wave, Q in ðtÞ, and a space-independent Windkessel pressure, p w ðtÞ, which approximates pressure waves throughout the computational domain. The parameters in these analytical equations are P out and the wall compliance, C T , and blood flow resistance, R T , of the entire systemic circulation, which in turn depend on the distributed compliances of the 1D model segments, and peripheral resistances and compliances of the 0D Windkessel models [24]. Appendix B provides further details of all the equations used in this section.

Total peripheral resistances
First, the resistance at the root of the 1D/0D arterial network, R T , was calculated as with P the mean blood pressure of the ensemble averaged carotid waveform (figure 3a) and Q in the mean flow rate at the aortic root (figure 2b). The outflow pressure P out was either set to reported values of capillary pressure, P cap ¼ 4:4 kPa [19], or calculated by fitting a curve of the form to the ensemble averaged carotid, brachial and radial pressures in diastole (figure 3, red lines), as described in [24]. Here, T 0 and p w ðT 0 Þ are the time and pressure, respectively, at the beginning of the fit, when pressure starts decaying exponentially, and t ¼ R T C T is the time constant of the exponential decay of pressure in diastole. The average values obtained from the three pressures in figure 3 are P out ¼ 9:2 + 0:2 kPa and t ¼ 0:32 + 0:03 s. Then, the total resistance at the outlet of each terminal vessel

Peripheral compliances
First, the total peripheral compliance, C p , was calculated using  Figure 5. Two different calculations of pulse wave velocity at diastolic pressure, c d , used in this study: (a) uniform c d in all 1D model arterial segments calculated using the foot-to-foot method applied to the flow waveforms at the Asc Ao (solid line) and Desc Ao 4 (dashed line); or (b) diameter-dependent c d calculated from estimates of c d at five aortic sites using the QA-loop method [22]. In (a), the transit time (TT) was calculated from the time of the foot of each wave, which was identified as the intersection (red circles) of a horizontal projection through the local minimum (black circles) and a tangential projection through the maximum gradient in systole (crosses), as detailed in [21]. The two sites are separated by Dx ¼ 28:5 cm measured along the aortic 3D centreline. In (b), D d is the luminal diameter at diastolic pressure. The curve is a least-squares fitting of the form c d ¼ jD À0:5 d (j ¼ 0:72 m 3=2 s À1 ) as described in [23].
rsif.royalsocietypublishing.org J. R. Soc. Interface 13: 20160073 with C T the compliance of the entire 1D/0D model network and C c the sum of the compliances of each 1D model segment, C seg , which were calculated using equation (B 6) in appendix B.
For the models in which the diastolic pressure decay was fitted to match the experimental pressure data (P out ¼ 9:2 kPa), C T was calculated from the time constant, t, and net resistance, R T , as [24] For the models in which the capillary pressure was assumed (P out ¼ P cap ¼ 4:4 kPa), we adopted [24] where Q max and Q min are the maximum and minimum flow rates calculated from the aortic root inflow Q in ðtÞ (figure 2b), Dt is the difference between the time at Q max and the time at Q min , and P s and P d are the systolic and diastolic pressures, respectively, of the ensemble average carotid waveform in figure 3a.
Once C p ¼ C T À C c was known, we then calculatedC j , j [ N t , by distributing C p to each terminal branch in proportion to the flow distribution used in equation (2.8) [26], that is The value of C j in the Windkessel model, which is different from the resistance-weighted complianceC j , is given by For the models in which C T was approximated by equation (2.11), we used the iterative procedure described in [27] to refine the values of C T and R T used to calculate R

Proximal resistances
The proximal resistance, R j 1 , in the three supra-aortic vessels (j ¼ 22, 24, 26) was assumed to be equal to the characteristic impedance, Z j 0 , of the end point of the terminal vessel j, that is where c j d and A j d are, respectively, the wave speed and area at diastolic pressure at the end point. This choice of R j 1 minimizes the magnitude of waves reflected at the outlet of the 1D model upper branches [25].
At the descending aorta outlet (j ¼ 20), a proximal resistance, R Ao 1 , was calculated to reproduce characteristic inflection points of the carotid pressure waveform (figure 3a). This modelling assumption implies that the main site of wave reflections in our model was located in the descending aorta, in agreement with clinical observations [28], whereas the upper branches were treated as non-reflective boundaries. The following equation was derived to calculate the value of R Ao 1 (see appendix C): Z Ao 0 is the characteristic impedance at the outlet of the descending aorta, computed using equation (2.14) with the parameters of the descending aorta (Vessel 20), P s ¼ 15:7 kPa is the systolic pressure of the ensemble average carotid waveform, P out ¼ 9:2 kPa, and P si ¼ 14:1 kPa is the pressure at the inflection point during the pressure decay in the second half of systole (figure 3a). The inflection point was calculated as the local maximum of the first derivative of the ensemble average carotid waveform. For our volunteer, we obtained R Ao 1 ¼ 0:60Z Ao 0 and a negative reflection coefficient of V ¼ À0:25 calculated using equation (B 7).

Reconciling data inconsistencies
We made every effort to acquire the in vivo data in the volunteer under similar physiological conditions, since we aimed to simulate waveforms for a single physiological state. However, we observed several data inconsistencies in (i) HR, (ii) MRI-derived vessel geometry, and (iii) mean volume flow rates. This section describes how these inconsistencies were addressed to mitigate their impact on the simulation workflow.

Heart rate
The average HR calculated from the flow (61.5 + 3.4 bpm) and area (65.2 + 1.5 bpm) measurements was greater than that calculated from the pressure (59.3 + 2.4 bpm) measurements, possibly because MRI acquisition required the volunteer to be constrained inside the scanner for almost 2 h and to perform several breath holds. The impact of this inconsistency was mitigated by calculating an average HR representative of all in vivo measurements (61.2 + 3.4 bpm). This average HR was then assigned to the aortic inflow waveform (figure 2b, colour line).

Magnetic resonance imaging-derived vessel geometry
Differences in diastolic areas were observed between 3D SSFP (used to define the 3D aortic geometry) and 2D SSFP MRI acquisitions (used to generate the dynamic area images). These can be explained by the different spatial resolutions of the 2D and 3D data used. To reduce scanning time, the 3D SSFP anatomy images were acquired with a lower spatial resolution than the 2D cine SSFP dynamic area images (1.5 versus 1.0 mm 3 , respectively). Therefore, we assumed that the 2D data capture better the true diastolic areas and corrected the 3D geometry by matching the diastolic cross-sectional areas given by the five aortic area waveforms in figure 2c. This was achieved by (i) computing the ratio of the measured diastolic areas using 2D cine SSFP at the five aortic locations (green contours in figure 1d) to the actual cross-sectional areas obtained from slicing the 3D mesh with the measurement planes, (ii) rescaling the 2D contours associated with the aortic 3D mesh, and (iii) regenerating the 3D geometry to match the measured areas in those five locations. The aortic 2D contours in-between measurement locations were scaled with a linearly interpolated scaling factor from the nearest two measurement locations. The 2D contours associated with the arch branches were re-scaled using an average of the scaling factors in the aorta.

Mean volume flow rates
Measures of flow rate by PC-MRI are less exposed to errors in contour segmentation than area measurements. This is because the flow is approximately zero near the arterial wall, and therefore errors in contour segmentation are likely to have a small impact on the overall flow rate. As a result, blood flow data acquired by MRI should be more accurate than luminal area data. However, mean volume flow rates computed from the aortic PC-MRI increased towards the periphery, from Desc Ao 1 to Desc Ao 2 or 3 (figure 2b), probably due to the smaller HR recorded in Desc  Ao 1 (figure 2b). This violates the principle that mean aortic flow decreases as we move down the aorta due to perfusion to small side branches. To mitigate this inconsistency, we calculated the percentage of cardiac output leaving the aorta in equation (2.8) using the mean value among the four aortic flow waves measured by PC-MRI at the locations labelled Desc 1 to Desc 4 in figure 2a.

Results and discussion
We investigated multiple modelling assumptions in which vessel stiffness and boundary condition parameters were defined via the different methods described in §2.3. Using the 1D/0D framework for its computational efficiency, rsif.royalsocietypublishing.org J. R. Soc. Interface 13: 20160073 multiple permutations of uniform and non-uniform pulse wave velocity, reflective and matched proximal outflow resistances, and capillary and data-fitted P out were tested against the available in vivo data. These are listed in appendix C. We obtained two different sets of modelling assumptions that produced the smallest relative errors 1 for aortic area and carotid pressure. Hereafter, we refer to these sets as the 'best area' ( §3.1) and 'best pressure' ( §3.2) models. We then studied the impact of uniform versus non-uniform pulse wave velocities ( §3.3), and matched versus reflective proximal outflow resistances ( §3.4) on these models. Lastly, 3D and 1D simulations were compared for the 'best area' and 'best pressure' models by using compatible material laws and identical inflow and outflow boundary conditions in both modelling frameworks ( §3.5).    Figure 6. (a -c) Waveform predictions by the 'best area' (colour solid lines) and 'best pressure' (colour dashed lines) models in the five aortic sites where MRI measurements were taken (first five rows) and at the outflow of the left common carotid artery (CCA, last row). In vivo waveforms (black lines) are shown for the aortic flows, aortic areas and carotid pressure, together with average (avg), maximum (max), systolic (sys) and diastolic (dias) relative error metrics (first column: 'best area' model; second column: 'best pressure' model).

Best aortic area waveform predictions
rsif.royalsocietypublishing.org J. R. Soc. Interface 13: 20160073 corresponding numerical predictions calculated by both the 'best area' and 'best pressure' 1D models. Aortic pressures and carotid flows and areas are also displayed for both models, though corresponding in vivo data were not acquired.
In this section, we focus on the results obtained by the 'best area' model (colour solid lines), for which all four terminal branches have an outflow pressure (P out ) equal to the capillary pressure given in [19] (P out ¼ P cap ¼ 4:4 kPa), and hence P out is not specific to our volunteer. The three supra-aortic vessels are coupled to matched Windkessel models (i.e. R 1 ¼ Z 0 ), while the descending aorta has a reflective R Ao 1 ¼ 0:60Z Ao 0 calculated using equation (2.15). Elastic mechanical properties were estimated assuming a uniform pulse wave velocity (c d ¼ 4:56 m s À1 ) calculated using the foot-to-foot method.
The 'best area' model is able to capture well the overall shape of in vivo aortic flow and area waves ( figure 6a,b). The following features are well described by the model: the time and magnitude of the feet of the waves, the amplitude of the flow and area waves, the skewed flow peak in systole, the relatively smaller flow in diastole and the decay in diastolic area. Predictions for aortic flow waveforms show oscillations in systole and early diastole not seen in the in vivo measurements. Despite these oscillations, average relative flow errors for the flow do not exceed 7%. Discrepancies in area predictions occur mainly in systole, leading to average relative area errors smaller than 4%.
Systolic and diastolic relative errors for left carotid pressures became smaller than 1% after one iteration of the Windkessel compliances and resistances (figure 6c). The simulated carotid pressure wave, however, overestimated the in vivo pressure from the systolic peak to the end of the cardiac cycle, leading to an average relative pressure error of almost 10%.

Best carotid pressure waveform prediction
The outflow pressure P out is the only different independent parameter between the 'best area' and 'best pressure' models: P out ¼ 9:2 kPa for the 'best pressure' model versus P out ¼ 4:4 kPa for the 'best area' model. Therefore, P out plays a very important role in shaping arterial pulse waveforms (figure 6). Unlike in the 'best area' model, P out ¼ 9:2 kPa is specific to our volunteer: it is based on an exponential fit of the form given by equation (2.7) to the decay in pressure during diastole in the carotid, brachial and radial arteries (figure 3). We have obtained very similar values of t and P out in these three vessels (figure 3) which suggests that late-diastole pressure is approximately uniform in space, in agreement with the derivation of equation (2.7) [24].
The predicted carotid pressure wave by the 'best pressure' model features the systolic and diastolic inflection points observed in vivo, captures the time and shape of the systolic peak better than the 'best area' model and reproduces well the decay in pressure in diastole (figure 6c). As a result, average and maximum relative errors for carotid pressure are considerably reduced with respect to the 'best area' model. However, aortic areas are not predicted as well as by the 'best area' model (figure 6b): area is underestimated in all five aortic sites during most of the cardiac cycle by around 10%. Given that average relative errors in in vivo aortic flow are smaller for the 'best pressure' model from Asc Ao to Desc Ao 2 and for the 'best area' model in Desc Ao 3 and 4 (figure 6a), we cannot determine the best choice for P out from the in vivo data that are available for this study. Invasive in vivo aortic pressure waveforms would be required for this purpose.
By correcting the 3D arterial geometry based on the five 2D SSFP area measurements ( §2.4.2), we obtained, in both the 'best area' and 'best pressure' models, relative errors for diastolic areas smaller than 1% in all aortic sites except for Desc Ao 1, where errors were up to 5% (figure 6b). This is because the PC-MRI plane selected to acquire the in vivo luminal contours at Desc Ao 1 was more oblique than the planes considered at the other four aortic locations (figure 1d): indeed, the angle between the MRI contours (shown in green) and the contours on planes perpendicular to the aortic axis used to generate 1D model areas (shown in black) was 178 at Desc Ao 1 and smaller than 98 at the other four aortic sites. Automatic calculation of planes perpendicular to the aortic axis by the MRI acquisition software may, therefore, improve the quality of the PC-MRI data required for subject-specific modelling.

Uniform versus non-uniform pulse wave velocity
A uniform pulse wave velocity, c d , calculated using the foot-tofoot method (figure 5a) was used in both 'best area' and 'best pressure' models to determine the elastic properties of all 1D model arterial segments. Using a non-uniform c d ðxÞ calculated by the QA-loop method (figure 5b) also predicted well the arrival time of the feet of in vivo aortic flow, aortic area and carotid pressure waveforms (figure 7 shows the results for the 'best area' model only). However, in both the 'best area' and 'best pressure' modelling assumptions, a distributed c d ðxÞ increased relative errors for aortic areas (figure 7b). These were consistently underestimated by a non-uniform c d ðxÞ, which suggests that this modelling assumption led to an overestimation of aortic stiffness. Indeed, in all arterial segments (except for those in the ascending aorta), non-uniform values of c d ðxÞ are greater than the uniform c d ¼ 4:56 m s À1 (see table 1), which leads to greater elastic moduli EðxÞ (stiffer walls) according to equation (2.5). Stiffer arterial walls reduced the total compliance compared with the uniform c d case, preventing the iterative process described in [27] from producing the target in vivo pulse pressure and, hence, increasing systolic and diastolic relative errors for the carotid pressure (figure 7c).
In the 'best pressure' model, we also observed overall larger relative errors for aortic flow when a non-uniform c d ðxÞ was used. In the 'best area' model, relative errors in aortic flow waveforms decreased in descending aorta sites if a nonuniform c d ðxÞ was used (figure 7a). In either model, however, the main changes introduced by a non-uniform c d ðxÞ were observed in area and pressure waves of the aorta and carotid artery, rather than in the flow waves. According to these results, the foot-to-foot method provides a better estimate of pulse wave velocity in the upper thoracic aorta than the QA-loop method. Indeed, considerable relative errors (over 30%) in the estimates of c obtained by the QA-loop method have been reported in [29,30].

Matched versus reflective proximal outflow resistances
An important result of our study is the fact that a reflective resistance R 1 was required at the outflow of the descending aorta (R Ao calculated by the 'best pressure' model, with a reflective aortic R Ao 1 computed using equation (2.15), with the corresponding waveforms simulated using a well-matched R 1 ¼ Z 0 in all four terminal branches. A reflective R Ao 1 had the opposite effect on aortic flows to aortic pressures/areas (figure 8): it increased systolic flows and decreased systolic pressures/areas, leading to smaller average relative errors than those produced by a matched R Ao 1 . The effect of the reflective R Ao 1 is in agreement with the linear 1D pulse wave theory used to obtain equation (2.15): the theory described in [25] shows that a negative reflection coefficient at a 1D model terminal end (V ¼ À0:25 in our case) will decrease the amplitude of reflected pressure wavefronts and will increase the amplitude of reflected flow wavefronts. Moreover, a reflective R Ao 1 was required to generate the diastolic pressure peak observed in vivo (figure 8c), which suggests that reflected waves in the lower body travelling upstream of the aorta play an important role in shaping aortic and carotid pulse waveforms. This result is consistent with clinical observations [28].
Similar results were obtained when comparing matched versus reflective R 1 for the 'best area' model (not shown). Lastly, we note that using equation (2.15) to calculate R 1 in the three supra-aortic vessels had an insignificant effect on the waveforms produced by the 'best area' and 'best pressure' models.

One-dimensional versus three-dimensional haemodynamics
We performed a comparison between 1D and 3D results for the 'best area' and 'best pressure' modelling assumptions. In general, 1D and 3D theories produced similar flow, area and pressure waveforms with comparable relative errors at sites where in vivo data were available. Figure 9 displays 1D and 3D results computed under the 'best area' assumptions, together with in vivo waveforms, at the Asc Ao, left CCA, Desc Ao 2 and Desc Ao 4 locations. Figure 10 shows the corresponding results under the 'best pressure' assumptions. Simulated aortic flow, aortic area and carotid pressure waveforms using 1D and 3D modelling were both able to capture the main features of corresponding in vivo waveforms. This is confirmed by average relative errors smaller than 7%, 4% and 11%, respectively, using the 'best area' assumptions and smaller than 9%, 8% and 4%, respectively, using the 'best pressure' assumptions. Such small errors suggest that subject-specific 1D/0D and 3D/0D models are both able to capture the main features of in vivo aortic flow, aortic area and carotid pressure waveforms under normal anatomical and physiological conditions.  In vivo measurements were not available to test the accuracy of simulated aortic pressures, carotid flow and carotid area waveforms. For these waveforms, we compared 1D model against 3D model predictions. Previous studies [5,27] had carried out 1D versus 3D comparisons in idealized geometries, but not in subject-specific geometries as we have done in this study. Our results showed that simulated aortic pressures, carotid flow and carotid area waves using the 1D formulation contained the main features of corresponding 3D model waves. Average relative errors were smaller than 2%, 5% and 3%, respectively, in both the 'best area' and 'best pressure' models, which are similar to corresponding errors obtained in idealized geometries [5,27]. The greatest discrepancies were observed in mid-systole, as we had previously observed in idealized geometries [27]: the 1D formulation underestimates 3D model systolic pressures along the aorta, with differences becoming smaller towards the periphery. A greater pressure gradient is therefore necessary between the root and outflow of the aorta in the 3D model during mid-systole. This indicates that more energy is dissipated along the 3D aortic model: part of the energy is consumed in driving secondary flow patterns that develop in the aortic arch due to the curvature of the vessel and which are not captured by the 1D formulation. The small relative errors between 1D and 3D modelling approaches suggest that our combined 1D/0D method for calculating mechanical properties and outflow boundary conditions for 3D aortic flow modelling-introduced in [27] for idealized aortic and carotid geometries-can also be applied to subject-specific geometries under normal physiological and anatomical conditions; for example in the absence of occlusive or aneurysmal disease.
We note that the overall shape of the flow waveform in the left common carotid artery obtained by the 'best area'   Figure 9. Blood flow (Q), luminal area (A) and blood pressure (P) waveforms computed using the 1D (red lines) and 3D (blue lines) 'best area' models in four arterial sites: (a) ascending aorta, (b) left common carotid artery (CCA) (one diameter away from the outlet), (c) descending aorta 2 and (d ) descending aorta 4. Available in vivo waveforms at these sites are shown in black lines. Errors are shown for the 1D model (first column) and the 3D model (second column) relative to the in vivo data, if available; else only for the 1D model relative to the 3D model. and 'best pressure' models is qualitatively similar to that measured with PC-MRI in [6] for a younger normal subject. The predicted flow is unidirectional for the entire cardiac cycle and contains a characteristic peak in early diastole (figures 9b and 10b). Previous studies had verified 1D and 3D model waveforms by comparison against in vivo data in more extensive vascular networks than the thoracic aorta used here [3,6,31]. While including additional arteries will make the computational domain more complete, it also makes it more challenging to personalize all modelling parameters to a given subject. In previous studies, numerous parameters had to be obtained from the literature or tuned; for example using optimization techniques [17]. As a result, only a small proportion of the total number of parameters could be determined from the available in vivo data. Here we have minimized the number of arbitrary parameters by focusing on a confined region of the systemic vasculature and acquiring a rich array of non-invasive MRI and tonometry data. Therefore, all parameters of the aortic 1D/0D models were directly calculated or estimated from the data, except for the blood density r f ¼ 1060 kg m À3 and viscosity m ¼ 3:5 mPa s [18], the polynomial order z ¼ 9 of the velocity profile [20], and, for the 'best area' model, the capillary pressure P cap ¼ 4:4 kPa [19]. For the 3D/0D models, in addition to these four parameters, we had to arbitrarily specify values of wall density and thickness, mechanical properties of external tissue support, and transverse shear factor.
Lastly, we note that blood viscosity had a small effect on the aortic waveforms computed by the 1D model: the inviscid 1D solution changed average relative errors for flow, area and pressure waveforms by less than 0.4% in both the 'best area' and 'best pressure' models (see appendix C). According to this result, the inviscid 1D solution could also be used to simulate pulse wave propagation in the upper thoracic aorta, under normal anatomical and physiological conditions. In this scenario, we avoided estimating blood viscosity and the shape of the velocity profile, reducing to two the number of 1D model parameters that are not subject specific for the 'best area'   Figure 10. Results for the 'best pressure' modelling assumptions (same format as in figure 9). rsif.royalsocietypublishing.org J. R. Soc. Interface 13: 20160073 model (the blood density, r f , and capillary pressure, P cap ) and to only one parameter (r f ) for the 'best pressure' model.

Clinical implications
The 1D/0D methodology described in this study enables calculation of pressure, flow and area waves in the upper aorta and supra-aortic vessels from non-invasive measurements in a discrete number of locations. It provides a means of calculating patient-specific estimates of haemodynamic quantities that are relevant to assess cardiovascular function, such as aortic (central) blood pressure [32], aortic pulse wave velocity [33] and wall shear stress [18]; the latter using the 3D/0D formulation with parameters determined by the 1D/0D methodology. However, the 1D/0D algorithm has only been verified in one subject and using non-invasive data. Further verification in a clinical cohort is, therefore, essential before our method could become a clinical tool for non-invasive detailed characterization of aortic haemodynamics. Verification against intra-vascular aortic pressures will allow us to test aortic pressure predictions and elucidate whether the 'best area' model provides more accurate results than the 'best pressure' model or vice versa. Clinical cohorts with different cardiovascular pathologies should also be investigated, for example, to determine the applicability limits of this 1D/0D method in anatomies such as aortic coarctation or aneurysms.
In §2.4, we have discussed the importance of reconciling several inconsistencies in the clinical data used for calibrating our 1D and 3D models. In addition, the impact of uncertainties affecting the measured data should be analysed before translating the method to the clinic. For example, uncertainties in the flow split at the supra-aortic vessels measured by PC-MRI have been shown to have a significant effect on wall shear stress-based indicators [34] and may also affect considerably the pattern of predicted pressure, flow and area waveforms.

Conclusion
We have shown that accurate, subject-specific, 1D/0D and 3D/0D models of pulse wave haemodynamics in the upper aorta of a young healthy volunteer can be obtained using non-invasive clinical data. By simulating blood flow in a confined region of the systemic vasculature and acquiring a substantial amount of in vivo measurements, we have minimized the number of arbitrary modelling assumptions and determined most of the model parameters from the in vivo data. We have provided a detailed workflow for calculating the geometrical and mechanical properties of the computational domains, as well as outflow boundary conditions, from non-invasive data acquired by MRI and tonometry. The following are key aspects of this workflow to minimize relative errors for aortic flows, aortic areas and carotid pressures computed using the 1D/0D formulation: (i) elimination of inconsistencies in the clinical data, (ii) a uniform pulse wave velocity calculated from PC-MRI flow waves by the foot-to-foot method, (iii) a reflective resistance at the outflow of the descending aorta calculated from the shape of the carotid pressure wave, and (iv) an outflow pressure either equal to the capillary pressure to better predict aortic area waves or calculated from the tonometry pressure waves to better predict the carotid pressure wave. We have also demonstrated good agreement between 1D and 3D model predictions, especially during early systole and diastole. This study supports the use of 1D and 3D models for subjectspecific modelling of aortic pulse wave haemodynamics, as well as the use of the 1D/0D formulation for an efficient calculation of parameters for 3D/0D modelling.
Data accessibility. The in vivo datasets supporting this article have been uploaded as part of the electronic supplementary material.

B.1. Total peripheral resistances
By assuming periodic flow and neglecting nonlinearities and viscous dissipation in the 1D model segments of an arterial network coupled to three-element Windkessel outlet models (figure 4), the time-average pressure over one cardiac cycle, P, measured in any 1D model segment (the aorta and supra-aortic vessels in this study) is given by [24] P P out is the pressure at the outflow of each Windkessel model, assumed to be the same at each outlet, R T is the resistance at the root of the 1D/0D arterial network, which depends on the individual resistances of the 0D Windkessel models, N t is the set of arterial segments that are terminal branches and Q in is the time-averaged inflow waveform. For the upper thoracic aorta considered in this study, Q in is the mean flow rate at the ascending aorta (figure 2b) and Table 3. Modelling assumptions and average relative errors for all tested permutations. The minimum average relative errors among all modelling assumptions are highlighted in bold. rsif.royalsocietypublishing.org J. R. Soc. Interface 13: 20160073 N t ¼ f20, 22, 24, 26g (figure 1e). Rearranging equation (B 1) and taking P to be the pressure wave at the carotid artery leads to equation (2.6). By applying equation (B 1) to a network consisting of a single 1D model vessel (each terminal branch j [ N t in this study), we can write P ¼ P out þ ðR with Q j out the mean flow at the outlet of the vessel (which is equal to the mean flow at the inlet for a single vessel). Equation (2.8) follows from combining equations (B 1) and (B 2) and solving for R j 1 þ R j 2 .
As also detailed in [24], pressure waves in any 1D model segment of the network can be approximated by a spaceindependent 0D Windkessel pressure, p w ðtÞ, which depends on the 1D/0D model parameters of the network. This is achieved by neglecting nonlinearities, blood flow inertia and viscous dissipation in 1D model segments. Under these assumptions, p w ðtÞ is given by where p w ðT 0 Þ is the pressure p w at a reference time t ¼ T 0 , C T is the compliance of the entire 1D/0D arterial network, Q in ðtÞ is the inflow waveform and Q j out ðtÞ is the outflow at each terminal branch. For the upper thoracic aorta considered in this study, Q in ðtÞ is the blood flow at the ascending aorta shown in figure 2b (colour line). Figure 3a shows the spaceindependent pressure, p w ðtÞ, calculated using the parameters of the 'best pressure' model, with T 0 corresponding to the time at the start of systole.
By taking T 0 to be the time when pressure starts decaying exponentially in diastole, and assuming Q in ðtÞ ¼ dQ j out ðtÞ=dt ¼ 0 for t ! T 0 , j [ N t -note that Q in ðtÞ is approximately zero for t ! T 0 (figure 2b)-equation (B 3) reduces to equation (2.7), with t ¼ R T C T .

B.2. Peripheral compliances
The total compliance C T consists of both 1D and 0D compliances [24]: where C c is the total compliance of the 1D model arterial network and C p is the total peripheral compliance. Equation (B4) can be rearranged in the form of equation (2.9). The compliances C c and C p depend on, respectively, the compliance of each 1D model segment (C seg ) and the parameters of the 0D Windkessel outlet models (R j 1 , C j and R j 2 , with N the total number of segments in the 1D model arterial network (N ¼ 26 in this study) andC the resistance-weighted compliance of each Windkessel outlet model. Rearranging the third expression in equation (B5) yields equation (2.13). The compliance C seg of each 1D model segment can be calculated as with L the segment length, and A d ðxÞ and c d ðxÞ, respectively, the area and wave speed at diastolic pressure.

B.3. Proximal resistances
The linear analysis presented in [25] shows that a pressure wavefront, dp, propagating towards the outlet of a straight 1D model arterial segment coupled to a single-resistance outlet, R 1 , produces a reflected pressure wavefront, dp, given by V is the terminal reflection coefficient and Z 0 is the characteristic impedance of the 1D model segment. According to equation (B 7), the incoming pressure wavefront is completely absorbed if R 1 ¼ Z 0 , leading to equation (2.14) for calculating R 1 in the three supra-aortic vessels. At the outlet of the descending aorta, a reflective R 1 was calculated by assuming that the amplitude at the inflection point, P si À P out , of the carotid pressure waveform (figure 3a) is caused by the reflection at the aortic outlet of a wave of amplitude P s À P out propagating from the ascending aorta in early systole. Thus, equation (2.15) follows from taking R 1 ¼ R Ao 1 , Z 0 ¼ Z Ao 0 , dp ¼ P si À P out and dp ¼ P s À P out in equation (B 7). Table 3 lists all the different permutations of pulse wave velocity, c, proximal aortic outflow resistance, R Ao 1 , and terminal outflow pressure, P out , that were tested against the available in vivo data using the 1D/0D framework. For each set of modelling assumptions, average relative errors between in vivo and predicted waveforms are provided at the same arterial sites as in figure 6. The last two rows show the results for the 'best pressure' and 'best area' models under the assumption of inviscid flow; that is m ¼ 0.