Sensitivity of hemodynamics in a patient specific cerebral aneurysm to vascular geometry and blood rheology. Mathematical Biosciences and Engineering

Newtonian and generalized Newtonian mathematical models for blood flow are compared in two different reconstructions of an anatomically realistic geometry of a saccular aneurysm, obtained from rotational CTA and differing to within image resolution. The sensitivity of the flow field is sought with respect to geometry reconstruction procedure and mathematical model choice in numerical simulations. Taking as example a patient specific intracranial aneurysm located on an outer bend under steady state simulations, it is found that the sensitivity to geometry variability is greater, but comparable, to the one of the rheological model. These sensitivities are not quantifiable a priori. The flow field exhibits a wide range of shear stresses and slow recirculation regions that emphasize the need for careful choice of constitutive models for the blood. On the other hand, the complex geometrical shape of the vessels is found to be sensitive to small scale perturbations within medical imaging resolution. The sensitivity to mathematical modeling and geometry definition are important when performing numerical simulations from in vivo data, and should be taken into account when discussing patient specific studies since differences in wall shear stress range from 3% to 18%.


1.
Introduction. An increasing demand from the medical community for scientifically rigorous and quantitative investigations of vascular diseases has recently given a major impulse to the development of mathematical models and numerical tools for the computer simulations of the human cardiovascular and cerebrovascular systems, in both healthy and pathological states. The cardiovascular system is not degrading steadily over time, but rather is capable of growth and remodeling in response to changes in mechanical and biochemical stimuli [9,13]. The initiation and progression of a number of threatening diseases, like intracranial aneurysms (IA) are believed to be intimately tied to the inability of the vasculature to respond satisfactorily to such complex changes.
IA occur commonly at the apices of bifurcations and outer bends of curved arterial segments. It has been postulated that IA formation, growth and rupture are associated with local hemodynamics as well as lumen structural mechanics. However the mechanisms and parameters are not fully understood and optimal treatment stratagies are still actively sought. Commonly, parameters associated with diseases in the cerebrovascular system are located at or near the lumen wall, including wall shear stress (WSS) and its spatial gradients (WSSG), the jet and impingement sizes [7,23,17]. Abnormal values of WSS and WSSG can be associated to disturbed flow, which is linked to both the local and large-scale conduit morphology.
The role of blood flow physiological parameters in regulating aneurysm morphology and natural history is poorly understood [10]. The mechanical response of blood largely depends on the state of its red blood cells, which can range from three dimensional structures to relatively uniformly dispersed individual cells. Throughout most of the arterial system of healthy individuals, the red blood cells are dispersed and it is considered sufficient to model blood as an inelastic, constant viscosity fluid (Newtonian) [21]. However, in some disease states the vascular geometry is locally altered and induces relatively stable regions of slow recirculation (e.g. in IA or downstream of a stenosis). In such cases, non-Newtonian behavior may become important and constitutive models with shear-thinning viscosity, viscoelasticity, thixotropy or yield-stress could be necessary (see [1,3,21] and [9,Chapter 6]).
The hemodynamic flow field strongly depends on boundary surface definition; it is therefore necessary to study aneurysms using realistic geometries. Experimental studies using idealized aneurysm models tend to simplify the complex flow patterns observed in aneurysms while still providing a wealth of knowledge to this challenging topic [12,24]. However, there is also a need to perform studies that are patient specific. Medical data acquisition has greatly progressed with recent development of high resolution imaging that is tendencially non-invasive, providing in vivo clinical data that can be readily used to perform both experimental and numerical tests, that yield fine detail and a greater set of measures. Correlation of cerebral vasculature shape to risk of aneurysm rupture or formation have brought to light several key parameters, principally associated with the upstream curvature and radius of the vessel [22] and the size, shape and characteristics of the branching vessels [14,15]. These geometric factors are insensitive to the smaller scales, however that are typical of uncertainty inherent from the medical image resolution.
Despite the feasibility and capability of performing detailed and sophisticated numerical simulations, the reproducibility of patient specific studies from clinical data is compromised by modeling uncertainty. This stems principally from both the errors associated with using in vivo data that cannot be benchmarked, and also from mathematical modeling of both the blood, vessel walls and biochemical effects and responses. Studies to quantify the sensitivity of specific parameter change with respect to hemodynamic measures have indicated so far a dominant influence of the variability in the geometry definition reconstructed from medical images [6] and parent-artery inflow geometry [5] in most cases studied.
The present work is devoted to the sensitivity of parameter changes in the numerical simulations with respect to the computed flow field. The sensitivity arises from modelling uncertainty in deriving numerical techniques and mathematical models that represent the in vivo physiological environment. The sensitivity is taken to be in the form of both mathematical modeling for the rheological behavior of blood and in the reconstruction procedure of an anatomically realistic geometry of a cerebral aneurysm obtained from in vivo rotational CTA scan.
2. Image reconstruction. The medical images used in this study were obtained in vivo from rotational CTA, with resulting voxel resolution of 0.4 mm on a 512 3 grid. The reconstruction of the virtual 3D geometry for numerical simulations consists of image segmentation, surface extraction, and finally surface smoothing. The initial surface triangulation is obtained using a constant threshold value for segmenting the image data, and a marching tetrahedra algorithm for extracting the 3D surface. Several other segmentation methods exist [6,2], however each methods yields a different geometry definition limited by the acquisition modality. The resulting virtual model of the vasculature is then prepared for the numerical simulations by identifying the regions of interest and removing secondary branches. Surface smoothing is required due to medical imaging noise and limited resolution that yields an initial surface definition that is not anatomically representative.
Smoothing is performed using the bi-Laplacian method [11]. Figure 1 shows the initial cerebral arterial geometry, with some secondary branches and a saccular aneurysm, as well as the region of interest that includes the aneurysm and is used in the numerical simulations.
It is apparent that there is a degree of uncertainty with respect to the model definition that stems from using in vivo data. This is due to the fact that different segmentation and reconstruction procedures will lead to virtual geometry variance, and thus to different numerical solutions of the flow field.
Here we focus the study of image reconstruction uncertainty by considering two levels of smoothing, with 300 and 20,000 iterations of the bi-Laplacian method, that corresponds to light and intense smoothing, respectively. For the rest of the paper, the two virtual model definitions are referred to as geometries G s (lightly smoothed) and G S (intensely smoothed). These intensities are representative of acceptable levels of smoothing, where the lightly smoothed case retains small-scale geometric features while the intensely smoothed geometry reduces these [11]. These geometry models, with the two different smoothing levels, are shown in Figure  2. The alterations to the surface are on average below 0.2 mm (1/2 voxel size) and within the uncertainty of the medical image resolution. The unconstrained smoothing is performed and the regions of highest curvature variation are the most affected. Here this corresponds to the distal region of the aneurysm neck, where the geometric difference is of the order of the voxel size.
Clearly, other approaches of performing a sensitivity analysis on the geometry can be sought, using different segmentation and reconstruction methods, as well as different levels of smoothing. Furthermore, different portions of the arterial system to include, for instance, the side branches, or a larger portion for the region of interest also should be considered. Nevertheless, the study presented here, using the light and intensely smoothed model definitions, as a means of introducing small virtual model definition perturbations, is a good compromise to the level of detail that a rigorous sensitivity analysis would require. The sensitivity of the geometry definition to within the medical image resolution is being considered, assuming that an optimal segmentation has been performed. It provides a preliminary indicative measure of the effect that geometric variability has on the resulting flow solution.
3. Fluid Equations. In this work, blood is modeled as an incompressible and isothermal fluid. The governing equations, corresponding to the conservation of mass and linear momentum, can be written as where ρ is the constant fluid density, u and p are the fluid velocity and pressure, and τ is the deviatoric or extra stress tensor. The rheological properties of blood are modeled through the specification of a constitutive relation for τ . In this study, the only non-Newtonian effect that is taken into account is the shear-thinning viscosity.
In this case, the extra stress is computed explicitly from the strain rate tensor through the relation Here, µ is the shear dependent viscosity andγ, the shear rate, is a scalar measure of the strain rate tensor, given bẏ The viscosity function, µ(γ), is fixed a priori by fitting experimental viscosity data and can be written in the general form where µ 0 and µ ∞ are the asymptotic viscosities at zero and infinite shear rate, respectively, and F (γ) is a continuous and monotonic function such that Viscosity model. The correct specification of the viscosity model is crucial to capture the correct rheological behavior of blood. When using models and material parameters coming from literature, it is important to know the specific conditions under which the viscosity was measured, namely, temperature, hematocrit and any specific pathological conditions of the blood donors. In the present work we use viscosity data obtained by Prof. M.V. Kameneva (Univ. Pittsburgh) for normal human blood at 23 o C, for an hematocrit of 40%. Using the estimates found in [16], we converted the data, obtaining realistic viscosity values at 37 o C) (body temperature). The model parameters shown in Table 1 were computed using a non-linear least squares fitting of the viscosity data. More specifically, for each viscosity model of the form µ(γ; Λ) shown in Table 1, we computed the parameter vector Λ = (λ 1 , · · · λ n ) solution of the minimization problem where K ⊂ R n is the set of admissible parameters and µ i are the measured values of viscosity, for eachγ i . In the numerical optimization procedure, we start by performing an unconstrained optimization, and then adding inequality constraints incrementally, until all parameters are within physiological ranges. For details on constrained numerical optimization see [19]. Oldroyd Table 1. Selected generalized Newtonian models for blood viscosity with the corresponding material constants. Constants were obtained by nonlinear least squares fit from viscosity data (Poi) with Ht = 40% and T = 37 o C.
The most striking fact that can be easily observed in Figure 3 is that the zero shear rate viscosity µ 0 can change considerably, depending on the viscosity model. In that figure we can observe that, although the Carreau and Cross models fit very well the experimental data, they present quite different values regarding the zero shear viscosity. This is related to the lack of viscosity data for very low shear rates which, in the present study, start at a minimum of 0.06 s −1 . In commonly used parameter sets for blood viscosity models, like the one proposed in [8], the asymptotic viscosities µ 0 and µ ∞ are fixed a priori and only the remaining parameters are considered in the numerical fitting. Although this technique removes the heterogeneity in the computed µ 0 , in absence of experimental data for very low shear rates there is no valid reason to fix the asymptotic viscosities.  Figure 3. Apparent viscosity as a function of shear rate for whole blood at Ht= 40 %, T = 37 o , obtained using a Contraves LS30 for low shear rates (γ ≤ 128s −1 ) and a CMSM (γ ≥ 300s −1 ).
In this study the Carreau viscosity model was adopted, with the set of parameters shown in Table 1. The reference constant viscosity considered for the Newtonian vs. non-Newtonian comparisons was µ = 0.04 Poi, corresponding to the average experimental viscosity in the rangeγ ∈ [6, 1000] s −1 . Another possible criteria to choose the Newtonian viscosity is to compute the viscosity value that would yield the same volume flow rate in an infinite cylinder of comparable radius. This criteria would give a reference value of µ = 0.0345 Poi.
Flow parameters. The fluid equations given by (1) must be endowed with proper initial and boundary conditions, corresponding to each particular physical problem. In the simulations presented in section 5, the vessel is assumed to be rigid and noslip conditions are imposed in the physical part of the boundary. The fluid starts from rest and is ramped up to a steady state, by imposing a velocity profile at the inflow section, that is parabolic with respect to the distance from the centroid of the section to its boundary. The velocity ramp is of the form: The chosen velocity profile corresponds, for t ≥ t ramp , to a physiological average value of the volume flow rate in the artery under study. More precisely, we use the reference value Q = 4cm 3 s −1 , taken from [20]. At the outflow section we are imposing a no-traction condition. Using as characteristic length and velocity, the average radius of the artery and the maximum velocity magnitude at the inflow section, respectively, the Reynolds number is Re = 250, in the Newtonian case, for a mean radius of 0.25 cm.
This choice of boundary conditions can be greatly improved, namely by considering a geometrical multiscale approach [9, Chapters 10 and 11], capable of taking into account the remaining circulatory system, as well as the cerebral auto-regulation mechanisms leading to changes in collateral flows in response to alterations in afferent blood pressure (see [18]).

Relevant hemodynamic indicators.
The importance of using mathematical models to obtain numerical simulations of blood flow in arteries is due to the link between hemodynamics and disease formation, such as atheroma and aneurysms. Although the correlation between the blood dynamics and disease is not fully comprehended, it is clear that a good understanding of the flow field helps in having a better knowledge on the disease. The correlators between hemodynamics and disease most commonly sought are the mechanical parameters on and near the wall, such as WSS and derived measures.
The WSS is the tangential stress that the fluid exherts on the wall, consisting on the tangential component of the stress tensor, which is given by σ = pI − τ , on the wall: W SS = σ n − (σ n · n)n = τ n − (τ n · n)n, where n is the outward normal to the wall surface, and σ n and τ n are the normal components of the stress and extra-stress tensors, respectively.
Amongst the derived measures of the WSS, are the WSS gradient (WSSG), the oscillatory shear index (OSI) and the gradient oscillatory number(GON). While the OSI and GON are wall measures related to unsteady flow, the WSSG is the spatial variation of the WSS: where τ ζ1 = µn · ∇(u · ζ 1 ) and τ ζ2 = µn · ∇(u · ζ 2 ), with ζ 1 and ζ 2 the perpendicular unit tangent vectors to the surface. It is commonly accepted that the regions of hight WSSG are related to disturbed flow patterns associated with the location of pathologies of the vessel wall, like myointimal hyperplasia or atheroma.
4. Numerical method. The mathematical model described in the previous section was discretized in space with stabilized P 1 − P 1 finite elements [4]. The time discretization was performed with a backward difference formula (BDF) of order 3, and all simulations where performed with an in house version of the freely available FEM software LifeV (see http://www.lifev.org). This requires the differential problem (1) to be written in the weak form. Being Ω the fluid domain, let us define the Hilbert spaces V = H 1 Γ d (Ω) and Q = L 2 (Ω), where Γ d denotes the portion of the boundary where Dirichlet conditions are imposed. Formally multiplying equations (1) by suitable test functions and integrating by parts, yields The finite element method consists in substituting the functional spaces V and Q by suitable sequences of finite dimensional subspaces {V h |h > 0} and {Q h |h > 0} in which we can compute an approximate solution (u h , p h ). A proper choice of these subspaces ensures that the problem is solvable in the finite dimensional spaces and the approximate solutions converge to the solution of (5) when h → 0. The convective term is linearized in a semi-implicit way, by considering that at time step t k+1 we have u k+1 · ∇u k+1 ≈ u k · ∇u k+1 . Finally, in order to fully linearize the discrete problem, we use the approximation The fully discrete problem then reads: and The numerical simulation of the previous equations can be carried out as in the constant viscosity case, having the same number of degrees of freedom. The computational effort is not exactly the same in both cases, since on one hand to solve (6) there is the need, at each time step, of computing the shear rate, and on the other hand the number of iterations to solve the linear system changes. Nevertheless the computational costs are of the same order. The BDF method can be memory consuming because the solutions in the previous three time steps need to be stored.

5.
Numerical results and discussion. A total of four cases were studied with respect to uncertainty in both the geometry and constitutive model of the blood. In particular, geometries G s and G S , obtained from different smoothing intensities were considered, each with two different rheological models, Newtonian and Carreau. From here onwards the superscripts N and C will be used to refer to the Newtonian and Carreau models, respectively.
Numerical simulations were performed using graded meshes with element size of 0.16 mm within the aneurysm. This amounts to using around 0.85 M tetrahedral elements in each geometry.
Steady state solutions were sought, however an instability associated to a small recirculation region at the distal portion of the aneurysm neck meant that this was in practice not achieved. The instability did not affect the flow upstream to this region, including within the aneurysm, while downstream a greater influence was noted. The simulations discussed here are taken at the same time step, at a total of 10 s after a steady inflow velocity was imposed, hence at time T = t ramp + 10 s.
The general flow patterns within the aneurysm indicate separated flow in the overhang regions at the proximal and distal portions. Since the aneurysm is located on the outer bend of the artery, part of the core flow enters the aneurysm naturally, impacting the upper side wall and forming two counter rotating vortical structures. These spread out to fill the separated regions and move to the outflow zones.
In Figure 4 the particle tracks for each case are depicted. These tracks are seeded at the neck of the aneurysm, at the same locations for all cases. It is evident that both rheological model and geometry have a considerable impact on the resulting velocity field. Moreover, from cross sections of the velocity magnitude shown in Figure 5 an emphatic difference in the flow field is perceived. We can also observe that for both solutions using G s (hence G C s and G N s ), the bulk of the velocity flux along the artery axis occurs in what the artery would have conceivably been. However, the aneurysm contains high velocity swirling motion that is driven from the neck of the aneurysm. For G S a greater portion of flow with a higher velocity sweeps into the aneurysm, however this is seen to exit more readily and the regions of swirling and separated flow tend to be more contained. Figure 6 displays the WSS distribution. Within the aneurysm the WSS is low, with values typically below 2 Pa and on average around 1.5 Pa. The highest values, around 8 Pa, are found at the distal portion of the neck. For this patient specific study, the geometry uncertainty given by the surface smoothing has an influence on the WSS comparable to that of different rheological models.
In order to quantify the differences in WSS given by these uncertainties, illustrated in Figure 7, the mean values are computed (for the aneurysm region shown only) and found to be While these values show that the geometrical uncertainty has a greater impact on the WSS with respect to the rheological modeling, both have a significant magnitude. However, within the aneurysm, the peak differences are found to be more pronounced when comparing different models, while the largest differences are found to be at the aneurysm neck due to the geometric variations.  Figure 7. WSS difference [dyn/cm 2 ] between (a) the Carreau and the Newtonian models on G s (b) the Carreau and the Newtonian models on G S (c) the Newtonian model on G s and G S (d) the Carreau model on G s and G S .
Differences in the WSS do not affect the location of high WSSG, and these are consistently located at the distal portion of the aneurysm neck, as reported in Figure 8. With respect to a physiological impact, it is evident that the WSSG is not emphatically affected by the modeling uncertainties for this case study, and may more generally constitute a more robust parameter. For the patient case studied, no high WSSG is seen within the aneurysm with an average value of 7 Pacm.
Regarding the viscosity distribution obtained for the Carreau model, the average values of µ S = 0.0337 Poi and µ s = 0.0334 Poi were observed in G S and G s , respectively. It should be noted that these values are smaller than the constant Newtonian viscosity µ = 0.04 Poi, chosen as an average high shear viscosity over the rangeγ ∈ [6, 1000] s −1 . This indicates that the alternative method of computing the equivalent Newtonian viscosity proposed in section 3, giving µ eq = 0.0345 Poi is in fact more accurate in this specific situation. However, the value µ = 0.04 Poi used for the Newtonian viscosity can be found, with different justifications, in several works concerning hemodynamical applications [11]. In fact, the choice of the Newtonian viscosity, not addressed in this work fully, is yet another source of uncertainty. It was also observed that the higher viscosity values, near µ = 0.04 Poi, are only found close to the inflow section, before the first bend, where the fluid starts exhibiting strong secondary flows.
Conclusion. The study presented is a preliminary work to examine the effects of uncertainty on numerical simulations of cerebral aneurysms from in vivo medical data. The sensitivity to the computational pipeline of performing hemodynamic simulations was introduced in the form of variations in the geometry and rheological model. Significant differences were observed in solutions for both cases, being of the order of 5% with respect to the rheological model, and 15% for the geometrical variations. This indicates that since blood rheology includes non-Newtonian effects, these should not be neglected in the mathematical modeling. Furthermore, sensitivity to the geometry has a greater, though comparable, influence on the computed WSS.
Naturally, the results obtained in this work are patient specific and do not necessarily extend to general situations. In order to be validated in a more general framework, a comprehensive study must be performed in a large representative medical imaging database. Moreover, the analysis of a steady state solution is not physiological and the fully unsteady problem should be considered, using a realistic waveform. The prominent differences in computed flow field detected for steady flow are expected to be present in the time dependent simulations, however the extent cannot be inferred from these results.
Finally, it must be noted that uncertainty can come from other sources. These include: the use of different viscosity models based on distinct fitting procedures of the experimental data; the use of different structural models for the vessel wall and the inclusion of fluid-structure interaction; the inclusion of smaller side branches that are usually cut-off, inflow and outflow boundary conditions. Kameneva from the University of Pittsburgh, USA, for providing experimental viscosity data.