The personalized parameter A = Aopt and the corresponding unloaded reference configuration were obtained at the minimum of fobjective,p. The suction problem, as previously used in the Living Heart Human model (Dassault Systèmes, 2018), was chosen to approximate the inverse mechanics problem within only one simulation, in order to reduce the computational complexity compared to an iterative backward displacement method (Sellier, 2011). A further advantage of this approach is that it can be used without changes to the FEM implementation. The forward simulation with active part included an initial loading with a linear pressure increase from the reference configuration to diastasis and consecutively a cardiac cycle. The initial active material parameters were taken from the Living Heart Human Model (Dassault Systèmes, 2018). The parameters m and b, that determine the relaxation time tr, and thus influence timing and shape of the relaxation (Eq. 14), were adapted manually to the pressure curve, in order to adapt the point in time when zero activation is reached. This adaptation was required to account for the longer duration of the pressure trace due to the initial loading step and a simplified synchronous electrical activation compared to the Living Heart Human Model (Dassault Systèmes, 2018). The manual fit was performed with the RBM fiber field and the resulting parameters were kept constant for all simulations (m = 300, b = −0.38). In accordance to previous works, the parameter Tmax, scaling the active tension (Peirlinck et al., 2019), and the scaling factor n, providing the active sheet contribution (Sack et al., 2018a), were personalized to each model with varying fiber definition. First, the parameter Tmax was adapted using parameter sweeps with manual refinement, such that the end-systolic volume evaluated during early iso-volumetric relaxation at a pressure value of p = 60 mmHg coincided with the value extracted from the data (VES). Next, the parameter n was optimized in a parameter sweep between n = 0.4 and n = 0.9 with a sampling distance of Δn = 0.05. The objective function, defined as the absolute error between the data and the simulation of the maximal left-ventricular shortening during the cardiac cycle with respect to the configuration in diastasis, was minimized. The left-ventricular length was defined as the Euclidean distance between the apex and the center of the base and extracted from each forward simulation. If the minimum was located between two sampling points of the parameter sweep, the optimal parameter nopt was approximated by linear interpolation. The choice of the objective function was motivated by the work of Sack et al. (2018a) that also included this global metric related to longitudinal strain into the fitting process of the active sheet contribution.

We studied the sensitivity of the simulation output depending on the underlying fiber field obtained from the same in-vivo cDTI input, but, computed with different interpolation and mapping methods with varying number of degrees of freedom and varying smoothing strength. To exclude the influence of the material parameters, for each model, simulations with the optimal parameter setting of all other models were conducted and compared. To evaluate the simulation output, the volume over the cardiac cycle was extracted and the pressure-volume relation (pV-loop) was analyzed. Global physiological parameters: end-diastolic volume (EDV), end-systolic volume (ESV), stroke volume (SV), defined as the difference between EDV and ESV, and ejection fraction (EF) were calculated. The circumferential, radial, longitudinal, and fiber strains relative to the initial state in diastasis were extracted. The median over the myocardium during the cardiac cycle and the distribution within the myocardium at peak systole were compared. The twist, defined as the difference in rotation of apex and base, was measured relative to the end-diastolic state.

In Figure 2, the interpolation errors were compared for the different in-vivo fiber models based on the leave-one-out technique with a mid-ventricular short-axis slice of the cDTI data as target (Section 2.2). Figure 2A shows the distributions of the interpolation error. The median interpolation error is the smallest for the HFC method (15.2°), followed by the PDG (18.9°), and POD models (24.2°). The highest error is observed for the RBM (34.0°). To assess the similarity of the fiber fields, the distributions of the mutual differences are shown in Figure 2B. The fiber field obtained from HFC is most similar to the fibers obtained from PGD and vice versa. Both show smaller median angular differences to the fiber field resulting from POD than compared to RBM. Remarkably the difference in angle distribution when comparing HFC to RBM is skewed with a peak at the 25th percentile and a tail of higher difference angles leading to a smaller overall spread. All other distributions have outliers at higher difference values. When comparing the fiber field of the POD model to the fibers from HFC and PGD, the distributions are similar. A slightly higher difference is observed when compared to RBM with higher median and 75th percentile. When comparing all fiber fields to RBM, the POD method shows the largest degree of similarity with the smallest median and 75th percentile.

FIGURE 2. Analysis of the interpolated fiber fields from in-vivo cDTI data obtained with the four interpolation methods: the direct tensor interpolation method using heat flux coordinates (HFC) in red, the low rank model estimation based on Proper Generalized Decomposition (PGD) in violet, the low rank model estimation based on Proper Orthogonal Decomposition (POD) in blue, and the rule-based method adapted to the data (RBM) in green. Subplot (A) shows the absolute angular error distribution and median value for the in-vivo interpolation experiment. The in-vivo cDTI data is interpolated to a previously excluded mid-ventricular short-axis slice and compared to the measured cDTI data. Subplot (B) analyses the distributions of the mutual angular difference between the resulting interpolated fiber fields. The reference fiber field is noted as label on the horizontal axis. Subplots (C,D) depict the distributions of the helix and transverse angle obtained from the data (gray) and the interpolated fiber fields.

Figures 2C,D show the histograms of the characteristic helix and transverse angles from the four interpolated fiber fields compared to the in-vivo data. The fibers obtained from the HFC interpolation show the most similar distribution compared to the data for negative helix angles. The positive helix angles, corresponding to the subendocardial layer, are underestimated by 41.7°/52.8% for the 95-percentile. The helix angle distributions of the fiber fields, resulting from the PGD, POD, and RBM models, underestimate high positive (for the 95th percentile by PGD: 47.3%/POD: 52.9%/RBM: 79.3%) and negative (for the fifth percentile by PGD: 25.8%/POD: 16.5%/RBM: 69.9%) helix angles. The fiber field from RBM has the smallest spread in helix angle. The fibers from the POD model have a bias towards negative helix angles with a negative median of −10.2°. For all other interpolated fiber fields, the negative bias is smaller, with medians between −1.6° and −3.6°. The data has a small imbalance towards positive values with a median of 4.3°. As depicted in Figure 2D, the spread of the transverse angle distribution is smaller for all interpolated fiber fields compared to the data, with the histogram of the RBM model showing the smallest spread. The shape of the distribution is most similar to the data for the HFC-model and PGD-model, and more similar for the POD-model than the RBM-model.

Figure 3 shows the objective function (yellow line) of the joint estimation of the passive scaling parameter A and the unloaded reference configuration (outlined in Section 2.3.4), evaluated in a parameter sweep of A, between Alb = 20 and a variable upper boundary between 35 and 50, adapted to the stiffness of the model. The two contributions to the objective function are shown: the volume error at the initial diastatic pressure during inflation (inflated/gray line) and the volume error in end-diastole (ED/black line). The volume error, corresponding to geometry at the initial diastatic pressure during inflation, decreases for all models with increasing stiffness (increasing A). The volume error for end-diastole shows a minimum in the range of investigated values for A. The resulting optimal passive stiffness scaling factors A are listed in Table 1.

FIGURE 3. Objective function (sum; yellow) of the joint estimation of the passive material parameter and unloaded configuration, together with the single volume error contributions as function of the passive scaling parameter A. The ventricular volume of the simulation results is compared to the volume obtained from cine MRI data at the configuration inflated to the initial pressure (inflated; gray) and at the end-diastolic configuration (ED; black). Each subplot shows the result of the parameter sweep for one interpolated fiber field obtained by the four compared methods: HFC, PGD, POD, and RBM. The vertical lines indicate the optimum.

TABLE 1. Personalized passive and active tissue parameters (A, Tmax) and the active contribution in sheet direction (n) for the four models with different fiber orientation obtained by the four interpolation methods: HFC-model, PGD-model, POD-model, and RBM-model.

The optimized scaling factor of the active tension Tmax is listed in Table 1. The mutual differences show the highest similarity in personalized Tmax between HFC-model and PGD-model, followed by POD-model, and RBM-model.

Figure 4 depicts the personalization of the active tension in sheet direction (n). The gray line depicts the maximal LV shortening obtained from the CINE data. The black line shows the maximal LV shortening in the simulation as a function of the parameter n. The yellow line corresponds to the objective function, defined as the absolute difference of the maximal long axis shortening between simulation and data. The simulated maximal LV shortening decreases monotonically with increasing n. Small values of active sheet contribution, underestimate absolute shortening, while high values lead to an overestimation of absolute shortening. The optimal values of n are listed in Table 1.

FIGURE 4. Objective function (yellow) of the fitting procedure of the parameter that scales the active contribution in sheet direction (n), which is given by the difference between data and simulation outcome of the maximum of the left ventricular shortening during the cardiac cycle with respect to the initial configuration in diastasis. The maximum of the left ventricular shortening as a function of the parameter n is shown, for the data (gray) and the simulation output (black), together with their absolute difference (yellow). The optimized parameter value resulting in the smallest absolute shortening error is indicated by the vertical line. Each subplot shows the results of the parameter estimation for one interpolated fiber field obtained by the four compared methods: HFC, PGD, POD, and RBM. The subtitles indicate the previously adjusted parameter values of the passive (A) and active (Tmax [MPa]) scaling factors.

Figure 5 shows the pressure-volume relations of the forward simulation starting in diastasis, after prior inflation from the unloaded reference configuration. The four subplots correspond to the four parameter configurations obtained from personalization to the four fiber models. For each parameter setting, simulations with all four models, were performed. No sharp isovolumetric contraction is present but rather an indentation of the pv-loop is found. This non-physiological behaviour is related to an imbalance of the increase in pressure and active tension and observed for all models and parameter configurations. The decrease in volume directly after the end-diastolic state is less pronounced for Settings three and four obtained from personalization to the POD-model and RBM-model. In each subplot, the pressure-volume loops are similar for all four models. Their mutual similarity is highest between the HFC-model and the PGD-model and between the POD-model and the RBM-model. The end-diastolic volume (EDV) shows only small differences (below 3 ml) between the models within each setting. The highest EDV is observed for the HFC-model, followed by PGD-model, POD-model, and RBM-model. The EDV for POD-model, and RBM-model are similar with differences <0.2 ml. The order corresponds to the decrease of the personalized stiffness scaling parameter A. For all parameter configurations, the smallest end-systolic volume at the beginning of isovolumetric relaxation is observed for the RBM-model, followed by slightly higher values for POD, HFC and PGD model, with differences ≤7 ml or <13%. The order by increasing end-systolic volume corresponds to the order by increasing active tension scaling factor Tmax, obtained from personalization. The same order applies to the stroke volume and ejection fraction (2). The global physiological parameters, minimal ESV, EDV, SV and EF are listed in Table 2. The maximal variation between the models for each parameter settings is <6 ml (equivalent 8%) for SV and ≤5 percentage points (equivalent 9%) for EF. A high similarity of the global physiological parameters, resulting from the models simulated with their individually optimized parameter setting, is observed with maximal variation of 1.3 ml (equivalent 1.8%) in SV and 0.9 percentage points (equivalent 1.6%) in EF.

FIGURE 5. Simulated pressure-volume relation for the four different parameter settings shown in the four subplots. Each parameter settings was optimized to one interpolated fiber field obtained with the four compared methods: 1) HFC, 2) PGD, 3) POD, 4) RBM. For each parameter configuration (subplot) the pV-loop for all four interpolated fiber fields are depicted. The marker indicates the diastatic state, i.e. the initial state of the simulation.

TABLE 2. Simulated global physiological parameters: minimal end-systolic volume (ESV), end-diastolic volume (EDV), stroke volume (SV), and ejection fraction (EF). Each was evaluated for all four models: HCF, PGD, POD and RBM obtained from simulations with all four parameter configurations.

Figure 6 shows an example time course of circumferential, radial, longitudinal, and fiber strain during the cardiac cycle, starting in diastasis. The strain values are relative to the diastatic state. The presented traces of all fiber models are exemplary for the parameter setting personalized to the HFC-model and are similar in shape for all other parameter settings. During diastole, all strains show only small changes within the first 400 ms after diastasis, corresponding to a plateau in the pressure curve. In late diastole, at the time point t = 400 ms the pressure starts to increase and reaches the end-diastolic pressure (EDP) at t = 514 ms. In parallel the circumferential, longitudinal, and fiber strain increase and the radial strain decreases. The end-diastolic median strains averaged over the models are: circumferential: 0.029, radial: 0.060, longitudinal: 0.038, fiber: 0.023, with a maximal variation between the models smaller than 0.005 for all strains. After end-diastole, during the steep pressure increase, a spike in all strain curves is observed and is most pronounced for the circumferential strain. Similar to the indentation of the pV-loop during iso-volumetric contraction, this behaviour is attributed to the imbalance of active tension and pressure increase (see Figure 5). This spike is followed by rapid decrease in circumferential, longitudinal, and fiber strain, and rapid increase in radial strain, corresponding to systolic contraction. During relaxation and subsequent early diastolic pressure increase, the circumferential, longitudinal, and fiber strain decrease and radial strain increases to the initial configuration with zero strain relative to diastasis.

FIGURE 6. Simulated strain curves over the cardiac cycle, starting from the initial state in diastasis for all models with the four different fiber fields, obtained by the compared interpolation methods (HFC: red, PGD: violet, POD: blue, and RBM: green) and simulated with one exemplary parameter setting optimized for the fiber field interpolated with the HFC method. The subplots show the median of the circumferential, radial, longitudinal, and fiber strain over the myocardium, relative to the inflated state in diastasis.

Figure 7 shows the peak systolic strain distributions, clipped at the 5th and 95th percentile, for the circumferential, radial, longitudinal, and fiber strains relative to the initial diastatic state in the four subplots. Each subplot depicts the simulation results obtained with all four personalized settings, indicated by the numbers on the horizontal axis. For each setting, the distributions of the simulated strains are shown for all models with the different fiber fields (HFC-model in red, PGD-model in violet, POD-model in blue, RBM-model in green). When comparing the results between the models, the same trends exist for all settings, however, an offset of the strain values between the settings is observed. Also the same variation between the models in the spread of the distributions is observed for all settings.

FIGURE 7. Distribution of the simulated strains in systole, cropped at the 5th and 95th percentiles. The subplots correspond to the circumferential (A), radial (B), longitudinal (C), and fiber (D) strain at the end-systolic time-point assigned to the peak radial strain. Each subplot shows four strain distributions (one for each interpolated fiber field) for four different parameter settings (labeled on the x-axis), respectively. The parameter settings are obtained by optimization with one interpolated fiber field: 1) HFC, 2) PGD, 3) POG, 4) RBM. Round markers indicate the median, horizontal lines indicate the fifth,25th,75th, and 95th percentile.

For the circumferential strain, shown in Figure 7A, the spread of the distribution is wider for the HFC-model and the PGD-model than for the POD-model and the RBM-model. The absolute median strain value is the smallest for the PGD-model, followed by smaller absolute strains for the HFC-model, the POD-model, and RBM-model. This order coincides with the order of decreasing value of Tmax obtained during personalization. The same order, with lowest absolute strain for the PGD-model and highest for the RBM-model, is present for the radial strains as shown in Subplot (B). The median strain values are more similar, when comparing the models with their individually personalized setting (HFC-model/PGD-model/POD-model/RBM-model: circumferential: 0.12/-0.12/-0.12/-0.13; radial: 0.57/0.58/0.59/0.56; Supplementary Figure S2). For the longitudinal strain, shown in Subplot (C), the absolute median value of the HFC-model is slightly higher than for the PGD-model. Lower absolute values are obtained with POD-model, and RBM-model. For the personalized settings for each model, the RBM-model also results in the smallest absolute median value (HFC-model/PGD-model/POD-model/RBM-model: longitudinal: 0.15/−0.16/−0.16/−0.14). For the median absolute fiber strain, shown in Subplot (D), the highest value is observed for the RBM-model, followed by the HFC-model, the POD-model, and the PGD-model. The spread between the 25th and 75th percentile is the smallest for the RBM-model, followed by a slightly higher spread for the HFC-model, higher spread for the POD-model, and largest spread for the PGD-model. The fiber strain with the personalized setting, respectively is: HFC-model/PGD-model/POD-model/RBM-model: 0.16/−0.14/−0.14/−0.15.

In Figure 8, the twist relative to end-diastole is shown. Figure 8B depicts the twist over the cardiac cycle starting from the initial diastatic time point. Figure 8A shows the corresponding peak twist in systole corresponding to the time point of 675 ms after diastasis. In Figure 8B, the four subplots correspond to the four parameter settings obtained by personalization to the four models. The vertical lines indicate end-systole and correspond to the settings illustrated in Subplot (A). For all parameter settings the HFC-model results in the highest end-systolic twist (at t = 675 ms), followed by the twist of the PGD-model. The POD-model and RBM-model result in similar twist as the PGD-model for Setting one and lower twist for the other settings at t = 675 ms. When comparing the twist of the models simulated with their personalized settings, the trend remains (twist: HFC-model: 7.4°; PGD-model: 6.4°; POD-model: 5.3°; RBM-model: 5.3°; Supplementary Figure 2S). A variation of the shape of the twist curve over the cardiac cycle and the value of the twist in end-systole is observed between the settings. When analyzing the change of twist over time in Subplot (B), for all settings, first untwist is observed between diastasis and end-diastole. The decrease in twist in end-diastole (the reference frame for twist calculation) occurs parallel to the EDP increase. This is followed by an oscillation at t = 535 ms. Subsequently, 575 ms after diastasis, during the pressure plateau and ejection, a steep increase in twist is observed with clockwise rotation of the base and counter-clockwise rotation of the apex. At the peak of the clockwise rotation of the base, a distinct spike in twist is present for the POD-model and RBM-model and less pronounced for the PGD-model and HFC-model at the t = 625 ms. Subsequently, twist further increases for the PGD-model and HFC-model, until the maximal apical counter-clockwise rotation is reached at the time point 675 ms after diastasis, corresponding to the end of the pressure plateau. This time point is marked by a vertical gray line and evaluated in Subplot (A). For the POD-model and RBM-model, this further end-systolic twist is not observed for Setting four and also not observed for the RBM-model with Setting 3.

FIGURE 8. Twist relative to end-diastole. Subplot (A) shows the twist in end-systole at 675 ms after the initial state in diastasis, for the four parameter settings listed on the horizontal axis and marked by vertical, gray lines. Each setting was obtained by parameter optimisation with one interpolated fiber field (with interpolation methods: HFC (red), PGD (violet), POD (blue), and RBM (green)). For each setting, simulations with each fiber field, obtained by the four interpolation methods, were performed. The resulting twist is indicated by the 4 markers for each parameter setting. In Subplot (B), the twist over the cardiac cycle is shown. Each Subplot (1–4) corresponds to one parameter setting (optimized to: 1: HFC, 2: PGD, 3: POG, 4: RBM), respectively. All plots show the twist for each interpolated fiber field. The gray vertical lines indicate the time point evaluated in Subplot (A).

The personalization of the material parameters, scaling the stiffness of the myocardium (A), maximal active tension (Tmax), and dispersion (n), showed a dependency on the choice of fiber interpolation method. The parameter values differ significantly between the microstructural models with a maximal increase of 31% for A, 53.5% for Tmax, and 11.8% for n. This dependence of stiffness and stress on the fiber orientation has been previously observed by Wang et al. (2013). The variation in parameter values for different interpolation methods hinders comparability of model parameters to physical quantities measured in laboratory experiments (e.g. stiffness or active tension) and moreover requires consideration when comparing modelling studies (e.g. active fiber stress). The parameters not only reflect the physical tissue property but also act as auxiliary parameters to compensate for inaccuracies in the underlying fiber representation. The similarity of the parameters A and Tmax between the HFC-model and the PGD-model is correlated to the highest similarity in the interpolated fiber field providing the smallest interpolation errors compared to the data. The POD-model and RBM-model, both result in higher difference angles compared to the data and lower parameter values. This trend and separation into two groups is the same as observed in the simulation outputs when keeping the parameter setting constant. ESV, median circumferential, and radial strain follow the trend of the resulting personalized Tmax. EDV follows the resulting personalized stiffness scaling A.

The analysis for all parameter settings reveals a systematic bias due to the underlying fiber interpolation method. However, in practice the parameters are optimized to the fiber field. This on the one hand results in differences in the personalized parameters and on the other hand in remaining effects on the simulation output. After personalization of each model individually, the differences in median strain values were reduced to: circumferential: ≤ 10.1%; radial: ≤ 3.9%; longitudinal: ≤ 13.8%; fiber: ≤ 17.7%. Differences in the distribution of the fiber strain were observed. A smaller spread of the distribution is present for the HFC-model and the RBM-model compared to the PGD-model and POD-model, resulting from the smoother microstructure representation. Due to the personalization with objective functions including the EDV and ESV (evaluated during relaxation at a pressure of 60 mmHg), the remaining variation in SV and EF are small (SV: 1.3ml/1.8%; EF: 0.9 percentage points/1.6%) and can be attributed to minor volume changes during relaxation (Supplementary Figure S2). Compared to literature values all models resulted in circumferential, longitudinal, and fiber strains within a physiological range. (The average median strains were: circumferential: 0.1475, radial: 0.575, longitudinal: 0.153, fiber: 0.148.) Literature values of global strains from porcine studies are: circumferential: 0.17/−0.14/−0.14/−0.15; longitudinal: 0.1675/−0.17/−0.11/−0.13; radial: 0.435/0.65/0.21/0.39 from: (Stoeck et al., 2021): (cine SSFP, mean value)/(Berberoğlu et al., 2022) (cine SSFP)/(Ferreira et al., 2018) (DENSE MRI, mean value)/(Verzhbinsky et al., 2020) (DENSE MRI, median value)), and fiber strain: 0.14 from (Verzhbinsky et al., 2020) (combined DENSE and cDTI). It is noted that a small underestimation compared to the literature values is expected due the difference in the reference configuration: The strains were referred to mid-diastole for this simulation study and end-diastole for the data. The simulated radial strain is higher than three out of four literature values, thus, it might be overestimated in the simulation. However, the literature values of radial strain are subject to a high variation, corresponding to high uncertainty in the data. Furthermore, different approaches for strain estimations, e.g., 2D vs. 3D, are found in literature.

All models underestimated twist (HFC-model/PGD-model/POD-model/RBM-model: twist = 7.4°/6.4°/5.3°/5.3°, corresponding to a maximal torsion of 0.15°/mm, obtained with the HFC-model), compared to literature values measured by 3D tagging with CSPAMM (Rutz et al., 2008) in pigs: twist: 11.05°(Berberoğlu et al., 2022); torsion: 0.27°/mm (Stoeck et al., 2021) averaged over all measurements and all healthy animals. In contrast to twist, torsion excludes the bias of heart size (Young and Cowan, 2012) and is calculated by the ratio of the twist and the end-diastolic ventricular length. The HFC-model led to the highest end-systolic twist, closest to physiological values, followed by the PGD-model and lower values for the POD-model and the RBM-model. This trend agrees well with the interpolation error of the models and unlike other physiological readouts is not compensated by personalization based on volume and LV length data. Similarly, higher twist with a more realistic representations of microstructure has been observed by Gil et al. (2019). Comparing the temporal course of twist during the cardiac cycle to a healthy physiological curve shape (Rüssel et al., 2009; Omar et al., 2015), the twist obtained with the HFC-model and the PGD-model have a more physiological shape, showing less pronounced non-physiological oscillations and spikes during isovolumetric contraction and early systole than the POD-model and the RBM-model. For both, HFC-model and the PGD-model twist increased during ejection after the maximum counter-clockwise rotation of the base and a peak in twist at end-systole is observed as found in physiology. Contrary, the POD-model and RBM-model showed a drop in twist after the early-diastolic spike.

The tensor interpolation method (HFC) corresponds to the smallest interpolation error. It exploits the high spatial coherence of the tensor field in shape adapted coordinates. For optimal performance this method requires the adaptation of the weighting matrix H for the anisotropic Gaussian kernel according to the spatial resolution and coverage of the input data. The two low-rank models (PGD and POD) are based on basis functions that were extracted from high-resolution data. This prior information results in a lower degree of flexibility to represent the input data compared to the direct tensor interpolation. The PGD model uses independent basis function in each spatial direction and thus better adapts to the input data compared to the POD model. This results in a smaller interpolation error with the PGD model than with the POD model. The rule-based method is a simple linear model with only four degrees of freedom. Thus, it over-smooths the fiber field and results in the highest interpolation error. The computational cost to generate the personalized fiber fields was: HFC: 21s, PGD: 34s, POD: 18s, RBM: 2s. Times were measured on an eight-core Intel Core i7-10700K, 3.79 GHz desktop computer. The times are three orders of magnitude smaller than the simulation time: HFC: 4 h 39 min, PGD: 4 h 42 min, POD: 4 h 32 min, RBM: 4 h 38 min. The standard deviation of 4.2 min between the simulations for the four models is negligible compared to the total simulation time and no differences in stability were present.

This study has shown the feasibility of personalization of microstructure in a biomechanical model based on a pre-clinical animal experiment with cDTI data including nine short-axis slices. In a more clinical setting, three slices are currently acquired (Khalique et al., 2020, 2018; Gotschy et al., 2021). While full coverage and isotropic spatial resolution are desirable for biomechanical modelling of the heart, cDTI in clinical practice suffers from long scan duration (Nguyen et al., 2020), which is being addressed as part of ongoing imaging research (Nguyen et al., 2021). Due to a steep, non-linear error increase for less than five short-axis slices for all interpolation methods Stimm et al. (2022) (average error increase from nine to three slices in an ex-vivo experiment: HFC: 4.6°/PGD: 10.6°/POD: 5.2°/RBM: 0.2°), further development in clinical data acquisition is required. When applying the interpolation approaches presented in this study to clinical data, an additional uncertainty in the simulation results would be present due to this error increase for less than five short-axis slices. However, the tensor interpolation method (HFC) remains the interpolation method with the smallest interpolation error also with three input slices (with an advantage of a 6.7°smaller error compared to the rule-based method (RBM) on average (Stimm et al., 2022)). In this study, we showed that a systematic bias in the simulation results was introduced by the interpolation method and that a correlation of the trend with the increase in interpolation error exists. Combining these two observations of 1) the smaller interpolation error that was observed for the HFC method and 2) the systematic bias of the simulation outputs, we expect a reduction in the bias of the simulation results when using the HFC method also with clinical data. However, the overall higher interpolation error increases the uncertainty compared to the pre-clinical setting and the difference in interpolation performance between the methods deceases in the clinical setting thus the bias is expected to be smaller. To reduce the uncertainty for patient-specific simulations based on clinical in-vivo cDTI data, advancements in clinical DTI that would enable a minimum of five input slices are required.

In-vivo cDTI studies have observed a reorientation of microstructure in pathology, such as dilated cardiomyophathy (DCM) (von Deuster et al., 2016a), hypertrophic cardiomyopathy (HCM) (Ferreira et al., 2014; Das et al., 2022) and aortic stenosis (AS) (Gotschy et al., 2021). Patient-specific modelling based on a patient-specific microstructure enables to perform modelling studies investigating the link between structure and function of the heart. This can improve the understanding of the mechanical conditions that lead to structural remodelling, and consequently disease progression. Remodelling might onset in an early disease state (Gotschy et al., 2021) and therefore early biomarkers can be revealed by investigation of remodelling. Together with future advances in both in-vivo cDTI and patient-specific modelling based on clinical data predictive modelling is a future goal.

Systematic trends of the simulation results were observed between the fiber models, when using the same model parameters. The correlation of this systematic bias with the interpolation performance suggests that a more realistic representation of in-vivo microstructure affects the simulation output and therefore reduces the uncertainty in the simulation results. This suggests that using the interpolation method with the lowest interpolation method (HFC tensor interpolation) leads to the simulation results with the lowest uncertainty. The higher twist observed with the HFC-model, remaining present after personalization, further supports that models profit from more realistic representation of microstructure.

In this study, we have investigated the dependence of model behaviour on interpolation errors and compared the results to cohort values. The study concentrated on the isolated bias introduced by the representation method to personalize the model using in-vivo cDTI data. The direct comparison to patient-specific physiological measurements is currently not feasible. Results would be biased by model inaccuracy, due to simplifications, such as neglecting the influence of the right ventricle and data uncertainty (e.g. for radial strain and consequently fiber strain estimation) on the individual basis. Further, the improvement compared to a generic fiber representation was not evaluated, which would require a cohort study with enough cases to represent the fiber variability across the population. In pathology, microstructural remodelling is observed (Ferreira et al., 2014; von Deuster et al., 2016a; Gotschy et al., 2021), consequently the variability of structure increases, rendering patient-specific microstructure more important.

The personalized microstructure, based on in-vivo cDTI data, is subject-specific, but both, the data acquisition and the interpolation techniques introduce errors. The error of the first eigenvector of the in-vivo diffusion tensor, induced by the measurement, can be estimated with the cone of uncertainty (Aliotta et al., 2018). Aliotta et al. (2018) estimated an uncertainty of 15.5°(human data). We estimated an uncertainty of 11.3°in our previous work (Stimm et al., 2022) for porcine data, which was equivalently used in this study. The interpolation error calculated on a mid-ventricular slice was: HFC: 15.2°, PDG: 18.9°, POD: 24.2°, RBM: 34.0°. Sensitivity studies by Geerts et al. (2003); Pluijmert et al. (2017) have found that a difference angle of 8°already influences the simulation output. Consequently, the error minimization by using a HFC representation compared to a RBM representation has an effect on the simulation results, as was confirmed in this study. The analyzed interpolation error reflects the difference angle to the measurement, that is subject to noise and imaging artifacts. Consequently, denoising by smoothing or by exploiting prior information included in the fiber models, also contributes to the measured angular difference. The interpolation performance of the models has been evaluated in detail in (Stimm et al., 2022), presenting an average interpolation error of: HFC: 13.7° ± 1.1°, PGD: 17.8° ± 5.2°, POD: 19.0° ± 3.3°, RBM: 24.8° ± 6.0°in-vivo. Thus, the interpolation error of the underlying case is slightly above the average error for all methods. The number of degrees-of-freedom of the data-based models (PGD and POD) was adapted to ex-vivo data Stimm et al. (2021). Optimization in an in-vivo study might increase denoising and result in smoother microstructure representation, thereby, potentially reducing the spread of the strain distributions. All methods were adapted based on all data points, however, the endo- and epicardial boundaries were subject to partial-volume effects, resulting in higher uncertainty of these data points. Further, segmentation errors, in the proximity of the papillary muscles at the endocardium, might have influenced the errors at the boundaries. This resulted in wider spread of the helix and transverse angles in the data. This spread is reduced by denoising effects of the interpolation. Although, over-smoothing has the same effect, it reduces the spread of the helix and transverse angle distributions of the interpolated fiber representations too radically. Adaptation of the fitting such that it only uses data in areas not prone to partial voluming and segmentation errors, might reduce the representation error, especially for the RBM-method with only four degrees of freedom. While the HFC-model resulted in a globally more heterogeneous microstructure than the RBM-model, both result in smoother fiber representations than the PGD-model and POD-model. This might be a reason for the more narrow distribution of the fiber strain found in this study.

The patient-specific model used in this work contains simplifications, influencing the accuracy of the simulations. While right-ventricular deformation and pressure influence left-ventricular mechanics (Palit et al., 2015), a LV only model was used in this study because in-vivo microstructure is currently not available for the thinner right-ventricular wall. Further, the surrounding of the heart was modeled by a boundary condition constraining the base as proposed by Peirlinck et al. (2019). Adding a ribcage and pericardial boundary condition can improve model accuracy (Pfaller et al., 2019; Ponnaluri et al., 2019; Strocchi et al., 2020; Hadjicharalambous et al., 2021). However, both, neglecting influences of the RV and the surrounding are expected to result in a consistent offset of the simulation results, not affecting the isolated influence of the microstructure representation, analyzed in this study. Due to a steeper increase in ventricular pressure compared to the active tension, during isovolumetric contraction, a non-physiological change in volume has been observed. The pressure trace was generated based on subject-specific in-vivo measurements and is inconsistent to the cosine-shape active tension resulting from the active model of the living heart framework. The strong hemodynamic constraint of a fixed pressure trace leads to a challenging, if even possible, active parameter fitting problem. An alternative, would be a lumped-parameter model, however, resulting in an idealized pressure trace not directly representing the measured pressure. Another data-driven alternative would be to incorporate the active parameter as Lagrange multipliers in the model implementation, while including both, pressure and volume measurements as model inputs (Asner et al., 2016; Miller et al., 2021). This option is, however, not part of the Living Heart Human model. Therefore, the results during isovolumetric contraction were not evaluated in this study and physiological parameters were calculated from the volume corresponding to the EDP and not from that corresponding to maximal volume.

The twist, given by the relative rotation of apex and base, was evaluated and follows a physiological behaviour. Nevertheless, a remaining net counter-clockwise rotation of the left ventricle was observed at the end of the simulation of one cardiac cycle, related to the missing constraint of the moments. Constraining the average moments of the base to zero on the contrary would prevent basal rotation and thus affect twist. A potential way to compensate the net rotation would be to constrain the average moments to the average moments of the data (Asner et al., 2017).

Linear tetrahedral elements were used in this study to simulate the nearly-incompressible material. Despite the small element size counteracting potential inaccuracies, this might introduce locking and thus an artificial stiffening and torsional rigidity (Maas et al., 2016). Therefore, a potential way to mitigate the underestimation of twist for all compared models could be the use of quadratic tetrahedral elements.

The data analyzed in this study is subject to the following licenses/restrictions: The data was collected from previous publications and is available upon reasonable request to the corresponding author. Requests to access these datasets should be directed to stoeck@biomed.ee.ethz.ch.

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2022.1042537/full#supplementary-material

Supplementary Figure S1 | Pressure trace of the forward simulation during one cardiac cycle, starting at the initial diastatic state. This pressure corresponds to the loading condition on the endocardial surface applied upon linear inflation.

Supplementary Figure S2 | Personalized model outputs for the HFC-model (with Setting 1) in red, PGD-model (with Setting 2) in violet, POD-model (with Setting 3) in blue, and RBM-model (with Setting 4) in green. (A) shows the circumferential, radial, longitudinal, and fiber strain distributions (similar to Figure 7), (B) shows the pV-loops, with a marker indicating the initial state (similar to Figure 5), and (C) shows the twist over the cardiac cycle. The vertical black line in (C) marks the end-systolic time-point where peak twist is evaluated (similar to Figure 8).

Aliotta E., Moulin K., Magrath P., Ennis D. B. (2018). Quantifying precision in cardiac diffusion tensor imaging with second-order motion-compensated convex optimized diffusion encoding. Magn. Reson. Med. 80, 1074–1087. doi:10.1002/mrm.27107

PubMed Abstract | CrossRef Full Text | Google Scholar

Asner L., Hadjicharalambous M., Chabiniok R., Peressutti D., Sammut E., Wong J., et al. (2017). Patient-specific modeling for left ventricular mechanics using data-driven boundary energies. Comput. Methods Appl. Mech. Eng. 314, 269–295. doi:10.1016/j.cma.2016.08.002

CrossRef Full Text | Google Scholar

Asner L., Hadjicharalambous M., Chabiniok R., Peresutti D., Sammut E., Wong J., et al. (2016). Estimation of passive and active properties in the human heart using 3D tagged MRI. Biomech. Model. Mechanobiol. 15, 1121–1139. doi:10.1007/s10237-015-0748-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Baillargeon B., Rebelo N., Fox D. D., Taylor R. L., Kuhl E. (2014). The living heart project: A robust and integrative simulator for human heart function. Eur. J. Mech. A Solids 48, 38–47. doi:10.1016/j.euromechsol.2014.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Barbarotta L., Bovendeerd P. H. M., Ennis D W. V., Perotti L. E. (2021). “A computational approach on sensitivity of left ventricular wall strains to fiber orientation,” in Functional imaging and modeling of the heart. FIMH 2021. Lecture notes in computer science (Cham: Springer), 296–304. doi:10.1007/978-3-030-78710-3_29

CrossRef Full Text | Google Scholar

Bayer J. D., Beaumont J., Krol A. (2005). Laplace–dirichlet energy field specification for deformable models. An FEM approach to active contour fitting. Ann. Biomed. Eng. 33, 1175–1186. doi:10.1007/s10439-005-5624-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Bayer J. D., Blake R. C., Plank G., Trayanova N. A. (2012). A novel rule-based algorithm for assigning myocardial fiber orientation to computational heart models. Ann. Biomed. Eng. 40, 2243–2254. doi:10.1007/s10439-012-0593-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Bayer J., Prassl A. J., Pashaei A., Gomez J. F., Frontera A., Neic A., et al. (2018). Universal ventricular coordinates: A generic framework for describing position within the heart and transferring data. Med. Image Anal. 45, 83–93. doi:10.1016/j.media.2018.01.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Berberoğlu E., Stoeck C., Kozerke S., Genet M. (2022). Quantification of left ventricular strain and torsion by joint analysis of 3D tagging and cine MR images. HAL pre-print. hal-03604931.

Google Scholar

Berberoğlu E., Stoeck C. T., Moireau P., Kozerke S., Genet M. (2021). In-silico study of accuracy and precision of left-ventricular strain quantification from 3D tagged MRI. PLOS ONE 16, e0258965. doi:10.1371/journal.pone.0258965

PubMed Abstract | CrossRef Full Text | Google Scholar

Beyar R., Sideman S. (1984). A computer study of the left ventricular performance based on fiber structure, sarcomere dynamics, and transmural electrical propagation velocity. Circ. Res. 55, 358–375. doi:10.1161/01.RES.55.3.358

PubMed Abstract | CrossRef Full Text | Google Scholar

Buljak V., Maier G. (2011). Proper orthogonal decomposition and radial basis functions in material characterization based on instrumented indentation. Eng. Struct. 33, 492–501. doi:10.1016/j.engstruct.2010.11.006

CrossRef Full Text | Google Scholar

Buoso S., Manzoni A., Alkadhi H., Plass A., Quarteroni A., Kurtcuoglu V. (2019). Reduced-order modeling of blood flow for noninvasive functional evaluation of coronary artery disease. Biomech. Model. Mechanobiol. 18, 1867–1881. doi:10.1007/s10237-019-01182-w

PubMed Abstract | CrossRef Full Text | Google Scholar