Blood Flow Simulations of the Renal Arteries - Effect of Segmentation and Stenosis Removal

Patient specific based simulation of blood flows in arteries has been proposed as a future approach for better diagnostics and treatment of arterial diseases. The outcome of theoretical simulations strongly depends on the accuracy in describing the problem (the geometry, material properties of the artery and of the blood, flow conditions and the boundary conditions). In this study, the uncertainties associated with the approach for a priori assessment of reconstructive surgery of stenoted arteries are investigated. It is shown that strong curvature in the reconstructed artery leads to large spatial- and temporal-peaks in the wall shear-stress. Such peaks can be removed by appropriate reconstruction that also handles the post-stenotic dilatation of the artery. Moreover, it is shown that the effects of the segmentation approach can be equally important as the effects of using advanced rheological models. This fact has not been recognized in the literature up to this point, making patient specific simulations potentially less reliable.


Introduction
Advances in medical imaging technologies have led to the possibility to study patient specific blood flow in more detail by utilizing computational fluid dynamics (CFD) models in geometries derived from medical images. The advantages of using patient specific data in order to assess the details of blood flow in the vessel of interest also enables assessment of the merits of different possible surgical measures. Computed Tomography Angiography (CTA) with improved spatial resolution has become an important tool in diagnosis of cardiovascular diseases. The diagnostic quality of the non-invasive CTA technique is now almost on par with invasive (catheter) angiography for detecting for example coronary artery diseases (c.f. [1,2]). CTA enables measurements of not only the shape of the lumen of the vessel but also the thickness of its wall. However, the accuracy of such measurements depends on the accuracy of the segmentation method used to delineate the arterial wall. Significant differences in results among different algorithms were demonstrated in a large study comparing different segmentation algorithms, as reported by Kirişli et al. [3]. In the literature, there are several review papers describing the most common approaches for vessel lumen segmentation (c.f. [4][5][6][7]).
The lumen segmentation of a blood vessel may act as a first step in several different applications such as better diagnosis of stenosis, evaluation of possible of surgery alternatives or, as in our case, reconstruction of the vessel to the shape prior to being affected by atherosclerosis. In all these cases, the required level of accuracy of the lumen segmentation may be different; the accuracy required for determining the presence of a stenosis is much lower as compared to the accuracy needed for a reconstruction that would resemble the shape of a healthy vessel.
Atherosclerosis is a disease that develops over a long period of time (decades). The process is slow and manifests in form of morphological changes in the arterial wall at several common locations (branches of certain arteries). A hypothesis that has been around for decades states that the blood flow in these locations is a major factor in the formation of stenotic lesions. It has also been hypothesized that the process is related to the Wall Shear Stress (WSS) [8]. Yet, it has been claimed that both high WSS [9] as well as low WSS [10] can play a role in the processes. In recent years, the temporal-and spatial gradients of the WSS has been recognized as more important as compared to the absolute value of the WSS itself [10]. Several measures for defining the accumulated effect of the WSS has been proposed (e.g. [11][12][13]). The most common criteria used nowadays are related to averages of the WSS. To this group of measures one may include the Time-Averaged WSS (TAWSS), the Oscillatory Shear Index (OSI) (c.f. [14][15][16]), and the Relative Residence Time (RRT) (cf. [17]). In order to assess such prognostic criteria for the development of atherosclerosis, the logical starting point is the affected artery as it looked like at an earlier time, i.e. at its healthy state. To achieve this, CTA data needs to be converted to a three-dimensional structure and thereafter "reconstructed" into the smooth shape that the vessel originally had. The main issues are how to reconstruct the vessel as well as the needed accuracy for the predictive purpose that we seek. The latter question can be expressed in terms of sensitivity of vessel pathology characterizing parameters, such as TAWSS, OSI and RRT to errors in the reconstruction. In recent years, patient specific geometries for numerical simulations has also been used where the effects of the blood rheology models were considered (c.f. [16]). Blood rheology depends, among other factors, primarily on the local hematocrit (Red Blood Cell, RBC, concentration) and the local shear-stress. The effects of rheological model are non-negligible and should be taken into account. In this work, it will be shown that the treatment of segmentation is equally important as the rheological approaches we applied.
In the following, the computational framework for arterial blood flow simulations is described where after methods for reconstructions are addressed. The sensitivity of flow parameters to the reconstruction is detailed in the result section.

Computational Methods
The computational methods used to analyze blood flow in the segmented arteries are presented below. First, the set of model equations are discussed, followed by the methods used to reconstruct the diseased artery into a healthy state.

Governing equations and boundary conditions
The computational model is based on the following assumptions: Blood is an incompressible, non-Newtonian fluid whose properties depend on the local Red Blood Cell (RBC) concentration. Arterial flows are by nature always time-dependent due to the heart pulsations, and in diseased arteries, the flow is transitional or intermittently turbulent. In this work, the conservation of momentum and mass are supplemented by a convection-diffusion (transport) equation for the RBC concentration (α).
where γ = (∇u + (∇u) T )/2 is the shear rate tensor, ρ = (1 − α)ρ plasma + αρ rbc the density of the mixture of plasma and RBCs, D the diffusivity of the RBC phase, and μ the mixture viscosity. Two different viscosity models are considered; a Newtonian blood analog fluid and a non-Newtonian Quemada model. For the Newtonian model, a constant value of μ = 3.38 · 10 −3 Pa·s is used. In the Quemada model [18,19], the viscosity of blood depends on the local RBC concentration and the local shear rate. This viscosity model has the following form: where γ = √ 2(γ : γ ), μ plasma = 1.32 · 10 −3 Pa·s. The viscosity coefficients, k 0 and k ∞ , as well as the critical shear rate, γ c , have been fitted to a large range of hematocrits and shear rates by Cokelet [20]: The geometry of the blood vessels, including the descending aorta and the renal arteries is obtained from CTA data from an anonymous patient. In the following, the focus is only on a segment of the abdominal aorta and its main branches, including the renal arteries, depicted in Fig. 1a. As noted by the black arrow, the left renal artery has a severe stenosis. At the inlet into the abdominal aorta, the volume flow is time-dependent as in the other branches of the aorta. At the inlet into the abdominal aorta, the velocity distribution is provided by previous simulations of the descending aorta [21]. All branching arteries have been extended in order to reduce the influence of the outlet boundary conditions. A constant (spatial) average gauge pressure (p = 0) is imposed at the left and right renal artery outlets, allowing the flow to redistribute between the arteries. For the other outlets, the flow rate is set to a fraction of the instantaneous inflow. The fraction is 1/3 for the descending aorta, and 1/9 each for the superior mesenteric artery and the celiac artery. The flow rate in the different artery branches are depicted in Fig. 1b. The RBC transport equation, Eq. 3, requires a condition on all boundary points. The hematocrit at the inlet is described by a hyperbolic tangent distribution along a radial line, shown in Fig. 2a. This profile implies that there is a lower RBC concentration near the wall of the inlet artery, a feature that has been physiologically verified [22].
where α min = 0.25, α = 0.25, a = 4000/m, d the wall distance and d min = 0.8 mm. The resulting distribution is shown in Fig. 2b. The walls of the arteries are assumed to be solid, hence the no-slip condition is imposed (u = 0). This assumption is justified for atherosclerotic arteries. For the RBC concentration, no flux of cells into the wall is assumed. Hence, the wall-normal gradient of the hematocrit is assumed to vanish (i.e. n wall · ∇α = 0). A zero-gradient condition is applied for the hematocrit at all outlets.

Flow solver
The governing equations are solved with an inhouse version of the OpenFOAM solver, with added support for hematocrit dependent viscosity models. A second order finite volume discretization is applied to the governing equations, Eqs. 1-3. Convective fluxes are evaluated with limited central differences (see [23]) . Diffusive fluxes are evaluated in a similar fashion, but with corrections applied for mesh non-orthogonality. Time derivatives are approximated to second order with an implicit backward time stepping scheme, where t is the timestep and φ is the considered field. The time step is adapted throughout the simulations to keep the Courant number close to, but below 1. The time step size ranges between 10 −4 and 2 · 10 −6 s, where the smaller timestep is applied during peak systole for the stenoted geometry, and the higher during diastole for the non-stenoted case. The discretized equations are solved sequentially, starting with the hematocrit (3), coupled to the flow equations through the density and viscosity fields. The momentum equation is solved component by component, after which three pressure correction steps with the PISO algorithm are performed in order for the velocity field to fulfill (1).
The initial value of the velocity and pressure fields are set to zero, whereas the hematocrit is α = 0.45. The simulations are carried out over a few heart cycles such that the effects of the initial conditions diminish. This is carried out by monitoring the cycle-to-cycle variation. Once the effect of the initial conditions is small, statistics of the relevant parameters is evaluated for one heart-cycle.

Vessel Reconstruction
A patient specific vessel geometry, including the descending aorta, the superior mesenteric artery, the celiac artery, and the renal arteries is considered. The left renal artery has a stenosis where the vessel area abruptly reduces by approximately 90%. Segmentation of CTA data was performed from which the vessel surface was extracted. The degree of smoothing applied to the geometry was varied in order to enable for the sensitivity to geometry perturbations to be assessed.

Stenosis removal
Starting from the patient geometry, the left renal artery without stenosis was reconstructed. Two different methods were applied; a Voronoi diagram based method and a method based on mesh sweeping. The general idea behind the methods can be summarized as: 1. Compute the centerline of the stenoted vessel 2. Cut out the stenosis from the surface and the centerline 3. Recompute the centerline in the removed part of the vessel 4. Interpolate the surface along the reconstructed centerline 5. Stitch together the meshes Each of these steps will be described in more detail below.

Centerline computation
Following the work presented in [24], the centerline of a vessel is defined as the path between the centerpoints of the inlet and the outlet section that maximizes the distance to the boundary. The curves are obtained with the vmtk library [25], according to the algorithm proposed in the aforementioned reference.

Voronoi diagram based reconstruction
The Voronoi diagram based reconstruction is a modification of a method originally developed for reconstruction of vessels with aneurysms [26]. The method operates on the Voronoi diagram of the vessel that can be described as a partition of the lumen volume into a set of minimal spheres. The spheres may overlap, but no sphere is fully contained within another (see Fig. 3a).
Two cutting points, placed upstream and downstream of the stenosis, are chosen, and the intermedial centerline points and Voronoi vertices are removed (Fig. 3b). The now missing part of the centerline is replaced by a cubic spline (Fig. 3c), used in order to guide an interpolation of the Voronoi diagram along the gap. For the interpolation, a set of Voronoi vertices close to each cutting point is chosen and vertex start-end pairs are formed. The associated sphere radius as well as distance to and rotation angle around the centerline is interpolated along the line (Fig. 3d). The interpolated Voronoi cells are mapped to a grid from which the surface is reconstructed through the Marching Cube algorithm [27] (Fig. 3e).

Mesh sweep reconstruction
The mesh sweep method is similar to the Voronoi method, but with the main difference that the geometry manipulations are performed directly on the surface geometry. The original surface and centerline (see Fig. 4a) are divided by two cutplanes whose normals are taken to be the local tangents of the centerline (Fig. 4b). The reconstructed centerline is approximated by two piecewise cubic splines, r 1 (s), r 2 (s), of the form, where s is 0 at the start point of the spline and 1 at the end, and a i , b i , c i and d i are constants. At each point along the spline, an orthogonal coordinate system can be defined (the Frenet Frame), consisting of the curve tangent, t F , the normal, n F , and the binormal, b F . These are defined as where the dot indicates differentiation with respect to s. The splines are required to pass through the centerpoints (p 0 and p 2 ) of the cutplanes used to remove the stenosis, and with tangent vectors coinciding with the plane normals (n 0 and n 2 ). This gives four conditions for the spline constants: An additional, intermediate guiding point is taken from the original centerline (p 1 ) and the last four conditions can be derived by requiring the position, as well as the Frenet tangent and normal to be continuous at the knot point: To prevent the curve from forming loops and folds, the scaling laws for the spline constants derived by Yong et al. [28] are applied. The resulting non-linear system of equations are solved iteratively with the Newton-Raphson method.
Once the corrected centerline has been computed (Fig. 4c), the non-stenoted geometry can be obtained by sweeping the mesh along the line. The only known information is the surface curve at the start and the end point. A transformation is thus required that smoothly morphs the surface between these curves. The surface curves are transformed into their local Frenet frame, and an optimal pairing of points is obtained by finding start-end point pairs minimizing the Frenet distance. The Frenet coordinates are interpolated linearly along the centerline, shown in Fig. 4d. A transformation to world coordinates gives the points on the surface (Fig. 4e). Triangular elements are created connecting the surface points, yielding the final reconstructed surface (Fig. 4f).

Vessel area evaluation
Using the centerline as starting point, the area distribution of a vessel is computed. At each point along this line, the position, p i , and the tangent vector, t i are known. The local area is defined as the area of the intersection between the vessel volume, V vessel and a cutplane having an origin at p i and normal t i . If this plane is denoted P i , the local area can be written as The volume-plane intersection is an expensive operation. Thus, it would be preferable to compute the area from the intersection of the vessel surface and the plane instead. From Stokes theorem it is known that, where n is the normal of the surface (here, the centerline tangent, t i ). Set F = 1 2 n × r. If the normal is constant, it is straightforward to show that ∇ ×F · n = 1, and Eq. 15 can be restated as, The procedure applied to compute the vessel area distribution can be summarized as: For each centerline point i: 1. Find the intersection between the vessel surface and a plane with origin p i and normal t i , producing an ordered list of points {r j } j making up the perimiter of the cutsurface.

Results
The effects of the methods used for the segmentation of the entire geometry and reconstruction of the left renal artery are considered by studying different combination of segmentation approaches and rheological models. Such a comparison provides insight into the relative importance and uncertainty associated with the treatment of the CTA data. Similarly, by comparing the effect of the segmentation/reconstruction and rheology on common parameters for characterizing risk for atherosclerosis, the possibility of using these parameters as risk for disease indicators can be assessed. The geometrical effect of handling the left renal artery can be shown by considering the cross-sectional area along the axis of the artery. Figure 5 depicts the cross-sectional area for the different cases. The original segmentation case ("Stenoted") curve shows the presence of a severe stenosis ranging from approximately 0.5 to 1.5 cm downstream of the renal artery bifurcation. Further downstream and up to about 3 cm, an increase in the area of the artery can be noted, commonly referred to as "post-stenotic dilatation". A minor area increase can also be observed, located just prior the stenosis due to flow separation at that location. Using the Voronoi approach for stenosis removal is limited to only the stenosis itself since the reconstructed centerline otherwise deviates significantly from the stenosed case. Hence, this approach has no effect on the post-stenotic dilated cross-sectional area. Using mesh sweeps, both the stenosis and the post-stenotic expansion are removed, leaving a smooth and monotonically decreasing cross-sectional area. The non-uniformities of the geometry may lead to unsteady flow separation, increased losses as well as larger and fluctuating WSS. The reconstruction methods outlined above require input of coordinates for the start and endpoints of the stenoted region. The effect of moving the points on the area distribution of the reconstructed vessel is shown in Fig. 6a and b. It can be concluded that both methods are sensitive to the choice of input coordinates. Both methods have been found to be fairly insensitive to mesh resolution, as long as a centerline can be computed for the stenosed case. This, in turn, requires the resolution to be finer than the minimal stenosis diameter.
As a pressure condition is applied at the outlet of the renal arteries, the flow may be redistributed (though within the limits imposed by the mass flow through the other vessels). The flow rate in the two arteries is given in Fig. 7a and b. As shown in these figures, the flow rate is identical for both reconstruction methods. Moreover, the most important parameter is the presence or absence of stenosis.
The most common way of considering losses due to a stenosis is by considering the pressure drop over an arterial segment. This parameter can be measured directly by catheterization or estimated from flow measurements. The pressure loss over the left renal artery as a function of time for different artery shape treatments is depicted in Fig. 8a. Another measure is the needed power to drive the flow. This parameter (the integral of the pressure and local normal velocity to the inlet and outlet boundaries) expresses more correctly the energy loss per time unit (J/s). However, this parameter is not commonly used in clinical situations since it requires the spatial and temporal distribution of both the pressure and velocity in two cross-sectional planes of the artery. The losses across the left renal artery, corresponding to the pressure drop, are depicted in Fig. 8b. As seen, the presence of stenosis results in significant differences in both pressure drop and losses as compared to the reconstructed cases. However, the loss of energy is a more sensitive parameter and a difference between  Fig. 7 The flow-rate vs time in the right a and the left renal arteries b before and after reconstruction using different segmentation procedures. Note that the curves for the mesh sweep and Voronoi stenosis removal methods coincide the Voronoi and the so called Sweep approaches is detected, although this difference is relatively small.

Atherosclerotic indicators
The wall shear stress has been considered to be a central factor when it comes to atherosclerosis. The dependence of this parameter on the segmentation method is essential for the assessment of the models intended to be used to predict risks for an arterial wall disease. Figure 9a and b depict the spatial mean and peak of the wall shear stress as function of time, respectively. As shown, the WSS is largest for the case with stenosis with the peak coinciding with the peak flow rate. The smoothed cases have smaller WSS values, in particular with respect to peak WSS.
The effect of segmentation on the commonly used Time Averaged WSS (TAWSS) and the Oscillatory Shear Index (OSI) is considered in Fig. 10. Here, the spatial distribution of the cross-sectional average of these time-averaged parameters is shown. The TAWSS has, as expected, a very strong peak at the stenosis and at the region where the artery cross-section experiences non-uniformity, just upstream of the reconstructed region. The OSI behaves similarly. This index is a dimensionless number with maximal value of 0.5. The OSI is strong not only at the stenosis, like the TAWSS, but also at locations where the geometry is non-smooth. Both the smoothed and the Voronoi approaches do have a large value at the edge of the reconstructed, non-smooth region. For the reconstructed arteries, a lower level of TAWSS is observed, which varies in space due to the spatial non-uniformity of the artery (as expressed in the cross-sectional area of the artery). The variations in OSI are much stronger at the stenosis and much weaker at the post-stenotic dilatation region. However, OSI is less sensitive to the type of reconstruction even though there are significant differences in terms of geometrical parameters, losses and WSS.
The Relative Residence Time (RRT) is defined commonly as function of the OSI and TAWSS [17]. The RRT using different segmentations is depicted in Fig. 11a. The Figure  also compares the effects of the rheological models. In Fig. 11b, the Voronoi method is The importance of WSS gradients as an indicator for potential arterial wall pathology has been suggested [10]. Obviously, as the WSS is not smooth in a vessel with surface non-smoothness, the spatial gradient of the WSS becomes a sensitive indicator for such non-uniformities. Figure 12a depict the sensitivity of the WSS gradient to identify the stenosis. The level of the WSS gradient is reduced by more than an order of magnitude for the reconstructed vessels, using sweeps or Voronoi's method. Moreover, it can be noted that the fluctuation in WSS gradients diminish for both methods downstream of the reconstructed region. For the Voronoi approach, an additional peak is observed around the end of the post-stenotic dilation (position D). The mesh sweep method yields a smoother gradient, and with a small peak around position (C), corresponding to the point of maximal centerline curvature.

Summary
In order to be able to carry out reliable patient specific blood flow simulations of artery reconstruction, one has to be aware of the effects due to the segmentation method used for the original CTA images. The level of uncertainty due to segmentation techniques is of the same order as uncertainties and errors due to using rheological models neglecting temporal-and spatial-variations in the concentration of RBC. Commonly used parameters for characterizing atherosclerosis, such as the level of WSS, its temporal and/or spatial averages, or parameters as the OSI, RRT and the WSS streamwise gradient, are all sensitive to different degree to the segmentation parameters. Thus, these parameters should be used with caution when different segmentation methods are used.