Flow residence time in intracranial aneurysms evaluated by in vitro 4D flow MRI

The process of an intracranial aneurysm development, growth, and rupture is multifaceted and complex. In addition, clinical observations have identified the potential of thrombus formation within such aneurysms. While the underlying mechanism is not fully understood, the thrombi represent a potential risk factor for ischemic stroke. Emerging studies indicate that blood residence time (RT) is a promising hemodynamic metric associated with the aneurysm rupture and formation of intra-aneurysmal thrombi. Here, we present a methodology to experimentally evaluate both trajectory-wise and local RT based on magnetic resonance imaging (MRI) veloc-imetry, and apply it to in vitro flow measurements in scaled-up replicas of 9 patient-specific intracranial aneurysms. Lagrangian tracks of massless tracers are integrated from the velocity fields and averaged to return the mean RT in the aneurysm sac. This is found to be closely approximated by a simple time scale based on the sac diameter and space – time average of the aneurysmal fluid velocity. The mean RT is also correlated with the inflow time scale at the parent artery. These results also provide a basis for the estimation of RT when high-resolution hemodynamic maps are not available. With the continuous increase in accuracy and resolution enabled by progress in MRI technology, the methodology described here may in the future be applicable to in vivo data.


Introduction
Although the mechanisms of formation and progression of intracranial aneurysms (IAs) are not fully understood, it is believed that morphology, biochemical processes, hemodynamics, and tissue mechanics play important roles in their evolution (Lasheras, 2007;Rayz and Cohen-Gadol, 2020;Sforza et al., 2009).Historically, the risk of intracranial aneurysm rupture has been predominantly on the aneurysm size and factors such as hypertension, gender, or family history (Greving et al., 2014).While such factors continue to be used when determining the risk of rupture, certain statistical studies (e.g.Backes et al. (2014)) suggest that their use is only weakly correlated with the actual aneurysm rupture.As such, certain intracranial aneurysms undergo intervention when the actual risk of rupture is low, while others are observed when they are at high risk.This prompts development of biomarkers that can be implemented in clinic using viable technologies.
Wall shear stress (WSS), oscillatory shear index (OSI), and residence time (RT) are among hemodynamic metrics that have been proposed to investigate the association between hemodynamics and the aneurysm initiation, growth, and rupture.While the research has mainly focused on predicting risk of lesion rupture, there are clinical observations of thrombus formation at the aneurysmal wall as a secondary complication (Sugiyama et al., 2013;Wiebers et al., 2003).On one hand, thrombus formation inside the aneurysm can lead to stroke.On the other hand, thrombosis is a key healing process after endovascular treatments like coiling and flow diversion (Szikora et al., 2010).
RT as a metric for the aneurysm rupture that also has association with thrombosis formation is the time during which a particle/fluid element is present in a certain region.As recently reviewed by Reza and Arzani (2019), several definitions of RT have been used, following Lagrangian approaches and Eulerian ones; see e.g., (Esmaily-Moghadam et al., 2013;Himburg et al., 2004;Rayz et al., 2010;Reza and Arzani, 2019;Shadden and Arzani, 2015;Suh et al., 2011).Prolonged RT promotes aggregation and binding of activated platelets to adhesive proteins (Jesty et al., 2003;Shadden and Arzani, 2015).The blood flow recirculation can then hinder the transport of thrombotic factors and adhesive proteins out of the aneurysm, thus facilitating platelet activation and deposition of thrombotic substances to the wall.
Using computational fluid dynamics (CFD), Rayz et al. (2008) showed a similarity between the locations of intra-aneurysmal regions of low shear stresses and regions of thrombus deposition observed by magnetic resonance imaging (MRI).Other computational studies also suggest that intra-aneurysmal thrombosis may be related to prolongation of RT (Sugiyama et al., 2013).Rayz et al. (2010) studied RT by simulating passive scalar advection by CFD and showed a correlation between thrombotic region and prolonged RT.While powerful, numerical simulations are impacted by assumptions for the governing equations and boundary conditions, prompting the need for further experimental validations.In vitro 4D flow MRI or PIV measurements in flow phantoms provide valuable insight as well as benchmarking data (Amili et al., 2018;Brindise et al., 2019;Discetti and Coletti, 2018;Markl et al., 2012).
This work presents a methodology to experimentally evaluate RT based on 4D flow MRI data, and applies it to the flow measured in scaled-up replicas of 9 patient-specific intracranial aneurysms.One of the main objectives is to investigate the correlation between RT (for which high-resolution data is required and involved post-processing) and other more immediate observables, such as aneurysm size and mean blood velocity in the aneurysm sac.Determining such correlation as presented in this study has important clinical implications.First, RT is an important hemodynamic descriptor that is difficult to accurately measure in vivo with present-day technology.Therefore, being able to estimate it from readily available quantities would enable its use as a biomarker in large-scale studies.Second, RT may inform our predictive understanding of the evolution of intracranial aneurysms.Therefore, deciphering its relationship with commonly used observables may pave the way to explain the contradictory data in the literature.

In vitro models and flow velocity measurements
4D flow MRI was used to obtain velocity data in 3D-printed replicas of 9 intracranial aneurysms presented in Fig. 1(a).In vivo aneurysm geometries were obtained from patients recruited to this study using the approved institutional review board (IRB) consent via computed tomography (CT) or MRI Time of Flight (TOF) scans.The methodology to reconstruct the vascular anatomy and realize the in vitro experiments was described in detail in Amili et al. (2018) and is only briefly  summarized here.Table 1 lists the aneurysm location and size.The diameter D sac was calculated as the cubic root of the sac volume, and ranges between 2 mm and 6 mm.Patient-specific models were manufactured by 3D printing at the W. M. Keck Center, University of Texas at El Paso.High-resolution stereo-lithography guarantees hydrodynamically smooth walls (Amili et al., 2019;Amili et al., 2020;Amili et al., 2018).The native geometries were scaled up by a factor of 3 to increase the effective flow measurement resolution.The models were placed into a flow circuit driven by a computercontrolled piston pump, generating a waveform such that the Reynolds number and Womersley number are matched with physiological values.To this end, the kinematic viscosity (ν) of the working fluid (water-glycerol mixture) is set to 50% lower than blood; the velocity is 6 times slower, and the period (T) of the waveform is 18 times longer than in vivo, according to the Reynolds number and Womersley number definitions: Here U cycle and D are the cycle-averaged bulk flow velocity and diameter of the parent artery at the inlet of the 3D-printed model, respectively.As localized in vivo flow information was not available for all cases, we applied to all cases a waveform (Fig. 1(b)) obtained by averaging the transcranial Doppler ultrasound from two of the patients at two sites each: right internal carotid artery and right middle cerebral artery.
Although the cerebral blood flow can vary between different sites (Zarrinkoob et al., 2015), for a given geometry, flow rate and waveform are not expected to strongly affect the hemodynamics (Cebral et al., 2005).While the flow in aneurysms with more than one inlet or outlet is likely to be affected by the flow ratio between branches in real applications, it is a common practice to divide the flow rate proportional to the cross-sectional area at the branching sites in the absence of such data.Re and Wo values are listed in Table 1.
The in vitro 4D flow MRI measurements were carried out on a Siemens MAGNETOM Prisma 3 T MRI scanner at the University of Minnesota Center for Magnetic Resonance Research (CMRR).The 4D flow MRI acquisition was based on the sequence described in Markl et al. (2012) with a GRAPPA acceleration factor of 2, reconstructing 3D velocity fields with a voxel size of 0.6 mm in all directions over 16 phases of the cardiac cycle.CuSO 4 with a concentration of approximately 0.06 mol/L was added to the working fluid in order to enhance the MRI signal-to-noise ratio (SNR) following Elkins and Alley (2007).The basic parameters of the MRI measurements are listed in Table 2.The VENC parameter was comparable to physiological velocity to avoid aliasing, and the scan was repeated to enhance SNR.Following Pelc et al. (1994) and Haacke et al. (1991), the SNR is calculated as: where 〈S〉 is the mean of MR signal inside the wetted volume, 〈B〉 is the mean of the MR signal in the background, and σ b is the standard deviation of MR in the background.We notice that the spatial variability of the SNR is moderate (typically within ±10% of the mean).Therefore, even the low-velocity regions are deemed to yield sufficiently accurate results to evaluate Lagrangian trajectories.Using the estimated SNR level, the corresponding velocity uncertainty, depending on the case, is in the range of 2-10% with respect to the cycle-average velocity at the parent vessel.For further details of the experimental methodology and MRI pre-and post-processing steps, see Amili et al. (2018).

Lagrangian particle tracking
The tracking methodology from MRI velocimetry was presented and validated in Amili et al. (2019) and Amili et al. (2020).Based on the 4D flow MRI data scaled back to the physiological condition (blood flow in the native size artery), Lagrangian particle tracking is implemented by solving the particle equation of motion in a general form of Maxey and Riley (1983): where ρ p is the particle density, ρ f is the fluid density, u p is the particle velocity, u is the fluid velocity interpolated at the particle location, g is the gravitational acceleration, μ is the fluid dynamic viscosity, and Dp is the particle diameter.The right-hand side includes the net force from gravity and buoyancy; the undisturbed flow stress; the added mass term; and the drag force.The solver is general and thus includes forces acting on inertial particles which may represent, e.g., blood clots (Amili et al., 2020).For the purpose of this study, we consider small neutrally buoyant particles that act as fluid tracers.The problem reduces to the kinematic tracking of fluid elements along pathlines.In the following, the term particle and tracer will be used interchangeably.
To avoid potential bias due to different release timings and locations, particles are released throughout the cardiac cycle over the 16 temporal bins and volumetrically on a Cartesian grid.The latter is congruent to the collocation points of the MRI data and refined by a factor of 2 in all directions, to obtain a sufficient number of trajectories.The number of released particles depends on the individual aneurysm geometry and volume, ranging from 80,000 to 300,000.Although volumetric seeding does not necessarily reflect the physiological concentration of blood components (Tambasco and Steinman, 2002), it can avoid the potential bias of planar seeding for the purpose of residence time calculation.A similar strategy has been reported by Reza and Arzani (2019), Shadden and Arzani (2015), and Suh et al. (2011).By solving the particle equation of motion, velocity and position vectors of every tracer are computed at every time step (~0.02 ms) with the 4th-order Runge-Kutta method.The computation is terminated when the majority of tracers are washed out of the flow domain.To complete Lagrangian tracking of all phases per aneurysm case, 16 cumulative cores were used at the Minnesota Supercomputing Institute for up to 12 h.

Lagrangian residence time
The Lagrangian RT is computed as the time a tracer spends along its trajectory within the sac before being washed out: where t in and t out indicate the time the tracer enters and leaves the sac, respectively.The aneurysm sac is identified as the bulged or dilated volume, past a neck plane that separates it from the parent artery (Fig. 1).The time steps between the tracer entrance and exit through the "neck plane" add up to the RT associated to each trajectory.By inspecting the trajectories (Fig. 2(a)), a large fraction of them appear to barely enter the aneurysm.Because we are interested in fluid elements representative of the hemodynamics within the aneurysm, we define a threshold plane parallel to the neck plane at a distance d from the sac apex (Fig. 2(b)).Fig. 2(c) shows the effect of varying d/D sac on the probability density function (PDF) of the Lagrangian RT for one aneurysm case.Other definitions of the sac size (e.g. the maximum dimension) lead to equivalent conclusions.As d/D sac decreases (i.e. the threshold plane moves towards the apex), the trajectories that marginally penetrate the sac are excluded.The PDF thus shifts to higher RT values, and the mode of the distribution becomes robust around d/D sac = 0.7.This threshold is applied to all cases; the results are robust to a ±10% variation of its value.

Local residence time
The recirculation regions, where the blood tends to reside for extended periods, are identified by evaluating the local residence time RT local (also known as the mean exposure time), calculated in each voxel of the aneurysm volume using the Lagrangian trajectories.For robust particle sampling, the voxels of the original data are volumetrically binned with a kernel size of 3 in all directions; i.e., RT local is calculated over a grid three times coarser than the native voxel resolution.Lagrangian trajectories of 5 evenly spaced bins of the 16 phases throughout the cardiac cycle are used which is verified to be sufficient to provide statistical convergence.The particle accumulation time in each spatial bin (of volume V bin ) is calculated by summing up the number of time steps each tracer is present in the bin over the cardiac cycle.This is then normalized dividing by the number of particles passing through each bin (N), its side length (Shadden and Arzani, 2015), and number of released phases: Integrating the local residence time [s/m] along a tracer trajectory within the aneurysm sac, returns within numerical error the trajectorywise RT [s] defined in the last section.

Mean residence time (volumetric time scale)
As Lagrangian tracking of particles is computationally expensive, a simple volumetric time scale based on the cycle-averaged velocity field is considered.This Eulerian metric represents an order-of-magnitude estimate of the eddy-turnover-time inside the aneurysm: where U sac is the spatio-temporal average of the fluid velocity within the sac at all measured points.Compared to the other measures, T vol lifts stringent requirements on the spatial and temporal resolution of 4D flow MRI, as it can be obtained from non-phase-resolved measurements with relatively coarse resolution achievable in many in vivo settings.

Results
PDFs of Lagrangian RT are presented in Fig. 3, displaying a broad variability of distributions that reflects the differences in the morphology and flow patterns between the nine cases.Here and in the following, time scales are reported after re-scaling from the laboratory to the physiological conditions.The RT amounts to less than 10% of the cardiac cycle duration for all cases.Fig. 4 presents iso-surfaces of RT local and the normalized velocity magnitude U mag /U cycle .The former is normalized by the characteristic length πD sac (affine to the sac perimeter) and RT mean (the mean Lagrangian RT averaged over all trajectories for each case); while the latter is normalized by the time average of the bulk velocity at the parent vessel, U cycle .RT local maps reveal the complex topology of the quasi-stagnation areas in the vasculature.As expected, the regions where fluid resides for extended periods (high RT local , in shades of green) are complementary to those where the velocity is high (in shades of purple).However, low velocity regions (not shown) are topologically different from the high-RT local regions, and largely overlap the near-wall envelope (where the velocity vanishes due to the no-slip boundary condition).In other words, the local residence time cannot be deduced immediately from the velocity fields.
The mean of the Lagrangian RT is plotted in Fig. 5 against the volumetric time scale T vol for all cases.A strong positive correlation exists between both metrics, and in fact T vol appears to be a robust estimator of RT.Without case 9 (the only case that qualifies as a large aneurysm), the correlation coefficients are larger than 0.9.Using the mode and median of RT instead of the mean leads to analogous results.The mean of the local residence time RT local (averaged within the sac volume and multiplied by πD sac to convert it to temporal units) is plotted versus the mean Lagrangian RT and versus T vol in Fig. 6(a) and 6(b), respectively.Overall, a strong correlation (with a p-value less than 0.05) and close quantitative agreement between the three metrics are found.Although we have a limited number of aneurysm cases, the correlation coefficient between the above parameters is statistically significant with the largest p-value of 0.03 that is below the threshold level of 0.05 for 95% confidence bounds.We highlight that the p-value calculation for the significance of the correlation coefficient is based on the"student ttest" that assumes a normal distribution for the sample population (Walpole et al., 2012).The sample size for the point-wise data as shown in Figs. 5 and 6 is small and its normality cannot be confirmed.However, as the sample size becomes larger, one can expect a normal sampling distribution as observed in the PDFs of residence time of the Lagrangian tracks in Fig. 3.
The global Eulerian time scale T vol is expected to be determined in part by the time scale of the flow in the parent artery.Thus, the latter is expected to be influential for and correlate with the residence time in the aneurysm.This is confirmed in Fig. 7, where the inflow time scale D in / U in is shown to be strongly correlated with the mean Lagrangian RT (as well as with mean RT local and T vol , not shown).Here U in and D in are the mean flow velocity and diameter of the parent artery, both measured approximately two vessel diameters upstream of the aneurysm location.
On the other hand, all the considered metrics of RT are not found to be significantly correlated with the sac diameter for the small aneurysms; for example, the mean Lagrangian RT and the diameter D sac exhibit a correlation coefficient r = 0.2.A similar trend is observed for the correlation between the mean Lagrangian RT and the aneurysm largest dimension.

Discussion
Based on in vitro 4D flow MRI carried out in 9 patient-specific models of intracranial aneurysms, virtual particle tracking is implemented, and three metrics of residence time are compared: Lagrangian RT, local RT, and volumetric time scale.The first metric is based on trajectories obtained by integrating high-resolution 4D velocity fields.The second metric leverages the same trajectories to generate a spatial distribution of RT, providing detailed topological maps.Appropriate averaging of these two metrics returns mean RT values, quantifying the time spent by fluid elements inside the aneurysmal sac.The third metric is a global time scale based on the size and the mean fluid velocity inside the sac.The correlation between these metrics is investigated and found to be strong.In particular, the correlation with T vol implies that RT can be robustly estimated from in vivo velocity measurements at low spatiotemporal resolutions using magnetic or ultrasonic imaging methods.This is significant because, despite remarkable progress in ultra-highfield MRI (Ertürk et al., 2017), the accuracy and resolution attained in phantom studies remains presently out of reach in in vivo settings.Moreover, evaluating T vol does not require the computationally expensive integration of the velocity field.
With respect to many existing studies predominantly focused on large/giant aneurysms, we have considered relatively small aneurysms (≤7.0 mm in maximum dimension).The choice was driven by the clinical interest in cases for which decisions concerning treatment are complex, while data on and knowledge of the hemodynamics are scarce.While the flow structure in larger aneurysms is likely to be more complex, with regions of low and high RT, our approach will still be applicable.In fact, the accuracy is likely to improve due to alleviated constraint in relative spatial resolution.
Given the relation with thrombus formation, the clinical relevance of RT as a predictive biomarker in the cycle of aneurysm progression is evident.Although a direct, mechanistic connection between RT and rupture risk is not known, a strong inverse correlation exists between RT and WSS, which is known to have a profound influence on the aneurysm growth and fate (Xiang et al., 2011).In the present cases (small aneurysms), however, there is no measurable correlation between RT and sac diameter.In this context, it shall be remarked that RT is an inherently multi-factorial observable dependent on the vascular anatomy and morphology, flow direction and velocity, among other parameters.Still, the present study shows RT can be robustly predicted by combining one appropriate length scale and one velocity scale.This lowers the demands for high resolution and is especially advantageous for small aneurysms.Moreover, a significant correlation also exists between RT and the inflow time scale D in /U in , which can be estimated a priori from detailed knowledge of the brain circulation.This suggests the possibility of predicting trends of RT in aneurysms at given locations, independently of the ability to image the lesion and measure its blood flow velocity, which in turn could offer explanations for the variation in rupture rates found between aneurysms of identical size at different locations in cerebrovascular circulation.
The present study provides novel information on the Lagrangian transport in intracranial aneurysms.For small aneurysms, RT in the sac is found to be much shorter than the cardiac cycle.Rayz et al. (2010) reported RT durations of many cardiac cycles, but it is important to remark multiple definitions for RT exist, emphasizing different aspects of the transport process (Reza and Arzani, 2019).The relatively small dimension of the considered aneurysms is also likely to play a major role in dictating the magnitude of RT.In the limit of small sac diameter, RT reduces to the time the fluid takes to travel across a normal vessel; that is, as D tends to D in , RT tends to D in /U in which can be regarded as a lower bound for RT.Indeed, Fig. 7 shows that the Lagrangian RT correlates with and is several times larger than D in /U in .
It is important to remark that, while mapping a detailed local RT is possible using scaled-up aneurysm models in vitro (e.g., Fig. 4), it may not be achievable in vivo yet.However, this may be well envisioned in the near future due to ever-increasing ultra-high-field MRI capabilities (e.g., the first human head scan at 10.5 T by Sadeghi-Tarakameh et al. ( 2020)) and the recent approval of U.S. Food and Drug Administration for the clinical use of 7 T MRI.Until the fruits of these advances are ripened, the correlation of local/Lagrangian RT with global/Eulerian RT offers attractive alternatives, given the wide availability of MRI and/or transcranial Doppler (TCD) ultrasound.
A limitation of the present study is that it does not directly address whether the considered aneurysms were eventually thrombotic, and whether thrombi were generated in regions with high local RT.We remark that recent reviews (e.g., Ngoepe et al. (2018)) reported the thrombosis formation in IAs and pointed to a link between residence time and thrombosis.In addition, several case studies such as Brownlee et al. (1995) and Gonçalves et al. (2014) reported spontaneous thrombosis of large unruptured aneurysms.In such cases, a long residence time in the aneurysm sac can be naturally conjectured.The number of cases is  also limited, and the correlation among the quantities need to be confirmed by more extensive studies.In addition, the material used for the in vitro experiments have idealized properties: the phantom walls are rigid, and the working fluid is Newtonian.The challenge of reproducing tissue and blood properties, along with the expected influence of such simplifications have dictated analogous choices in previous studies (Rayz and Cohen-Gadol, 2020).Finally, although the scaled-up models allow for the effective spatial resolution to be finer than typical in vivo flow measurements, this may not be sufficient to resolve all the details of the near-wall fluid dynamics, which can be consequential for RT.The present smaller aneurysms also pose a greater challenge in terms of the relative resolution.In general, the spatio-temporal resolution is not sufficient to rule out the possibility that instantaneous Lagrangian trajectories be significantly different from those based on the timeaveraged flow fields reported by Valen-Sendstad et al. (2011).However, it is likely to reduce such discrepancies by averaging the instantaneous trajectories.
Importantly, with the continuous increase in accuracy and resolution warranted by progress in MRI technology (Ertürk et al., 2017;Sadeghi-Tarakameh et al., 2020;Schnell et al., 2017), the presented methodology may eventually be applied to in vivo data.Moreover, the generality of the method allows to immediately incorporate finite inertia associated to density and/or size of the particles, making it suitable to investigate the fate and transport of blood cells and clots.This will be the focus of future studies.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.(a) The nine considered aneurysms with sacs highlighted and flow directions indicated by the arrows.(b) The inflow waveform used for all cases.

Fig. 2 .
Fig. 2. (a) Example tracer trajectories colored by Lagrangian RT, (b) Definition of the threshold plane, (c) PDFs of Lagrangian RT at different d/D sac values for case 5.

Fig. 4 .
Fig. 4. Iso-surfaces of the scaled local residence time (green) and normalized velocity magnitude (magenta) with flow directions indicated by the arrows.The regions with velocity higher than the threshold (indicated by the color) are spatially complementary to those with prolonged local RT.

Fig. 6 .
Fig. 6.Comparison of mean local RT versus (a) mean Lagrangian RT and (b) volumetric time scale.

Fig. 7 .
Fig. 7. Comparison of the inflow time scale versus mean Lagrangian RT.

Table 1
Location, largest dimension, diameter, Reynolds and Womersley number for the considered aneurysms.

Table 2
Parameters of the in vitro 4D-flow MRI measurements.The same velocity encoding is applied along all directions.