Search tips
Search criteria 


Logo of springeropenLink to Publisher's site
Annals of Biomedical Engineering
Ann Biomed Eng. 2010 October; 38(10): 3058–3069.
Published online 2010 May 25. doi:  10.1007/s10439-010-0065-8
PMCID: PMC2940011

Flow Residence Time and Regions of Intraluminal Thrombus Deposition in Intracranial Aneurysms


Thrombus formation in intracranial aneurysms, while sometimes stabilizing lesion growth, can present additional risk of thrombo-embolism. The role of hemodynamics in the progression of aneurysmal disease can be elucidated by patient-specific computational modeling. In our previous work, patient-specific computational fluid dynamics (CFD) models were constructed from MRI data for three patients who had fusiform basilar aneurysms that were thrombus-free and then proceeded to develop intraluminal thrombus. In this study, we investigated the effect of increased flow residence time (RT) by modeling passive scalar advection in the same aneurysmal geometries. Non-Newtonian pulsatile flow simulations were carried out in base-line geometries and a new postprocessing technique, referred to as “virtual ink” and based on the passive scalar distribution maps, was used to visualize the flow and estimate the flow RT. The virtual ink technique clearly depicted regions of flow separation. The flow RT at different locations adjacent to aneurysmal walls was calculated as the time the virtual ink scalar remained above a threshold value. The RT values obtained in different areas were then correlated with the location of intra-aneurysmal thrombus observed at a follow-up MR study. For each patient, the wall shear stress (WSS) distribution was also obtained from CFD simulations and correlated with thrombus location. The correlation analysis determined a significant relationship between regions where CFD predicted either an increased RT or low WSS and the regions where thrombus deposition was observed to occur in vivo. A model including both low WSS and increased RT predicted thrombus-prone regions significantly better than the models with RT or WSS alone.

Keywords: Intra-aneurysmal thrombus, Computational fluid dynamics, Flow separation regions, Magnetic resonance imaging


In addition to various biological factors, local hemodynamics is thought to affect the progression of aneurysms. State of the art CFD methods are capable of predicting the flow fields in patient-specific vascular geometries2,9,23,24 and the actual role of hemodynamics in aneurysmal disease progression can then be determined by correlating hemodynamic parameters obtained from patient-specific models (either numerical or experimental) with aneurysmal changes observed in vivo. Hemodynamic forces acting on the arterial wall may trigger mechanoregulated wall remodeling resulting in lesion enlargement.9 Recent studies have reported a correlation between aneurysm growth and regions of low wall shear stress (WSS).1,12,22 Low WSS has also been proposed as a predictor of aneurysm rupture.11 In the course of aneurysmal disease progression, intracranial aneurysms may not only grow in size but, in some cases, can also develop intra-aneurysmal thrombus. Although thrombus formation can sometimes stabilize lesion growth, it can also present an additional risk of thrombo-embolism.27 In our previous work, a strong similarity was found between the location of intra-aneurysmal regions predicted by CFD to have low velocities and shear stresses and regions where thrombus deposition was observed to occur in vivo in the follow-up MR studies.19 In this study, we extended our previous work to investigate the effect of the flow residence time (RT) on thrombus formation.

The complicated process of thrombus formation is initiated by aggregation of platelets. A few studies have suggested that platelet aggregation is determined by a combination of factors, which include abnormally high shear stresses activating the platelets and the time of exposure to such stresses.3,10,13,17 Three patients, considered in our study, had severe reduction of the flow through one of the supplying vertebral arteries, either due to stenosis or surgical occlusion. This resulted in increased flow through collateral vertebral with consequent increase of the shear stress proximal to the aneurysm. Activated platelets are carried by the flow and may accumulate in recirculation and stagnation zones. Low WSS observed in these regions decrease anti-thrombotic factor production from the endothelial cells.16 Increased RT in the near-wall regions with low WSS is also likely to facilitate platelet aggregation and adhesion of formed blood clots to the vessel wall.

In order to model the transport of platelets, it is necessary to solve for the background velocity field and then calculate particle trajectories.14,15 Particle residence time (PRT) is typically considered in the Lagrangian approach, where a large number of particles are seeded in the flow, their individual trajectories are computed, and the time spent by a particle in the domain is calculated. In addition to providing a useful method for visualization of the flow patterns, PRT can be used to assess the transport and mixing for a given flow. A powerful new method to characterize flow kinematics, as demonstrated by Shadden et al.,20 is to compute a Finite-Time Lyapunov Exponent (FTLE) field that reveals Lagrangian Coherent Structures (LCS). An LCS analysis of the flow in a patient-specific abdominal aortic aneurysm is presented in Shadden and Taylor.21 The disadvantage of the Lagrangian approach is the high computational cost of tracking a large number of particles in unsteady flows. In order to estimate PRT in complex flows, it is important to ensure that the particles enter flow separation regions. This may require very large number of seeded particles, which makes their tracking computationally expensive. An additional concern is the computational mesh resolution. As reported in the study of Tambasco and Steinman,25 “impracticably high mesh resolutions are required to resolve the complex flow fields for Lagrangian studies in these more realistic models,” since the path dependent quantities are “very sensitive to the quality of a velocity field.” The Eulerian approach offers an opportunity to estimate PRT with a reduced computation time. In order to obtain an Eulerian map of the flow RT, it is necessary to estimate the average time particles spend at a given site. A volumetric residence time (VRT) concept was proposed in Kunov et al.14 and later extended in Tambasco and Steinman25 for three-dimensions. In this approach, small Lagrangian fluid volumes, each assumed to contain a uniform concentration of particles, were tracked as they advected by a steady background flow field. The VRT for each element of the computational domain was then computed by adding the total time each particle resided in the element and dividing the obtained sum by the concentration of particles in the volume (within a normalization factor). This method, while quantifying the time each particle stays in each location, still requires computationally expensive Lagrangian particle tracking and a very high resolution mesh. An alternative method is to obtain an Eulerian residence time map by solving an additional RT equation, as proposed by Haimes.7 This equation models the advection of a passive scalar that is injected into the computational domain at the inlets. This method, while being simpler than the Lagrangian approach, provides RT information for every cell of the computational mesh.

In the current study, we have studied the RT distribution in three patient-specific basilar aneurysms with known regions of thrombus deposition using an Eulerian approach. The luminal geometries at base-line and following thrombus formation were obtained from MR angiography. Numerical modeling of the flow in these aneurysms was carried out and flow RT maps were obtained by simulating passive scalar advection. The RT was computed for a number of small sampling volumes placed regularly over the aneurysmal walls. In addition, WSS values were computed at the wall surfaces of these sampling volumes. The WSS values were then averaged over the cardiac cycle. The RT and WSS values computed for sampling volumes placed in regions that remained thrombus-free were compared with RT and WSS for sampling volumes placed in regions that were later observed to clot.


Patients with intracranial aneurysms were recruited to this study using institutionally approved patient consent. Numerical simulations were carried out to model the flow in three basilar aneurysms where intra-aneurysmal thrombus had formed. At base-line, the aneurysms were thrombus-free and flow was modeled in those geometries. The CFD results were correlated with regions where thrombus deposition was observed to occur on MR images at the follow-up studies. Patient-specific computational models constructed previously and described in Rayz et al.19 were used in this study to obtain the flow fields and calculate the shear stresses and RT. For the sake of completeness, we again describe here the methodology used for model construction and CFD simulation.

Construction of Patient-Specific CFD Models

Three patients had fusiform basilar artery aneurysms that were thrombus-free at initial presentation and that then subsequently proceeded to form intraluminal thrombus. Medical histories of the patients were described in Rayz et al.19 In two patients the thrombus formed during the natural progression of the disease; in the third case, intra-aneurysmal thrombus formed following surgical occlusion of one vertebral artery. High resolution contrast-enhanced MR angiography images (CE-MRA) of cerebral blood vessels obtained at base-line (thrombus-free) were used to obtain the luminal geometries for the CFD modeling. The maximum intensity projection (MIP) images of the base-line and follow-up aneurysmal geometries for these three patients are presented in Fig. 1.

Maximum intensity projection images from MRA studies of patients with basilar aneurysms—Top row: images of aneurysms at baseline; bottom row: images acquired after thrombus deposition

In-house software was used to form three-dimensional iso-surfaces from the MRA data. In each case, the 3D surface was co-registered with the two-dimensional MR slices and a threshold intensity value was selected to ensure that the segmented surface coincided with the luminal boundaries of the gray-scale MR image. The surface was then transferred into a 3D modeling software package, Rapidform (INUS Technology, Seoul, South Korea), where the volume of interest, which included the aneurysm with the major proximal and distal vessels, was selected. Minor vessels, including the small perforators branching off the aneurysm, and background noise were removed. The follow-up luminal geometries, obtained with CE-MRA following thrombus formation, were co-registered with the base-line geometries to determine the location of the regions of thrombus deposition. In order to ensure that differences were due to actual physiological changes rather than to differences in image processing, settings were selected such that reference cerebral vessels unaffected by aneurysmal disease appeared unchanged. The co-registration results showing the clotted-off regions for all three patients are presented in Fig. 2. All three models shown in the figure are oriented to match the MIP images obtained from in vivo MRA data and shown in Fig. 1. The surfaces constructed with Rapidform were divided into a number of rectangular patches for import into the CFD pre-processor Gambit (Fluent Inc., Lebanon, NH), to be used as the boundaries for the computational grid.

Co-registration of the base-line (gray) and follow-up (red) luminal geometries showing the regions occupied by thrombus. The models are oriented to match the maximum intensity projection (MIP) images obtained from in vivo MRA data (see Fig. 1 ...

Patient-specific flow boundary conditions were obtained with phase-contrast MR-velocimetry (PC-MRI) measurements. Flow velocities in the vertebral arteries proximal to the aneurysms were measured in planes transverse to the vessels. An electrocardiogram-triggered through-plane PC-MRI was acquired at 40-ms intervals through the cardiac cycle, providing approximately 15 time points through the cardiac cycle. These measurements were used to prescribe time-dependent inlet boundary conditions. Equal flow in the distal cerebral arteries was assumed for the outlet conditions, since the two hemispheres in the posterior circulation require comparable total volume flow. In spite of some simplifying assumptions made here, we have shown in our previous work that the flow fields obtained in these models correspond to in vivo measurements conducted with PC-MRI velocimetry.18

The governing incompressible Navier–Stokes equations were solved numerically with a finite-volume package, Fluent (Fluent Inc., Lebanon, NH). The modeling assumptions and numerical schemes used for discretization and solution of the flow equations were described in Rayz et al.19 For each patient, pulsatile flow simulations were conducted in the base-line geometry using the measured, time-dependent flow rates. Since our previous results showed that accounting for non-Newtonian blood behavior improved the agreement between the low velocity and thrombus deposition regions in two patients with larger aneurysms,19 all simulations in the current study were carried out for non-Newtonian flow. The non-Newtonian blood viscosity was modeled with the Herschel–Bulkley formula, which accounts for yield stress and shear thinning properties of blood. Pulsatile WSS distributions were obtained from the numerical solutions and used to analyze a relation between the WSS and thrombus deposition regions. A number of sampling patches were distributed over the aneurysmal surfaces, excluding the proximal and distal vessels, and for each patch a WSS value averaged over the cardiac cycle was calculated.

Residence Time Estimation: Virtual Ink Technique

In this work, we estimated the RT using a “virtual ink” technique. We virtually injected a passive scalar (“ink”) at the inlets and monitored the transport of the “ink” throughout the flow domain by solving an advection–diffusion equation:

equation M1

where equation M2 is the velocity vector, D is the diffusion coefficient, and C is the local concentration of the virtual ink. In this study, we set the diffusion coefficient D to zero, which ensures a pure advection of the virtual ink. This pure advection equation is identical to the level-set equation, which is widely used in numerical simulations to track the interface between different materials (liquid–solid, liquid–gas, and liquid–liquid).

The Eulerian concept of virtual ink is different from the Lagrangian particle tracking. There is a continuous field of ink values rather than a discrete distribution of particles. A propagation of the interface or front between the ink-painted and ink-free flow, however, is quite similar to an advection of a single injection of neutrally buoyant particles. A two-dimensional flow in a tube around a cylinder was used to qualitatively validate the virtual ink technique by comparing the positions of injected particles with the ink interface at the same time points (Fig. 3). It can be observed that the particles follow the same flow patterns as the ink. While the particles show the front propagation at discrete points, the ink technique highlights a continuous interface.

Propagation of virtual ink front along with advection of neutrally buoyant particles. All particles are injected at time zero, while the ink injection is started at time zero and continued for the whole time of the simulation. At different time steps ...

For the patient-specific models, the value of the scalar at the inlets was set to 1 during the injection of virtual ink and zero otherwise; therefore, the regions where the scalar value is 1 are filled with ink, while the regions where its value is zero are ink-free. It should be noted that, due to the Eulerian nature (which represents the concentration C as the volume average over the control volume) of this method and certain amount of numerical diffusion, the value of the scalar obtained in the numerical solution is not always 0 or 1 but can assume values in between. The flow regions with the scalar value above a certain threshold can be readily visualized with a color coding scheme. For all models, the flow domain was considered to be free of ink at time 0 (i.e., C = 0). The virtual ink was injected at the inlets starting from time 0, and carried by the flow through the aneurysmal geometries. The injection at both vertebral arteries ensured that no “clean” or ink-free flow mixed with the ink. It takes a number of cardiac cycles (depending on aneurysm size) to achieve a complete filling of the aneurysmal geometries. After the aneurysms were filled with the virtual ink, the simulations were continued until the aneurysms were completely flushed with “clean” or ink-free blood.

It should be noted that there was a one-way coupling between the Navier–Stokes and advection–diffusion equations, i.e., the distribution of ink at each time step was computed after velocity field solution had been obtained and therefore ink concentration had no effect on the velocity field. A second-order discretization schemes were used for numerical solution of the advection–diffusion equation. The values of ink at the inlets were provided at each time step and zero concentration of ink was maintained at the outlets. In addition, zero flux of ink was specified at the arterial walls.

The RT map could be obtained by calculating the difference between the times of the virtual ink arrival and departure at a given point. Starting the RT measurement from the time of uniform ink distribution avoids the unnecessary complications associated with waiting for ink to reach all regions of the model. The RT maps in the current study were computed as the time it takes for ink to wash out from the flow geometry, after establishing a fully ink-laden model and subsequently injecting ink-free blood at the model inlets for a number of cardiac cycles. Thus, flow RT calculations were started at the beginning of the ink-free blood injection, and stopped once the ink value in all regions dropped below the 0.02 threshold. This threshold value was chosen from practical considerations: the threshold should be close to zero, which corresponds to an ink-free domain, while not equal to zero as it would erroneously highlight the outlets as regions with increased flow RT. In order to quantify the RT near the walls, small three-dimensional sampling volumes were distributed over the aneurysm surface, excluding the proximal and distal vessels. Each of the sampling volumes had a side coinciding with the luminal surface and dimensions of 0.6 × 0.6 × 0.6 mm3 (Fig. 4). For each sampling volume, the RT is updated at every time step: if the virtual ink value remains higher than the threshold, the RT is increased by the time step interval. The total RT for a sampling volume is therefore equal to the time step interval multiplied by the number of time steps at which the scalar ink value was above the threshold.

Surface sampling with small boxes (0.6 × 0.6 × 0.6 mm3). The sampling volumes were distributed over the aneurysm surface and used to calculate flow residence time and WSS

Mesh Independence

In our previous study with the same computational models, we have verified mesh independence of the flow field solution. This, however, may not be sufficient to claim that this mesh resolution is also adequate for the scalar advection equation, particularly for the boundary layers near the walls where the flow RT was calculated. In order to ensure that the RT results are not affected by the mesh, we have carried out pulsatile simulations with the ink injection on three different meshes for patient 1: a crude mesh with a nominal density of 0.625 mm, a medium-density mesh with a nominal density of 0.313 mm, and a refined mesh with a density of 0.156 mm. The RT calculated on different meshes was compared for each sampling volume and the mean differences were found to be within 3%. While there was no observed improvement between the differences calculated between crude and medium, and medium and refined meshes, respectively, only five sampling volumes out of 280 had a difference higher than 10% (with none exceeding 18%). Therefore, the medium-density mesh used in the previous study was found adequate for the RT estimations and each simulation was made using a similar nominal cell size. Because of the different size of the aneurysms, this nominal density resulted in roughly 752, 315, and 893 thousand cells for patients 1, 2, and 3, respectively.

Statistical Analysis

For each patch, the mean value of the RT and the presence or absence of thrombus (binary variable) were recorded. A Mann–Whitney test was performed on RT between the thrombosed and nonthrombosed areas. A logistic regression with a clustering by patient method was performed to test the statistical significance of the association between RT and thrombus location. A similar analysis was performed on mean values of WSS. Finally, a Likelihood-ratio test was performed to compare the model RT alone and WSS alone vs. WSS and RT for the prediction of thrombus formation.


Flow Visualization

The propagation of the virtual ink through the vascular geometries provides important information about the flow fields and can serve as an excellent method for flow visualization. The images showing the washout of the virtual ink by the ink-free flow in all models are presented in Fig. 5. Each row of images shows temporal progression of the ink-free flow in a different aneurysm model. In the washout phase of the simulation, ink-free blood is injected into the inlets of the ink-filled models. In all three cases, the ink-free flow quickly propagates through the proximal vessels and is carried through the lesion by the high-velocity jet. The ink is quickly washed from the healthy part of the vascular geometries, while remaining significantly longer in the aneurismal bulges, where the flow velocities are particularly low. The total washout time for ink from the aneurysmal geometries was 40 s (approximately 33 cardiac cycles) for Patient 1; 14.7 s (17 cycles) for Patient 2; and more than 30 s (20 cycles) for Patient 3. In the last case, although the simulations were stopped before complete washout of the ink occurred, we were still able to determine that the total washout times for Patients 1 and 3 with larger aneurysms were roughly twice as long as the washout time for the smaller aneurysm of Patient 2. The images shown in Fig. 5 may create a misleading impression that the regions of increased RT extend throughout the aneurysmal bulges. In actuality, the high RT zones are confined to the regions adjacent to aneurysmal walls, as demonstrated in Fig. 6 showing cross-sectional views of the models. The images in Fig. 6 show an advanced stage of the ink-free flow injection, roughly corresponding to the images shown on the right side of Fig. 5, when most of the lesions are ink-free and the ink remains only in the near-wall areas of aneurysmal bulges. The sectional views in Fig. 6 show RT contours corresponding to the ink distribution 12 s after the start of ink-free flow injection for Patients 1 and 3 and only 2.5 s after the start of ink-free flow for Patient 2. The difference in time is due to the smaller size of the Patient 2 aneurysm. For all patients, there is a similarity between the regions of increased RT that are shown in the last column of Fig. 5 and the regions later observed to clot (shown in Fig. 2).

Tracer washout in all patients at different time instances. First row: Patient 1; second row: Patient 2; third row: Patient 3. The regions filled with virtual ink are shown in red. The ink is entrapped in the recirculating flow observed in the aneurysmal ...
Cross-sectional planes showing intra-aneurysmal regions near arterial walls with increased flow residence times (red). Patient 1: RT map 12 s after the start of ink-free flow; Patient 2: RT map 2.5 s after start of ink-free ...

Wall Shear Stress Distribution

The WSS distributions for all patients at peak systole are shown in Fig. 7. Low WSS values were observed everywhere at the aneurysmal walls, except the areas adjacent to high-velocity jets entering the lesions. In order to show the variation of WSS on the aneurysmal walls, the scale was magnified with the maximum value chosen to be 1.5 Pa (15 dynes/cm2), which is a normal physiological value observed in healthy arteries5,6 The WSS at the healthy walls of the proximal and distal vessels is, therefore, shown in red, although its values are actually within the normal range, rather than elevated.

Wall shear stress at systole. The scale is adjusted to show regions of abnormally low WSS; red color corresponds to a normal, physiological range (1.5 Pa)

The WSS is particularly low in the regions of flow recirculation, where slowly rotating vortices were observed throughout the cardiac cycle. The WSS in these regions remained abnormally low at all cardiac cycle phases. For all patients, the aneurysmal wall areas that were later observed to be covered by intraluminal thrombus are similar to the areas with CFD-predicted low WSS values.

Correlation Between Regions of Increased Residence Time, Regions of Low WSS, and Regions of Thrombus Deposition

Mean RT was 3.87 ± 4.48 s (range 0.40–28.51 s) in non-thrombosed area (335 patches) and 18.22 ± 11.21 s (range 0.63–40.13 s) in thrombosed area (456 patches) (p < 0.0001 on the Mann–Whitney test, Fig. 8). The logistic regression demonstrated that the odds of a location developing a thrombus increased by over 20-fold with each standard deviation increase in RT (95% CI: 13.3–25.9, p < 0.0001, r2 = 0.41).

Box and whisker plots of change in RT from area with and without thrombus. The boxes show the 25th and 75th percentile (interquartile) ranges. Median values (respectively, 2.2 and 16.9 for the nonthrombus and the thrombus groups) are shown as a horizontal ...

Mean WSS was 0.77 ± 0.57 Pa (range 0.03–3.13 Pa) in non-thrombosed area and 0.13 ± 0.14 Pa (range 0.03–1.56 s) in thrombosed area (p < 0.001). The logistic regression demonstrated an estimated change in thrombus location odd ratio of 81.7 per standard deviation change in WSS (95% CI: 40.7–163.8, p < 0.0001, r2 = 0.47).

The logistic regression model including both RT and WSS demonstrated a significantly higher prediction value than the model with RT or WSS alone (p < 0.0001, r2 = 0.57).


This study investigated the relation of two hemodynamic factors to intra-aneurysmal thrombus deposition. The WSSs and flow RTs were computed with CFD in base-line geometries for three intracranial aneurysm patients. A significant relation was determined between areas with low WSS and increased RT and regions that were observed to clot at the follow-up studies. In our previous study, we have shown that thrombus is likely to form in regions of both low velocities and low maximum shear stress—a parameter equal to the difference between the maximum and minimum eigenvalues of the shear stress tensor.19 In this study, we obtained a quantitative comparison between WSS and thrombus location and analyzed the effect of another important parameter—flow RT. The results show that either one of these two flow variables considered alone can be used to predict aneurysmal regions prone to thrombus formation, however a model including both WSS and RT demonstrated a significantly higher prediction value.

It is possible that a combination of several hemodynamic factors may contribute to the thrombus formation process. Slow flow in itself, like that observed in wall boundary layers in a healthy vessel, is not likely to cause thrombus deposition. Likewise, although low shear is observed in the core flow of a normal vessel it does not lead to thrombus deposition. An increased flow RT in a region with low shear stresses may promote clot formation and adhesion to arterial wall, since the forming clot is not washed away by the flow. The increased RT may also indicate compromised mass transfer between the flow and the vessel which may locally impact biological processes in the artery wall. A low WSS and high flow RT are likely to be related. In the work of Himburg et al.,8 the concept of “relative residence time” (RRT) is considered and expressed in terms of the oscillating shear index (OSI) and time-averaged shear stress magnitude. It is shown that near-wall RT at a particular site is inversely proportional to wall shear. Our results, based on the study of three patients, show that both low WSS and increased RT may facilitate thrombus deposition; however, it is not yet clear which flow factor is a better predictor of thrombus-prone regions. More patient-specific intra-aneurysmal thrombus cases should be considered in order to answer this question with confidence. Such cases are not easy to come by; however, we are currently monitoring 30 cerebral aneurysm patients and will expand our analysis if more thrombus deposition cases are observed.

Increased flow RT and low WSS, of course, cannot be considered as a necessary and sufficient condition for thrombus deposition. Two of the patients considered in the study (patients 1 and 2) were monitored with MR for several years prior to thrombus formation and during this time still remained thrombus-free, even though the observed flow fields were for the most part unchanged. Additionally, the flow RT in veins is likely to be much higher than that predicted in aneurysmal arteries, yet venous thrombus formation does not occur under normal conditions. Initiation of clotting within an aneurysm is likely influenced by biochemical processes such as inflammation, but once the clotting process is triggered, the thrombus is likely to form in the regions of flow separation characterized by low shear stresses and increased RT. The virtual ink technique used in the study is capable of showing flow separation regions. Thus, we are not proposing an entirely new hemodynamic mechanism linked to the regions of thrombus deposition, but rather a convenient method for identifying these regions in aneurysmal geometries.

In this study, we simulated a pure advection of a virtual ink. This technique is useful for visualization of flow separation zones, since these are the last to be reached by the ink as well as the last to become ink-free, once the ink is trapped in recirculating flow. Arterial flow, however, is pulsatile and three-dimensional which facilitates mixing and allows the virtual ink to eventually escape from the flow separation regions. It should be noted that ink trapped in the recirculation regions will eventually exit through the distal vessels, thus artificially elevating the RT downstream of these regions. In order to remove this effect from the RT calculations, the regions with ink concentration below a certain threshold were considered ink-free. The mixing of the ink still exiting from the recirculation regions with the large amount of “clean” flow coming from the inlets ensures that the ink concentration downstream of the aneurysms is below the threshold and does not affect the RT map. The value of the ink scalar in any given point can vary between 0 and 1, and it is an arbitrary decision on what value is taken as a threshold separating ink-free and ink-carrying flow. In this study, the threshold defining virtual ink washout was set to 0.02. This threshold value was chosen to be very close to zero, which corresponds to an ink-free domain. In practice, the threshold cannot be set equal to zero, since even after most of the ink is flushed and only a few pockets remain in the near-wall regions, there is still a very small amount of ink constantly leaving these pockets and exiting through the outlets. Thus, a zero threshold would erroneously highlight the outlets as regions with increased flow RT. We have found that the value of 0.02 was sufficient to ensure that the highlighted regions are only those where ink remains for a long time, and not just passes through on the way to exit. The contours of the regions with increased RT depend on the threshold values as well as on the time that has passed from the beginning of the washout. The delay time that provides contours of the virtual ink that are very similar to the region of thrombus deposition in one aneurysm model may be much longer than is needed to completely washout ink from a smaller geometry. Because of this difference in washout times, the images in Figs. 5 and and6,6, showing virtual ink contours for different aneurysms, were not obtained at the same time intervals from the start of the ink-free injection. As demonstrated in Fig. 6, the virtual ink remains in the near-wall regions of aneurysmal bulges, while the thrombus can extend somewhat deeper into the core of the flow. It is plausible that the thrombus forms in layers, with the initial layer adhering to the arterial wall in regions of increased flow RT and then gradually expanding into the aneurysmal bulge.

It should be also noted that the distribution of ink could be affected by flow changes caused by current patient conditions, such as psychological stress, exercise, anti-platelet drugs, etc. In order to take these into account, one may need to model the influence of the proximal and distal circulation. This can be incorporated into the model by various methods where a relatively simple model of the downstream circulation is linked to a three-dimensional, patient-specific model.4,26 While this is a very important topic, we feel this is beyond the scope of this study. We have shown in our previous work that the flow fields obtained in these models correspond to in vivo measurements conducted with PC-MRI velocimetry.18 In this study, we used these flow fields for an initial attempt to determine if there is a link between regions of increased RT and thrombus deposition.

The virtual ink technique provides a powerful tool for flow visualization and can be used for estimation of the flow RT distribution. While this method has been previously used in other CFD studies, e.g., in geophysical fluid dynamics, this is, to our knowledge, the first application of this concept to the flow in aneurismal arteries. The time difference between the arrival and washout of the virtual ink provides a “poor man’s” method for obtaining an Eulerian RT map. The main difficulty we found in applying this method is in the case when the virtual ink cycles in and out of the same area, e.g., a recirculating region where ink-free and ink-painted flow can replace each other several times.

In spite of these deficiencies of the method, it is a useful tool to visualize general flow patterns and to determine regions of increased RT for a given patient. The correlation determined between the CFD-predicted regions with increased RT and the regions of thrombus deposition observed in vivo with MR angiography demonstrates that along with the shear stress, flow RT can be a clinically relevant hemodynamic parameter. Thus, the flow RT and WSS distributions computed with CFD models can provide additional information to clinicians and help them in the prediction of aneurysm evolution for a given patient.


The results show that aneurysmal wall regions with increased flow RT and low WSSs may be prone to thrombus deposition. Patient-specific CFD models were constructed from MRI data for three basilar aneurysms that were thrombus-free and then proceeded to develop intra-aneurismal thrombus. A correlation was determined between the regions with increased RT computed with CFD simulations in baseline geometries and regions observed to clot at the follow-up MR studies. The relation between low WSS and thrombus location was also found significant. A model including both increased RT and low WSS was found to have a significantly higher thrombus prediction value than models considering either RT or WSS alone. The “virtual ink” technique, based on the advection of a passive scalar, provides flow visualization as well as an easy estimation of flow RT. The study demonstrates that hemodynamic factors can provide important additional information on aneurysm progression on a patient-specific basis.


We acknowledge grant support from the NIH (K25NS059891 (VR) and R01NS059944 (DS), and a VA Merit Review Grant (DS). Vitaliy Rayz would also like to thank David A. Steinman for insights into the issue of mesh resolution effects on the flow residence time calculations.

Open Access This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.


1. Boussel L, Rayz V, McCulloch C, Martin A, Acevedo-Bolton G, Lawton M, Higashida R, Smith WS, Young WL, Saloner D. Aneurysm growth occurs at region of low wall shear stress: patient-specific correlation of hemodynamics and growth in a longitudinal study. Stroke. 2008;39:2997–3002. doi: 10.1161/STROKEAHA.108.521617. [PMC free article] [PubMed] [Cross Ref]
2. Cebral JR, Castro MA, Appanaboyina S, Putman CM, Millan D, Frangi AF. Efficient pipeline for image-based patient-specific analysis of cerebral aneurysm hemodynamics: technique and sensitivity. IEEE Trans. Med. Imag. 2005;24:457–467. doi: 10.1109/TMI.2005.844159. [PubMed] [Cross Ref]
3. Einav S, Bluestein D. Dynamics of blood flow and platelet transport in pathological vessels. Ann. NY Acad. Sci. 2004;1015:351–366. doi: 10.1196/annals.1302.031. [PubMed] [Cross Ref]
4. Formaggia, L., J. Gerbeau, F. Nobile, and A. Quarteroni. On the coupling of 3D and 1D Navier–Stokes equations for flow problems in compliant vessels. In: Rapport de Recherche, INRIA, Editor, 2000.
5. Giddens DP, Zarins CK, Glagov S. The role of fluid mechanics in the localization and detection of atherosclerosis. J. Biomech. Eng. 1993;115:588–594. doi: 10.1115/1.2895545. [PubMed] [Cross Ref]
6. Glagov S, Zarins C, Giddens DP, Ku DN. Hemodynamics and atherosclerosis. Insights and perspectives gained from studies of human arteries. Arch. Pathol. Lab. Med. 1988;112:1018–1031. [PubMed]
7. Haimes, R. Using residence time for the extraction of recirculation regions. In: 14th Computational Fluid Dynamics Conference, Norfolk, VA, USA, 1999, pp. 345–354
8. Himburg HA, Grzybowski DM, Hazel AL, LaMack JA, Li XM, Friedman MH. Spatial comparison between wall shear stress measures and porcine arterial endothelial permeability. Am. J. Physiol. Heart. Circ. Physiol. 2004;286:H1916–H1922. doi: 10.1152/ajpheart.00897.2003. [PubMed] [Cross Ref]
9. Humphrey JD, Taylor CA. Intracranial and abdominal aortic aneurysms: similarities, differences, and need for a new class of computational models. Annu. Rev. Biomed. Eng. 2008;10:221–246. doi: 10.1146/annurev.bioeng.10.061807.160439. [PMC free article] [PubMed] [Cross Ref]
10. Jesty J, Yin W, Perrotta P, Bluestein D. Platelet activation in a circulating flow loop: combined effects of shear stress and exposure time. Platelets. 2003;14:143–149. doi: 10.1080/0953710031000092839. [PubMed] [Cross Ref]
11. Jou LD, Lee DH, Morsi H, Mawad ME. Wall shear stress on ruptured and unruptured intracranial aneurysms at the internal carotid artery. Am. J. Neuroradiol. 2008;29:1761–1767. doi: 10.3174/ajnr.A1180. [PubMed] [Cross Ref]
12. Jou LD, Wong G, Dispensa B, Lawton MT, Higashida RT, Young WL, Saloner D. Correlation between lumenal geometry changes and hemodynamics in fusiform intracranial aneurysms. Am. J. Neuroradiol. 2005;26:2357–2363. [PubMed]
13. Kroll MH, Hellums JD, McIntire LV, Schafer AI, Moake JL. Platelets and shear stress. Blood. 1996;88:1525–1541. [PubMed]
14. Kunov MJ, Steinman DA, Ethier CR. Particle volumetric residence time calculations in arterial geometries. J. Biomech. Eng. 1996;118:158–164. doi: 10.1115/1.2795954. [PubMed] [Cross Ref]
15. Longest PW, Kleinstreuer C. Numerical simulation of wall shear stress conditions and platelet localization in realistic end-to-side arterial anastomoses. J. Biomech. Eng. 2003;125:671–681. doi: 10.1115/1.1613298. [PubMed] [Cross Ref]
16. Malek AM, Izumo S, Alper SL. Modulation by pathophysiological stimuli of the shear stress-induced up-regulation of endothelial nitric oxide synthase expression in endothelial cells. Neurosurgery. 1999;45:334–344. doi: 10.1097/00006123-199908000-00028. [PubMed] [Cross Ref]
17. Ramstack JM, Zuckerman L, Mockros LF. Shear-induced activation of platelets. J. Biomech. 1979;12:113–125. doi: 10.1016/0021-9290(79)90150-7. [PubMed] [Cross Ref]
18. Rayz VL, Boussel L, Acevedo-Bolton G, Martin AJ, Young WL, Lawton MT, Higashida R, Saloner D. Numerical simulations of flow in cerebral aneurysms: comparison of CFD results and in vivo MRI measurements. J. Biomech. Eng. 2008;130:051011. doi: 10.1115/1.2970056. [PubMed] [Cross Ref]
19. Rayz VL, Boussel L, Lawton MT, Acevedo-Bolton G, Ge L, Young WL, Higashida RT, Saloner D. Numerical modeling of the flow in intracranial aneurysms: prediction of regions prone to thrombus formation. Ann. Biomed. Eng. 2008;36:1793–1804. doi: 10.1007/s10439-008-9561-5. [PMC free article] [PubMed] [Cross Ref]
20. Shadden SC, Lekien F, Marsden JE. Definition and properties of Lagrangian Coherent Structures from finite-time Lyapunov exponents in two-dimensional aperiodic flows. Physica D. 2005;212:271–304. doi: 10.1016/j.physd.2005.10.007. [Cross Ref]
21. Shadden SC, Taylor CA. Characterization of coherent structures in the cardiovascular system. Ann. Biomed. Eng. 2008;36:1152–1162. doi: 10.1007/s10439-008-9502-3. [PMC free article] [PubMed] [Cross Ref]
22. Shojima M, Oshima M, Takagi K, Torii R, Hayakawa M, Katada K, Morita A, Kirino T. Magnitude and role of wall shear stress on cerebral aneurysm: computational fluid dynamic study of 20 middle cerebral artery aneurysms. Stroke. 2004;35:2500–2505. doi: 10.1161/01.STR.0000144648.89172.0f. [PubMed] [Cross Ref]
23. Steinman DA. Image-based computational fluid dynamics modeling in realistic arterial geometries. Ann. Biomed. Eng. 2002;30:483–497. doi: 10.1114/1.1467679. [PubMed] [Cross Ref]
24. Steinman DA, Milner JS, Norley CJ, Lownie SP, Holdsworth DW. Image-based computational simulation of flow dynamics in a giant intracranial aneurysm. Am. J. Neuroadiol. 2003;24:559–566. [PubMed]
25. Tambasco M, Steinman DA. On assessing the quality of particle tracking through computational fluid dynamic models. J. Biomech. Eng. 2002;124:166–175. doi: 10.1115/1.1449489. [PubMed] [Cross Ref]
26. Vignon-Clementel IE, Figueroa CA, Jansen KE, Taylor CA. Outflow boundary conditions for three-dimensional finite element modeling of blood flow and pressure in arteries. Comput. Methods Appl. Mech. Eng. 2006;195:3776–3796. doi: 10.1016/j.cma.2005.04.014. [Cross Ref]
27. Wiebers DO, Whisnant JP, Huston J, III, Meissner I, Brown RD, Jr, Piepgras DG, Forbes GS, Thielen K, Nichols D, O’Fallon WM, Peacock J, Jaeger L, Kassell NF, Kongable-Beckman GL, Torner JC. Unruptured intracranial aneurysms: natural history, clinical outcome, and risks of surgical and endovascular treatment. Lancet. 2003;362:103–110. doi: 10.1016/S0140-6736(03)13860-3. [PubMed] [Cross Ref]

Articles from Springer Open Choice are provided here courtesy of Springer