Fluid–Structure Interaction Models of Bioprosthetic Heart Valve Dynamics in an Experimental Pulse Duplicator

Computer modeling and simulation is a powerful tool for assessing the performance of medical devices such as bioprosthetic heart valves (BHVs) that promises to accelerate device design and regulation. This study describes work to develop dynamic computer models of BHVs in the aortic test section of an experimental pulse-duplicator platform that is used in academia, industry, and regulatory agencies to assess BHV performance. These computational models are based on a hyperelastic finite element extension of the immersed boundary method for fluid–structure interaction (FSI). We focus on porcine tissue and bovine pericardial BHVs, which are commonly used in surgical valve replacement. We compare our numerical simulations to experimental data from two similar pulse duplicators, including a commercial ViVitro system and a custom platform related to the ViVitro pulse duplicator. Excellent agreement is demonstrated between the computational and experimental results for bulk flow rates, pressures, valve open areas, and the timing of valve opening and closure in conditions commonly used to assess BHV performance. In addition, reasonable agreement is demonstrated for quantitative measures of leaflet kinematics under these same conditions. This work represents a step towards the experimental validation of this FSI modeling platform for evaluating BHVs. Electronic supplementary material The online version of this article (10.1007/s10439-020-02466-4) contains supplementary material, which is available to authorized users.


B Parameter Fitting for the Valve Material Models
We perform model fits for the constitutive model of the valves detailed in the "Leaflet Mechanics" section in the main text using experimental data from Billiar and Sacks [1,2] for the porcine tissue BHV and from Kim et al. [3] for the bovine pericardial BHV. In the biaxial tests of Billiar and Sacks, the specimens were mounted onto the testing device so that the radial and circumferential directions were aligned with the direction of the applied forces, as shown in Figure S1a. Kim et al. used a stress-control biaxial testing method on a gluteraldehyde-treated bovine pericardial tissue sample. As shown in Figure S1c, the specimen was aligned at 45 • to the direction of the applied forces, mimicking the fiber alignment of the pericardial valve leaflets. To determine constitutive model parameters, we assume that the solid is incompressible, compute the second Piola-Kirchhoff stress (S), and compare it to the experimental values for given values of the Green-Lagrange strain E = 1 2 (C − I). The stress is computed by in which W = W iso + W aniso is the Holzapfel-Gasser-Ogden model described by (4) and (5) in the "Leaflet Mechanics" section of the main text, S dev is the deviatoric part of the second Piola-Kirchhoff stress, C = F T F is the Cauchy-Green strain, F is the deformation gradient tensor, and p is the pressure assuming that the specimen is acted on only by in-plane loads. We use lsqcurvefit in MATLAB (The MathWorks, Inc., Natick, MA, USA), a nonlinear least-squares solver, to determine the model parameters. Figures S1b and S1d show constitutive model fits obtained using this approach.  Figure S1: (a) Schematic of the biaxial tensile test of Billiar and Sacks [1,2] for porcine aortic valve tissue specimens to study their material response. X C is the circumferential axis, and X R is the radial axis. The dashed lines are the directions in which forces were applied. (b) Parameter fitting for the porcine aortic valve using the experimental tensile test data from Billiar and Sacks [1,2] for glutaraldehyde-fixed porcine aortic valves. (c) Schematic of the biaxial tensile tests of Kim et al. [3] for bovine pericardium tissue specimens to study their material response. X 1 is the preferred mean fiber direction (PD) also shown by gray lines, X 2 is the cross-preferred fiber direction (XD), and X 1 and X 2 are the directions in which forces were applied. (d) Parameter fitting for the bovine pericardial valve using the equibiaxial data from Kim et al. [3] compared to the plot using parameters determined by Kim et al. using a finite element model of the biaxial test.

C Parameter Fitting for the Windkessel Model
Pulse-duplicator system components upstream of the aortic test section, including the resistance and compliance of the pump, the viscoelastic impedance adapter (VIA) subsystem, and the left ventricular chamber of the pulse duplicator, are described by a three-element Windkessel (R-C-R) model for the porcine aortic valve, whereas a more detailed model is used for the bovine pericardial valve because additional experimental data are available for the bovine pericardial BHV cases. See Figure 3 of the main text. We remark that the VIA subsystem includes tunable compliance chambers (see Figure 1b of the main text) that help to ameliorate unphysiological pressure oscillations in the left ventricular pressure [4]. The upstream model used in the porcine aortic valve simulations is in which C VIA , R 1 , and R 2 characterize the VIA system, and Q LV and P LV are the volumetric flow rate and pressure, respectively, at the inlet of the aortic test section model. P pump is derived from (S4) and (S5) using experimental measurements of Q LV and P LV . The upstream model used in the bovine pericardial valve simulations is in which C VIA1 , C VIA2 , R VIA , and R out characterize the VIA system, R MV characterizes the resistance at the mitral position, P LA is the left atrial pressure, Q pump is the prescribed volumetric flow rate of the pump, and Q LV and P LV are the volumetric flow rate and pressure, respectively, at the inlet of the aortic test section model. The downstream loading conditions of the pulse duplicator are also described by a three-element Windkessel model, in which C is the compliance, R c is the characteristic resistance, R p is the peripheral resistance, P Wk is the Windkessel pressure, and Q Ao and P Ao are the volumetric flow rate and pressure, respectively, at the outlet of the test section. The Windkessel parameters R c , R p , and C are calibrated using experimental measurements of Q Ao and P Ao . We use the nonlinear optimization tool fmincon in MATLAB to determine the model parameters by comparing P Ao from (S9) and (S10) with that of experimental system ( Figure S2). When calibrating the reduced-order models, the experimental measurements of Q Ao are used as inputs to the Windkessel models. A first-order (Godunov) time step splitting is used to couple the three-dimensional FSI model to the reduced-order models, as described by Griffith et al. [5]. At the outlet of the aortic test section, for example, the mean flow rate through the outlet surface computed by the three-dimensional FSI model is provided as an input to the Windkessel model. The pressure generated in the Windkessel model is used, along with the flow rate, to determine the pressure along the outlet surface, which is used as a boundary condition for the three-dimensional model.

D.1 Continuum Equations
The IB formulation uses an Eulerian description of the momentum, viscosity, and incompressibility of the coupled fluid-structure system, and it uses a Lagrangian description of the deformations, stresses, and resul-tant forces of the immersed structure. In particular, we model the immersed structure as a viscoelastic solid, in which the viscous stresses in the solid are small compared to elastic stresses, as in previous work by us and others [6][7][8][9][10][11]. The coupling between Eulerian and Lagrangian variables is mediated by integral transforms with delta function kernels. In the continuous IB formulation, we consider a fixed three-dimensional Eulerian computational domain Ω that is divided into the time-dependent solid (Ω s t ) and fluid (Ω f t ) subdomains, so that Ω = Ω s t ∪ Ω f t . Here, x = (x 1 , x 2 , x 3 ) ∈ Ω are physical coordinates, X = (X 1 , X 2 , X 3 ) ∈ Ω s 0 are reference coordinates attached to the structure, N (X) ∈ ∂Ω s 0 is the outward unit normal to the reference configuration of the solid region at material point X, and χ(X, t) ∈ Ω s t is the physical position of material point X at time t. ρ and µ are the mass density and viscosity, u(x, t) and p(x, t) are the Eulerian velocity and pressure fields, f (x, t) is the Eulerian elastic force density, P(X, t) is the first Piola-Kirchhoff elastic stress of the immersed structure, and δ(x) = 3 i=1 δ(x i ) is the three-dimensional Dirac delta function. The IB form of the equations of motion for the fluid-structure system is: in which D Dt = ∂ ∂t + u · ∇ is the material derivative. Notice that because ∂χ/∂t(X, t) = u(χ(X, t), t) and ∇ · u(x, t) = 0, the immersed structure is described in the continuous IB framework as an exactly incompressible material.

D.2 Numerical Approximations
The computational domain Ω, which includes both the solid and fluid subregions, is described using a blockstructured locally refined Cartesian grid consisting of multiple nested levels of Cartesian grid patches. High spatial resolution is used dynamically near fluid-structure interfaces ( Figure S3) and near flow features, such as vortices shed from the valve leaflets, that are identified by feature detection criteria that select regions of high vorticity for enhanced spatial resolution. A staggered-grid discretization is used for the incompressible Navier-Stokes that includes a version [12,13] of the xsPPM7 variant [14] of the piecewise parabolic method (PPM) [15] to approximate the convective term.
In our simulations, we use a regularized delta function δ h (x) to approximate the singular Dirac delta function δ(x). The three-dimensional regularized delta function is the tensor product of one-dimensional delta functions, In our computations, different regularized delta functions are used for the rigid aortic test section and flexible valve leaflets. These choices are based on preliminary benchmarking studies extending the tests reported by Griffith and Luo [11]. For the aortic test section, we use a simple piecewise-linear (PWL) kernel that takes the form, and we use a three-point B-spline (BS3) kernel for the valve leaflets, (S16) Figure S3: The aortic test section in the computational domain. The computational domain is a blockstructured adaptively refined Cartesian grid. A snapshot of the locally refined grid is shown on a plane that bisects the aortic test section.
Notice that the gap between the valve leaflets during closure reflects the size of regularized delta function, and in our simulations, the valve is effectively closed during the diastolic phase of the simulated cardiac cycle. This is also shown by Figures 4b and 4d in the main text, in which there is no leakage once the valve is closed. An advantage of the IB formulation of such problems is that it does not require an explicit contact model to simulate contact between the valve leaflets.

E Stabilization Method for the Hyperelastic Immersed Boundary Method
In the continuum equations of the IB framework (Supplemental Materials Section D.1), the Cauchy stress on the full domain can be seen to take the form on the Lagrangian finite element mesh. Consequently, the method only approximately accounts for the incompressibility of the immersed structure, although the accuracy of this approximation improves under grid refinement. Vadala-Roth et al. [16] introduced a volumetric energy in the solid region to stabilize and improve the accuracy of the numerical scheme substantially, in which the stabilization term p stab = −∂U (J)/∂J acts as a pressure that reinforces incompressibility in the solid region in the numerical method. We refer to this as a stabilization method because its effect vanishes under grid refinement. Specifically, in the continuum limit, the structure is automatically incompressible in the IB framework, and so J ≡ 1 and p stab ≡ 0. This stabilization method has been tested using standard quasi-static solid mechanics benchmark problems.

E.1 Compression Benchmark
This benchmark involves a rectangular block with a downward traction of 200 dyn cm −2 applied at the center of the top side along with zero horizontal displacement. Zero vertical displacement is imposed at the bottom side of the block, and zero traction is applied at all sides other than the top side [17], as shown in Figure S4a. The block is a neo-Hookean material, in whichĪ 1 = tr(C) is the first invariant of the modified right Cauchy-Green strain tensor Here, c 1 is the shear modulus, and κ stab is the numerical bulk modulus controlled by the shear modulus and the numerical Poisson ratio, ν stab , via in which µ s is the linearized (F = I) shear modulus. For the neo-Hookean model, µ s = c 1 . This test uses c 1 = 80.194 dyn cm −2 and ν stab = 0.4. The density and viscosity are ρ = 1.0 g cm −3 and µ = 0.01 dyn s cm −2 , respectively. We again emphasize that the numerical bulk modulus, κ stab , and the numerical Poisson ratio, ν stab , act as stabilization parameters, and not as physical parameters of the model, and that their effects vanish under grid refinement. This test shows that using the stabilization method yields superior volume conservation and smooth deformations for the mesh elements ( Figures S4b and S4d).

E.2 Torsion Benchmark
This benchmark involves applying torsion to an elastic beam [18] as shown in Figure S5a. Torsion is applied to the top face of the beam while zero displacement is imposed at the opposite face. Zero traction boundary conditions were used on all the other faces. The elastic beam is a Mooney-Rivlin material, This test uses material parameters c 1 = 9000 dyn cm −2 , c 2 = 9000 dyn cm −2 , and ν stab = 0.4. The density and viscosity are ρ = 1.0 g cm −3 and µ = 0.04 dyn s cm −2 , respectively. This test also shows that the stabilization method gives reasonable deformation ( Figure S5c) as well as yielding superior volume conservation ( Figures S5b and S5d).

F Solution Verification
We perform grid convergence studies to verify the level of consistency in our numerical results. We consider both dynamic and steady flow conditions at three different spatial resolutions. For the dynamic cases, we consider the consistency under grid refinement of all of the measurements reported in the manuscript. For the steady flow conditions, we observe the time-averaged flow rates and the turbulence kinetic energy (TKE) between different resolutions, in addition to bulk measurements. We use locally refined Cartesian grids with an effective fine-grid resolution that corresponds to a uniform  Figure S6 compares the bulk flow properties under grid refinement, indicating that our simulations are consistent under grid refinement. This suggests that the results reported in the main text are sufficiently resolved for these measurements by the resolutions considered therein. Tables S1, S2, S3, and S4 quantify the discrepancies relative to the experimental data under grid refinement, as well as the relative difference between resolutions (N = 192, N = 256, and N = 320). The relative differences are calculated by in which (0, T ) indicates an integral over time, and q = 2 or ∞ for the L 2 -and L ∞ -norms, respectively. This shows that the discrepancies in bulk flow measurements compared to the experimental data are within 8.7% at all resolutions. The relative differences between successive resolutions are within 1.6%, indicating consistency under grid refinement. This suggests that the bulk properties are reasonably resolved by our FSI model of the BHVs at the higher spatial resolutions (N = 256 and N = 320).  The comparison of PDVA results is similar to that of other bulk measurements under grid refinement, as * Comparisons of L 2 -norms of the discrepancies in the bulk measurements from the simulation and experiment shown in Figure S6 under the finer (N = 256) and finest (N = 320) resolutions, normalized by the L 2 -norms of the measurements from the experiment. Relative indicates the L 2 -norms of the differences between the measurements from the two resolutions, normalized by the L 2 -norms of the measurements from the finer (N = 320) case. * Comparisons of L ∞ -norms of the discrepancies in the bulk measurements from the simulation and experiment shown in Figure S6 between the coarser (N = 192) and finer (N = 256) resolutions, normalized by the L ∞ -norms of the measurements from the experiment. Relative indicates the L ∞ -norms of the differences between the measurements from the two resolutions, normalized by the L ∞ -norms of the measurements from the finer (N = 256) case. * Comparisons of L ∞ -norms of the discrepancies in the bulk measurements from the simulation and experiment shown in Figure S6 between the finer (N = 256) and finest (N = 320) resolutions, normalized by the L ∞ -norms of the measurements from the experiment. Relative indicates the L ∞ -norms of the differences between the measurements from the two resolutions, normalized by the L ∞ -norms of the measurements from the finer (N = 320) case. observed in Figure S6. Figures S8 and S9 also give a qualitative sense of model consistency under grid refinement. Although Figures S8 and S9 indicate that additional small-scale flow features are generated by the model at higher spatial resolutions, the large scale flow features are consistent between grid resolutions.
We also observe that the time-averaged flow fields (Eq. S26) as well as the flow profiles (Figures S10 and S11) are in reasonable agreement between different spatial resolutions. The mean flow is determined via We consider steady inflow conditions at two different flow rates: 100% and 50% of the maximum flow rate from the fully dynamic simulations. In particular, we find that the difference between N = 256 and N = 320  Figure S6d.
cases are small, and we conclude that N = 256 results are reasonably well resolved for large-scale flow measures. We observe that the turbulence kinetic energy (Eq. S27) as well as its profiles for each bovine pericardial valve case (Figures S12 and S13) still show some differences between different spatial resolutions. Turbulence kinetic energy is defined by in which u i = u i − u i is the fluctuating velocity. We again consider steady inflow conditions at two different flow rates: 100% and 50% of the maximum flow rate. At the maximum flow rate, the axial flow fields agree to approximately 1.6% in the L 2 norm for N = 256 and N = 320, with peak velocities differing by approximately 7%. At 50% of maximum flow, the flow fields agree to within approximately 1.1% in the L 2 norm, with peak velocities agreeing to within approximately 4%. The differences are larger for the TKE, with L 2 differences on the order of 10% and maximum TKE values differing by 17.5% at the maximum flow rate but to within approximately 1.5% at 50% of the maximum flow rate. The mean flow and TKE profiles are in reasonable agreement and appear to be converging under grid refinement. Although we find that the differences between N = 256 and N = 320 are not large, these results suggest that we are able to resolve the large-scale flow features, but higher spatial resolution is clearly needed to fully resolve the flow fields.
With steady inflow conditions, as in the pulsatile cases, the computed values of PDVA yielded by the model are in good agreement at the different grid spacings (data not shown). It is apparent, however, that the peak flow velocities for the coarser cases are larger than those of the finer cases. The observed differences in the peak flow velocities are a consequence of the numerical scheme that is used to couple the Lagrangian and Eulerian variables in the present IB method. Specifically, the use of regularized delta functions in the Lagrangian-Eulerian coupling operators causes the "effective" sizes of the Lagrangian structures to be larger than their "actual" sizes, with a length scale that is related to the regularization width of the delta function. In practice, the width of the regularized delta function is tied to the Cartesian grid spacing. This results in an effective narrowing of the flow path that is on the order of a few Cartesian grid cell widths, which is more pronounced on coarser grids. The effect of this narrowing is apparent in Figures S10 and S11 and results in higher flow velocities at the same volumetric flow rates. This effect also diminishes under grid refinement, as observed in the grid convergence study.
Figures S14 and S15 qualitatively compare the leaflet kinematics for N = 192, N = 256, and N = 320 during valve opening and closure for the porcine aortic valve and bovine pericardial valve, respectively. These results indicate that the computed structural deformations are consistent under grid refinement. Comparing the stress distributions of the leaflets in two grid resolutions (Figures S16 and S17) also confirms that the stress distributions are consistent under grid refinement.

G Comparison Between Saline and Newtonian Blood Analogue Test Fluids
The goal of this study is to present our FSI models of bioprosthetic heart valves in the in vitro pulse duplicator system in the same setting as the completed experiments, which use saline as the test fluid. Even so, we performed additional simulations with a Newtonian blood analogue fluid with density 1.0 g cm −3 and viscosity 3.5 cP to show that our FSI models can be used in such flow regime without numerical instabilities ( Figures S18 and S19). We are also interested in understanding the influence of using saline as compared to a more realistic blood analogue fluid in in vitro testing. To compare the predictions of the model for different test fluids, we perform simulations at the finer resolution (N = 256) for both pulsatile and steady mean simulations. For the comparison of the steady mean simulation, we only show the results for the bovine pericardial valve case. The porcine aortic valve model yields similar results and is not shown.
Our simulations indicate that the large scale flow structures are very similar between saline and the blood analogue fluid (Figures S18, S19, and S20). However, we observed differences in the small-scale turbulent flow structures present between saline and blood analogue fluid case ( Figure S21), especially in the maximum steady flow case ( Figure S21a). This is expected because of the larger Reynolds number associated with saline compared to the blood analogue.

H Simulation Environment
For N = 256, FSI simulations take about six days to compute a full cardiac cycle using 64 processor cores on Dogwood cluster managed by University of North Carolina Information Technology Services Research Computing. In this study, we focus on modeling aspects rather than parallel performance analysis. There is currently an ongoing effort to improve the parallel performance of the IBAMR software library, and we aim to report on these efforts, and their impacts on models like those described herein, in the future.