Determination of a Strain Energy Density Function for the Tricuspid Valve Leaflets Using Constant Invariant-Based Mechanical Characterizations

The tricuspid valve (TV) is located between the right atrium and the right ventricle of the heart. Anatomically, the TV consists of three leaflets that are attached to the surrounding heart tissue by the ring-like TV annulus and anchored to the papillary muscles on the right ventricular wall by the chordae tendineae. These subvalvular components work harmoniously to regulate the blood flow of the heart for the millions of cardiac cycles in a lifespan—the leaflets will relax to allow blood flow from the right atrium into the right ventricle during diastole, while closing to seal the right atrial-ventricular orifice during systole. If the geometry or function of these subvavular components are altered, the intricate valvular function will likely become compromised and the patient may develop tricuspid regurgitation [1–5]. Historically, tricuspid regurgitation was often overlooked due to the common sentiment among cardiologists that it would naturally regress after treatment of left-sided lesions [6]. However, the hallmark clinical investigation by Dreyfus et al. [7] revealed that this may not be accurate, and tricuspid regurgitation should receive more deliberate clinical attention.

The TV has received notably more attention within the clinical and scientific communities after the change of clinical perspectives by Dreyfus et al. [7]. One substantial area that is highly relevant to the present study is the in vitro mechanical characterizations of the TV leaflets. Since 2016, carefully designed biaxial testing experiments have been used to characterize the mechanical behaviors of porcine [8,9], ovine [10], and human [11] TV leaflets. Interested readers are referred to recent extensive reviews [12,13]. Some of these investigations had specific objectives, such as characterizing how short-term freezer storage alters the TV leaflet mechanical properties [14,15], whereas others aimed to establish databases of biaxial tissue mechanics for the TV leaflets that can inform constitutive models for in-silico simulations of the TV function [8–11]. Khoiy et al. [16] provided a constitutive modeling-focused investigation of the TV leaflets. Although their study answered several important questions, they did not determine which strain energy density function (SEDF) is most suitable for the TV leaflets. Thus, there is a pressing need to determine a suitable SEDF for the TV leaflets that can be used in computational biomechanical analyses of the TV function.

Three approaches have been conventionally employed to determine a satisfactory SEDF form for a soft biological tissue. In the first structurally guided approach, considerations for the underlying microstructures of the tissue are utilized to theoretically construct the SEDF form. Holzapfel et al. [17] used this method to establish a well-received SEDF for the arterial wall tissue [18,19]. The second structurally informed approach is quite similar and directly employs the tissue's microstructural information, such as the collagen fiber orientation density function, to construct more complex constitutive models. This has been employed for skin [20], the mitral valve [21], and the aortic valve [22]. In contrast, the third phenomenological approach relies on the empirical assessments of the SEDF's applicability to the soft biological tissue. Experimental data from traditional planar biaxial tension experiments are fit with a series of candidate SEDF forms, and the model, which provides the best fit, is taken as the most suitable. Since a plethora of SEDF forms exist, there is likely at least one acceptable representation of the tissue's mechanical behaviors. However, it can be difficult to delineate which SEDF form is the best for the tissue, especially when several models yield the fits to the experimental data with similar goodness of fit. Therefore, an alternative phenomenological approach uses the experimental data from a series of carefully designed constant invariant-based mechanical characterizations rather than the traditional planar biaxial tension experiments [23]. These unique protocols enable the contributions of the invariants to be decoupled so that a more representative SEDF form can be realized without needing extensive microstructural information like the above-mentioned structurally guided and structurally informed approaches. Thus, the constant invariant-based approach is excellent for emerging areas in soft tissue biomechanics, such as in the infancy of the left ventricle [23] or the mitral valve biomechanics [24], where the extensive microstructural data or information is often times lacking.

The objective of this study is to determine a suitable phenomenological SEDF form for the TV leaflets using the data from constant invariant-based mechanical characterizations. First, specimens for each TV leaflet (n = 6) are characterized using a series of constant I₁ and constant I₄ protocols, and the experimental data is transformed into the partial derivatives of the SEDF. Three candidate SEDF forms are fit to the experimental data using an in-house differential evolution optimization algorithm to estimate the unknown model parameters. The results from these fits are analyzed to determine the most suitable SEDF form based on: (i) model fitting capability, (ii) model parameter correlations, and (iii) the well-established D-optimality and condition optimality criteria. The advantages and disadvantages of coupling the invariant terms within the SEDF are also discussed to provide our recommendation of the most representative SEDF form for the TV leaflets.

Herein, $I1=tr(C), I2=12[tr(C)2−tr(C2)]$⁠, and $I3=det(C)$ are the first three invariants of $C=FTF, I4=N·(CN)$ and $I5=N·(C2N)$ are the pseudo-invariants of C, and F is the deformation gradient. Note that the incompressibility condition, i.e., $I3=det(C)=1$⁠, reduces Eq. (1) to a function of two invariants I₁, I₂ and two pseudo-invariants I₄, I₅. Ideally, one could determine an appropriate SEDF of this form for the TV leaflets; however, as shown in Sec. 2.1.3, it is infeasible to employ planar biaxial testing to fully determine a SEDF that contains invariants more than I₁ and I₄. Thus, we further assume that the SEDF is only a function of I₁ and I₄, resulting in $W=W(I1,I4)$⁠. This agrees with the prior assumptions for the constitutive modeling of the heart valve leaflets [24,26] and is necessary to establish the constant invariant-based mechanical characterization experiment, which will be described next.

Herein, I is the identity tensor, N is the fiber orientation at the undeformed (reference) configuration, $B=FFT$ is the left Cauchy–Green tensor, $W,1=∂W/∂I1$ and $W,4=∂W/∂I4$ are the partial derivatives of W with respect to I₁ and I₄, respectively, and $p=2W,1(λ1λ2)−2$ is the penalty parameter to enforce the incompressiblity condition than can be analytically derived from $det(F)=1$⁠, considering plane-stress condition.

In the above relationships, $ξ1=λ12−(λ1λ2)−2$ and $ξ2=λ22−(λ1λ2)−2$ are the combinations of the tissue stretches. These relationships enable us to directly use the experimental data to compute the partial derivatives of W without needing a predefined SEDF form.

Six porcine hearts (n = 6) were obtained from a regional tissue vendor (Animal Technologies, Inc., Tyler, TX) overnight on dry ice. On the day of experiments, the hearts were thawed and dissected to acquire the three TV leaflets—the anterior leaflet (TVAL), the posterior leaflet (TVPL), and the septal leaflet (TVSL)—which were then frozen for brief storage prior to the mechanical testing within 24 h [14].

An hour prior to the biaxial mechanical testing, the TV leaflets were thawed in an in-house room-temperature phosphate-buffered saline solution. The entire leaflet was promptly mounted to a commercial biaxial testing system—Bio-Tester (CellScale, Canada)—with an effective specimen size of 13 × 13 mm to first characterize the underlying collagen fiber architecture of the leaflet via an in-house polarized spatial frequency imaging (pSFDI) device [27] (Fig. 1(a)). In brief, the pSFDI results provided the pixel-wise collagen fiber orientation for the effective region of the TV leaflet, which was transformed into a histogram to determine the probability density distribution of the collagen fiber orientation. Since these histograms did not qualitatively follow a normal distribution, the median fiber orientation was used to represent the specimen's preferred fiber orientation (e.g., $θfiber≈−16$ deg for a representative TV leaflet specimen in Fig. 1(a)).

Once the preferred collagen fiber orientation was determined, the leaflet was unmounted from the biaxial testing device, and an 8 × 8 mm specimen was excised from the leaflet such that the $X′$-axis of the specimen approximately aligned with the quantified median collagen fiber orientation (Fig. 1(b)). Finally, the specimen was remounted to the biaxial testing system using a series of five-tined Bio-Rakes, resulting in a 6 × 6 mm effective testing size for the subsequent biaxial mechanical characterizations.

Each of the constant I₁ protocols and constant I₄ protocols was repeated five times at 1.4%/s with a 3-minute rest period between any two protocols to minimize the potential viscoelastic effects. The linear actuator displacements and load cell force values were continuously recorded at 10 Hz throughout the experiment.

The experimental data from the last cycle were used to compute the tissue stretches and the normal components of the first Piola–Kirchhoff stresses using: $λi=(Li+di)/Li$ and $[P]ii=fi/(Lit)$⁠, respectively. Here, subscript i indicates the direction of interest (1: circumferential or $X′$ direction; 2: radial or $Y′$ direction), L_i are the specimen edge lengths following preconditioning, d_i are the actuator displacements, f_i are the load cell force readings, t is the specimen thickness following postpreconditioning, i.e., $t=t0/(λ1PCλ2PC)$⁠, and t₀ is the specimen thickness at the undeformed (reference) configuration as measured via a noncontact laser displacement sensor (Keyence IL-030, Itaska, IL). The first Piola–Kirchhoff stresses were transformed to the Cauchy stresses using $σ=PFT$⁠, and then the computed Cauchy stresses, together with the tissue stretches, were used to calculate the partial derivatives of W in Eq. (5).

The computed values of $W,1$ and $W,4$ were plotted against the experimental I₁ protocols and I₄ for the constant I₄ and constant I₁ protocols, respectively (black dots in Figs. 2(a)–2(c)). These four plots of the experimental raw data revealed noticeable noise for some specimens and loading protocols, which may negatively impact the quality of the fitting results. Hence, quantile regression was used to fit a third-order polynomial function (the initial surrogate function) to the 50% quantile of the experimental data for each constant invariant protocol. Each fit was qualitatively assessed to observe any potential problems associated with the surrogate function fits, such as undesired oscillations in nearly linearly trends, and to determine if an alternative polynomial order is required. The final surrogate functions were resampled at the experimental values of I₁ or I₄ to construct the smooth curves that represent the experimental data (red solid lines in Figs. 2(a)–2(c)) while reducing the noise effects.

We considered two study scenarios to investigate two separate sets of SEDF forms. The first study scenario sought to examine the predictivity for a series of candidate SEDF forms, whereas Study Scenario 2 aimed to evaluate the coupling between the I₁ and the I₄ terms in the SEDF forms.

These three SEDF forms were considered in the baseline comparisons (i.e., study scenario 1) to assess their applicability to the TV leaflets.

where the superscript denotes either the experimental data or the model predictions, and j denotes the current model parameter set. The 1000 parameter sets with the lowest residuals were selected for use in the model fitting procedure (⁠$npop=1000$⁠). In each iteration of the DEO algorithm, the n_pop residuals and the Euclidian distances to the best parameter set were used to update the parameter sets toward the global minimum. Although global search algorithms, such as the adopted DEO algorithm, are robust and widely used in minimization, they may have difficulty finding a unique minimum when the residual search space near the global minimum is flat with small changes in the residual value [33]. The minimization processes were, therefore, repeated n = 50 times for each SEDF form to mitigate the potential influence of this inherent limitation of the DEO global optimization algorithm.

Three assessments were performed for each of the six studied SEDF forms. First, the distributions of the coefficients of determination R² were computed for the four plots of the experimental data to evaluate the model fitting capabilities. Second, the n = 900 sets of model parameters for each SEDF form were used to determine the Pearson's correlation matrix using the corr() function in matlab (MathWorks, Natick, MA), which was subsequently plotted using the imagesc() function to visually assess the potential linear relationships between the model parameters. The presence of potential linear relationships may be indicative of more well-defined search space. Finally, the D-optimality and condition optimality criteria were computed [34,35] for each SEDF form with each of the three TV leaflets. Within the realm of soft tissue biomechanics, these optimality criteria have most notably been employed to improve constitutive model fitting [34] and determine appropriate protocols for the mechanical characterization of myocardium tissue [36] In brief, an in-house matlab program was used to compute the Hessian maxtrix H of Eq. (13) for each leaflet specimen considering the parameter set with the smallest residual value. The D-optimality criterion was then computed as the $det(H)$⁠, whereas the condition criterion was computed as $cond(H)$⁠. The D-optimality criterion should be maximized to improve the model fitting as $det(H)$ represents the inverse of the volume of the confidence region for the model parameters. On the other hand, the condition optimality criterion provides insight into the shape of the confidence region for the model parameters—larger relative values indicate more eccentric (or more nonuniform) confidence regions and smaller relative values for more uniform confidence regions.

The experimental results for the three representative specimens (one per TV leaflet) are shown in Figs. 2(a)–2(c) with the remaining five specimens for each leaflet provided in Appendix (Figs. 9–11). All the three TV leaflets shared similar trends in the experimental data: (i) $W,1$ increasing nonlinearly as I₁ increased; (ii) $W,4$ decreasing nonlinearly with I₁; (iii) $W,1$ decreasing almost linearly with I₄; (iv) $W,4$ increasing almost linearly as an increasing I₄. These observed trends were generally consistent for all the values considered in the constant I₁ or constant I₄ protocols; however, the experimental trends sometimes deviated from these overall observations at the larger values of the constant invariants. For example, the maximum constant I₁ protocol, i.e., $I1=3.4$⁠, produced more nonlinear behaviors for the representative TVAL and TVPL specimens (Figs. 2(a) and 2(b)) and the nonlinear increasing response was found to be more profound for the representative TVAL specimen (Fig. 2(a)), rather than the linearly decreasing trend in $W,1$ versus I₄ observed for the other values of constant I₁. The reasoning for this departure is unclear, but, nevertheless, it underscores the complex relationships between the partial derivatives of W and the invariants that can be determined from the constant invariant-based mechanical characterizations. In addition to these experimental trends, Fig. 2 demonstrates the good-quality fit of the third-order polynomial function to the experimental data.

Three candidate SEDFs were considered in this study scenario—the transversely isotropic version of the polynomial Mooney–Rivlin form (Eq. (7)), an exponential version of the polynomial form (Eq. (8)), and the combination of the polynomial and exponential form (Eq. (9)). The fitting capabilities of these three SEDF forms were first evaluated by using the distributions of the R² values from the nonlinear least-squares fits to the experimental data. Then, the model parameter correlations and the two optimality criteria (i.e., D-optimality and condition optimality) were assessed to examine the model parameter search space and confidence regions.

Representative best fits of the three SEDF forms in Eqs. (7)–(9) for the TVPL specimen in Fig. 2(b) are shown in Figs. 3(a)–3(c), respectively. Of the three SEDF forms, the exponential SEDF in Eq. (8) yielded the best fit to the experimental data of this representative TVPL specimen Fig. 3(b), whereas the polynomial SEDF and the combined polynomial-exponential SEDF had difficulties with fitting to the constant I₄ protocols Figs. 3(a) and 3(c). Moreover, these representative best fits show the inability of the polynomial SEDF form to accurately capture the relationship between $W,1$ and I₁ at the small values of $W,1$ (i.e., the toe region). In contrast, the exponential SEDF, however, can successfully capture this behavior but produces poor fitting at the larger values of $W,1$ for the $I1=3.4$ protocol, which is shared among all three considered SEDF forms.

We further examined the distributions of the R² values as shown in Fig. 4, which provide further insight into the fitting capabilities of the three SEDF forms across different specimens and the corresponding fitting results. These distributions echo the superb performance of the exponential SEDF form for capturing all of the relationships between the partial derivatives of W and the invariants (⁠$R2=0.92$–0.99). Further, they also elucidate the disparity between the polynomial SEDF form (⁠$R2=0.91$–0.97) and the combined polynomial-exponential SEDF form (⁠$R2=0.91$–0.97 versus 0.65–0.95), which is not directly discernible from the representative fits in Fig. 3. Much of the difference in the fitting capabilities of the considered SEDF forms stems from the $W,1$ versus I₄ relationship under the constant I₁ protocols. This relationship, interestingly, holds the most qualitative interprotocol differences as outlined in Sec. 3.1, i.e., starkly different relationships between $W,1$ and I₄ under the larger values of constant I₁. Perhaps the unique coupling of the I₁ and I₄ terms in the exponential SEDF form enables this model to more accurately capture these protocol-specific behaviors. This is further explored in Study Scenario 2 (Sec. 3.3), where the decoupled forms of the exponential SEDF are considered.

The fits for the three SEDF forms were further analyzed to determine the Pearson's correlation matrix (Fig. 5). These color maps, which represent the correlation matrix values, provide further insight into the parameter search space near the global minimum. For example, the polynomial SEDF form had a notably larger mean magnitude for the Pearson's correlation coefficients than the exponential SEDF form (0.51 versus 0.30). A similar trend was shared for the polynomial and exponential components of the combined polynomial-exponential SEDF form (0.51 versus 0.21), underscoring the intrinsic differences in model parameter correlations between the polynomial form and the exponential form. Nevertheless, the larger mean absolute values suggest that the parameter search space is more well-defined for the polynomial SEDF form. These interpretations are primarily qualitative, and the following assessment of the optimality criteria can elucidate more quantitative information of the model parameter's confidence regions.

The D-optimality and condition optimality criteria of the three SEDFs for all of the specimens are summarized in Tables 1 and 2, respectively. These criteria provide insight into the parameter confidence region, such as the volume of the confidence region (D-optimality criterion) or the eccentricity of the confidence region (condition optimality criterion). The D-optimality criterion, which is the most commonly used criterion for design of experiments [34–36], showed that the mixed polynomial-exponential SEDF form had the smallest parameter confidence volume (i.e., the largest D-optimality criterion) for 12 out of 18 specimens and the second-smallest for the remaining 6 specimens. Of the remaining two SEDF forms, the exponential SEDF form had the smallest parameter confidence volume for 6 out of 18 specimens and the second-smallest for the other 12 specimens, while the polynomial SEDF form consistently had the worst D-optimality criterion values (Table 1). On the other hand, the condition optimality criterion, which is used to determine the eccentricity of the parameter confidence region, shows that the polynomial SEDF form often had the least-eccentric confidence region of the three SEDF forms (Table 2). This observation may agree with the uniformly larger Pearson's correlation coefficients highlighted in Sec. 3.2.2 and Fig. 5, and it also follows the prior observation that the exponential SEDF form have an undesired long valley in the model parameter search space near the global minimum.

We further considered the three decoupled forms of the exponential SEDF (Eqs. (10)–(12)) to understand how decoupling the invariant terms within the exponential term influences the model's fitting capabilities, parameter correlations, and optimality criteria.

The exponential SEDF in Eq. (8) provided the best fit to the representative TVPL specimen as demonstrated in Fig. 2(b) (and the same in Fig. 6(a)), whereas the decoupled SEDF forms had remarked difficulties in capturing the experimental relationships (Figs. 6(b)–6(d)). For example, all the three decoupled SEDF forms show the lack of protocol-specific predictions for $W,1$ versus I₁ in the constant I₄ protocols and $W,4$ versus I₄ in the constant I₁ protocols, despite a seemingly good representation of the exponential relationships. On the other hand, the predictions for $W,4$ versus I₁ and $W,1$ versus I₄ kept their protocol-specific capabilities but became linear and unable to capture the nonlinear relationships.

These observed limitations of the three decoupled SEDF forms to fit the experimental data were also noticeable in the distributions of the R² values (Fig. 7). Specifically, the three decoupled SEDF forms showed considerably worse fitting performances for the experimental data for all the TV leaflet specimens (DCPL I–III: $R2=0.79$–0.97 versus 0.74–0.97 versus 0.75–0.97) compared to the exponential SEDF (⁠$R2=0.92$–0.99). It is apparent from these observations that the coupling between the I₁ and I₄ terms is likely critical to the representation of the experimentally quantified behaviors of the TV leaflets. This was alluded to in study scenario 1 as one potential reason for the exponential SEDF form out-performing the polynomial SEDF form when predicting protocol-specific behaviors—the disparity is more profound for the decoupled exponential SEDF forms. Nevertheless, the results show that the decoupled forms may still somewhat capture the experimental trends, which could be beneficial when using traditional planar biaxial mechanical data where the decoupling of the invariants is infeasible.

The color maps of the Pearson's correlation coefficient matrices shown in Fig. 8 qualitatively highlight the reduction in the model parameter correlations as the invariant terms were decoupled. Quantitatively, the median magnitude of the Pearson's correlation coefficients reduced from 0.30 for the exponential SEDF form to 0.15, 0.09, and 0.15 for the progressively decoupled exponential SEDF forms (i.e., DCPL I, II, III), respectively. This observation follows the expectation that the decoupled SEDF forms should have fewer model parameter correlations; however, it is unclear how the reduction in the model parameter correlations alters the residual search space (i.e., the confidence regions), which will be examined next.

The D-optimality and condition optimality criteria for the decoupled SEDF forms provide further quantitative insight into the parameter search space near the global minimum (Tables 3 and 4). Although it was unclear in study scenario 1 whether the exponential SEDF or the combined polynomial-exponential SEDF had the smallest confidence region (i.e., largest D-optimality criterion), the exponential SEDF showed the best criterion values for all the 18 TV leaflet specimens in study scenario 2. Interestingly, the most decoupled SEDF form (DCPL III) held the most second-best D-optimality criterion for 12 specimens followed by DCPL I for 4 specimens and DCPL II for only two specimens. These observations show no apparent trend in the parameter search space as the SEDF is progressively decoupled—the notably worse criterion values underscore the value of the coupled I₁ and I₄ terms in the exponential SEDF. Considering the condition optimality criterion, the exponential SEDF also outperforms the decoupled SEDF forms by having the least eccentric search space (i.e., the smallest condition optimality criterion) for 12 of the 18 specimens and the second-least eccentric for the remaining six specimens. In contrast to the D-optimality, there is a clearly second-best SEDF form candidate: DCPL III having the least eccentric search space for six specimens and the second-least eccentric for the remaining 12 specimens. These observations indicate the superb performance of the exponential SEDF and highlight the benefits of the coupled I₁ and I₄ terms in the exponential SEDF that reduce the confidence region size and eccentricity of the model parameter search space.

Our constant invariant-based mechanical characterization provided the first-of-its-kind experimental quantification of the relationships between the partial derivatives of W and the invariants I₁ and I₄ for the TV leaflets. Similarities in the trends of the experimental data for the n = 18 tested TV specimens suggest that only one SEDF form must be found to capture the mechanical behaviors of all the three TV leaflets (i.e., no leaflet-specific SEDF forms). Interestingly, the results from the present study have stark differences from the prior constant invariant-based characterization by May-Newman and Yin [24] for the mitral valve leaflets, highlighting the tissue-specific nature of the constant invariant characterizations. The mitral valve leaflets showed consistent nonlinearly increasing relationships between the partial derivatives of W and the invariants, whereas our results for the TV leaflets had a mixture of nonlinear/linear increasing and decreasing relationships. Some of these differences could be attributed to how the specimens were prepared in May-Newman and Yin compared to the our investigation. For example, we used our in-house pSFDI device to quantify the initial preferred orientation of the collagan fiber architecture, which then informed the orientation of the excised specimen from the central region of the TV leaflets (see Fig. 1). In contrast, May-Newman and Yin assumed the collagen fiber orientation following the native leaflet's fiber architecture and maintained the entire leaflet entity in their testing (i.e., leaflet, annulus, and chordae insertions). The presence of the annulus and chordae insertions may also influence the observed mechanical behaviors of the mitral valve leaflet. Further, assuming a preferred collagen fiber orientation for the mitral valve leaflets is likely reasonable; however, their post-testing histology showed relatively larger interspecimen and spatially heterogeneous variations in the preferred fiber orientation.

On the other hand, our results herein are more similar to the constant invariant-based mechanical characterization of the left ventricle myocardial tissue by Humphrey et al. [23]. The primary difference is that the TV leaflets appear to have more nonlinear relationships between I₁ and the partial derivatives of W for constant I₄ protocols than the left ventricle myocardial tissue. This may be a consequence of different microstructural properties between the two soft tissues—myocardial tissue is well-known to have highly aligned collagen fibers that vary transmurally [36,37], whereas the TV leaflet collagen fiber architectures potentially have more splay and spatial heterogeneity [38]. Thus, when the fiber stretch is held constant in the constant I₄ protocols, the cross-fiber stress–stretch relationship may be more linear for the myocardial tissue, leading to more apparently linear relationships between the partial derivatives of W and I₁. Future computer simulation-based investigations of the constant invariant-based mechanical characterizations could help elucidate the underlying mechanisms of these differences.

Prior constant invariant-based mechanical characterization studies made qualitative assessment of the experimental results to determine a more suitable SEDF form for the tissue of interest. For example, Humphrey et al. [23] began with a transversely isotropic version of the Mooney–Rivlin SEDF and employed the observed experimental trends (e.g., linearly increasing relationship between $W,1$ and I₁) to tailor the model form for the myocardial tissue. This approach can successfully lead to an SEDF for the tissue of interest; however, it remains to be answered whether another SEDF form would be more favorable in capturing the mechanical behaviors. We, therefore, adopted an alternative approach, in which we considered three candidate SEDF forms and evaluated their respective performance in fitting to the constant invariant-based experimental data. The distributions of R² values in Fig. 4 show the superior performance of the polynomial and exponential SEDFs compared to the combined polynomial-exponential SEDF form. Further analyses of the polynomial and exponential SEDF forms based on the well-established optimality criteria [34,35] revealed that the exponential SEDF had fewer model parameter correlations and larger values of the D-optimality and the condition criteria than the polynomial SEDF form. These observations indicate that the exponential SEDF has a smaller but more eccentric model parameter confidence region. The indication of an eccentric model parameter confidence region is supported by a prior investigation highlighted a long, narrow valley in the model parameter search space near the global minimum for SEDFs with an exponential term [33]. Nevertheless, the D-optimality criterion is commonly taken as the primary consideration to inform experimental design despite both optimality criteria having implications on how the global search algorithms (e.g., DEO) will perform for fitting the experimental data [35].

The fundamental investigations of Humphrey and Yin [23,39] provided the two techniques to determine a SEDF form using experimental data acquired from biaxial mechanical characterizations or constant invariant-based mechanical characterizations. As Humphrey et al. [23] pointed out in the latter constant invariant-based investigation, these two approaches can arrive at a unique SEDF form due to the inability of typical biaxial mechanical characterizations to decouple the contributions from the invariants. This inadequacy of the traditional biaxial mechanical characterizations prevents any coupling of the invariants in calibrating the SEDF form, while the constant invariant-based characterizations allow the invariant coupling being considered. We, thus, sought to explore how the coupling of the I₁ and I₄ invariants in the exponential SEDF may alter the model's fitting performance and the confidence region of the model parameters.

Decoupling of the invariants in the exponential SEDF led to a notable reduction in the fitting performance (DCPL I–III: $R2=0.79$–0.97 versus 0.74–0.97 versus 0.75–0.97) compared to the exponential SEDF (⁠$R2=0.92−0.99$⁠). Looking closer at the representative best fits of the decoupled exponential SEDFs (Fig. 6), one will notice that some of the protocol-specific behaviors (⁠$W,1$ versus I₁ and $W,4$ versus I₄) and the nonlinearity (⁠$W,4$ versus I₁ and $W,1$ versus I₄) cannot be accurately captured by the decoupled SEDF forms. Hence, although the decoupled forms can apparently provide reasonable fits to the experimental data as shown by the median R² values, they provide least accurate predictions of the constant invariant-based protocol features that may be more critical to the representation of the tissue's mechanical behaviors. These observations regarding the decoupled forms are unique to the TV leaflets; however, the changes in nonlinearity and protocol-specific features arising from the decoupling can guide future SEDF selections for other soft biological tissues.

This study is not without limitations. First, although the collagen fiber architecture was characterized using our in-house pSFDI device and careful attention to detail was given to excising the specimen at the median fiber orientation, the collagen fiber architecture may not be perfectly aligned with the Bio-Tester's X-axis. This could be due to the spatial heterogeneity in the collagen fiber architecture when excising the specimen from the leaflet or the small rotational errors when mounting the specimen to the CellScale BioTester. Nevertheless, our experimental setup allows us to ensure that all specimens were nearly aligned with the X-axis, which is key to our assumptions when determining the analytical form of $W,1$ and $W,4$ and facilitating the constant I₄ protocols. Second, the constant invariant technique employed herein is only valid for tissues with minimal or no transmural variations in the collagen fiber architecture. Although previous studies have demonstrated minimal transmural variations for the TV leaflets [10], this precludes the applicability of this method to other biological tissues with more pronounced transmural variations (e.g., myocardial tissue). Third, the tissue stretches were computed by means of the linear actuator displacements, which may not be as accurate as the deformation in the center of the specimen where the boundary or loading effects are minimal. All the tested specimens also had an array of four fiducial markers that were used to compute the tissue stretches following our well-established procedure [9,40]. Our internal study (not presented) that compared the computed stretches between the method of using the linear actuator displacements and the method of the fiducial markers showed negligible differences in the relationships of the partial derivatives of W and the invariants. From this internal study, we elected to use the linear actuator displacements for the tissue stretches, since the constant invariant values are enforced while there may be slight deviations when using the tissue stretches from the fiducial markers. Finally, the DEO algorithm used in this work [32] may possess difficulties finding a unique global minimum for SEDFs with a long valley [33] (e.g., the exponential SEDF form) in the model parameter's search space. We attempted to mitigate this potential limitation by performing numerous minimization procedures (n = 50) for each SEDF form, but the influence of this uncertainty in the global minimum on our results was not systematically evaluated.

One primary future extension of this study is fitting the proposed exponential SEDF (Eq. (8)) to available biaxial mechanical data for the TV leaflets. This can be done in the conventional manner, such as the prior investigation of Khoiy et al. [16], where the SEDF form is fit to the averaged biaxial mechanical data for several specimens. On the other hand, nonparametric bootstrapping [41] can be employed to construct the model parameter likelihoods and compute other meaningful metrics (e.g., the parameter confidence intervals). In addition to the present study using I₁ and I₄, additional constant invariant characterizations can be performed using alternative invariants, such as those proposed by Criscione et al. [42]. These alternative invariants are not covariant unlike the traditional invariants of C [43], which we expect would improve the robustness of the DEO algorithm-based minimization and potentially the model parameter confidence regions or eccentricity. The exponential phenomenological SEDF form can also be enhanced to consider more rich microstructural information, such as the collagen fiber splay about the preferred direction. Its performance can then be compared to the more complex structurally informed constitutive models, such as those proposed by Zhang et al. [21] or Rego and Sacks [22], to determine the further benefits or limitations of this new SEDF form. Ultimately, the exponential SEDF form, together with the optimal model parameters determined in the future extensions, can be employed in in-silico simulations of TV function to provide a computationally tractable representation of the TV leaflet's mechanical behaviors.

This study has provided the first-of-its-kind constant invariant-based mechanical characterization of the TV leaflets. Our experimental results showed the three TV leaflets shared similar relationships between the partial derivatives of W and the invariants I₁ and I₄, suggesting one SEDF form may suffice to represent all three TV leaflets. Three types of SEDFs—a polynomial SEDF, an exponential SEDF, and a combined polynomial-exponential SEDF—were fit to the experimental data to evaluate their potential use for the TV leaflets. The exponential SEDF form displays a superior performance by providing excellent fits to all of the experimental data (⁠$R2=0.92−0.99$⁠) and having a smaller model parameter confidence region, as determined by the D-optimality criterion, than the polynomial SEDF with a similar fitting performance (⁠$R2=0.91−0.97$⁠). Motivated by the prior investigations of Humphrey and Yin that separately led to candidate coupled and decoupled SEDFs for the myocardial tissue, we further explored how decoupling the invariants in the exponential SEDF could alter the model performance. Interestingly, our results showed that the three decoupled exponential SEDF forms showed a worsened fitting performance than the original exponential SEDF by failing to capture many of the protocol-specific and nonlinear relationships in the experimental data. The further analysis revealed that there was no notable improvement in the confidence region shape or eccentricity of the model parameter space near the global minimum. Thus, although the decoupled SEDF forms could possibly represent the TV leaflet's mechanical responses, they have intrinsic limitations that prevent them from accurately representing the leaflet mechanical behaviors, especially under constant invariant-based mechanical characterizations. It is, thus, recommended to use the generalized exponential SEDF for the TV leaflets. Future extensions of this work will determine the optimal model parameters for the TV leaflets, further enabling in-silico simulations of functioning TVs.

• American Heart Association Scientist Development Grant (Grant No. 16SDG27760143; Funder ID: 10.13039/100000968).

• National Science Foundation Graduate Research Fellowship (Grant No. GRF 2019254233; Funder ID: 10.13039/100000001).

• Presbyterian Health Foundation Team Science Grants (Years 1–3) (Funder ID: 10.13039/100001298).

• American Heart Association (Award No. #821298; Funder ID: 10.13039/100000968).

The authors of this paper have no financial or personal relationships with other people or organizations that could inappropriately influence (bias) our work.

The remaining experimental data (black dots) that are not shown in Fig. 2 are presented in Figs. 9–11 for the TV anterior leaflet, posterior leaflet, and septal leaflet, respectively. Note the superb fitting capability of the smooth surrogate function (solid line) to reduce the noise effects in the SEDF form data fitting.