A novel framework for fluid/structure interaction in rapid subject-specific simulations of blood flow in coronary artery bifurcations

Background/Aim. Practical difficulties, particularly long model development time, have limited the types and appli-cability of computational fluid dynamics simulations in numerical modeling of blood flow in serial manner. In these simulations, the most revealing flow parameters are the endothelial shear stress distribution and oscillatory shear index. The aim of this study was analyze their role in the di-agnosis of the occurrence and prognosis of plaque development in coronary artery bifurcations. Methods. We developed a novel modeling technique for rapid cardiovascular hemodynamic simulations taking into account interactions between fluid domain (blood) and solid domain (artery wall). Two numerical models that represent the observed subdomains of an arbitrary patient-specific coronary artery bifurcation were created using multi-slice computed tomography (MSCT) coronagraphy and ultrasound measurements of blood velocity. Coronary flow using an in-house finite element solver PAK-FS was solved. Results. Overall behavior of coronary artery bifurcation during one cardiac cycle is described by: velocity, pressure, endothelial shear stress, oscillatory shear index, stress in arterial wall and nodal displacements. The places where (a) endothelial shear stress is less than 1.5, and (b) oscillatory shear index is very small (close or equal to 0) are prone to plaque genesis. Conclusion. Finite element simulation of (cid:193) uid–structure interaction was used to investigate patient-specific (cid:193) ow dynamics and wall mechanics at coronary artery bifurcations. Simulation model revealed that lateral walls of the main branch and lateral walls distal to the carina are exposed to low endothelial shear stress which is a predilection site for development of atherosclerosis. This conclusion is confirmed by the low values (cid:3)(cid:3) of oscillatory shear index in those places.


Introduction
Bifurcation sites of human arteries are among the most frequent locations affected by atherosclerosis, being involved in up to 20% of percutaneous interventions. Several studies on the distribution of atherosclerotic plaques in human arterial systems have shown that atherosclerosis occurs predominantly at branching locations of the vascular tree where the arteries have relatively complex geometry that result in disturbed blood flow behavior [1][2] . In these regions, complex hemodynamic conditions prevail and local flow disturbances dictate the localization and progression of atheroma. Endothelial shear stress (ESS), as the main local flow-related factor, was investigated for the first time on computational fluid dynamics models of carotid and coronary bifurcation 3 and distal abdominal aortas 4 based on autopsy material. Those studies demonstrated that areas with low ESS correlate with atherosclerotic plaque development confirmed at autopsy 5 . In addition, several in vivo experiments on animal models confirmed proatherogenic effects of low ESS 6 . Recently human in vivo studies of arterial models derived from multislice computed tomography (MSCT) or intravascular ultrasound demonstrated the role of low ESS in the initiation and progression of atherosclerosis and arterial remodeling 7 . According to ESS spatial distribution plaques are more prevalent in low ESS areas, lateral walls of the main vessel and side branches, while they are less common in the flow divider or carina, which is characterized by high ESS 8 . Recent basic research studies have shed light on the precise pathways by which low ESS leads to atherosclerosis 5 .
Over the years, mathematical modeling has become complementary to experimental testing in predicting the biomechanical behavior and investigating clinical problems 9,10 . Finite element analysis allows us to repeat numerical experiments changing some parameters, thus analyzing the effect and influence of a single component within the observed phenomenon 11,12 . While multiphysics problems, particularly fluid-structure interaction (FSI) problems, are too complex to solve analytically, computational simulations have become a useful tool in studying blood flow in the cardiovascular system. Coupled problems require software to analyze the fluid domain (blood) and the solid domains (artery wall).
There are many reports in which the arterial wall is observed as rigid. Some papers report that wall deformation caused by interaction between the blood flow and artery wall changes the hemodynamic forces acting on the arterial wall compared to the rigid wall simulation 13,14 . Rigid wall simulations overvalue the magnitude and the distribution of ESS on the arterial wall 14,15 .
The aim of this study was to present a new computational framework for rapid modeling, considering a fluid structure interaction, and calculation of ESS and oscillatory shear index (OSI) of atherosclerotic plaque free coronary artery bifurcation obtained from MSCT, as well as to prove if specific bifurcation segments (lateral walls of main vessel, side branch lateral walls) correlate with low ESS and discuss the influence of this local hemodynamics conditions for potential plaque formation.

Definition of ESS and flow characteristics
ESS is the component of stress coplanar with an endothelial surface caused by the friction of the flowing blood and the arterial wall. ESS is expressed in units of force per unit area (1 N/m2 = 1 Pascal [Pa] =10 dyne/cm2) 16,17 .
In relatively straight segments of the artery, ESS is pulsatile and unidirectional and varies within a range of 1.5 to 3 Pa (physiological ESS) throughout the cardiac cycle and yields a positive time-average 18,19 . In contrast, in arteries with geometrically irregular regions, and disturbed laminar flow, pulsatile flow causes low and/or oscillatory ESS. Low ESS may be unidirectional at any point, but has a periodically fluctuating magnitude in a considerably low timeaverage (less than 1-1.5 Pa) [18][19][20] . Low ESS usually appears at the inner segments of curvatures and upstream of stenosis 16 . Oscillatory ESS is described by changes in both directions (bidirectional) and magnitude between systole and diastole, with a result of very low time-average, commonly close to 0 18,19 . Oscillatory ESS qualifies predominantly downstream of stenosis, at the bifurcation lateral walls, and in the area of branch points 3,19,[21][22] .
The OSI is a dimensionless index that describes the variation of ESS at the arterial wall during the cardiac cycle taking into account only the changes of ESS direction 23 . At the places where ESS have extreme changes the value of this index is 0.5. As opposed, at the places where ESS have no changes the value of this index is 0.
Numerical simulation was performed on one model of plaque-free coronary artery bifurcation based on MSCT coronary angiography.

Fluid-structure interaction governing equations
Fluid and structure occupy different subdomains, thus their system of equations can be set-up individually within their subdomains. The governing equations for the fluid domain are based on the Navier-Stokes equations and the equation of continuity 24 . These equations can be expressed in the matrix form as: where M is the constant mass matrix, C is the constant damping matrix, K is the stiffness matrix which may depend on displacement, s F is the vector of external forces , , U U U are the vectors of acceleration, velocity and displacement, respectively. A detailed explanation of variables in the equations 1-3 is available in [24][25][26] .
Among others, we calculate the distribution of OSI and ESS in the desired number of discrete time instants within the cardiac cycle. ESS was calculated based on the equation: (4) where ij are the components of endothelial shear stress, is the dynamic viscosity, t v is the tangential velocity and n is the direction of a unit vector normal to the wall at the moment t .
The OSI was calculated based on equation: where s t is the surface shear vector, defined as: and t is the shear vector, defined as: t n In the equation (7) is the stress tensor.

Simulation code
The two main approaches for the numerical simulation of fluid-structure interaction problems are: (a) monolithic approach (strong coupling), and (b) partitioned approach (weak coupling). The monolithic approach involves the development of new software, while the partitioned approach uses the existing software to solve new classes of problems. The software PAK-FS is developed by coupling existing computational solid dynamics (CSD) software PAK-S and computational fluid dynamics (CFD) software PAK-F the using partitioned approach 27 . The algorithm of developed code is shown in Figure 1. This algorithm partially separates the used solvers in space. As boundary conditions the CFD program uses the current geometry of the interface and the velocity of the interface nodes. This information is obtained from the CSD program for solving the solids. On the other hand, loads from the fluid domain obtained using the CFD program is transmitted to CSD program.

Model generation
The finite element method requires the physical domain in which the problem is discretized. To date, the methodologies to characterize local shear stress distribution and other hemodynamic parameters in circulation have been used for research purposes only, because no tools are available to provide this information in a routine manner. Artery bifurcations may have a very complicated configuration 28 . Lately, the trend and need is patient-specific modeling 29,30 . It is very challenging to generate quality hexahedral meshes for complicated structures 31 . The hexahedral mesh quality has a significant role in the simulation by the finite element method 32 . This is very important in CFD, where numerical errors become visible in the flow field. During generating procedures and developing software for the creation of the finite element model the authors have sought to satisfy all the listed constraints.
The original images, comprising MSCT image data of the human coronary artery bifurcation, were acquired from patients of Clinical Center of Serbia using a CT750HD scanner (Discovery, General Electric) 33 . Image segmentation and geometry identification of the blood vessels were made in the software Mimics ( Figure 2).
Structured mesh generation is the most appropriate method for discretizing domains. It produces a body fitted mesh, exactly describing domain's boundaries and hence boundary conditions can be accurately modeled. Application of structured mesh for complex arterial geometries is very difficult or infeasible. Hence, the entire computational model  is subdivided into the number of small blocks 34 . Blocks are generated around a volumetric model (Figure 3a). This topology of building blocks is used to generate the finite element mesh of real patient-specific bifurcation of coronary artery using structured mesh techniques (Figures 3b and 3c). Blocks are then linked together to produce finite element mesh describing arterial lumen. The described method is implemented in in-house rapid mesh generator STL2FEM. The two obtained computational models are generated in less than an hour.

Fig. 3 -Multiblock approach implemented in a rapid mesh generator stl2fem: (a) topology of building blocks, (b) computational model before mapping to the real model, and (c) computational model. Blocks [bold lines in (a)] are represented by 8-node hexahedrons, described by vertices and edges.
Quadratic 8-node isoparametric elements are generated using transfinite interpolation without smoothing 35 . The solid domain is modeled by linear 4-node quadrilateral (shell) elements. The solutions were obtained using 20,237 elements and 21,848 nodes for lumen while arterial wall contains 2,303 elements and 2,344 nodes ( Figure 4).

Initial and boundary conditions
According to Perktold et al. 36,37 , blood is approximated as an incompressible Newtonian fluid with the density of 1,050 kg/m 3 and a dynamic viscosity of 0.003675 Pas. Blood vessel walls are modeled with a linearly elastic material with the Poisson's ratio of 0.49, the wall density of 1100 kg/m 3 , and the uniform wall thickness of 0.8 mm. One cardiac cycle (0.89 s) is simulated in 20 steps with the time step size of 44.49 ms. The maximal diastolic flow velocity in the inlet is 23 cm/s. The velocity profile used for simulating pulsatile flow is shown in Figure 5. For fluid-structure interface, small displacements of the vessel wall are assumed, while the lumen mesh is fixed 12 .

Fig. 5 -Maximal velocity profile in one human's cardiac cycle. When viewed in cross-section plane, velocity is increased gradually. Velocity profile is parabolic: the value of velocity near to the arterial wall is equal to 0, and value at the center of the vessel's cross section has a maximum value. The obtained ultrasound measurement of the flow rate gives the amount of blood that elapses through the vessel's cross section in the time unit. The maximal velocity value is determined on the basis of equality of volumetric flows through vessel's cross section calculated based on the maximal and mean velocities.
The mesh discretization of the fluid and the structure on the interface are compatible. The nodes of the solid model in the inlet section of the artery are constrained to have no displacement in the direction of fluid moving. Velocities of nodes on the contact surfaces are equal. We calculate ESS, OSI index, velocity, pressure and equivalent stress in coronary artery wall throughout the cardiac cycle. Figure 6. The calculated values of ESS through the cardiac cycle were from the lowest 0.03 Pa to 3.3 Pa. Low ESS was present in specific bifurcation segments: main branch lateral wall, lateral wall distal to carina, and lateral carina wall. In this segments ESS values were in the range from 0.03 Pa (in lateral wall distal to carina at the beginning of the cycle) to 0.8 (at the main branch lateral wall). Consequently, for specific segments: ESS of main branch lateral wall 0.07-0.8 Pa, carina lateral wall was 0.04-0.37 Pa and lateral wall distal to carina 0.03-0.3) (Figure 7). In contrast, the highest ESS was at the main branch anterior wall. A more complete picture of hemodynamic conditions inside the lumen gives the OSI index ( Figure 8). In specific low ESS bifurcation segments OSI index is 0 or close to 0. Areas with high values (up to 0.5) of this index are characterized by changes in the ESS direction. On the other hand, the places where this index is very small (close or equal to 0) are disposed to plaque genesis.

Endothelial shear stress in some steps is shown in
The overall picture of hemodynamic conditions prevailing coronary artery complements streamlines and pressure field in the blood (Figure 9), and the equivalent stress in the arterial wall ( Figure 10).

Fig. 8 -Oscillatory shear index: segments with low endothelial share stress (ESS) correlate with oscillatory share index (OSI) index O or close to O. In the areas with the OSI index
value of 0 or close to 0 ESS is unchanged in the direction within the cardiac cycle.

Discussion
The main findings of our study are as follows: simulations of fluid-structure interaction can be used in serial manner to investigate patient-specific hemodynamics conditions at the coronary artery bifurcation; CFD coronary artery model reveals that the lateral walls of the main branch and lateral walls distal to the carina are exposed to low ESS which is a predilection site for development of atherosclerosis; thanks to a novel algorithm that we implement, hardware-and time-consuming procedures for the surface reconstruction and computed aided design (CAD) modeling are now completely skipped.
The benefits of applying this procedure in clinical practice are presented on the basis of the numerical results. CFD is currently being used as a means of enhancing our understanding fluid flow in arteries, while CSD is used to calculate associated vessel wall forces. In this paper we presented a computational framework for FSI simulation.
Making new solvers to simultaneously solve FSI is very difficult and time-consuming task. To successfully solve such interdisciplinary problems, CFD and CSD software that al-ready exists are coupled. The multiblock approach embedded in in-house software STL2FEM is used for mesh generation. It is a very effective tool for generating computational meshes for complex geometries such as artery bifurcations in biomechanics and engineering. This approach dramatically speeds up the process of generating patient-specific finite element bifurcation models with minimum errors in the nu-  merical solution. The methodology was applied to a vascular patient-specific model of coronary artery bifurcation including elastic wall modeling in simulations. Previous considerations indicate that developed computational frame gives useful inputs to cardiologists. Thanks to the obtained results, cardiologists are set in the role of decision makers. They have a clear view of blood flow through a coronary artery bifurcation, so they can suggest optimal treatment strategies. ESS is one of the most essential factors influencing endothelial structure and function, and it is the main local flow related factor responsible for coronary atherosclerotic plaque formation and progression 2, 5 . A CFD model of plaque free coronary bifurcation reveals that the lateral walls of the main branch, and the lateral walls distal to the carina are exposed to low ESS. Thanks to the clear picture of this parameter behavior during the cardiac cycle, it is possible to predict the positions that are suspectible to plaque with a high probability. Using this facts, formation and progression of atherosclerotic, a plaque can be predicted. This low ESS specific bifurcation segment has higher probability for plaque development. In contrast ESS in the carina region is higher and this area is less likely to develop plaque.
According to several recent studies, low ESS promote atherosclerosis with several mechanisms: through nitric oxide (NO)-dependent athero-protection, increased low-density lipoprotein cholesterol uptake, sustain increase in oxidative stress, inflammation and smooth muscle cell proliferation and migration 5,20 .
Better understanding a complex dynamic interplay between local flow conditions and the plaque formation and progression in bifurcations may express areas of angiographic interest 8 . Plaque distribution may provide a useful information for appropriate technique selection in treatment of complex bifurcation lesions. It gives the operator a chance consider to protection of the side branch with additional wire or to decide in advance between one and two stent techniques. It is also valuable knowledge for future stent designs. We presented a model of plaque-free coronary artery, but this methodology may be applicable in different clinical and angiographic circumstances. This analysis takes into account only kinematic and dynamic effects. Calcifications on the wall and other disorders are not included in the model because based on current knowledge it was not possible to predict the extent of calcification deposition in relationship with fluid behavior. However, the presence of calcification is possible to simulate in future studies by changing characteristics of the arterial wall at calcification place.
Examining ESS in a diseased coronary artery may enable identification of the initial stages of plaques that potentially evolve a high-risk lesion.
Thanks to rapid modeling, this methodology not only allows characterization of arterial plaque at a single point in time but, also allows prognostic insight into how plaques evolve over time 38 . If it were possible to recognize in vivo plaques that were expected to develop features of rupture prone lesions, specific interventions could be performed to avoid adverse cardiac event.

Conclusion
Rapid modeling of patient-specific artery bifurcation opens new avenues in the investigation of the role of ESS in the natural history of atherosclerosis and in future we hope it will be a powerful tool in everyday medical practice in the decision making process for every single patient.