Copyright c 2008 Tech Science Press CMES, vol.28, no.2, pp.95-107, 2008 Meshless Generalized Finite Difference Method and Human Carotid Atherosclerotic Plaque Progression Simulation Using Multi-Year MRI Patient-Tracking Data Chun Yang 1, Dalin Tang 2, Chun Yuan 3, William Kerwin 2,FeiLiu 3, Gador Canton 3 Thomas S. Hatsukami 3,4, Satya Atluri 5 Abstract: Atherosclerotic plaque rupture and progression have been the focus of intensive investigations in recent years. Plaque rupture is closely related to most severe cardiovascular syndromes such as heart attack and stroke. A computational procedure based on 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 (T 1, T 2,andT 3,atintervals 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, hyperelastic, isotropic and nearly incompressible. The linear elastic model was used. The 2D plaque model was discretized and solved using a meshless generalized finite difference (GFD) method. Starting from the T 2 plaque geometry, plaque progression was simulated by solving the solid model and adjusting wall thickness using plaque growth functions iteratively until T 3 is reached. Numerically simulated plaque progression agreed very well with actual 1 School of Math, Beijing Normal University, Beijing, China 2 Corresponding author, dtang@wpi.edu, Math Dept, WPI, Worcester, MA 01609 USA 3 Deparment of Radiology, University of Washington, Seattle, WA 98195 4 Division of Vascular Surgery, University of Washington, Seattle, WA 98195 5 Center for Aerospace Research & Education, University of California, Irvine, CA 92612, USA. plaque geometry at T 3 given by MRI data. We believe this is the first time 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. Keyword: meshless, generalized finite difference, artery, plaque progression, plaque rupture, atherosclerosis. 1 Introduction Cardiovascular disease (CVD) has remained the number one cause of death in America since 1902. 36% of 45 year olds and 80% of those 75 and older have CVD to some degree. 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. Our long-term goal is to develop non-invasive methods to assess plaque vulnerability and predict possible rupture before the fatal event actually happens. 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); Scotti et al. (2005); Yuan, Mitsumori, Beach,
96 Copyright c 2008 Tech Science Press CMES, vol.28, no.2, pp.95-107, 2008 and Maravilla (2001)]. Most efforts 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 (200-700 data points/patient) 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. However, computational models using patient-specific plaque progression data to simulate plaque growth and predict future plaque rupture risk are lacking in the current literature. In this paper, a computational procedure based on meshless generalized finite difference (MGFD) method and serial MRI data is introduced to quantify patient-specific carotid atherosclerotic plaque growth functions and simulate plaque progression. By adding time dimension into our play, plaque vulnerability assessment and clinical decisions can be based on multi-time MRI scans and simulated virtual plaque progression. Seeing is believing. With validation, our procedure can be implemented in clinical applications 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, a meshless GFD method is introduced in this paper to avoid frequent mesh updates and adjustments. Computational modeling for engineering applications with meshless methods have made considerable advances in recent years [Ahrem, Beckert and Wendland (2006); Atluri (2004, 2005); Atluri, Yagawa, and Cruse (1995); Bathe (1996, 2002); Ling, Atluri (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)]. A new meshless interpolation scheme for MLPG method was developed by Ma [Ma (2008)]. 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)]. 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)]. 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 userselected 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)]. Details are given in the following sections. 2 Models and methods Due to the complexity of the problem, we start from 2D models in this paper to get some insight for further full 3D investigations. Patientspecific plaque progression data was acquired by serial MRI (scanning patients multiple times with time span at about 18 months). Corresponding slices were matched using carotid bifurcation as the registration point. 2D models were constructed for selected slices and solved by mesh-
Finite Difference Method and Plaque Progression Simulation 97 less GFD method to obtain the stress distributions in the plaque. Point-wise plaque growth functions were determined based on three time point stress and vessel wall thickness data. 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. 2.1 In vivo Serial MRI Data Acquisition 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 shows 5 (selected from 12) MRI slices with 5 different weightings obtained from a human carotid plaque sample. Figure 2 gives the re-constructed 3D geometries of the plaque at three time points showing plaque progression. Figure 3 gives three-time point segmented contour plots of 5 selected slices. These slices were used for model construction, quantification of growth function, and validation of simulated plaque progressions. 2.2 The structure model Since there was insufficient data to quantify individual plaque component growth, the plaque was treated as a uniform material, which was assumed to be hyperelastic, isotropic, incompressible and homogeneous. The governing equations and cor- (a) MRI T1-Weighting (b) T2 (c) Proton Density (PD) (d) TOF (Time of Flight) (e) CET1 (Contrast Enhanced T1) Figure 1: Selected multi-weighting MRI slices (5 of 12 slices) of carotid plaque from a participating patient. Multi-weighting MRI techniques can better differentiate plaque components and provide more accurate plaque vulnerability assessments [Yuan and Kerwin (2004)]. responding initial and boundary conditions are given below [Fung (1994)]: ρu i,tt = σ ij, j, i, j = 1,2; sum over j, (1) ε ij =(u i, j +u j,i )/2, i, j = 1,2,3, (2) σ ij n j out_wall = 0, (3) σ ij n j Γ = p in (t) Γ, (4) u i t=0 = u i0, (5) u i,t t=0 = u i0 (6) where ρ is material density, u =(u 1, u 2 ) is the displacement vector, σ isstresstensor, ε is strain tensor, P in is the specified lumen pressure, Γ is vessel inner boundary, f, j stands for derivative of f with respect to the jth variable. Using a linear model, the strain-stress relationship
98 Copyright c 2008 Tech Science Press CMES, vol.28, no.2, pp.95-107, 2008 (a) T 1, Front (c) T 2, Front (e) T 3, Front (b) T 1, Back (d) T 2, Back (f) T 3, Back Figure 2: Re-constructed 3D geometry of a carotid plaque based on in vivo serial MRI data. Three time point data are shown. T 1, T 2 and T 3 refer to time points from here on, unless otherwise indicated. T 1 T 2 : 525 days; T 2 T 3 : 651 days. Red lumen; Cyan- outer wall; Yellow - necrotic core; Fresh red - hemorrhage in necrotic core); Purple - loose matrix; Dark blue calcification; Green - fibrous cap. (a) Time 1 (b) Time 2 (c) Time 3 1.15 cm Figure 3: Segmented contour plots using CAS- CADE showing plaque components and plaque progression. 5 slices were selected with the bifurcation serving as the registration point. Data from three time points are shown. Red lipid; blue calcification (Ca). is [Fung (1994)]: σ 11 ε 11 σ 22 = C ε 22 (7) σ 12 γ 12 C = E 1 μ 0 0 0 1 μ0 2 μ 0 1 0 0 0 C 1 C 2 0 = C 2 C 1 0 0 0 C 3 1 μ 0 2 (8) where E 0 is the Young s Modulus, μ 0 is the Poisson ratio, C 1, C 2, and C 3 are coefficients defined by (8) for convenience. Substituting (2) and (7)- (8) into (1), we have the displacement equations in scalar form: ρ 2 u 1 t 2 = C 2 u 1 2 u 1 1 x 2 +C 3 1 x 2 +(C 2 +C 3 ) 2 u 2, 2 x 1 x 2 (9) ρ 2 u 2 t 2 = C 2 u 2 2 u 2 3 x 2 +C 1 1 x 2 +(C 2 +C 3 ) 2 u 1, 2 x 1 x 2 (10) The boundary conditions are: ( u 1 u 2 u1 C 1 n 1 +C 2 n 1 +C 3 n 2 + u ) 2 = t 1, x 1 x 2 x 2 x 1 (11) ( u 1 u 2 u1 C 2 n 2 +C 1 n 2 +C 3 n 1 + u ) 2 = t 2, x 1 x 2 x 2 x 1 (12) where (n 1, n 2 ) is the normal direction of the boundary and (t 1,t 2 )=(P in n 1,P in n 2 ) is the fluid force applied at the inner boundary (lumen), and the outer boundary is treated as a free boundary, with (t 1,t 2 )=0. 2.3 The meshless GFD method The advantage of MGFD method is that generalized finite difference schemes can be derived using arbitrarily distributed points (see Fig. 4). With
Finite Difference Method and Plaque Progression Simulation 99 MGFD, 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 is explained by the following example. Fig. 4(a) gives a selected nodal point Pi, its round support and all surrounding points Z j, j = 1,...,Ni, which are used to derive the GFD scheme. To derive the second order GFD schemes for the first and second order derivatives of the unknown function f (x 1,x 2 ), we use the Taylor expansion of f at Pi. Omitting higher order terms, we have, f f j = f i +h j f x +k j 1 i x + h2 j 2 i 2 x 2 1 i +h j k j x 1 x + k2 j 2 i 2 +o(h 2 j +k2 j ) (13) i x 2 2 where the subscript i indicates the selected node Pi at which derivative GFD schemes are being derived, j = 1...Ni are the other nodes in the Ninode-star of Pi, f j is the function value at Zj, Zj are the neighboring points of Pi within the support, h j = x 1 j x 1i,k j = x 2 j x 2i. Ignore the higher order terms in (13), we have:. h 1 k 1 h 2 1 2 h 1 k 1 k 2 1 2.. h h Ni k 2 Ni kni Ni 2 h Ni k 2 Ni 2.. = f x 1 f x 2 x 2 1 x 1 x 2 x 2 2 f 1 f i. f Ni f i i (14) Rewrite (14) as A f = g, and use the least-squares method to get the formula for f : f =(A T WA) 1 A T W g (15) (a) A Nodal-Star with Round Support x 2 Z1 Z6 Pi Z7 Z10 Z11 Z12 Z14 Z15 x 1 (b) Nodal Points for Slice 4, Time 2. Z2 Z3 Z4 Z5 Z8 Z9 Z13 Z14 Z16 (c) Nodal Points for Slice 4, Time 2, with Components. Lipid Figure 4: Meshless GFD scheme derivation and nodal point distributions. (a) A schematic plot illustrating the derivation process of meshless GFD schemes with round support; (b)-(c) Nodal points distributions on Slice 4, Time 2, with and without plaque components. where W is the weight matrix with w 1 0 W =... 0 w Ni Let w j = B =(b kl ) 5 Ni =(A T WA) 1 A T W Ca 1 (h 2 j +k2 j ) j = 1,2,...Ni, (16)
100 Copyright c 2008 Tech Science Press CMES, vol.28, no.2, pp.95-107, 2008 we have, f x 1 f x 2 x 2 1 x 1 x 2 x 2 2 = i ni j=1 ni j=1 ni j=1 ni j=1 ni j=1 b 1 j ( f j f i ) b 2 j ( f j f i ) b 3 j ( f j f i ) b 4 j ( f j f i ) b 5 j ( f j f i ) (17) Substituting (17) into (9)-(10), we get the discrete approximation for inner nodes. For boundary nodes, we do a similar work using 1 st order GFDM method. That is: f f j = f i +h j f x +k j 1 i x +o(h j +k j ) (18) 2 i h 1 k 1.. h Ni k Ni ( ) f x 1 f = x 2 i f 1 f i. (19) f Ni f i b 1 j ( f j f i ) (20) b 2 j ( f j f i ) ( ) f x 1 f = x 2 i ni j=1 ni j=1 Substituting (20) into (11)-(12), we get the discrete approximation for boundary nodes. Second order scheme for boundary nodes were also tested. However, it did not give good result. One order lower schemes for boundary nodes are common in computational schemes. Time derivatives are discretized using the 2 nd order center difference scheme: t 2 = f n+1 2 f n + f n 1 Δt 2 +O(Δt 2 ), (21) Assembling all the discretized equations and boundary conditions, we get the final linear system: K u = f, (22) The vector ( ) u1 u = u 2 is the displacement solution at time step (n+1). Equation (22) is solved by a sparse solver under Matlab environment. 2.4 The shrink-pressurize process and MGFD model validation When constructing in vivo MRI-based plaque models, a shrink-pressurize process needs to be used to (a) shrink the original in vivo MRI geometry to get the numerical starting geometry with zero lumen pressure and then (b) pressurize the reduced starting geometry to recover the original in vivo geometry with specified lumen pressure. Using Slice 4 (S4) at Time 2 (T 2 )as an example, Fig. 5 shows its original in vivo MRI geometry, the numerical starting geometry with 16.75% inner boundary shrinkage and 6.5% outer-boundary shrinkage, and the pressurized geometry obtained by solving the GFD plaque model. The outer boundary was reduced less so that the conservation law of mass (area for 2D models) is enforced. Young s modulus was set at E 0 =2.3 10 6 dyn/cm 2, 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)]. Patient-specific pressure = 136 mmhg (by arm) was used as the lumen pressure. A commercial finite element software package ADINA (ADINA R & D, Inc., Watertown, MA) was used to validate our MGFD 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 10 years to construct and solve 2D/3D artery models which were validated by experimental measurements [Tang et al. (2005,2008)]. A finite element ADINA model was constructed using S4 geometry and following the same procedures using in [Tang et al. (2005)]. Figure 6 compares maximum principal stress (Stress-P 1 ) from MGFD and ADINA models and shows that results from both models had very good agreement (error < 3%).
Finite Difference Method and Plaque Progression Simulation 101 a) Original MRI Geometry (b) Numerical Starting Geometry c) In vivo MRI Compared with Starting Geometry (d) In vivo MRI Compared with Pressurized Geometry Figure 5: Comparison of (a) original in vivo MRI plaque geometry with (b) numerical starting geometry, inner boundary shrinkage, 16.75%, outerboundary shrinkage 6.5%, and (d) the pressurized geometry obtained by solving the GFD plaque model, Slice 4 at T 2 was used. 3.1 A piecewise equal-step method to define wall thickness Vessel wall thickness was selected as the measure for plaque progression. In our previous paper, the shortest distance method was used to determine vessel thickness, i.e., for a selected nodal point on the inner boundary (lumen), the shortest distance between that point and the out-boundary was defined at the vessel thickness at that lumen point. That led to uneven selection of nodal points from the out boundary as shown by Fig. 7(a). A piecewise equal-step method is introduced to fix the problem. The vessel is divided into several pieces (segments) according to its geometry (4 in the S4 case). For each piece, 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. 7). This method is sufficient for the cases covered in this paper. (a) Stress-P 1 by GFD (b) Stress-P 1 by ADINA * Min=0.00 KPa * (a) Shortest Distance Method (b) Piecewise Equal-Step Method Δ Max=105.2 KPa Δ Max=108.6 KPa Min Universal Scale Max Figure 6: Comparison of maximum principal stress (Stress-P 1 ) by MGFD and ADINA for model validation. Figure 7: Piecewise equal-step method for determination of vessel wall thickness. Slice 4 at T2 is used for demonstration. The vessel wall was divided into 4 segments. 25 points were equally distributed on each segment. 3 Plaque growth function and progression simulation With the MGFD model validated, we can use it to derive plaque growth function and simulate plaque progression. 3.2 Quantifying plaque growth function (PGF) using serial MRI data Plaque growth function (PGF) is a function we use to determine the vessel wall thickness or nodal point displacement for every numerical time step based on current and past plaque geometry and
102 Copyright c 2008 Tech Science Press CMES, vol.28, no.2, pp.95-107, 2008 mechanical conditions. The following assumptions were made when deriving the plaque growth functions: a) Plaque growth depends on local plaque morphological and mechanical conditions; b) For simplicity, a single-point correspondence approach is used, i.e., plaque growth at a given nodal point is determined by information from the past and current time steps at the same point. In reality, influence from neighboring points should already be factored in by the plaque growth at the chosen nodal point; c) The first order time derivatives of the vessel wall thickness and stress conditions (which are functions of time, for every nodal point under consideration) should be included in PGF; d) The information at most current time step is most directly related to plaque growth in the nearest future. Two different plaque growth function forms were used in our derivation and simulation process. Three-time-point function. Data from three time steps T 1, T 2 and T 3 are used to fit wall thickness or displacement at T 3. The fitting function is: f T3 fit(i)=a i f T2 (i)+b i df dt (i) T2 +c i σ T2 (i)+d i dσ dt (i) (23) T2 where f is the wall-thickness or displacement function, i is the nodal point numbering index (100 points were chosen from the inner boundary), σ is the Stress-P 1 function, a i, b i, c i,and d i are coefficients determined by the least squares method to fit T 3 wall thickness or displacement data, and the time derivatives are calculated by dg dt (i)= T2 w g T 3 (i) g T2 (i) +(1 w) g T 2 (i) g T1 (i) T 3 T 2 T 2 T 1 (24) where g can be either f or σ, w is a weight function to be adjusted for better agreement when the growth function is used in progression simulation. When the growth function (23) is used in the plaque progression simulation code to adjust the inner and outer boundary at each numerical time step, the time derivatives are evaluated with T3 replaced by the current time corresponding to the numerical step: dg dt (i)= tn w g t n (i) g T2 (i) +(1 w) g T 2 (i) g T1 (i) t n T 2 T 2 T 1 (25) where t n indicates the current time in the numerical simulation. Two-time-point function Data from two time steps T 1 and T 2 are used to fit wall thickness or displacement at T 2. The fitting function is the same as that given by (23) except that the derivatives are calculated using two time points by: dg dt (i)= g T 2 (i) g T1 (i). (26) T2 T 2 T 1 In the simulation code, (25) will be used to calculate time derivatives, the same way as it is done for the three-point growth function. So the main difference between the three-point formula and the two-point formula is that the coefficients a i, b i, c i, and d i are determined using different time (two or three) data points. The formulas for the two growth functions in the simulation code are almost the same. 3.3 Plaque progression simulation Using S4 at T 2 and the plaque growth function determined in 3.2, the following procedure is used to simulate plaque progression: Step 1. Start from the original in vivo MRI geometry (S4 at T 2 ), use proper shrinkage to get the zero-pressure numerical starting geometry; Step 2. Discretize the geometry using the meshless GFD method, solve the model to get
Finite Difference Method and Plaque Progression Simulation 103 plaque geometry and stress/strain distributions under specified pressure conditions; (a) S4 at T 2 under 136 mmhg (b) S4 at T 3 under 136 mmhg (c) Overlapping (a) and (b) Step 3. Use the growth function to determine the plaque geometry for next numerical time step by adjusting the nodal points on inner and outer boundaries; Step 4. Adjust internal nodal points as needed; (d) Simulated S4 at Day 220 (blue). (e) Simulated S4 at Day 400 (blue). (f) Simulated S4 at Day 651 (T3). Step 5. Solve the plaque model using the updated plaque geometry; Step 6. Repeat Steps 3-5 till numerical time reaches T 3 (MRI scan time). Results obtained from the simulation code are presented in next section. 4 Results Slice 4 was used as the sample to demonstrate the simulation process. The procedure described in 3.3 was followed. Fig. 8 gives the starting, target and three simulated geometries, showing good agreement between simulated and actual plaque progression. Define the absolute and relative errors as: Absolute Error = ( WT_num(i) WT_T 3(i) ), (27) 100 Figure 8: Simulated plaque growth has good agreement with MRI data (error = 4.52%). Progression code starting time: T 2 ; ending time: T 3. S4 was used as the sample slice. The three-point growth function was used. S4 geometry at T 3 determined from MRI data was used as the benchmark data for validation. Lumen pressure: 136 mmhg. T 2 T 3 time span was 651 days. the two-point growth function (26). The errors are: Absolute Error = 0.0329cm, (31) Relative Error = 11.4%, (32) It is clear that the three-point growth function gave more accurate predictions. Relative Error = [( WT_num(i) WT_T 3(i) )/(WT_T 3(i) ], 100 (28) (a) Simulated S4 at Day 220 (blue). (b) Simulated S4 at Day 440 (blue). (c) Simulated S4 at Day 651 (T 3 ). where 100=total number of nodal points selected from inner boundary, we have: Absolute Error = 0.0126cm, (29) Relative Error = 4.52%, (30) for the case simulated (S4, 3-point formula). Figure 9 gives the simulated plaque geometries at Days 220, 440, and 651 (ending time, T 3 )using Figure 9: Simulated plaque geometries at three numerical time steps showing the prediction accuracy (error = 11.4%) using the two-point growth function was not as good as that from the threepoint growth function. Simulations were also conducted for S2, S3 and
104 Copyright c 2008 Tech Science Press CMES, vol.28, no.2, pp.95-107, 2008 S5. Results from S3 were similar to that of S4. S2 and S5 did not give good results because the original MRI geometries differ too drastically from T 1 to T 2 and from T 2 to T 3 and no clear trend could be quantified. 5 Discussion 5.1 The computational approach 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, cell activities, blood conditions such as cholesterol level, diabetes, changes caused by medication such as statin, and other chemical conditions, inflammation and lumen surface condition. 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. In stead of trying to identify the individual factors which contribute to plaque growth, our computational simulation approach takes the actual multiyear 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 trend would continue; 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 growth function. Our results from the 4 slices simulated indicated that our model worked well when the trend from T 2 to T 3 was more or less similar to that from T 1 and T 2. And the model did not work very well when T 2 geometry seemed to be out of place between T 1 and T 3. Our model was applicable to 50% of the cases considered. 5.2 Model limitations Clearly the current model is very limited and serves as initial demonstrations of the MGFD method and the potential significant contributions from the progression simulation model. The model needs to be extended to 3D. And the growth functions need to be adjusted to include 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. Fluid forces and blood conditions (cholesterol, lipid lowering medication factor) can be included for better accuracy of predictions. 5.3 Progression simulation and plaque vulnerability assessment The long term goal of the current study is that patient-specific quantitative plaque growth functions can be determined based on multiple annual MRI scans and used to simulate plaque progression. Future plaque vulnerability can be assessed based on predicted future plaque morphologies by the progression simulation models. We are adding the time dimension into plaque assessment technology to improve the predicting power and accuracy. If our studies are successful, annual MRI scans would be recommended to patients who are in their early-to-middle stages of atherosclerosis. Simulated plaque growth would be generated for their early diagnosis and proper treatment for prevention of serious or even fatal clinical cardiovascular events. 6 Conclusion We believe that this is the first time that human carotid atherosclerotic plaque progression was simulated based on patient-specific plaque morphology and point-wise 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 meshless GFD method worked well for the progression model. The predicted progression by the three-time-point growth function was considerably more accurate than that given by the two-time-point (4.52% vs. 11.4%). The current 2D model can be extended to 3D model with more terms added to the growth functions for better predictions. More case studies are needed to validate our findings. Accurate plaque
Finite Difference Method and Plaque Progression Simulation 105 progression simulation adds the time dimension to plaque vulnerability assessment strategies and should improve our predicting accuracies. Acknowledgement: This research was supported in part by NIH grant NIH/NIBIB, R01 EB004759 as part of the NSF/NIH Collaborative Research in Computational Neuroscience Program, and in part by NSF grant DMS-0540684. Drs Vasily Yarnykh, Baocheng Chu, and Fei Liu contributed in the in vivo MRI data acquisition and image processing and their efforts are happily acknowledged. References Ahrem, R.; Beckert, A.; Wendland, H. (2006): A meshless spatial coupling scheme for largescale fluid-structure-interaction problems, CMES: Computer Modeling in Engineering & Sciences, 12(2):121-136. Atluri,S.N.(2005): Methods of Computer Modeling in Engineering & the Sciences-Part I, Tech Science Press, Forsyth, GA. Atluri,S.N.(2004): The Meshless Local-Petrov- Galerkin Method for Domain & BIE Discretizations, Tech Science Press, Forsyth, GA. Atluri, S. N.; Han, Z. D.; Rajendran, A. M. (2004): A new implementation of the meshless finite volume method, through the MLPG Mixed approach, CMES: Computer Modeling in Engineering & Sciences, 6 (6): 491-513. Atluri, S. N.; Liu, H. T.; Han, Z. D. (2006a): Meshless local Petrov-Galerkin (MLPG) mixed collocation method for elasticity problems, CMES: Computer Modeling in Engineering & Sciences, 14(3):141-152. Atluri,S.N.;Liu,H.T.;Han,Z.D.(2006b): Meshless local Petrov-Galerkin (MLPG) mixed finite difference method for solid mechanics, CMES: Computer Modeling in Engineering & Sciences, 15(1): 1-16. Atluri, S. N.; Yagawa, G.; Cruse, T. A. (Editors), (1995): Computational Mechanics 95, Vols. I & II, Springer-Verlag, New York. Bathe, K. J. (1996): Finite Element Procedures, Prentice Hall, Inc., New Jersey. Bathe, K. J. (Ed.) (2002): Theory and Modeling Guide, Vol I & II: ADINA and ADINA-F, ADINA R & D, Inc., Watertown, MA. Friedman, M. H.; Bargeron, C. B.; Deters, O. J.;Hutchins,G.M.;Mark,F.F.(1987): Correlation between wall shear and intimal thickness at a coronary artery branch, Atherosclerosis, 68:27-33. Friedman, M. H.; Giddens, D. P. (2005): Blood flow in major blood vessels - modeling and experiments, Annals of Biomedical Engineering, 33(12):1710 1713. Fung, Y. C. (1994): A First Course in Continuum Mechanics, Third Edition, Englewood Cliffs, New Jersey. Fung, Y. C. (1993): Biomechanics: Mechanical properties of Living Tissues, Springer-Verlag, New York, pp. 68. Fuster, V. (1998): The Vulnerable Atherosclerotic Plaque: Understanding, Identification, and Modification, Editor: Fuster, V. Co-Editors: Cornhill, J. F.; Dinsmore, R. E.; Fallon, J. T.; Insull, W.; Libby, P.; Nissen, S.; Rosenfeld, M. E.; Wagner, W.D. AHA Monograph series, Futura Publishing, Armonk NY. Fuster V.; Stein, B.; Ambrose, J. A.; Badimon, L.; Badimon, J. J.; Chesebro, J. H. (1990): Atherosclerotic plaque rupture and thrombosis, evolving concept. Circulation. 82 Supplement II:II-47 II-59, 1990. Giddens, D. P.; Zarins, C. K.; Glagov, S. (1993): The role of fluid mechanics in the localization and detection of atherosclerosis. Journal of Biomechanical Engineering, 115:588-594. Han, Z. D.; Atluri, S. N. (2004a): Meshless local Petrov-Galerkin (MLPG) approaches for solving 3D problems in elasto-statics, CMES: Computer Modeling in Engineering & Sciences, 6(2): 169-188. Han, Z.D.; Atluri S. N. (2004b): A meshless local Petrov-Galerkin (MLPG) approach for 3- dimensional elasto-dynamics, CMC: Computers, Materials & Continua, 1(2):129-140. Han, Z. D.; Atluri, S. N. (2007): A System-
106 Copyright c 2008 Tech Science Press CMES, vol.28, no.2, pp.95-107, 2008 atic Approach for the Development of Weakly Singular BIEs, CMES: Computer Modeling in Engineering & Sciences, 21(1):41-52. Han, Z. D.; Liu, H. T.; Rajendran, A. H.; Atluri,S.N.(2006): The Applications of Meshless Local Petrov-Galerkin (MLPG) Approaches in High-Speed Impact, Penetration and Perforation Problems, CMES: Computer Modeling in Engineering & Sciences, 14(2):119-128. Han, Z. D.; Rajendran, A. M.; Atluri, S. N. (2005): Meshless local Petrov-Galerkin (MLPG) approaches for solving nonlinear problems with large deformations and rotations, CMES: Computer Modeling in Engineering & Sciences, 10(1):1-12. Hu, S. P.; Young, D. L.; Fan, C. M. (2008): FDMFS for Diffusion Equation with Unsteady Forcing Function, CMES: Computer Modeling in Engineering & Sciences, 24(1):1-20. Humphrey, J. D. (2002): Cardiovascular Solid Mechanics, Springer-Verlag, New York. Kerwin, W.; Hooker, A.; Spilker, M.; Vicini, P.; Ferguson, M.; Hatsukami, T.; and Yuan, C. (2003): Quantitative magnetic resonance imaging analysis of neovasculature volume in carotid atherosclerotic plaque. Circulation. 107(6), 851-856. Kleiber, M. (1998): Handbook of Computational Solid Mechanics, Springer-Verlag, New York. Kobayashi, S.; Tsunoda, D.; Fukuzawa, Y.; Morikawa, H.; Tang, D.; Ku, D. N. (2003): Flow and compression in arterial models of stenosis with lipid core, Proceedings of 2003 ASME Summer Bioengineering Conference, Miami,FL, 497-498. Ku, D. N.; Giddens, D.P.; Zarins, C.K.; Glagov, S. (1985): Pulsatile flow and atherosclerosis in the human carotid bifurcation: positive correlation between plaque location and low and oscillating shear stress. Arteriosclerosis, 5:293-302. Ling,X.W.;Atluri,S.N.(2006): Stability analysis for inverse heat conduction problems, CMES: Computer Modeling in Engineering & Sciences, 13(3):219-228. Liszka, T.; Orkisz, J. (1980): The finite difference method at arbitrary irregular grids and its application in applied mechanics, Computers and Structures, Vol. 11, 83-95. Ma,Q.W.(2008): A New Meshless Interpolation Scheme for MLPG_R Method, CMES: Computer Modeling in Engineering & Sciences, 23(2):75-90. Mai-Duy, N.; Tran-Cong, T. (2004): Boundary integral-based domain decomposition technique for solution of Navier-Stokes equations, CMES: Computer Modeling in Engineering & Sciences, 6(1): 59-76. Masuda, S; Noguchi, H. (2006): Analysis of Structure with Material Interface by Meshfree Method, CMES: Computer Modeling in Engineering & Sciences, 11(3):131-144. Naghavi, M.; Libby, P.; Falk, E.; Casscells, S. W.; Litovsky, S.; Rumberger, J.; Badimon, J. J.; Stefanadis, C.; Moreno, P.; Pasterkamp, G.; Fayad, Z.; Stone, P. H.; Waxman, S.; Raggi P.; Madjid, M.; Zarrabi, A.; Burke, A.; Yuan, C.; Fitzgerald, P. J.; Siscovick, D. S.; de Korte, C. L.; Aikawa, M.; Juhani, K. E.; Airaksinen, K. E.; Assmann, G.; Becker, C. R.; Chesebro, J. H.; Farb A.; Galis, Z. S.; Jackson, C.; Jang, I. K.; Koenig, W.; Lodder, R. A.; March, K.; Demirovic, J.; Navab, M.; Priori, S. G.; Rekhter, M. D.; Bahr, R.; Grundy, S. M.; Mehran, R.; Colombo, A.; Boerwinkle, E.; Ballantyne, C.; Insull, W. Jr.; Schwartz, R. S.; Vogel, R.; Serruys, P. W.; Hansson, G. K.; Faxon, D. P.; Kaul, S.; Drexler, H.; Greenland, P.; Muller, J. E.; Virmani, R.; Ridker, P. M.; Zipes, D. P.; Shah, P. K.; Willerson, J. T. (2003a): From vulnerable plaque to vulnerable patient: a call for new definitions and risk assessment strategies: Part I. Circulation. 108(14):1664-72. Naghavi, M.; (co-authors are the same as above). (2003b): From vulnerable plaque to vulnerable patient: a call for new definitions and risk assessment strategies: Part II. Circulation, 108(15):1772-8. Perko, J.; Sarler, B. (2007): Weight Function Shape Parameter Optimization in Meshless Methods for Non-uniform Grids, CMES: Computer
Finite Difference Method and Plaque Progression Simulation 107 Modeling in Engineering & Sciences, 19(1):55-68. Scotti, C. M.; Shkolnik, A. D.; Muluk, S. C.; Finol, E. A. (2005): Fluid-structure interaction in abdominal aortic aneurysms: effects of asymmetry and wall thickness. Biomed Eng Online, vol. 4:64-70. Shu, C.; Ding, H.; Yeo, K. S. (2005): Computation of incompressible Navier-Stokes equations by local RBF-based differential quadrature method, CMES: Computer Modeling in Engineering & Sciences, 7(2):195-206. Tang, D.; Chen, X. K.; Yang, C.; Kobayashi, K. and Ku, D. N. (2002): A viscoelastic model and meshless GFD method for Blood Flow in collapsible Stenotic Arteries, Advances in Computational Engineering & Sciences, Chap. 11, Tech Science Press, International Conference on Computational Engineering and Sciences. Tang, D.; Yang, C.; Kobayashi, K. and Ku, D. N. (2001): A Generalized Finite Difference Method for 3-D Viscous Flow in Stenotic Tubes with Large Wall Deformation and Collapse, Applied Num. Math., 38: pp. 49-68. Tang, D.; Yang, C.; Mondal, S; Liu, F.; Canton,G.;Hatsukami,T.S.andYuan,C.(2008): A Negative Correlation between Human Carotid Atherosclerotic Plaque Progression and Plaque Wall Stress: In Vivo MRI-Based 2D/3D FSI Models, J. Biomechanics, 41(4) :727-736. Tang, D.; Yang, C.; Zheng, J.; Woodard, P. K.; Saffitz, J. E.; Petruccelli, J. D.; Sicard, G.A.; Yuan, C. (2005), Local maximal stress hypothesis and computational plaque vulnerability index for atherosclerotic plaque assessment, Annals of Biomedical Engineering, 33(12):1789-1801. Wacher, A.; Givoli, D. (2006): Remeshing and Refining with Moving Finite Elements. Application to Nonlinear Wave Problems, CMES: Computer Modeling in Engineering & Sciences, 15(3):147-164. Wu, X. H.; Shen, S. P.; Tao, W. Q. (2007): Meshless Local Petrov-Galerkin Collocation Method for Two-dimensional Heat Conduction Problems, CMES: Computer Modeling in Engineering & Sciences, 22(1):65-76. Yang, C.; Tang, D.; Yuan, C.; Hatsukami, T. S.; Zheng, J. and Woodard, P. K. (2007): In Vivo/Ex Vivo MRI-Based 3D Models with Fluid- Structure Interactions for Human Atherosclerotic Plaques Compared with Fluid/Wall-Only Models, CMES: Computer Modeling in Engineering and Sciences, 19(3):233-245. Yuan, C.; Kerwin, W. S. (2004): MRI of atherosclerosis. Journal of Magnetic Resonance Imaging. 19(6):710-719. Yuan, C.; Mitsumori, L. M.; Beach, K. W.; Maravilla, K. R. (2001): Special review: carotid atherosclerotic plaque: noninvasive MR characterization and identification of vulnerable lesions, Radiology, 221, 285-299. Zohouri S.; Pirooz, M. D.; and Esmaeily, A. (2004): Predicting wave run-up using full ALE finite element approach considering moving boundary, CMES: Computer Modeling in Engineering & Sciences, 7(1):107-118.