Computer modeling and simulation (CM&S) 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 ﬁnite element extension of the immersed boundary method for ﬂuid–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 ﬂow 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 leaﬂet kinematics under these same conditions. This work represents a step towards the experimental validation of this FSI modeling platform for evaluating BHVs. Keywords: immersed boundary method, ﬁnite element method, porcine aortic valve, bovine pericardial valve 1 1 Introduction Worldwide, nearly 300,000 aortic valve replacements are performed annually, and the rate of heart valve replacement is projected to exceed 850,000/year by 2050 [1]. Bioprosthetic heart valves (BHVs) are commonly constructed from ﬁxed porcine heart valves or bovine or porcine pericardial tissues. BHVs generate ﬂow patterns that mimic those of the normal human aortic valve and typically allow patients to be managed without chronic anticoagulation [2]. Unfortunately, BHVs have a limited durability, and they typically fail 10–15 years post implantation, primarily from tissue degeneration or calciﬁcation [3]. Computer modeling and simulation (CM&S) can be used throughout the life cycle of prosthetic valve design and regulatory approval. CM&S is a cost- and time-ecient complement to traditional bench testing that can assess device performance under a broader range of conditions than those listed in the instructions for use, including patient-speciﬁc conditions [4]. Simulations may be used in the design phase to optimize device design. Credible simulation data may also be leveraged in regulatory applications to support claims of device safety and e↵ectiveness. Indeed, using modeling and simulation to support regulatory decision- making is a strategic priority area for the U.S. FDA Center for Devices and Radiological Health [5]. CM&S is well-suited for performing root cause analyses to understand impaired device function. For example, CM&S provides information about valve performance that is dicult to acquire on the bench, including assessing the impact of non-circular conﬁgurations on transcatheter BHVs [6]. To realize its full impact in device regulation, however, CM&S results must be shown to be credible through veriﬁcation and validation (V&V) [7]. A ﬂuid–structure interaction (FSI) approach is necessary to model heart valves across the full cardiac cycle [8]. Accounting for coupling between the ﬂexible valve leaﬂets and the ﬂuid ﬂow is crucial in studying the e↵ect of vortices in the aortic sinuses, predicting ﬂuid-induced shear stress on the leaﬂets, and assessing valve performance by quantifying the valve oriﬁce area and regurgitation [9]. A widely used approach to simulating cardiovascular FSI is the arbitrary Lagrangian–Eulerian (ALE) method [10,11], which uses body-conforming meshes for the ﬂuid and solid. ALE methods have realized limited success in simulating the dynamics of heart valves to date, however, because of the substantial challenges posed by dynamically generating geometrically conforming discretizations of thin structures that undergo substantial motion [10,11]. Non-body conforming discretizations, which avoid the diculties of body-ﬁtted grids, are now widely used to model heart valve dynamics. One of the earliest of these types of approaches is the immersed boundary (IB) method [12], which was introduced by Peskin to simulate heart valves [13, 14]. The IB formulation allows the structural discretization to be independent of the ﬂuid grid and thereby facilitates models with very large structural deformations [15]. Extensions to the IB method have also been used to simulate heart valves [16, 17], but in most cases, these prior simulations were not fully resolved, and the valve leaﬂets were described using only simple structural models based on linear elasticity. Other studies [18, 19] used anisotropic models of native and bioprosthetic valves but do not address the diastolic phase of the cardiac cycle. Flamini et al. [20] used an IB approach to simulate aortic valve FSI across multiple cardiac cycles, but they used a simpliﬁed description of the aortic valve mechanics. Hasan et al. [21] also used an IB approach to simulate FSI in a subject-speciﬁc aortic root model, and they used realistic hyperelastic constitutive models to describe the valve leaﬂets, but their model has not been validated. In addition, neither study used subject-speciﬁc driving or loading conditions. Methods also have been developed that combine features of ALE and IB-like approaches, including the hybrid ﬁctitious domain/ALE method [22] and the immersogeometric (IMGA) method [9, 23–26]. These methods also seek to relax the need to use body-conforming discretizations. Prior models of aortic valves using a ﬁctitious domain/ALE method [22] showed instabilities when simulated under physiological Reynolds numbers and transvalvular pressures. The IMGA method has been used in several studies on pericardial BHVs that include experimentally derived constitutive models of the valve leaﬂets. Hsu et al. [23] validated the leaﬂet kinematics of their model by comparing the cross-sectional proﬁles of the leaﬂets to those from dynamic in vitro experimental measurements [27]. Xu et al. [25] used in vivo imaging data to drive their simulations and compared the ﬂuid ﬂow patterns with those from magnetic resonance imaging. These models both employed isotropic descriptions of the pericardial BHV leaﬂets. Wu et al. [26] used an anisotropic model that incorporated the ﬁber structure of a bovine pericardial valve that was developed and validated with experimental data from Sun and Sacks [28]. These prior studies did not use boundary condition models that established ﬂow conditions directly comparable to available experimental or clinical data. 2 Other methods that avoid body-conforming discretizations include smoothed particle hydrodynamics (SPH) [29]. The SPH approach [29], however, cannot fully impose incompressibility, and prescribing bound- ary conditions is dicult. In addition, several studies have used cut-cell-like methods [30–32] to simulate valve FSI, but these studies also used linearly elastic leaﬂet models and have not yet been validated. To date, there have been relatively few experimental validation studies of FSI models of native or bio- prosthetic aortic valves [33, 34]. The two studies by Tang et al. [33] and Sig¨ uenza et al. [34] use idealized isotropic leaﬂet models that do not account for leaﬂet anisotropy and, as a result, their simulations show discrepancies in the dynamics of the valve compared with experimental data. These studies also do not use ﬂow domains that are long enough to allow the complex ﬂow patterns downstream of the valve to develop fully. The model of Sig¨ uenza et al. exhibits incomplete closure and experiences substantial regurgitation dur- ing diastole, leading to an underestimation of the transvalvular pressure gradient. The study of Sig¨ uenza et al. also uses ﬂow rate boundary conditions throughout the cycle, which introduces bias in the valve dynamics. This study develops a computational FSI model based on the IB method of an experimental pulse- duplicator platform for simulating BHV dynamics. The model is calibrated using relatively limited ex- perimental data, and this study describes initial work towards the V&V of this model for porcine tissue and bovine pericardial BHVs. Our models of the leaﬂet mechanics are based on experimental tensile test data [35–37] of ﬁxed tissues that are similar to the biomaterials used to construct porcine aortic and bovine pericardial BHVs. To provide realistic driving and loading conditions across the full cardiac cycle, we use reduced-order models that are calibrated using pressure and ﬂow data acquired from the pulse-duplicator systems. Simulation results are compared to experimental data obtained from a commercial ViVitro pulse duplicator and a customized experimental apparatus based on the ViVitro pulse duplicator. 2 Materials and Methods 2.1 Experimental Pulse Duplicator In vitro experiments are performed using two di↵erent experimental platforms, including a ViVitro Pulse Duplicator (ViVitro Labs, Inc., Victoria, BC, Canada) available through the FDA Cardiac Device Flow Lab and a customized pulse duplicator developed by Scotten that is similar to the commercial ViVitro system. Figure 1a details Scotten’s customized system. It includes a prototype Leonardo electro-optical sub- system [38] to assess projected dynamic valve area (PDVA). The commercial ViVitro system is similar but uses high-speed videography to assess valve kinematics, from which we reconstruct PDVA data via automatic image analysis in DataTank (Visual Data Tools, Inc., Chapel Hill, NC, USA). We also measure ﬂow rates and pressures, as indicated in Figure 1b. Experiments in the customized pulse duplicator use a 25 mm Labcor TLBP A Supra (Labcor Laborat´ orios Ltd., Belo Horizonte, Brazil) porcine aortic valve. Flow and pressure signals are ﬁltered at 100 Hz, and PDVA signals are not ﬁltered. Experiments at the FDA use a Model 2800 25 mm Carpentier–Edwards PERIMOUNT RSR (Edwards Lifesciences, Irvine, CA, USA) bovine pericardial aortic valve. Flow signals are ﬁltered at 100 Hz, and pressure signals are not ﬁltered. Both experiments use saline as the test ﬂuid for this initial study because the viscosities of common blood analogues are more sensitive to temperature than saline. In addition, saline is widely used for hydrodynamic assessments of heart valves [39]. For instance, the ISO 5840-3 standard allows saline to be used as a test ﬂuid [40], and the most recent inter-laboratory study by Wu et al. [41] also used saline as the test ﬂuid. Results presented in Supplemental Materials Section G examine the e↵ect of using a Newtonian blood analogue instead of saline in the computational model. 3 (a) (b) Figure 1: (a) A customized pulse duplicator based on the commercial ViVitro pulse-duplicator system adapted with prototype electro-optical subsystem for measuring aortic valve projected dynamic valve area (PDVA), or alternate conﬁguration for measuring mitral valve PDVA. (b) A schematic diagram of the custom pulse duplicator based on the commercial ViVitro pulse-duplicator system adapted with prototype electro- optical subsystem for measuring aortic valve PDVA, or alternate conﬁguration for measuring mitral valve PDVA. The commercial ViVitro system is similar but lacks the back light and the photo sensor for acquiring PDVA. 2.2 Solid Mechanics Models 2.2.1 Leaﬂet Mechanics We describe the biomechanics of the valve leaﬂets using the framework of nonlinear solid mechanics [42]. Brieﬂy, leaﬂet deformations are described by the mapping =(X,t) between reference coordinates Xand 4 current coordinates xat time t. The valve leaﬂets are treated as anisotropic, incompressible, hyperelastic materials. For a hyperelastic material, the ﬁrst Piola–Kirchho↵stress Pis related to a strain-energy functional (F) via [21] P=@ @F, (1) in which F=@/@Xis the deformation gradient tensor. We split the strain energy functional into isochoric and volumetric parts, (F)=W(F)+U(J), (2) in which F=J1/3FandJ=d e t F. We use a modiﬁed version of the Holzapfel–Gasser–Ogden model from Murdock et al. [43] with the addition of angle dispersion [44]. The isochoric part of the model, W(F), includes an isotropic contribution from the extracellular matrix and an anisotropic contribution from the collagen ﬁbers embedded in the leaﬂets, W(F)=Wiso(F)+Waniso(F). (3) This model describes the extracellular matrix as an exponential neo-Hookean material along with a collagen ﬁber reinforcement model accounting for angle dispersion, Wiso(F)=C10{exp⇥ C01(¯I13)⇤ 1}, (4) Waniso(F)=k1 2k2{exp⇥ k2(¯I1+( 1 3)¯I? 41)2⇤ 1}, (5) in which ¯I1=t r ( C) is the ﬁrst invariant of the modiﬁed right Cauchy–Green strain tensor C=FTF, ¯I? 4= max( ¯I4,1) = max( eT 0¯Ce0,1), and e0is a unit vector aligned with the mean ﬁber direction in the reference conﬁguration. By construction, ¯I? 4is nonzero in extension but not in compression. The parameter 2[0,1 3] describes collagen ﬁber angle dispersion. If = 0, the ﬁbers are perfectly aligned and the constitutive model becomes the same as that of Murdock et al. [43] By contrast, if =1 3,t h e nt h em o d e l describes an isotropic distribution of ﬁbers, and we obtain an isotropic model with no preferred direction of ﬁber reinforcement. The parameters for the porcine aortic valve are ﬁt to experimental tensile test data from Billiar and Sacks [35,36] for glutaraldehyde-ﬁxed porcine aortic valves (Fig. 2a). We obtain C10=0.302 kPa, C01=3.25,k1=0.197 MPa, k2=0.001, and =0.0. The parameters for the bovine pericardial valve are ﬁt to biaxial data from Kim et al. [37], and we obtain C10=0.119 kPa, C01= 22 .59,k1=2.38 MPa, k2= 149 .8, and =0.292. For further details, see Supplemental Materials Section B. The mathematical framework used in this study treats the leaﬂets as exactly incompressible. Thus, within our numerical framework, the volumetric part of the strain energy, U(J)=(JlnJJ+ 1), (6) can be viewed as a stabilization term. For further details, see the summary of the study of Vadala-Roth et al. [45] in Supplemental Materials Section E. We use = 14.1 MPa in our simulations. The geometry of the Labcor TLBP A Supra porcine aortic valve is constructed based on literature values [46]. A model collagen ﬁber architecture is created using Poisson interpolation [47] (Fig. 2a). The pericardial valve geometry is reconstructed from micro-CT images of a Carpentier–Edwards PERIMOUNT RSR Model 2800 surgical aortic heart valve (Fig. 2b). The mean ﬁber orientation is chosen to be 45(Fig. 2b), following the small angle light scattering data of Sun et al. [28] 2.2.2 Aortic Test Section The wall of the aortic test section is glass. To avoid the expensive linear solvers required by exactly imposing the rigidity constraint within the present computational framework [48], we instead use a penalty method that models the test section as a sti↵neo-Hookean material with [21] Wwall=cwall 2(¯I13). (7) 5 (a) (b) Figure 2: (a) Model porcine bioprosthetic valve geometry and ﬁber architecture. The idealized geometry of the Labcor TLBP A Supra porcine aortic valve is reconstructed based on literature values [46]. The model ﬁber structure is generated using Poisson interpolation [47]. (b) Model bovine pericardial bioprosthetic valve geometry and ﬁber architecture. This valve geometry is obtained from micro-CT imaging of a Carpentier- Edwards PERIMOUNT RSR Model 2800 surgical aortic heart valve. The model ﬁber structure is generated based on the small angle light scattering (SALS) data of Sun et al. [28] The SALS data show that the mean ﬁber orientation of a bovine pericardial valve leaﬂet is 45. Additional structural forces are included in the rigid test section model, F(X,t)=wall(X(X,t)). (8) In the limit as cwallandwallbecome large, the test section becomes e↵ectively rigid and stationary. We use cwall= 33.1 kPa and wall=8.52⇥105kPa cm2. 2.3 Fluid Model and Boundary Conditions We use the incompressible Navier–Stokes equations to model the test ﬂuid in the aortic test section of the pulse duplicator as a viscous incompressible ﬂuid. We model the saline solution with a uniform density ⇢=1.0gc m3and a uniform dynamic viscosity µ=1.0 cP (saline at 25C) for both cases. Results from additional simulations using a blood analogue ﬂuid with density 1 .0gc m3and viscosity 3 .5 cP are provided in Supplemental Materials Section G. These simulations indicate that the large-scale ﬂow structures and leaﬂet kinematics are similar between saline and the blood analogue. 6 (a) (b) Figure 3: Three-dimensional aortic test section models for the porcine (a) and bovine pericardial (b) BHV simulations along with the reduced-order models that provide driving and loading conditions. Three-element Windkessel models are used at the downstream (outlet) for both cases. (a) A three-element Windkessel model is used at the upstream (inlet) for the porcine aortic valve simulations. The pump pressure is derived from pressure and ﬂow data from the ventricular outﬂow tract of the pulse duplicator. (b) Because the pump ﬂow waveform and atrial pressure data are available for the bovine pericardial valve experiments, a more detailed pump model is used upstream for the bovine pericardial valve simulations. Three-element Windkessel (R–C–R) models establish downstream loading conditions for the aortic test section for both cases (Fig. 3). A three-element Windkessel model is also used for the porcine BHV simu- lations to capture the upstream driving conditions for the aortic test section (Fig. 3a). Additional data are 7 available for the bovine pericardial case, including the experimental pump ﬂow and atrial pressure wave- forms, which allow for a more complete description of the upstream components of the system (Fig. 3b). In both cases, we impose a combination of normal traction and zero tangential velocity boundary conditions at the inlet and outlet of the FSI model to couple the reduced-order models to the detailed description of the ﬂow within the aortic test section. The values of the resistances and compliance for the upstream model are CVIA=0.1 mmHg mL1,R1=0.15 mmHg mL1s, and R2=0.15 mmHg mL1s for the porcine aortic valve case, which characterize compliance and resistance of the VIA system. The values for the bovine pericardial valve case are CVIA 1=0.0275 mmHg mL1,CVIA 2=0.0347 mmHg mL1, RVIA=0.15 mmHg mL1s, and Rout=0.0898 mmHg mL1s. The mitral valve is modeled as a diode, with a resistance of RMV=0.0280 mmHg mL1s when the valve is open. The values for the downstream model are Rc=0.0218 mmHg mL1s,Rp=1.31 mmHg mL1s, and C=0.915 mmHg mL1for the porcine BHV, and Rc=0.0282 mmHg mL1s,Rp=1.22 mmHg mL1s, and C=1.27 mmHg mL1for the bovine pericardial BHV. Supplemental Materials Section C provides additional details. We use solid wall (zero velocity) boundary conditions on the remainder of the boundaries in the computational domain of the FSI model. 2.4 Fluid–Structure Interaction We use a ﬁnite element (FE) extension of the IB method [49] to simulate FSI. Supplemental Materials Section D provides details on the continuum formulation and numerical approximations. Reduced-order boundary condition models are coupled to the FSI model following the approach detailed by Grith et al. [16] 2.5 Numerical Discretizations Simulations are performed using a geometrically realistic model of the aortic test section of the pulse du- plicator with length 10.1 cm and diameter 28 mm (see Fig. 3). The aortic test section is embedded in a square computational domain with side lengths of 10.1 cm. In our IB formulation, Eulerian variables are approximated using a block-structured adaptively reﬁned Cartesian grid, and Lagrangian variables are ap- proximated using an unstructured FE mesh that conforms to the geometry of the structure [49]. The e↵ective ﬁne-grid resolution of the Cartesian grid is approximately 0 .39 mm. The structural meshes use second-order hexahedral (27-point) elements for the BHV leaﬂets and ﬁrst-order tetrahedral (4-point) elements for the aortic test section. The aortic test section and porcine BHV meshes were generated using Trelis (Compu- tational Simulation Software, LLC, American Fork, UT, USA), and the bovine pericardial BHV mesh was generated using Pointwise (Pointwise, Inc., Ft. Worth, TX, USA). The average grid-spacings of the meshes are 0 .4 mm for the aortic test section, 0 .59 mm for the porcine aortic valve, and 0 .75 mm for the bovine pericardial valve. We use a piecewise-linear kernel for the aortic test section and a three-point B-spline kernel for the valve leaﬂets as regularized delta functions (see Supplemental Materials Section D). The time step size starts at  t=7.5⇥106s, and it is systematically reduced if needed to avoid instabilities related to our time stepping scheme. The penalty parameters cwallandwallin Eqs. (7) and (8) detailed in the “Aortic Test Section” section are empirically determined as approximately the largest values allowed by our explicit time stepping algorithm at the time step sizes used in the simulations. Supplemental Materials Section F shows results from a grid convergence study to quantify the level of consistency in our numerical results and to justify the chosen level of grid resolution for the ﬁnal results reported here. 2.6 Software Infrastructure FSI simulations use the IBAMR software infrastructure, which is a distributed-memory parallel implemen- tation of the IB method with adaptive mesh reﬁnement (AMR) [50, 51]. IBAMR uses SAMRAI [52] for Cartesian grid discretization management, libMesh [53] for ﬁnite element discretization management, and PETSc [54] for linear solver infrastructure. 8 0 0.1 0.2 0.3 0.4 0.5 0.6time (s)050100150pressure (mmHg)PLV (simulation)PLV (experiment)PAo (simulation)PAo (experiment) (a)Porcine aortic valve0 0.1 0.2 0.3 0.4 0.5 0.6time (s)-5051015202530flow rate (l/min)QAo (simulation)QAo (experiment) (b)Porcine aortic valve 0 0.1 0.2 0.3 0.4 0.5 0.6time (s)050100150pressure (mmHg)PLV (simulation)PLV (experiment)PAo (simulation)PAo (experiment) (c)Bovine pericardial valve0 0.1 0.2 0.3 0.4 0.5 0.6time (s)-5051015202530flow rate (l/min)QAo (simulation)QAo (experiment) (d)Bovine pericardial valve Figure 4: Comparisons between simulated and experimental pressure and ﬂow rate waveforms for the porcine aortic valve (a)-(b) and bovine pericardial valve (c)-(d). The experimental waveforms shown are the average waveforms over 10 consecutive cycles of data, with shaded regions showing where 95% of the data fall. The experimental and computational stroke volumes for the porcine aortic valve are 69 .4±0.4m L and 72 .7 mL, respectively, and 71 .6±0.7 mL and 72 .1 mL for the bovine pericardial valve. The maximum experimental pressure di↵erences during forward ﬂow for the porcine aortic and bovine pericardial valve are 22.8±0.2 mmHg and 19 .7±0.5 mmHg, respectively. The maximum computational pressure di↵erences during forward ﬂow are 22 .4 mmHg and 16 .4 mmHg, respectively. 3 Results We perform corresponding experiments and simulations using a pulse rate of 70 beats per minute in all cases. We use 10 consecutive cycles of experimental pressure and ﬂow rate waveforms for both valves to characterize the reduced-order models that provide driving and loading conditions for the three-dimensional FSI models. The stroke volumes of the average ﬂow waveforms are 69 .4±0.4 mL and 71 .6±0.7 mL for the porcine and bovine pericardial BHVs, respectively. We assess the computational results by comparisons to available experimental data, including ﬂow rates, upstream and downstream pressures, and leaﬂet kinematics. The computational models also provide detailed ﬂow patterns and leaﬂet stress distributions, which are not readily available in the present experimental models. Figure 4 shows comparisons between simulated and experimental pressures and ﬂow rates for both valves, which are in good agreement. To quantify this, we calculate the discrepancy between the simulation data 9 Table 1: Comparisons of normalized L2- and L1-norms of discrepancies in the bulk measurements between simulation and experiment* porcine ( L2) pericardial ( L2) porcine ( L1) pericardial ( L1) PLV 3.4% 4.6% 8.6% 14.4% PAo 1.9% 1.4% 4.5% 3.6% QAo 4.7% 8.6% 9.5% 20.2% *Comparisons of L2-a n d L1-norms of the discrepancies in the bulk measurements from simulation and experiment shown in Figure 4, normalized by the L2-a n d L1-norms of the measurements from the experiment. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5time (s)00.511.522.5PDVA (cm2)simulationexperiment (a)Porcine aortic valve0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5time (s)00.511.522.5PDVA (cm2)simulationexperiment (b)Bovine pericardial valve Figure 5: Comparisons between simulated and experimental projected dynamic valve area (PDVA) for the porcine aortic (a) and bovine pericardial (b) valves. (a) The experimental data (acquired using custom apparatus depicted in Fig. 1b) are manually aligned with the beginning of the valve opening with the simulation data. The experimental PDVA measurement shown is the average PDVA over 10 consecutive cycles of data with shaded region showing where 95% of the data fall. (b) The experimental data are acquired using a high-speed videographic method, from which we reconstruct PDVA data via automatic image analysis via DataTank. Videographic data are available for only a single cycle for the pericardial BHV. and the mean experimental data by Mq=MsimulationMexperiment Lq(0,T) kMexperiment kLq(0,T),q =2,1, (9) in which (0 ,T) indicates an integral over time and Mexperimentis the experimental data averaged over 10 consecutive cycles. Comparisons are shown in Table 1. The computational stroke volumes are 72.7 mL and 72.1 mL for the porcine and bovine pericardial BHVs, which are 4 .68% and 0 .71% larger than the mean experimental stroke volumes, respectively. The maximum experimental pressure di↵erences during forward ﬂow for the porcine aortic and bovine pericardial valves are 22 .8±0.2 mmHg and 19 .7±0.5 mmHg, respectively. The maximum computational pressure di↵erences during forward ﬂow are 22 .4 mmHg and 16 .4 mmHg, which correspond to di↵erences of 2 .0% and 16 .5% compared to the experimental pressure di↵erence in the mean pressure waveforms. Figure 5 uses PDVA to compare the leaﬂet kinematics in the experimental and computational models. The experimental data are systematically shifted in time to align with the simulation data at the beginning of valve opening. The timing of opening and closing are in excellent agreement, and the open area is also in 10 Figure 6: Detailed comparison of the bovine pericardial valve leaﬂet kinematics during closure from simu- lation (top) and experiment (bottom). The simulation captures the behavior of each of the leaﬂets closing one at a time (order: bottom !right !left leaﬂet) as observed in the experiment. The time increment between frames for simulation is 9.6 ms. excellent agreement. There are some di↵erences in the ﬂuttering frequency, and we discuss potential sources of these discrepancies in the “Discussion” section. Figure 6 compares the leaﬂet kinematics of the bovine pericardial valve during closure in our simulation to images acquired via high-speed videography. We observe that each of the leaﬂets close one at a time in the experiment as well as in the simulation. Figure 7 provides detailed ﬂow patterns generated in the FSI simulations of the porcine aortic (Fig. 7a) and bovine pericardial (Fig. 7b) BHVs. We observe that the large-scale ﬂow features are similar in both cases. We classify the ﬂow regime of our simulations by computing the peak Reynolds number Repeak, Repeak=⇢QpeakD µA, (10) in which ⇢is the density, Qpeakis the peak ﬂow rate, Dis the diameter of the aortic test section, µis the dynamic viscosity of the ﬂuid, and Ais the cross-sectional area of the aortic test section. Repeakis 20,576 and 19,330 for the porcine tissue and bovine pericardial BHV, respectively. We also compare the leaﬂet kinematics of the porcine tissue (Figs. 8a and 9a) and bovine pericardial (Figs. 8b and 9b) valves. One important di↵erence between the two valves is the symmetry breaking in the pericardial leaﬂets, especially during closure. Because of the asymmetric ﬁber architecture in each of its leaﬂets, the bovine pericardial valve shows a swirling motion during closure, as also observed in real pericardial BHVs [55]. Di↵erences between the porcine valve and bovine pericardial valve resulting from di↵erent ﬁber architec- tures are also clear in the stress distributions (Fig. 10). The von Mises stresses on the valves are distributed according to the anisotropic material responses. The porcine valve shows a stress distribution that is aligned from commissure to commissure when the valve is loaded in diastole (Fig. 10a). The stress is also distributed symmetrically on each porcine leaﬂet during valve opening (Fig. 10c), which is again concentrated around the commissures. In contrast, the bovine pericardial valve shows an asymmetric stress distribution, and the von Mises stresses are concentrated at one of each leaﬂet’s commissures. The location where stress is the highest on the leaﬂet occurs where the ﬁbers collect at the commissures (Fig. 10b). The opposite is true during valve opening, when the highest leaﬂet stresses occur near the free edges of the leaﬂets (Fig. 10d). These patterns appear because the parts of the leaﬂets with less commissural support experience larger deformations. 4 Discussion This study has developed FSI models of BHVs in an experimental pulse-duplicator platform. Our leaﬂet models include anisotropic descriptions of the leaﬂet biomechanics that are based on experimental tensile test data. In addition, we use experimental pressure and ﬂow data obtained from the pulse-duplicator systems 11 (a)Porcine aortic valve (b)Bovine pericardial valveFigure 7:Cross-section view of simulated ﬂow patterns using the porcine aortic (a) and bovine pericar-dial (b) valve models. The color shows the axial velocity through the aortic test section at the center plane,with red indicating forward ﬂow and blue indicating reverse ﬂow. (a)Repeak= 20,576. (b)Repeak= 19,330.The time increment between frames is 57.6 ms.to establish realistic boundary models for the detailed FSI models. These boundary models are calibrated inisolation from the rest of the system, independent of the FSI model of the valves. The simulated pressures,ﬂow rates, and leaﬂet kinematics all emerge from integrating these three model components, and the motionof the leaﬂets, including the timing of valve opening and closing, is not prescribed. Further, because the ﬂowrate is not imposed in the model, and because the time-dependent conﬁguration of the valve determines theresistance of the aortic test section, achieving both pressures and ﬂow rates that are in good agreement withthe experimental data represents a nontrivial test of the model.The numerical results detailed in Figure 4 show excellent agreement in ﬂow rates, pressures, valve openareas, and the timing of valve opening and closure between the simulations and experiments. For instance, theexperimental ﬂow rate oscillations that are present during valve closure in the porcine aortic valve (Fig. 4b)are physical because they come from the interaction between the momentum of the ﬂuid and the complianceof the system, and these oscillations are captured by the computational model. We observe from the bovinepericardial BHV case that a more comprehensive reduced-order model of the upstream components of thepulse-duplicator system, including the pump, VIA, left ventricle, left atrium, and mitral valve, allows usto capture system dynamics, including closure (Fig. 4d), more completely. Table 1 shows theL2- andL1-12 (a)Porcine aortic valve (b)Bovine pericardial valveFigure 8:Leaﬂet kinematics of the porcine aortic (a) and bovine pericardial (b) valves during opening. (a)The time increment between frames is 1.92 ms. (b) The time increment between frames is 3.84 ms.norms of the discrepancies in the simulation results relative to the experimental measurements. In particular,theL2-di↵erences demonstrate good quantitative agreement, in which the relative discrepancies are within4.7% for the porcine aortic valve and 8.6% for the bovine pericardial valve. These results show that ourvalve models generate leaﬂet kinematics that yield realistic pressures and ﬂow rates in the fully coupled FSImodels.The experimental and computational leaﬂet kinematics are assessed using PDVA (Figs. 5a and 5b). ThePDVA measurements for the porcine aortic valve indicate reasonable agreement between simulations andexperiments. There is a key di↵erence in the techniques used to acquire these data. An electro-optical system(Fig. 1b) was used to measure PDVA directly for the porcine aortic valve. For the bovine pericardial valve,PDVA was calculated indirectly using automatic image analysis via DataTank from videographic images.An advantage of the electro-optical subsystem is that it is calibrated with a reference area at the samelocation as the valve in the test section, but high-speed video is not available for this system. In contrast,although high-speed video is available for the pericardial BHV as shown in Figure 6, the image analysismethod uses an estimated area-to-pixel scaling. In addition, the frame rate (400 fps) of the high-speedvideo of bovine pericardial valve is not sucient to capture the full dynamic response. Consequently, weare missing data needed to resolve the full dynamic waveform. We plan to use higher speed (5000 fps)videography in future work to quantify leaﬂet motion more completely. An important limitation of the presentmodel is that we do observe discrepancies between the computational and experimental leaﬂet ﬂutteringfrequencies and amplitudes. BHVs are known to be viscoelastic [56], and fully viscoelastic models may beneeded to achieve better agreement between simulation and experiment. At present, however, experimentallyconstrained viscoelastic models of BHV biomaterials suitable for three-dimensional mechanical analysesappear to be lacking. The structural models could be improved by incorporating additional experimentaldata that characterize the ﬂexural properties of the valve leaﬂets.Detailed comparisons of the leaﬂet kinematics for the bovine pericardial BHV showed that the leaﬂets13 (a)Porcine aortic valve (b)Bovine pericardial valveFigure 9:Leaﬂet kinematics of the porcine aortic (a) and bovine pericardial (b) valves during closure. (a)The time increment between frames is 9.6 ms. (b) The time increment between frames is 19.2 ms.closed one at a time both in the experiment and simulation (Fig. 6). This result indicates that the overallkinematics are in reasonable agreement between the experiment and our model. We speculate that this maybe because each of the leaﬂets has slightly di↵erent size and geometric features, along with the contributionfrom its ﬁber structure. Alternatively, ﬂow instabilities may induce the sequential closure of the valve leaﬂets.The large-scale ﬂow features for the porcine aortic and bovine pericardial valves are similar (Fig. 7).We also observe additional small-scale turbulent ﬂow features at current spatial resolutions, which are alsopresent in physiological ﬂow regimes [57].Repeakfor both cases is approximately 20,000, which clearlymotivates the need for further studies on the treatment of turbulence in these types of models. In thepresent study, we perform implicit large-eddy simulation (ILES) using high-resolution slope limiters, basedon the piecewise parabolic method (PPM) [58–61], to model the ﬂow ﬁeld. Explicit LES methods havenot yet been completely developed for the present IB approach to FSI. We plan to compare ILES [58–61]and explicit LES [62, 63] models for cardiovascular ﬂows in future work. The similarity in the large-scaleﬂow features for the two BHVs suggests that it is important to consider the leaﬂet kinematics in additionto the ﬂow patterns in comparing di↵erent BHVs. Studying leaﬂet kinematics could also be important inidentifying the factors that a↵ect the durability of di↵erent BHVs.Di↵erences in the porcine and bovine pericardial leaﬂet kinematics are particularly prominent duringleaﬂet closure (Figs. 8 and 9). Speciﬁcally, the twisting motion during closure that we see in the bovinepericardial valve results from the asymmetric alignment of its ﬁbers. This is also evident in stress analysesof the leaﬂets. Di↵erences in stress distributions (Fig. 10) are caused by the reduced strain at one of twocommissure points of each of the bovine pericardial leaﬂets, whereas the stress distributions in the porcineBHV model reﬂects the symmetric ﬁber architecture of the porcine valve cusps. Also, unlike the porcineaortic valve, the pericardial valve lacks nodes of Arantius, which are intrinsic anatomic structures in nativeaortic valve leaﬂets. It has been hypothesized that the nodes may play a role in distributing and equalizingleaﬂet closing stresses [64]. We also observe that the porcine BHV leaﬂets experience higher stresses, which14 (a)Porcine aortic valve (diastole) (b)Bovine pericardial valve (diastole) (c)Porcine aortic valve (systole) (d)Bovine pericardial valve (systole) Figure 10: von Mises stress (kPa) on the porcine aortic (a), (c) and bovine pericardial (b), (d) valves during diastole (a)-(b) and systole (c)-(d). The time increment between frames in panels (a)-(b) is 30.72 ms, and the time increment between frames in panels (c)-(d) is 11.52 ms. may also explain why porcine BHVs have a greater tendency to develop leaﬂet tears with regurgitation compared to bovine pericardial BHV [65, 66]. The stress concentrations at the commissures of the bovine pericardial BHV also agrees with the known failure regions for bovine pericardial BHVs [67,68]. The present study has several limitations. Although we report 10 consecutive cycles of experimental data, one future aim is to better quantify the uncertainty associated with cycle-to-cycle variability in the experi- ments, and to integrate this uncertainty into our computational models. We shall also quantify variability in the leaﬂet kinematics along with variability in the pressure and ﬂow rate measurements by collecting high-speed video and data over multiple cycles. We also have not yet systematically performed sensitivity analysis or uncertainty quantiﬁcation of the model. Performing such analyses is crucial for establishing the credibility of the model but is challenging because of the substantial computational requirements of such analyses. We also plan to validate the ﬂow ﬁelds by particle image velocimetry (PIV) under di↵erent con- ditions, as well as using a blood analogue (dynamic viscosity of 3–4 cP) as our test ﬂuid. With PIV, we also plan to perform direct comparison between computational and experimental ﬂows in steady state cases in which we expect to obtain converged mean ﬂows. This will help us quantify the accuracy our model by comparing time-averaged ﬂows as well as turbulence kinetic energies to compare both large- and small-scale ﬂow features to those of the experiments. In this study, we also use homogeneous structural models that 15 omit descriptions of the discrete layers of the valve leaﬂets, which will a↵ect how the leaﬂets deform [69,70]. Therefore we shall also further validate speciﬁcally the leaﬂet kinematics by focusing on quasi-static BHV deformations, as done by Sun et al. [28], or by comparing the leaﬂet kinematics with reconstructions of leaﬂet deformations from in vitro experiments, as done by Iyengar et al. [27] or Sugimoto et al. [71] Other potential validation tests include assessing BHV performance using additional hydrodynamics performance measures (pressure drop, e↵ective oriﬁce area, etc.) [72–74]. In summary, this study describes work to model BHVs in an experimental pulse-duplicator platform that is used in various settings to assess their performance. Ultimately, fully validated FSI models of BHVs may be used to develop a high-ﬁdelity model for predicting dynamic performance amongst di↵erent valve designs and addressing persistent challenges posed by current BHV designs. They could also facilitate matching patient requirements with valve performance speciﬁcs and the development of regulatory guidelines for evaluation of novel designs. This methodology may also be extended and further validated to be applicable to study the e↵ectiveness of TAVR devices with varying degrees of intra- and paravalvular leak, reduced leaﬂet mobility, subclinical valve thrombosis, and potential pannus formation [75]. Acknowledgment J.H.L. acknowledges funding from the NIH Integrative Vascular Biology Training Fellowship (5T32HL069768- 17) at the University of North Carolina School of Medicine. B.E.G. acknowledges funding from NSF (Awards OAC 1450327, OAC 1652541, and CBET 1757193) and NIH (Awards R01HL117063 and U01HL143336). Thanks to B. Vadala-Roth for providing data for Supplemental Materials Section E, and to K. Aycock and T. Morrison for reviewing the manuscript. The ﬁndings and conclusions in this article have not been formally disseminated by the U.S. FDA and should not be construed to represent any agency determination or policy. The mention of commercial products, their sources, or their use in connection with material reported herein is not to be construed as either an actual or implied endorsement of such products by the Department of Health and Human Services. Conﬂict of interest No conﬂict of interest.