|Home | About | Journals | Submit | Contact Us | Français|
Cardiovascular disease (CVD) is becoming the number one cause of death worldwide. Atherosclerotic plaque rupture and progression are closely related to most severe cardiovascular syndromes such as heart attack and stroke. Mechanisms governing plaque rupture and progression are not well understood. A computational procedure based on three-dimensional meshless generalized finite difference (MGFD) method and serial magnetic resonance imaging (MRI) data was introduced to quantify patient-specific carotid atherosclerotic plaque growth functions and simulate plaque progression. Participating patients were scanned three times (T1, T2, and T3, at intervals of about 18 months) to obtain plaque progression data. Vessel wall thickness (WT) changes were used as the measure for plaque progression. Since there was insufficient data with the current technology to quantify individual plaque component growth, the whole plaque was assumed to be uniform, homogeneous, isotropic, linear, and nearly incompressible. The linear elastic model was used. The 3D plaque model was discretized and solved using a meshless generalized finite difference (GFD) method. Four growth functions with different combinations of wall thickness, stress, and neighboring point terms were introduced to predict future plaque growth based on previous time point data. Starting from the T2 plaque geometry, plaque progression was simulated by solving the solid model and adjusting wall thickness using plaque growth functions iteratively until T3 is reached. Numerically simulated plaque progression agreed very well with the target T3 plaque geometry with errors ranging from 11.56%, 6.39%, 8.24%, to 4.45%, given by the four growth functions. We believe this is the first time 3D plaque progression simulation based on multi-year patient-tracking data was reported. Serial MRI-based progression simulation adds time dimension to plaque vulnerability assessment and will improve prediction accuracy for potential plaque rupture risk.
Cardiovascular disease (CVD) is becoming the number one cause of death world-wide. There has been considerable effort investigating mechanisms governing atherosclerotic plaque progression and rupture [Friedman, Bargeron, Deters, Hutchins and Mark (1987); Friedman and Giddens (2005); Giddens, Zarins, Glagov, S. (1993); Ku, Giddens, Zarins and Glagov (1985); Gibson et al. (1993); Scotti et al. (2005); Yang, Tang, Atluri et al. (2008); Yuan, Mitsumori, Beach, and Maravilla (2001)]. A large number of the fatal clinical events are caused by rupture of a vulnerable atherosclerotic plaque [Fuster (1998); Fuster et al. (1990); Naghavi et al. (2003a, 2003b)]. Many victims of the disease who are apparently healthy die suddenly without prior symptoms. While major advancements in the treatment of CVD continue, progress has been very limited in early detection and treatment of at-risk individuals. Computational models and methods based on real patient data to simulate plaque progression and predict possible future rupture are lacking in the literature.
Most efforts for plaque progression research were focused on fluid dynamics side since it has been well accepted that atherosclerosis initiation and progression correlate positively with low and oscillating flow wall shear stresses. However, this “low and oscillating shear stress hypothesis” cannot explain why moderate and advanced plaques continue to grow under elevated flow shear stress conditions [Tang et al. (2005)]. Our recent results using serial MRI patient-tracking data and computational models indicated that 18 out of 21 patients studied showed significant negative correlation between plaque progression measured by wall thickness increase and plaque wall (structure) stress [Tang et al. (2008)]. Two-dimensional meshless generalized finite difference (GFD) computational models were developed using patient-specific plaque progression data to simulate plaque growth and predict future plaque rupture risk [Yang, Tang, Atluri et al. (2008)].
In this paper, we extend our previous 2D meshless GFD models to three-dimensional carotid plaque models with bifurcation. A computational procedure based on 3D meshless generalized finite difference method and serial MRI data was introduced to quantify patient-specific carotid atherosclerotic plaque growth functions and simulate plaque progression. 3D models are much closer to real physical plaques compared to 2D models, leading to more accurate and reliable predictions for plaque growth and mechanical stress conditions. By using plaque progression simulation, plaque vulnerability assessment and clinical decisions can be based on multi-time MRI scans and simulated “virtual” plaque progression. With validation, our procedure can be implemented in clinical applications. Simulated future plaque morphologies with stress/strain distributions can be displayed to physicians and patients and will lead to considerable improvement in prediction accuracy for potential plaque rupture risk and possible prevention of fatal cardiovascular events.
Because of the complexity of plaque geometry and structure, meshless GFD method has the advantage that it is meshless, thus avoids frequent re-mesh process that finite element methods require. Computational modeling for engineering applications with meshless methods have made considerable advances in recent years [Atluri (2004, 2005); Atluri, Yagawa, and Cruse (1995); Bathe (1996, 2002); Ling, Atluri (2006); Ahrem, Beckert and Wendland (2006); Shu, Ding, and Yeo (2005); Wu, Shen, Tao (2007)]. A series of meshless local Petrov-Galerkin (MLPG) methods were introduced to solved 3-dimensional elasto-static and dynamical problems [Han and Atluti (2004a, 2004b, 2007); Han, Liu, Rajendran, Atluri (2006)] and nonlinear problems with large deformation and rotations [Han, Rajendran and Atluri (2004)]. A “mixed” approach was introduced to improve the MLPG method using finite volume method [Atluri, Han and Rajendran (2004)] and finite difference method [Atluri, Liu, and Han, (2006a, 2006b); Hu, Young, Fan (2008)]. Sladek et al. applied the MLPG method to models of intelligent material systems where the mechanical property of the material being studied changes over time, and models of 3D anisotropic functionally graded solids [Sladek, Sladek, Solek, and Atluri, (2008); Sladek, Sladek, Tan, and Atluri, (2008)]. Li and Atluri further developed the MLPG mixed collocation method for material orientation and topology optimization [Li and Atluri, (2008)]. Cai and Zhu introduced a local meshless Shepard and least square interpolation method based on local weak form [Cai and Zhu, 2008]. A new meshless interpolation scheme for MLPG method was developed by Ma [Ma (2008)]. Young et al. developed hypersingular meshless method for solving 3D potential problems with arbitrary domain [Young et al. (2009)]. A finite volume meshless local Petrov-Galerkin method for topology optimization design of the continuum structures was introduced by Zheng et al. (2009). Analysis of structure with material interfaces was performed by Masuda and Noguchi [Masuda, Noguchi (2006)]. Perko and Sarler studied weight function shape parameter optimization in meshless methods for non-uniform grids [Perko and Sarler (2007)]. A meshless local natural neighbour interpolation method was applied to structural dynamic analysis by Liu et al. [Liu, Chen, Li, and Cen, 2008]. Remeshing and refining with moving finite elements were investigated by Wacher and Givoli [Wacher, Givoli (2006)]. Numerical methods were also developed to solve problems with free and moving boundaries [Zohouri, Pirooz, and Esmaeily (2005); Mai-Duy and Tran-Cong, (2004)]. Zhang et al. developed high-fidelity tetrahedral mesh generation techniques for models based on medical imaging data with fluid-structure interactions [Zhang et al., (2009)]. Regularized meshless method for irregular domain problems was investigated by Song and Chen [Song and Chen, (2009)]. Sellountos, Sequeira and Polyzos performed elastic transient analysis with MLPG(LBIE) method and local RBFs [Sellountos, Sequeira and Polyzos, (2009)], Lucas et al. proposed an automated approach for mesh adaptation to enhance accuracy of algorithms [Lucas, van Zuijlen and Bijl, (2009)]. A boundary element formulation for incremental nonlinear elastic deformation of compressible solids was proposed by Colli, Gei and Bigoni [Colli, Gei and Bigoni, (2009)]. GFD methods have been used in many engineering applications and in our previous papers where irregular geometries and free-moving boundaries are involved [Kleiber (1998); Liszka and Orkisz (1980); Tang, Chen, Yang, Kobayashi and Ku (2002); Tang, Yang, Kobayashi and Ku (2001)]. One advantage of using GFD is that generalized finite difference schemes can be derived for user-selected irregular grid points which can be freely adjusted to accommodate plaque deformation and growth. The MGFD method introduced in this paper uses grid points from the local support of each nodal point so that theoretical MLPG framework can be applied [Atluri (2004)].
Due to the complexity of the problem, we start from a linear 3D structure model for a human carotid plaque with bifurcation in this paper to get some insight for further full 3D fluid-structure interaction investigations. Patient-specific plaque progression data was acquired by magnetic resonance imaging (MRI) at multiple time points (scan interval set at about 18 months). Corresponding slices were matched using carotid bifurcation as the registration point. 3D plaque bifurcation models were constructed and solved by meshless GFD method to obtain the stress distributions in the plaque. The meshless GFD results were validated by ADINA (ADINA R & D, Watertown, MA), a commercial finite element software package that we have used to solve the plaque models in the past 12 years. Vessel wall thickness (WT) was used to measure plaque progression. 100 points on the lumen and their corresponding points on the out boundary for each matched slice (total matched slices: 8) were selected to define the inner and out boundaries of the plaque which will be adjusted in numerical simulations using plaque growth functions. Pointwise plaque growth functions were determined based on three time point stress and vessel wall thickness data on the selected points from 8 matched sliced. The growth functions were then used to simulate plaque progression starting from Time 2 plaque morphology. The simulated Time 3 plaque morphology was compared with actual Time 3 plaque geometry to validate our modeling method. Details are given below.
Serial MRI data from one patient was provided by the University of Washington (UW) group using protocols approved by the University of Washington Institutional Review Board with informed consent obtained. Scan time intervals were about 18 months, subject to scheduling variations. MRI scans were conducted on a GE SIGNA 1.5-T whole body scanner using an established protocol outlined in the papers by Yuan and Kerwin et al. [Kerwin, Hooker, Spilker, Vicini, Ferguson, Hatsukami, and Yuan (2003); Yuan and Kerwin (2004)]. Upon completion of a review, an extensive report was generated and segmented contour lines for different plaque components for each slice were sent to Tang’s group for model construction and further computational mechanical analysis. Details of the model construction process can be found from [Yang et al. (2007); Tang et al. (2008)]. Figure 1 gives the re-constructed 3D geometries of the plaque at three time points showing plaque progression. Figure 2 gives segmented contour plots of the plaque at three time points. 5 selected slices. These slices were used for model construction. The matched 8 slices were used to determine plaque growth functions. Plaque geometry at Time 3 was also used to validate simulated plaque progression.
Since there was insufficient data to quantify individual plaque component growth, the plaque was treated as a uniform material, which was assumed to be linear, isotropic, incompressible and homogeneous. The governing equations and corresponding initial and boundary conditions are given below [Fung (1993;1994)]:
where ρ is vessel material density, (u1, u2, u3) is the displacement vector corresponding to x,, and z directions, with z being the axial (longitudinal) direction. σ is stress tensor, ε is strain tensor, Pin is the specified lumen pressure, Γ is vessel inner boundary, u30(t), u3L(t) are the pre-stretch functions for the z-displacement at the two ends of the tube, f, j stands for derivative of f with respect to the jth variable. For the 3D linear model, the strain-stress relationship is:
The boundary conditions are:
where (n1, n2, n3) is the normal direction of the vessel surface. In our model, the inner boundary (lumen surface) is with (1, 2, 3) = (Pinn1, Pinn2, Pinn3), Pin is the specified inner pressure. The outer boundary was set as a free boundary, with (1, 2, 3) = 0.
Meshless generalized finite difference (GFD) methods have been used in many engineering applications where irregular geometries and free-moving boundaries are involved [Kleiber (1998); Liszka and Orkisz (1980)]. The advantage of GFD method is that generalized finite difference schemes can be derived using arbitrarily distributed points (see Fig. 4). With GFD, we will be able to use denser nodal point distributions where plaque has higher stress/strain concentration or critical morphological features. We will also be able to adjust, move, add or drop nodal points as needed. This leads to the desired flexibility in handling the complex geometry and plaque growth where not only the geometry changes, the total plaque area (volume if 3D) also changes.
The GFD concept and derivation of the generalized finite difference schemes are explained by the following example. Fig. 3(a) gives the carotid plaque with distributed nodal points. Fig. 3(b) gives 4 leader-nodes (Pi) with their neighbor nodes (Zj). One leader node was a boundary node. The other 3 were inner nodes. The neighbor nodes were selected using a sphere support with radius: R = s*max (dx, dy, dz). The support size control data s = 3.0 was determined numerically to reach best agreement with solutions obtained by ADINA. Fig. 3(c) shows a simplified 2D plot of a leader node with its neighboring points, illustrating the derivation process of GFD scheme.
For inner nodes, we use 2nd order GFD scheme. For boundary nodes, we use 1st order GFD scheme. For each leader node Pi with Ni neighbor nodes Zj, j = 1, …, Ni (see Fig. 3(c)), the 2nd order Taylor expansion of f (x, y, z) at Pi is given by:
where the subscript j = 1… Ni is the subscript for the neighbor nodes Zj of Pi. f is u1, u2, u3 in displacement formula, fj = f (Zj). Dropping the last term, we have:
where fj = f (Zj), . GFD schemes for all the first and second order derivatives are determined from (19) using function values at the Ni neighbor nodes and least-square approximation techniques.
The first-order Taylor expansion is given by:
The first order GFD schemes can be obtained similarly from the following system:
2nd order center difference scheme was used for the time derivative term:
The vector = (u1, u2, u3)T is the displacement solution at time step (n + 1).
Under the in vivo condition, the artery is axially stretched and pressurized, thus axial and circumferential shrinking was needed a priori to generate the starting plaque geometry for the computational model. The shrinkage in axial direction was 9% so that the vessel would regain its in vivo length with a 10% axial stretch. Circumferential shrinkage for lumen and outer wall was determined so that: 1) total mass volume was conserved; 2) plaque geometry after 10% axial stretch and pressurization had the best match with the original in vivo geometry. Steady pressure condition Pin = 100 mm Hg was used in the model. The Young’s modulus was set to E = 3.6e + 6 = 360 KPa based on our experimental data [Kobayashi, Tsunoda, Fukuzawa, Morikawa, Tang, Ku (2003); Tang et al. (2008); Tang, Yang, Zheng, Woodard, Saffitz, Petruccelli, Sicard and Yuan (2005)] and current literature [Fung (1993); Humphrey (2002)]. Fig. 4 compares the contours obtained from the GFD model with MRI contours. While some variation can be observed due to the stretch and pressurization, the net wall thickness errors (averaged over all 1200 data points) were less than 2%.
A commercial finite element software package ADINA (ADINA R & D, Inc., Watertown, MA) was used to validate our GFD model. ADINA has been validated by hundreds of realistic engineering and real life applications and is well accepted in the industry and research communities [Bathe (1996); Bathe (2002)]. We have been using ADINA in the past 12 years to construct and solve 2D/3D artery models which were validated by experimental measurements [Tang et al. (2005, 2008)]. Finite element ADINA models were constructed for all three time points T1, T2 and T3 following the same procedures using in [Tang et al. (2005)]. Comparisons of results from the two models are given in Section 3.
Plaque GFD models were constructed based on in vivo MRI data at Time 1 (T1), Time 2 (T2) and Time 3 (T3) and solved to get 3D plaque geometry (which matched in vivo MRI geometry with errors less than 2%) and stress/strain data. Slices from the three scans (T1, T2 and T3) were matched using the carotid bifurcation as the registration point (see Fig. 2). In this paper, vessel wall thickness (WT) was selected as the measure for plaque progression. For each matched slice, 100 evenly-spaced points from the lumen were selected and a piecewise equal-step method was introduced to calculate wall thickness (Fig. 5). Each slice was divided into several segments according to its geometry (4 for this plaque sample). For each segment, equal step is used for inner and outer boundaries respectively to choose equal number of nodal points. The corresponding points on the inner and out boundaries are paired and the distance between the paired points are defined as vessel wall thickness at the given lumen point (Fig. 5). This method is sufficient for the cases covered in this paper.
Using the plaque geometry and stress data at T1 and T2, four plaque growth functions (GF1, …, GF4) were introduced to predict “next time step” plaque geometry:
where ak(j) are coefficients of the growth functions to be determined using T3 data, f can be displacement variables or wall thickness, σ is the maximum principal stress at the nodal point, which could also be calculated using numerical steps in the simulation, Δt = time step in simulation, j = 1, 2, … 8 is the slice number, i is the index for the points on each slice. GF1 uses wall thickness at one nodal point xi only. GF2 added neighboring points xi−1 and xi + 1 to GF1. GF3 included both WT and stress. GF4 added neighboring points xi−1 and xi + 1 to GF3.
Nodal positions of the selected 1200 nodal point pairs (100 per lumen including CCA, ECA and ICA) on plaque inner and outer boundaries obtained from the GFD model at T3 were used to determined the coefficients ak(j) in GF1-GF4. Least square approximations were used to reach the best fit. R2 values of the fitting results are given in Table 1.
Starting from the plaque geometry at T2 and using the plaque growth functions determined in 2.5, the following procedure was used to simulate plaque progression and try to reach best agreement with plaque geometry at T3 obtained from GFD model (called the target T3 geometry from here on):
The neighboring points and stress terms in GF2–GF4 were handled in similar ways. The w-value in the simulation formula was determined numerically to have the best matching with plaque T3 geometry. Results obtained from the simulation code are presented in next section.
One patient data with three time point MRI scans was used to construct the GFD models and demonstrate the simulation process. Results from GFD models were compared with that from ADINA models for model and method validation. Simulation errors from the four selected growth functions were compared. Details are given below.
Figure 6 compares plaque maximum principal stress (Stress-P1) obtained from both GFD and ADINA models on slice S8 and on a bifurcation cut surface. GFD results agreed well with ADINA results (error < 3%).
Four growth functions were used to simulate plaque progression from T2 to T3 following the procedure described in 3.3. Fig. 7 gives the overlapping contour plots of the simulated results at t = T2 + (T3 − T2)/4, T2 + (T3 − T2)/2, T2 + 3*(T3 − T2)/4, and T3, showing good agreement between simulated and actual plaque progression. Fig. 8 gives the 3D stacked view of the simulated plaque geometries at 4 time steps.
Fig. 9 shows 3D stack views of simulated plaque progression using GF2, GF3, and GF4. It is clear that GF2 and GF4 (with neighboring points) provided better agreement with the target T3 geometry, compared to GF1 and GF3, respectively. Similarly, GF3 and GF4 using both wall thickness and stress data led to better match with target T3 geometry than that given by GF1 and GF2.
To quantitatively compare the errors associated the four growth functions, we define the absolute and relative errors as:
where 1200 = total number of nodal points selected from inner boundary, target WT_T3 is the plaque wall thickness at T3 obtained from GFD model with axial stretch and pressure applied to reach best match with MRI geometry. The errors for the simulated plaque geometries using GF1–GF4 are given by Table 2. GF4 gave the best match with the target T3 geometry, which is 4.45%, only about 1/3 of the error given by GF1.
It is well-known that atherosclerotic plaque progression is a multi-faceted process involving not only mechanical factors, but also plaque type, component size and location, inflammation and lumen surface condition, biological processes including many different cell activities, blood conditions such as cholesterol level, diabetes, changes caused by medication, diet, life styles, exercise, emotions, etc. To reproduce the realistic plaque progression process accurately, we need to know the real mechanisms governing the progression and real-time information about all the factors involved in the process. To be honest, little is known about real human plaque progression mechanisms based on patient-tracking data. And it is absolutely impossible to have complete real-time patient-tracking data for the duration of plaque progression to its maturity (up to 50 years). Investigations and findings from all the channels, modalities and disciplines could be integrated together to obtain better and more thorough understanding of the complicated atherosclerotic progression process. Our computational simulation approach takes the actual multi-year plaque progression data to simulate future growth, based on the assumptions that the trend that led to the current state would continue into the future. In other words, the governing mechanisms (whether we know them or not) would remain the same: the contributing factors would continue to contribute the same way as they did in the past, with the changes taken into consideration by the terms included in the plaque growth function.
Our results indicated that the growth function with neighboring points and structural stress included provided better matching with the target geometry, with an error of 4.45%, compared to an error of 11.56% given by a growth function without neighboring points and structural stress terms. It is reasonable to expect that growth functions including flow shear stress, blood chemical conditions, cholesterol level and other related factors may lead to better progression simulation results.
Human atherosclerotic plaques have complex geometries. Mesh generation will take considerable effort to build a finite element plaque model. Plaque progression simulation requires that plaque geometry be adjusted for each numerical time step. One major limitation of finite element models or other similar models which involve mesh generation is that re-meshing of the structure domain is extremely time consuming and hard to be automated. This makes mesh-based simulation code for progression simulation near impossible. By using the meshless GFD method, we can adjust computational nodes in the plaque anyway we want to grow the plaque as the growth function dictates. And this can be written into the computational code so that it is done automatically. The method and model developed in this paper can be applied to many other biological and engineering processes where the computational domain changes in the process being investigated.
Clearly the current model is very limited and serves as initial demonstrations of the GFD method and the potential significant contributions from the progression simulation model. The material model needs to be extended to 3D nonlinear models to reflect artery stiffening behaviors. Plaque components should be included in the model. Fluid-structure interactions should be added. And the growth functions need to be adjusted to include fluid shear stress, blood chemical conditions, and the proper terms that can represent the plaque growth trend. Better understanding of the biological and mechanical factors will help us to better formulate the growth function.
We believe that this is the first time that 3D human carotid atherosclerotic plaque progression was simulated based on patient-specific plaque morphology and pointwise plaque growth functions derived from multi-year MRI data. Our results indicated that our proposed progression simulation process was able to accurately predict future plaque morphology if the current progression trend was continued. The 3D meshless GFD method worked well for the progression model. The predicted progression by the growth function including neighboring points and structural stress terms was considerably more accurate than that given by other growth functions. The current 3D structure model can be extended to 3D model with fluid-structure interactions and more terms can be added to the growth functions for better predictions. More case studies are needed to validate our findings. Accurate plaque progression simulation adds the time dimension to plaque vulnerability assessment strategies and should improve our predicting accuracies.
This research was supported in part by NIH grant NIH/NIBIB, R01 EB004759, and in part by NSF grant DMS-0540684. Professor Chun Yang’s research was partially supported by National Sciences Foundation of China grant 10871028 and priority discipline of Beijing Normal University. Drs Chun Yuan, Gador Canton, Thomas S. Hatsukami, Vasily Yarnykh, Baocheng Chu, and Fei Liu from University of Washington provided the in vivo MRI data and their efforts are happily acknowledged.