Controlled Comparison of Simulated Hemodynamics across Tricuspid and Bicuspid Aortic Valves

The model aortic valves were constructed using a design-based approach to elasticity, as previously introduced in [32]. We specified that the tension in the leaflets must support a pressure, then derived an associated system of partial differential equations. The solution of this system specified the loaded model geometry and tensions needed to support the specified pressure load. Using this information, a reference configuration and constitutive law were derived. This configuration automatically included local fiber structure and heterogeneous material properties that would be challenging to measure experimentally. By tuning free parameters in these differential equations, this construction allowed us to design a model that is consistent with known anatomy and material properties. This method was previously applied to create realistic and effective mitral valve models for FSI simulations [30, 29].

We created four model aortic valves. The first is trileaflet and serves as a control, and is identical to the model derived in [32]. Three are bicuspid, representing LC/RC, RC/NC and NC/LC fusion. Aside from fusion at the free edge, all four models have identical material properties and leaflet reference geometry.

We represented the valve leaflets as an unknown parametric surface in $\mathbb{R}^{3}$,

The parametrization was assumed to conform to the fiber and cross fiber directions, meaning that curves of constant $v$ ran circumferentially in the leaflets and curves of constant $u$ ran radially. Let subscripts denote partial derivatives. The unit tangents in the circumferential and radial directions in the leaflets were defined as

The leaflet exerted tensions $S$ and $T$ in the circumferential and radial directions, respectively. Shear tensions were assumed to be identically zero. Let $p$ denote the pressure supported by the leaflets. Consider the equilibrium of pressure and tension on an arbitrary patch of leaflet corresponding to $[u_{1},u_{2}]\times[v_{1},v_{2}]$. This dictated that tension integrated over on the boundary of the leaflet and pressure integrated with respect to area sum to zero, or

We then applied the fundamental theorem of calculus, differentiating and integrating to convert tensions integrated on the boundary of the patch to derivatives of tension integrated on the entire patch. The order of integration was swapped formally and the integrals were combined to obtain

Since the patch of leaflet is arbitrary, the integrals can be dropped. This process gave the following system of partial differential equations for equilibrium of the leaflets:

The equations (5) have three components and five unknowns, and so are not closed. To close them we temporarily specified that

The values $\alpha,\beta,a,b$ are tunable free parameters that are not required to be constants. The parameters $\alpha,\beta$ specify the maximum tension that each fiber can achieve. The parameters $a,b$ can be tuned to influence spacing of fibers and achieve realistic gross morphology in the closed configuration of the valve. Equations (6) are not meant to represent a physical constitutive law, rather they provided a way to arrive at a closed configuration for the valve and the heterogeneous tensions it supported in its loaded state without specifying either a reference configuration or physical constitutive law. Based on preliminary experiments in prior work [30], we expect the leaflet geometry to be fairly insensitive to the specific functional form of equations (6), and that other smooth functions varying from zero to the maximum allowed tension would work similarly. The emergent leaflet geometry, however, depends strongly on the selected coefficients.

Incorporating the tensions (6), the equilibrium equations (5) were discretized with centered finite differences. Let $\mathbf{X}^{j,k}$ denote an arbitrary node in the discretized leaflets. The nonlinear system of equations that corresponded to this node was given

The annulus and commissure positions were prescribed as Dirichlet boundary conditions. The free edges were treated with homogeneous Neumann (zero tension) boundary conditions. The nonlinear system of equations was solved with Newton’s method with line search.

Next, we used the solution to equations (5), which represents the closed, loaded configuration of the valve, to specify a reference configuration and constitutive law. We prescribed uniform, experimentally measured circumferential strains of $E_{u}$ = 0.15 and radial strains of $E_{v}$ = 0.54 radially [55]. (In-membrane shear strain was not reported in [55]. We do not include it here nor expect it to contribute substantially were it to be included.) The equation $E=L/R-1$ was solved for the reference length $R$, where $E$ denotes engineering strain and $L$ is the current length. (Note that stretch ratio, $L/R$, could be used and is equivalent; use of engineering strains is not meant to denote that deformations are small.) The tension/strain relation for each link was taken to be exponential through the origin with exponential rates of $57.46$ circumferentially and $22.40$ radially [41]. We solved for the stiffness coefficient $\kappa$ in the equation

to scale stiffness for each link to match tension $t$ at the relevant strain $E$.

This process produces material properties that are comparable to those measured in in vitro testing. The tension in the solution to equations (5) is heterogeneous, so each link in the model has its own stiffness coefficient. This process thus created heterogeneous material properties. We estimated the tangent modulus at the loaded circumferential strains of $E_{u}$ and radial strains of $E_{v}$, and compared to an experimental study on human data that tested a small patch on the central or belly region of the leaflet [45]. The circumferential tangent modulus has minimum, maximum and mean values of $6.6\cdot 10^{7}$, $1.7\cdot 10^{9}$ and $1.4\cdot 10^{8}$ dynes/cm², respectively, compared to an experimental estimate of $9.9\pm 1.8\cdot 10^{7}$ dynes/cm² [45]. The radial tangent modulus has minimum, maximum and mean values of $2.4\cdot 10^{5}$, $3.2\cdot 10^{7}$ and $5.7\cdot 10^{6}$ dynes/cm², respectively, compared to an experimental estimate of $2.3\pm 0.4\cdot 10^{7}$ dynes/cm²[45]. At the belly of the leaflet, we estimated the pointwise tangent modulus to be $9.4\cdot 10^{7}$ dynes/cm² circumferentially and $8.5\cdot 10^{6}$ dynes/cm² radially. Thus, at the belly of the leaflet, our model has circumferential tangent modulus that agrees with this experiment to within standard error. Further, along the midline of the leaflet, the radial tangent modulus varies from a maximum of $3.0\cdot 10^{7}$ dynes/cm² at the attachment of the leaflet and annulus to a minimum of $4.0\cdot 10^{5}$ dynes/cm² at the free edge. The radial stiffness of our model on the centerline varies from just above their experimentally reported range at the annulus, through the range for a portion of leaflet, then is lower at the free edge, where the leaflet is expected to be more compliant radially. Thus, the tangent moduli in our model agree with measured data in certain specific, relevant locations. However, the tangent moduli are heterogenous and thus not equal everywhere. Studies have found substantial variety in the material properties of leaflet tissue, especially near the commissures [9]. We observe local increases in the circumferential tangent modulus near the commissures (see [30]), which explains slightly elevated mean values compared to those of [45].

The resulting fiber orientations, which run circumferentially, are comparable to experimental data, especially in the belly region of the leaflet, as discussed further in [30]. The model exerts tension only in the circumferential and radial directions, and we do not include in-plane coupling terms. We are aware of little data that addresses the material response to in-plane shear, but one study showed that the changes to the closing kinematics of the leaflet were minimal with or without such terms on a similar material model [25]. In future work, these radial and circumferential tensions and strains could be fit to a variety of constitutive laws including those with in-plane coupling or other general hyperelastic laws. One minor limitation of these models is that they lack bending rigidity (beyond what emerges from the thickening process described below), and hence may not capture fluttering behavior in the leaflets accurately. We do not alter material properties in the bicuspid leaflets from those of the tricuspid case, as to avoid confounding factors associated with pathological remodeling of bicuspid valves and further isolate the effects of changes to leaflet geometry alone. Thus, the material properties prescribed and derived for our model agree with experimental data in terms of tangent modulus, shape of the tension/strain curves and fiber orientation.

To obtain initial conditions for FSI, we solved for open configurations of the tri- and bicuspid valve as follows. Using this newly generated constitutive law for tensions, we again solved the generic equations of equilibrium (5) with zero pressure. The positions of the free edges were prescribed as a Dirichlet boundary condition to ensure an open configuration. Evenly spaced points were computed on a chord between the two commissure points on either side of this leaflet. Each free edge point was placed on a line between the annulus and the corresponding point on the chord. The distance from the annulus was set to be the total reference length in the radial direction associated with this radially-oriented curve multiplied by $E_{v}$, the radial strain in the predicted loaded configuration. This distance reduced overall residual tension in the leaflets compared to putting the points directly on the chord or at the rest length from the annulus. To construct a bicuspid valve, two leaflets were prescribed to coincide exactly on half of their free edges. For the fused half of the leaflets, the point was determined to lie at the reference length circumferentially from the commissure, and at the reference length multiplied by $E_{v}$ from each corresponding point on the annulus. The free half of the fused leaflets was set with a spline interpolant from the center point to the relevant commissure. The third leaflet, which has no fusion, is set identically to the tricuspid case. This process created a raphe in the center of the fused cusps. The reference configuration and constitutive law of the bicuspid valve are identical to those of the tricuspid valve; the models differ only in their initial configuration and coinciding points on the free edge. For all models, this process produced an initial configuration that is not the reference configuration; thus the initial condition included pre-strain and residual tension. Further, a configuration that is tension-free may not exist. Experiments have shown that the aortic valve has pre-strain and residual tension throughout the cardiac cycle, including when the valve is open, both in vitro [55] and in vivo [1]. The resultant configurations are shown in Figure 1. Until this point we worked with a zero-thickness membrane and membrane tensions; next we thickened the model to a physiological thickness of $0.044$ cm [47]. Two additional layers were placed normal to the membrane, and membrane stiffnesses were divided by three to obtain the relevant stiffness for each layer. Linear springs attached the layers together in the normal directions; these prevent the leaflets from separating from each other. Their stiffnesses were tuned empirically such that the layers undergo minimal movement from each other without causing further time step restrictions. This also served to mitigate the “grid-aligned artifact” that occurs with pressure discontinuities across thin, tensioned membranes in the IB method [30, 29].

Image segmentations to construct a patient-specific model aortic geometry were constructed from CT scan data using SimVascular [35]. The patient was a 50 year old male with no known complications in the aorta or aortic valve. A pathline for the approximate centerline of the vessel was created. Two-dimensional contours were segmented manually on slice views normal to the pathline. These contours were then lofted into a surface, and a triangular mesh of this surface was exported. The aorta model was held in an approximately constant position with target points and linear springs on each edge of the triangular mesh (Section 2.3). Keeping the aorta approximately rigid serves to further isolate the effects of changes to valve morphology alone, because in the bicuspid cases it is unclear whether to use healthy material properties to remove effects of material property changes due to disease, or material properties from adult bicuspid patients, which potentially be influenced by pathological remodeling. Because the fluid domain was a Cartesian box, flow extenders of length 1 cm were added to the inlets and outlets. This ensured that the normal to the vessel was aligned with the normal to the fluid domain at the inlet and outlet. An interpolant connected the aorta from immediately below the aortic annulus to the flow extender at the edge of the fluid domain. As with the valve, two additional layers were added by extruding the model in the normal direction at one and two times the structure mesh width, or $0.025$ and $0.05$ cm.

We positioned the valve manually in the aorta model as follows. We imported the aorta geometry into Paraview [2], then applied the slice operator and interactively manipulated its position to determine the origin and normal of the minimum plane of the annulus. A rigid body rotation and translation were applied to align the axis and origin of the valve with the axis and normal of the minimum annulus plane. We applied an additional rotation, computed by trial and error, around the axis of the valve to ensure that the commissures align with the aortic sinus. The model has exactly a 120 degree angle between each of the commissures, and this aligned well with the patient’s aortic sinus after segmentation. Expert clinicians confirmed the anatomical accuracy of this placement by visual inspection of the models.

Simulations were run using the immersed boundary (IB) method for FSI [44]. In the IB method, fluid quantities are represented with respect to the lab or Eulerian frame. Structure quantities, such as the aorta and valve, are represented with respect to a material or Lagrangian frame. The fluid occupies the entire domain of interest. The structure occupies some subset of this domain and is treated as neutrally buoyant, meaning that it has no additional mass beyond the mass of the local fluid. The two frames are coupled via convolutions with the Dirac-$\delta$ function, as described below. The structure influences the fluid through a body force computed via these convolutions, which is one distinctive feature of the IB method.

The fluid velocity and pressure are denoted $\mathbf{u},p$ respectively. The Eulerian frame force $\mathbf{f}$ represents force exerted by the structure onto the fluid. These functions take arguments of fixed spatial position $\mathbf{x}$ and time $t$. The constants $\rho$ and $\mu$ denote fluid density and dynamic viscosity. Let $\mathbf{X}(\mathbf{s},t)$ denote the current configuration of the structure, where $\mathbf{s}$ labels a material point. The Lagrangian frame force $\mathbf{F}$ denotes the force exerted by the structure in the Lagrangian frame. (We previously used $u,v$ to label material points, but switch to $\mathbf{s}$ to avoid confusion with fluid velocity here.) The Dirac delta function is denoted as $\delta$.

The governing equations of the IB method are

Equations (9) and (10) are the Navier Stokes equations for conservation of momentum and volume of an incompressible Newtonian fluid, plus the Eulerian frame force $\mathbf{f}$. Equation (11) is a mapping from the configuration of the structure to the force it exerts on the fluid. The omitted argument indicates that the mapping $\mathcal{F}$ takes the entire configuration of the structure at time $t$ and produces the Lagrangian frame force $\mathbf{F}$. Equations (12) and (13) are interaction equations that couple the fluid and structure. The interpolation equation (12) specifies that the structure moves at the local fluid velocity. The force spreading equation (13) determines the Eulerian frame force $\mathbf{f}$ from the Lagrangian frame force $\mathbf{F}$.

Simulations were conducted with the software library IBAMR (Immersed Boundary Adaptive Mesh Refinement) using a staggered Cartesian grid for the fluid [22, 21]. Simulations were run for three cardiac cycles at 75 beats per minute, or 0.8 s per beat. The fluid density and viscosity were set to $\rho=1.0$ g/cm³ and 0.04 Poise, respectively. The velocity was initialized to zero and the simulation began at early diastole. The second and third cardiac cycles appeared qualitatively similar in all simulations and sufficient to wash out initialization effects; no further cycles were run due to to computational expense and similarity of further cycles in previous validation tests. Simulations were run on Stanford University’s Sherlock cluster on 48 Intel Xeon Gold 5118 cores across two nodes with a 2.30 GHz clock speed.

The fluid mesh width was set to $\Delta x=$ 0.05 cm. The fluid domain was taken to be $7.2\times 4.8\times 11.2$ cm, corresponding to $144\times 96\times 224$ points. The structure mesh width was targeted to half that of the fluid mesh width, or $\Delta s=\Delta x/2=0.025$ cm. The lengths of individual links were determined as described in Section 2.1 and evolved according to the FSI dynamics. A scaffold was placed around the valve to ensure that there are no holes between the model valve and aorta. The aorta, scaffold and the annular edge of the leaflets were held in place by target points, a penalty method that approximately enforces a fixed position. For a point $\mathbf{X}$ and its desired location $\mathbf{X}_{target}$, this is a force of the form $\mathbf{F}=-k(\mathbf{X}-\mathbf{X}_{target})$, representing a linear spring of zero rest length connecting the current and desired positions with $k=5.8\cdot 10^{4}$ dynes/cm. A 5-point discrete $\delta$ function was used [4]. At the outlet, open-boundary stabilization was applied with penalty coefficient $\eta=\rho/(4\Delta t)$ dynes s/cm⁴ [8]. Open boundary stabilization was necessary to prevent spurious, non-physical flows from emerging from the intersection of the aortic domain with the boundary of the fluid domain. At the inlet, a flow-averaging force $\mathbf{f}=-\eta(u-\bar{u},v,w)$ was applied, where the $x$ direction is normal to the inlet and $\bar{u}$ is the mean inflow and $\eta=\rho/(4\Delta t)$ dynes s/cm⁴. This approximately enforced a uniform or plug profile to the inflow. Empirically, this was necessary to avoid numerical instabilities in the form of rapidly increasing spurious flows at the boundary of the aortic inlet. This was the only stabilization applied at the inlet, but in the even of a local flow reversal, this stabilization applies a similar force to that of open boundary stabilization. Outside of the aorta, a force $\mathbf{f}=-\eta\mathbf{u}$ was applied to damp flow to zero with $\eta=\rho/(4\Delta t)$ dynes /cm⁴. In the first 0.1 s, this force was applied on the entire domain to damp out initial transients associated with residual stress in the initial conditions.

Pressure boundary conditions created a pressure gradient across the valve and aorta to drive the simulations. (Here the term “pressure gradient” denotes the difference in two local values of pressure in mmHg, as is standard in clinical literature.) To prescribe the pressure boundary conditions, we prescribed zero tangential slip and applied a normal traction with value equal to the negative of the desired pressure. Applying the incompressibility constraint, zero tangential slip implies that this normal traction is then equal to the pressure locally. We targeted aortic systolic, diastolic and mean pressures of 120, 80 and 96 mmHg respectively and a flow rate of 5.6 L/min. The aortic pressure was initialized to be equal to the initial ventricular pressure of 23 mmHg at $t=0$ s, meaning that the initial pressure gradient was zero. The aortic pressure then ramped linearly to 94 mmHg at $t=0.1$. For the remainder of the simulation, the pressure at the aorta side outlet was determined by an RCR lumped parameter network [33]. The total resistance was computed as the ratio of the targeted mean pressure over the targeted mean flow. The ratio of the distal to proximal resistors was set to $0.064$ [37]. Capacitance was set using an exponential exact solution to the governing differential equations in the case of assumed zero flow in diastole. This procedure gave values of proximal resistance of $R_{p}=83.67$ s dynes cm^-5, distal resistance of $R_{d}=1287.65$ s dynes cm^-5 and capacitance of $C=0.0017\text{ dynes}^{-1}\text{ cm}^{5}$. The pressure on the left ventricular side was prescribed based on experimental measurements [56]. During peak systole, a pressure loss was observed across the location of the flow averaging force, which removes energy from the flow locally. Thus, we added 15 mmHg to the systolic pressure upstream of the flow extender to obtain the desired physiological ventricular systolic pressure of 120-130 mmHg immediately proximal to the valve. In results, we report ventricular pressures taken from the pointwise fluid pressure proximal to the valve and downstream of the flow extender. Outside of the aorta, no slip boundary conditions were prescribed on faces that include the inlet and outlet, keeping these structures stationary through the simulation. Zero pressure boundary conditions were prescribed on all other faces.

To validate these methods, we conducted a physical experiment with a 3D-printed vascular geometry and prosthetic aortic valve and measured its flow field via 4D MRI, then, with the same methods used in this work, simulated the physical experiment [31]. We achieved excellent qualitative agreement on the hemodynamics and valve kinematics between simulation and experiment. We computed the $L^{1}$ relative error of streamwise velocity and velocity magnitude on the entire flow domain and relevant two-dimensional slices, and achieved reasonable qualitative agreement of order 30% during forward flow on all quantities. Other studies on flows through heart valves have performed qualitative comparisons or quantitative validation of leaflet kinematics, but we not are aware of previous studies have directly compared the three-dimensional velocity fields and achieved agreement in relative error. We also compared varying levels of resolution for the bicuspid aortic simulations shown in the current paper, which produced qualitatively and quantitatively similar results (see Appendix). Thus, we consider the hemodynamics that result from our method to be well-validated, and these comparisons ensure that downstream hemodynamics shown in this work are realistic.

To quantify the nature of the flows, we computed a number of integral metrics on two-dimensional cross sections of the flow. These metrics were previously used to study flow in the lungs [3] and through the pulmonary valve [50]. The two-dimensional cross sections are denoted $S_{1},S_{2}\dots S_{5}$, numbered streamwise, and located at the aortic annulus, the sinotubular junction and three additional cross sections in the ascending aorta.

The metric $I_{1}$ represents the nondimensional streamwise momentum flux and is given as

where $A$ is the area of the cross section, $U_{T}$ is a constant velocity scale, and $\mathbf{n}$ is the unit normal to the cross section. This is the $L^{2}$ norm of the normal component of velocity, made non-dimensional by inclusion of the velocity scale and area. In the case of a constant or plug velocity profile with velocity equal to $U_{T}$, the velocity scale, $I_{1}=1$.

The metric $I_{2}$ represents secondary or tangential flow strength and is given as

Similarly, this is the $L^{2}$ norm of the tangential component of velocity on the cross section, nondimensionalized by the inclusion of the velocity scale and area.

For both $I_{1}$ and $I_{2}$, the velocity scale was set to the spatial average of the maximum instantaneous flow rate,

where we explicitly notate the space and time variables for clarity. This value was computed for each case individually.

Finally, we compute $I_{R}$, the fraction of reverse flow on the slice. The indicator or characteristic function of the set on which flow moves in the reverse direction is defined as

and the metric is defined as

This value represents the portion of cross sectional area in which flow moves towards the left ventricle, rather than up the aorta.

To evaluate the flow fields and integrals on the selected cross sections, the flow field was first restricted to the interior to the aorta model using PyVista [52]. Using Paraview filters [2], data were projected to the normal and tangential components using Gram Scimidt orthogonalization, two dimensional slices were extracted from these data, then integrated.

Simulations were run with a tricuspid valve and three bicuspid bicuspid valves: LC/RC fusion, RC/NC fusion and NC/LC fusion. In all four cases, the simulated valves functioned as expected over multiple cardiac cycles, opening to allow forward flow and closing under back pressure. The flow fields and forward jets showed substantial differences between all cases. Visualizations of velocity are shown in Figure 2, depicting the vertical component of velocity on a slice through the annulus, and the normal component of velocity on fives slices approximately normal to the local axis. These slices are located the aortic annulus, the sinotubular junction, and three locations in the ascending aorta (numbered $1\dots 5$ streamwise). Movies of these flows are included in the supplemental information. Vector plots of the tangential component of velocity are shown in Figure 3, which illustrate the locations and orientation of secondary flows.

In the tricuspid case, the flow appeared plug-like immediately after opening. As the flow accelerated during peak systole, a jet developed through the center of the valve orifice, then separated slightly from the inner wall at mid systole. At peak and mid systole, there were regions of slight recirculation back towards the valve orifice. We refer to the outer side of the aorta with greater radius of curvature, also called the convexity, as the greater curvature and the inner side with lesser radius of curvature, also called the concavity, as the lesser curvature. As the valve began to close, the jet velocity decreased and the flow reversed along the lesser curvature first, where the flow was already decelerating due to the jet angle towards the greater curvature.

In each of the bicuspid cases, the valve orifice was far narrower than in the tricuspid case, and a more complex flow pattern with greater recirculation occurred. Due to the orifice orientation, the jet developed in close approximation to the vessel wall, in contrast to the tricuspid case where the jet remained approximately centered in the vessel lumen. This difference in location and increased jet intensity produced substantially different flow features in all four cases.

In the bicuspid case with LC/RC fusion, the valve orifice was oriented along the greater curvature. At peak systole, a concentrated jet developed immediately distal to the valve orifice and remained adherent to the greater curvature. As the jet progressed downstream, it spread along the greater curvature normal to the flow direction while maintaining intensity. There was a region of reverse flow distal to the fused LC and RC leaflets. The strongest reverse flows were present at peak systole, generating recirculation in the vertical component of velocity (Figure 2). Downstream of the aortic sinus, lower velocity fluid moved away from the center of the jet along the vessel wall towards the lesser curvature. This generated a counter-rotating vortex pair at peak systole in the middle of the ascending aorta, which began to shift and lose strength as the jet spread along the greater curvature at mid systole (Figure 3). This vortex pair structure (slices 3 and 4) persisted with an additional secondary flow rotating counterclockwise that extended along almost the entire circumference of aorta. Further downstream (slice 5), the larger counterclockwise rotation had formed into a vortex and the counter-rotating vortex pair did not form.

In the bicuspid case with RC/NC fusion, the valve orifice was oriented along the lesser curvature, and correspondingly at the sinotubular junction, the jet was centered along the lesser curvature. As the flow propagated, the jet traveled along the wall from the lesser curvature to the greater curvature and began to spread. Further downstream, (slice 5) there is reverse flow at the lesser curvature of the aorta, directly distal of the valve orifice. The forward flow extended into the center of the vessel, especially at mid systole, separating from the wall more than in the LC/RC fusion case. Strong secondary flows accompanied the movement of the jet. An asymmetric counter-rotating vortex pair was generated that was eventually subsumed by a counterclockwise rotation around the vessel wall.

With NC/LC fusion, the valve orifice was oriented along the greater curvature. The streamwise flow behaved similarly to the LC/RC fusion case. The jet was concentrated immediately downstream of the valve orifice and then spread as flow progressed through the aorta, with reverse flow distal to the fused leaflets. The secondary flows have rotation along the vessel wall that traveled clockwise, particularly in the ascending aorta, and no distinct counter-rotating vortex pair appeared. Therefore, the secondary flows are substantially different from the two other bicuspid cases.

Thus, in all cases, the jet ultimately traveled towards the greater curvature of the aorta. The LC/RC and RC/NC fusion cases generated a counterclockwise rotation, whereas the NC/LC case generated a clockwise rotation. In addition, RC/NC fusion produced more substantial rotating flow, as the jet traveled from the lesser to greater curvature.

Pressure and flow waveforms for the tricuspid and LC/RC fusion cases are shown in Figure 4. The second and third cardiac cycles are shown; the first cardiac cycle contains initialization effects and is omitted. The aortic pressure decreased gradually through diastole, while the ventricular pressure remained low. Next, the ventricular pressure increased as systole began. A transient forward pressure gradient across the valve occurred, with peak values as high as 16.3 mmHg. Then, the pressure gradient decreased, leaving minimal pressure gradient across the valve. At times, the pressures even crossed, while forward flow continued due to inertia. Thus, the model valve offered minimal resistance to flow. When the ventricular pressure began to drop, a prominent dicrotic notch occurred in the aortic pressure, followed by an oscillation associated with valve closure. This oscillation is similar to that seen in experimental pressure traces [51, 36] and is caused by a sudden pressure loading of the elastic valve. The simulation oscillation has possibly higher amplitude and longer duration due to the rigidity in the model aorta than those of experiments, but noise and sampling rates in experimental traces makes this difficult to assess. An oscillation also occurs in the ventricular pressure immediately following valve closure, again due to sudden pressure loading on the valve (recall that the ventricular pressure reported is the pointwise fluid pressure immediately proximal to the valve). The cycle then repeated, with pressure gradually decreasing through diastole. On the third cycle, trace regurgitation occurred and a corresponding lack of smoothness is seen in the pressure curve.

The flow waveform showed a flow rate of approximately zero during diastole, as expected. At the beginning of systole, the flow rate increased rapidly under the transient forward pressure gradient. Next, the ventricular pressure began to fall below the aortic pressure, the flow rate decreased and subsequently became negative. An oscillation was seen in the flow rate that then decayed and was followed by a flow rate of approximately zero for the remainder of diastole. A similar oscillation occurs in experimental measurements of flow rate, though the simulation oscillation is of possibly longer duration due to rigidity in the model aorta, as in the oscillation in pressure [51, 36]. The cycle then repeated. The cumulative flow had an average value of 75.24 ml/cycle, or 5.64 L/min, corresponding to the nominal cardiac output of 5.6 L/min. The peak flow rate was $Q_{max}$ = 522.2 ml/s. We estimate the peak Reynolds number as $\rho(Q_{max}/A)2r/\mu\approx 6600$, where $r=1.25$ cm is the valve radius and $A=\pi r^{2}$ is the annulus area.

With LC/RC fusion, a forward pressure gradient across the domain was sustained through systole, with peak pressure gradient of 30.6 mmHg and sustained pressure gradients over 18 mmHg. These pressure gradients represent clinically significant, stage B, progressive, mild-to-moderate aortic stenosis [42]. This model does not include calcification or other material property changes. The stenosis, therefore, was exclusively caused by fusion of the free edges. The occurrence of stenosis is in stark contrast to the tricuspid case, in which a near zero pressure gradient occurs throughout forward flow, as expected for normal aortic valve morphology.

The bicuspid flow cases had lower forward flow during systole, with a cumulative flow rate of 56.90 ml/cycle, or 4.27 L/min. Thus, without any changes to boundary conditions or modeling any compensatory measures in the circulatory system, the stenosis associated with the bicuspid valves resulted in a loss of about 1.37 L/min (or 24%) of cardiac output. The valve sealed reliably for all three cycles. Pressure and flow rates for the RC/NC and NC/LC fusion cases are similar to those of the LC/RC fusion case.

Line plots of velocity for all cases are shown in Figure 5, depicting the velocity normal to the slices along three lines through the center of each slice, enabling direct comparison. In the tricuspid case, at the annulus the flow was plug-like. Immediately downstream of the valve, the center of the jet had nearly constant velocity. Some recirculation occurred along the edge of the jet, particularly along the greater curvature at the sinotubular junction (slice 2). Further downstream, trace local reverse flow appeared, especially on the lesser curvature. The normal velocity later became somewhat uniform in the ascending aorta, as the wide jet spread throughout the vessel.

For all three bicuspid valves, flow proximal to the valve was less symmetric. Distal to the valve, this asymmetry became pronounced. There was a localized plug-like jet immediately distal to the valve, corresponding to the unrestricted leaflet in each case. The line plots for downstream slices capture how the jet dispersed and spread along the outer wall in each case, creating secondary and reverse flows.

We further analyzed the flow with three non-dimensional integral metrics: $I_{1}$, representing the excess streamwise momentum, equation (14), $I_{2}$, representing secondary flow strength, equation (15), and $I_{R}$, the fraction of reverse flow, equation (18).

Figure 6 shows the integral metrics in the final 0.4 s of simulation, including all of systole in the final cycle. For the tricuspid case, the slices at the annulus and immediately downstream of the valve had $I_{1}$ values greater than 1, indicating higher momentum than uniform flow. As the flow moved downstream, the streamwise momentum decreased, with the slices in the ascending aorta having $I_{1}$ values of approximately 1, indicating a return to uniform flow.

For the bicuspid cases, the streamwise momentum at slice 1, the annulus slice, is similar to the tricuspid case throughout the cardiac cycle. At slices 2 and 3, the streamwise momentum drastically increased, with values of $I_{1}$ reaching nearly 2.5, thereafter decreasing at slice 4 and more still at slice 5. The exception was in the LC/RC fusion case, where $I_{1}$ did not substantially decline until slice 5, as opposed to slice 4 in the other two bicuspid cases. Thus, higher momentum flow persisted further downstream with LC/RC fusion.

The metric $I_{2}$ represents the non-dimensional secondary flow strength and quantifies the vortices and rotational motion seen in Figure 3. For the tricuspid case, aortic flow exhibited low $I_{2}$ throughout the cardiac cycle. For the bicuspid cases at slice 1, though $I_{2}$ remained low, its values are still greater than those of the tricuspid case. Immediately downstream of the valve and through the ascending aorta, $I_{2}$ in the bicuspid cases is substantially higher than in the tricuspid case. For all three bicuspid cases, $I_{2}$ was higher at slice 3 than slice 2, as vortex structures developed downstream. With LC/RC fusion and NC/LC fusion, $I_{2}$ increased again at slice 4, indicating strong rotating flow along the vessel wall. In the LC/RC case, the counter-rotating vortex pairs also contributed to $I_{2}$. The value of $I_{2}$ decreased at slice 5 as the jet and secondary flow structures broke up. In RC/NC fusion, $I_{2}$ is substantially higher at slices 4 and 5, as the jet moved from the lesser to greater curvature. Thus, differences in leaflet morphology cause higher secondary flow strength further downstream of the valve.

At a given aortic cross section, $I_{R}$ defines the fraction of area in which flow moves in the reverse direction. During systole, the tricuspid case generally had low amounts of reverse flow. Notably, there is the most reverse flow at the sinotubular junction (slice 2). As the flow decelerated into early diastole, the value of $I_{R}$ increased to nearly one.

For all bicuspid cases, $I_{R}>0.5$ was seen in slice 2 at early systole. Since the jet was localized to the area above the freely moving cusp in each case, there was more area for reverse flow to develop and flow back towards the fused cusps. This reverse flow generated large areas of recirculation in the aorta (Figure 2). The value of $I_{R}$ remained above 0.3 for almost the entirely of systole in all three bicuspid cases. The time at which each slice first experienced a high percentage of reverse flow increased monotonically downstream, illustrating that reverse flow was not generated until the jet reached that slice.