Numerical computation of blood flow for a patient-specific hemodialysis shunt model

Hemodialysis procedure is usually advisable for end-stage renal disease patients. This study is aimed at computational investigation of hemodynamical characteristics in three-dimensional arteriovenous shunt for hemodialysis, for which computed tomography scanning and phase-contrast magnetic resonance imaging are used. Several hemodynamical characteristics are presented and discussed depending on the patient-specific morphology and flow conditions including regurgitating flow from the distal artery caused by the construction of the arteriovenous shunt. A simple backflow prevention technique at an outflow boundary is presented, with stabilized finite element approaches for incompressible Navier–Stokes equations.


Introduction
Over the last decades, end-stage renal disease (ESRD) has emerged as a severe clinical difficulty leading to high rates of morbidity [1,2]. The condition of ESRD is associated with impaired kidney function for which hemodialysis procedures are strongly advisable [3]. Such procedures require well-functioning vascular access placed between an artery and a vein, which is known as arteriovenous shunt. Brescia et al. [4] surgically introduced the use of arteriovenous shunt for hemodialysis in 1966. Morphologies of anastomotic sites play a vital role in efficient shunt operation. They are also responsible for altering hemodynamic forces, as discussed in the literature [5][6][7]. Indeed, with these clinical advancements, some underlying issues remain, such as the existence of intimal hyperplasia to the artery, stenosis development, and insufficient flow occurrence after remodeling, which are related to high longevity of vascular access, as described in clinical studies [8][9][10][11]. These phenomena occur because strong arterial blood flow dynamics affect vein flow directly, which alters the vein hemodynamics dramatically. To maintain arteriovenous shunt functions, elucidation of blood flow fields through them and fluid dynamical forces exerted on the vessel walls are strongly required.
Numerous experimental and computational studies have elucidated the flow dynamics within shunt models [12][13][14][15]. However, because of the large variety of actual shunt morphologies among individual patients, physical experiments using phantoms have presented some difficulties. Patient-specific computational simulation is expected to play an important role in assessment and planning for surgical shunt placements by elucidating complex flow dynamics and wall shear stress (WSS) distributions using medical imaging technologies such as magnetic resonance imaging (MRI) and computed tomography (CT). Browne et al. [16] presented a comparative study designed to observe the mean pressure drop across arteriovenous shunts in physical experiments and computational simulations. Other studies [17][18][19] have presented extensive information about arteriovenous shunt related blood flows.
This study was conducted to investigate a physiological flow in an arteriovenous shunt using patient-specific CT scan data for vessel morphology and corresponding phase-contrast MRI (PC-MRI) data for velocity boundary conditions. The patient case targeted in this study has two characteristics from physiological points of view. The first is that regurgitation occurs in a certain time period from the distal artery, which strongly affects flow behavior near the anastomotic site. Second is that the anastomotic site has three-dimensionally complex morphology because of constraints in surgery. These characteristics engender several difficulties for computation, but provide new insights into medical practices. Although several studies [12][13][14][15][16][17][18][19] have presented extensive information about blood flows for arteriovenous shunts, a class of patient cases with complex morphologies and flow conditions has apparently been only insufficiently investigated. Therefore, the novelty of this paper is presentation of detailed flow behaviors, especially for three-dimensionally complex morphology of anastomotic site and a flow condition with direction change caused by physiological regurgitation.
In simulating blood flows in complicated vessel geometries, unphysical backflows frequently occur numerically, as reported from several studies [20,21]. These phenomena arise in several situations such as occurrences of vortical structures, local flow separations, and widening of the cross-sectional area near outlets. Bazilevs et al. [22] and Moghadam et al. [23] presented outflow boundary conditions to prevent backflows and performed a comparative study of various types of outflow boundary conditions. Because such outflow boundary conditions require consideration from mathematical analysis perspectives, Saito et al. [24], and Zhou and Saito [25] have explained well-posedness and energy inequality of the outflow boundary conditions of Signorini's type for incompressible Navier-Stokes equations. They also proved the existence of a unique solution for Galerkin method with outflow boundary conditions. As described in this paper, we present another type of backflow prevention technique by giving an external force near the outflow boundary without affecting the inner region. It is noteworthy that the unphysical backflow from the outflow boundary is attributable to numerical instability, whereas the regurgitating flow from the distal artery represents a physiological phenomenon.
The remainder of the paper is presented as follows. Section 2 presents governing equations, model construction using medical imaging data, stabilized finite element formulation and backflow prevention technique. Section 3 explains the results and presents some discussion. Section 4 concludes this study.

Patient data
This study, which was conducted in accordance with the Declaration of Helsinki, was approved by the institutional review boards (IRBs) at Akashi Medical Center Hospital (#2020-12) and Tohoku University School of Medicine (#2020-1-511). The solid surface of shunt geometry is extracted from CT scan data of a hemodialysis patient in digital imaging and communications in medicine (DICOM) format. Figure 1 portrays the configuration of an arteriovenous shunt: the target of this study. During the surgery conducted for this patient case, artery-side-to-end-vein anastomosis [26] was performed between the side of the radial artery and the end of the cephalic vein in the wrist area. This configuration consists of three parts connected at the anastomotic site presented in Fig. 1: the proximal artery (PA), distal artery (DA), and distal vein (DV). Therefore, three boundaries, PA , DA , and DV , are included in our configuration in addition to the vessel wall boundary wall . Here, the direction vectors of centerlines of three vessels PA, DA, and DV cannot be laid

Governing equations and model configurations
Let denotes the domain boundary, and n sd = 3 represents the number of space dimensions, we assume that blood behaves as a Newtonian fluid, which is governed by unsteady, incompressible Navier-Stokes equations.
where , , and respectively denote the density, velocity vector, and the unit outward normal vector on the boundary , and ( , p) = − p + 2 ( ) is the stress tensor. Here p, , and are the pressure, identity tensor, and blood viscosity respectively. Also, ( ) = 1 2 ∇ + (∇ ) T is the strain-rate tensor. Stress in the matrix form is expressed as ij = −p ij + Here, ij is symmetric, which implies that ij = ji . Boundary conditions (3) and (4) respectively represent the typical form of the essential and natural boundary conditions defined on g and h with the given functions g and h. Divergence free velocity vector 0 is given as an initial condition. In this study, g = wall ∪ PA ∪ DA is a disjoint union of the vessel wall and the inlet or outlet boundaries. Also, DV is set as h defined in (4). The vessel wall is assumed to be rigid. Therefore, a no-slip condition = is applied on wall in (3).

Stabilized finite element formulation
Numerical solutions of (1)-(5) with complex geometries are obtainable using Galerkin finite element methods. However, such formulations are often associated with potential numerical instabilities [27]. Tezduyar et al. [28,29] proposed stabilized finite element formulations for incompressible fluid flows. This formulation also allows for equal-order polynomial approximation for velocity-pressure variables.
For discretizing (1)-(5), we used the streamline upwind Petrov-Galerkin (SUPG) and pressure stabilizing Petrov-Galerkin (PSPG) methods [28,29]. Spatial domain is discretized into the elements e , e = 1, 2, … , n el , where n el represents the where P 1 stands for the set of first-order polynomials in space. The stabilized finite element formulation of (1)- (5) with the stabilization parameter is given below.
The advection velocity is represented by in (10) and (11). These equations are rewritten in the form of The space discretization of (12) and (13) leads to the following system of ordinary differential equations.
Therein, matrices M, C, D, and G respectively correspond to acceleration, convection, diffusion, and pressure gradient terms. Vectors F and E are related to the boundary conditions. Subscripts S and P respectively stand for the SUPG and PSPG contributions. Crank-Nicolson method is used for time-discretization of differential equations (14) and (15). Then the advection velocity is approximated by the Adams-Bashforth method, which leads to the system of linear equations. The stabilization parameter for each element is given as [28] where t and || e || respectively represent the time-step and norm of the advection velocity. The element length h e of e computed as where is the unit vector in the direction of flow, n en denotes the number of nodes in an element, and N a is the basis function corresponds to the node a.

Method for preventing backflows
In this section, we examine a backflow presence at the outflow boundary and a prevention technique for it. Computational blood flow simulation for complex vessel geometries extracted from realistic medical data mostly yield undesirable instabilities at outflow boundaries. This situation might occur because of a lack of accurate velocity profiles or necessity for setting relevant boundary conditions for these boundaries. We consider a stabilized weak formulation of the governing equation given in Sect. 2.3 with an artificial external force term as shown below.
where denotes the artificial external force term to account the backflow instability at the outflow boundary h . For this study, we express this term as a function of the velocity vector, Here is a positive constant, ̃h h represents a region of backflow near h , and ̃h ( h ) denotes the indicator function, We choose a union of elements in which the backflow occurs at some boundary node to be ̃h h ;

3
Numerical computation of blood flow for arteriovenous shunt Consequently, ( h ;̃� h ) is linear in the velocity vector h if the region of backflow, ̃′ h , is fixed. The last term of (18) is expected to be applicable in the situation of flow reversal. After implementation of this external force, the system of equations becomes the following.
for h ∈ X h n sd and p h ∈ (X h ) n sd . We discretize acceleration and artificial external force terms M + M S Zhou and Saito [25] proved that boundary conditions of Signorini's type satisfy the energy inequality, by which backflows are prevented by introducing a conditional vector, similarly to our technique. The energy inequality is not assured, in general, because of a possibly negative boundary integral term in the weak formulation when we use the do-nothing boundary condition. Although our prevention technique is combined with the do-nothing boundary condition, we might expect that the possibly negative term corresponding to the backflow is canceled by adding artificial force .

Boundary conditions using PC-MRI
For this study, boundary conditions are constructed using PC-MRI measurements. As presented in Fig. 1, our target configuration has three inlet or outlet boundaries: PA , DA , and DV . Flow rates Q PA = ∫ PA ⋅ ds and Q DV = ∫ DV ⋅ ds at 36 discrete time-steps during one heart period are obtained by PC-MRI. Cubic spline interpolation is used to obtain flow rate profiles with sufficient time resolutions. Because we assume incompressible fluid flow, total incoming and outgoing flow (22) rates must be conserved, i.e., Q PA + Q DA + Q DV = 0 , where flow directions are taken in outward normal directions. Therefore, we calculate Q DA = − Q PA + Q DV . We set DV as h with = 0, and DA as g because this boundary plays an important role in reproducing the flow dynamics for this patient's case with regurgitation. Actually, the velocity direction at this boundary changes from outlet to inlet and vice-versa during one heartbeat. Another choice for boundary conditions is to set DV as g and DA as h . However, this setting is too restrictive because we only have integrated flow rate data on DV boundary but do not know the velocity profile there, and our target is the complex flow near the anastomosis site, especially in the vein part.
According to the special feature of this patient case with flow direction change at DA because of regurgitation from DA, we define splitting and merging phases. In the splitting phase, incoming flow from PA splits into DA and DV from t = 0 to t = 0.18 and t = 0.49 to t = 1.1, as presented in Fig. 2. In the merging phase, DA acts as an inflow boundary, which means that two incoming flows from PA and DA merge at the anastomotic site. The merging phase lasts from t = 0.18 to t = 0.49 as presented in Fig. 2. The flow rate at DA becomes 0 at t = 0.18 and t = 0.49, at which flow behaviors at DA change respectively from outlet to inlet and vice-versa. Detailed characteristics of flow behaviors in these phases are discussed in Sect. 3.

Mesh generation and numerical computation
Vascular modeling toolkit (VMTK) [30] was used for drawing the centerline in the extracted solid model. Then, a tetrahedral finite element mesh was generated using Gmsh [31] with a thin layered mesh near the vessel wall to resolve a large velocity gradient in boundary layers. The finite element mesh is depicted in Fig. 3, where the total numbers of nodes n nd and elements n el are, respectively, 26292 and 140052.
A sparse non-symmetric system of linear equations is solved using GPBi-CG method, as discussed in [32,33]. Density and blood viscosity are set respectively as 1060.0 kg/m 3 and 2.66 × 10 −3 N s/m 2 . Time step length of t = 2 × 10 −4 s is adopted for a cardiac cycle T = 1.1 s. We have heuristically taken = 0.5 in this study for efficient prevention of backflow from outflow boundary DV .

Figures 4 and 5 portray velocity vectors on the outflow boundary DV with outward normal vectors
DV in black color. In Fig. 4a at t=0.09, backflow appears in  Fig. 4b at t=0.14 and immediately diverges. Figure 4c shows the magnitude of velocity on the plane perpendicular to DV . When applying the external force presented in Sect. 2.4, the instability is well suppressed. Figure 5 portrays velocity vectors and their magnitude for the same time steps as those shown in Fig. 4. This comparison portrays that our backflow prevention technique works well in the region near the outflow boundary without exerting too strong an effect on upstream regions.

Periodicity check
We have no appropriate initial conditions. Therefore, we performed a periodicity check by computing flow fields for five cardiac cycles starting from an initial condition 0 = . Spatially-averaged WSS magnitudes are depicted in Fig. 6 for different cardiac cycles. Based on the periodicity check, we adopted the computational results from third to fourth periods (from t = 2.8 to t = 3.9 ) as portrayed in Fig. 7 for discussing blood flow characteristics in the following subsections. The four yellow lines of (a), (b), (c), and (d) in Fig. 7 respectively correspond to time-instants t = 2.88, 3.48, 3.63, and 3.78. The period of (b)-(d) is a merging phase, whereas the other is a splitting phase. Lines (a) and (c) respectively correspond to the instants at which maximum flow comes from PA and DA in the splitting and merging phases. Lines (b) and (d) respectively correspond to time instants at which phase changes from splitting to merging, and from merging to splitting.

Characteristics of flow fields in the arteriovenous shunt model
As described in Sect. 2.5, blood enters from PA and exits from both DA and DV during the splitting phase, whereas blood enters from both PA and DA, and exits from DV during the merging phase. Therefore, DA plays two roles in the present  Figure 8a shows the splitting phase in which the flow from PA splits into DA and DV in a complex way affected by the threedimensional morphology of anastomotic site. Then, flow rate in DA begins to decrease and change directions, which causes strong circulation in DA, as presented in Fig. 8b. In the merging phase presented in Fig 8c, flows coming in from both PA and DA are comparable. Then, because of the morphology around the anastomotic site, it merges to form strong spirals in DV. Figure 9 presents the same timings with Fig. 8 from a different viewing angle. In Fig. 9c, because of the patient-specific morphology of this patient case, formation of spiral flow in DV driven by two incoming flows from PA and DA is readily apparent, which implies that the flow in merging phase is much smoother than that in the splitting phase. It is noteworthy that if an arteriovenous shunt is not surgically constructed, then no direct connection exists between the radial artery and cephalic vein. Blood coming in from PA simply passes out of DA. Therefore, the patient case analyzed in this paper represents a specific problem in blood flows with arteriovenous shunt.

Characteristics of this morphology from physiological viewpoints
Wall shear stress (WSS) exerted on the endothelial surface of blood vessels is defined as where is a traction vector computed from the stress tensor and surface normal vector, as = . We first compute the elementwise WSS distribution. Then, we get the nodal approximated distribution as the local arithmetic mean of elementwise one. Figures 10 and 11 respectively present the WSS distribution for splitting and merging phases from two viewing angles. At t=2.88 during the splitting phase in Fig. 10, maximum flow coming in from PA, colliding to the anastomotic site wall, and splitting into DA and DV causes high WSS around the anastomotic site. At t=3.63 during the merging phase in Fig. 11, two incoming flows from PA and DA merge smoothly and generate swirling flow in DV, causing the slightly high WSS near the anastomotic site. Other high-WSS regions are apparently because of the local bumps in this patient's morphology.
High OSI denotes areas with frequent changes in WSS vectors directions. Figure 12 portrays the streamlines in splitting and merging phases and the OSI distribution. Around the root of DV near the anastomotic site, circulations rotating in opposite directions in splitting and merging phases are apparent, which brings about high OSI there. OSI is known to be related strongly to atherosclerosis [34]. Therefore, this information can be important for planning surgical procedures considering shunt morphologies. Fig. 12 Circulations at the anastomotic site and OSI distribution

Concluding remarks
This study was conducted to investigate details of blood flow dynamics in an arteriovenous shunt model extracted from patient-specific CT scan and PC-MRI data. We have presented a simple backflow prevention technique to avoid computational instability. Characteristics of hemodynamics within the arteriovenous shunt model depend greatly on patient-specific morphologies and flow conditions, which differ widely among individual patients. In the patient case examined for this study, complex three-dimensional morphology around the anastomotic site and the existence of regurgitating flow from the distal artery (DA) bring about characteristic flow structures and WSS/OSI distributions in different phases and in transitions between them, which provides useful information for clinical practice.