Numerical simulations of flow patterns in the human left ventricle model with a novel dynamic mesh morphing approach based on radial basis function

We present a new numerical simulation framework for prediction of flow patterns in the human left ventricle model. In this study, a radial basis function (RBF) mesh morphing method is developed and applied within the finite-volume computational fluid dynamics (CFD) approach. The numerical simulations are designed to closely mimic details of recent tomographic particle image velocimetry (TomoPIV) experiments. The numerically simulated dynamic motions of the left ventricle and tri-leaflet biological mitral valve are emulated through the RBF morphing method. The arbitrary Lagrangian-Eulerian (ALE) based CFD is performed with the RBF-defined deforming wall boundaries. The results obtained show a good agreement with experiments, confirming the reliability and accuracy of the developed simulation framework.


Introduction
Cardiovascular disease remains a leading cause of deaths globally. The intracardiac flow analysis has been attracting a lot of interest in the field of cardiology in the last years. Among various disease conditions of the human cardiovascular system, an early diagnosis of congestive heart disease plays a crucial role in treating this potentially fatal condition. The blood flow patterns in the left ventricle can be used as one of the clinical indicators of heart dysfunction. The focus of many investigations was the formation of the vortices, which play an essential role in fluid dynamics.
The process of the vortex ring formation in a model of the left ventricle was presented in the pioneering work of [1]. This is followed by a series of studies which investigated momentum and energy transfer of blood flow in simplified ventricle models, [2][3][4][5]. More recently, in extensive state-of-the-art reviews of [6][7][8][9], it is postulated that dynamics of vortices can immediately reflect important physiological changes of the left ventricle, and as such, could provide early signs of heart disease. The remaining research question is to reveal the exact mechanisms behind the pathophysiological implications of vortices in the left ventricle.
To investigate blood flow patterns and surrounding tissue dynamics, various imaging techniques are applied in clinical settings such as Doppler ultrasound echocardiography [6,8,[10][11][12][13] and magnetic resonance imaging (MRI) [9,[14][15][16][17]. There is a continuing need for further improvements of these techniques -especially to increase spatial and/or temporal resolutions, which will bring even more details on early changes of the characteristic blood flow patterns in the left ventricle.
Besides the experimental techniques, computer simulation methods such as computational fluid dynamics (CFD) have been also applied to study the idealized, subject-and patient-specific blood flow in the left ventricle. The advantages of CFD are in providing a direct non-invasive local variation of the pressure, improved accuracy in estimating the wall shear stress (WSS) along the left ventricle wall, and a significantly increased spatial-and temporal-resolution which can reveal detailed insights into blood flow patterns. As stated in Ref. [18], the leading challenge in this field is to have a credible simulation approach which can capture the transient flow features in a geometrically complex deforming domain.
Starting from the pioneering work of [19], a wide variety of computational studies addressing the blood flow in the left ventricle, have been performed. Idealized geometries of the left ventricle with a prescribed wall motion were studied in Ref. [20,21], and vortex structures generated were analyzed. The fully coupled flow/wall structure interactions (FSI) simulation approaches include the surrounding tissue dynamics, [22][23][24][25]. However, the complex composition of the myocardial muscles makes this approach particularly challenging to accurately reproduce the dynamic movement of the ventricle. In addition to the movement of the left ventricle wall, the dynamic behavior of the mitral and aortic valves impose other challenges for the numerical simulations, [26][27][28][29][30]. Along with the FSI simulations, image-based methods with a prescribed ventricle motion have been reported in several studies. Among them, two distinct approaches in modeling the LV moving boundary can be identified: the immersed boundary method (IBM) [31][32][33], and the arbitrary Lagrangian-Eulerian (ALE) method, [34][35][36][37][38]. Each of the proposed approaches has certain advantages, e.g., the IBM in using efficient numerical solvers and the ALE in a more accurate treatment of the wall-boundary conditions, respectively.
In the present study, we have developed a computational workflow that translates the clinical/medical (in vivo) or experimental (in vitro) images to a series of deforming wall boundary conditions of the left ventricle for the numerical simulations. The dynamic behavior of the left ventricle and the biological valve is based on the novel radial basis function (RBF) method for the morphing of the deforming numerical mesh used for CFD computations. There is a constant increase in applications of the morphing mesh techniques in automotive [39,40] and aerospace [41][42][43] industries. In contrast to these industrial applications, some initial attempts of the radial basis function (RBF) based CFD in biomedical applications are shown in Ref. [44,45], which simulated blood flow in the patient-specific carotid aorta and ascending thoracic aorta aneurysm, respectively. To the authors best knowledge, we present here the first generation of the fully integrated left ventricle simulation based on the RBF morphing approach.
Our choice is based on the great flexibility of the radial basis function (RBF) to interpolate scattered data, [46,47]. By defining so-called control nodes at the surface of geometry, even complex movements of the surface can be easily mimicked. Similarly, the RBF can be used to move the nodes of a fully unstructured numerical mesh without changing the internal addressing of the surrounding nodes, leading to a numerically efficient simulation (no re-meshing needed). The more detailed mathematical rationale of the RBF method used is provided in Appendix (B).

Development of the image based dynamic mesh morphing approach
We developed a novel image-based dynamic mesh morphing approach, which provides an efficient numerical simulation framework for predictions of blood flow patterns in the patient-or subject-specific left ventricle. The entire flow-chart of the developed procedure is given in Fig. 1. The basic idea behind the developed framework is to utilize clinical-and/or laboratory-obtained 4D image data (4D = 3D spatial coordinates + time instants) as the input to define dynamic boundary conditions of a moving left ventricle wall, as well as corresponding mitral and aortic valves, to finally provide detailed insights into instantaneous blood flow patterns and local pressure distribution. The entire workflow can be summarized as follows: 1. The clinical/experimental data are pre-processed first (this includes repairing, smoothing, cropping, etc., of particular geometrical segments). 2. The finite number of the control points at the surface of the geometrical segments, at particular time instants, is introduced. Typically, the number of particular time instants is not sufficient to have an adequate smooth variation of the geometry. To achieve this, the cubic spline interpolation of the control points from available time instants (frames) is used to achieve the desired temporal variation. 3. By implementing the radial basis function (RBF) method (i.e. geometries are mapped by control points on each frame), old surfaces are replaced by the newly morphed ones. The morphed surfaces are then calibrated by comparing some of the most important global parameters defined in clinical or laboratory settings. These parameters include the total flow rate variation during the pulsating cycle, the opening area of the valves, characteristic angle of valve leaflets, etc. 4. Next, the volume mesh is created from the morphed surfaces at the reference time instant (e.g. at the peak systole). Performing the RBF mapping, the isomorphic volume meshes of each particular time instant are generated. 5. Initiating the CFD simulation from the basic mesh and specified initial boundary conditions, the surface node speed can be derived at each time step (i.e. v b = ds/dt, where s is the node displacement). The interior volume mesh nodes are smoothed regarding the predefined skewness of the control volumes. Additional local mesh refinement is performed if this predefined skewness criterion cannot be fulfilled. 6. After completion of the CFD simulations, the postprocessing of the results obtained is performed and compared with available experimental results.

Geometry reconstruction and registration
To be able to numerically mimic as closely as possible the experimental conditions, special attention is devoted to geometric reconstruction of the most important elements: (i) the left ventricle, (ii) the biological valve, (iii) the casing containing simplified left atrium and aortic root. The available information contains the CAD model of the casing, the video fragments of a top view of the valve movement recorded during the experiments, and finally, the shape of the left ventricle resulting from a combination of the CT-derived statistically shape model and reconstructed tomographic PIV images. Next, we will provide details of this process for each of the essential geometrical segments.

Left ventricle
The geometry reconstruction of the left ventricle starts from the available CT scans of 150 patients. From these scans, by applying the statistical shape modeling (SSM) approach, a representative average geometry of the left ventricle is generated for 20 characteristic time instants of the cardiac cycle [49], as shown in Fig. 2 (a). This statistically  averaged shape of the left ventricle was used to construct the optically accessible geometrical model, which was then experimentally investigated by combined particle image velocimetry (PIV) [37,48] and ultrasound (US) [12,50] techniques, respectively.
In the present study, we focus on mimicking experimental conditions of our recent tomographic PIV measurements, [48]. The statistically averaged geometry of the left ventricle (at the peak systole) is shown in Fig. 2(b). The 3D rapid prototyping method was used to print the basic mold of the left ventricle on which a silicon membrane with a thickness of 1 mm was made. Note that this silicone membrane also includes a part of the casing, i.e. it also covers parts of the simplified left atrium and aortic root. The dynamics of the left ventricle membrane is experimentally imposed through the pressure variation of the surrounding fluid, which in turn generates corresponding compression and expansion of the basic shape of the left ventricle.
An example of the instantaneous TomoPIV measurements is given in Fig. 3, where characteristic flow patterns ( Fig. 3(b),(c)-left) and shape of the left ventricle ( Fig. 3(b),(c)-right) in the central vertical plane at two characteristic time-instants (at the beginning of diastole, t/t 0 = 0.38, Fig. 3(b), and at the end of diastole, t/t 0 = 0.82, Fig. 3(c), as indicated in Fig. 3(a)) are shown.
From the time-resolved TomoPIV data, we extracted the left ventricle surfaces at 30 equally distributed time-instants of the entire cardiaccycle, Fig. 4 (Data from PIV). Due to limitations of the experimental technique, we observe some unphysical surface deformations containing randomly distributed holes and bumps, Fig. 4 (Smoothing and repairing). The very first step in our geometrical reconstruction of the left ventricle shape is to filter these imperfections and to generate a smooth surface. This geometry repairing is done by applying the smoothing algorithm of [51]. In order to morph the CT based model to the repaired TomoPIV model, approximately 400 control nodes are defined on both geometries, Fig. 4 (Define control nodes). This number of control nodes was sufficient to properly capture the typical local surface curvature and to mimic the spatial resolution of the TomoPIV. Then, by applying the Fig. 4. Illustration of the entire process of the geometry conversion from the TomoPIV measurements to the morphed STL format geometry used for numerical simulations.  RBF morphing approach, in total 30 isomorphic geometries are generated and registered for the particular instants of the cardiac-cycle, Fig. 4 (Morphed STL). From these 30 isomorphic geometries, by applying the cubic spline interpolation, particular geometries for each time-instant can be easily generated. It should be noted that in addition to the geometry registration of the left ventricle volume (from the reconstructed surrounding surface), the characteristics flow rates and corresponding inlet and outlet boundary conditions must also match the TomoPIV experiments. In the current numerical simulations, the numerical mesh is generated only once (matching the peak systole phase from TomoPIV measurements). The coordinates of the surface nodes are extracted from the entire volume mesh. Then, by applying the RBF morphing method, coordinates of every node at each time instant can be easily calculated.

Biological valve
In contrast to the above described left ventricle geometry reconstruction based on the combined CT/TomoPIV imaging, due to limited view of the used PIV measurements (i.e. because of the non-transparent material of the biological valve), the geometry of the biological valve needs to be generated differently. The three-dimensional geometry of the biological valve, which is used in the previous experiments [12,48], is shown in Fig. 5(a) (a side view) and (b) (a top view). From the 3D CAD drawing, geometries for a fully open (the orange leaflet) and fully closed (red leaflet) valve are constructed first, as indicated in Fig. 5(c). To mimic the characteristic opening extension during the pulsating cycle, the horizontal cross-sections are recorded from an additional camera installed for a top-view, as shown in Fig. 6, where white lines indicate extent of the opening. These video recordings (taken with the time-resolution of 120 frames per second) revealed quite important non-linear geometric deformations of the leaflets of the biological valve. Based on these recordings and calculation of the effective opening area (the opening area as calculated from the contrast-image postprocessing MATLAB module), jointly with synchronization and registration of the measured flow rates, characteristic key-frames of the 3D valve geometry were identified, Fig. 7. In addition to the already mentioned fully open and fully closed valve state, an intermediate key frame obtained by these top-view camera recordings is shown in Fig. 5(c) (green leaflet). Then, the control nodes at the valve surface were generated and the identical RBF morphing method (as used for left ventricle geometry) was applied at these characteristic key-frames. To generate valve geometries of all intermediate time-instants, a cubic interpolation between the registered key-frames was used (typically about 200 intermediate geometries are generated).

The final assembly
In the final step, the complete geometry containing the left ventricle, the biological valve, and fixed casings, is assembled. Note that the casings are fixed as well as the positioning of the biological valve within it. There it is of crucial importance to have fully synchronized registrations of the left ventricle geometries and characteristic positions of the leaflets of the biological valve at characteristic time-instants of the heart cycle. The above mentioned RBF morphing method is implemented in the PyGeM source code package (Python Geometrical Morphing) [52]. The detailed mathematics of the RBF method is provided in detail in Appendix (A). The interface between the PyGeM results and the CFD package ANSYS/Fluent is established through in-house-developed user-defined functions (UDF). Based on the package the geometric morphing of the dynamic left ventricle as well as the biological valve is achieved. The complete assembly is shown in Fig. 2 (c).

Governing equations
In our computer simulations, we apply the ALE (Arbitrary Lagrangian-Eulerian) formulation for conservation of mass and momentum. The conservation of mass (continuity equation) can be written as: where v represents velocity vector, v b is the velocity vector on the surface boundary of the model including the valve, n is the normal vector on the surface boundary, ρ is the density of the working fluid. The momentum equation is:   (2) where p is the pressure, I is the unit tensor, τ = μ(∇v +∇v T ) is the tensor of viscous stress.

Simulated experimental working conditions
In the present simulation, the fluid properties (density (ρ) and dynamic viscosity (μ)) are assumed to be constant and consistent with the properties of the working fluid (i.e. a mixture of the 60% glycerine and 40% water to assure the transparency and a good refractive index matching with the silicone membrane) as in the reference experiments of [48]. All relevant fluid properties and flow parameters of the experiment which is numerically simulated are given in Table 1.

Numerical details
The discretised equations are solved by the finite-volume based ANSYS Fluent (version 19.1) solver. A summary of the numerical schemes used is provided in Table 2. The simulation starts from the beginning of the diastole with zero-velocity as initial condition. The entire fluid flow is generated by the prescribed motions of the left ventricle surface and the valve. The outlet is closed during the diastole and then opened during the systole. The entire simulation domain is meshed by non-uniformly distributed tetrahedral control volumes (CVs). The typical number of total control volumes varied between approximately 4 × 10 6 and 12 × 10 6 CVs. To ensure a good mesh quality and to reduce the cell skewness, the extensively deformed volume cells were locally smoothed or additionally locally re-meshed. The local refinement and smoothing of numerical mesh was mostly needed during the valve opening phase.

Results and discussion
The simulations results showing characteristic transient ventricular flow patterns are plotted in Fig. 8. We simulated in total 10 cardiac cycles to eliminate possible effects of the initial transients in the flow field. The results shown are for the last cardiac cycle. We select three time instants, which include fully closed (a), start of the opening (b), and finally, fully open (c) biological valve states -identical to those shown in Fig. 7. The streamlines colored by the velocity magnitude in the central vertical plane are superimposed with iso-surface of the velocity magnitude (|V| = 0.3 m/s) for all time-instants. The first time-instant is for a fully closed mitral and fully open aortic plane, Fig. 7 (a). The flow is oriented towards the ejection plane and two low intensity recirculation regions are generated in the proximity of the closed valve. The initial phase of the mitral valve opening is shown in Fig. 7 (b). The formation of the central jet is clearly visible in the proximity of the valve, whereas the complex multi-roll patterns are generated in the ventricle interior. The following time-instant is at the peak of the generated flow-rate, with a fully open tri-leaflet mitral valve, Fig. 7(c). The distinct central jet can be observed followed by generation of several local circulation regions. The simulated flow patterns are in very good qualitative agreement with experimental observation described in Ref. [48] or similar studies in the literature, [3,7,14].
Next, we perform a quantitative direct comparison with the Tomo-PIV measurements by selecting characteristic profiles at various instants of the cardiac cycle, Fig. 9. At the early diastole, early formation of the central jet can be seen, Fig. 9(b-c). Both, the angle and the penetrative length of the jet, show a good agreement between simulations and experiments. This is confirmed by extracting the horizontal profiles across the left ventricle, as shown in Fig. 9(a). A good agreement is obtained across the entire profile, with proper capturing of the central peak, and the outlet value at the aortic plane. This good agreement is also confirmed by calculating the Pearson correlation coefficient between CFD and experimental velocity profiles of r 1D = 0.96. For the mid diastole, the jet reaches its maximum, Fig. 9(d-e-f). The jet angle is again well predicted (Fig. 9(e-f), while the CFD-RBF simulations indicate a slightly stronger central jet velocity, Fig. 9(d). The correlation coefficient between measured and simulated velocity profiles of r 1D = 0.98 confirms good agreement. Finally, for the end diastole phase, when the valve is closed and jet inertia governs formation of the complex flow pattern interactions inside the left ventricle, results show a reasonable   agreement, Fig. 9(h-i). The remaining coherent central jet pattern and after-impinging anti-clockwise flow along the right side of the ventricle are well captured by simulations. The direct comparison of the velocity magnitude along the characteristic horizontal profile, indicated in Fig. 9 (h-i), is shown in Fig. 9(g). The CFD produced a double-peak profile, with a good match with peaks location, but the values of the characteristic peaks deviate from the experiments. The correlation coefficient between the CFD and TomoPIV profiles is poor at this location, i.e. r 1D = 0.1. This low value of the correlation coefficient is mainly due to the stronger backflow (as predicted by CFD) located in the right part of the left ventricle. Note that at this particular time instant, the flow is entirely driven by remaining jet inertia and is very sensitive to the surrounding flow variations. Even relatively small deviations between imposed working conditions of CFD and experiments will be augmented at this point of the cardiac cycle.
In addition to the characteristic velocity magnitude contours and profiles, we also present TomoPIV versus CFD correlation (Bland-Altman) plots based on the entire two-dimensional planes (the central vertical and central horizontal plane, where the CFD results are downsized to match the spatial resolution of TomoPIV) at different time instants, and we calculate corresponding integral correlation coefficients (i.e. Pearson (r 2D ) and Spearman (r s2D ), Appendix (A). In both (the central vertical and horizontal) planes, the correlation coefficients are in the 0.62 ≤ r 2D ≤ 0.91 and 0.69 ≤ r s2D ≤ 0.92 range. The obtained correlation coefficients confirm a good agreement between the CFD and experiments for analyzed cross-sections and time-instants. We conclude that an overall good qualitative and quantitative agreement is obtained between present simulations and experiments, confirming the suitability and accuracy of the developed dynamic mesh morphing approach.
The mesh dependency was checked by re-running simulations on two additional mesh levels (M1-M3) that have 4 M (coarse) and 12 M (fine) CVs, respectively. For the purposes of a short analysis, we select the horizontal profile at the early diastole (shown in Fig. 9(a)) and compare results for a different mesh in Fig. C.11. It can be seen that a good agreement is obtained between all numerical meshes, confirming good numerical accuracy and robustness of the present approach.

Limitations of the present study
A majority of the limitations of the present study can be closely associated with limitations of the referent experiment, [48]. The geometry of the left ventricle is generated from a statistically-shaped model built from 150 individuals. The present computational framework can easily deal with the patient-and/or subject-specific left-ventricles too. There is a lack of TomoPIV data in the upper part of the left ventricle and in the proximity of the biological valve, which imposed some uncertainties in reconstructing the entire shape of the left ventricle. Furthermore, the movement of tri-leaflets biological valve was recorded with a relatively low-resolution camera from the top, which can result in some deviations in calculating the exact cross-section during opening and closing phases of the cardiac cycle. In experiments, the LV configuration has the second biological valve (aortic valve), which was not modeled in the present simulations (a simple outflow boundary condition or a no-slip wall are assumed depending from particular time-instants of the cardiac cycle). The working fluid which was used in the experiments has higher dynamic viscosity compared to the real blood. This was an experimental necessity in order to get a proper matching refractive index needed for the optical measurements. Consequently, the generated flow is in a laminar flow regime, where some short transient flow regimes can be generated in the patient-specific geometries.

Conclusion
In the present study, we have numerically simulated dynamic flow patterns in a patient-specific model of the left ventricle. The numerical simulations are based on the identical working conditions as presented in recent experimental studies which combined the state-of-art tomographic laser-optics based particle-image velocimetry and echocardiography (ultrasound) measurements. The motion of the left ventricle and the biological tri-leaflet mitral valve are modeled and numerically simulated by applying a novel dynamic mesh morphing approach based on the radial basis function (RBF), within the arbitrary Lagrangian-Eulerian (ALE) discretization framework. The numerical simulations revealed all experimentally observed flow patterns inside the left ventricle. Furthermore, a direct comparison of the numerically predicted and measured velocity profiles at characteristic locations in the left ventricle at different instants of the cardiac cycle, demonstrated very good agreement, confirming the validity and accuracy of the present simulations. The image-based simulation workflow presented here can be easily adapted for similar patient-specific bio-medical applications where deforming geometries play a crucial role in generating characteristic blood flow patterns. The potentials of the RBF based CFD approach can be particularly attractive to mimic various geometrical variations caused by a disease progression (e.g., development of the stenotic regions, aneurysm growth, heart valve disease, etc.) or for an efficient geometry and numerical mesh generation for virtual surgery applications.

Declaration of competing interest
The authors declare that they have no known competing interests or personal relationships that could have influenced the work reported here.

Declaration of competing interest
None Declared.

Appendix B. Mathematical rationale of the radial basis function (RBF)
In the present work, a radial basis function (RBF) interpolation method is used to define the displacement of the internal fluid nodes from predefined control nodes at the surface of a computational domain. The interpolation function (d i ) describes the displacements in the entire domain and is approximated by a sum of basis functions as: are the coordinates of the control nodes where the displacements are known (in our case the surface nodes), p i is an arbitrary polynomial which is dependent on the kernel basis function (ϕ), and, N is the total number of control nodes. The α ij and the polynomial p i coefficients are determined by the interpolation condition: where g C is the known displacements of the control nodes, and to complete the linear system an orthogonality condition needs to be satisfied: for all polynomials q with a degree less than, or equal to degree of polynomial p i . The degree of polynomial p i depends on the basis kernel function [47]. In the present study we adopt the multi-quadratic basis kernel function where r = ||x − x Cj ||. This choice of the basis kernel function proved to be both numerically accurate and robust for the simulated left ventricle configuration with a dynamic biological valve. In combination with this radial basis function, a simple linear polynomial can be used: The unknown coefficients α and β can be obtained by solving the system where M C is the interpolation matrix defined by calculating all the radial interactions between source points M Cij = ϕ||x Ci − x Cj ||, and P C is the coordinate matrix in which the j-th row is where M m is evaluated by basis functions M mij = ϕ||x mi − x Cj || and P m is the coordinate matrix with its j-th row given as [ 1 x mj y mj z mj ] . In summary, in the mesh-morphing approach presented here, with predefined RBF, a vector field of displacement is interpolated as an independent scalar field as: