The Effect of Arterial Curvature on Blood Flow in Arterio-Venous Fistulae: Realistic Geometries and Pulsatile Flow

Arterio-Venous Fistulae (AVF) are regarded as the “gold standard” method of vascular access for patients with End-Stage Renal Disease (ESRD) who require haemodialysis. However, up to 60% of AVF do not mature, and hence fail, as a result of Intimal Hyperplasia (IH). Unphysiological flow and oxygen transport patterns, associated with the unnatural and often complex geometries of AVF, are believed to be implicated in the development of IH. Previous studies have investigated the effect of arterial curvature on blood flow in AVF using idealized planar AVF configurations and non-pulsatile inflow conditions. The present study takes an important step forwards by extending this work to more realistic non-planar brachiocephalic AVF configurations with pulsatile inflow conditions. Results show that forming an AVF by connecting a vein onto the outer curvature of an arterial bend does not, necessarily, suppress unsteady flow in the artery. This finding is converse to results from a previous more idealized study. However, results also show that forming an AVF by connecting a vein onto the inner curvature of an arterial bend can suppress exposure to regions of low wall shear stress and hypoxia in the artery. This finding is in agreement with results from a previous more idealized study. Finally, results show that forming an AVF by connecting a vein onto the inner curvature of an arterial bend can significantly reduce exposure to high WSS in the vein. The results are important, as they demonstrate that in realistic scenarios arterial curvature can be leveraged to reduce exposure to excessively low/high levels of WSS and regions of hypoxia in AVF. This may in turn reduce rates of IH and hence AVF failure.


INTRODUCTION
End-Stage Renal Disease (ESRD) is characterized by an irreversible loss of kidney function, and is treated predominantly with haemodialysis or renal transplantation. 25,27,61 Efficient haemodialysis is dependent upon high-quality vascular access, allowing the rapid removal of blood (at up to 450 mL min À1 ), which passes through a dialyser removing metabolic waste and water, before being returned to the body. The preferred method of vascular access 60 is via an established Arterio-Venous Fistula (AVF) created surgically by connecting a patient's own artery and vein, usually in the arm. The connection or ''anastomosis'' can take a number of configurations. However, an end-to-side arrangement, where the end of the vein is connected to the side of the artery, has emerged as the clinically preferred configuration. 36 Following creation of an AVF, there is a rapid increase in blood flow within the artery and vein. In an ideal situation the venous lumen proceeds to enlarge, along with concurrent remodeling of the vessel wall. 20 After a period of six to twelve weeks the vein should have ''matured'' such that it can support insertion of two large gauge needles several times a week -providing access for high quality haemodialysis over many years. Unfortunately, up to 60% of AVF do not mature as anticipated 1 resulting in unacceptable patient morbidity and a significant burden on healthcare expenditure. Failure occurs predominantly due to the development of Intimal Hyperplasia (IH). Material obtained from immature dialysis access reveals a hyperplastic lesion arising from the intima that progressively occludes the lumen of the vessels, resulting in slow flow, and ultimately thrombosis and irreversible loss of the AVF. 2,50 Specifically, IH is found in the venous section of the AVF close to the anastomosis 2,8,56 and more distally, 8,56 as well as in the arterial section of the AVF close to the anastomosis [6][7][8]56 and more proximally. 6,8 Whilst the exact mechanisms underlying development of IH are unknown, there is considerable evidence to suggest that the inherently unphysiological flow patterns within AVF play an important role. In particular regions of highly oscillatory flow, abnormally low/high Wall Shear Stress (WSS), and low Lumen-to-Wall Normal Oxygen Flux (LWNOF)-leading to wall hypoxia, have all been implicated in the initiation of IH. 10,13,14,16,19,22,23,26,30,38,39,45,46,49,58,59 In 2015 Iori et al. 32 investigated the effect of planar arterial curvature on blood flow and oxygen transport in AVF using an idealized model. Specifically, idealized straight venous sections were connected to either the inside or outside of a curved planar arterial bend, with nonpulsatile inflow conditions. It was found that inner configurations suppressed potentially pathological low WSS and low LWNOF but, conversely, outer configurations suppressed potentially pathological high-frequency flow unsteadiness. The present study extends this work to more realistic non-planar AVF configurations with pulsatile inflow conditions.

Geometries Native Vessels
Native arterial geometries from two ''patients,'' henceforth referred to as P1 and P2, were recon-structed using Computed Tomography (CT) scan data of the left arm and hand. The anonymized CT datasets were obtained retrospectively from adult subjects who had undergone clinically-indicated upper-body CT angiography in the setting of major trauma and were found to have normal vasculature. Specifically, the P1 dataset was collected using a Philips Brilliance 64 Helical CT scanner (841 slices in total with a voxel size of 0.91 mm 9 0.91 mm 9 1 mm), and the P2 dataset was collected using a Philips ICT 256 Helical CT scanner (372 slices in total with a voxel size of 0.60 mm 9 0.60 mm 9 2 mm). The data collection had local research governance approval.
Segmentation of the CT datasets was undertaken with ITK Snap 2.4.0, 63 using a semi-automatic approach for reconstructing the brachial artery and the proximal sections of the radial and ulnar arteries, and manual segmentation for reconstructing the distal sections of the radial and ulnar arteries. Vessel surfaces were then extracted using VMTK, 3 and smoothed using a combination of VMTK and PTC-Creo Parametric. Specifically, VMTK was used to perform an initial smoothing, from which vessel centrelines were extracted. Variation of maximum inscribed sphere radius was then calculated and averaged in five equisized sections of centreline. Finally, PTC-Creo Parametric swept a circle along each centreline to form the final geometries, where the circle radius was varied smoothly between the averaged maximum inscribed sphere radii in each section of the centreline using the ''loft'' command, which interpolates using a NURBS surface. Images of the resulting arterial geometries for P1 and P2 are shown in Fig. 1  bones in the arm and hand are included for reference in Fig. 1.

Forming an Arterio-Venous Fistulae
Straight veins were digitally anastomosed onto the inner/outer curvature of an arterial bend in P1/P2 using the CAD software PTC Creo Parametric to create brachiocephalic fistulae. Specifically, following Iori et al., 32 inner configurations were designed such that the anastomosis was located on the opposite side of the artery to the bulk of any skewed arterial flow (verified to be the case a posteriori), and outer configurations were designed such that the anastomosis was located on the same side of the artery as the bulk of any skewed arterial flow (also verified to be the case a posteriori). In all cases the vein was anastomosed proximal to the radial/ulnar bifurcation following Chen et al. 18 Moreover, in line with various previous studies, 11,24,31,32,57 the vein had a diameter equal to that of the BAI and made an angle of 35°with the local centre-line of the native artery. In all cases the native arterial geometry was cropped approximately 6 Â 10 À2 m upstream and downstream of the anastomosis. The resulting brachiocephalic AVF configurations for P1 and P2 are shown in Fig. 2.

Blood Flow
For all configurations blood was treated as an incompressible Newtonian fluid. Specifically, blood flow was modeled using the time-dependent incompressible 3D Navier-Stokes equations for a fluid with constant viscosity, which can be written as where l is the viscosity of human blood, q is the density of human blood, u is the three-dimensional blood velocity (vector) field, and p is the pressure. Values of l ¼ 3:5 Â 10 À3 Pa s and q ¼ 1060 kg m À317 were employed.
The assumption of Newtonian rheology ignores the well known ''shear-thinning'' property of blood, which begins to occur at shear rates below 50-100 s À140, 44,62 and is particularly significant below shear rates of 10 s À1 . 40 An a posteriori analysis of results from P1-OUT revealed up to 100, 99 and 93% of the total blood volume was exposed to a shear rate above 10, 100 and 250 s À1 , respectively, at some stage during a single pulse, with analogous figures of 100, 98 and 83%; respectively, obtained for P2-OUT. These significant exposures to high shear, at some stage during a single pulse, suggest widespread suppression of rouleaux formation, which is responsible for shear-thinning and occurs on a time-scale of order 1 s. 51 Consequently, the assumption of Newtonian rheology is considered reasonable. We also note that the assumption has been widely used in previous studies of blood flow in AVF and vascular grafts. 9,11,32,34,38,43,54 Oxygen Transport For all configurations oxygen was treated as a passive scalar dissolved in blood plasma. Specifically, oxygen transport was modeled using the time-dependent 3D advection-diffusion equation, which can be written as where j is the diffusivity of oxygen in human plasma, u is the three-dimensional blood velocity (vector) field, and C is the oxygen concentration. A value of j ¼ 1:2 Â 10 À9 m 2 s À119 was used in this study.
The assumption that oxygen is a passive scalar dissolved in blood plasma ignores the effect of haemoglobin, to which oxygen can bind. However, previous studies suggest that haemoglobin simply acts to augment oxygen transport patterns by a spatially constant factor of approximately two, 37,45 which can be accounted for a posteriori when the results are analysed (see Section ''Wall Shear Stress and Wall Normal Oxygen Flux'').

Blood Flow
A no-slip boundary condition was applied at the AVF wall, which was considered to be rigid. A spacetime varying velocity boundary condition was imposed at the BAI of each AVF configuration. The inflow rate waveform, which had a period of 1 s, was taken from Sigovan et al. 54 Specifically the first 15 Fourier modes were extracted, and the signal was scaled such that the peak Reynolds number at the BAI was 1300, in agreement with Sigovan et al. 54 The resulting timeaveraged Reynolds number at the BAI, Re BAI ; was $750, similar to the time-constant inflow Reynolds number used in Iori et al. 32 All flow was prescribed normal to the BAI inflow plane, and following previous studies, 28,34,43 a Womersley profile was imposed in space. The inflow rate waveforms at the BAI Q BAI for P1 and P2 are shown in Fig. 3. An RCR Windkessel model was applied at the RAO, UAO and VO. Specifically where P RAO ; P UAO and P VO are the spatially averaged pressure waveforms at the RAO, UAO, and VO respectively, Q RAO ; Q UAO and Q VO are the outflow rate waveforms at the RAO, UAO, and VO respectively, and R 1RAO ; Values for the Windkessel parameters were selected using an approach similar to Pant et al. 47 Specifically, they were obtained for P1-IN and P2-IN configurations using an iterative approach, that aimed to minimise where P BAI is the pressure waveform at the BAI, P RS ¼ 130 mmHg and P RD ¼ 80 mmHg are reference systolic and diastolic pressures, respectively, at the BAI taken from Sigovan et al., 54 Q R is a reference outflow rate waveform at the VO taken from Sigovan et al., 54 and the integrals are taken to be over a single pulse, when the solution has become period independent. The iterative process itself employed a 0D lumped parameter model, detailed in Fig. 4, combined with a low-resolution time-dependent 3D model of flow in each configuration. The 0D model represented the lowresolution 3D model in terms of three inductors, with inductances L 1A ; L 2A ; and L 1V as per Fig. 4 as well as three quadratically non-linear resistors with resistances R 1A ; R 2A ; and R 1V as per Fig. 4. The 3D model solved the governing flow equations as per Section ''Blood Flow'' using StarCCM+ v9.06.9 as per Section ''Computational Method'', but with $1=11 of the spatial resolution and 1 / 5 of the temporal resolution. It was assumed that the UAO and RAO has the same peripheral behavior i.e., R 1UAO ¼ R 1RAO ; R 2UAO ¼ R 2RAO and C UAO ¼ C RAO : Hence, both outlets were merged in the 0D model as per Fig. 4, and were parameterized in terms of R 1UAO ; R 2UAO and C UAO alone.
To begin the iterative process, L 1A ; L 2A ; L 1V ; R 1A ; R 2A ; and R 1V were estimated, and the 0D model was used to identify values of the Windkessel parameters that minimized U: These parameters were then used to define boundary conditions for the low-resolution 3D model, which was run until the solution became period independent, and from which improved estimates for L 1A ; L 2A ; L 1V ; R 1A ; R 2A ; and R 1V were then extracted. These were fed back into the 0D model, and the process was repeated until the Windkessel parameters were seen to converge.
Values of the Windkessel parameters obtained for the P1-IN and P2-IN simulations are shown in Table 1.
These parameters were also used for P1-OUT and P2-OUT simulations, respectively.
The resulting average flow splits between the VO and the distal arteries in a temporal window spanning the full pulse period (henceforth referred to as WT), systole from 1.1 to 1.4 s (henceforth referred to as WS), and part of diastole from 1.4 to 1.7 s (henceforth referred to as WD), were $66:34; $60:40; and $78:22 respectively.
The resulting pressure at the BAI varied in the range 82-123 mmHg for P1-IN, and 80-132 mmHg for P2-IN, with similar ranges for P1-OUT and P2-OUT, respectively. The average pressure-drop between the BAI and the VO was $4 mmHg for both P1-IN and P2-IN, with a similar value for both P1-OUT and P2-OUT. The maximum pressure-drop between the BAI and the VO was $11 mmHg for P1-IN, and $9 mmHg for P2-IN, with similar values for P1-OUT and P2-OUT, respectively.
. Schematic illustration of the 0D model used in the Windkessel parameter estimation process. Grey components were used to represent the low-resolution 3D model. We note that a Dean number De defined as ; where j A and D A are the local arterial curvature and arterial diameter respectively at a given point in the artery, never exceeds 320 anywhere in the proximal arterial section of any configuration i.e., it is always well below the critical Dean number of $900 at which multiple unsteady vortices develop in a planar curved tube. 21 Oxygen Transport Following Iori et al., 32 a steady-state spatially-constant oxygen concentration of C NBAI ¼ 1:305Â 10 À1 mol m À3 was applied at the BAI (based on an oxygen partial pressure of 75 mmHg 12 and a Henry's law constant of 1:74 Â 10 À3 mol m À3 mmHg À135 ). Zero boundary-normal oxygen concentration gradients were applied at the VO, RAO, and UAO. Also, for all configurations a steady-state spatially-constant oxygen concentration of C W ¼ 1:044 Â 10 À1 mol m À3 was applied at the arterial and venous walls (based on an oxygen partial pressure of 60 mmHg 12 and a Henry's law constant of 1:74 Â 10 À3 mol m À3 mmHg À135 ).
The imposition of a spatially-constant oxygen concentration at the proximal arterial inlet is based on the assumption that oxygen is ''well mixed'' upstream of the domain. The imposition of a spatially-constant oxygen concentration at the arterial and venous walls is based on the assumption that the walls act as an oxygen sink, readily consuming oxygen. 58

Computational Method
Polyhedral unstructured volume meshes, with prismatic boundary layer meshes adjacent to the arterial wall, were constructed for each configuration. The meshes were refined near the anastomosis. Specifically, elements near the anastomosis had an average size of $6 Â 10 À5 m, expanding progressively to $1:5Â 10 À4 m beyond a distance of $0:045 m from the anastomosis. The prismatic boundary layer meshes were 25 elements thick, with the first element having a thickness of less than 3 Â 10 À6 m; in line with the mesh resolution employed in Coppola and Caro. 19 Each mesh had $11 Â 10 6 elements in total. Results of a mesh independence study are presented in Appendix A.
Solutions for the blood velocity field and the oxygen concentration were obtained using Star-CCM+ v9.06.9 via the following procedure: Each simulation was initialized with zero velocity, zero oxygen concentration, and a pressure of 78 mmHg, and run with the coupled steady solver until every momentum and mass continuity residual fell below 10 À10 : Each steady-state solution was then used as the initial condition for the coupled implicit unsteady solver, which advanced solutions 1 s in time (1 pulse period). For all simulations a timestep of 1 Â 10 À4 s was used.
Each simulation was then advanced a further 1 s (1 pulse period) using the coupled implicit unsteady solver, during which time data were exported for analysis. Results of a period independence study are presented in Appendix B. For all simulations a timestep of 1 Â 10 À4 s was used.
Each simulation was carried out on 64 cores of a Dell AMD Opteron 64-core server with 512 GB RAM, and required approximately 4 weeks to complete.

Unsteady Analysis
Qualitative Insight

Quantitative Insight
Temporal variation of r l ; the WSS magnitude averaged over a ring of arterial section (see Fig. 6), was extracted from each simulation. Plots of r l for each AVF configuration are shown in Fig. 6 together with the Power Spectral Density (PSD) extracted from WT, WS, and WD.
It is clear that the WSS in arterial sections of all AVF configurations exhibit high-frequency fluctuations. This unsteadiness appears to be irrespective of whether the venous section of the AVF is connected to the inner or outer curvature of an arterial bend; converse to the idealized model of Iori et al. 32 which predicted that connecting the venous section of the AVF to the outer curvature of an arterial bend would suppress unsteadiness in the artery. Further observations can also be made. Firstly, within WS during systole, P1-OUT and P2-OUT appear to have relatively more energy in higher frequency modes than P1-IN and P2-IN respectively. This results is in line with the qualitative observations of Section ''Qualitative Insight.'' However, within window WD during diastole, the picture is less clear. Converse to the results from WS, it appears that P1-IN and P2-IN have relatively more energy in higher frequency modes than P1-OUT and P2-OUT respectively, but the differences are less distinct.
Further quantitative assessment of unsteadiness was undertaken using snapshot Proper Orthogonal Decomposition (POD) 28,41,55 of r À r t ; where r is the space-time WSS vector field and r t is the temporal average of the space-time WSS vector field, in an arterial section of each AVF configuration near the anastomosis (see Fig. 6). Specifically, snapshot POD was undertaken in WT, WS, and WD. For WT, 1000 temporal snapshots with a uniform spacing of 0.001 s were used. For WS and WD, 1500 temporal snapshots with a uniform spacing of 0.0002 s were used.
The fraction of POD modes required to reconstruct 96% of the energy content of each r À r t is given in Table 2. Taking this fraction as a proxy for the range of energized frequencies and overall unsteadiness, various observations can be made. Firstly, within WS, during systole, P1-OUT and P2-OUT appear to be relatively more unsteady than P1-IN and P2-IN respectively. This results is in line with the qualitative observations of Section ''Qualitative Insight'' and the quantitative analysis of Section ''Quantitative Insight.'' However, within window WD, during diastole, the picture is less clear. Converse to the results from WS, it appears that P1-IN is relatively more unsteady than P1-OUT. However, P2-IN is relatively only slightly more unsteady than P2-OUT.  Fig. 3), and at 1.60 s, which is during diastole (marked with vertical dashed-dotted lines in Fig. 3).

Summary
The main finding of the above unsteady analysis is that, converse to results from previous more idealized studies, 32 forming an AVF by connecting a vein onto the outer curvature of an arterial bend does not, necessarily, suppress unsteady flow in the artery. Specifically, during systole, outer configurations are actually more unstable than inner configurations, and whilst this appears somewhat reversed in diastole, the situation is less clear, and all configurations exhibit a degree of unsteadiness.
We note that during WS, in systole, the average flow split between the vein and the distal arteries is $60:40: This is significantly more even than the fixed flow split of 80 : 20 prescribed by Iori et al. 32 We can speculate that it is this, more even flow split during WS, that causes the outer configurations to become unstable; noting in particular from Fig. 5 that during WS the outer configurations become unstable in arterial sections distal to the anastomosis, when/where the current simulations exhibit a relatively higher flow rate cf. Iori . Views of P1-OUT (a) and P2-OUT (b) with rings of arterial sections from which r l were extracted for spectral analysis marked in black, and regions of arterial sections from which r À r t were extracted for POD analysis shaded in grey. Plots of r l for each configuration (c), PSD from WT (spanning the full pulse period) of r l for each configuration (d), PSD from WS (spanning systole) of r l for each configuration (e), and PSD from WD (spanning part of diastole) of r l for each configuration (f).

Time-Averaged Analysis Wall Shear Stress and Wall Normal Oxygen Flux
Quad-color maps of WSS À (''pathologically low'' timeaveraged WSS magnitude <0.5 Pa), WSS + (''pathologically high'' WSS magnitude >30 Pa), LWNOF À (''pathologically low'' time-averaged LWNOF <4:275Â 10 À7 mol m À2 s À1 ), and the intersection between regions of WSS À and LWNOF À are shown in Fig. 7, and the corresponding excised and flattened arterial and venous sections from each AVF configuration are shown in Fig. 8. Bar charts showing the percentage area of WSS À , WSS + , LWNOF À , and the intersection between regions of WSS À and LWNOF À in the excised and flattened arterial and venous sections of each configuration are shown in Fig. 9. Note that LWNOF is defined here as À2jn Á $c; where n is the outward facing wall normal and the factor of two accounts for the role of hemoglobin. 37,45 Following Iori et al., 32 the pathologically low WSS threshold of 0.5 Pa was chosen to be at the lowerbound of estimates made by various authors, including Masuda et al. 42 22 who suggest up to 30 Pa. The pathologically low LWNOF threshold was obtained by assuming that the innermost region of the vascular wall (i.e., the intima and media) receive oxygen solely from luminal blood 12,29 and not from the adventitial vasa vasorum. With the further assumptions that this innermost region requires oxygen at a rate of at least 8:55 Â 10 À3 mol m À3 s À1 in order to avoid hypoxia (based on measurements of oxygen consumption in smooth muscle cells of dog femoral arteries 48 ), and that it has a thickness of 5 Â 10 À5 m (based on measurements in dog femoral arteries 15 ), one obtains the threshold of 4:275 Â 10 À7 mol m À2 s À1 .
It can be seen that for the P1-OUT and P2-OUT configurations the artery is exposed to significantly more WSS À than for the P1-IN and P2-IN configurations, respectively. Also, for the P1-OUT and P2-OUT configurations, a triangular shaped region of LWNOF À is present in the artery opposite the anastomosis, where stenosis in AVF have previously been observed, 56    results of previous more idealized studies, 32 and add further support to the assertion that connecting a vein to the outer curvature of an arterial bend may suppress IH if it is caused by low WSS and/or low LNWOF (leading to wall hypoxia). Finally, for the P1-OUT and P2-OUT configuration the vein is exposed to significantly more WSS + than for the P1-IN and P2-IN configurations. Previous more idealized studies, where the vein and artery had the same diameter, 32 exhibited the same trend. However, the difference in this study between inner and outer curvature configurations was far more marked.

Summary
The main finding of the above time-averaged analysis is that, in agreement with results from previous more idealized studies, forming an AVF by connecting a vein onto the inner curvature of an arterial bend will suppress exposure to regions of low WSS and hypoxia in the artery. Results also show that forming an AVF by connecting a vein onto the inner curvature of an arterial bend will significantly reduce exposure to high WSS in the vein.

CONCLUSIONS AND FUTURE WORK
We have extended previous work by Iori et al., 32 which investigated the effect of planar arterial curvature on flow and oxygen transport patterns in idealized AVF with non-pulsatile inflow conditions, to more realistic non-planar brachiocephalic AVF configurations with pulsatile inflow conditions. Results from our more realistic simulations show that forming an AVF by connecting a vein onto the outer curvature of an arterial bend does not, necessarily, suppress unsteady flow in the artery. Specifically, during systole, outer configurations are actually more unstable than inner configurations, and whilst this is somewhat reversed in diastole, the situation is less clear, and all configurations exhibit a degree of unsteadiness. This finding is converse to results from the previous more idealized study of Iori et al. 32 However, results also show that forming an AVF by connecting a vein onto the inner curvature of an arterial bend will suppress exposure to regions of low WSS and hypoxia in the artery. This finding is in agreement with results from the previous more idealized study of Iori et al. 32 Finally, results show that forming an AVF by connecting a vein onto the inner curvature of an arterial bend will significantly reduce exposure to high WSS in the vein. The results are important, as they suggest that in realistic scenarios arterial curvature can be leveraged to reduce exposure to excessively low/high levels of WSS and regions of hypoxia in AVF. This may in turn reduce rates of IH and hence AVF failure.
Future work should address limitations of the current modeling approach. In particular, the effect of wall compliance should be considered. Previous studies have shown time-averaged WSS patterns to be qualitatively similar between models with rigid/compliant walls. 43 However, the effect of compliance on unsteady dynamics has yet to be elucidated. Future work should also investigate the effect of geometric variations across an extended cohort of patient specific arterial datasets, as well as the effect of anastomotic connection angle, vein-to-artery diameter ratio, venous curvature, and variations in morphology induced by arm movement and bending. Previous studies suggest that neck movement can significantly alter vascular morphology in a patient specific fashion. 4,5 It is therefore possible that the morphology of an AVF changes when a patient moves or bends their arm. Finally, future  work should development both bench-top and in-vivo pre-clinical experiments, in order to validate the computational results.

Polyhedral Unstructured Volume Mesh
Grid convergence of the polyhedral unstructured volume mesh was assessed by running simulations on seven successively finer discretizations of the P1-IN configuration using Star-CCM+ v9.06.9. All discretizations had a fixed prismatic boundary layer mesh that was 25 elements thick, with the first element having a thickness of $2:6 Â 10 À6 m; in line with the boundary layer mesh resolution employed in Iori et al. 32 For each mesh, the element size near the anastomosis D , and the overall number of mesh elements N e ; are given in Table A.3. Each simulation had steady inflow boundary conditions. Specifically, the inflow Reynolds number at the BAI was fixed at 1, 300, the flow split between the VO and the distal arteries was 45 : 55, and the flow split between the UAO and the RAO was 50 : 50. These conditions correspond to peak inflow in the associated pulsatile simulation. The oxygen boundary conditions were the same as those used for the pulsatile simulation. Each simulation was initialized by converging to a steady-state solution such that all the residuals had an absolute value smaller than 1 Â 10 À9 : Each steadystate solution was then used as the initial condition for the coupled implicit unsteady solver, which advanced  solutions 0.2 s, and then a further 0.07 s over which data was extracted for analysis. All unsteady simulation used a time step of 2:2 Â 10 À4 s.
Time-averaged WSS magnitude r t and time-averaged LWNOF f t were extracted from each simulation. Distributions of r t and f t are shown in Figs. A.1 and A.2 respectively. Space-time-averaged WSS magnitude r ts and space-time-averaged LWNOF f ts were also calculated for each simulation, and used to obtain percentage errors E r and E f defined as and where r 6ts and f 6ts are values of r ts and r ts respectively from mesh M6 (the finest discretization). Plots showing how E r and E f vary with N e are shown in Fig. A.3.

Prismatic Boundary Layer Mesh
Grid convergence of the prismatic boundary layer mesh was assessed by running simulations on versions of Mesh 5 with either increased or decreased boundary layer resolution. Specifically two meshes, henceforth referred to as Mesh 5þ and Mesh 5À; were made. Each mesh had the same overall boundary layer mesh thickness as Mesh 5, but a different number of boundary layer mesh elements N l in the wall normal direction (see Table A

APPENDIX B: PERIOD INDEPENDENCE
The P1-IN simulation was started from a steadystate solution as per Section ''Computational Method'' and advanced 3 s in time (3 pulse periods) using the coupled implicit unsteady solver and a timestep of 1 Â 10 À4 s. Period convergence was assessed by comparing solution data from t 2 ½1; 2 s, hence forth referred to as T1, with solution data from t 2 ½2; 3 s, hence forth referred to as T2.
Distributions of r t over T1 and T2 are shown in Fig. B.1. Visually the distributions appear very similar. Quantitatively, spatial averages of r t over T1 and T2 were found to be within 0:08%: Distributions of f t over T1 and T2 are shown in Fig. B.2. Once again, visually the distributions appear very similar. Quantitatively, spatial averages of f t over T1 and T2 were found to be within 0:03%:   Unsteady content was compared between T1 and T2 using snapshot POD of r À r t in an arterial section of each AVF configuration near the anastomosis (see Fig. 6). Specifically, 1000 temporal snapshots with a uniform spacing of 0.001 s were used from each of T1 and T2. Plots of the first, fourth, and seventh temporal POD modes, denoted a 1 , a 4 , and a 7 respectively, extracted from T1 and T2, are shown in Fig. B.3. Visually the modes appear very similar.
Based on the results presented above, data from t 2 ½1; 2 s were exported for analysis from simulations in P1-IN, P1-OUT, P2-IN, and P2-OUT.