Validation of Patient-Specific Cerebral Blood Flow Simulation Using Transcranial Doppler Measurements

We present a validation study comparing results from a patient-specific lattice-Boltzmann simulation to transcranial Doppler (TCD) velocity measurements in four different planes of the middle cerebral artery (MCA). As part of the study, we compared simulations using a Newtonian and a Carreau-Yasuda rheology model. We also investigated the viability of using downscaled velocities to reduce the required resolution. Simulations with unscaled velocities predict the maximum flow velocity with an error of less than 9%, independent of the rheology model chosen. The accuracy of the simulation predictions worsens considerably when simulations are run at reduced velocity, as is for example the case when inflow velocities from healthy individuals are used on a vascular model of a stroke patient. Our results demonstrate the importance of using directly measured and patient-specific inflow velocities when simulating blood flow in MCAs. We conclude that localized TCD measurements together with predictive simulations can be used to obtain flow estimates with high fidelity over a larger region, and reduce the need for more invasive flow measurement procedures.

1. Introduction

Computational fluid dynamics (CFD) has been widely applied by researchers to model blood flow in cerebral arteries and specifically within aneurysms (Cebral et al., 2011; Miura et al., 2013; Mountrakis et al., 2013; Byrne et al., 2014; Ouared et al., 2016). There is considerable interest in exploring the correlation between the dynamical properties of blood flow and clinical outcomes, with the long-term aim to provide a personalized, predictive simulation approach for aneurysm formation, growth, and/or rupture (Jou et al., 2008; Bernabeu et al., 2013; Xiang et al., 2014). When performing such simulations it is essential that computational models are able to deliver a realistic prediction of patient-specific flow velocities.

A range of simulation studies have been performed using patient-specific flow measurements derived from phase contrast magnetic resonance angiography (pc-MRA, see e.g., Boussel et al., 2008). However, Marzo et al. (2011) found that using this type of measurement provides limited accuracy benefits in comparison with modeled boundary conditions. The use of CFD in combination with transcranial Doppler (TCD) velocity measurements has been less extensively researched (see e.g., Hassan et al., 2004), primarily because reliable TCD measurements can only be made in a limited subset of the cerebral arteries. In addition, TCD measurements with hand-held devices may contain errors if held at an incorrect angle (e.g., an underprediction of approximately 1.6% if the angle is off by 10 degrees). However, the excellent time resolution of TCD allows for a more reliable detection of peak velocities, and helps to establish more precise pulsatile flow profiles. Indeed, the maximum velocity values found by TCD are frequently around 30% higher than those found through pc-MRA (Chang et al., 2011; Meckel et al., 2013). In addition, TCD is non-invasive, unlike pc-MRA, and both are widely applied in clinical settings today.

Blood consists of blood cells which reside within a liquid medium known as blood plasma. Blood has a viscosity which decreases under shear flow (shear-thinning), unlike water which exhibits a constant Newtonian viscosity regardless of shear strain rate. Many well-known CFD studies of cerebrovascular blood flow are performed using a Newtonian fluid model (e.g., Cebral et al., 2011; Miura et al., 2013; Byrne et al., 2014), though recent research has found that such an assumption could lead to over-estimation of wall shear stresses (WSS) in cerebral arteries and aneurysms (Bernsdorf and Wang, 2009; Xiang et al., 2011; Khan et al., 2016). As a result, it can also alter the outcome of related diagnostic techniques such as three-band diagram analysis (Bernabeu et al., 2013), a technique proposed by Chatzizisis et al. (2008) to compare WSS at a specific location, over a period of time, to a set of pathological threshold values.

Existing CFD studies of cerebrovascular flow frequently derive inflow velocities not from the specific patient of interest, but from healthy subjects (e.g., Miura et al., 2013; Byrne et al., 2014) or idealized pulsatile profiles (Womersley flow, e.g., Castro et al., 2006; Alnæs et al., 2007; Cebral et al., 2011). However, blood flow velocities in middle cerebral arteries (MCA) from healthy subjects are typically much lower than those from stroke patients or patients suffering from hypertension. In this context Venugopal et al. (2007) found that mean WSS properties of simulations at Reynolds numbers (Re) below 200 do not correspond in any linear way to WSS properties of simulations at Re = 340–675. Itani et al. (2015) investigated how the mean, maximum, and minimum wall shear stress changes when a patient is subject to different levels of exercise intensity. They also found a non-linear relation between maximum inflow velocity and extracted WSS.

In this work, we simulate blood flow in a patient-specific MCA model using patient-specific TCD measurements as inflow boundary conditions, and compare our predictions against clinical measurements at four locations. Our simulations employ the lattice-Boltzmann method at high resolution, a technique which has been shown by Jain et al. among others, to be particularly well-suited for simulating cerebrovascular and aneurysmal blood flow Jain et al. (2016). We perform simulations imposing the measured velocity from the individual patients at the inlet, and investigate how the choice of rheology model affects the predicted flow velocities throughout the MCA. In addition, we report on the accuracy of velocity predictions when running simulations with downscaled inlet velocities, and rescaling the velocities obtained from the measurement planes.

2. Materials and Methods

To perform our simulations, we use the HemeLB software (Groen et al., 2013; Nash et al., 2014) for lattice Boltzmann simulations of blood flow in cerebral arteries. The lattice Boltzmann method (LBM) is a highly scalable simulation approach which uses a discretized kinetic model on a regular lattice to reproduce the dynamics of incompressible fluid flow. The LBM can be interpreted as a numerical solver for the Navier-Stokes equation with the advantage that it algorithmically separates the non-linearity from the non-locality. Specific boundary conditions are applied to create accurate representations of fluid flow near vessel walls, as well as inflow and outflow boundaries. In our case, we adopt a 3-dimensional LBM which propagates fluid flow in 19 directions per grid point (D3Q19) using a BGK collision operator (see e.g., Succi, 2001 for details). For the boundary conditions, we used the Bouzidi (Bouzidi et al., 2001) model to represent flow interactions with the vessel walls. Patient-specific inflow conditions were obtained from TCD measurements performed at the National Hospital for Neurology and Neurosurgery (NHNN) using the Doppler BoxX (with a handheld device) from the DWL company, and used rotational angiography data from NHNN to obtain imaging data from the same patient. TCD measurements were recorded for at least six cardiac cycles each in the right MCA, consecutively at depths of 49, 54, 57, 59, and 63 mm away from the temple area (see Figure 1 for the location of the TCD validation planes in the 3D model, Table 1 for the velocity measurements, and Figure 2 for the TCD image measurement at the inflow boundary). The Doppler BoxX provides a flow direction indication at all depths whenever a measurement is made. In our case, this feature enabled us to hold the TCD device such that the flow was observable in the right MCA, as well as the right Anterior Communicating Artery (ACA). This is important, because retaining such a tight orientation minimizes TCD measurement errors caused by holding the device at a wrong angle. In addition, to align the TCD measurements precisely with the corresponding planes of flow direction in the simulation domain, we performed a triangulation and an angle correction with respect to the perpendicular flow direction (see Table 2 for our triangulation results). The maximum velocity at the inflow boundary, extracted from the TCD data, was 1.50 m/s.

FIGURE 1

Figure 1. Overview of the overall patient vasculature in the medical images (A), and the patient-specific 3D model used in our simulations (B). As part of the simulation model overview, we indicate the four TCD measurement planes used for validation. Both the inlet and the 63 mm TCD measurement plane are at the right hand side of the image.

TABLE 1

Table 1. Overview of measured and simulated flow characteristics in the MCA, as well as relative differences between measurement and simulation.

FIGURE 2

Figure 2. Raw TCD input image of the measured velocity at a depth of 63 mm (inflow boundary plane). The measured velocity at the selected depth (63 mm) is given at the top, while the general flow direction at all depths is given at the bottom, either toward the device (red) or away from it (blue). We observe a change in flow direction around a depth of 67 mm, which is at the junction between the right MCA and the right ACA.

TABLE 2

Table 2. Triangulation points, input, output, and measurement plane locations (and orientations where applicable) in the simulation domain, used to calculate the angle correction.

Extracted cardiac cycle lengths vary for each cardiac cycle and each measurement. The patient is known to have an existing aneurysm in the MCA on the opposite (left) side, within which the velocity magnitudes could not be clearly resolved using TCD due to its unfavourable orientation. We segmented the images using VMTKlab (vmtklab.orobix.com), and voxelized the 3D model using the HemeLB setup tool. The resulting geometry has one inflow region and five outflow regions—two small ones at the top near the inflow boundary, two larger ones at the bottom, and the largest one left of the 49 mm plane (see Figure 1B).

The 2D inflow profiles were reconstructed from the 1D TCD data by mapping a parabolic profile to the non-circular inlets. This parabolic inlet profile has the original velocity from the 1D TCD data mapped to the centre of the inlet (the lattice site which is furthest from any wall), and 0 velocity values mapped to wall boundary sites. The velocity magnitude of a given lattice site is then calculated using a parabolic equation, which depends on the distance of the lattice site to the nearest vessel wall site in the inlet plane (0 for wall sites, 1 for the site in the centre, and values in between for other sites).

The boundary conditions in the lattice Boltzmann method were implemented as follows. To set the reconstructed velocity profile ${\overrightarrow{u}}_{\mathrm{\text{TCD}}}({\overrightarrow{x}}_{\mathrm{\text{in}}},t)$ at the inlet, we use a method introduced by Ladd (1994). A simple bounce-back boundary condition is augmented with a momentum term that results in a time-dependent Dirichlet condition for the velocity

$\begin{array}{rc}\overrightarrow{u}({\overrightarrow{x}}_{\text{in}},t)={\overrightarrow{u}}_{\text{TCD}}({\overrightarrow{x}}_{\text{in}},t).& (1)\end{array}$

At the outlet, we employed an open boundary condition in terms of a mixed Dirichlet-Neumann boundary condition (Nash et al., 2014)

$\begin{array}{rc}{\overrightarrow{u}}_p({\overrightarrow{x}}_{\text{out}},t)=0,& (2)\end{array}$

$\begin{array}{rc}\widehat{n}·\nabla {\overrightarrow{u}}_n({\overrightarrow{x}}_{\text{out}},t)=0,& (3)\end{array}$

where $\widehat{n}$ is the normal vector of the outlet plane, and ${\overrightarrow{u}}_p$ and ${\overrightarrow{u}}_n$ are the in-plane and normal components of the outlet velocity, respectively. The gradient in Equation (3) is taken as the first-order finite difference approximation on the lattice Boltzmann grid. In the implementation by Nash et al. (2014), the density $\rho ({\overrightarrow{x}}_{\mathrm{\text{out}}},t)$=ρ₀ at the outlet is prescribed in order to determine the unknown fluid variables. It is worth noting that prescribing the density at the outlet fixes the static pressure through the ideal gas equation of state. However, this does not constrain the dynamic pressure which varies over a cardiac cycle as shown in Figure 3.

FIGURE 3

Figure 3. Differential pressure at the main outlet plane, relative to the ideal gas pressure for average density in the simulation. The maximal pressure found in the plane is given by the red dashed line, while the minimal pressure found in the plane is given by the blue dotted line. The average pressure in the plane is given by the black line.

The shear-thinning behavior of blood is modeled using the Carreau-Yasuda (CY) model which employs the expression (Boyd et al., 2007; Bernabeu et al., 2013)

$\begin{array}{rc}\frac{\eta (\dot{\gamma })-{\eta }_{\infty }}{{\eta }_0-{\eta }_{\infty }}={\left(1+{(\lambda \dot{\gamma })}^a\right)}^{\frac{n-1}{a}}& (4)\end{array}$

to account for the dependence of the dynamic viscosity η on the shear rate $\stackrel{°}{\gamma }$. Here, η₀ and η_∞ are the asymptotic values at zero and infinite shear rate, and a, n, λ are empirical materials parameters that describe the shear-thinning curve. The CY model represents a smooth transition between Newtonian behavior at η₀ and η_∞.

The HemeLB simulations were performed on the ARCHER supercomputer at EPCC in Edinburgh, United Kingdom, and the SuperMUC supercomputer at LRZ in Garching, Germany. We used between 1,536 and 24,768 cores, depending on the chosen resolution.

2.1. Choice of Lattice Boltzmann Parameters

Our lattice Boltzmann model uses a D3Q19 lattice with the Lattice Bhatnagar-Gross-Krook (LBGK) collision model (Bhatnagar et al., 1954). The relaxation parameters are set to yield the dynamic viscosity of blood of η = 0.004 Pa·s. The parameters used in the CY model are η₀ = 0.16 Pa·s, η_∞ = 0.0035 Pa·s, λ = 8.2 s, a = 0.64 and n = 0.2128 as given by Boyd et al. (2007) and previously adopted by Bernabeu et al. (2013). In our full-resolution, full-velocity simulations, we used a voxel size of 10 μm, a time step size of 0.28 μs, and a geometry consisting of 174,738,326 fluid sites. The simulations ran for 21.43 million time steps, which corresponds to 5 s of simulated time following a one-second “warmup” period (during which the inlet flow speed is increased gradually from rest in order to avoid flow instability or shockwaves resulting from a step change). The Reynolds number of our full-velocity simulation is approximately 966, based on a characteristic diameter of 24 mm with the highest measured peak velocity of 1.61 m/s.

We also performed simulations at reduced velocity and resolution, multiplying the velocities by 50 or 25%, as well as with increased voxel sizes of 20 and 40 μm. We discuss the implications of using this type of velocity scaling in detail in the next subsection.

2.2. Velocity Scaling

The LBM is valid in the incompressible regime and introduces compressibility errors that scale quadratically in the Mach number Ma = U/c_s, where U is the flow velocity and c_s the speed of sound. The cardiac flow is characterized by the Reynolds number Re = UD/ν and the Womersley number α = (ω^D2/ν)1/2, where D is the vessel diameter, ν = η/ρ is the kinematic viscosity, and ω is the angular frequency of the oscillating flow due to the cardiac cycle. In terms of the simulation parameters, the kinematic viscosity of the lattice BGK model and the speed of sound are given by

$\begin{array}{rc}\nu =\frac{1}{3}\left(\widehat{\tau }-\frac{1}{2}\right)\frac{{(\Delta x)}^2}{\Delta t},& (5)\end{array}$

$\begin{array}{rc}{c}_s=\frac{1}{\sqrt{3}}\frac{\Delta x}{\Delta t},& (6)\end{array}$

where $\widehat{\tau }$ is the dimensionless relaxation parameter of the BGK model, and Δx and Δt are the discrete lattice spacing and time step, respectively. Based on the Reynolds and Mach numbers, we have the following relation for the dimensionless relaxation parameter

$\begin{array}{rc}\widehat{\tau }-\frac{1}{2}=\sqrt{3}\frac{D}{\Delta x}\frac{Ma}{Re}.& (7)\end{array}$

Linear stability requires $\widehat{\tau }>0.5$ which guarantees a positive viscosity. However, it is mandatory to keep the Mach number small in order to reduce compressibility errors and make the system less prone to instabilities due to density fluctuations. In the standard diffusive scaling, convergence is achieved by reducing the Mach number while keeping the Reynolds number constant. This implies (Δx)²~Δt. Thus, reducing the Mach number by means of increasing resolution results in an increase in computational costs due the cubic scaling of volume.

Therefore, some authors have been tempted to use lower flow velocities, e.g., from healthy subjects (Miura et al., 2013; Byrne et al., 2014), in order to maintain stable simulations at a larger voxel size Δx. The ratio of the reduced velocity U′ and the original velocity U is denoted by a scaling factor s. This leads to a scaling relation

$\begin{array}{rc}s=\frac{U^{\prime }}{U}=\frac{{\nu }^{\prime }R{e}^{\prime }D}{\nu Re{D}^{\prime }}=\frac{{\alpha }^2{\omega }^{\prime }{D}^{\prime }R{e}^{\prime }}{{{\alpha }^{\prime }}^2\omega DRe},& (8)\end{array}$

where the prime denotes the quantities associated with the scaled velocity U′. If one insists on a fixed system size D′ = D and cardiac cycle length ω′ = ω, it is not possible to fix both the Womersley number and the Reynolds number at the same time such that the simulation is performed in a flow regime different to that of the full velocity simulation. In section 3.2, we demonstrate that this can significantly impact the simulated flow patterns.

3. Results and Discussion

We present results from three types of simulation. First, we compare our full velocity and full resolution (10 micron voxel size) simulations against the TCD measurements on the same patient. Second, we present the results from simulations at reduced velocity and reduced resolution, and compare them both with results from our full-scale simulations and with the TCD measurements. Third, we compare the results of simulations using a Newtonian rheology model to simulations using a non-Newtonian (Carreau-Yasuda) rheology model (Abraham et al., 2005).

3.1. Validating Full Velocity Haemodynamics Predictions Against Measurements

In Table 1 we present the maximum velocity v_max as measured with TCD and the simulation results for all four measurement planes. Our simulations predict the flow velocity with a relative error of less than 9% in all cases. In Figure 4 we present a direct comparison of our TCD measurements in the four planes over time, and our velocity predictions derived using HemeLB at the same locations. We observe good agreement between the simulation results and the measured TCD profile. The differences can be ascribed to the limitations of our approach (see section 3.4) and uncertainties in the measurements, including phase misalignments due to the sequential nature of the TCD measurement.

FIGURE 4

Figure 4. Comparison of the maximum velocity using both TCD (blue line) HemeLB (green line), given for the planes at 49 mm (Top left), 54 mm (Top right), 57 mm (Bottom left), and 59 mm (Bottom right). HemeLB results presented here are for the run at 100% velocity and with Newtonian rheology. The phase has been shifted to align both results with the start of the first cardiac cycle.

In Figures 5A–D, we present the two-dimensional velocity profiles extracted from the simulation at the four measurement planes. These profiles were extracted at the peak systole of the second cycle, corresponding to a velocity at the inlet of approximately 1.42 m/s. The figures show how the profile changes along the flow through the MCA. Compared to the inflow profile, the velocity profile at 59 mm is already substantially different, as a high velocity region is visible on the left side of the artery. The profiles at 57 and 54 mm show a strong concentration of (high) velocity near the top, which is presumably due to the bend present in that region of the artery, while a bend in the opposite direction just before the 49 mm plane is the likely cause of the more evenly distributed velocities there at peak systole (Figure 5E). In Figures 5E,F we show the calculated wall shear stress (WSS) across the MCA at peak systole and diastole (at 2.18 s). The front in Figures 5E,F corresponds to the left side in Figures 5A–D. We observe a WSS of >40 Pa during the systole in at least three locations. The WSS at the subsequent diastole (Figure 5F) remains relatively high at the location near the second outlet at the top, which indicates that this location could be susceptible to the formation of a new aneurysm.

FIGURE 5

Figure 5. Calculated flow velocity magnitude, in the direction along the vessel center lines, at the second peak systole (at 1.44 s) in m/s for each of the four TCD validation planes. We show the velocity profiles in (A–D) for planes at a depth at 49, 54, 57, and 59 mm, respectively (run at 100% velocity, Newtonian rheology). We present the calculated wall shear stress (WSS) at peak systole in (E), and at diastole (at 2.18 s) in (F) (using the same scale).

3.2. Full vs. Reduced Velocity Simulations

In this section we compare the velocity profiles at peak systole from simulations at 10 μm voxel size and full velocity with those at reduced velocity and/or increased voxel size. Reduced velocity and resolution runs are attractive because they are cheaper, faster to run, and more likely to become computationally tractable in a clinical setting. For example, at time of writing, a full velocity run across five cardiac cycles costs approximately £4200 on the ARCHER supercomputer (EPCC, 2017), whereas a run at 50% velocity and the same resolution costs £2100 and a run at 50% velocity and 20 μm voxel size costs £500 to perform. However, reduced velocity simulations have a lower Reynolds number which affects a wide range of flow properties. In this study we have performed runs at 50% velocity (Re ~ 483) and 25% velocity (Re ~ 242).

We compare our simulation results at full velocity and resolution with those at reduced velocity and resolution in Figure 6 and Table 3. When we reduce the inflow velocity by 50%, the maximum inflow velocity at the inlet is 0.75 m/s (not an uncommon value for healthy volunteers) (Bishop et al., 1986) instead of 1.50 m/s (not an uncommon value for stroke patients) (Manno et al., 1998). We multiply the extracted velocities from our reduced velocity runs by two for simulations at 50% inflow velocity, and by four for simulations at 25% velocity. When comparing the runs with full inflow velocity runs with those at 50%, we already observe major differences in the extracted velocities. Here the comparisons at all four locations feature velocity differences of more than 0.4 m/s, and more than 30% of the maximum absolute flow velocity extracted in the corresponding plane. For the planes at 49 and 57 mm we see very large velocity differences near the vessel wall. This is likely due to the much higher Reynolds number of the flow in the full velocity run. When we compare the rescaled 50% velocity runs to the TCD measurements, the velocities differ by up to 15.5%, which is almost twice as large as the 8.8% maximum difference between TCD measurements and full velocity runs.

FIGURE 6

Figure 6. Absolute difference in flow velocity, between the run with Newtonian rheology at 10 μm resolution and 100% velocity and other runs for each of the four validation planes. Comparisons are made with runs at 10 μm and 50% velocity (Left column), 20 μm and 50% velocity (Middle), and 20 μm, and 25% velocity (Right) respectively. The velocities in reduced velocity runs are multiplied by 2 (for the 50% velocity runs) or 4 (for the 25% velocity runs). The snapshots were made at the second peak systole (at 1.44 s), with differences in m/s.

TABLE 3

Table 3. Comparison of full and reduced velocity simulations against TCD velocity measurements.

The results of the 50% velocity run with 20 μm voxel size are almost identical to the one with 10 μm voxel size, with only very small differences in all the velocity planes. However, the run with 25% velocity is considerably less accurate, with absolute velocity differences up to 0.75 m/s, in particular close to the vessel walls. These errors are still smaller close to the inflow boundary at 59 mm, but dominate the overall result in the validation planes that are beyond the bifurcation with lenticulostriate arteries.

We conclude that simulations with reduced velocities affect the accuracy of the results significantly. This is particularly important because realistic velocities close to the wall are essential to obtain accurate wall shear stress estimates. We find that no such estimates can be reliably made for half velocity runs.

3.3. Comparing Rheology Models

To compare different rheology models, we performed simulations on our MCA geometry using a Carreau-Yasuda (CY) rheology model (Abraham et al., 2005). When the CY model was adopted, Bernabeu et al. found important differences in the WSS and Three-Band-Diagram analysis outcome for the MCA under “healthy human” flow conditions. Here we focus on differences in velocities obtained from the two rheology models, as we are interested in comparing simulation predictions to TCD measurements.

The difference in flow velocity between the Newtonian rheology model and the CY rheology model at peak systole is shown in Figure 7. We observe differences in velocity of up to 0.12 m/s in three of the four validation planes, and a difference of up to 0.21 m/s in a highly concentrated central region in the 54 mm measurement plane. In all cases the velocity differences are largest in regions where the absolute velocity is relatively small in the Newtonian rheology results, cf. Figure 5, while only smaller differences exist in regions where the velocity is relatively large. These results suggest that the choice of using either a CY or Newtonian rheology model has little effect on $v_{\mathrm{\text{max}}}^{\mathrm{\text{sim}}}$ in all our comparisons (see Table 3).

FIGURE 7

Figure 7. Absolute difference in flow velocity, between the run with Newtonian rheology and run with CY rheology. Both of these runs were performed at 100% of the full velocity. The snapshots were made at the second peak systole (at 1.44 s), with differences in m/s, for each of the four TCD validation planes.

The difference between the Newtonian and the CY rheology model for 50% reduced velocity is shown in Figure 8 at peak systole. As noted above, velocity extractions from runs at 50% velocity are multiplied by 2 to enable a direct comparison with full velocity runs. The difference in velocities between the 50% runs is considerably smaller than for 100% velocity runs, reaching at most 0.05 m/s in any of the measurement planes. The velocity difference is largest close to the arterial wall, but is in all cases much smaller than the velocity mismatch introduced by the velocity reduction (see Figure 6, left row). This is in line with the finding that the choice of the rheology model has a small effect, and in the reduced velocity runs the impact of scaling down the velocity on accuracy is the dominating factor.

FIGURE 8

Figure 8. Absolute difference in flow velocity, between the run with Newtonian rheology and run with CY rheology. Both of these runs were performed at 50% of the full velocity, with the differences rescaled by a factor 2. The snapshots were made at the second peak systole (at 1.44 s), with differences in m/s, for each of the four TCD validation planes.

3.4. Limitations of our Study

The main limitations of our validation study are related to data acquisition, model construction and simulation constraints.

Regarding TCD measurements, phase misalignments are common when directly comparing simulation results to these measurements, due to differences in apparent cardiac cycle length between the consecutively measured TCD planes (see Figure 4). Furthermore, due to the proprietary nature of the TCD numerical data, numerical velocity values were extracted semi-automatically from JPEG images obtained with the Doppler BoxX software, which may introduce small transcription errors of up to 0.0064 m/s due to the resolution of the images. The measurement quality and level of background noise can vary with different measurements, as different depths are subject to varying levels of occlusion and wave propagation interference.

In the area of segmentation it is particularly challenging to accurately reproduce the small lenticulostriate arteries originating near the origin of the MCA (Kang et al., 2012). These arteries are not always clearly captured in the medical imaging data, and many existing haemodynamics models of MCAs do not include them, while our geometry contains two of these arteries. However, omitting them altogether can lead to velocity overestimations in the remainder of the MCA. In our model we were able to resolve the lenticulostriate arteries to a limited extent after extensive segmentation efforts.

Due to the one-dimensional nature of the TCD measurement, we used a parabolic inflow velocity profile and fitted it to the non-circular shapes of the inflow boundaries (see section 2). Real inflow velocity profiles can vary depending on the morphology of the arterial network, as shown for example by Takeuchi and Karino (2010). Regarding the outlets, a more physiologically accurate choice of boundary condition would take into account the downstream peripheral resistance. However, such an approach introduces additional patient-specific parameters. For the purposes of the validation conducted in this study we intentionally limit the complexity of the model and thus use a simple mixed Neumann-Dirichlet boundary condition.

Furthermore, our simulation model is based on a rigid geometry and does not include elastic deformations of the vessel. In the case of blood flow in cerebral aneurysms, Dempere-Marco et al. (2006) found that incorporating wall motion has relatively little effect on the WSS. Understanding the dynamical response of arterial walls in the MCA, on a patient-specific level, is a particularly challenging area of research. However, recent studies show promising results that should soon allow us to examine these processes (Oubel et al., 2010; Vanrossomme et al., 2015).

3.5. Future Work

There are a range of factors that we seek to incorporate in our model as part of our future research. Firstly, we aim to develop techniques to create more realistic inflow profiles by using simulation data of arteries upstream from the patients MCA. Secondly, we seek to enhance our model by incorporating mechanisms for arterial wall deformations and damage. Such mechanisms are highly complex and very difficult to measure experimentally, and therefore modelling them is a particularly challenging area of research. Thirdly, we seek to provide more realistic outflow properties by extending our geometry to arteries further downstream. This could be accomplished for example by investigating how existing (1D) resistance models could be accurately applied within the context of complex 3D simulation models, or by attempting to simulate the full human brain in 3D over realistic time scales, and using patient-specific flow conditions.

4. Conclusions

We have conducted a validation study comparing flow velocities from patient-specific lattice-Boltzmann simulations to clinical TCD measurements in the MCA. As part of the study, we analyzed simulation results obtained at reduced velocities and variable resolution. Moreover, we investigated the impact of using the Carrueau-Yasuda rheology model compared to a Newtonian rheology model.

We achieved very good agreement of the maximum velocity between our full patient-specific velocity simulation results and TCD measurements, with an error of less than 9% independent of the choice of rheology model. Simulating blood flow at reduced velocities, for example by scaling down the velocity or using velocity measurements from healthy subjects, is attractive because the simulation runs are computationally cheaper and deliver results faster. However, we found that scaling down the flow velocities leads to substantially larger errors, and an accurate comparison between simulations and TCD measurements is no longer achieved. Adopting a CY rheology model instead of a Newtonian one results in small changes in maximum velocities in the planes and in our validation, whereas substantial flow velocity differences are observed near the arterial wall and in the resulting WSS. However, the CY rheology model does not enable a significant improvement when the velocity is already scaled down (e.g., by using inflow profiles of healthy volunteers or reduced velocity Womersley profiles), as errors caused by this velocity scaling then dominate the overall accuracy. Figures 7, 8 suggest that a Newtonian rheology model may be a justifiable approximation for MCA simulations at lower (i.e., < 0.75 m/s) peak flow, but that this could quickly become problematic for the higher flows typically recorded in unhealthy patients (in which 1.5 m/s is not unusual).

Computational haemodynamics predictions that accurately match patient-specific TCD measurements are likely to become an important asset in clinical settings and pave the way to using computer models in the process of clinical decision making (Fenner et al., 2008; Sadiq et al., 2008). Compared to clinical measurements alone, patient-specific simulations allow us to obtain information about a much wider range of flow properties, such as detailed flow velocity characteristics and wall shear stress estimates. In addition, simulations can help predict flow velocity in areas that have not been directly measured, and thereby help reduce the number and intensity of invasive measurements that need to be performed. Here we have shown that a combination of non-invasive TCD measurements with haemodynamics simulations can lead to accurate predictions of blood flow velocity throughout the MCA. The ability to make these accurate predictions constitutes an important step in making computational haemodynamics a viable approach for assessing intracranial blood flow.

Ethics Statement

The patient-specific data (3D angiography, TCD measurements) used in this study was available on the shelf and did not contain identifiable information. The segmented geometry is already published Itani et al. (2015). The present study involved only secondary analysis of de-identified data that is not linked to the subjects from whom it was originally collected.

Author Contributions

DG conceived the study, while DG and RR carried out the simulations, performed the validation comparison, and wrote the manuscript. US advised on the choice of simulation parameters and contributed to writing the manuscript. RC segmented the medical images and extracted the TCD velocity profiles from the measurement images, with help from DG, US, and PC. FR obtained the original angiography images, while HC performed the TCD measurements. Both FR and HC advised on the medical aspects of the manuscript. PC coordinated the study and helped draft the manuscript. All authors gave final approval for publication.

Funding

This work received funding from the EU FP7 CRESTA project (grant number, 287703), the EU H2020 projects ComPat (grant no. 671564) and CompBioMed (grant no. 675451), the EPSRC-funded 2020 Science Programme (EP/I017909/1), and the Qatar National Research Fund (NPRP), project No. 5-792-2-238. RC is supported through doctoral training grant SP/08/004 from the British Heart Foundation (BHF), through the UCL CoMPLEX doctoral training programme.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Aditya Jitta for his efforts in analyzing a preliminary version of our MCA simulation model. We are grateful to Rupert Nash, Miguel Bernabeu, and Timm Krüger for useful discussions pertaining the simulation parameters. We thank Ann Warner for her help in obtaining angiographic data, and Sebastian Schmieschek in assisting with the supervision of RC. Access to the ARCHER supercomputer at EPCC in the UK was provided by the UK Consortium on Mesoscopic Engineering Sciences (EP/L00030X/1). In addition, we are grateful to the Leibniz Rechenzentrum (LRZ) in Germany for providing access to the SuperMUC supercomputer.

Supplemental Data

We provide the source data for our publication via Figshare under doi: 10.17633/rd.brunel.5001962.

References

Abraham, F., Behr, M., and Heinkenschloss, M. (2005). Shape optimization in steady blood flow: a numerical study of non-newtonian effects. Comput. Methods Biomech. Biomed. Eng. 8, 127–137. doi: 10.1080/10255840512331388588

PubMed Abstract | CrossRef Full Text | Google Scholar

Alnæs, M. S., Isaksen, J., Mardal, K. A., Romner, B., Morgan, M. K., and Ingebrigtsen, T. (2007). Computation of hemodynamics in the circle of Willis. Stroke 38, 2500–2505. doi: 10.1161/STROKEAHA.107.482471

PubMed Abstract | CrossRef Full Text

Bernabeu, M. O., Nash, R. W., Groen, D., Carver, H. B., Hetherington, J., Krüger, T., et al. (2013). Impact of blood rheology on wall shear stress in a model of the middle cerebral artery. J. R. Soc. Interface Focus 3:20120094. doi: 10.1098/rsfs.2012.0094

PubMed Abstract | CrossRef Full Text | Google Scholar

Bernsdorf, J., and Wang, D. (2009). Non-Newtonian blood flow simulation in cerebral aneurysms. Comput. Math. Appl. 58, 1024–1029. doi: 10.1016/j.camwa.2009.02.019

CrossRef Full Text | Google Scholar

Bhatnagar, P. L., Gross, E. P., and Krook, M. (1954). A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Phys. Rev. 94, 511–525.

Google Scholar

Bishop, C. C., Powell, S., Rutt, D., and Browse, N. L. (1986). Transcranial Doppler measurement of middle cerebral artery blood flow velocity: a validation study. Stroke 17, 913–915.

PubMed Abstract | Google Scholar

Boussel, L., Rayz, V., McCulloch, C., Martin, A., Acevedo-Bolton, G., Lawton, M., et al. (2008). Aneurysm growth occurs at region of low wall shear stress. Stroke 39, 2997–3002. doi: 10.1161/STROKEAHA.108.521617

PubMed Abstract | CrossRef Full Text

Bouzidi, M., Firdaouss, M., and Lallemand, P. (2001). Momentum transfer of a Boltzmann-lattice fluid with boundaries. Phys. Fluids 13, 3452–3459. doi: 10.1063/1.1399290

CrossRef Full Text | Google Scholar

Boyd, J., Buick, J. M., and Green, S. (2007). Analysis of the Casson and Carreau-Yasuda non-Newtonian blood models in steady and oscillatory flows using the lattice Boltzmann method. Phys. Fluids 19:093103. doi: 10.1063/1.2772250

CrossRef Full Text | Google Scholar

Byrne, G., Mut, F., and Cebral, J. (2014). Quantifying the large-scale hemodynamics of intracranial aneurysms. Am. J. Neuroradiol. 35, 333–338. doi: 10.3174/ajnr.A3678

PubMed Abstract | CrossRef Full Text | Google Scholar

Castro, M. A., Putman, C. M., and Cebral, J. R. (2006). Computational fluid dynamics modeling of intracranial aneurysms: effects of parent artery segmentation on intra-aneurysmal hemodynamics. Am. J Neuroradiol. 27, 1703–1709.

PubMed Abstract | Google Scholar

Cebral, J. R., Mut, F., Weir, J., and Putman, C. (2011). Quantitative characterization of the hemodynamic environment in ruptured and unruptured brain aneurysms. Am. J. Neuroradiol. 32, 145–151. doi: 10.3174/ajnr.A2419

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, W., Landgraf, B., Johnson, K. M., Kecskemeti, S., Wu, Y., Velikina, J., et al. (2011). Velocity measurements in the middle cerebral arteries of healthy volunteers using 3d radial phase-contrast HYPRFlow: comparison with transcranial Doppler sonography and 2d phase-contrast MR imaging. Am. J. Neuroradiol. 32, 54–59. doi: 10.3174/ajnr.A2240

PubMed Abstract | CrossRef Full Text | Google Scholar

Chatzizisis, Y. S., Jonas, M., Coskun, A. U., Beigel, R., Stone, B. V., Maynard, C., et al. (2008). Prediction of the localization of high-risk coronary atherosclerotic plaques on the basis of low endothelial shear stress:: an intravascular ultrasound and histopathology natural history study. Circulation 117, 993–1002. doi: 10.1161/CIRCULATIONAHA.107.695254

PubMed Abstract | CrossRef Full Text | Google Scholar

Dempere-Marco, L., Oubel, E., Castro, M., Putman, C., Frangi, A., and Cebral, J. (2006). “CFD analysis incorporating the influence of wall motion: application to intracranial aneurysms,” in International Conference on Medical Image Computing and Computer-Assisted Intervention (Copenhagen: Springer), 438–445.

Google Scholar

EPCC (2017). ARCHER kAU Calculator. Available online at: http://www.archer.ac.uk/access/au-calculator/

Fenner, J. W., Brook, B., Clapworthy, G., Coveney, P. V., Feipel, V., Gregersen, H., et al. (2008). The EuroPhysiome, STEP and a roadmap for the virtual physiological human. Philos. Trans. R. Soc. Lond. A Math. Phys. Eng. Sci. 366, 2979–2999. doi: 10.1098/rsta.2008.0089

PubMed Abstract | CrossRef Full Text | Google Scholar

Groen, D., Hetherington, J., Carver, H. B., Nash, R. W., Bernabeu, M. O., and Coveney, P. V. (2013). Analyzing and modeling the performance of the HemeLB lattice-Boltzmann simulation environment. J. Comput. Sci. 4, 412–422. doi: 10.1016/j.jocs.2013.03.002

CrossRef Full Text | Google Scholar

Hassan, T., Ezura, M., Timofeev, E. V., Tominaga, T., Saito, T., Takahashi, A., et al. (2004). Computational simulation of therapeutic parent artery occlusion to treat giant vertebrobasilar aneurysm. Am. J. Neuroradiol. 25, 63–68.

PubMed Abstract | Google Scholar

Itani, M. A., Schiller, U. D., Schmieschek, S., Hetherington, J., Bernabeu, M. O., Chandrashekar, H., et al. (2015). An automated multiscale ensemble simulation approach for vascular blood flow. J. Comput. Sci. 9, 150–155. doi: 10.1016/j.jocs.2015.04.008

CrossRef Full Text | Google Scholar

Jain, K., Jiang, J., Strother, C., and Mardal, K.-A. (2016). Transitional hemodynamics in intracranial aneurysmscomparative velocity investigations with high resolution lattice boltzmann simulations, normal resolution ansys simulations, and MR imaging. Med. Phys. 43, 6186–6198. doi: 10.1118/1.4964793

PubMed Abstract | CrossRef Full Text | Google Scholar

Jou, L.-D., Lee, D. H., Morsi, H., and Mawad, M. E. (2008). Wall shear stress on ruptured and unruptured intracranial aneurysms at the internal carotid artery. Am. J. Neuroradiol. 29, 1761–1767. doi: 10.3174/ajnr.A1180

PubMed Abstract | CrossRef Full Text | Google Scholar

Kang, C.-K., Woerz, S., Liao, W., Park, C.-A., Kim, Y.-B., Park, C.-W., et al. (2012). Three dimensional model-based analysis of the lenticulostriate arteries and identification of the vessels correlated to the infarct area: Preliminary results. Int. J. Stroke 7, 558–563. doi: 10.1111/j.1747-4949.2011.00611.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Khan, M., Steinman, D., and Valen-Sendstad, K. (2016). Non-Newtonian versus numerical rheology: Practical impact of shear-thinning on the prediction of stable and unstable flows in intracranial aneurysms. Int. J. Numer. Method Biomed. Eng. 33:e2836. doi: 10.1002/cnm.2836

PubMed Abstract | CrossRef Full Text | Google Scholar

Ladd, A. J. C. (1994). Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation. J. Fluid Mech. 271, 285–309.

Google Scholar

Manno, E. M., Gress, D. R., Schwamm, L. H., Diringer, M. N., and Ogilvy, C. S. (1998). Effects of induced hypertension on transcranial Doppler ultrasound velocities in patients after subarachnoid hemorrhage. Stroke 29, 422–428.

PubMed Abstract | Google Scholar

Marzo, A., Singh, P., Larrabide, I., Radaelli, A., Coley, S., Gwilliam, M., et al. (2011). Computational hemodynamics in cerebral aneurysms: the effects of modeled versus measured boundary conditions. Ann. Biomed. Eng. 39, 884–896. doi: 10.1007/s10439-010-0187-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Meckel, S., Leitner, L., Bonati, L. H., Santini, F., Schubert, T., Stalder, A. F., et al. (2013). Intracranial artery velocity measurement using 4D PC MRI at 3 T: comparison with transcranial ultrasound techniques and 2D PC MRI. Neuroradiology 55, 389–398. doi: 10.1007/s00234-012-1103-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Miura, Y., Ishida, F., Umeda, Y., Tanemura, H., Suzuki, H., Matsushima, S., et al. (2013). Low wall shear stress is independently associated with the rupture status of middle cerebral artery aneurysms. Stroke 44, 519–521. doi: 10.1161/STROKEAHA.112.675306

PubMed Abstract | CrossRef Full Text | Google Scholar

Mountrakis, L., Lorenz, E., and Hoekstra, A. G. (2013). Where do the platelets go? a simulation study of fully resolved blood flow through aneurysmal vessels. Interface Focus 3:20120089. doi: 10.1098/rsfs.2012.0089

PubMed Abstract | CrossRef Full Text | Google Scholar

Nash, R. W., Carver, H. B., Bernabeu, M. O., Hetherington, J., Groen, D., Krüger, T., et al. (2014). Choice of boundary condition for lattice-boltzmann simulation of moderate-Reynolds-number flow in complex domains. Phys. Rev. E 89:023303. doi: 10.1103/PhysRevE.89.023303

PubMed Abstract | CrossRef Full Text | Google Scholar

Ouared, R., Larrabide, I., Brina, O., Bouillot, P., Erceg, G., Yilmaz, H., et al. (2016). Computational fluid dynamics analysis of flow reduction induced by flow-diverting stents in intracranial aneurysms: a patient-unspecific hemodynamics change perspective. J. Neurointerv. Surg. 8, 1288–1293. doi: 10.1136/neurintsurg-2015-012154

CrossRef Full Text | Google Scholar

Oubel, E., Cebral, J., De Craene, M., Blanc, R., Blasco, J., Macho, J., et al. (2010). Wall motion estimation in intracranial aneurysms. Physiol. Meas. 31:1119. doi: 10.1088/0967-3334/31/9/004

PubMed Abstract | CrossRef Full Text | Google Scholar

Sadiq, S. K., Mazzeo, M. D., Zasada, S. J., Manos, S., Stoica, I., Gale, C. V., et al. (2008). Patient-specific simulation as a basis for clinical decision-making. Philos. Trans. A Math. Phys. Eng. Sci. 366, 3199–3219. doi: 10.1098/rsta.2008.0100

PubMed Abstract | CrossRef Full Text | Google Scholar

Succi, S. (2001). The Lattice Boltzmann Equation: For Fluid Dynamics and Beyond. Oxford: Oxford University Press.

Google Scholar

Takeuchi, S., and Karino, T. (2010). Flow patterns and distributions of fluid velocity and wall shear stress in the human internal carotid and middle cerebral arteries. World Neurosurg. 73, 174–185. doi: 10.1016/j.surneu.2009.03.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Vanrossomme, A. E., Eker, O. F., Thiran, J.-P., Courbebaisse, G. P., and Zouaoui Boudjeltia, K. (2015). Intracranial aneurysms: wall motion analysis for prediction of rupture. Am. J. Neuroradiol. 36, 1796–1802. doi: 10.3174/ajnr.A4310

PubMed Abstract | CrossRef Full Text | Google Scholar

Venugopal, P., Valentino, D., Schmitt, H., Villablanca, J. P., Vinuela, F., and Duckwiler, G. (2007). Sensitivity of patient-specific numerical simulation of cerebal aneurysm hemodynamics to inflow boundary conditions. J. Neurosurg. 106, 1051–1060. doi: 10.3171/jns.2007.106.6.1051

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiang, J., Tremmel, M., Kolega, J., Levy, E. I., Natarajan, S. K., and Meng, H. (2011). Newtonian viscosity model could overestimate wall shear stress in intracranial aneurysm domes and underestimate rupture risk. J. Neurointerv. Surg. 4, 351–357. doi: 10.1136/neurintsurg-2011-010089

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiang, J., Tutino, V. M., Snyder, K. V., and Meng, H. (2014). CFD: computational fluid dynamics or confounding factor dissemination? The role of hemodynamics in intracranial aneurysm rupture risk assessment. Am. J. Neuroradiol. 35, 1849–1857. doi: 10.3174/ajnr.A3710

PubMed Abstract | CrossRef Full Text | Google Scholar