Principal Investigator/Program Director (Last, first, middle): Yoganathan, Ajit, P. GRANT NUMBER PROGRESS REPORT SUMMARY HL70262-08 PERIOD COVERED BY THIS REPORT PRINCIPAL INVESTIGATOR OR PROGRAM DIRECTOR FROM THROUGH Ajit P. Yoganathan 05/16/11 06/30/13 APPLICANT ORGANIZATION Georgia Tech Research Corp. TITLE OF PROJECT (Repeat title shown in Item 1 on first page) Computational Modeling of Mechanical Heart Valves A. Human Subjects (Complete Item 6 on the Face Page) Involvement of Human Subjects No Change Since Previous Submission Change B. Vertebrate Animals (Complete Item 7 on the Face Page) Use of Vertebrate Animals No Change Since Previous Submission Change SEE PHS 2590 INSTRUCTIONS. WOMEN AND MINORITY INCLUSION: See PHS 398 Instructions. Use Inclusion Enrollment Report Format Page and, if necessary, Targeted/Planned Enrollment Format Page. a. Specific Aims There are no changes to the original specific aims. b. Studies and Results OVERALL SUMMARY The overall progress of the grant during the last 5 years is described below. We have integrated and tested several high-resolution experimental and computational methods to develop a translational research platform for the verification and validation of medical imaging-based computational models. These models are aimed to provide predictive assessment of the performance of cardiac devices post-implantation (including but not limited to prosthetic heart valves) prior to in vivo intervention. For the first time, we have validated a computational fluid-structure interaction (FSI) solver specifically developed to simulate internal flow induced by moving boundaries as seen in the human left ventricle (LV). Two types of experimental in vitro test beds have been developed and tested for this purpose: a tetrahedral chamber with a single moving wall and an anatomical flexible-walled LV physical model. Multiple modalities of state-of-the-art experimental diagnostics, including high-speed stereo photogrammetry, particle image velocimetry (PIV) and cardiac magnetic resonance imaging (MRI) have been used to obtain high spatiotemporal resolution data of the wall motion kinematics and flow fields to use for validating FSI simulations. In addition, a separate computational method has been developed that can model pulsatile flows through prosthetic cardiovascular devices with the presence of realistic platelets. This method employs optimal parallel processing with very high spatiotemporal simulations to quantify blood damage in cardiovascular flows. This generic method can also be used for a variety of other flow conditions and medical devices as well. We have developed advanced image processing methods to reconstruct the in vivo wall motion and flow fields. We have validated the accuracies of these methods using our in vitro experimental platforms using benchmark data obtained using higher spatiotemporal experimental modalities. Finally, several sub- studies examining the interaction between anatomical components of the LV, including the mitral annulus and mitral leaflet length, on filling fluid dynamics have been investigated. We have also conducted detailed studies on the effect of varying heart rate on diastolic flow patterns to investigate the impact of mild and moderate exercise. These fundamental studies conducted under controlled conditions in vitro and in silico provide a previously unavailable physical insight into the fluid flow experienced within the LV under normal and altered functional conditions. The results of these studies can be used to optimize the design of cardiac devices such as prosthetic heart valves and mitral annuloplasty rings. Future directions of this work involve optimizing well known current designs, but also evaluating newer devices such as percutaneous valves. This can be performed as a pre-clinical assessment of novel devices. In the long-term, we intend to further our clinical translational framework to employ computational modeling for patient specific geometries for use in ____1___ PHS 2590 (Rev. 05/01) Page Form Page 5 Principal Investigator/Program Director (Last, first, middle): Yoganathan, Ajit, P. surgical planning and intervention. SPECIFIC AIM I: Development of CFD tools for left ventricle/aorta configuration 1. The Development of an efficient solver for left heart hemodynamics simulations The left heart system is a complex system consisting of several organs. Each organ has its own physiological functionality and constraints. The left heart system is depicted in Figure 1. It consists of two main parts: i) the left ventricle (LV) and ii) the aorta. We simplify the system by assuming that the atrium and the descending aorta are truncated from the whole cardiovascular systems. Appropriate boundary conditions are applied in the domain interface. For example, the inlet and outlet boundaries are specified using measured quantities. The motion of the left ventricle can be obtained via non-invasive measurements of patients or results of cell- based activation models. In order to model this complex system, we develop an overset-curvilinear immersed boundary (overset- CURVIB) method in general non-inertial frames of reference to deal with this situation. Our numerical method incorporates overset-curvilinear grids to efficiently handle multi-connected geometries and increase the resolution locally near the immersed boundaries. Complex bodies undergoing arbitrarily large deformation may be embedded within the overset-curvilinear background grid and treated as sharp interfaces using the curvilinear immersed boundary (CURVIB- (Ge and Sotiropoulos, 2007)). The blood flow (incompressible fluid) equations are formulated in a general non-inertial frame of reference to enhance the overall versatility and efficiency of the numerical approach. Efficient search algorithms to identify areas requiring blanking donor cells and interpolation coefficients for constructing the boundary conditions at grid interfaces of the overset grid are developed and implemented using efficient parallel computing communication strategies to transfer information among sub-domains. The governing equations are discretized using a second order accurate finite-volume approach and integrated in time via an efficient fractional step method. Various strategies for ensuring globally conservative interpolation at grid interfaces suitable for incompressible flow fractional step methods are implemented and evaluated. The method is verified and validated against experimental data and its capabilities are demonstrated by simulating the systolic phase in an anatomic left ventricle with an implanted mechanical heart valve in the aortic position as the final check. Figure 1. The sketch depicts the computational framework and the partition between the fluid and the solid domains. Γ and Γ are the inlet and outlet of the computational domain. Γ is the aortic inlet outlet aorta portion of the domain where the no-slip boundary condition is applied. Γ is the interface between FSI leaflets Ω and the blood flow Ω , which is simulated via the fluid-structure interaction methodology. s F The Γ represents the endocardium surface where the left ventricle beats. The kinematics of Γ can LV LV be obtained via non-invasive measurement in patient or simulation results of the cell-based model. 2. Kinematic models for patient-specific left ventricle In order to simulate the motion of the left ventricle, we reconstruct the anatomic left heart geometry from MRI data of a healthy subject. The aortic and mitral valve, coronary and carotid arteries as well as one part of the left atrium are removed from the original anatomy. The final geometry is shown in Figure 2. The values for the ____2___ PHS 2590 (Rev. 05/01) Page Form Page 5 Principal Investigator/Program Director (Last, first, middle): Yoganathan, Ajit, P. various geometrical parameters of our LV model are as follows: the LV long and short axes lengths are L=80 mm and D =47 mm, respectively; the mitral annulus diameter is D =37 mm; L M Figure 2. The left heart model reconstructed from MRI images. (r,θ, z) is the cylindrical coordinate system defined for the LV with corresponding unit vectors i, i , i . L and D are the lengths of the long r r z L and short LV axes, respectively. To model the LV kinematics we use the cell-based electric excitation methodology. The approach follows lumped-type parameter models (Beyar and Sideman, 1986; Beyar and Sideman, 1984,; Chadwick, 1982), which, nevertheless, is able to reproduce the physiological features of the LV kinematics. The model is based on the following assumptions: 1) Only the LV moves. The atrium and the LV base remain stationary in the cardiac cycle; 2) Only the endocardium is modeled as the wall surface. The endocardium surface is discretized with an unstructured grid with material nodes as shown in Figure 3. Each material node is assumed to represent one endocardium cell; and 3) The response of the endocardium cell to the cardiac electrical stimulation is assumed to be a function of a time-dependent potential p (t). Figure 3. The moving LV model, discretized with the unstructured grid is immersed in a background stationary curvilinear mesh as required by the CURVIB method. The blood flow is driven by the LV wall motion resulting from the cell-activation model. The "red" material point denotes one material point on the LV surface. To demonstrate that the emerging, large-scale LV wall kinematics resulting from the proposed model is physiologic, major global LV kinematics parameters calculated from the model are compared with available in vivo data. The most important parameter during diastole is the E-wave velocity at the mitral orifice. In our case, the E-wave peak velocity, which is estimated as the bulk velocity from the mitral inflow wave from, ____3___ PHS 2590 (Rev. 05/01) Page Form Page 5 Principal Investigator/Program Director (Last, first, middle): Yoganathan, Ajit, P. reaches 70 cm/s. This value is in the range of measurements obtained from in vivo MRI, which range from 50 to 70 cm/s (Saber et al., 2001; Schenkel et al., 2009b) . The instantaneous direction of the calculated wall velocity field is visualized in Figure 4 in terms of instantaneous LV wall streamlines at several instants during the cardiac cycle. The cell-activation model produces wall surface motion that exhibits complex twisting motion as it relaxes in the clockwise direction (looking from the base) during diastole and contracts in the counter-clockwise direction during systole. These twisting patterns are consistent with in vivo observations (Hunter et al., 2003,), in which the rotation of the apex relative to the base is in the counter-clockwise direction during systole and clockwise during diastolic untwist (as viewed from the base). Figure 4. The deformation of the endocardium surface during one cardiac cycle from diastole to systole. The instantaneous wall velocity field is indicated by the projected streamlines on the endocardium surface. The red dot in the inset shows the time instance in the cardiac cycle. The calculated LV volume rate of change dV/dt during one cardiac cycle is shown in Figure 5. During diastole, the distinct early diastolic filling peak (E-wave) is separated from the passive filling peak (A-wave) by a phase of very slow volume expansion, referred to as diastasis. An important physiologic parameter is the ratio of the E- and A-waves dV/dt, peaks (E/A ratio). As shown in Figure 5, the E/A ratio resulting from our model is approximately 1.1, which is in good agreement with physiologic values (Kim et al., 1995; Thomas and AE, 1991). In the systolic phase, the dV/dt resulting from the model becomes as low as -19 liter/min is also well within the physiologic in the order of -20 liter/min (Yoganathan et al., 2004). The details of the derivation have been published in the Journal of Computational Physics (Le and Sotiropoulos, 2012a). Figure 5. The left ventricle volume rate of change during one cardiac cycle. There are two distinct positive E-wave and A-wave peaks separated by the diastasis during diastole. The negative peak is the systolic peak. ____4___ PHS 2590 (Rev. 05/01) Page Form Page 5 Principal Investigator/Program Director (Last, first, middle): Yoganathan, Ajit, P. 3. Simulations of Pulsatile BMHV Suspension Flow using Parallelized LBM-EBF Method: An alternative numerical method is employed to study platelet damage in BMHV flows. These computational studies employ a lattice Boltzmann with external boundary force (LBM-EBF) method for platelet flow simulations. In this method, the fluid flow is solved using the lattice Boltzmann method, and each platelet is mapped onto a Lagrangian frame moving through the fluid domain. Platelets are realistically shaped as ellipsoidal particles, and not as point particles as in other limited computational studies. The two-way fluid- solid interactions are enforced with no-slip conditions by the external boundary force method. The motion and orientation of the platelet are obtained from Newtonian dynamics equations. The use of the lattice Boltzmann method for particle suspension fluid flows has been extensively validated. The LBM-EBF FSI solver has the ability to track the interaction forces (collisions and shear) of platelets flowing through BMHVs and can thus be used to evaluate blood damage. The advantages of this method are in its parallel computation (shown to scale to 96.6% of ideal efficiency at 16,384 computational cores), multiscaling (modeling particles with subgrid resolution), and realistic particle modeling (platelets with mass and volume and meshed surfaces). The first two studies using LBM-EBF have been published (Wu et al. 2011, Yun et al. 2012), where the platelet damage quantification method was validated. This was performed against experimental blood damage results and for flows through actual BMHV hinges. These papers have validated the numerical method for modeling platelet damage, particularly for BMHV flows. Figure 6. A) The computational model of the experimental setup, B) Side and angled transparent views of the SJM valve, C) 5 time points from acceleration to peak systole to valve closure and diastole for pulsatile BMHV flow. The LBM fluid solver has been validated by simulating single-phase fluid pulsatile flow through a SJM regent valve for an entire cardiac cycle and comparing the results with experiments of pulsatile flow of a blood analog fluid through a SJM valve in an in vitro loop with the same design (same geometry, leaflet motion, boundary conditions and blood kinematic viscosity). The instantaneous time snapshots of out-of-plane vorticity in the perpendicular plane are shown in 2D snapshots (Figure 6) with the highest spatiotemporal resolution to date (80µm, 2.4µs). This includes reverse flow and the ability to capture the smallest scale eddies in this complex, high Reynolds number flow. Particularly evident are vortex shedding past the leaflet tips and past the sinus expansion step. Once the leaflets close, leakage jets are formed in the central orifice and through the leaflet-valve gaps. These results were all obtained after simulation of two full cardiac cycles ____5___ PHS 2590 (Rev. 05/01) Page Form Page 5 Principal Investigator/Program Director (Last, first, middle): Yoganathan, Ajit, P. to remove the effect of initialization to peak flow. Figure 7. 5400 platelets throughout BMHV after end of systolic flow modeling, with leaflet closure. Simulations are then modeled with platelets released during pulsatile flow through the BMHV. 5400 platelets are released, with 300 platelets released every 20ms during the systolic flow phase (Figure 7). Platelets are released upstream of the valve at the same axial position but with random cross-sectional positions and orientations. Highest damaged platelets are found residing in the valve and sinus regions, despite 64% of platelets advecting downstream of the sinus at the simulation end timepoint. The highest damaged platelets are found near the leaflet surfaces, valve walls, and near the walls of the sinus expansion. Figure 8. Contour plot of platelet BDI (dynes s/cm2) – Eulerian view For Eulerian views (figure 8), the flow domain can be partitioned into small regions, and average damage in each region can be determined. From this, a platelet damage contour plot is created showing the evolution of platelet damage in space, over time. The largest number of high damage regions is found in the sinus expansion near the walls. The flow in the sinus expansion is characterized by strong vortices and recirculation regions. This is particularly dangerous as increased platelet residence time in the recirculation region could lead to platelet activation, consistent with previous results based on a model geometry. For a Lagrangian view (Figure 9), platelets with highest accumulated damage are isolated and pathlines are mapped. This can be viewed in multiple angles and connected with flow fields, showing how instantaneous platelet damage and instantaneous flow features are connected. This can also be used as a guideline for future single-phase flow simulations to interpret flow results and connect results with potential for blood damage complications. Figure 9. Particle tracking showing highlighted platelet damage path and associated flow field This baseline case is performed for a SJM valve with physiologic adult conditions, and has demonstrated the capabilities of LBM-EBF as a numerical tool to evaluate the blood damage performance of medical devices. BMHVs are a suitable initial case due to the simplified geometry and flow, and the extensive ____6___ PHS 2590 (Rev. 05/01) Page Form Page 5 Principal Investigator/Program Director (Last, first, middle): Yoganathan, Ajit, P. previous data available for comparison. The numerical method can be applied to model various cardiovascular flows and devices. This can be applied to pediatric applications for BMHVs, to be used as a predictive assessment of the fluid mechanics performance and blood damage performance of pediatric sized BMHVs. This is also useful for new devices that have not undergone significant analysis. For example, percutaneous valves or left ventricular assist devices require accurate assessment of both generic device performance and blood damage performance, which can be provided by the numerical method. SPECIFIC AIM II: In vitro experiments for validation in a phantom left ventricle/aorta configuration with moving boundaries 1. Pulsatile Flow Loop: One of the Georgia Tech Left Heart Simulators used for the experimental validation is shown in Figure 10. The model left ventricle is an acrylic truncated tetrahedron with a single flexible wall made of latex. The motion of the flexible membrane is driven by a Vivtro Superpump System pump (Model SPS3891, Vivitro Systems; British Columbia). There is a 23mm St. Jude Medical Regent BMHV is housed just below the aortic sinus in the simulator. Downstream of the aortic section (red) there is an ultrasonic flow probe (Model ME- PXN, Transonic Sonic Systems; Ithaca, NY). The compliance chamber and point resistance are used to control the hemodynamic waveforms. The mitral inflow comes from the reservoir and through another ultrasonic flow probe, and no mitral valve is positioned at the mitral orifice at this stage. Directionality of the flow was maintained using a unidirectional valve 15 diameters upstream of the mitral orifice. Flow conditioners are included 10 valve diameters upstream of the mitral orifice and are used to reduce the influence of the non-streamwise components to obtain a uniform velocity field as the inflow to the ventricle. 2D digital particle image velocimetry measurements can be performed at the mitral inflow tract, ventricle, and the flow through the BMHV in the aortic section. 50 phase locked measurements are ensemble averaged to construct velocity fields in each region. Stereo-photogrammetry is used to track the kinematics of the flexible membrane. The hemodynamic data acquisition system consists of two flow probes Q1 and Q2 (Model ME-PXN, Transonic Sonic Systems; Ithaca, NY) and three pressure transducers (TruWave Disposable Pressure Transducers, Edwards Lifesciences; Irvine, CA). The inline flow probes are used to measure the instantaneous flow upstream of the mitral orifice to the ventricle and the downstream flow through the aorta. The instantaneous absolute pressure measurements are measured upstream of the mitral valve, downstream of the aortic valve, and near the apex of the ventricle as shown in Figure 6. A trigger signal activates the acquisition of the flow field measurements using 2D digital particle image velocimetry (DPIV), flexible membrane position tracking using high-speed stereo photogrammetry, and the instantaneous hemodynamic variables. The trigger is an analog pulse created in LabView and output by a National Instruments DAQ system. The DAQ output waveform to the piston can be modified to change the contraction/relaxation timings and flow magnitude to simulate disease. The changes to the kinematic of flexible membrane and the fluid mechanics within the system can be recorded. The model described above has since been redesigned and improved both structurally and functionally. The new Georgia Tech Left Heart Simulator is shown in Figure 11. The left ventricular geometry was generated in SolidworksTM (Dassault Systèmes Solidworks Corporation, Waltham, MA, USA) by constructing a series of concentric elliptical sketches on parallel planes connected by spline curves. A 125° cut was made at the base of the ventricle such that the mitral and aortic annular planes matched in vivo conditions. The ventricle is encased in an acrylic chamber to which a pulsatile linear actuator is connected. The linear actuator causes a change in pressure in the outer chamber, which in turn forces the ventricle to contract and relax. The system is then placed in a similar flow loop as the simplified model (Figure 12). Pressure measurements can also be taken at the atrium, ventricle, and aorta. Flow measurements are taken at the inflow and outflow of the ventricle. The advantage of this simulator is in its MRI and echocardiography compatibility, which provides a direct link between bench top and clinical modalities. This allows for a comparison between planar flow fields from PCMR and PIV, wall motion from stereo-photogrammetry and CMR, mitral and aortic flow rates from PCMR and flow probes, etc. ____7___ PHS 2590 (Rev. 05/01) Page Form Page 5 Principal Investigator/Program Director (Last, first, middle): Yoganathan, Ajit, P. Figure 10. Left: The simplified model of a left ventricle with mitral inlet (blue) and the aortic outlet (red). The BMHV is implanted at the aortic position just below the aortic sinus. Right: schematic of the ventricular chamber with membrane (blue), inflow and out flow tracts. Figure 11. Left: The realistic model of a left ventricle with mitral inlet and the aortic outlet. BMHVs are implanted at the mitral annulus and aortic annulus just below the aortic sinus. Right: schematic of the realistic ventricular bag membrane. Figure 12. Schematic of new flow loop for the realistic Georgia Tech Left Heart Simulator. ____8___ PHS 2590 (Rev. 05/01) Page Form Page 5 Principal Investigator/Program Director (Last, first, middle): Yoganathan, Ajit, P. For the validation of FSI models of the left heart complex, three sets of experiments were conducted on phantom left ventricle models. The design of each experiment was selected to gradually include more complexity in the system. The first set of experiments used the idealized left ventricle model consisting of a tetrahedral chamber. The fluid dynamics of diastolic filling through inclined inflow nozzles, coarsely representing the asymmetry of mitral leaflets, was investigated as a part of this sub-study. Numerical simulations for this setup were conducted in tandem and the flow field results were compared to experiments for validation. The second and third sets of experiments were conducted on the anatomical left ventricle model, to examine the effects of varying the mitral annulus shape (between D- and O-orifices) and heart rate (100 bpm and 120 bpm). 2. Validation of vortex formation phenomenon in complex geometries 2.1 Vortex ring formation from inclined nozzles In order to validate our computational framework to simulate complex flows in the left ventricle, we validate our numerical procedures with well-controlled experimental data. The first validation is the investigation of vortex dynamics in laminar, impulsively driven flows through inclined nozzles in a piston-cylinder apparatus. Our simulations are motivated by the need to provide a complete description of the intricate vortical structures and governing mechanisms emerging in such flows as documented in the experiments of (Troolin and Longmire, 2010). Our computational result agrees well with the experimental data as shown in Figure . Figure 13. Comparison between a) measured (Webster and Longmire) (left) and b) computed (right) instantaneous non-dimensional, out-of-plane vorticity contours for the inclined nozzle on the symmetry plane. The first contour is 5 U/D and the increment is 2.5 U/D. Dash lines indicate negative vorticity values. Ring 1 and 2 are the primary and secondary (stopping) rings. We show that the flow is dominated by the interaction of two main vortical structures: the primary inclined vortex ring at the nozzle exit and the secondary stopping ring that arises due to the entrainment of the flow into the cylinder when the piston stops moving. These two structures are connected together with pairs of vortex tubes, which evolve from the continuous vortex sheet initially connecting the primary vortex ring with the interior cylinder wall as seen in Figure . In the exterior of the nozzle the key mechanism responsible for the breakup of the vortical structure is the interaction of the stronger inclined primary ring with the weaker stopping ring near the longest lip of the nozzle. In the interior of the nozzle the dynamics are governed by the axial stretching of the secondary ring and the ultimate impingement of this ring on the cylinder wall. Our simulations also clarify the kinematics of the azimuthal flow along the core of the primary vortex ring. We show that the azimuthal flow is characterized by a pair of two spiral saddle foci at the long and short lips of the nozzle through which ambient flow enters and exits the primary vortex core. The details of the comparison and computational results can be found in the Journal of Fluid Mechanics (Le et al., 2011). ____9___ PHS 2590 (Rev. 05/01) Page Form Page 5 Principal Investigator/Program Director (Last, first, middle): Yoganathan, Ajit, P. Figure 14. The three-dimensional topology and evolution of the vortical structure for the inclined nozzle visualized by plotting the iso-surface of non-dimensional vorticity magnitude colored with contours of non-dimensional helicity density. The vorticity dynamics during this early stage is characterized by the formation and interaction of the primary ring (R1) and the stopping ring. 2.2 Vortex ring formation in idealized left ventricle models In order to examine the vortex formation process of the left ventricle in a well-controlled condition, we carry out an experiment with an idealized model of the LV. The idealized LV consists of a tetrahedral chamber representing the LV. The mitral and aortic orifices are connected to a circular pipe as shown in Figure . Particle Imaging Velocimetry (PIV) measurements are carried out on representative planes to provide sufficient data in order to validate the numerical model. The computational model closely follows the experimental setup. The diameter of the mitral orifice is D = 1 inch (25.4 mm). The viscosity of the working fluid is set as the experimental value of 3.5x10-6 m2/s. The inflow waveform at the mitral orifice is shown in Figure with the peak flow rate of Q = 16 L/min and the bulk velocity corresponds to this flux is U = 0.52 m/s.. The peak Reynolds number is therefore Re = UD/υ= 3943. A structured grid with size of 45 million grid points is used to simulate this case. The computational grid is a box covering the size of the LV chamber. The aortic outlet is specified using homogeneous Neumann boundary condition (no resistance). The simulation results are compared with the experimental data on the measured plane shown in Figure . Good agreement between the experiment and numerical simulation are shown at early diastolic filling as seen in Figure at t=300 ms. All important features of the flow are captured accurately by the computational model. The main structure is an asymmetrical vortex ring, which is indicated by the negative and positive vorticity cores as shown in Figure . At early diastole, this vortex ring is formed at the mitral orifice tip and propagates toward the apex of the LV chamber. Note that due to the asymmetrical configuration of the geometry, the vortex core near the aortic orifice (negative) grows significantly faster than the other. The negative core also impacts directly with the LV solid wall right after formation. The computational model also simulates correctly the vorticity extraction created by the impingement of this vortex ring on the inclined wall surface. The three- dimensional structure of this vortex ring is shown in Figure . The vortex core thickness varies continuously along the ring’s circumference indicating a deviation from the perfectly circular ring. ____10___ PHS 2590 (Rev. 05/01) Page Form Page 5