Suppressing Unsteady Flow in Arterio-Venous Fistulae

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, a large proportion of AVF do not mature, and 1 hence fail, as a result of various pathologies such as Intimal Hyperplasia (IH). Unphysiological flow patterns, including high-frequency flow unsteadiness, associated with the unnatural and often complex geometries of AVF are believed to be implicated in the development of IH. In the present study we employ a Mesh Adaptive Direct Search optimisation framework, computational fluid dynamics simulations, and a new cost-function, to design a novel non-planar AVF configuration that can suppress high-frequency unsteady flow. A prototype device for holding an AVF in the optimal configuration is then fabricated, and proof-of-concept is demonstrated in a porcine model. Results constitute the first use of numerical optimisation to design a device for suppressing potentially pathological high-frequency flow unsteadiness in AVF.


Introduction
Haemodialysis is a common form of renal replacement therapy for patients who suffer with End-Stage Renal Disease (ESRD) [1,2]. Efficient haemodialysis requires highquality vascular access, allowing for the rapid removal of a patient's blood, which passes through a dialyser before being returned to the body. The 'gold-standard' method of vascular access is via an established Arterio-Venous Fistula (AVF) [3] created by surgically connecting an artery and vein, usually in the patient's arm.
The connection or 'anastomosis' can take a number of configurations. However, an end-to-side layout, where the end of the vein is connected to the side of the artery, has emerged as the clinically preferred configuration [4].
After an AVF is formed, a low-resistance return pathway is created from the arterial 2 system to the venous system, with a consequential increase in blood flow through the venous section of the AVF. In an ideal scenario these changes stimulate vascular dilation and re-modeling, such that the AVF 'matures', and can be used for haemodialysis [5]. Unfortunately, however, up to 60% of AVF do not mature as anticipated [6] resulting in unacceptable patient morbidity and a significant burden on healthcare expenditure. A major cause of AVF failure is the formation of stenotic regions due to development of Intimal Hyperplasia (IH) in the juxta-anastomotic area and the vein [7,8,5]. Whilst the exact mechanisms underlying development of IH in AVF are unknown, there is considerable evidence to suggest that their inherently unphysiological flow patterns play an important role [9,8,5]. In particular regions of low Wall Shear Stress (WSS) [10,11,12,9,13], low lumen-to-wall oxygen flux (leading to wall hypoxia) [14,15,16,17], and high-frequency flow unsteadiness [18,19,20,21,22] have all been individually implicated in the initiation and development of IH, although it is of course possible that these factors are related, and a combination is responsible.
The present study is based on the hypothesis that high-frequency flow unsteadiness plays an important role in development of IH. The genesis and nature of unsteadiness in pipe flow has been studied for over a century. The seminal work of Reynolds identified a critical Reynolds number for flow in a straight pipe, beyond which flow transitions to turbulence [23] -although only now is the entire process of transition being fully understood [24]. Subsequently, the seminal work of Dean identified a critical Dean number beyond which multiple unsteady vortices develop in a planar curved pipe [25]. More recent work, undertaken specifically in the context of vascular fluid mechanics, includes studies of stenotic flows [26], pulsatile flow in curved pipes [27,28,29,30], flow in helical pipes [31,32,33,34], and flow at junctions [35,36,37].
In terms of simulating flow in AVF configurations, there have been a range of previous Computational Fluid Dynamics (CFD) studies, including 1D simulations [38,39,40], 2D simulations [41], 3D simulations with rigid walls in idealised geometries [42,43,44,45,46,47], 3D simulations with rigid walls in more realistic geometries [48,49,50,51,52,53,54,55], and simulations with distensible walls [56]. Of these studies, several investigated flow unsteadiness. In particular, Ene-Iordache et al. [46,45] undertook CFD simulations of pulsatile blood flow in idealised end-to-end and end-to-side AVF configurations with various anastomotic angles. According to their results, AVF formed by connecting a straight vein to a straight artery with an acute anastomotic angle of 30 degrees are optimal, because they reduce exposure to 'disturbed flow', the definition of which included an unsteady component. This finding appears to have inspired development of an externally mounted device to help maintain an acute anastomotic angle [57]. More recently, CFD simulations in idealised configurations with an acute anastomotic angle of 35 degrees were undertaken by Browne et al. [43,44]. Results exhibited unsteady flow in both the artery and vein, as did results from simulations in four more realistic configurations undertaken by Bozzetto et al. [48]. Finally, Iori et al. [42] used CFD simulations in idealised configurations to demonstrate that forming an AVF by connecting a vein onto the outside of an arterial bend could suppress high-frequency unsteadiness inside the artery. However, flow within the venous section, where pathology is also known to develop, typically remained unstable, and a recent follow-on study by Grechy et al. [55], in more realistic configurations, suggested that unsteady flow is present throughout such configurations at certain points in the pulse cycle.
In the present study we employ a Mesh Adaptive Direct Search (MADS) optimisation framework, CFD simulations, and a new cost function, to design a novel non-planar 4 AVF configuration that suppresses high-frequency flow unsteadiness throughout the AVF, in both the arterial and venous sections. A prototype device for holding an AVF in the optimal configuration is then fabricated, and proof-of-concept is demonstrated in a porcine model. Results constitute the first use of numerical optimisation to design a device for suppressing potentially pathological high-frequency flow unsteadiness in AVF.

Overview
In this section we present the optimisation framework, which is based on the MADS method proposed by Audet and Dennis [58,59], and uses a Kriging based surrogate model [60,61] with a zeroth-order regression method, and a Gaussian correlation function, to predict values of the objective function. Despite limitations around convergence rates cf. gradient based methods [62], and no strict guarantee of a global optimum, similar methods have been used to design a range medical devices and surgical procedures [63,64,65,66], including Y-shaped grafts for the Fontan procedure [64,65], and cardiovascular stents [66,67]. We note that the cost function used here is new, as is application of the MADS method to design of a non-planar vascular configuration. 5

Parameterisation
The artery was constructed by sweeping a circular section with diameter D a = 5 × 10 −3 m along a plane-arc centreline x a defined as which runs from the Proximal Arterial Inlet (PAI) to the Distal Arterial Outlet 30,30] is the arterial arc length. A priori prescription of a curved arterial section is based on the recent findings of Iori et al. [42], who showed that arterial curvature can have an significant impact on unsteadiness in AVF. The vein was constructed by sweeping a circular section with diameter D a along a twopiece cubic hermitian spline x v defined as where s v ∈ [0, 1] is the venous arc length, is where the arterial and venous centrelines intersect, is a mid-point on the venous centreline, 6 is the location of the Distal Venous Outlet (DVO), is a tangent where the arterial and venous centrelines intersect, is a tangent at the DVO, and finally is a tangent at a mid-point on the venous centreline, defined such that the venous centreline is C2 continuous.
The resulting AVF configuration is described by a set of six parameters: x p1 , y p1 and z p1 , which define the location of a mid-point on the venous centreline, θ and φ, which define the angle between the vein and artery, and ξ, which defines the location of the DVO. The parameters were constrained such that x p 1 ∈ [−7D a , 7D a ], Additionally, the parameters were constrained such that the configuration was not self-intersecting, the maximum venous curvature was less than 1.8/D a , the venous length was less than 20D a , the whole AVF configuration fitted within an infinitely long cylinder, aligned with the x-axis, and centred on the origin, with a radius of 7.5D a . The parameterisation and constraints were chosen to ensure that the resulting geometry had clinical/anatomical relevance. Specifically, the artery to vein diameter ratio of 1 : 1 is in line with various measurements made post AVF surgery [38,68], 7 the maximum venous curvature of 1.8/D a is in line with measurements of maximum curvature in aneurysmatic internal carotid arteries [69], and the restriction to fit within a cylinder of radius 7.5D a is in line with measurements of average human arm radius [70].
A schematic illustration of the arterial and venous centrelines is shown in Fig. 1.
We note that in general the parameterised AVF configuration is non-planar. Figure 1: Schematic illustration of the arterial centreline x a and the venous centreline x v , along with associated parameterisation in terms of θ, φ, x p 1 , y p 1 , z p 1 and ξ. Flow enters at the PAI and then divides at the anastomosis, exiting both the DVO and DAO. In general the parameterised AVF configuration is non-planar. 8

Definition
We introduce a new cost function that indicates whether a particular AVF configuration is likely to induce unsteady flow patterns. To evaluate the cost function for a particular configuration we first find a steady-state solution to the continuity equation ∇ · u = 0, (9) and the steady-state incompressible Navier-Stokes equations for a Newtonian fluid where µ is the viscosity of human blood, ρ is the density of human blood, u is the three-dimensional blood velocity (vector) field, and p is the pressure. For all configurations, values of µ = 3.5 × 10 −3 Pa · s and ρ = 1060kg · m −3 were employed.
Also, for all configurations, steady-state parabolic velocity profiles were applied at the PAI with an average velocity of 0.528m · s −1 into the configuration (equivalent to ∼ 622 ml/min), and at the DAO with an average velocity of 0.106m · s −1 out of the configuration. This resulted in a Reynolds number of 800 at the PAI, and an 80:20 flow split between the DVO and DAO for each configuration [38,71]. Finally, a constant (and arbitrary) pressure was applied at the DVO of each configuration, and a zero-velocity no-slip condition was applied at the arterial and venous walls, which were assumed to be rigid.
The resulting steady-state solution is then used as the initial condition for Eq. 9 and the unsteady incompressible Navier-Stokes equations for a Newtonian fluid with the same parameter values and boundary conditions employed for the steadystate problem. Any subsequent divergence of the solution from the steady-state initial condition can be used as a proxy for unsteadiness in the particular configuration.
Specifically, based on qualitative prior experience of when and how fundamentally unsteady results begin to diverge from a steady-state initial condition, we define our cost function F as the maximum percentage surface area of the configuration that has at least a 10% deviation in WSS magnitude from the steady-state solution during the interval 0.04732s ≤ t ≤ 0.09464s, thus where Γ is the surface area of the given configuration and where σ S (x) is the WSS magnitude from the steady-state solution, and σ U (x, t) is the WSS magnitude from the unsteady solution.
We note that the cost function F is defined in the context of non-pulsatile boundary conditions. This is clearly an approximation, since in reality blood flow is pulsatile. nificant high-frequency unsteadiness even with non-pulsatile inflow conditions; and it is this high-frequency unsteadiness that we are looking to suppress. Moreover, performance of the optimal configuration will be assessed a posteriori with a fully pulsatile simulation in Section 3. We also note that the cost function F , as defined, is relatively cheap to evaluate cf. one based on a full unsteady simulation. This makes it feasible to optimise in the six-dimensional parameter space.

Computational Method
Star-CCM+ v9.06.9 (CD-Adapco, Melville, NY USA) was used to obtain F for a given AVF configuration. Specifically, 1000 iterations of the segregated steady-state solver were used to solve Eq. 9 with Eq. 10, and the segregated implicit unsteady solver with a time-step of 4.735 × 10 −4 s was used to solve Eq. 9 with Eq. 11. Polyhedral unstructured volume meshes, with prismatic boundary layer meshes adjacent to the wall, were constructed for each configuration. The meshes were refined near the anastomosis. Specifically, elements near the anastomosis had an average size of 1.5×10 −4 m, expanding progressively to 3×10 −4 m beyond a distance of ∼ 2.5×10 −2 m from the anastomosis. The prismatic boundary layer meshes were 10 elements thick, with the first element having a thickness of 3.23 × 10 −5 m. Each mesh had ∼ 1 × 10 6 elements in total. We note that this mesh resolution is lower than that employed by Iori et al. [42]. However, it is similar to that employed in other studies of unsteady flow in AVF [43,51]. Moreover, performance of the optimal configuration will be assessed with a higher-resolution simulation in section 3.

Initial Seeding
The Kriging-based surrogate model of F was initially seeded with data at a Latin Hypercube Sampling (LHS) of the six-dimensional parameter space [60,61]. Specifically, 600 initial configurations that satisfied whereσ D (x) is the time-averaged σ D (x, t) over the period 0.04732 − 0.09464s. The richness of the parameter space is apparent. The initial Kriging-based surrogate model of F was calculated using the Matlab Dace toolbox [72].
Initial Sampling Configuration

Mesh Adaptive Direct Search Method
Having initialised the Kriging-based surrogate model of F , the MADS method was used to search for a feasible configuration that minimised F . Values of x p1 , y p1 , z p1 , θ, φ, ξ, and F for the 29 feasible configurations at which F was evaluated as the MADS method proceeded through so-called SEARCH and POLL steps [58,59] are shown in Fig. 4. The Matlab Dace toolbox [72] was employed to re-calculate the Kriging-based surrogate model of F , as required, as the MADS method proceeded.
The MADS method self-terminated when the mesh size parameter associated with the SEARCH step was refined more than three times, as per Yang et al. [65].  Figure 4: Values of x p1 , y p1 , z p1 , θ, φ,ξ, and F for the 29 feasible configurations at which F was evaluated as the MADS method proceeded through so-called SEARCH and POLL steps [58,59]. The solid black bars correspond to feasible configurations predicted to have the lowest values of F by the Kriging-based surrogate model during the SEARCH steps. The solid grey bars correspond to feasible configurations for which the Kriging-based surrogate model has the highest predicted error during the SEARCH steps. The diagonal striped bars correspond to feasible configurations at which F was evaluated during the POLL steps. Vertical dashed lines correspond to the optimal configuration. Note that it was identified during a POLL step. 16

Results and Discussion
The lowest value of F obtained from the optimisation process was 0.0062. The associated parameter values are reported in Table 1. The resulting AVF configuration is shown in Fig 5, and an STL file is provided as Electronic Supplementary Material.
We note that the optimal configuration features an anastomosis on the outer curvature of the arterial bend, inline with the results of Iori et al. [42] who showed that connecting the vein onto the outside of an arterial bend could suppress unsteady flow in the artery. We also note that the configuration features an inherently non-planar looped venous section, and a relatively obtuse anastomotic connection angle -in contrast to previous designs for suppression of 'disturbed flow', which were planar and had a relatively acute 30 degree anastomotic connection angle [45,46]. Finally we note that, qualitatively, our optimal non-planar configuration bears a resemblance to the 'curved' (as opposed to the 'straight') AVF configuration of Krishnamoorthy et al. [53], which increased AVF patency rates in a pig model.

Parameter
Optimal Value x y z Figure 5: Optimal configuration coloured byf (x) (white whenf (x) = 0 and black whenf (x) = 1). The black areas are too small to be visible. Salient aspects of the optimal configuration include an anastomosis on the outer curvature of the arterial bend, a looped non-planar venous section, and a relatively obtuse anastomotic connection angle. An .stl file definition of the optimal configuration is included as electronic supplementary material.
Having identified an optimal configuration, it is useful to assess the sensitivity of F to each of the six geometric parameters. Fig. 4 clearly shows that θ and φ vary little during the MADS process, suggesting they play an important role, and must remain close to their optimal values. However, ξ varies significantly, indicating it is less critical. Further insight can be gained from Fig. 6, which plots the median, interquartile range, and range, of parameters associated with configurations that achieve θ and φ appear to be the most important parameters, with ξ in particular having little effect. In summary, the optimal configuration appears particularly sensitive to perturbations of θ and φ, which define the angle made between the vein and the curved artery. Non-planarity also appears to be important. However, the location of the distal end of the vein appears to be less important.  Figure 6: Box plots of the median (horizontal black line), inter-quartile range (grey box), and range (whiskers), of parameters associated with configurations that achieve F < 0.1. Black stars indicate the parameter values of the optimal set. It is clear that θ and φ play an important role, and must remain close to their optimal values. However, ξ varies significantly, indicating it is less critical. The plots also suggest that z p1 should remain close to its optimal (non-zero) value, signalling the importance of non-planarity.

Overview
In this section we undertake high-resolution unsteady pulsatile simulations of flow within the optimal AVF configuration to ascertain whether high-frequency unsteadiness remains suppressed when the inflow is pulsatile.
In all cases a space-time varying velocity boundary condition was imposed at the PAI.
The inflow waveform, which had a period of 1s, was obtained from [49]. Specifically the first 15 Fourier modes were extracted, and the signal was scaled such that the peak Reynolds number at the PAI was 1300, in agreement with [49]. The resulting time-averaged Reynolds number at the PAI was ∼ 750, similar to the time-constant PAI Reynolds number of 800 used during the optimisation process. The resulting Womersley number was 6.9. All flow was prescribed normal to the PAI inflow plane, and in line with previous studies [56,54,77], a Womersley profile was imposed in space. The PAI inflow rate Q P AI as a function of time is shown in Fig. 8.

22
In all cases an RCR Windkessel model was applied at the DAO and DVO. Specifically (16) and where where P P AI is the pressure waveform at the PAI, P RS = 130mmHg and P RD = 80mmHg are reference systolic and diastolic pressures, respectively, at the PAI, Q R is a reference outflow rate waveform at the VO taken from Sigovan et al. [49], and the integrals are taken to be over a single pulse, when the solution has become period independent.
The iterative process itself combined a 0D lumped parameter model, with a lowresolution time-dependent 3D flow model. The 0D model represented the lowresolution 3D model in terms of three inductors, as well as three quadratically nonlinear resistors. The 3D model solved Eq. 9 and Eq. 11 using StarCCM+ v9.06.9.
To begin the iterative process, the 0D lumped parameters were estimated using the Nelder-Mead simplex algorithm implemented in Matlab, and the 0D model was used to identify values of the Windkessel parameters that minimised Φ. 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 the 0D parameters were then extracted. These were fed back into the 0D model, and the process was repeated until the Windkessel parameters were seen to converge.
The values of the resulting RCR parameters are given in Tab. 2.  Fig. 8.

24
Finally, in all cases a no-slip zero-velocity boundary condition was applied at the AVF wall, which was assumed to be rigid.

Computational Method
Star-CCM+ v9.06.9 (CD-Adapco, Melville, NY USA) was used to perform all simulations. To begin, zero flow and a reference pressure of 80mmHg were set as initial conditions for the coupled steady-state solver, which was used to solve Eq. 9 with Eq. 10, converging every momentum and mass continuity residual below 10 −10 . The resulting steady-state solution was then set as the initial condition for the coupled implicit unsteady solver, which was used to solve Eq. 9 with Eq. 11. Specifically, the had ∼ 12 × 10 6 elements in total. This mesh resolution is similar to that used in [42]. Fig. 9 shows snapshots of velocity magnitude u(x) on the arterial symmetry plane, and on a three-dimensional surface that tracks the venous centreline, at four points in T 1 for the Newtonian simulation. Fig. 10 shows snapshots of plane-normal vorticity ω(x) on the arterial symmetry plane, and on a three-dimensional surface that tracks the venous centreline, at four points in T 1 for the Newtonian simulation. Fig. 11 shows streamlines, originating from the PAI, coloured by velocity magnitude, at four high-frequency flow unsteadiness in previous studies [42,43,44,51,55], are absent. u(x) (ms −1 ) Figure 11: Streamlines, originating from the PAI, coloured by u(x), at four points in T 1 for the Newtonian simulation.

Snapshot Proper Orthogonal Decomposition
A quantitative assessment of unsteadiness was made via snapshot Proper Orthogonal Decomposition (POD) [79,80,77] of WSS fluctuations σ f (x, t) = σ(x, t) −σ(x) for the Newtonian simulation, where σ(x, t) is the WSS field, andσ(x) is the timeaveraged WSS field. Specifically, snapshot POD of σ f (x, t) was undertaken in T 1 using 1000 temporal snapshots with a uniform spacing of 0.001s. Fig. 12 shows 29 the Power Spectral Density (PSD) of a 1 (t) and a 5 (t), the first and fifth temporal POD modes respectively for the Newtonian simulation. The second through fourth modes exhibited similar behaviour and are hence omitted for brevity. All significant energy is below 10Hz, and on comparison with Fig. 13, which shows the PSD of Q P AI , it is reasonable to attribute the majority of this low-frequency content to the pulse waveform. There is certainly no significant energy content around 60Hz as per Iori et al. [42] and Grechy et al. [55], nor above 100Hz as per Bozzetto et al. [48], Browne et al. [43,44] and Iori et al. [42]. Consequently, the results demonstrate that the optimal configuration can suppress potentially pathological high-frequency flow unsteadiness even when the inflow is pulsatile.   Figure 13: PSD of Q PAI on a semi-log scale (a) and a log-log scale (b) for the Newtonian simulation.

Pressure Drops
The time-averaged pressure drop between the DVO and PAI for the Newtonian simulation was 2.78mmHg, while the time-averaged pressure drop between DAO and PAI for the Newtonian simulation was 1.013mmHg. These values are substantially lower than those obtained in previous studies that observed unsteady flow in the vein [51,43,44]. This is inline with the expectation that more stable flow is associated with a reduced pressure drop. In this section we present results from a proof-of-concept study. Specifically, a prototype device is developed to hold an AVF in the optimal configuration, and deployed in a porcine model.

Device Design and Manufacturing
A CAD model of an extravascular device for holding an AVF in the optimal configuration was developed (see Fig. 14). Specifically, the device consisted of a solid central section for structural integrity, with attached 'guttering' sections to accommodate both the arterial and venous segments of an AVF. The device was designed to be implanted ('slotted' onto the the AVF) after forming the anastomosis between the artery and the vein.
Going forwards, it is envisaged that the device be fabricated from a bio-compatible silicone such as Nusil Med-4830, or a bioresorbable hydrogel. However, for this non-recovery proof-of-concept the device was fabricated in Objet TangoBlackPlus FLX980 using a Connex 350 3D printer. Specifically, a device was fabricated with a 5 × 10 −3 mm internal guttering diameter, along with a counterpart device that was mirror-symmetric in the plane of the arterial curve (see Fig. 15).  A mid-line incision was made in the neck and the jugular and carotid vessels were dissected out on both the left and right sides of the neck (with left-right laterality defined in respect to the animal). Bilaterally, the jugular vein was then cut distally 35 and mobilised, with all visible tributaries ligated, and a slit arteriotomy was made in the common carotid artery before a microsurgical anastomosis (10/0 nylon) was made between the mobilised jugular vein and the common carotid artery to create an endto-side AVF. Following creation of the AVF, the right-handed and left-handed devices were implanted to shape the left-side and right-side AVF, respectively. Given the anatomical limitations of the surgical field a short (∼ 5×10 −3 m) portion of the device in the region of the DVO had to be trimmed on both sides to permit implantation.
Personal and project licences for the study were granted by the United Kingdom Home Office.

Stethoscope Recording
Stethoscope spectra with high frequency peaks have been successfully correlated with the presence of stenosis in vascular connections [81,82,83,84]. It is well known that stenotic vessels produce highly unsteady flow [77,26], therefore the stethoscope spectra were adopted as a proxy for unsteadiness. Digital stethoscope recordings were made directly over the anastomosis both before (left-side AVF only) and after (leftside and right-side AVF) device implantation. Specifically, recordings were made using using an electronic 3M Littmann 3200 series stethoscope for approximately 20s with a sample rate of 4,000 Hz, and three independent but otherwise identical recordings were made at each location.

Results and Discussion
Figs. 16 17, 18 show the left-side AVF before implanting the device, the left-side AVF after implanting the device, and the right-side AVF after implanting the device, 36 respective. On both the left-and right-hand sides the device held itself in place via the guttering, and visually formed the AVF into the desired configuration. PAI DVO DAO Figure 16: Left-side AVF before implanting the device. The PAI is the proximal common carotid artery, the DAO is the distal common carotid artery, and the DVO is the jugular vein. 37 DVO PAI DAO Figure 17: Left-side AVF after implanting the device. The PAI is the proximal common carotid artery, the DAO is the distal common carotid artery, and the DVO is the jugular vein. DVO PAI DAO Figure 18: Right-side AVF after implanting the device. The PAI is the proximal common carotid artery, the DAO is the distal common carotid artery, and the DVO is the jugular vein.
Figs. 19,20,21 show PSD produced via a Burg Autoregressive Model of the stethoscope recordings from the left-side AVF before implanting the device, the left-side AVF after implanting the device, and the right-side AVF after implanting the device, respectively. From Fig. 19 it is clear that that all three recordings for the left-side AVF, before implanting the device, exhibit a peak at ∼ 600Hz. This is in a similar range to peaks found by previous computational [42] and experimental [84] studies of flow in AVF configurations, and is presumed to be caused by the AVF junction (although to be sure one should compare with PSD of the native vasculature, prior to AVF creation, which were not available). From Fig. 20 it is clear that that these peaks are suppressed after the device is implanted. Moreover, from Fig. 21 it is clear that the PSD for the right-side AVF, after implanting the device, look visually similar to their left-side counterparts; without a significant peak at ∼ 600Hz. However, PSD for right-side AVF, before implanting the device, were not available for direct comparison. The results indicate that the device is acting to suppress high-frequency unsteady flow in the AVF.  Figure 19: PSD produced via a Burg Autoregressive Model of three independent but otherwise identical stethoscope recordings from the left-side AVF before implanting the device.  Figure 21: PSD produced via a Burg Autoregressive Model of three independent but otherwise identical stethoscope recordings from the right-side AVF after implanting the device.
In summary, results from the proof-of-concept indicate that the device can suppress high-frequency unsteady flow in an AVF. However, various questions remain, which should be addressed by future studies. In particular, i.) does the device itself play an important viscoelastic structural role in suppression of high-frequency flow unsteadiness, ii.) is suppression of high-frequency flow unsteadiness sustained across a series of porcine models, and in time over a period of months if a bio-compatible device is left in-place inside a porcine model, iii.) will the external nature of the device interfere with vascular mass transport, or hinder favourable remodelling, and thus unintentionally cause AVF pathology, and finally, iv.) can such a device actually reduce the prevalence of IH, and hence AVF failure rates.

Conclusions
High-frequency flow unsteadiness has been associated with development of vascular pathology in AVF [18,19,20,21,22]. In this study we employed a MADS optimisation framework, and CFD simulations, to design a non-planar AVF configuration that suppresses high-frequency flow unsteadiness. Salient aspects of the novel configuration include an anastomosis on the outer curvature of an arterial bend, and a looped non-planar venous section with a relatively obtuse anastomotic connection angle. A prototype device for holding an AVF in the optimal configuration was then fabricated, and proof-of-concept was demonstrated in a porcine model. Results constitute the first use of numerical optimisation to design a device for suppressing potentially pathological high-frequency flow unsteadiness in AVF.

Electronic Supplementary Material
An STL file of the optimal AVF configuration is provided as Electronic Supplementary Material.