Comparison of fractional flow reserve value of patient-specific left anterior descending artery using 1D and 3D CFD analysis

Int J Adv Eng Sci Appl Math https://doi.org/10.1007/s12572-020-00257-7 IIT, Madras Comparison of fractional flow reserve value of patient-specific left anterior descending artery using 1D and 3D CFD analysis Supratim Saha1 • T. Purushotham1 • K. Arul Prakash1 Ó Indian Institute of Technology Madras 2020 Abstract Coronary heart disease is a primary source of mortality in India and around the globe. There are different types of disease state, but the primary focus of this study is on the disease where there is a narrowing down of the lumen area in the blood vessel. The left anterior descending artery is under investigation, which is a part of the left coronary artery, and it delivers the heart muscles with blood. The stenosed coronary artery impedes the blood flow and triggers heart failure. In a mild situation, medicine can be prescribed, but in severe cases, surgery is required. The cases which fall between mild and severe are a dilemma to the doctor for taking clinical decisions. The fractional flow reserve (FFR) tells us about the functional acuteness of stenosed coronary artery in this situation. The patient-specific left anterior descending artery of a human arterial tree is numerically investigated based on computational fluid dynamics approach for quantifying the functional acuteness of stenosis in terms of FFR. The CAD models of patient-specific geometries are generated from multislice computed tomographic scan data obtained from various Indian patients. The simulation of the 3D model is done using finite volume-based solver in OpenFOAM. Similarly, the 1D stenosed artery is generated using an analytical equation, and the solution is obtained using locally conservative Galerkin method. A comparative study of FFR between 3D and 1D models is carried out in this study, and a correlation is proposed among 3D FFR, 1D FFR and different geometric attributes. & K. Arul Prakash arulk@iitm.ac.in 1 Department of Applied Mechanics, Indian Institute of Technology Madras, Chennai, India Keywords Coronary heart disease
Coronary lesion
Computational fluid dynamics
Fractional flow reserve 1 Introduction Atherosclerotic heart disease is a condition that causes about 20 million loss of life worldwide. The number of deaths associated with CHD in India went up from 17% in 2001–2003 (26 per cent in adults) to 23% in 2010–2013 (32% in adults) [1]. Any quick and precise approach of diagnosing defects in the coronary artery will, therefore, boost its succession planning, including clinical therapy anywhere in the world. For very critical coronary occlusions, it is often straightforward to make angioplasty or more effective surgical repair processes. Medication is a pervasive way to treat the mild situation of coronary artery disease. Nevertheless, most of the CHD cases fall between very mild and serious cases. It creates a quandary situation for clinicians to determine the safest option in this scenario. Fractional flow reserve (FFR) is the coronary stenosis pressure ratio used as a diagnostic measure to assess the severity. Pd is the distal pressure of the stenotic lesion, and Pp is the proximal pressure of the stenotic lesion [2] FFR ¼ Pd : Pp ð1Þ Patient’s diagnostic technique for determining stenosis seriousness is intrusive. The process requires insertion of a wire probe in a patient’s blood vessel and sending it for examination to the stenosis site. Ultimately, pressure values are assessed across the stenosis to quantify the seriousness of the blockage. However, because of its an intrusive process, FFR measurement can be involved with clinical complications, such as blood vessel rupturing, and 123 Int J Adv Eng Sci Appl Math is not only unsuitable for monitoring medically treated stenoses, but also for cost implications [3]. A computed tomography angiogram (CTA) acquires 3D images of the beating heart and larger blood vessel of patients. This 3D picture is replicated, and computational calculations are carried out to predict the FFR value on the created 3D artery model. There is no risk associated with this treatment as it is non-intrusive. Although 3D modeling predicts FFR performance close to that achieved employing the intrusive wire insertion method, the time needed for computations is massive because it incorporates for the realistic stenosed configuration and also considers the integration of terminal vessels as a lumped model into the 3D computational domain, which is responsible for the overall computing expense. The numerical simulation is done to conduct complete arterial network with the stenosed condition in one of the blood vessels, and with reasonably good precision 1D simulation can model such accuracy at a low computational cost. The guiding equations which govern flow of blood through the artery were provided by Sherwin et al. [4]. A complete 1D numerical model of arterial flow was considered by taking into account of the connections between the systemic circulatory system and coronary flow developed by Mynard and Nithiarasu [5] and also incorporated a model of the heart in the arterial circulation for numerical simulation. Almost all the earlier research includes a 3D computational model to calculate FFR. However, 3D numerical study require a great deal of computational cost, while 1D simulations not only comparatively require less time and but also less cost. Iguchi et al. [6] clinically investigated the effect of lesion size on functional severeness for intermediate coronary stenosis and found that lesion size has clinical importance in intermediate-grade coronary stenosis. Different authors have clinically investigated the importance of lesion size on the severeness of the stenosed artery, but there has been no significant computational study reported. Due to the effect of the size of the lesion, there is a dearth of comprehension about intermediate level stenosis. A comprehensive study of the hemodynamics parameter is needed due to the concurrent effect of blockage percentage and the length of lesion. Several studies were conducted using either 1D simulations or 3D simulations for FFR computing, but there are very few comparative studies that take together 1D and 3D models. Etienne et al. [7] perform a comparative work of 1D model and 3D model for arterial blood flow for hemodynamic parameters. Due to the enormous simulation time for full-scale 3D model, a comparison between 1D model and 3D model for patientspecific cases is necessary to determine the functional seriousness of stenosis using the FFR values. For different geometric attributes of coronary stenosis, hemodynamic 123 parameters are evaluated in this study. The comparative study of 1D FFR and 3D FFR values in the left anterior descending arteries is also performed. Error investigation is being studied between the two models. A correlation between 1D and 3D result will be proposed in this study that may be used instead of carrying out 3D simulations. 2 Numerical methodology 2.1 1D computational details The artery is perceived to be a cylindrical-shaped vessel with an elastic boundary condition at the wall. The 1D equations for continuity and momentum conservation are given in Sherwin et al. [4] and are oA oðAuÞ þ ¼0 ot ox ou ou 1 op þu þ ot ox q ox ð2Þ f ¼0 qA ð3Þ where A considered to be the cross-sectional area, u is the average velocity across flow area, p is considered as the pressure within the artery, q & 1060 kg/m3 is the blood density and f is considered as the friction force per unit length. For the modeling of the friction term, a steady, laminar and Poiseuille flow is assumed. The closing of the system of equations is achieved by providing additional constraint to relate the pressure with cross-sectional area based on vessel elasticity, Poisson’s ratio, and thickness of the wall. The relationship is defined in Formaggia et al. [8] and Olufsen et al. [9]. The pressure in the artery and area dependence is used as: pffiffiffi pffiffiffiffiffi p ¼ pext þ b A A0 ð4Þ where pext is the transmural pressure, A0 is the cross section where the transmural pressure is considered zero (i.e., p ¼ pext ) and material properties of the vessel are b. The unknown parameters (A) and (u) are gleaned from the expression of forward (w1 ) characteristics and the backward (w2 ) characteristics as:   ðw1 w2 Þ2 q 2 A¼ ð5Þ b 1024 1 u ¼ ðw1 þ w2 Þ 2 ð6Þ The numerical simulation is done using the process (Mynard et al. [5]) developed by the locally conservative Galerkin (LCG). Each element is treated with its own boundaries as a sub-domain in this method. Int J Adv Eng Sci Appl Math 2.2 Modeling of the arterial tree The entire arterial network is modeled with both systemic and coronary circulation and is shown in Fig. 1. The arteries and their simulation properties were taken from Mynard et al. [5] and implemented in the current study. The left endocardial artery is modeled as stenotic one and all the other parts of the arterial tree are modeled to be normal. 2.3 Geometry of the stenosed segment The equation for a straight duct representing the stenosed region is given by Sonu et al. [10]. Some improvement in the equation is made by incorporating the tapering nature of the duct which is shown in Fig. 2 and the modified expression is   e Lseg x y¼ Lseg    x Lm L2s R S0 1 cos 2p ð7Þ 2 Ls to the maximum cross-sectional blockage location, Ls is the length of the stenosis region of the artery, Di is the diameter at inlet, Do is the exit diameter and the severity is percent blockage of the lumen area through which blood flows in the artery. 2.4 1D validation study The simulated (present) results obtained are compared with the numerical results of Low et al. [11] as shown in Figs. 3 and 4. The artery selected for validation study is the right carotid artery. For normal conditions and heart function, the pressure and flow waveform pattern are found to be almost similar. 2.5 Boundary conditions Using a sigmoid function, the pressure profile is generated as an inlet to the arterial tree as shown in Fig. 5. The arterial tree outlet is given a zero coefficient of reflection. where Lm S0 ¼ Ls  x  Lm 2 severity ; 100 e¼ Ls ; 2 ðD i 2 ð8Þ Do Þ ð9Þ Lseg is the arterial segment length, Rs is the radius of the stenosed zone, Lm is the distance from the start of the artery 2.6 Grid independent study The numerical investigation using three mesh sizes of 50 nodes, 100 nodes and 150 nodes is performed. For the three mesh sizes, the axial velocity magnitude along the arterial length (Left Endocardial artery) is compared as shown in Fig. 6. It is observed that the mesh sizes of both 100 and 150 yield similar results. Therefore, for the simulation, 100 nodes per artery are considered. 2.7 3D computational details Blood is perceived as a Newtonian fluid in this study, and it is assumed that its movement through the arterial vessel is laminar and incompressible. The mass and momentum conservation is given by, oui ¼0 oxi   oui o ui uj þ ¼ oxj ot ð10Þ 1 oq o2 ui þ# q oxi oxi oxj ð11Þ The above equations are in conservative form. A OpenFOAM solver based on finite volume method is used to solve the CFD problem being studied. The pimpleFoam solver of OpenFOAM is used with secondorder accurate scheme for time, gradient and divergence. Fig. 1 1D branched network of arteries having important blood vessels taken from Mynard et al. [5] 123 Int J Adv Eng Sci Appl Math Fig. 2 Geometry of the stenosed segment with label Fig. 5 The pressure profile which is provided as a inlet to arterial tree Fig. 3 Comparison of right carotid artery pressure waveforms with the results of Low et al. [11] Fig. 6 Variation of axial velocity along the arterial length for grid independence test of a 70% severity 2.8 3D computational domain Fig. 4 Comparison of right carotid artery flow waveforms with the results of Low et al. [11] 123 The flowchart for the generation of computational domain is shown in Fig. 7. The multislice computed tomographic Int J Adv Eng Sci Appl Math scan data collected of Indian patients are segmented using the snake algorithm and transformed into a CAD model. It is further meshed in an ICEM-CFD, which is a commercial software package where the aspect ratio of the cells is kept close to one. 2.9 3D validation study The flow through a stenosed artery was validated with Sonu et al. [10] under steady-state condition. Reynolds number (Re) based on characteristics length (D) and reference velocity (Uavg ) at inlet is defined as Re = qUavgD/l. The hydrodynamic properties considered in the simulation are density (q = 1000 kg/m3), viscosity (l = 1 Pa s), velocity (Uavg = .5 m/s) and diameter (D = 1 m). Based on the following relationships, velocity and diameter are nondimensionalized u y U¼ ð12Þ ; U¼ : Uavg D The dimensionless velocity here is U, the dimensional speed is u, the average inlet velocity is Uavg . The inlet diameter is D, the cross-stream and the dimensionless length is y and Y. At the inlet of the computational domain, a parabolic velocity profile is imposed. The wall is given noslip boundary condition along with pressure (p = 0) at the outlet for simulation. Velocity at different locations along the computational domain is computed, as shown in Fig. 8. It has been observed to be in good accordance with the findings of the literature. The numerical results of FFR are compared with the experimental results by Saha et al. [12] using the 3D numerical model used in the present study and showed a similar trend with experimental results. 2.10 Boundary condition The Womersley velocity variation is achieved from the patient-specific left anterior descending artery’s volume Fig. 8 Comparison of velocity results at different downstream distance to stenosis flow rate waveform using the 1D numerical simulation. Figure 9 indicates the patient-specific flow waveform obtained from 1D simulation. The profiles of velocity are determined from the volume flow rate waveform based on the FFT operation. The velocity profile of Womersley is applied at the inlet of the left anterior descending artery which is the computational domain and constant pressure (p = 0) condition for outlet of the computational domain. 2.11 Grid independence test The analysis of grid independence is carried out for all cases, but for only one case it is shown here. The error is measured to the finest mesh that the simulation considers. The simulation is conducted in the computational domain using 300,232 elements. Fig. 7 Flowchart for generation of the computational domain 123 Int J Adv Eng Sci Appl Math Fig. 10 Location where primitive variables are computed in the computational domain and distal location of the stenosis, respectively. For the computation of FFR, the pressure values at probe a and probe c will be required. 3 Results and discussion Fig. 9 Flow waveform for 1D patient-specific case 3.1 1D numerical results No. of elements FFR Error (%) 105,433 .949 1.05 200,543 .954 .89 300,232 .956 .74 450,232 .963 – 2.12 Modeling viscosity The different viscosity models are incorporated in 3D simulation for computation of FFR and comparison are done among them. Non-Newtonian Models like Carreau Yasuda and Power law are compared with the Newtonian model for computation of FFR. The model constants of Carreau Yesuda and Power law are taken from Shibeshi et al. [13]. As the shear rate in our patient-specific models are more than 100 s-1, so it is a fair good approximation of Newtonian fluid [14]. It also seen that with different non-Newtonian model, the FFR value marginally changes with different viscosity models as shown in Table 1. The flow waveform indicates the amount of blood flow supplied to one of the coronary arteries. The reduction of flow due to any disease state means that the supply of all essential components such as oxygen and vital nutrients to the artery’s neighbouring cells are impeded. This may lead to the failure of cardiac function. Thus, flow waveform gives us more understanding of the disease state condition of the cardiovascular system. The flow waveform for different geometric attributes like severity and lesion length of stenosis is investigated. Figure 11 shows that with higher severity or blockage percent, the flow through the artery is hampered and indicative of more severe disease state. Another geometric 2.13 Probe locations The location where the computation of primitive variables like pressure and velocity is done is shown in Fig. 10. The probe location at a, b and c indicates the proximal, throat Table 1 FFR comparison for different Non-Newtonian model with the Newtonian model Model Newtonian Carreau Yasuda Power law FFR .68 .72 .71 123 Fig. 11 Flow waveform at different severity for corresponding lesion length Int J Adv Eng Sci Appl Math Fig. 12 Pressure waveform at different lesion length for corresponding severity parameter which plays a vital role in the functional severity of stenosis is lesion length. The flow investigation through various stenosis models having different lesion length for the same percent blockage is done. Figures 12 and 13 clearly show that the values of pressure and flow are significantly different for various lesion lengths. It is observed that with the increment in the length of lesion, the pressure and flow are reducing for given severity. The possible explanation of these above results can be given with the help of the two characteristics which dominate in 1D pulse propagation study. The two characteristics are forward traveling characteristics (w1 ) and backward Fig. 13 Flow waveform at different lesion length for 70% severity traveling characteristics (w2 ). Whenever there is a blockage or obstruction in the flow path, the backward traveling component dominates due to wave reflection. The wave reflection increases with the increase in blockage percent of stenosis. Because of this, w2 will increase the backward traveling characteristics. As a result, the velocity value ðu / ðw1 þ w2 ÞÞ that depends on the wave characteristics will decrease with the increment of percent blockage. Thus, the blood flow in the left anterior descending artery reduces with the increment in severity. For the same blockage percent, the same trend is found with the increasing lesion length as shown in Fig. 13. These observations suggest that not only blockage percent is a significant hemodynamic parameter but also the length of the lesion plays an vital role in the intermediate-grade coronary artery. The pressure waveforms for 70 severity of stenosis are plotted in Fig. 14 for one heart cycle, and the magnitude of pressure at the probes is observed in order as pa [ pc [ pb . The velocity waveform for 70 severity of stenosis is plotted in Fig. 15 for one heart cycle, and the magnitude of velocity at the probes is observed in order as vb [ vc [ va . 3.2 3D computational results The comparison is made between a diseased state artery and an artery without occlusion at a particular instant of the cardiac cycle. The typical pressure contours obtained from the patient-specific CFD simulation are shown in Fig. 16. The diseased artery is reconstructed into normal one by recovering the reduced lumen area due to build-up of Fig. 14 Variation of pressure waveform at different probe location 123 Int J Adv Eng Sci Appl Math Fig. 17 Streamlines along with velocity magnitude of the patientspecific case, a normal artery, b stenosed artery Fig. 15 Variation of velocity at different probe location models, which are consistent with different literature findings as shown in Fig. 18. This finding is obvious as the pressure head loss occurs due to the separation of the flow after the lowest cross-sectional region where the flow of blood has to take place against the gradient of the adverse pressure. The error between the models is shown in Table 2. Various geometric parameters may contribute to error. Severity, lesion length and coronary tortuosity may play a significant role in error. The wall distensibility can also be responsible for the error but while comparing with the silicon model, it qualitatively followed the same trend [12]. It is also seen that artery having almost the same severity, the error differs due to the tortuous nature of the artery. Regression testing suggests a correlation exists based on Fig. 16 Pressure contour of the patient-specific case, a normal artery, b stenosed artery plaque. The change in pressure is significant in the stenosed artery near the throat location than in normal artery. The streamlines in diseased and normal artery are shown in Fig. 17 for the same patient-specific case. The corresponding velocity magnitude is used to color each streamline. The velocity is maximum at the throat location in the diseased artery. 3.3 Comparison between FFR3D and FFR1D results In this study, the comparison of three patient-specific cases is performed for 1D and 3D FFR value. It is observed that with increasing severity, the FFR value decreases for all 123 Fig. 18 Comparison of FFR1D and FFR3D for patient-specific cases Int J Adv Eng Sci Appl Math Table 2 Error between the two models for three patient-specific cases 4 Conclusion Data FFR3D FFR1D Error (%) Case 1 .956 .939 1.8 Case 2 .889 .84 5.5 Case 3 .67 .626 6.5 The blood flow within the coronary artery is influenced by the forward and backward traveling wave characteristics. The flow and pressure waveform depends not only on the severity of coronary lesion, but also on the span of the lesion. The size of the lesion seems to have a clinical significance on the functional seriousness of stenosed coronary artery. FFR comparison between 1D and 3D model is done for different patients, and error is quantified between the models. The comparison of FFR is also performed based on statistical analysis, and a correlation is proposed that can be used instead of the real 3D numerical computation. A good estimation of the FFR values can be achieved by Windkessel component which can be implemented at the terminal of vessels instead of the model which has resistance and it may boost the 1D model. The effort is made to a clinical relevant problem to hasten the decision taken by practitioners using fractional flow reserve (FFR) parameter using a 1D model employing locally conservative Galerkin method, and the study predicts FFR errors with full-scale 3D numerical simulation done in an open-source software using CT scan data from patients. This may be stepping stone for generating a more feasible and accurate correlation which will predict the value of FFR instead of carrying out the full-scale 3D simulation. Fig. 19 Confrontation between numerical simulation and correlation for calculating FFR computational results. The results obtained by means of numerical methods and correlation are in good concurrence. The correlation for FFR calculation can be used instead of 3D full-scale simulation FFR3D ¼ aFFR1D þ bB þ cL þ d where B is severity, L is size of lesion in mm, a = .81626, b = 1.2677 9 106, c = 4.5419 9 106, d = .2599. It is found that at 99 confidence, the value of R2 , for FFR is .9939 for the correlation model is sufficiently significant. Figure 19 shows the effectiveness of the FFR regression model, where the contrast between the expected values of the correlation model and the numerical results is shown. The errors have been determined between the computational results for simulation and those provided by correlation. The variation from the simulated values is within ± 2%. The correlation may not be very feasible, but it may be a stepping stone for developing a more reliable and accurate correlation which will incorporate more intricate details of the 3D model for using it clinically. Acknowledgements The authors appreciate exploratory research support (APM1718844RFERKARU) provided by IIT Madras and Global Challenges Research Fund (GCRF) Grant (RB1819APM003SWA-NKARU) partial funding from Swansea University, UK. The authors also appreciate the arrangement for Computed Tomographic scan results by Sri Ramachandra Medical Center Institute, Chennai. The authors want to thank Prof. P. Nithiarasu from Swansea University for his scientific deliberation. References 1. Gupta, R., Mohan, I., Narula, J.: Trends in coronary heart disease epidemiology in India. Ann. Glob. Health 82(2), 307–315 (2016) 2. Tonino, P.A., De Bruyne, B., Pijls, N.H., Siebert, U., Ikeno, F., vant Veer, M., Klauss, V., Manoharan, G., Engstrøm, T., Oldroyd, K.G., Ver Lee, P.N.: Fractional flow reserve versus angiography for guiding percutaneous coronary intervention. N. Engl. J. Med. 360(3), 21324 (2009) 3. Koo, B.K., Erglis, A., Doh, J.H., Daniels, D.V., Jegere, S., Kim, H.S., Dunning, A., DeFrance, T., Lansky, A., Leipsic, J., Min, J.K.: Diagnosis of ischemia-causing coronary stenoses by noninvasive fractional flow reserve computed from coronary computed tomographic angiograms: results from the prospective multicenter DISCOVER-FLOW (Diagnosis of Ischemia-Causing Stenoses Obtained Via Noninvasive Fractional Flow Reserve) study. J. Am. Coll. Cardiol. 58(19), 1989–1997 (2011) 4. Sherwin, S.J., Formaggia, L., Peiro, J., Franke, V.: Computational modelling of 1D blood flow with variable mechanical properties and its application to the simulation of wave propagation in the 123 Int J Adv Eng Sci Appl Math 5. 6. 7. 8. 9. human arterial system. Int. J. Numer. Methods Fluids 43(67), 673–700 (2003) Mynard, J.P., Nithiarasu, P.: A 1D arterial blood flow model incorporating ventricular pressure, aortic valve and regional coronary flow using the locally conservative Galerkin (LCG) method. Commun. Numer. Methods Eng. 24(5), 367–417 (2008) Iguchi, T., Hasegawa, T., Nishimura, S., Nakata, S., Kataoka, T., Ehara, S., Hanatani, A., Shimada, K., Yoshiyama, M.: Impact of lesion length on functional significance in intermediate coronary lesions. Clin. Cardiol. 36(3), 172–177 (2013) Boileau, E., Nithiarasu, P., Blanco, P.J., Müller, L.O., Fossan, F.E., Hellevik, L.R., Donders, W.P., Huberts, W., Willemet, M., Alastruey, J.: A benchmark study of numerical schemes for onedimensional arterial blood flow modelling. Int. J. Numer. Methods Biomed. Eng. 31(10), e02732 (2015) Formaggia, L., Nobile, F., Quarteroni, A., Veneziani, A.: Multiscale modelling of the circulatory system: a preliminary analysis. Comput. Vis. Sci. 2(2–3), 75–83 (1999) Olufsen, M.S., Peskin, C.S., Kim, W.Y., Pedersen, E.M., Nadim, A., Larsen, J.: Numerical simulation and experimental validation of blood flow in arteries with structured-tree outflow conditions. Ann. Biomed. Eng. 28(11), 1281–1299 (2000) 123 10. Varghese, S.S., Frankel, S.H., Fischer, P.F.: Direct numerical simulation of stenotic flows. Part 1. Steady flow. J. Fluid Mech. 582, 253–280 (2007) 11. Low, K., van Loon, R., Sazonov, I., Bevan, R.L., Nithiarasu, P.: An improved baseline model for a human arterial network to study the impact of aneurysms on pressure-flow waveforms. Int. J. Numer. Methods Biomed. Eng. 28(12), 1224–1246 (2012) 12. Saha, S., Purushotham, T., Prakash, K.A.: Numerical and experimental investigations of fractional flow reserve (FFR) in a stenosed coronary artery. In: E3S Web of Conferences, EDP Sciences, vol. 128, p. 02006 (2019) 13. Shibeshi, S.S., Collins, W.E.: The rheology of blood flow in a branched arterial system. Appl. Rheol. 15(6), 398–405 (2005) 14. Jhunjhunwala, P., Padole, P.M., Thombre, S.B., He, C., Wang, J., Huang, Y., Zhu, T., Miao, Y., Li, Z., Fan, L., Yao, J.: NonNewtonian blood flow in left coronary arteries with varying stenosis: a comparative study. Mol. Cell. Biomech. 1313(11), 1–261 (2016) Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.