Arterio-Venous Fistulae (AVF) are regarded as the “gold standard” method of vascular access for patients with end-stage renal disease who require haemodialysis. However, a large proportion of AVF do not mature, and 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.
I. 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 high-quality 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 a 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 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–12,9,13 low lumen-to-wall oxygen flux (leading to wall hypoxia)14–17 and high-frequency flow unsteadiness18–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 have 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 turbulence23—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–30 flow in helical pipes,31–34 and flow at junctions.35–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–40 2D simulations,41 3D simulations with rigid walls in idealised geometries,42–47 3D simulations with rigid walls in more realistic geometries,48–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° 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° 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 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.
II. OPTIMISATION
A. Overview
In this section, we present the optimisation framework, which is based on the MADS method proposed by Audet and Dennis58,59 and uses a Kriging-based surrogate model60,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 of medical devices and surgical procedures,63–66 including Y-shaped grafts for the Fontan procedure64,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 a non-planar vascular configuration.
B. Parameterisation
The artery was constructed by sweeping a circular section with diameter along a plane-arc centreline xa, defined as
which runs from the Proximal Arterial Inlet (PAI) to the Distal Arterial Outlet (DAO), where 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 a significant impact on unsteadiness in AVF. The vein was constructed by sweeping a circular section with diameter Da along a two-piece cubic hermitian spline defined as
where is the venous arc length,
is where the arterial and venous centrelines intersect,
is the mid-point on the venous centreline,
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: xp1, yp1, and zp1, which define the location of the 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 , , , , , and . Additionally, the parameters were constrained such that the configuration was not self-intersecting, the maximum venous curvature was less than 1.8/Da, the venous length was less than 20Da, and the whole AVF configuration fitted within an infinitely long cylinder, aligned with the x-axis, and centered on the origin, with a radius of 7.5Da. 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 the maximum venous curvature of 1.8/Da 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.5Da 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.
Schematic illustration of the arterial centreline xa and the venous centreline , along with associated parameterisation in terms of , ϕ, , , , 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.
Schematic illustration of the arterial centreline xa and the venous centreline , along with associated parameterisation in terms of , ϕ, , , , 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.
C. Cost function
1. 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
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 Pa s and were employed. Also, for all configurations, steady-state parabolic velocity profiles were applied at the PAI with an average velocity of into the configuration (equivalent to ) and at the DAO with an average velocity of 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 steady-state 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 , 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. However, Iori et al.42 have previously demonstrated that AVF can exhibit significant 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 Sec. III. 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.
2. 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 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 , expanding progressively to beyond a distance of from the anastomosis. The prismatic boundary layer meshes were 10 elements thick, with the first element having a thickness of . Each mesh had 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 Sec. III.
D. Optimisation framework
1. 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 , , , , , and were identified via LHS. Of these 600 configurations, 509 were rejected because they were self-intersecting, or they violated the other constraints on maximum venous curvature, venous length, and overall AVF size. Values of xp1, yp1, zp1, , ϕ, and ξ for the remaining 91 valid AVF configurations, along with calculated values of F, are shown in Fig. 2. A sample of the 91 valid AVF configurations is shown in Fig. 3, where the surfaces are coloured by , defined as
where is the time-averaged σD(x, t) over the period 0.04 732–0.09 464 s. The richness of the parameter space is apparent. The initial Kriging-based surrogate model of F was calculated using the Matlab Dace toolbox.72
Values of xp1, yp1, zp1, , ϕ, ξ, and F for the 91 AVF configurations used to seed the Kriging-based surrogate model of F.
Values of xp1, yp1, zp1, , ϕ, ξ, and F for the 91 AVF configurations used to seed the Kriging-based surrogate model of F.
A sample of the 91 valid AVF configurations used to seed the Kriging-based surrogate model of F. Surfaces are coloured by [white when and black when ]. The richness of the parameter space is apparent.
A sample of the 91 valid AVF configurations used to seed the Kriging-based surrogate model of F. Surfaces are coloured by [white when and black when ]. The richness of the parameter space is apparent.
2. 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 xp1, yp1, zp1, , ϕ, ξ, and F for the 29 feasible configurations at which F was evaluated as the MADS method proceeded through so-called SEARCH and POLL steps58,59 are shown in Fig. 4. The Matlab Dace toolbox72 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
Values of xp1, yp1, zp1, , ϕ, ξ, 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.
Values of xp1, yp1, zp1, , ϕ, ξ, 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.
E. Results and discussion
The lowest value of F obtained from the optimisation process was 0.0062. The associated parameter values are reported in Table I. The resulting AVF configuration is shown in Fig. 5, and an STL file is provided in the supplementary material. We note that the optimal configuration features an anastomosis on the outer curvature of the arterial bend, in line 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° 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.
Optimal parameter values.
Parameter . | Optimal value . |
---|---|
xp1 | m |
yp1 | |
zp1 | |
ϕ | |
ξ |
Parameter . | Optimal value . |
---|---|
xp1 | m |
yp1 | |
zp1 | |
ϕ | |
ξ |
Optimal configuration coloured by [white when and black when ]. 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 in the supplementary material.
Optimal configuration coloured by [white when and black when ]. 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 in the supplementary material.
Having identified an optimal configuration, it is useful to assess the sensitivity of F to each of the six geometric parameters. Figure 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 that it is less critical. Further insight can be gained from Fig. 6, which plots the median, inter-quartile range, and range of parameters associated with configurations that achieve , along with the optimal set of parameter values. Once again, it is clear that and ϕ play an important role and must remain close to their optimal values. However, ξ varies significantly, indicating that it is less critical. The plots also suggest that zp1 should remain close to its optimal (non-zero) value, signalling the importance of non-planarity. Finally, Fig. 7 plots variation of Fe, the Kriging-based surrogate model of F after completion of the MADS process, with xp1, yp1, zp1, , ϕ, and ξ, when other parameters are fixed to have values associated with the four most stable configurations specifically, the optimal configuration with F = 0.0062 and three other configurations with F = 0.0078, F = 0.0108, and F = 0.0113. Once again, 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.
Box plots of the median (horizontal black line), inter-quartile range (grey box), and range (whiskers) of parameters associated with configurations that achieve . 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 that it is less critical. The plots also suggest that zp1 should remain close to its optimal (non-zero) value, signalling the importance of non-planarity.
Box plots of the median (horizontal black line), inter-quartile range (grey box), and range (whiskers) of parameters associated with configurations that achieve . 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 that it is less critical. The plots also suggest that zp1 should remain close to its optimal (non-zero) value, signalling the importance of non-planarity.
Variation of Fe with xp1, yp1, zp1, , ϕ, and ξ, when other parameters are fixed to have values associated with the four most stable configurations. Vertical dashed lines denote the optimal parameter values.
Variation of Fe with xp1, yp1, zp1, , ϕ, and ξ, when other parameters are fixed to have values associated with the four most stable configurations. Vertical dashed lines denote the optimal parameter values.
III. PULSATILE SIMULATIONS
A. 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.
B. Governing equations and boundary conditions
Newtonian simulations were undertaken by solving Eqs. (9) and (11), with and . In addition, analogous non-Newtonian simulations were also undertaken using the Carreau-Yasuda model73,74 and the Casson model,75 parameterised as per Boyd et al.76
In all cases, a space-time varying velocity boundary condition was imposed at the PAI. The inflow waveform, which had a period of 1 s, was obtained from Ref. 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 Ref. 49. The resulting time-averaged Reynolds number at the PAI was , 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 QPAI as a function of time is shown in Fig. 8.
Flow waveform imposed at the PAI, QPAI, and the outflow waveforms at the DAO, QDAO, and DVO, QDVO, for the Newtonian simulation.
Flow waveform imposed at the PAI, QPAI, and the outflow waveforms at the DAO, QDAO, and DVO, QDVO, for the Newtonian simulation.
In all cases, a Resistor-Capacitor-Resistor (RCR) Windkessel model was applied at the DAO and DVO. Specifically
and
where PDAO and PDVO are the spatially averaged pressures at the DAO and DVO, respectively, QDAO and QDVO are the outflow rates at the DAO and DVO, respectively, and R1DAO, R2DAO, CDAO, R1DVO, R2DVO and CDVO are relevant Windkessel parameters. These were selected using an approach similar to Pant et al.78 such that PDAO had a physiological range 80–130 mm Hg, and flow splits were similar to those of Sigovan et al.,49 with an average flow split between the DVO and DAO of 66:34, a minimum flow split between the DVO and DAO of 37:63 at the beginning of the acceleration phase, a maximum flow split between the DVO and DAO of 99:1 at the end of the deceleration phase, and a flow split between the DVO and DAO at peak inflow of 45:55. Specifically, the RCR values were obtained using an iterative approach that aimed to minimise
where PPAI is the pressure waveform at the PAI, PRS = 130 mm Hg, and PRD = 80 mm Hg are reference systolic and diastolic pressures, respectively, at the PAI, QR 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 low-resolution time-dependent 3D flow model. The 0D model represented the low-resolution 3D model in terms of three inductors as well as three quadratically non-linear resistors. The 3D model solved Eqs. (9) and (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 Table II.
Windkessel parameters used for the unsteady pulsatile simulations.
Parameter . | Value . |
---|---|
R1DAO (mm Hg ml−1 s) | 0.78 |
R2DAO (mm Hg ml−1 s) | 27.0 |
CDAO [ml(mm Hg)−1] | 0.018 |
R1DVO (mm Hg ml−1 s) | 7.33 |
R2DVO (mm Hg ml−1 s) | 8.46 |
CDVO [ml(mm Hg)−1] | 3.48 |
Parameter . | Value . |
---|---|
R1DAO (mm Hg ml−1 s) | 0.78 |
R2DAO (mm Hg ml−1 s) | 27.0 |
CDAO [ml(mm Hg)−1] | 0.018 |
R1DVO (mm Hg ml−1 s) | 7.33 |
R2DVO (mm Hg ml−1 s) | 8.46 |
CDVO [ml(mm Hg)−1] | 3.48 |
The resulting outflow waveforms QDAO and QDVO at the DAO and DVO, respectively, for the Newtonian simulation are shown in Fig. 8.
Finally, in all cases, a no-slip zero-velocity boundary condition was applied at the AVF wall, which was assumed to be rigid.
1. 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 80 mm Hg 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 solver was run from 0 to 1 s for one pulse period in order to allow transient “start-up” phenomena to pass. The Newtonian simulation [solving Eq. (9) with Eq. (11)], the Carreau-Yasuda model simulation, and the Casson model simulation were then all initiated from this point and run from 1 to 2 s for another pulse period (hence forth referred to as T1), during which time data were exported for analysis, and then from 2 to 3 s for a final pulse period (hence forth referred to as T2), during which time data were exported to verify period convergence. The period convergence results are presented in Appendix B. In all cases, the coupled implicit unsteady solver used a fixed time step of .
A polyhedral unstructured volume mesh, with a prismatic boundary layer mesh adjacent to the wall, was used for each simulation. The mesh was refined near the anastomosis. Specifically, elements near the anastomosis had an average size of , expanding progressively to beyond a distance of from the anastomosis. The prismatic boundary layer meshes were 25 elements thick, with the first element having a thickness of . The mesh had elements in total. This mesh resolution is similar to that used in Ref. 42.
C. Results and discussion
1. Velocity and vorticity
Figure 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 T1 for the Newtonian simulation. Figure 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 T1 for the Newtonian simulation. Figure 11 shows streamlines, originating from the PAI, coloured by velocity magnitude, at four points in T1 for the Newtonian simulation. It is clear that the majority of arterial flow is shifted to the outside of the arterial curvature throughout the pulse cycle, as per Iori et al. It is also apparent that the majority of venous flow is shifted to the outside of the venous curvature throughout the pulse cycle. Finally, it can be observed that a region of re-circulating flow is present adjacent to the inside of the venous curvature, proximal to the anastomosis. The nature of this re-circulating region changes through the pulse cycle. However, fine-scale separated vortical structures, associated with high-frequency flow unsteadiness in previous studies,42–44,51,55 are absent.
Snapshots of u(x) on the arterial symmetry plane, and on a three-dimensional surface that tracks the venous centreline, at four points in T1 for the Newtonian simulation. It is clear that the majority of arterial flow is shifted to the outside of the arterial curvature and that the majority of venous flow is shifted to the outside of the venous curvature throughout the pulse cycle. A region of re-circulation can also be observed to the inside of the venous curvature, proximal to the anastomosis. The nature of this region of re-circulation changes through the pulse cycle; however, disturbed flow associated with high-frequency flow unsteadiness is absent.
Snapshots of u(x) on the arterial symmetry plane, and on a three-dimensional surface that tracks the venous centreline, at four points in T1 for the Newtonian simulation. It is clear that the majority of arterial flow is shifted to the outside of the arterial curvature and that the majority of venous flow is shifted to the outside of the venous curvature throughout the pulse cycle. A region of re-circulation can also be observed to the inside of the venous curvature, proximal to the anastomosis. The nature of this region of re-circulation changes through the pulse cycle; however, disturbed flow associated with high-frequency flow unsteadiness is absent.
Snapshots of ω(x) on the arterial symmetry plane and on a three-dimensional surface that tracks the venous centreline, at four points in T1 for the Newtonian simulation. It is clear that the region of re-circulation observed to the inside of the venous curvature, proximal to the anastomosis varies through the pulse cycle, together with the arterial distal flow. However, fine-scale separated vortical structures, associated with high-frequency flow unsteadiness, are absent.
Snapshots of ω(x) on the arterial symmetry plane and on a three-dimensional surface that tracks the venous centreline, at four points in T1 for the Newtonian simulation. It is clear that the region of re-circulation observed to the inside of the venous curvature, proximal to the anastomosis varies through the pulse cycle, together with the arterial distal flow. However, fine-scale separated vortical structures, associated with high-frequency flow unsteadiness, are absent.
Streamlines, originating from the PAI, coloured by u(x), at four points in T1 for the Newtonian simulation.
Streamlines, originating from the PAI, coloured by u(x), at four points in T1 for the Newtonian simulation.
2. Snapshot proper orthogonal decomposition
A quantitative assessment of unsteadiness was made via snapshot Proper Orthogonal Decomposition (POD)79,80,77 of WSS fluctuations for the Newtonian simulation, where is the WSS field, and is the time-averaged WSS field. Specifically, snapshot POD of was undertaken in T1 using 1000 temporal snapshots with a uniform spacing of 0.001 s. Figure 12 shows the Power Spectral Density (PSD) of a1(t) and a5(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 10 Hz, and in comparison with Fig. 13, which shows the PSD of QPAI, it is reasonable to attribute the majority of this low-frequency content to the pulse waveform. There is certainly no significant energy content around 60 Hz as per Iori et al.42 and Grechy et al.,55 nor above 100 Hz 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.
PSD of a1(t) on a semi-log scale (a) and a log-log scale (b), and PSD of a5 on a semi-log scale (c) and a log-log scale (d) for the Newtonian simulation.
PSD of a1(t) on a semi-log scale (a) and a log-log scale (b), and PSD of a5 on a semi-log scale (c) and a log-log scale (d) for the Newtonian simulation.
PSD of QPAI on a semi-log scale (a) and a log-log scale (b) for the Newtonian simulation.
PSD of QPAI on a semi-log scale (a) and a log-log scale (b) for the Newtonian simulation.
3. Pressure drops
The time-averaged pressure drop between the DVO and PAI for the Newtonian simulation was 2.78 mm Hg, while the time-averaged pressure drop between DAO and PAI for the Newtonian simulation was 1.013 mm Hg. These values are substantially lower than those obtained in previous studies that observed unsteady flow in the vein.51,43,44 This is in line with the expectation that more stable flow is associated with a reduced pressure drop.
4. Non-Newtonian simulation results
Results from the Carreau-Yasuda model simulation and the Casson model simulation were very similar to those obtained via the Newtonian simulations. Quantitative comparisons are provided in Appendix A.
IV. PROOF-OF-CONCEPT
A. Overview
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.
B. 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 AVF) after forming the anastomosis between the artery and the vein.
CAD model of an extravascular device for holding an AVF in the optimal configuration. The device consists of a solid central section, with attached “guttering” sections to accommodate both the arterial and venous segments of an AVF.
CAD model of an extravascular device for holding an AVF in the optimal configuration. The device consists of a solid central section, with attached “guttering” sections to accommodate both the arterial and venous segments of an AVF.
Going forward, 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 mm internal guttering diameter, along with a counterpart device that was mirror-symmetric in the plane of the arterial curve (see Fig. 15).
Devices fabricated in Objet TangoBlackPlus FLX980 using a Connex 350 3D printer. Each device had a mm internal guttering diameter. The left-handed device (a) was mirror-symmetric with the right-handed device (b) in the plane of the arterial curve.
Devices fabricated in Objet TangoBlackPlus FLX980 using a Connex 350 3D printer. Each device had a mm internal guttering diameter. The left-handed device (a) was mirror-symmetric with the right-handed device (b) in the plane of the arterial curve.
C. Porcine model
1. Surgery
A 60 kg Large White/Landrace swine was kept in standard animal care facilities at the Northwick Park Institute for Medical Research (Harrow, UK). A non-recovery experiment was performed with the animal intubated and anaesthetised using a combination of xylazine, ketamine, and isoflurane. AVF were formed in the neck between the jugular and carotid vessels. This is a widely used site for creation of vascular models in the pig since it offers significantly increased access compared to the fore- or hind-limbs.
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 with respect to the animal). Bilaterally, the jugular vein was then cut distally 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 end-to-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 ( 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.
2. Stethoscope recording
Stethoscope spectra with high-frequency peaks have been successfully correlated with the presence of stenosis in vascular connections.81–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 (left-side and right-side AVF) device implantation. Specifically, recordings were made using an electronic 3M Littmann 3200 series stethoscope for approximately 20 s with a sample rate of 4000 Hz, and three independent but otherwise identical recordings were made at each location.
D. Results and discussion
Figures 16–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, respectively. 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.
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.
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.
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.
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.
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.
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.
Figures 19–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 all three recordings for the left-side AVF, before implanting the device, exhibit a peak at Hz. This is in a similar range to peaks found by previous computational42 and experimental84 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 these peaks are suppressed after the device is implanted. Moreover, from Fig. 21 it is clear that the PSDs for the right-side AVF, after implanting the device, look visually similar to their left-side counterparts; without a significant peak at Hz. However, PSDs 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.
PSD produced via a Burg autoregressive model of three independent but otherwise identical stethoscope recordings from the left-side AVF before implanting the device.
PSD produced via a Burg autoregressive model of three independent but otherwise identical stethoscope recordings from the left-side AVF before implanting the device.
PSD produced via a Burg autoregressive model of three independent but otherwise identical stethoscope recordings from the left-side AVF after implanting the device.
PSD produced via a Burg autoregressive model of three independent but otherwise identical stethoscope recordings from the left-side AVF after implanting the device.
PSD produced via a Burg autoregressive model of three independent but otherwise identical stethoscope recordings from the right-side AVF after implanting the device.
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 remodeling, and thus unintentionally cause AVF pathology? And finally (iv) can such a device actually reduce the prevalence of IH and hence AVF failure rates?
V. CONCLUSIONS
High-frequency flow unsteadiness has been associated with development of vascular pathology in AVF.18–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.
SUPPLEMENTARY MATERIAL
See supplementary material for an STL file of the optimal AVF configuration.
ACKNOWLEDGMENTS
The authors are grateful for support from the National Institute for Health Research Imperial Biomedical Research Centre, CD-adapco, the British Heart Foundation, the Aldama Foundation, and the Engineering and Physical Sciences Research Council.
L.G., P.V., and R.C. have a pending patent application for “A device for maintaining vascular connections” (PCT Application No. PCT/EP2017/054420 and UK Application No. 1603501.6). P.V. is also a partner in Quadrature Solutions LLP.
APPENDIX A: NON-NEWTONIAN RESULTS
Snapshot POD of for both of the non-Newtonian simulations was undertaken in T1 using 1000 temporal snapshots with a uniform spacing of 0.001 s. Figure 22 shows the PSD of a1(t) and a5(t), the first and fifth temporal POD modes, for the Carreau-Yasuda model simulation and the Casson model simulation, as well as analogous results from the Newtonian simulation for reference.
PSD of a1(t) on a semi-log scale (a) and a log-log scale (b), and PSD of a5 on a semi-log scale (c) and a log-log scale (d) for the Newtonian simulation (solid line), the Carreau-Yasuda model simulation (dashed-dotted line), and the Casson model simulation (dashed line).
PSD of a1(t) on a semi-log scale (a) and a log-log scale (b), and PSD of a5 on a semi-log scale (c) and a log-log scale (d) for the Newtonian simulation (solid line), the Carreau-Yasuda model simulation (dashed-dotted line), and the Casson model simulation (dashed line).
Results obtained from the Newtonian simulation, the Carreau-Yasuda model simulation, and the Casson model simulation are almost identical.
APPENDIX B: PERIOD CONVERGENCE
Snapshot POD of for all simulations was undertaken in T1 and T2. Figures 23–25 show the PSD of a1(t) and a5(t), the first and fifth temporal POD modes, for the Newtonian simulation, Carreau-Yasuda model simulation, and Casson model simulation, respectively. Results obtained in T2 and T3 are visually identical, indicating period convergence.
PSD of a1(t) on a semi-log scale (a) and a log-log scale (b), and PSD of a5 on a semi-log scale (c) and a log-log scale (d) for the Newtonian simulation in T1 (solid line), and T2 (dashed-dotted line).
PSD of a1(t) on a semi-log scale (a) and a log-log scale (b), and PSD of a5 on a semi-log scale (c) and a log-log scale (d) for the Newtonian simulation in T1 (solid line), and T2 (dashed-dotted line).
PSD of a1(t) on a semi-log scale (a) and a log-log scale (b), and PSD of a5 on a semi-log scale (c) and a log-log scale (d) for the Carreau-Yasuda model simulation in T1 (solid line), and T2 (dashed-dotted line).
PSD of a1(t) on a semi-log scale (a) and a log-log scale (b), and PSD of a5 on a semi-log scale (c) and a log-log scale (d) for the Carreau-Yasuda model simulation in T1 (solid line), and T2 (dashed-dotted line).
PSD of a1(t) on a semi-log scale (a) and a log-log scale (b), and PSD of a5 on a semi-log scale (c) and a log-log scale (d) for the Casson model simulation in T1 (solid line), and T2 (dashed-dotted line).
PSD of a1(t) on a semi-log scale (a) and a log-log scale (b), and PSD of a5 on a semi-log scale (c) and a log-log scale (d) for the Casson model simulation in T1 (solid line), and T2 (dashed-dotted line).