Comparison of Immersed Boundary Simulations of Heart Valve Hemodynamics against In Vitro 4D Flow MRI Data

The experimental data in this study were obtained using the 4D flow MRI sequence developed by Markl et al. [28], which produced full 3D, three component, phase-averaged velocity fields in a physical model of the right ventricular outflow tract (RVOT) and pulmonary arteries. The experiments were performed in a whole body, 3 Tesla General Electric MRI scanner at the Richard M. Lucas Center for Imaging at Stanford University (Figure 1), and these data were discussed in a previous publication [39]. The gating signal for the 4D flow sequence was provided by the data acquisition system that was also used to control the ventricle pump. The cardiac cycle duration was 0.8325 s or approximately 72 beats per minute. Based on a retrospective gating scheme, the velocity data were acquired with a temporal resolution of approximately 83 ms, providing 10 time phases for the final data. The velocity at a point in the flow domain for a given phase represents the mean velocity during the time interval of that phase. Therefore, the velocity for each phase represents the phase-average of the velocity during an 83 ms interval during the cardiac cycle. This phase average is constructed from the flow fields of over hundreds of cardiac cycles and is therefore expected to smooth out small variations between cycles. The 4D flow MRI scans used 0.9 mm thick sagittal slices and spatial resolutions of 0.894 mm and 0.897 mm in the coronal and axial directions, respectively.

Multiple 4D flow MRI scans were conducted and averaged to increase signal strength and reduce noise. Additionally, pump-off scans were averaged and subtracted from the pump-on scans to reduce various sources of error, including eddy current effects [7]. To account for drift in the eddy currents and other signal properties, pump-on and pump-off scans were alternated. All of the pump-on scans, with the pump-off data already subtracted, were then averaged to comprise the final dataset. Four pump-on scans and three pump-off scans were measured, with each individual scan taking approximately 20 minutes. This 4D flow MRI procedure results in time-resolved, three-component, phase-averaged velocity data over the entire 3D test section volume. The following references have further background on MRI [8, 33, 32].

The expected uncertainty of 4D flow MRI velocity measurements is given by

$$\delta_{u}=\frac{\sqrt{2}}{\pi}\frac{\mathrm{VENC}}{\mathrm{SNR}}\,,$$

(1)

where SNR is the signal to noise ratio and VENC is the maximum velocity that can be measured by the 4D flow MRI sequence, known as the velocity encoding value [8]. The VENC was set at 250 cm/s for all runs. The signal to noise ratio was calculated as the ratio of the signal in the flow region to the signal in the solid wall, as done by Banko et al. [3]. The SNR was 23.7, resulting in an uncertainty of ±4.7 cm/s throughout the domain.

Note that this uncertainty was calculated using the SNR for the entire experimental dataset, which was comprised of averages of multiple 4D flow MRI scans and many cardiac cycles over approximately four hours of total scan time. There were sources of variation during the experiment that could produce different flow conditions cycle to cycle, including drift in the kinematics of the centrifugal and ventricular pumps in the flow loop and variations in the air pressure pneumatically driving the flow. The flow loop was observed during the MRI sequence and small manual adjustments were made to maintain the flow rate of 3.5 L/min. We observed that the cycle-to-cycle variations in the volumetric flow rate were relatively small and did not substantially contribute to the overall uncertainty in the velocity.

The anatomic physical model used in this work was designed to analyze the flow fields local to bioprosthetic valves in a healthy RVOT anatomy. The geometry modeled a 11-13 year-old pediatric patient with cardiac output of 3.5 L/min [39]. The model extended from the outlet of the right ventricle (RV), through the main pulmonary artery (MPA) to the left pulmonary artery (LPA) and right pulmonary artery (RPA) branches after the first bifurcation. The healthy geometry design was based on measurements made from MRI data of the right heart and pulmonary arteries of six healthy subjects between ages 11 and 13. MRI data were collected for clinical purposes and used for modeling under an IRB-approved protocol. For each patient, we measured the diameters of the MPA, LPA, and RPA in the sagittal and axial planes. The measurement location was at the valve annulus while the locations for the LPA and RPA were immediately downstream of the bifurcation. We measured the turning angle, radius of curvature, and arc length of the RVOT and MPA in both sagittal and axial planes for each subject. The other key parameters were the three angles at the bifurcation: LPA to RPA, MPA to LPA, and MPA to RPA. Measurements were done in ImageJ (Bethesda, MD).

Median values for all subjects of the diameters, lengths, radii of curvature, and bifurcation angles were used to construct the healthy pulmonary geometry. Median values were chosen to reduce the effect of outliers in the patient measurements, though differences between median and mean were generally small. The model design was scaled to fit the 25mm bioprosthetic valve used in this experiment, resulting in an MPA diameter of 25mm, LPA diameter of 12mm, and RPA diameter of 14mm. All scaled values fell within the range of measured values from the patient cohort, resulting in a scaled model that accurately represented the 11-13 year-old patient population.

The model was designed in SolidWorks (Waltham, MA) and was manufactured using stereolithography (SLA) with DSM Somos WaterShed XC 11122 resin at the W.M. Keck Center at the University of Texas, El Paso. All components of the model and the valve itself are fully MR compatible. The 3D printed healthy RVOT model is shown in Figure 2a, with all components fully assembled.

The valve used in this experiment was a 25mm Epic valve (Abbott Cardiovascular, Plymouth, MN), which is a porcine bioprosthetic trileaflet aortic surgical valve, which we placed in the pulmonary position. The valve includes support scaffolding at the commissures and a sewing ring that is used to suture the valve to the RVOT vessel in the patient. The valve was implanted into the model using a two-piece component valve-holder piece with a groove that matched the dimensions of the sewing ring on the 25mm Epic valve, which has an estimated internal diameter of 20 mm, as measured via a caliper. This groove was clamped around the sewing ring with the two pieces screwed together to secure the valve (Figure 2b).

We developed an in vitro flow loop experimental setup that replicates physiological conditions over the cardiac cycle for the pediatric pulmonary system (Figure 3). The pulsatile flow loop was driven by a custom designed ventricle box that was inspired by pulsatile ventricular assist devices, such as the Berlin Heart EXCOR (The Woodlands, TX). A thin membrane in the ventricle box was pneumatically driven by a pulsatile air supply to create systole and diastole in the model and the working fluid. The digital trigger signal governing the pneumatic ventricle box was connected to the electrocardiogram (ECG) converter and trigger on the MRI system to signal the start of each cardiac cycle, allowing us to collect gated, phase-locked data.

The ventricle box was the key component of the RVOT model that drives the full flow loop to replicate right heart circulation and pulmonary physiology (Figure 3). The flow loop included an upstream inlet tank fed by a steady centrifugal pump from a supply reservoir that connected to a flexible bag that acted as a right atrium. This atrium proxy fed the ventricle box through a 27mm bioprosthetic tricuspid valve. The ventricle box connected to the main printed components of the healthy RVOT geometry, including the component that housed the bioprosthetic valve. The outlets of the 3D printed model connected to flexible tubing, which in turn connected to capacitor elements. Each capacitor element exited into another segment of flexible tubing which ran to the outlet tank of the flow loop.

All of the experimental flow loop components were tuned so that the pressure and flow rate waveforms reasonably matched physiological conditions. The pressures were measured with Millar (Houston, Texas) SPR-350 pressure transducer catheters inserted in the model through sealed ports. The RV pressure was measured directly upstream of the valve while the MPA pressure was taken immediately downstream of the valve; both catheters were centered in the vessel. The branch flow measurement was collected using a Transonic Systems (Ithaca, NY) ultrasonic flow probe placed around each branch PA. The locations of these flow probes, and the one placed upstream of the model to monitor the cardiac output are seen in Figure 3. The pressure waveform was averaged over approximately ten cardiac cycles.

The working fluid for all experiments was a blood analog consisting of 60% water and 40% glycerin, with a density of 1.09 g/cm³ and a viscosity of 3.9 centipoise at room temperature in the magnet room. Gadolinium was added to the blood-analog fluid to increase the signal contrast during the 4D flow MRI scan. We conducted the experiment at a heart rate of 72 beats per minute and targeted a 50-50 split for flow in the LPA and RPA. Additional details about the development of the physiological flow loop and experimental model and methods can be found in Schiavone et al. [39].

We constructed the model valves using a design-based approach to elasticity, which we introduced in our previous studies [16, 17]. With this method, we specified that the valve leaflets must, via tension in the leaflets, support a pressure. From this specification, we derived and solved a nonlinear partial differential equation for mechanical equilibrium of the leaflets under the prescribed load. By tuning free parameters and boundary conditions for these differential equations, we designed a closed configuration of the valve, including the tensions in the loaded configuration. From this closed configuration, we then derived the reference configuration and material properties of the model valve. This model was fiber based, in that the structure was discretized as a system of one-dimensional curves that fill a region in space. Alternatively, the model may be viewed as a system of nonlinear springs that includes a continuum limit. These methods were found to be highly effective under physiological pressures, opening and closing repeatedly over multiple cardiac cycles and achieving realistic flow rates with physiological driving pressures.

To design the model valve, we represented the leaflet as an unknown parametric surface

$\displaystyle\mathbf{X}(u,v):\Omega\subset\mathbb{R}^{2}\to\mathbb{R}^{3}.$

(2)

We assumed that the parameters conform to the fiber and cross fiber directions of the leaflets, meaning that for curves of constant $v$ on which $u$ varies represent the fibers, and curves of constant $u$ on which $v$ varies represent curves in the cross-fiber direction. Let single bars $|\cdot|$ denote the Euclidean norm. We denoted the unit tangents to the surface as

$\displaystyle\frac{\mathbf{X}_{u}}{|\mathbf{X}_{u}|}\quad\text{ and }\quad\frac{\mathbf{X}_{v}}{|\mathbf{X}_{v}|},$

(3)

which represented the local fiber and cross-fiber directions, respectively. These directions were not prescribed or required to be orthogonal. We denoted the magnitude of local membrane tensions in the circumferential and radial directions as $S$ and $T$, respectively.

We then prescribed a uniform pressure load of $p=30$ mmHg, the approximate diastolic pressure load in the experimental data immediately following the closing transient. We then considered the mechanical equilibrium on an arbitrary patch of leaflet $[u_{1},u_{2}]\times[v_{1},v_{2}]$. This requirement gives the following integral balance of pressure and leaflet tensions:

	$\displaystyle 0=\int_{v_{1}}^{v_{2}}\int_{u_{1}}^{u_{2}}p\left(\mathbf{X}_{u}(u,v)\times\mathbf{X}_{v}(u,v)\right)dudv$		(4)
	$\displaystyle+\int_{v_{1}}^{v_{2}}\left(S(u_{2},v)\frac{\mathbf{X}_{u}(u_{2},v)}{|\mathbf{X}_{u}(u_{2},v)|}-S(u_{1},v)\frac{\mathbf{X}_{u}(u_{1},v)}{|\mathbf{X}_{u}(u_{1},v)|}\right)dv$
	$\displaystyle+\int_{u_{1}}^{u_{2}}\left(T(u,v_{2})\frac{\mathbf{X}_{v}(u,v_{2})}{|\mathbf{X}_{v}(u,v_{2})|}-T(u,v_{1})\frac{\mathbf{X}_{v}(u,v_{1})}{|\mathbf{X}_{v}(u,v_{1})|}\right)du.$

Then, we applied the fundamental theorem of calculus to convert all single integrals to double integrals over the patch and swapped the order of integration formally as needed to obtain

$\displaystyle 0=$

$\displaystyle\int_{v_{1}}^{v_{2}}\int_{u_{1}}^{u_{2}}\bigg{(}p(\mathbf{X}_{u}\times\mathbf{X}_{v})+\frac{\partial}{\partial u}\left(S\frac{\mathbf{X}_{u}}{|\mathbf{X}_{u}|}\right)+\frac{\partial}{\partial v}\left(T\frac{\mathbf{X}_{v}}{|\mathbf{X}_{v}|}\right)\bigg{)}\;dudv.$

(5)

This equation represents the integrated form of the equations of equilibrium. The patch is arbitrary, and so, formally assuming sufficient smoothness, the integrand must be zero. This gave the following partial differential equation form of the equations of equilibrium:

$\displaystyle 0=p(\mathbf{X}_{u}\times\mathbf{X}_{v})+\frac{\partial}{\partial u}\left(S\frac{\mathbf{X}_{u}}{|\mathbf{X}_{u}|}\right)+\frac{\partial}{\partial v}\left(T\frac{\mathbf{X}_{v}}{|\mathbf{X}_{v}|}\right).$

(6)

As written, this system of differential equations is not closed, as it has five unknowns, the three components of $\mathbf{X}$ and the two tensions $S,T$ and only three equations, one for each component. To close the system, we temporarily specified that

$\displaystyle S(u,v)=\alpha\left(1-\frac{1}{1+|\mathbf{X}_{u}|^{2}/a^{2}}\right),\quad T(u,v)=\beta\left(1-\frac{1}{1+|\mathbf{X}_{v}|^{2}/b^{2}}\right).$

(7)

This functional form allowed the solution of the differential equation to find heterogeneous tensions to support the prescribed pressure load, while preventing extreme local tensions from occurring. The parameters $\alpha,\beta,a,b$ are tunable free parameters that we set to match the gross morphology of the prosthetic valve. The parameters $\alpha,\beta$ represent the maximum allowable tension in the circumferential and radial directions. The parameters $a,b$, which have units of cm, were tuned to control the gross morphology of the leaflets.

To tune these parameters and design a valve that corresponded to the prosthetic valve, we measured its gross morphology with a ruler and caliper. The free edge length of each leaflet was approximately 2.2 cm. We estimate the leaflet height at 1.3 cm, but this was highly dependent on loading; even pulling the leaflets flat along the ruler changed lengths substantially. The internal radius of leaflet attachment to the annulus was 1.0 cm and annular height was 1.1 cm. To match these values, we set the parameters $\alpha,\beta$ to $1.3\cdot 10^{7}$ and $4.6\cdot 10^{5}$ dynes, respectively. The value of $a$ varied linearly from $42$ cm at the annulus to $113$ cm at the free edge. Taking a variable value of $a$ was necessary to avoid excessive billowing of the leaflet near the annulus while maintaining adequate free edge length. The value of $b$ was set to $140$ cm.

We discretized the system of equations (6) including the tension formulas (7) with centered finite differences. The position of the valve at the annulus and commissures was prescribed as a Dirichlet boundary condition. At the free edge, we prescribe zero-tension (homogeneous Neumann) boundary conditions. The resulting nonlinear system of equations was solved with Newton’s method with line search. Additionally, we applied the method of continuations on pressure, in which pressure is adaptively increased from an initial value of zero to the prescribed value of $p=30$ mmHg, and the solution with the previous pressure is used as the initial guess for Newton’s method with the subsequent pressure.

We then used the solution to equations (6), which represents the predicted loaded configuration of the valve, to derive a reference configuration and constitutive law. For each link in the discretized model, the solution included the loaded length and the tension needed to support the prescribed pressure load. Based on the experiments of Yap et al. [46], we prescribed uniform stretch ratios of $\lambda_{c}=1.15$ circumferentially and $\lambda_{r}$ = 1.54 radially, then used this information to compute the reference length for each link from its loaded length. For a link with length $L$, we then solved $\lambda=L/R$ for the rest length $R$ at the stretch $\lambda$ corresponding to the direction of the link. We took the tension/stretch relationship to be exponential with a zero at $\lambda=1$, and based the shape, but not the local stiffness, on the strip biaxial tests of May-Newman et al. [29]. Based on a nonlinear least squares fit to their data, we took the exponential rate to be $\eta_{c}=57.46$ circumferentially and $\eta_{r}=22.40$ radially. For each link in the discretized, predicted loaded configuration with tension $\tau$ and the appropriate exponential rate and stretch for the circumferential or radial direction, we solved

$\displaystyle\tau=\kappa(e^{\eta(\lambda-1)}-1)$

(8)

for the local stiffness coefficient $\kappa$. Since the solution to equations (6) includes heterogeneous tensions, this resulted in heterogeneous material properties in the leaflets.

To generate an initial configuration suitable for FSI simulations, we sought a configuration of the leaflets that is open and free from external loading. Using the constitutive law we just set to determine tensions, we again solved the equilibrium equations (6) with $p=0$ mmHg. On the physical valve, a small portion of the leaflets point was attached radially inward from the commissures. To model this geometric feature, we force one eighth of the free edge points starting at each commissure of each leaflet to coincide (approximately 1/4 of each leaflet total), thus pinching the leaflets together. The remainder of the free edge was fixed as a Dirichlet boundary condition to ensure the leaflets do not self-intersect. This model has pre- or residual stretch and tension in this configuration, though substantially less than the predicted loaded configuration. Further, a configuration with zero tension, given the reference lengths that are computed locally for each link in the discretized model, may not exist.

Until this point, we used a membrane formulation of the leaflets. Next, we moved to a thickened formulation, while still using a fiber-based material model. To achieve a realistic thickness for the valve, we extruded the membrane for the valve normally in each direction, to form three adjacent membranes with a total thickness of $0.44$ mm, as reported in Sahasakul et al. [37]. This also served to mitigate the “grid aligned artifact” associated with large pressure differences across zero-thickness membranes in the IB method [16]. The stiffness of coefficients for each layer were set to be one third of the membrane stiffness coefficients computed previously. Linear springs of rest length $0.22$ mm were placed to keep the layers together through the simulation.

The gross morphology of the model valve that emerged from this process is shown in Figure 4. The free edge was 2.87 cm, corresponding to 3.3 cm in the predicted loaded configuration. After the pinching the leaflets together at the commissures, this left approximately 2.1 cm of free edge rest length per leaflet free to move independently of the other leaflets, within measurement error of ±0.1 cm from the free edge length of 2.2 cm measured on the prostheses. The leaflet rest height was 0.94 cm corresponding to a predicted loaded height of 1.44 cm. The measured leaflet height of 1.3 cm is nearly the predicted loaded height of 1.44 cm, which may be because the leaflets are so compliant in the radial direction, that pulling them flat to measure them achieved substantial stretches. The fiber orientation of the model runs from commissure to commissure and qualitatively matches experimental observations [38], though direct quantitative comparison is beyond the scope of this work. One minor limitation is that we do not add bending rigidity to the leaflets, beyond what emerges from the thickening process described above, and thus may not accurately capture leaflet flutter or other similar behaviors. Based on the thickness of $0.44$ mm, we estimated the mean tangent modulus at the predicted loaded stretches as $6.7\cdot 10^{7}$ dynes/cm² circumferentially and $8.5\cdot 10^{4}$ dynes/cm² radially. The prosthetic valve tissue is fixed in glutaraldehyde, and literature values for the fully-recruited circumferential tangent modulus of fixed porcine aortic valve tissue vary widely. Based on the experimental measurements of Billiar and Sacks and their constitutive law for valves fixed under 4 mmHg of pressure, we evaluated their constitutive law at the relevant stretches $\lambda_{r}$ and $\lambda_{c}$ and estimated the circumferential tangent modulus to be $4.0\cdot 10^{7}$ dynes/cm² [6]. Rousseau et al. reported moduli ranging from $5.1\cdot 10^{7}$ to $2.3\cdot 10^{8}$ dynes/cm², depending on the applied preload during fixation [36]. Sung et al. reported moduli ranging from $3.0\cdot 10^{8}$ to $5.2\cdot 10^{8}$ dynes/cm², depending on fixation pressure [43]. Thus our estimated tangent modulus falls within the range of existing studies, so we considered our resultant modulus in good agreement given the complexity of the steps involved, phenomenological nature of the constitutive law and uncertainties in experiments. We do not have access to the precise material properties of the prosthetic valve, and further, the only literature we could find on the material properties of a similar prostheses reported the tangent modulus at one particular loading, which did not appear to be at a relevant stretch for comparisons with our model [19]. Thus, our model has material properties in a reasonable range for a fixed aortic valve prostheses (placed in the pulmonary position in our simulations), but it does not directly model the material properties of the prostheses.

These steps completed the derivation and design of the model valve. We then used the model valve and constitutive law described throughout this section as initial conditions for fluid-structure interaction simulations.

We constructed the model vessel for FSI simulations from data from the MRI scans (Figure 4). The signal magnitude of 3D printed model material is distinct from the signal of the fluid in the scans, and we applied a thresholding operation to generate a three-dimensional model of the printed vessel surface. Using the MRI data ensured that the MRI and simulation coordinates were consistent in space and there were no alignment or registration errors. While using the files that generated the 3D printed model would have offered more spatial fidelity, the potential error in flow fields due to any mis-registration would have likely been much more substantial. Using Meshmixer (San Rafael, CA), we smoothed the mesh to remove stair-step effects and removed artifacts from the valve scaffold. We then remeshed to the desired edge length of 0.25 mm and extruded the model 0.25 mm and 0.5 mm to create a three-layer structure. As in the valve, this serves to eliminate the “grid aligned artifact” that can occur with pressure differences across thin membranes in the IB method [16]. Flow extenders of length 1 cm were added to the vessel at the inlet and both outlets to ensure that the normal to the vessel was aligned with the normal of the fluid box at the inlets and outlets. In FSI simulations, the vessel was held in place using target points, stiff springs of zero rest length that connect the current position of each model node to its desired position (Section 3.3). Additional linear springs are placed on each edge in the triangulated model. These springs are not meant to model a particular material and only serve to keep the vessel rigid and stationary throughout the simulation.

The valve leaflets are mounted to a scaffold, which we measured with a ruler and caliper as with the leaflets. The bottom of the scaffold has variable height with amplitude 0.2 cm, with the highest points aligned with the commissures, and we prescribed the bottom curve to take the functional form $0.2\left(\frac{1}{2}(\cos(3\theta)+1\right)$. The top of the scaffold has minimum height $h_{min}=0.65$ cm and maximum height $h_{max}=$1.5 cm aligned with the commissures, measured with the valve resting on a table. We prescribed the top of the scaffold to take the form $h_{min}+h_{max}|\cos\left(\frac{3}{2}\theta\right)|^{p}$, where $p$ is a power to be determined. We measured the scaffold height at two additional locations, then estimated $p=3.1$ with a least squares fit. From its minimum radius of 1 cm, the radius of leaflet attachment, we extended the scaffold 0.3 cm radially to form a subset of a cylinder. This thickness, slightly larger than the real scaffold thickness of 0.25 cm, ensured the scaffold overlapped with the vessel and prevented leaks between the vessel and scaffold. Finally, a small nub of thickness 0.15 cm and height 0.2 cm was placed inward radially at the commissure tip to model the bulging of the scaffold at that location. As with the vessel, the position was held nearly constant with target points and stiff linear springs connecting nodes, here placed with a cylindrical topology. We do not model the texture or irregularities of the real scaffold structure, but believe these measurements are accurate to approximately $\pm$1 mm throughout the model.

We used the immersed boundary method for fluid-structure interaction simulations [35]. The immersed boundary method uses two reference frames, an Eulerian or laboratory frame for fluid quantities and a Lagrangian or material frame for structure quantities. In one distinctive feature of the IB method, the fluid interacts with the structure through an Eulerian-frame body force, which in turn is computed from the Lagrangian frame structure force. Let $\rho,\mu$ denote the fluid density and viscosity, respectively. Let $\mathbf{x}$ represent a fixed location in space, and $t$ time. The fields $\mathbf{u},\mathbf{p}$ represent fluid density and pressure, and $\mathbf{f}$ denotes the IB body force in the Eulerian frame. Let $\mathbf{s}$ label a material point of the structure, which we previously denoted $u,v$ in Section 3.1 but switch to avoid confusion with fluid velocity, and $\mathbf{X}(\mathbf{s},t)$ denote the position of the structure point with label $\mathbf{s}$. The field $\mathbf{F}$ represents the structure force in the Lagrangian frame. Let $\delta$ denote the Dirac $\delta$-function.

The governing equations of the IB method are

	$\displaystyle\rho\left(\frac{\partial\mathbf{u}(\mathbf{x},t)}{\partial t}+\mathbf{u}(\mathbf{x},t)\cdot\nabla\mathbf{u}(\mathbf{x},t)\right)$	$\displaystyle=-\nabla p(\mathbf{x},t)+\mu\Delta\mathbf{u}(\mathbf{x},t)+\mathbf{f}(\mathbf{x},t)$		(9)
	$\displaystyle\nabla\cdot\mathbf{u}(\mathbf{x},t)$	$\displaystyle=0$		(10)
	$\displaystyle\mathbf{F}(\,\cdot\,,t)$	$\displaystyle=\mathcal{F}(\mathbf{X}(\,\cdot\,,t))$		(11)
	$\displaystyle\frac{\partial\mathbf{X}(\mathbf{s},t)}{\partial t}$	$\displaystyle=\mathbf{u}(\mathbf{X}(\mathbf{s},t),t)$		(12)
		$\displaystyle=\int\mathbf{u}(\mathbf{x},t)\delta(\mathbf{x}-\mathbf{X}(\mathbf{s},t))\;d\mathbf{x}$
	$\displaystyle\mathbf{f}(\mathbf{x},t)$	$\displaystyle=\int\mathbf{F}(\mathbf{s},t)\delta(\mathbf{x}-\mathbf{X}(\mathbf{s},t))\;d\mathbf{s}.$		(13)

The equations (9) and (10) are the Navier Stokes equations governing the dynamics of a viscous, incompressible fluid, plus the extra IB body force $\mathbf{f}$. Equation (11) represents the Lagrangian frame force via the mapping $\mathcal{F}$, which determines the field $\mathbf{F}$ as a function of the entire configuration of the structure. Equations (12) and (13) are the interaction equations, which couple the two frames. Equation (12) performs velocity interpolation and specifies that the structure moves with the local fluid velocity. The Lagrangian frame force is spread to the Eulerian grid via equation (13).

We ran simulations using the open-source solver IBAMR (Immersed Boundary Adaptive Mesh Refinement) using a staggered grid discretization of the fluid domain [12, 13]. Simulations were run on the Shrelock Cluster at Stanford University on four nodes or 96 Intel Xeon Gold 5118 cores with a 2.30GHz clock speed. The discrete $\delta$-function used was the five-point kernel, which we selected for its translational invariance [4] and lack of even-odd condition, which may create artifacts in the flow field when using a staggered-grid flow solver [15].

To drive the simulations, we prescribed pressures at the RV inlet and PA outlets based on the experimental measurements coupled to time-varying resistance boundary conditions. The motivation for this setup is as follows. In the physical experiment, an approximately equal flow split was achieved by tightening a clamp on downstream tubing. In the simulation, we wished to similarly achieve an approximately equal flow split, but the RPA and LPA portions of the domain had differing resistances due to their geometry, including curvature and length within the domain of interest. Further, the experimental pressure was measured approximately at the junction of these two vessels and pressures downstream at the simulation outlets were not available to use directly. To achieve an approximately even flow split on the two outlets, the total resistance of each outlet must be equal. Thus, we added a small resistance during systole to the less resistive outlet to make the sum of the resistance due to geometry and the resistance added to be equal on both outlets. During diastole, the valve became heavily loaded in a short duration of time, which created an oscillation in the flow waveform in simulations. On preliminary simulations without resistance this ringing persisted for an excessively long duration. The experimental pressure waveform includes an oscillation at the corresponding time, but if this oscillation was prescribed directly at the outlets extensive ringing occurred because the frequency of the experimental oscillation was not precisely equal to that of the simulation setup. (The time resolution on flow rates was ten frames per cardiac cycle, which was not enough to capture oscillations that may have occurred in the experimental flow rate.) Thus, we prescribed a resistance to both outlets, which allows the pressure and flow to oscillate in phase, which in turn caused the amplitude of the oscillations in both pressure and flow to decay over time. The experiment includes damping downstream of the region of interest that would be very complex to model. We selected the value of diastolic resistance by trial and error to introduce numerical damping in the system and reduce ringing in the flow rate. We expect that this resistance has minimal influence on the opening dynamics of the valve and almost none on the fully closed state, rather the resistance primarily influences the amplitude, decay rate and duration of the closing transient seen in the flow rate as the valve oscillates and settles into its closed state. We applied pressures, rather than prescribed flow rates, because pressures are necessary to load and properly close the valve during diastole. If we had applied a flow rate, during diastole the valve would not be loaded by the appropriate pressure difference to close in a physiological manner and would instead remain open under near-zero flow conditions. These strategies created reasonable agreement with the experimental pressures and flow rates throughout the cardiac cycle (minus small differences in the pressure oscillation frequency), while at the same time achieving a more even flow split and better damping of the closing transient than prescribing the experimental pressures alone.

To implement these strategies, we smoothed the experimental pressures measured on the RV and PA sides of the valve by taking a weighted average in time via convolution with a normalized cosine bump of radius 0.01 s. We then represented each curve as a finite Fourier series with 600 frequencies, which ensured the data were smooth and periodic in time. At the outlets, we applied a time-varying resistance boundary condition of the form

$\displaystyle R(t)Q=P_{outlet}-P_{exp}$

(14)

where $R(t)$ is the resistance, $Q$ is the instantaneous flow rate during the simulation, $P_{outlet}$ is the outlet pressure and $P_{exp}$ is the experimental pressures as processed above. To achieve an approximately even flow split during systole, we estimated the systolic resistances as follows. We ran a steady simulation for 0.2 s with the maximum experimental pressure difference of $p=15.66$ mmHg prescribed at the RV and zero pressure at both the RPA and LPA. This resulted in flow rates of $Q_{RPA}$ = 93.4 ml/s and $Q_{LPA}$ = 171.2 ml/s. We integrated the simulation pressure field at a slice in the PA junction to estimate a mean pressure of 3.53 mmHg. We then estimated that the resistance of the vessels themselves was 27.49 and 47.85 dynes/cm⁵ for the left and right PAs, respectively. The resistance at systole is set to $20.36$ s dynes/cm⁵ on the LPA and $0$ on the RPA, which makes the sum of the resistance in the vessel and the added outlet resistance approximately equal for both outlets. When the sign of the pressure difference across the vessel changed sign, indicating diastole, $R(t)$ smoothly increased to $1000.0$ s dynes/cm⁵. The resistance value remained constant through the remainder of diastole, then decreased smoothly to the systolic resistance after the pressure again changed sign, indicating systole.

The cardiac cycle duration was 0.8325 s, and simulations were run for two cycles. Results from the second cycle were analyzed and results from the first cycle was discarded due to initialization effects. (We subsequently report times in the second cardiac cycle modulo the cycle duration.) The time step was set to $\Delta t=7.5\cdot 10^{-6}$ s, which was the largest value we found to be numerically stable given the explicit time discretization of the structure and fluid-structure coupling. At peak flow, this time step resulted in an advective Courant Friedrichs Lewy (CFL) number of $\max\|\mathbf{u}\|\Delta t/\Delta x=0.036$. The observed stability limit is thus much less than that expected by the fluid velocity alone, and we do not consider lower time steps further. We output 1200 frames per second of the simulation for post processing. The fluid resolution was $\Delta x=0.45$ mm, exactly half that of the MRI resolution in the sagittal direction, and approximately half that of the MRI in the other two directions. (A convergence study was performed, see Appendix.) The structure resolution was targeted to half that of the fluid resolution, or 0.225 mm, but the precise lengths of each link in the model were determined by the steps in Section 3.1. We set $\rho=1.09$ g/cm³ and $\mu=$ 3.9 centipoise. The vessel was placed into a $9.36\times 6.48\times 5.76$ cm box, which corresponds to $208\times 144\times 128$ points. No local mesh refinement is used, which facilitates post-processing the velocity field (Section 3.4). The valve was placed at the origin with its axis aligned with the $x$ direction. The minimum height of the leaflet attachment is placed at -0.2 cm, just below the origin. Open-boundary stabilization was applied at the inlets and outlets [5]. Outside of the inlet and outlets, on faces that include the inlet and outlets, we prescribed no slip boundary conditions. On the remaining three faces, we prescribed zero pressure boundary conditions. No contact forces were needed, as the IB method automatically prevents contact between structures via the interaction equations and continuity of the velocity field [25].

We post-processed the simulation data to obtain a phase-averaged, resampled velocity field $\bar{\mathbf{u}}$ in a manner that mimicked experimental data acquisition. The scanner reported the phase-averaged velocity field at ten time steps per cardiac cycle, recorded from approximately ten cycles. To process simulation data, we computed the mean of the velocity field over the second cycle binned into ten time intervals of duration 0.083 s. (We also computed phase averages over three cycles but found no qualitative difference in the resulting flow fields and little quantitative difference between such results, see Appendix.) We converted phase-averaged velocity field from cell to point data, then linearly interpolated onto the nodes of the MRI mesh using PyVista [42].

We analyzed the flow quantitatively with the following metrics. We computed the $L^{p}$ relative error for $p=1,2$ of phase-averaged, resampled velocity magnitude,

$\displaystyle\frac{\|\bar{\mathbf{u}}-\mathbf{u}_{exp}\|_{L^{p}}}{\|\mathbf{u}_{exp}\|_{L^{p}}}=\frac{\left(\int|\bar{\mathbf{u}}-\mathbf{u}_{exp}|^{p}\;d\mathbf{x}\right)^{1/p}}{\left(\int|\mathbf{u}_{exp}|^{p}\;d\mathbf{x}\right)^{1/p}},$

(15)

and the $x$ (axial) component of the phase-averaged, resampled velocity only, which we denote as $\bar{u}$,

$\displaystyle\frac{\|\bar{u}-u_{exp}\|_{L^{p}}}{\|u_{exp}\|_{L^{p}}}=\frac{\left(\int|\bar{u}-u_{exp}|^{p}\;d\mathbf{x}\right)^{1/p}}{\left(\int|u_{exp}|^{p}\;d\mathbf{x}\right)^{1/p}}.$

(16)

Both relative errors were evaluated on the entire flow domain interior to the vessel and on three two-dimensional slices, $x=$ 0, 0.625 and 1.25 cm. We computed these metrics during the systolic phase only, because during diastole the flow fields are near zero in both simulation and experiment.

We also computed the integral metric $I_{1}$, which represents the nondimensional streamwise momentum and is computed on two dimensional slices of the domain [2, 39]. The metric is defined as

$\displaystyle I_{1}$

$\displaystyle=\left(\frac{1}{U_{T}^{2}A}\iint(\bar{\mathbf{u}}\cdot\mathbf{n})^{2}\;dA\right)^{1/2}=\left\|\frac{\bar{\mathbf{u}}\cdot\mathbf{n}}{U_{T}A^{1/2}}\right\|_{L^{2}}$

(17)

where $A$ is the area of the slice, $\mathbf{n}$ is the unit normal and $U_{T}$ is a velocity scale. The velocity scale is set to the maximum of the spatially averaged velocity over time,

$\displaystyle U_{T}=\max_{t\in[0,0.8325]}\left(\frac{1}{A}\iint\left(\bar{\mathbf{u}}(\mathbf{x},t)\cdot\mathbf{n}\right)\;dA\right).$

(18)

The velocity scale for the simulation and experiment were computed individually. This metric is 1 for a flow with a uniform velocity profile equal to $U_{T}$.