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

Abstract

The immersed boundary (IB) method is a mathematical framework for fluid–structure interaction problems (FSI) that was originally developed to simulate flows around heart valves. Direct comparison of FSI simulations around heart valves against experimental data is challenging, however, due to the difficulty of performing robust and effective simulations, the complications of modeling a specific physical experiment, and the need to acquire experimental data that is directly comparable to simulation data. Such comparators are a necessary precursor for further formal validation studies of FSI simulations involving heart valves. In this work, we performed physical experiments of flow through a pulmonary valve in an in vitro pulse duplicator, and measured the corresponding velocity field using 4D flow MRI (4-dimensional flow magnetic resonance imaging). We constructed a computer model of this pulmonary artery setup, including modeling valve geometry and material properties via a technique called design-based elasticity, and simulated flow through it with the IB method. The simulated flow fields showed excellent qualitative agreement with experiments, excellent agreement on integral metrics, and reasonable relative error in the entire flow domain and on slices of interest. These results illustrate how to construct a computational model of a physical experiment for use as a comparator.

Appendix

We performed a convergence study to evaluate the sensitivity of the phase-averaged velocity fields and flow rates to changes in simulation resolution. The fluid resolutions were set to $\Delta x =$ 0.9, 0.68, 0.45, 0.34 and 0.28 mm. We ran each simulation ran for two cardiac cycles, then performed phase averaging on the second cycle and sampled onto the MRI mesh for comparisons, as described in “Post Processing” section. For $\Delta x =$ 0.9 and 0.68 mm, we used a coarsened structure mesh with a target edge length of $\Delta s = 0.45$ mm, for $\Delta x =$ 0.45 and 0.34 mm we used the target edge length of $\Delta s = 0.225$ mm, as used in the remainder of this work, and for $\Delta x =$ 0.28 mm we use the leaflet mesh targeted to $\Delta s = 0.225$ mm, and a slightly finer mesh for the vessel targeted to $\Delta s = 0.2$ mm.

Exact convergence in IB method simulations is challenging to achieve due to diffuse coupling of the fluid–structure interface and the physical instability of the flow. The radius of the support of the discrete $\delta$-function is $2.5\Delta x$, which depends on the fluid mesh width, meaning that the effective thickness of all objects in the structure depends on the fluid resolution. Thus, IB method coupling reduces the effective radius of the valve and vasculature proportional to $\Delta x$, meaning that the effective radius converges to the true radius in a first order manner as $\Delta x \rightarrow 0$. Resistance to forward flow depends more than linearly on radius (to the fourth power in the case of Poiseuille flow), and thus slow convergence of the flow rate under a specified pressure difference is expected. Alternative schemes could include compensatory measures to adjust the boundary conditions to achieve a more constant flow rate, but this would not reduce numerical error at the fluid–structure interfaces. Additionally, the Reynolds number of the flow is much greater than one, the flow is inertially dominated and flow structures are unstable in time. This effect is partially mitigated by phase averaging the flow.

The phase-averaged, resampled velocity fields during peak systole and flow rates at each resolution are shown in Fig. 10. Despite the limitations discussed above, we observe similar qualitative trends in the flow field at all resolutions. At all resolutions, a jet formed and angled up downstream of the valve orifice, as shown in the sagittal view. The jets showed a triangle-like cross section at $x = 0$ with points aligned with the commissures. At $x = 0.625$ cm, the jet appears like a rounded triangle in the opposing orientation, with its points aligned with the center of the leaflets. At $x = 1.25$ cm, the jet is narrower downstream of the commissures, and wider downstream of the leaflets, again with a triangle-like cross section. The area of the jet increased with resolution, as expected given the IB method thickening of the valve structure. The narrowed jets at the two more coarse resolutions show locally elevated velocities relative to the two more fine resolutions. Figure 11 shows the instantaneous velocity fields at each resolution in the same axial and sagittal views. At $\Delta x =$ 0.9 mm, the sagittal view shows a qualitatively different jet than at finer resolutions, with regions of lower velocity farther from the vessel wall, indicating insufficient resolution. At $\Delta x =$ 0.9 and 0.68 mm, the jet is visibly narrowed compared to higher resolutions. While some features are similar at these two coarse resolutions, we conclude that the narrower jets indicate these simulations are under-resolved. Flows in the three finest resolutions, $\Delta x =$ 0.45, 0.34 and 0.28 mm. appeared qualitatively similar, with slightly more fine structure detail in both the axial and sagittal views present at the edges of the jet. The jets in the axial views all showed a similar triangle-like cross section, slightly narrower downstream of the commissures, as in the phase-averaged fields. In both the phase-averaged and instantaneous fields, the three finest resolutions appear sufficiently similar that the conclusions of this study would be identical with any of these resolutions.

Fig. 10

Comparison of three axial view slices and one sagittal slice of phase-averaged velocity across multiple resolutions

Fig. 11

Comparison of three axial view slices and one sagittal slice of instantaneous velocity across multiple resolutions

The flow rates, depicted in Fig. 12, increase with increasing resolution, also as expected given the increase in effective valve orifice area and radius. The maximum flow rates at $\Delta x =$ 0.45 and 0.34 mm were 248.6 and 281.4 ml/s, respectively, a relative difference of 11.7%. The mean flow rates at $\Delta x =$ 0.45 and 0.34 were 55.5 and 63.2 ml/s, respectively, a relative difference of 12.4%. The maximum flow rates at $\Delta x =$ 0.34 and 0.28 mm were 281.3 and 299.8 ml/s, respectively, a relative difference of 6.1%. The mean flow rates at $\Delta x =$ 0.34 and 0.28 mm were 63.2 and 67.4 ml/s, respectively, a relative difference of 6.2%. Thus the flow rates show signs of converging with increasing resolution, though some minor differences remain.

Fig. 12

Flow rates through the right ventricular inlet with various resolutions

Figure 13 shows the integral metric $I_{1}$ evaluated at $x = 0$ cm with various resolutions, where each resolution uses its own velocity scale (see Eq. (18)). Values at $\Delta x =$ 0.9, 0.68 were elevated relative to finer resolutions, indicating under resolution. Values of $I_{1}$ at $\Delta x =$ 0.45, 0.34 and 0.28 mm are extremely similar, with slight decrease as resolution increases. Values of $I_{1}$ at $x =$ 0.625 and 1.25 cm showed similar trends and are not shown.

Fig. 13

The integral metric $I_{1}$ evaluated at the origin with various resolutions

Given the overall qualitative similarity and relative changes in flow rates, we believed all conclusions in this work would be similar with $\Delta x =$ 0.45, 0.34 or 0.28 mm. We therefore selected $\Delta x =$ 0.45 mm, twice the MRI resolution, for the main portion of this study.

Additionally, we compared averaging the second cycle with phase averaging the second, third and fourth cycles with a resolution of $\Delta x = 0.45$ mm. We found the fields appeared qualitatively similar in both cases. Averaging the second, third and fourth cycles resulted in slightly worse relative error compared to experiment than the second alone. Comparing the two averaging methods at peak systole, relative error on velocity magnitude, the x-component of velocity was 0.07, and 0.09, respectively in $L^{1}$, and 0.10 and 0.11 in $L^{2}$. As these differences were substantially lower than relative differences between the simulation results and experiments, we ended all further simulations after the second cycle and use only the second cycle for phase averaging.