The present invention relates to a method for calculating pressures in a fluid stream through a tube section, especially a blood vessel with atherosclerotic plaque.
Atherothrombosis is the leading cause of death and severe disability worldwide (Rosamond et al., 2007). The disease generates fatty deposits inside vessel walls (lipid cores) covered by a protective fibrous cap—the atherosclerotic plaque. If the cap ruptures, blood clots are formed which may be carried downstream to lodge in small-diameter vessels. Disrupted blood flow results in causing heart attacks or strokes. Currently, the risk of cap rupture is assessed using the degree of luminal narrowing. This fails to take the morphology of the plaque into account. Indeed, unstable or vulnerable plaques are known to possess large lipid cores and thin fibrous caps. Morphology such as this generates severe internal stresses in the fibrous cap. In vitro studies have shown that cap rupture predominantly occur when static stresses exceed 300 kPa. The ability to estimate stress and strain magnitudes in the fibrous cap is thus expected to improve risk assessment.
Two large studies, the North American Symptomatic Carotis Endarterectomy Trial (NASCET, 1991) and the European Carotid Surgery Trial (ECST, 1991) have proved carotid endarterectomy to be highly beneficial in symptomatic patients with luminal narrowing of 70% or greater (Rothwell et al., 2003), henceforth designated severe stenosis. The Asymptomatic Carotid Atherosclerosis Study (ACAS, 1995) and the Asymptomatic Carotid Surgery Trial (Halliday et al., 2004) demonstrated that even asymptomatic patients with severe degrees of stenosis may benefit from surgical treatment. Therefore, current risk assessment of carotid atherosclerotic plaques is done by evaluating the degree of luminal stenosis, as measured by conventional intra-arterial angiography, duplex ultrasound, magnetic resonance angiography, or computational tomography.
While increasing degrees of stenosis carries inherent greater risks of atherothrombotic events, most ruptured coronary plaques are less than 60% stenosed (Casscells et al., 2003). As atherosclerosis develops, arteries often dilate in an attempt to normalize elevated wall shear stresses, a process known as expansive remodelling (Glagov et al., 1987). Consequently, risk assessment based purely on the degree of luminal stenosis will tend to underestimate the volume and thereby the severity of the atherosclerotic lesion. Therefore, there is a need for more sensitive methods to risk stratify patients with carotid atherosclerotic disease.
Thrombosis-prone (vulnerable) atherosclerotic plaques have been determined by histological examinations to possess large lipid pools and thin protective fibrous caps (Naghavi et al., 2003). Since plaque rupture, the most frequent cause of thrombosis, by definition represents a structural failure of the protective fibrous cap, it seems reasonable to assume that plaque morphology and biomechanical properties of the atherosclerotic lesion may influence the vulnerability. An in-vitro model study (Li et al., 2006b) showed that thin fibrous caps in combination with large lipid pools generate severe internal fibrous cap stresses. In vitro studies of coronary arteries have shown markedly elevated fibrous cap stresses in ruptured compared to stable lesions (Cheng et al., 1993). A recent publication comparing fibrous cap stress levels in asymptomatic versus symptomatic carotid patients found stresses in symptomatic patients to be nearly twice those of asymptomatic patients (Li et al., 2007).
To estimate stress levels in the fibrous cap, fluid structure interaction (FSI) analysis has emerged as a promising tool combining blood flow simulation through computational fluid dynamics (CFD) with finite element analysis (FEA) of the corresponding stress levels in the surrounding tissues (Tang et al., 2004). Magnetic Resonance Imaging (MRI) is a comprehensive imaging modality for providing geometric models for FSI simulations given its unique in vivo soft tissue imaging capabilities (Saam et al., 2005) combined with the ability to determine blood velocities (Chai and Mohiaddin, 2005) and vessel wall compliance (Metafratzi et al., 2002) non-invasively.
In clinical tests for risk assessment, it is highly desired to have evaluation methods that are fast and reliable. Performing 3D simulations of fibrous cap stresses is very time consuming, and therefore not applicable for rapid characterization of a patient's risk for thrombosis. Therefore, simulations in 2D cross-sections corresponding to either histological data (Cheng et al., 1993) or MRI scans (Li et al., 2006a) have been suggested. Even though the use of cross-sectional data matches the orientation of the available morphologic data, this approach precludes examinations of the longitudinal stress distribution. A previous study (Lovett and Rothwell, 2003) determined carotid plaque ruptures to be asymmetrically distributed longitudinally with a majority of ruptures occurring on the upstream side of the plaques. Furthermore, histological examinations of the longitudinal distribution of cell-types indicate differences in smooth muscle cell and macrophage content between upstream and downstream parts of the plaque (Dirksen et al., 1998).
A faster evaluation method has been disclosed in International patent application WO2007/002821. The method uses images of cross-sections of the artery or other vessel of interest to identify and locate components of the atherosclerotic deposit, including any hemorrhage, necrotic core, and calcification, and to determine the status and composition of the fibrous cap. A simple scoring system is applied that accounts for the presence of these components and more heavily weights the presence of these components in the juxtaluminal portion of the deposit. This method suffers from the fact that each image has to be evaluated by a person in order to assign values to the degree of occurrence of hemorrhage, necrotic core, and calcification. This human evaluation introduces uncertainties because different persons may judge differently form identical images, which is a great disadvantage.
International patent application WO 2006/069379 discloses a technology to identify structural parameters in order to determine distances, diameters, surface area, volume, etc. This includes the degree of stenosis. However, there is not given any advice as how to calculate pressure variations in the blood vessels for use in simulations of stress distribution around surface irregularities.
In U.S. Pat. No. 6,377,832 by Bergmann et al. discloses a method for neural networks in a calculation of the severity of a stenosis in a blood vessel.
In International patent application WO 2006/002353 by Chen and Dwyer, a method is disclosed for determining the risk of rupture in a blood vessel using 2-D slice images by scanning a blood vessel and generating a 3-D mesh model of the blood vessel. This method is complicated and requires high computing capacity due to the evaluation of the 3-D model. It would be desirable to provide a simpler model requiring less computing capacity.
It is the object of the invention to provide a simple, automated method and system for calculation of pressure variations in fluid tubes for use in simulations of longitudinal and radial stress distribution around surface irregularities. Especially, it is the object of the invention to provide specific methods and systems for simulation of stress distribution and blood pressure in blood vessels in order to facilitate predictions of plaque rupture risk.
This purpose is achieved with a method according to the invention for calculating pressures in a fluid streaming through a tube section from an upstream end to a downstream end of the tube section. The tube section comprises at least one flow divider dividing the fluid stream from an inlet to at least two arms with two outlets. The method comprises a first step of scanning the tube section with a scanner and providing a plurality of cross-sectional 2D scanning images along the tube section. By a computer program on the basis of these 2D images, a 3D image of the tube section is automatically calculated in a computer by using interpolation between the 2D images. Also, by using a computer, a 2D sectional image cut through the 3D image is performed, wherein the image cut follows the fluid stream in the tube section with the arms. Then, in the sectional image cut, a fluid pressure distribution is calculated along multiple locations along the tube section on the basis of given boundary conditions, the boundary conditions including fluid velocity or fluid pressure at the upstream end. Optionally, local fluid pressure extracted from the sectional cut is imposed on individual cross-section in order to calculate radial parameters.
By the method according to the invention, the pressure distribution in blood vessels can be calculated in a computerised manner. Once having the pressure distribution, methods can be applied for how to simulate stress distribution around surface irregularities. The method is relatively simple and rapid and, therefore, suitable for clinical tests.
In comparison to the aforementioned International Patent Application WO 2006/002353, the method according to the invention is simplified by the 2D sectional image cut following the fluid stream in the tube section with the arms. Thus, the 3-D image is reduced to a 2-D basis, for which calculations are performed fast, simple and reliable even with low computer capacity.
The term cross-sectional is to be understood as a 2D slicing of the tube section along the flow direction in the tube but no necessarily normal to the fluid direction. A certain deviation from a lateral direction, for example a scanning at an angle of between 70° and 89° or between 91° and 110° instead of scanning at 90° works equally well and can be taken account of in the calculations.
Though the method is highly suited for simulation of the stress and strain distribution and blood pressure in blood vessels in order to facilitate predictions of plaque cap rupture risk, the method is of general character and may be used for fluid flow in tubes in general, for example water pipes and oil pipes. Especially, it may be used in tubes with material deposits on the inner tube walls.
In a first embodiment, the method implies measuring the fluid pressure at the upstream end and measure the fluid velocity at the downstream end and include these measurements in the boundary conditions. For example, the method implies measuring the fluid velocity at the upstream end, and calculating on the basis of the pressure distribution in the tube section the expected fluid velocity at the upstream end, and comparing the measured fluid velocity and the expected fluid velocity at the upstream end as part in an evaluation procedure for the validity of the pressure calculation.
In a second embodiment, the method implies measuring the fluid velocity at the upstream end and wherein the boundary conditions include the measured fluid velocity at the upstream end and an estimated pressure at the downstream end. For example, the method implies measuring the fluid velocity at the downstream end, and calculating on the basis on the pressure distribution in the tube section the expected fluid velocity at the downstream end, and comparing the measured fluid velocity and the expected fluid velocity at the downstream end as part in an evaluation procedure for the validity of the pressure calculation.
In a third embodiment, the method implies using the derived longitudinal pressure distribution achieved by either the first or second embodiment as a boundary condition for calculating stresses and deformations in cross-sectional 2D models corresponding to the available image scans, for example MRI scans. Though both the first and second embodiments may be used for longitudinal pressure calculations it would be preferable for the second embodiment to be used to this method's inherent greater accuracy.
Advantageously, the method implies a comparison of the calculated deformations in cross-sectional 2D models with provided cross-sectional scans of the tube section, where the scans depict actual deformations in different cross sectional planes. By comparing the calculated cross sectional deformation with the actual deformation, validation of the chosen tissue parameters can be performed and the accuracy of the simulation can be checked. Preferably, the scans are dynamical scans, however, another possibility is pulsed scans at maximum blood pressure for maximum deformation. In addition, for blood vessels, MRI scans are preferred, however, other scanning methods are also possible, for example ultrasound or X-ray.
As in the case of atherosclerotic plaques in blood vessels or in cases where other kinds of material deposit lumps are found on the inner wall of the tube section obstructing part of a flow path of the fluid in the tube section, the method comprises calculating the forces on the surface of the lump on the basis of the pressure in the tube section at the location of the lump and on the basis of the fluid stream through the section. The calculated forces are then used to perform computer calculations for the deformation of the tube wall due to pressure on the lump. As part in an evaluation procedure for the validity of the applied mechanical parameters this deformation can be correlated with in situ measurements in cross-sectional planes.
As in the case of atherosclerotic plaques in blood vessels, the lump has a cap, and in a further step, the method implies calculating the stress on the cap on the basis of the pressure and the flow of the fluid. The calculated stresses can be assigned to both longitudinal sections and/or cross-sectional planes.
In a practical embodiment, the method comprises the finding of centre points for the fluid stream in the 3D image and calculating a 2D sectional cut through these centre points. Such a cut may not be a straight plane in the case where the tube is bent, as the plane in this case is bent in order to follow the stream. In order to facilitate the calculation, this can be overcome by projecting the image cut onto a straight plane before the calculation of the fluid pressure distribution.
The finding of the centre points may be performed automatically by centroid evaluation functions, however, very good results have been obtained by finding the centre points in the tube section in a manual selection process. In the case that the tube section, for example the blood vessel has a flow divider, which is a usual case, the centre points are found in all arms of the tube section, and a 2D sectional cut is calculated through all centre points.
In the case of blood vessels, the preferred imaging modality is magnetic resonance imaging (MRI), and the method involves measuring the velocity of the fluid in the upstream end or in the downstream end or in both by phase contrast MRI scans. In general the approach applies to any technique capable of acquiring 3D images.
In the following, a semi-automated method of creating longitudinal 2D Fluid Structure Interaction (FSI) models from MRI scans is presented. As shown, this allows simulations of longitudinal stress distributions and blood pressure levels, enables predictions of plaque rupture risk and illustrates the correlations between local stress variation and morphology.
The method according to the invention can also be used for pressure calculation in the heart, especially for evaluating the risk connected to plaque deposits in coronary arteries. In order to image the heart arteries, a method called Virtual Histology can be used, wherein an ultrasonic catheter is scanning the vessel walls from the inside for generating an image of the vessel walls.
It is fundamental to the approach to apply technology capable of separating tissues (morphologies) in the scanned region. The distinction between different cell phenotypes, the extra cellular matrix that cells produce, or the structural regions of different chemical components deposited in the tissue are central in order to generate a computational model of the mechanical stress distribution. The segmentation enables us to assign different properties to distinct regions and hereby model the tissue as a heterogeneous materials rather than a homogeneous. It is the heterogeneity of the structure that allows local peak stress to build up as blood pressure is increased.
The invention will be explained in more detail with reference to the Drawing with:
In the following three different methods will be explained in connection with the invention. The first method is a more simple method, which delivers results superior to prior art methods. However, even better results were obtained by the second method. The third method is applying results from method 1 or 2 in order to obtain new results in sections coincident with the original 2D scans.
Patients awaiting operation for severe carotid plaque can be scanned by a well-validated MRI protocol (Yuan and Kerwin, 2004) using a Philips 1.5 T scanner (Philips Intera Achieva 1.5T R1.5.4, Philips Inc., Best, The Netherlands). The plaque is scanned using a T1-weighted scan (T1W), T2-weighted scan (T2W), proton-density weighted scan (PDW) and time-of-flight (TOF) scans. Individual plaque components such as lipid, fibrous cap, and calcifications exhibits different image intensities on each scan allowing segmentation approaching histological definitions employed by AHA with a large degree of accuracy (Saam et al., 2005). Sixteen transverse slices centered on the carotid flow divider are acquired using cardiac-gated turbo spin echo sequences: T1W (TR/TE/Inv: 1RR/8/650 ms), T2W and PDW scans with TR=3RR and TE=40/20 ms, and a TOF sequence (TR/TE 34.9/2.4 ms). All sequences cover a field of view of 16×12 cm using a 256×256 matrix resulting in a raw resolution of 0.61×0.61 mm. The final images are interpolated using zero-filled Fourier transforms to a 512×512 matrix and final resolution of 0.3125×0.3125 mm with 2 mm slice thickness.
Flow velocities are measured in the transverse plane 2 cm up- and downstream from the flow divider using a phase-contrast turbo field echo sequence (TR/TENenc: 4.84 ms/2.90 ms/150 cm/s). To resolve the motion of the vessel walls, dynamic balanced steady-state free precession (SSFP) scans with 20 cardiac phases may be performed using a single longitudinal slice and five transverse slices (TR/TE/slice thickness/inplane resolution: 5.7 ms/2.8 ms/8 mm/0.74×0.76 mm and 7.7 ms/3.9 ms/2 mm/0.64×0.52 mm, respectively). To reduce artifacts stemming from vascular motion all scans should be acquired using ECG-triggering.
The image intensities of each scan compared to the adjacent sternocleid muscle can be analyzed using Cascade, a dedicated semi-automated segmentation tool, allowing segmentation into lipid core, fibrous cap, vessel wall, and blood stream (Liu et al., 2006) (
To create a longitudinal 2D model, the 3D isosurface model is sectioned along a non-uniform rational B-spline (NURBS) surface (Piegl and Tiller, 1997) (
In order to create a surface encompassing the entire model, these lines are extended a user-defined distance outwards. All four lines are then used to define a NURBS surface sectioning the central part of the blood vessel throughout the model (
Due to the large slice thickness, the flow divider will often tend to have an unphysiologically flat geometry. To correct for this abnormality, longitudinal MRI scans can be performed (
Though the demonstrated method uses a semi-automated approach with operator intervention, in many cases, an automated method for finding the vessel centre points would be sufficient. This automated method could be performed automatically by a computer program evaluating the density distribution and orientation of the identified centre points.
From the obtained lines a finite element model was developed employing smoothed splines describing the structural components. From this a coupled fluid-structure interaction (FSI) model was synthesised using COMSOL, a commercially available finite element solver (COMSOL 3.3a, COMSOL Inc, Stockholm, Sweden). The FSI simulation was performed using the FSI plane strain application mode available in the Solid Mechanics module of COMSOL. The FSI model employs a combination of flow simulation in the fluid domain described within the chemical engineering module with the solid mechanics module that simulates the mechanical response of the structure in the solid domain. The coupling is established by obtaining a continuation of the deformation and forces over the boundary that exists between the fluid and the solid domains. Details exists in the COMSOL Multiphysics User's Guide.
Blood flow is simulated in a Navier-Stokes model, and treated as an incompressible, homogeneous, Newtonian, viscous fluid with a density of 1050 kg/m3 and dynamic viscosity of 0.0035 N·s/m2. A parabolic inlet flow profile is applied with the maximal velocity Vin measured from the phase-contrast MRI scans (
Plaque components are assumed incompressible and to account for the non-linear stress/strain dependency of human tissues, a Neo-Hookean hyper-elastic model is used to specify the material properties of surrounding tissue (μ=6.20e6, κ=1.24e8, ρ=960) and vessel wall (μ=7.20e5, κ=1.44e7, ρ=1200). Lipid is treated as an isotropic materiel (Tang et al., 2004) with Young's modulus set to 1/100 of that of the equivalent Young's modulus of the vessel wall (E=1e5, ν=0.45, ρ=900).
To reduce computational demands and simulation times, different mesh densities for the individual sub-domains are applied: Vessel wall excluding the fibrous cap and lipid pool are meshed using a maximum element size of 0.3 mm. The blood stream is meshed with a maximum element size of 0.2 mm while the surrounding tissue receives a coarser mesh with a maximum element size of 1 mm. In order to ensure adequate numbers of cells across narrow regions, the fibrous cap separating lipid from bloodflow is meshed with a maximum element size of 0.1 mm.
A stationary analysis is employed to simulate the model at the time of maximal inlet velocity. Fluid velocities and corresponding plaque deformation as well as internal principal stresses are calculated and analyzed. To facilitate convergence, artificial isotropic diffusion is applied initially and stepped down logarithmically until finally being removed altogether.
First principal stress distribution and flow velocities, which are different for the two blood vessel branches, are presented in
Despite the indisputable clinical significance of high degrees of stenosis in carotid risk assessment it no longer seems appropriate to base risk stratification solely on this parameter. The invention is a novel solution describing a rapid semi-automated or automated method of determining longitudinal stress levels and blood pressures which can be used in their own right or to improve accuracy of cross-sectional simulations. The technique uses morphologic features derived from in vivo MRI scans and is as such unplagued by geometric problems resulting from histologic preparation or surgical excision.
Results indicate maximal stresses occur in the inflow shoulder region, the preferential site of plaque ulceration (Lovett and Rothwell, 2003). Interestingly, though stresses vary longitudinally, maximal stresses are neither found in the area of maximal constriction, nor in the area with the thinnest fibrous cap. Instead, stress levels were maximal at the transitional zone of high to low deformation. Performing simulations with physical characteristics of the lipid pool set to those of normal vessel wall, i.e. fibrous plaque (
The spatial resolutions of the MRI scans leave room for improvement, especially in the longitudinal directions with a slice-thickness of 2 mm. Improvement in this area could facilitate more accurate results. However, enhancing the resolution necessitates prolonged MRI scan times with added costs in terms of scanner time and patient discomfort. Going to higher field strengths (3 tesla or above) may prove beneficial trading off the greater signal to noise ratio (SNR) achieved through higher field strength with the lower SNR from thinner slices. It would have been preferable to use a dynamic solver instead of obtaining a static solution, since dynamic loading is hypothesized to be an additional risk factor. Since blood pressure is used as a boundary condition in the method above, dynamic solving requires knowledge of the blood pressure over time, which we did not have. Noninvasive measurements of dynamic blood pressure are possible using applanation tonometry (Zhao et al., 2002), originally developed for intraocular pressure measurements. It could prove of interest to include dynamically measured blood pressures in future simulations.
In the second method, the parameters were very much alike the first method as described above. However, a number of differences as explained in the following led to an improvement of the method.
Whereas in Method 1, the boundary conditions were
In both cases, the blood pressure was calculated around plaque near the flow divider, however, Method 2 yielded more accurate results. Results obtained by Method 2 are illustrated in
This method provides stress levels and deformations in 2D model planes coincident with MRI cross-sectional scans and histological sections. By using cross-sectional models, the geometry in the simulations matches the available real tissue morphology as realized by MRI or histology data, facilitating validation of the method by correlations of stress levels and known markers of plaque vulnerability.
Furthermore, whereas in method 1, geometrical errors are introduced due to the Gaussian smoothing when generating the smooth 3D shell models, this is not the case in method 3, since no smoothing needs to be performed when using method 3.
Using either cross-sectional MRI scans or histological data, cross-sectional 2D structural models can be constructed and rapidly simulated, since no flow parameters are calculated using these models. To calculate 2D cross-sectional deformation and stresses, blood pressure is applied along the fluid/solid domain interface, see
Method 3 follows the procedures outlined in method 1: model generation. The segmented grayscale images shown in
From method 1 and 2 the longitudinal pressure distribution can be derived and an average pressure calculated along horizontal cross-sections or selected from at the surface describing the fluid/solid domain interface, see
The model may be embedded in an artificial domain mimicking the surrounding tissue as depicted by the cross-sectional MRI scans. The boundaries surrounding the centrally located spine are fully constrained with respect to movements leaving the external interface free to deform according to tissue parameters.
A Neo-Hookean hyper-elastic model is used to specify the material properties of surrounding tissue (μ=6.20e6, κ=1.24e8, ρ=960) and vessel wall (μ=7.20e5, κ=1.44e7, ρ=1200). Lipid is treated as an isotropic materiel (Tang et al., 2004) with Young's modulus set to 1/100 of that of the equivalent Young's modulus of the vessel wall (E=1e5, ν=0.45, ρ=900).
To assure adequate numbers of elements across narrow regions, the resolution of narrow regions is set to 5.
The plane strain solid mechanics module in Comsol is used to simulate the stresses and deformation in the cross-sectional models resulting from the applied blood pressure derived from the longitudinal calculations of either method 1 or 2.
The deformation resulting from the mechanical cross-sectional simulations may be compared to dynamic cross-sectional in vivo MRI scans depicting the actual deformation in different planes. It is thus possible to validate the chosen tissue parameters and examine the accuracy of the simulation.
We have presented a novel technique of obtaining rapid longitudinal and radial information of plaque stress levels from 2D transverse MRI scans. The results may be used to investigate longitudinal and radial stress distributions and resulting morphologic and histologic changes as well as provide improved blood pressure boundary specifications for transverse 2D simulations. Initial results are promising conforming to established plaque physiology findings with maximal first principal stress levels found in the inflow region of the fibrous cap approaching established criteria for caps at risk of rupture (Cheng et al., 1993).
Number | Date | Country | Kind |
---|---|---|---|
PA200701555 | Oct 2007 | DK | national |
Filing Document | Filing Date | Country | Kind | 371c Date |
---|---|---|---|---|
PCT/DK2008/050266 | 10/31/2008 | WO | 00 | 7/12/2010 |
Number | Date | Country | |
---|---|---|---|
61001282 | Oct 2007 | US |