We discuss the optimal evaluation of endothelial shear stress for real-life case studies based on anatomic data acquisition. The fluid dynamic simulations require smoothing of the geometric dataset to avoid major artefacts in the flow patterns, especially in the proximity of bifurcations. A systematic series of simulations at different corrugation levels shows that, below a smoothing length of about 0.5 mm, the numerical data are insensitive to further smoothing.
Atherosclerosis is the most common cardiovascular disease and is responsible for about 35 per cent of annual deaths in developed countries . It primarily affects the arterial coronary arteries, and predictions of where and how the disease is likely to develop can be obtained by fluid dynamic simulations as a routine methodology to study blood flow patterns in human coronary arteries. A direct benefit of the simulation approach is to understand the connection between fluid-mechanical flow patterns and plaque formation and evolution, with important implications for predicting the course of atherosclerosis and possibly preventing or mitigating its effects.
The combined use of simulation and imaging techniques allows one to screen large numbers of patients non-invasively for incipient coronary disease. One option is to obtain the coronary artery wall, plaque morphology and lumen anatomy from the non-invasive multi-detector computed tomography (MDCT) imaging technique. In particular, the latest MDCT systems with 320-detector rows  enable three-dimensional acquisition of the entire coronary arterial tree in a single heart beat with the highest accuracy currently available from MDCT technology. In spite of recent technological progress, this resolution is still insufficient and the inherently noisy geometrical data pose a problem in the evaluation of endothelial shear stress (ESS), a quantity that proves extremely sensitive to the details of the wall morphology.
Although the ultimate goal of computational haemodynamics is a one-to-one comparison with medical data, such a task is still beyond reach in extended coronary systems. The reason is that the actual instrumental resolution of imaging techniques such as MDCT does not allow one to access the 0.1 mm scale and below. Previous work that compared magnetic resonance imaging (MRI)-reconstructed carotid vessels with computational fluid dynamics (CFD) results has established the lack of an unambiguous criterion to judge the match between computational and actual medical data [3,4]. More specifically, it was found that, in the case of MRI-reconstructed geometries, presenting an actual resolution of 0.1 mm, the computed shear stress shows a mean inaccuracy of approximately 15 per cent when comparing phantom data with CFD data in idealized geometries.
The main purpose of the present work is to analyse a specific aspect of this general problem, namely the convergence of computational data versus the corrugation level, for the specific case of MDCT-reconstructed data and lattice Boltzmann simulations.
2. Haemodynamics simulations
Our simulations of blood flow are based on the lattice Boltzmann method [5–7], combined with the acquisition of 320-detector MDCT data. In MDCT vascular profiling, the coronary arteries are imaged axially using a 320×0.5 mm detector configuration scanner with a rotation time of 350 ms. The geometry of each individual vessel of the coronary tree is reconstructed from the MDCT acquisition by using research-built post-processing software (VITREA 4.0, Vital Images, Minnetonka, MN, USA).
For a typical coronary artery system, the procedure to build the mesh from the MDCT raw data starts from a single vessel, formatted as stacked bi-dimensional contours (slices), with a nominal resolution of 0.1 mm. The contour path of each slice is irregular and shaped as a collection of 256 points. Each vessel is organized with the same format in a different file.
The slices are quasi-parallel and mostly transverse to the path connecting different centrelines. Moreover, the sequence of contour points is mostly aligned along the stacking sequence, that is, the contour index does not present any major twisting when moving between contiguous centrelines. Given that data relative to each vessel in the multi-branched geometry have the same format, each slice of a vessel is scanned to locate the set of mesh points that fall in its proximity and enclosed within the vessel surface. The mesh points arising from the different vessels are merged as the final step of the procedure.
As a conservative measure of the geometrical error, one can take the actual MDCT resolution to be of the order of approximately 0.4 mm . From these figures, a mild oversampling along the longitudinal direction is apparent. It also appears that the azimuthal resolution for a vessel of approximately 1 mm circumference is of the order of approximately 10−2 mm, which shows a more pronounced oversampling along the azimuthal direction.
Raw MDCT data present a level of geometric irregularities that can dramatically affect the quality of the lattice Boltzmann simulations. In order to improve the regularity of the initial geometry, we have smoothed the sequence of surface points via a linear filter along the longitudinal direction. This is achieved by averaging each point along the contour with K backward and forward slice neighbours, by using a simple weighting formula . For simplicity, we have taken wij=1/(2K+1) and studied the haemodynamics resulting from K=0 (unsmoothed geometry or level-0 smoothing), K=2 (level-2 smoothing) and K=3 (level-3 smoothing).
One could also smooth surface points along the azimuthal contour. However, the azimuthal corrugation is perceived as less strategic than the longitudinal ones, because most geometrical irregularities, such as stenosis or aneurysms, develop mostly along the longitudinal direction. Therefore, in this paper, we have studied the effect of the longitudinal smoothing procedure. We will discuss the resulting effects on the fluid dynamic results, when compared with the original, unsmoothed dataset.
An important point regards the choice of the sampling points xw used to compute the ESS. We take those points to coincide with the smoothed, off-lattice, points forming the slice contours. Each contour point is surrounded by a small number of fluid nodes and a linear extrapolation scheme is used to evaluate the ESS. Clearly, the xw points fall between the external fluid and wall mesh nodes and, as the mesh resolution increases, the set converges towards the exact no-slip hydrodynamic surface.
We have investigated a smooth human coronary artery system. According to standard naming convention, the coronary system is composed of a left main/left anterior descending (LM/LAD) vessel branching into two other primary vessels, LCX (left circumflex) and D1 (diagonal), and eight minor vessels. LM/LAD represents the longest vessel, named LM in the proximal region and LAD after the bifurcation with LCX. Overall, the test case is composed of 11 vessels, that is, one inlet vessel at the LM level and 10 outlets for all terminals. An image of this arterial system, obtained from the result of the simulation, is shown in figure 1. We have considered a mesh resolution of 0.05 mm corresponding to approximately 18 million voxels. At this resolution, the staircase approximation does not introduce major artefacts in the wall geometry, as analysed in our previous publication . In fact, the asperities in the vessel wall due to the staircase approximation account for about 5 per cent or less of the vessel diameter, while those present in the original imaging data can be as large as 50 per cent. More importantly, our previous study demonstrated convergence of the flow profiles and ESS maps at this resolution.
The effect of the wall smoothness is first inspected by looking at the sectional area and the volumetric flow profile along the full arterial system. In figure 2a,b we report the two profiles for the most representative vessel, the LM/LAD. Interestingly, the sectional area does not show major differences among the three datasets, only a higher level of noise is well visible. As a consequence, we expect to observe roughly the same flow profile for the three cases.
However, the flow profiles show marked differences between the unsmoothed and smoothed geometries. This highlights the critical role of smoothing near geometrically sensitive regions such as next to bifurcations. In fact, owing to mass conservation, the volumetric flow profile is quasi-uniform away from bifurcations, while bifurcations give rise to abrupt jumps in the volumetric flow. These jumps show large non-systematic differences as the vessel bifurcates (figure 2b). Apart from noting that the smooth geometry gives rise to a smoother flow profile, a full explanation for the sign of the difference between unsmoothed and smoothed geometries is still missing.
Overall, the level-0 profile presents a larger pressure drop between the inlet (left-most part of the profile) and outlet (right-most part).
It is reasonable to expect that wall smoothness strongly affects the ESS profile since this quantity is proportional to the gradient of the velocity field and is significantly more sensitive to resolution than basic hydrodynamic fields. Moreover, ESS is proportional to the fluid dissipation due to the fluid–wall contacts. The convergence of the flow profiles to a common profile is visible in the results of the level-2 and level-3 smoothing procedures. This is confirmed by the contour-averaged ESS, shown in figure 2c for the LM/LAD vessel and figure 2d for the LCX vessel. Again, the unsmoothed geometry shows dramatic departures from the smoothed ones. Such departures are of the order of 50 per cent for the case of a bifurcation-free vessel (figure 2d) in line with the observation of 15 per cent error with MRI-reconstructed vessels having about four times the resolution of MDCT [3,4]. The discrepancy between smoothed and unsmoothed vessels can be much larger in the case of bifurcating vessels.
By looking at volumetric and contour-averaged quantities, we conclude that our smoothing procedure provides reliable geometries for stationary flows. The complete superficial distribution of ESS, however, may still exhibit strong differences at various degrees of smoothness. This is a crucial test, since the capability of simulations to predict regions prone to atherogenesis relies on the quality of the computed ESS data. The simulation results again show strong differences between the unsmoothed and smoothed versions, while the comparison of the two smoothed geometries shows a small dispersion in ESS. A closer visual inspection confirms this result, in particular, in the proximity of bifurcations and stenosis.
To further point out the necessity of our smoothing procedure, we investigate the value of the traction vector on the surface of the artery. The traction vector t=σ⋅n−(n⋅σ⋅n)n, where σ is the stress tensor, is defined as the component of the wall force orthogonal to the wall normal n. It is a quantity of major importance in the study of haemodynamics, as it can be related to pathologies of the blood vessel, and therefore it needs to be properly reproduced in simulation.
Figure 3 displays the norm of the traction vector on the surface of a section of the artery, obtained from a simulation, once in an unsmoothed dataset and once in an artery with level-3 smoothing. The chosen artery section, marked by the rectangle on figure 1, consists of two subsequent loops. At the relatively low Reynolds number in this artery, the flow velocity of the blood is expected to be larger at the inner wall of the loop, leading to a larger wall shear stress on this wall. This property is well reproduced in the simulation at smoothing level 3. In the unsmoothed artery, on the other hand, the wall shear stress fluctuates due to the corrugation of the geometry, and fails to systematically reproduce the expected flow pattern. As a further test, a passive scalar particle was injected into the artery, and the smoothness of its trajectory used as a measure of the quality of the simulation. An extract of the particle’s trajectory is represented in figure 4 along a relatively straight section of the artery. It can be seen that, while the trajectory of the particle is perturbed in the unsmoothed artery, it is properly straightened in the smoothed case. More quantitatively, the particle trajectory xp(t) is approximated by an order-seven polynomial xp7(t), and the relative difference between the original trajectory and the smoothened polynomial is interpreted as an indication of the quality of the representation of the arterial wall. The polynomial represents a global fit over the time span of the particle over the section represented on figure 4. The order of the polynomial was chosen high enough to reproduce the overall shape of the artery on the chosen section. Consistently with previous observations, the unsmoothed case exhibits a worse result, with μp=0.0127, than the simulation with smoothing level 3, with μp=0.0081.
Summarizing, once sufficient smoothing is attained, the fine details of the ESS map show good convergence, rendering the method usable for quantitative analysis.
We have next analysed the case of a pulsatile flow. The analysis of pulsatile haemodynamics requires the investigation of the flow pattern in the time domain. Pulsatile flow is simulated by employing a time-dependent influx derived from physiological data . At the inlet, the time-dependent signal is shown in figure 5, while at the outlets we retain the constant pressure condition. We consider two representative time instants T1 and T2 and compare the flow profiles and the ESS distributions. These instants correspond to maximum Reynolds numbers of 200 and 75 at times T1 and T2, respectively. At both times, the flow profiles in figure 5 display similar features for the stationary flow, with a good convergence of profiles for the level-2 and level-3 degrees of smoothness.
The ESS data show that the unsmoothed geometry has a corrugated landscape, indicative of increased fluid–wall contacts. Interestingly, for this vessel the higher corrugation does not correspond to larger dissipation, that is, systematically larger ESS values. Overall, the ESS map bears a qualitative resemblance to that for the smoothed geometry, while a more insightful analysis of putative plaque-forming regions requires the use of the smoothed geometry.
In conclusion, the simulation of large-scale arterial systems relies on sophisticated imaging techniques and image reconstruction procedures. The imaging data from MDCT are intrinsically noisy and of little direct use in the fluid dynamic simulation of blood flow. Therefore, a smoothing procedure is needed in order to eliminate spurious corrugations in the vessel walls. This preliminary step must precede the construction of the Cartesian mesh needed for lattice Boltzmann simulations. We have shown that flow profiles and endothelial shear stress distributions contain major artefacts if the original geometry is not properly smoothed. Once the wall corrugations have been filtered out, the flow patterns show profiles that are converged and insensitive to further smoothing of the surface.
This work was supported by the Initiative in Innovative Computing at Harvard. We wish to thank Jonas Latt, Amanda Peters, Michelle Borkin, Frank Rybicky and Charles Feldman for useful discussions.
One contribution of 26 to a Theme Issue ‘Discrete simulation of fluid dynamics: methods’.
- This journal is © 2011 The Royal Society