Shear-thinning effects of hemodynamics in patient-specific cerebral aneurysms. Mathematical Biosciences and Engineering

. Two diﬀerent generalized Newtonian mathematical models for blood ﬂow, derived for the same experimental data, are compared, together with the Newtonian model, in three diﬀerent anatomically realistic geometries of saccular cerebral aneurysms, obtained from rotational CTA. The geometries diﬀer in size of the aneurysm and the existence or not of side branches within the aneurysm. Results show that the diﬀerences between the two generalized New- tonian mathematical models are smaller than the diﬀerences between these and the Newtonian solution, in both steady and unsteady simulations.

1. Introduction. Cerebral aneurysms are localized pathological dilations of the brain vessels, due to the weakening of the wall media. They potentially result in the sudden death or morbidity of the patient, through rupture and massive hemorrhage, or by pressing adjacent brain tissue. In the majority of the cases intracranial aneurysms remain asymptomatic until rupture. Based on autopsy and angiography data it is estimated that up to 6% of the population harbors one or more intracranial aneurysm [23,16]. They often occur at the apices of bifurcations and outer bends of curved arterial segments [19,30,8,12] suggesting, together with the numerous studies found in literature [23,5,7,22,9,3,24,26,20], their strong correlation with local hemodynamics, among other factors that include genetic predisposition, lifestyle and extra-vascular structure. Consequently, the medical community is increasingly demanding for scientifically rigorous and quantitive studies of blood flow dynamics within patient-specific aneurysms, driving the development and use of mathematical models and numerical methods for the computational hemodynamic simulations [23,22,7,27,13,6]. Nowadays, computer simulations already play a decisive role in better understanding the mechanisms of initiation and progression of aneurysms [23,22,7,27,13,6], providing valuable information on the details of hemodynamic characteristics believed to be related with the weakening of the wall [1,16,15], such as wall shear stress (WSS), oscillatory shear index (OSI), elevated pressure or area of flow impingement. Nevertheless, the numerical simulations are prone to a number of uncertainties [28,11,23,22,9] that should be quantified when using their outcome to analyze aneurysm formation and development, or to draw clinical conclusions such as the risk of aneurysm rupture or the success of surgical treatment. Indeed, despite the increasing feasibility and capability of performing numerical simulations, these hold intrinsic errors associated to in vivo medical image data and techniques for geometry reconstruction, and translation of physical observations to mathematical models.
The choice of the fluid mathematical model is mostly related with red blood cell (RBC) behavior. The mechanical response of blood, hence its rheological characteristics, largely depend on the state of its RBC, that can range from three-dimensional microstructures, at low shear rates, to uniformly dispersed individual cells, aligned with the flow field, at high shear rates [25]. This is relevant because RBCs have a high volume concentration in whole blood (hematocrit H t ≈ 45%). Besides RBC, blood is formed by white blood cells and platelets, suspended in an aqueous polymer solution, the plasma [25]. While plasma is nearly Newtonian, whole blood exhibits non-Newtonian properties, such as a varying shear-thinning viscosity. However, throughout most of the arterial system of healthy individuals, the RBC are dispersed, and the viscosity of blood can be considered effectively constant, hence it is sufficient to model blood as a Newtonian fluid [25,17,7,20]. Many authors advocate this argument, nevertheless this simplifying assumption may not be valid in some diseased states where the vascular geometry is altered and induces regions of slow recirculation, such as in aneurysms or downstream of stenoses. In these cases the non-Newtonian behavior may become relevant and constitutive models with shear-thinning viscosity or yield-stress should be applied [25,21].
Several studies can be found in literature regarding the comparison between Newtonian and non-Newtonian blood models, both in idealized and patient-specific scenarios [22,23,9,5]. Although the error bounds due to the fluid model might have a lower impact with respect to the one's related with the image reconstruction [17], it can be shown that the variability of the blood viscosity can lead to differences in the wall shear stress and other important hemodynamic indexes with clinical significance [23,9,5], indicating that non-Newtonian models should be analyzed. However, very little can be found in the literature regarding the comparison between different non-Newtonian models. In [21] the authors compare the Casson and Bingham's models in a specific pathological situation, with no direct relation with aneurysms. While the uncertainty quantification in the fluid mathematical model is mostly carried out through Newtonian vs non-Newtonian studies [23,9,5], there exist several non-Newtonian constitutive models considered appropriate to describe blood flow [9,5,21,25]. Hence, the choice of a certain non-Newtonian constitutive model also constitutes a source of uncertainty with respect to the computed flow field, and should be further analyzed. The present work focuses on the comparison of two non-Newtonian fluid models commonly used in blood flow simulations: the Carreau and the Cross shear-thinning models [25,9], together with the standard Newtonian one.
The fluid dynamics strongly depends on the geometry surface definition, making it essential to use anatomically realistic geometries to study aneurysms. Indeed, patient-specific geometries can be highly complex, with several side branches, non-planarity and sharp bends, leading to a disturbed flow field that is usually over simplified in idealized geometries [10]. The size of the aneurysm also plays an important role in the existence of recirculation areas, where the use of proper non-Newtonian fluid models might be crucial. Hence, three different anatomically realistic geometries, reconstructed from data obtained from in vivo rotational CTA scan, are studied in this work: a small aneurysm, a big aneurysm and an aneurysm with a side branch within the saccular region.
The paper is organized as follows. Section 2 is devoted to detailing the patientspecific geometries and their reconstruction. In Section 3 the mathematical models, including the non-Newtonian constitutive relations, are described in detail. In Section 4 the numerical methods and setup are introduced. The numerical simulations are carried out by means of the Finite Volume Method, using OpenFoam (www.openfoam.com). Steady and unsteady simulations are performed and analyzed. In Section 5 the numerical results are presented and discussed. Finally, the conclusions are drawn in Section 6.
2. Patient data sets. The medical images used in this study were obtained from in vivo rotational CTA as part of clinical procedures of cerebral aneurysm exclusion by coiling. The patient data sets are used and termed in this study 'big', 'branch' and 'small', due to their different characteristics (see figure 1). The chosen data sets reflect a variety of shape, location and presence of large side-branches within the aneurysm, that can be considered typical cases. The portion of the aneurysm location along the cerebral arterial system varies in each case. The voxel resolution was 0.4122 mm on a 512 3 grid for the 'big' and 'branch' cases, while 'small' was of 0.8244 mm on a 256 3 grid. The reconstruction of the Three-dimensional geometry for numerical simulations consists of image segmentation through a constant threshold approach, surface extraction using linear interpolation from the voxel grey-scale values, and finally surface smoothing based on the mesh connectivity [9,11]. The resulting surface models of the arteries are prepared for the numerical simulations by identifying the regions of interest, removing secondary branches, and creating an unstructured tetrahedral volumetric mesh. The meshes consisted in approximately 500K to 800K tetrahedral elements for the different cases studied. A mesh sensitivity analysis was performed in previous work [9]. 3. Blood rheology. The circulation of blood in the human cardiovascular system depends on the driving force of the heart, on the geometry and mechanical properties of the vascular system, and also on the mechanical properties of blood. Here we will focus on hemorheology, the study of flow properties of blood and its elements. For a comprehensive description of blood components and properties, we refer to Robertson, Sequeira and Kameneva [25] and references therein.
It is widely accepted that hemodynamic factors like flow separation, recirculation, low and oscillatory wall shear stress, all play a role in the localization and development of vascular diseases. For this reason, the development of mathematical models and numerical simulations in the vascular system can lead to a better understanding of the factors that trigger such diseases as well as to improve diagnosis and therapeutic planning.
In order to derive proper constitutive models for blood we must explore the connection between the composition of blood and its physical properties. As already mentioned, whole blood roughly consists of a suspension of cellular elements like red blood cells, white blood cells and platelets in plasma. While plasma is usually considered a Newtonian fluid, whole blood exhibits non-Newtonian properties like shear-thinning, as well as viscoelasticity, thixotropy and yield-stress, mainly due to RBCs behavior and to their high volume concentration in whole blood.
Shear-thinning. The main non-Newtonian property of blood, confirmed by very simple rheometric experiments, is its non constant viscosity. Blood viscosity is a decreasing function of shear rate. When blood is at rest or at low shear rates the RBCs aggregate in rod-shaped structures called rouleaux. These structures create an increased resistance to motion and cause blood to exhibit a higher apparent viscosity (Figure 2  As shear rate increases, these large aggregates breakdown into smaller aggregates and eventually into individualized RBCs (Figure 2 (b)), leading to a decrease in viscosity. As shear rate increases even further, the RBCs become deformed and align with the main flow, decreasing the viscosity. It has been shown that this process is reversible, i.e. if the shear rate is decreased the RBC will return to an undeformed state and eventually the rouleaux are formed again. The dynamics of viscosity variation is in fact more complex because viscosity does nor change instantly with shear rate. When shear rate is stepped up or down, it takes some time for the formation and breakdown of RBCs aggregates, and consequently for the value of viscosity, to stabilize. This effect is known as thixotropy and is commonly neglected in most rheological models. In fact, the amount of time necessary for the rouleaux to form, together with the small dimension of regions of consistently low shear rate, prevents their existence in most parts of the arterial system of healthy individuals and, for this reason, it is acceptable to consider a constant value for viscosity. However, if in some part of the arterial system the stability of the RBCs aggregates is augmented, for example in the case of saccular aneurysms where a large recirculation area is created, the viscosity variation and consequently shearthinning models must be taken into consideration.
Most shear-thinning models just amount to prescribe a functional dependence between shear rate and viscosity. The viscosity function µ(γ) is fixed a priori by fitting experimental 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 Notice that, given the properties of F (γ), the viscosity function (1) is bounded by µ 0 from above, and µ ∞ from below. These limit viscosities have solely a mathematical meaning since for very high shear rates the RBCs would be destroyed and for very low shear rates it is experimentally impossible to accurately measure the viscosity. The experimental viscosity data used in this work were obtained by Prof. M.V. Kameneva (Univ. Pittsburgh) for normal human blood at 23 o C and for an hematocrit (Ht) of 40%. The measurements were carried out with a Contraves LS30 machine for the values of shear rate under 128 s −1 , and with a Contraves MSM machine for the values of shear rate over 300 s −1 . Using the estimates found in [18], these data were converted to viscosity values at body temperature, i.e. 37 o C (see [9] for details).
In Table 1 are listed the well known Carreau and Cross viscosity functions, together with the correspondent parameter values, which were obtained by means of a nonlinear least squares fit from the viscosity data [9].

Model
Viscosity Model Model constants for blood Table 1. Selected generalized Newtonian models for blood viscosity with the corresponding material constants. Parameter constants were obtained by nonlinear least squares fit from viscosity data, with Ht = 40% and T = 37 o C.
In Figure 3 the viscosity models of Table 1 are plotted. It can be seen that the zero shear rate viscosity µ 0 changes considerably between models, for the same viscosity data. In fact, even if when different viscosity models fit very well the experimental data, they can 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 .  A possible fix for this increased uncertainty in the choice of the zero-shear viscosity is to use additional experimental or numerical data. More precisely, it is possible to construct semi-analytical solutions for the flow of pseudo-plastic fluids in pipes, under a constant pressure gradient (see [14]). This allows to choose µ 0 in such a way that a certain experimentally measured peak velocity is yielded by our viscosity model. In this way, the value of µ 0 , initially calculated from the least squares fit of the viscosity data, can be tuned to yield the correct velocity for a certain range of pressure gradients. This procedure is therefore dependent on the expected flow parameters. The maximum velocity for flow in a straight tube, for different pressure gradients, is displayed in Figure 4, for the Carreau model with the initial parameters presented in Table 1.
Fluid equations. Given the considerations on blood rheology presented so far, we model blood flow as an incompressible and isothermal fluid with varying sheardependent viscosity. These fluids are commonly known as generalized Newtonian fluids and the conservation of mass and linear momentum equations can be written as where ρ is the constant fluid density, the variables u, p are the fluid velocity and pressure, and τ is the extra stress tensor which, for this class of models is given explicitly by τ = µγ(∇u + ∇u T ). The variable viscosity is a function of the shear rate magnitude, a scalar measure of the strain rate tensor given by   This particular form of the extra stress tensor, being symmetric, is also compatible with the conservation of angular momentum. This set of equations must be endowed with appropriate boundary and initial conditions, depending on the particular geometry and situation to be modeled. Finally, thermodynamical effects are neglected since body temperature can be considered constant.

Numerical methods and flow parameters.
In the scope of this study both steady-state and unsteady problems were considered. For the steady-state simulations the set of equations described in the previous section was discretized using a standard finite volume method with linear Gauss integration. In the case of the unsteady simulations a pressure implicit operator splitting was used to discretize in time, while the space discretization was similar to the steady case. The simulations were carried out with the OpenFoam CFD software. In the steady-state, the boundary conditions consisted on imposing a flat velocity profile at the inlet of each geometry, in such a way that the Reynolds number Re = 250 is based on the hydraulic diameter. The average diameters in the undamaged part of the vessels are 4.8 mm, 5.16 mm and 3.35 mm for the geometries named as 'big', 'branch' and 'small', respectively. The no-slip condition was imposed at the vessel wall, and zero pressure at the outlets. The remaining conditions were set as homogeneous Neumann. In the unsteady simulations the inlet velocity was modulated in order to obtain a physiological volume flow rate.
Although much more elaborate options could be considered for the boundary conditions, namely using reduced mathematical models [22] to account for the global circulation, we decided to simplify this particular aspect of the work. The inlet sections are chosen to be far from the aneurysm, so that the flow can develop and the effect of choice of inflow condition is reduced.

5.
Numerical results and discussion.  It is possible to conclude that the three rheological models are qualitatively similar in the three geometries, although relevant quantitative differences are recorded. Figures 5-7 display the WSS distributions. The WSS is on the whole lower within the aneurysm than in the artery. The 'small' geometry exhibits the smallest velocity magnitude within the aneurysm, while the 'big' geometry has the largest due to its location at the outer bend of the artery and its size that facilitates the entrance of the flow, while in the 'branch' geometry the larger velocity is due to flow circulating around the dome and exiting through the secondary branch. Quantitatively, an average WSS within the aneurysm is approximately 3P a and the maximum WSS greater than 10P a is located at the distal part of the neck. The observed relative difference with respect to the Carreau model are as follows:     It is to be noted that even though the difference between the shear-thinning models is lower than between them and the Newtonian one, the choice of viscosity model impacts the results. This difference is noticeable mostly inside the aneurysm, where the low average shear rate enhances the role of µ 0 , which can be quite different from one shear-thinning model to another.
In order to observe the flow region within the aneurysm, a slice perpendicular to the artery is extracted. In this analysis we focus on the 'branch' geometry. Figure 8 shows the velocity magnitude and Figure 9 shows the dynamic viscosity magnitude in the slice. It can be noted that the differences in velocity and viscosity are less prominent than for the WSS, however the greatest differences are clearly identified to be at regions of both higher and lower velocity gradients. This indicates that in the regions of both low and high shear rates the computed solution will differ more noticeably. The presence of the secondary branch within the aneurysm causes the flow to impinge on the distal region of the aneurysm neck and in part flow along the top portion of the dome and exit in the secondary branch. This causes a region of As a final comparison we observe the integrated effects of the different rheological models for blood by computing the pressure drop between the inlet and outlet sections along the main artery. These pressure drops are given as follows: It can be seen that the Newtonian model has the highest pressure drop consistently, while the Carreau and Cross models behave similarly but with moderately higher pressure drop in the latter case. This indicates that the Newtonian model tends to have higher viscous forces than the shear-thinning models used here. 5.2. Unsteady simulations. The same nine cases reported for the steady simulations were considered in an unsteady physiological regime. The simulations were run for a total of twenty cardiac cycles, which was considered sufficient to reach a periodic behavior of the relevant hemodynamic parameters. The reported results correspond only to the last cardiac cycle.
All results are qualitatively similar to those obtained in the previous section but, depending on the geometry, the agreement between the steady solution and the averaged unsteady solution can vary considerably. In Figure 10 we plot, for each geometry and with a chosen model, the difference between steady and timeaveraged solutions, as well as the standard deviation for the unsteady solution. For each geometry, similar differences between steady and time-averaged solutions are observed. Results for only one model in each geometry are therefore shown. Considering the WSS magnitude obtained for the steady-state simulations, we observe that, for the 'branch' and 'small' geometries, there is a considerable gap between the WSS values obtained from the steady and unsteady simulations, while for the 'big' geometry the steady solution correctly represents the time-averaged unsteady solution.
The analysis of standard deviation in the time evolution of WSS gives useful information that allows to identify regions of higher WSS oscillation. In all cases under consideration this measure attains its lower values inside the aneurysm and higher values in or near the neck. In the case of the 'small' geometry the region of higher WSS standard deviation is most likely related to local surface curvature and not to the presence of the aneurysm.
With the purpose of illustrating the differences between the three models in a synthetic manner we present the time evolution of WSS in selected locations of the three geometries, namely in the region of the neck distal to inflow sections, see Figure 1 for the location of these points. The time evolution of velocity in a selected cell inside the aneurysm, in the case of the geometry labeled as 'small', is also shown in figures 11. 12. Here we represent the time evolution of the computed velocity field for the three rheological models, in the case of the 'small' aneurysm. The differences are significant, ranging from 6% to 13% when comparing the Newtonian with the non-Newtonian models, and around 2% when comparing the two non-Newtonian models. These results are compatible with the steady-state analysis. Figure 13 shows the time evolution of the WSS in a selected location in each geometry, as well as the relative difference between models. It is clear from these images that there is a significant difference between the WSS values reported by each model in absolute an relative terms along the cardiac cycle, reaching a maximum in peak systole and diastole. As expected, the difference between the Carreau and Cross models is smaller than between either of them and the Newtonian model. The magnitude of these differences varies between geometries and during the cardiac cycle, it is almost constant in the 'branch' geometry and presents higher variation in the other geometries. 6. Conclusions. Three patient-specific geometries of saccular cerebral aneurysms have been used to study the effects of shear-thinning viscosity models for blood on the computed flow field. The Carreau and Cross non-Newtonian models are used since these are commonly adopted to describe blood flow behaviour and the SHEAR-THINNING EFFECTS IN PATIENT-SPECIFIC CEREBRAL ANEURYSMS 13 Figure 10.
WSS difference [Pa] between steady and timeaveraged solutions (left) and standard deviation of the computed unsteady WSS (right) for the 'branch' geometry with the Cross model (top row), 'big' geometry with the Carreau model (middle row) and for the 'small' geometry with the Newtonian model. coefficients to these models are obtained from a common experimental data set. The major source of differences between these non-Newtonian models resides in relatively large discrepancies in the asymptotic zero-shear viscosity value, µ 0 , caused by the lack of experimental data for very low shear rates. This difference does not  have a physical meaning but rather is a consequence of data fitting. Both steady state and unsteady flow simulations highlight noticeable differences between the Newtonian and the shear-thinning models, however, differences between the non-Newtonian solutions are also visible. Overall, the wall shear stress indicates that approximately 10% variation can be present between the Carreau and Newtonian models, while approximately 4% is seen between the shear-thinning models. These differences are to be understood as average values, as they are not uniform in time and space. Cross-sectional slice through the aneurysm indicates similar trends and identifies the regions of greatest differences between the models as due to both the high and low shear rates. Pressure drops across the main arterial vessel show that the Newtonian model exhibits the greatest viscous drag. The three geometries considered present great differences in shape, size, location and characteristics. Even if this can be considered a limited set of geometries, it is nevertheless varied while the results indicate similarities in the effects of changing  the rheological model. The obtained results can help to understand the impact of the choice of the shear-thinning model to capture the rheological behavior of blood, as well as the possible uncertainty and variability when fitting different shear-thinning models to the same experimental data set, or opting for a Newtonian assumption. Our conclusions are limited by the lack of in vivo or in vitro measurements in these geometries, that could indicate which model is preferable under these flow conditions.