Causes of altered ventricular mechanics in hypertrophic cardiomyopathy: an in-silico study

Background Hypertrophic cardiomyopathy (HCM) is typically caused by mutations in sarcomeric genes leading to cardiomyocyte disarray, replacement fibrosis, impaired contractility, and elevated filling pressures. These varying tissue properties are associated with certain strain patterns that may allow to establish a diagnosis by means of non-invasive imaging without the necessity of harmful myocardial biopsies or contrast agent application. With a numerical study, we aim to answer: how the variability in each of these mechanisms contributes to altered mechanics of the left ventricle (LV) and if the deformation obtained in in-silico experiments is comparable to values reported from clinical measurements. Methods We conducted an in-silico sensitivity study on physiological and pathological mechanisms potentially underlying the clinical HCM phenotype. The deformation of the four-chamber heart models was simulated using a finite-element mechanical solver with a sliding boundary condition to mimic the tissue surrounding the heart. Furthermore, a closed-loop circulatory model delivered the pressure values acting on the endocardium. Deformation measures and mechanical behavior of the heart models were evaluated globally and regionally. Results Hypertrophy of the LV affected the course of strain, strain rate, and wall thickening—the root-mean-squared difference of the wall thickening between control (mean thickness 10 mm) and hypertrophic geometries (17 mm) was >10%. A reduction of active force development by 40% led to less overall deformation: maximal radial strain reduced from 26 to 21%. A fivefold increase in tissue stiffness caused a more homogeneous distribution of the strain values among 17 heart segments. Fiber disarray led to minor changes in the circumferential and radial strain. A combination of pathological mechanisms led to reduced and slower deformation of the LV and halved the longitudinal shortening of the LA. Conclusions This study uses a computer model to determine the changes in LV deformation caused by pathological mechanisms that are presumed to underlay HCM. This knowledge can complement imaging-derived information to obtain a more accurate diagnosis of HCM. Supplementary Information The online version contains supplementary material available at 10.1186/s12938-021-00900-9.


Background
Hypertrophic cardiomyopathy (HCM) is a relatively common inherited disorder, with a prevalence of 1:500, which develops in the absence of an identifiable cause [1,2]. There are several phenotypes of HCM, depending on the localization and distribution of hypertrophy in the heart -asymmetric, symmetric/concentric, apical or mid-ventricular obstruction [3]. HCM results in an increased ratio of wall to lumen volume, which can be diagnosed by echocardiographic or magnetic resonance assessment of left ventricular anatomy [1]. Besides this morphological modification, further abnormalities are underlying the HCM phenotype. In HCM hearts, fibrosis and myocardial cell disarray can be present and might have evolved for years before the onset of symptoms [1]. The disarray of the cells can be quantified by fractional anisotropy (FA)-a measure obtained by diffusion tensor magnetic resonance imaging (DT-MRI) or shear wave imaging (SWI). Ariga et al. [4] reported reduced FA in HCM patients compared to control subjects, measured by DT-MRI. Villemain et al. [5] reported similar findings in pediatric HCM patients using SWI. Further structural abnormalities in HCM were detected by SWI on an organ level: passive ventricular stiffness was significantly higher in HCM compared to the control group [5]. In HCM, increased stiffness on the organ level could not be explained by an alteration in the viscoelastic properties of the cardiac myocytes, since the passive stiffness of prepared HCM cardiac myocytes was measured to be the same as healthy donor myocytes [6]. The stiffer tissue behaviour might be due to further factors such as cell disarray or tissue fibrosis. Furthermore, the maximal active force was markedly lower in HCM myocytes than in donor myocytes [6]. In clinical routine, it is not possible to measure active force development of the myocytes in-vivo. Furthermore, the application of SWI for stiffness measurement is limited to the entire ventricle and might be not applicable for all patients [7]. A reconstruction of the myocardial cell orientation with DT-MRI is very time consuming and delivers limited anatomical coverage of the ventricle [8]. These limitations of the available imaging modalities make it impossible or at least cumbersome to identify abnormalities underlying the HCM phenotype in clinical routine. Nevertheless, the consequences of these structural changes can be measured and quantified to provide a basis to diagnose HCM. This diagnosis is often based on echocardiographic assessment of the systolic and diastolic function of the left ventricle (LV) [1] and parameters derived from tissue imaging (strain and strain rate) [9]. Furthermore, MRI can quantify heart motion and function by cine imaging, which enables LV wall thickness calculation. Tissue phase mapping and feature tracking provide LV radial, circumferential and longitudinal myocardial velocity time courses, as well as global and segmental systolic and diastolic peak velocities [10]. The left atrial longitudinal strain is measured as well to inform HCM diagnosis [11]. Such precise assessment of the cardiac function enables the quantification of altered mechanics in HCM patients compared to healthy volunteers. In this work, we conduct in-silico experiments to identify potential underlying causes of altered ventricular mechanics observed in HCM patients. The numerical heart simulator includes models of active force development, passive stiffness, the circulatory system and appropriate boundary conditions to conduct a sensitivity study. We define models capturing different pathological mechanisms in a virtual heart and explore different combinations of these to analyse their effect on ventricular mechanics. The identification of distinct underlying abnormalities leading to the altered mechanical behaviour of HCM hearts complementary to imaging could be valuable information for clinicians on the way to clearer and faster diagnoses. It can provide directions to differentiate HCM from other cardiac conditions, in which thickened walls are present. Moreover, it could help to separate healthy hearts from HCM genotype-positive but phenotype-negative hearts. This study aims at establishing cause-effect relationships between pathological mechanisms in HCM hearts and their effect on ventricular mechanics.

Results
Altered Mechanics Due to the Wall Thickness of LV We quantified isolated changes of the LV WT by comparing the deformation for different geometries-control geometry (Case 1), HCM 1 (Case 5) and HCM 2 (Case 6). During the diastolic period, the RMSD of the global wall thickening between each two out of the three cases was less than 1.3 %. During the systolic period, the RMSD of the global wall thickening between Case 1 and both hypertrophic cases was higher than 10.8 % and the RMSD between Case 5 and Case 6 was 2.6 %. The RMSD for all cases and all metrics are provided in the Supplementary Material ( Figure 1 and 2). The increase of the regional WT for the hypertrophic geometries (Case 5 and Case 6) during systole, for all 17 segments, was faster compared to the one of the control geometry (Case 1). Between Case 5 and Case 6, no marked difference was observed ( Figure 1). Wall thickening at end-systole (ES) was in the same range in all three cases (Case 1: between 18.1 % and 50.0 %, Case 5: between 18.4 % and 48.1 %, and Case 6: between 18.1 % and 51.3 %). Nevertheless, the distribution among the segments changed-the thickening of the basal segments (1-6) increased as the initial WT increased. Figure 2 shows the ES distribution of the wall thickening.   Figure 1 The time courses of the regional wall thickening for Case 1, 5 and 6 (left, middle, right, respectively). In each plot, the first vertical line (at 0.17 s) indicates ED and the second line (at 0.5 s) ES. Case 1: control case (initial WT = 10 mm); Case 5: hypertrophic geometry (15 mm); Case 6: hypertrophic geometry (17 mm).
The circumferential and radial strain and the strain rate (both global and segmental) for the control geometry differed from those of the hypertrophic geometries. During the systole, the RMSD of the global strain between Case 1 and each of the cases of hypertrophic geometries was between 1.4 % and 3.6 %. An increase of radial strains and a decrease of circumferential strains for Case 1 was observed during the entire systole, while those for the hypertrophic geometries occurred during the first half of the systolic period ( Figure 3 and Figure 4). This can be derived also from the strain rate-they are higher for the radial direction   Figure 3 The time courses of the regional circumferential strain for Case 1, 5 and 6 (left, middle, right, respectively). In each plot, the first vertical line (at 0.17 s) indicates ED and the second line (at 0.5 s) ES. Case 1: control case (initial WT = 10 mm); Case 5: hypertrophic geometry (15 mm); Case 6: hypertrophic geometry (17 mm).  The time courses of the regional radial strain for Case 1, 5 and 6 (left, middle, right, respectively). In each plot, the first vertical line (at 0.17 s) indicates ED and the second line (at 0.5 s) ES. Case 1: control case (initial WT = 10 mm); Case 5: hypertrophic geometry (15 mm); Case 6: hypertrophic geometry (17 mm).
and lower for the circumferential direction in the mid systole for the hypertrophic geometries compared to those in the control case. An initially similar strain rate is available for all three cases during the systole, with amplitude 200 %/s for the radial strain rate and around 90 %/s for the circumferential strain rate ( Figure 5, left).
The end-diastolic (ED) regional strains in all three directions were comparable in all three cases.  The ES longitudinal strain indicated less shortening of the tissue since its minimum over the segments became higher as the WT increased-from -11.0 % for Case 1 to around -6.5 % for Case 5 and Case 6. Furthermore, the circumferential and radial strain at ES indicated slightly less shortening of the tissue (Figure 3 and Figure 4). In Figure 6 the distribution of the strains at ES are visualized for each local direction (longitudinal, circumferential and radial). The velocities in all local directions for the control case were differing from those of the hypertrophic geometries during the entire systole (RMSD ranged from 0.01 m/s to 0.013 m/s). The velocities for the control case were between 0.04 m/s and 0.05 m/s during the beginning and the middle of the systole, while the velocities for the hypertrophic geometries were high at the beginning of the systole (around 0.06 m/s) and decreased quickly to 0.02 m/s until the middle of the systole. The ES velocities were comparable in all three cases. During the diastolic period, the RSMD between the control case and each of the cases of the hypertrophic geometries ranged from 0.002 m/s to 0.003 m/s. The time course of the velocities in all three local directions for Case 1 and Case 6 are visualized in Figure 5, on the right. The visualization of Case 5 was omitted, since it was comparable to Case 6. The longitudinal strain of the left atria increased slower during the systolic period for Case 1 compared to Case 5 and Case 6. The maximal longitudinal strain of the LA occurred at ES and it was 20 % for all three cases.
Summary: During the systolic period, the time course of the wall thickening and strain is altered between the control geometry and the hypertrophic geometries in all AHA segments. At ES, the circumferential and radial strain had only minor differences between the control geometry and the hypertrophic geometries, while the longitudinal strain indicated less shortening as the WT increased. The myocardial velocities for both hypertrophic geometries are increased during the systolic period compared to the control case. We did not observe any major differences in the measures between both hypertrophic geometries and therefore, in the following, we omit the less hypertrophic geometry (HCM 1) for further comparisons.  Figure 6 Bull's-eye displays for Case 1, 2, 6, 8, 10 and 13 showing the longitudinal, circumferential and radial strain at ES (first, second and third column, respectively). Each row corresponds to one case. Case 1: control case; Case 2: increased stiffness; Case 6: hypertrophic geometry (17 mm); Case 8: hypertrophic geometry (17 mm), decreased active force; Case 10: hypertrophic geometry, fiber disarray; Case 13: virtual HCM case (all pathologies included).

Altered Mechanics Due to the Stiffness of LV
We quantified isolated changes in the passive force by comparing the deformation of the simulation with control stiffness (Case 1) to the one with increased stiffness (Case 2), both with the control geometry. At the ED, the regional wall thickening was comparable between Case 1 and Case 2.
At the ES, it differed-the extent of the thickening values was reduced (Case 1: 18.7 %-50.0 % and Case 2: 16.9 %-35.8 % ), since the maximum decreased more than the minimum did, when the stiffness of the tissue was increased. In particular, the wall thickening of the segments in the free wall of the LV (5, 6, 11 and 12) were diminished compared to the one in the septal segments (2, 3, 8 and 9). This led to more homogeneous distribution of the wall thickening among all segments, when the stiffness of the tissue was increased (Figure 2, Case 1 vs. Case 2). Similar to the wall thickening, the regional strain in all directions (longitudinal, circumferential and radial) at the ES had reduced extent of the values when the stiffness was increased. Therefore, a more homogeneous distribution of the strain among the all segments at ES was present for Case 2 compared to Case 1 ( Figure 6). The global radial strain rate at the beginning of the systolic period was halved when the stiffness was increased (Figure 7, left). During the entire systolic period, the global velocities in all directions reduced when the stiffness was increased. During the diastolic period, the maximal velocities in all local directions slightly increased as the stiffness increased ( Figure 7, right). Furthermore, the velocities decreased quicker as the stiffness was increased. While for the control case (Case 1) the maximal velocities during the diastole are similar compared to the velocities during the systole, increased tissue stiffness (Case 2) led to higher velocities during the diastole compared to the velocities during the systole. The longitudinal strain of the left atria strongly decreased for the case of increased stiffness: 20 % for Case 1 to 10 % for Case 2.
A comparison the deformation of the simulation with control stiffness (Case 6) and the one with increased stiffness (Case 7), both with the hypertrophic geometry (HCM 2), confirmed these results. Summary: An increased stiffness of the tissue of the LV led to more homogeneous wall thickening among the AHA segments as well as equalized the strains at the ES.
The strain rates and velocities were reduced during the systole, especially visible at the beginning of the systole. In contrast, the maximal velocities during the diastole increased when the stiffness was increased and reduced quicker in the stiffer tissue. The longitudinal strain of the left atria was halved when the stiffness of the tissue was increased.

Altered Mechanics Due to the Active Force Development of LV
We quantified isolated changes in the maximal active force developed in the tissue by comparing the deformation of the simulation with control active force (Case 6) to the one with decreased active force (Case 8), both with the hypertrophic geometry (HCM 2). At the ED, the regional wall thickening was comparable between Case 6 and Case 8. At ES, the regional wall thickening differed-the extent of the thickening values remained similar, while the maximum decreased when the force was decreased. At the same time, the distribution of the thickening values among the AHA segments was retained when the active force was decreased ( Figure 2). The regional circumferential and radial strain at the ES had similar extent of the values (≈13.5 %) for Case 6 and Case 8. The ranges at ES were shifted-the circumferential strain indicated less shorting of the tissue (values are higher, since they are negative) and the radial strain indicated less elongation of the tissue (values are lower, since they are positive) ( Figure 6). The maximal global circumferential and radial strain rates during the entire heart cycle differed between between Case 6 and Case 8. The circumferential strain rate indicated slower decrease of the strain, while the radial strain rate indicated slower increase of the strain when the active force was reduced ( Figure 8, on the left).
The velocities in all directions (longitudinal, circumferential and radial) were reduced during the entire heart cycle when the active force was reduced ( Figure 8, on the right). The longitudinal strain of the left atria was lower for the case of reduced active force-20 % for Case 6 to 15 % for Case 8. A comparison of the deformation of the simulation with control active force (Case 1) and the one with decreased active force (Case 3), both with the control geometry, confirmed these results. Summary: A reduced maximal active force development in the tissue of the LV results in reduced wall thickening as well as reduced circumferential and radial strain at the ES in all AHA segments. It also results in reduced strain rates and velocities in the entire heart cycle. In total, less deformation is available in the entire LV.

Altered Mechanics Due to Fiber Disarray of LV
We quantified isolated changes in the FO of the LV by comparing the deformation of the simulation with control FO (Case 6) to the one with disarrayed FO (Case 10). The RMSD of the global wall thickening was less than 3.2 % for the entire heart cycle (Case 6 vs. Case 10: 1.1 % (ED) and 3.2 % (ES)). Similarly, the regional wall thickening showed minor differences. At the ES, the range of the values for Case 6 was 18.1 %-51.3 % and for Case 10 it was 19.7 %-49.4 %, while differences occurred mainly in the basal segments (1-6). The regional longitudinal strain at the ES was not influenced by the disarrayed FO.
The regional circumferential and radial strain at the ES had smaller extend of the values for the case with disarrayed FO. The difference was more pronounced for the circumferential strain-the range of the values for Case 6 was -15.5 %--3.4 % and for Case 10, it was -12.9 %--2.7 %. Therefore, the circumferential strain indicated less shorting of the tissue (values are higher, since they are negative) and the radial strain indicated slightly less elongation of the tissue (values are lower, since they are positive) ( Figure 6). The regional circumferential and radial strain rate indicated a slower change in the strains at the beginning of the systolic period. The values of the longitudinal strain rate were similar-the systolic RMSD was 2.2 %. The velocities in all three directions were comparable-the RMSD was 0.003 m/s during the systole and 0.001 m/s during the diastole. The longitudinal strain of the left atria slightly decreased when the FO was disarrayed-from 20 % to around 18 %. A comparison the deformation of the simulation with control FO (Case 7) to the one with disarrayed FO (Case 11), both cases with increased stiffness, confirmed these results. Summary: Disarrayed FO in the mid-wall of the LV led to less deformation in circumferential direction-the strain indicated less circumferential shortening of the tissue. Furthermore, it led to slightly less deformation in radial direction. The strain rates and the velocities were not considerably changed.

Altered Mechanics Due to Combination of Pathological Model Components
We compared the deformation of the control case (Case 1) to the virtual HCM heart (Case 13), which was the combination of hypertrophic geometry, stiffened passive behaviour, decreased active force development and disarrayed FO. The RMSD of the global wall thickening was 7.6 % for the systole and 7.7 % for the diastole. The maximum of the regional wall thickening decreased-from 50.0 % for the control case to 27.7 % for virtual HCM heart. As previously described, a decrease in the maximum was observed when the active force was decreased (e.g. Case 6 vs. Case 8). Furthermore, the extent of the wall thickening values reduced since the maximum decreased more than the minimum (at ES, the wall thickening ranged from 18.7 % to 50.0 % for Case 1 and from 10.4 % to 27.7 % for Case 13).
As previously described, a decrease in the range of the wall thickening, which led to a more homogeneous distribution of the wall thickening, was observed when the stiffness of the tissue was increased (e.g. Case 1 vs. Case 2). Similar to the wall thickening, the regional strain values in all directions (longitudinal, circumferential and radial) were closer to zero for the HCM case compared to control case ( Figure 9). The extent of the strain values reduced as well-the ES longitudinal strain was between -11.0 % and 7.8 % for Case 1 and between -3.5 % and 5.5 % for Case 13. At ES, the circumferential strain for Case 1 ranged from -17.8 % to -4.5 % and for Case 13: from -8.0 % to -3.5 %. The radial strain for Case 1 ranged from 10.6 % to 25.7 % and for Case 13: from 4.8 % to 14.8 % ( Figure 6). The global strain rates reduced for the HCM heart compared to control case for the entire heart cycle ( Figure 10, left). For each direction, the RMSD of the strain rates were similar for the diastolic and systolic period-around 12 % in longitudinal direction, 23 % in circumferential direction and 41 % in radial direction. During the systolic period, the major difference was observed during the first half of the systole. The global velocity in all three directions differed as well during the systolic and diastolic period-the RMSD was between 9.0 % and 10.5 % during the diastole and between 11.3 % and 14.4 % during the systole (Figure 10, right  case compared to control case. In total, the deformation and its rate of the LV was diminished for the virtual HCM heart compared to the control case. Table 2 provides an overview of the results. Note, that the measures, which are negative, indicate less or slower deformation as they increase-e.g. an increase in the diastolic radial strain rate indicates slower relaxation.

Altered Mechanics in HCM Patients Reported in the Literature
In HCM patients, global strains and strain rates are reported to be significantly lower [1,11,12], while Ito et al. [9] measured preserved (or even increased) circumferential shortening. Furthermore, altered myocardial velocities are detected in HCM patients compared to controls-global and segmental diastolic velocities are decreased and systolic longitudinal velocities were reduced in HCM [10]. In the left atria, a higher minimum volume and a lower peak atrial longitudinal strain was measured in HCM compared to controls [11].

Causes of Altered Mechanics in HCM Patients Reported in the Literature
Previous studies have examined potential origins of altered cardiac mechanics in HCM patients in clinical studies. To the best of our knowledge, no computational study was conducted on this topic. In a cohort of fifty-nine HCM patients, Urbano-Moral et al. [2] demonstrated the relation of a reduction in longitudinal shortening of the LV and the extent of hypertrophy. Furthermore, the reduction of the global strain and strain rate was correlated with the mean WT [13]. In a clinical study, Villemain et al. [5] suggested that altered LV relaxation might result from increased myocardial stiffness. Hoskins et al. [6] hypothesized that reduced active force might contribute to systolic dysfunction in HCM patients.
In a HCM patients' heart, distinct underlying phenomena are present simultaneously, but are differently pronounced. Thus, the effects of these phenomena can not be clearly separated from measurements obtained in clinical studies. In contrast, in the presented numerical study, we could relate the observed alterations of mechanics to their underlying causes.

Identification of Underlying Pathophysiology
An alteration of the wall thickening and the strain time course during the systole can be related to hypertrophic LV walls. Longitudinal strain, which indicates less shortening is as well related to hypertrophic LV walls. Reduced ES wall thickening and strain values (circumferential and radial) in all segments can be related to less active force development in the LV. A homogeneous distribution of the ES wall thickening and strains among all AHA segments can be related to stiffer tissue. A reduction of the circumferential strain can be attributed to the fiber disarray in the mid-wall of the LV, but also to less active force development or increased stiffness of the tissue. Strain rates that are reduced during the entire systole and strongly reduced during the beginning of the systole can be traced back to increased tissue stiffness. Strain rates that are reduced during the entire heart cycle are caused by reduced active tension developed in the tissue. An increase in the myocardial velocities in all directions during the systolic period can be related to hypertrophic walls of the LV. In contrast, these velocities decreased during the systolic period when the stiffness of the tissue was increased.
Additionally, an increase in the velocities during the diastolic period combined with a rapid decay of the velocities can be as well related to an increased stiffness of LV tissue. A reduction of the active tension development did not affect the velocities during the diastolic period. A reduced longitudinal strain of the left atria at ES was present in case of increased stiffness of the tissue or reduced active tension development.

Comparison between Simulated and Clinically Measured Deformation
Wall thickening: A segmentation of the endocardial and epicardial surfaces of the LV for the entire heart cycle is time consuming and user dependent. Therefore, in clinical routine, the wall thickening is not commonly evaluated, but instead, the radial strain is used. On the virtual hearts, we applied two distinct methods to calculate the radial strain and the wall thickening. For both measures, we obtained comparable time course morphology and ES distribution for each case, but different maximal values (compare for Case 1 the wall thickening (max 50 %) in Figure 1, left and the radial strain (max 26 %) in Figure 9, left). Kato et al. [14] showed that the radial thickening of the tissue arises due to the circumferential fiber shortening. This implies that fiber disarray in the mid-wall (fiber oriented in circumferential direction) will change the radial thickening and strain. In contrast, we observed a reduction of the deformation in circumferential direction, which was more pronounced than the reduction of the deformation in radial direction. Strain and strain rate: The strain values obtained for the virtual control case (Case 1) indicated less shortening compared to healthy volunteers in all three directions [15]. Nevertheless, the regional values of the circumferential and radial strain of the virtual control heart are inside the ranges (mean±SD) provided in Table  4 in [15] but deviate from the mean values (e.g. the mean values of the regional circumferential strain were between -26 % and -17 % and the mean values of the regional radial strain were between 12 % and 39 %in [15]). Additionally, we obtained heterogeneous distributions among the AHA segments and heterogeneous values were as well measured-e.g. radial strain was 39±21 % in the anterior basal segment and it was 12±8 % in the septum basal segment [15]. The heterogeneity in the circumferential strain was significant (p <0.05) [16]. The regional longitudinal strain in the virtual control heart strongly deviated from the ranges provided for the healthy volunteers (from -24±11 % to -13±7 % [15]). For the virtual HCM heart, the longitudinal strain also deviated from literature values, which were around -11 % for HCM hearts [11]. For each virtual heart, the values in the free wall were negative as well, but closer to zero (Figure 6), while the positive values of the longitudinal strain were mainly in the septal segments. The torsion of the ventricle, which is pronounced in the apical region, leads to an elongation of the septal segments-the apex pulls the septal segments as it rotates. At the same time, the free wall shortened.
For the virtual HCM heart, the strain values indicated less (and slower) shortening in the circumferential direction and less (and slower) thickening in the radial direction of the tissue compared to the virtual control case, in agreement with [1,11,12]. In contrast to Ito et al. [9], the circumferential strain of the virtual HCM heart indicated clearly reduced shortening. The stiff heart syndrome (cardiac amyloidosis) can be identified by the longitudinal strain, which shows a strong gradient between the apex and base regions (apical sparing) [17,18]. In the stiffened virtual hearts (e.g. Case 2 and Case 7), we observed an opposite effect when the stiffness was increased, the strain became more homogeneous. Velocity: Li et al. [10] measured reduced global and segmental diastolic radial and longitudinal peak velocities in patients with HCM vs. controls. In contrast, we did not observe a reduction of any peak velocity during the diastolic period in our virtual HCM heart compared to the control heart. Nevertheless, we observed a decrease in these velocities when the active force was reduced (Case 3) and an increase when the stiffness was increased (Case 2). Therefore, we showed that the effect of these two pathologies cancel out in our virtual HCM heart to obtain the same peak diastolic velocity as in the control case ( Figure 10, right). If we would either further reduce the active force in the LV or less deviate from the control stiffness (or both), we would obtain reduced diastolic velocities in line with [10]. When the stiffness of the tissue was increased, the amplitude of the diastolic velocity increased and also the morphology of the diastolic velocity course changed compared to the control case. Therefore, we confirm the suggestion of Villemain et al. [5] that modified LV relaxation might result from increased myocardial stiffness. Furthermore, Li et al. [10] measured reduced peak velocities during the systole in patients with HCM in the longitudinal direction (radial peak velocities were comparable). We observed a reduction of the velocities during the systolic period as well but for all directions (longitudinal, circumferential and radial). We did not compare the absolute values of velocities, since different imaging modalities provide different values for the velocities-feature tracking peak velocities are lower than directly measured tissue phase mapping velocities [10]. Hoskins et al. [6] hypothesized that reduced active force might contribute to systolic dysfunction in HCM patients. In agreement with [6], we observed reduced systolic function when the force was decreased-the velocities during the systole were decreased and the strain at ES was diminished. In a cohort of healthy volunteers, similar global longitudinal and radial velocities during the systole were measured (longitudinal was 2.6±0.55 cm/s and radial was 2.5±0.36 cm/s) [15]. Likewise, we obtained similar global longitudinal and radial velocities. Nevertheless, the circumferential velocity was also in the same range, while Lin et al. [19] measured negative circumferential velocities since the direction of deformation along the circumference was considered. We calculated the absolute values of these velocities which equal the speed of deformation in the circumferential direction. However, the direction of the deformation can be deduced from the strain values.
Left atria (LA) strain: The LV of the virtual HCM heart deformed less compared to the control case, which led to reduced maximal longitudinal strain of the LA. Similar behaviour was reported by Aly et al. [11]. Additionally, they discovered that LA dysfunction is present in HCM patients before global LV dysfunction can be measured. We could not reproduce this behaviour in our model-in each case, in which LA longitudinal strain was reduced, the deformation of the LV was as well reduced. This reported statement demonstrates the importance of using a whole heart model, in which the deformation of all chambers is related.

Limitations
In the following, we describe limitations of our study and provide directions for potential improvements. For the control geometry, the finite element mesh of the ventricle was coarse-up to two elements in the transmural direction. The linear course of the FO is represented by two or three fibers in the transmural direction. Thus, we could not create a combination of control geometry and disarrayed FO. Furthermore, the linearity of the solution to Laplace's equation depends on the width of the domain between the boundaries [20]. Therefore, the course of the transmural coordinates used to define the LV mid-wall region close to the apex is not linear. In future, a linear course can be obtained by solving a trajectory distance equation in the LV [20]. The measure of FA was calculated on the fine geometry (also used to create the FO, Figure 12) and delivered values close to one (0.95±0.11) in the mid-wall ring of the control FO and 0.81±0.25 for the disarrayed FO. Ariga et al. [4] measured for a control case a FA of 0.52±0.03 and for HCM patients 0.49±0.05, which are considerably lower values compared to the values in our virtual hearts. This is the result of a rule based algorithm, since it creates idealized FO. Nevertheless, the difference between the mean FA of the control and the one of the disarrayed FO in the virtual heart is 0.14, which is much higher compared to this difference measured by Ariga et al. [4] (0.03). Therefore, we consider the cases with disarrayed FO to be representative for severely disarrayed FO. For all numerical simulations, an identical input parameter set for the circulatory model was applied. It led to healthy systolic pressure (120 mmHg) for the virtual control heart, but to lower systolic pressures in the virtual HCM heart (70 mmHg). The pressure, applied on the endocardial surface, influences the deformation and therefore, the evaluated measures. In both control and HCM cases, the ejection faction (EF) was lower compared to literature values-44 % vs. 72 % [21]. An increase of the active force will increase the EF but also the systolic pressure in the control case. A diminished deformation results in a reduced EF and therefore in reduced strains, as discussed previously. The active force model provided the force based on a predefined curve (force over time) and a maximum value of the force (T max ). A more complex model will adjust the intervals of the increase and the decrease of the force based on the velocities (e.g. Land et al. [22]) and therefore, the duration of the contraction and relaxation of the LV. Ito et al. [9] reported that regional LV filling for HCM hearts was prolonged compared to control hearts and that the impairment of the diastolic relaxation is a major sign of HCM. We observed the major changes in the measures during the systole rather than diastole with exception of the strain rates-the RSMD are higher during the systole for the strains, velocities and wall thickening Additionally, we did not evaluate the diastolic relaxation time extent. In general, the values of the metrics in the apex segment might be misleading, since the local longitudinal directions in the apex strongly deviate from the global longitudinal direction ( Figure 15). This is a result of the algorithm which creates the local directions. Instead we could use the global longitudinal direction as a local one, but then the three local directions (longitudinal, circumferential and radial) will not build an orthogonal system in each volume element and could not provide linearly independent information. We calculated of the strain based on the deformation tensor F . Werys et al. [23] used the Green-Lagrangian strain tensor E = 0.5(F T F − I) to derive the strain directly from the motion of the myocardium from cine MRI images. Santiago et al. [24] projected the components of the Green-Lagrangian strain tensor on a global longitudinal direction to obtain the strain.
There is no agreement in the literature how to apply the deformation tensor to obtain the strain and therefore, we would obtain different strain values depending on the definition of strain.

Conclusions
We conducted an in-silico study on virtual human whole hearts to identify causes of altered mechanics in hypertrophic cardiomyopathy (HCM) hearts. We simulated the deformation resulting from combinations of physiological and pathological models and evaluated the mechanical behaviour by local and global measures (wall thickening, strain, strain rate and deformation velocities). The presented study shows which pathological mechanisms are required to be present in the LV to obtain altered mechanics and how they affect the deformation measures. An increased wall thickness leads to deformation alteration during the systole, while the ES values are comparable to control case. Stiffer tissue equalizes the strains at ES, while reduced active force development reduces the deformation of LV. Disarrayed FO in the mid wall did not influence the deformation of the LV. An inversion of these arguments allows to identify present pathological mechanisms in the tissue which cause an altered mechanical behaviour.
In the clinical routine, it is cumbersome to directly measure underlying pathological mechanisms and therefore, those, derived from a numerical simulation, might be a valuable information for clinicians and can contribute to a more accurate diagnosis in HCM patients.

Methods
In the following, we introduce the physiological and pathological mechanisms which were included in the sensitivity analysis. Then, we describe the heart geometry and the numerical solver used to obtain the deformation of the entire heart for several heart beats. Finally, we present the metrics which measure alteration of left ventricular mechanics.

Sensitivity Analysis
We created geometries with increased wall thickness (WT) of the left ventricle, varied internal forces (passive and active) and distorted fiber orientation (FO) in the LV. For all cases, the external forces were defined using the same parameterization of circulatory and pericardial models.

Wall Thickness
The WT (mean ± std) of the LV of the initial geometrical model (described in The Geometrical Model of the Heart), was 10 ± 2.3 mm and its cavity volume was 193 ml. We added tissue on the endocardial surface to increase the thickness of the LV to 1) 15 ± 3.3 mm and 2) 17 ± 4.1 mm with a concomitant decrease of the LV volume (118 ml and 94 ml, respectively). In both cases, the added tissue was distributed concentrically. Figure 11 shows the three geometries with distinct WT. WT of the right ventricle and both atria were not modified compared to the initial geometrical model. Figure 11 Three LV geometries clipped through their long axis. Left: initial geometry (mean wall thickness 10 mm). Middle and Right: hypertrophic geometries (15 mm and 17 mm, respectively).

Passive Forces
We varied the input parameters of the passive force model (described in The Numerical Solver ) which determines the tissue stiffness. To identify the parameters of the passive force model for the control case, we used a method based on the pressure volume relation of LV as described previously [25] and obtained the following parameters for the Guccione model: C = 309 Pa, b f = 17.8, b t = 7.1 and b f t = 12.4. In HCM myocardium, Villemain et al. [5] measured a five fold increase of the stiffness compared to controls. Thus, we increased the parameter C which determines the global stiffness to 1545 Pa for the entire myocardium of all four chambers to capture increased tissue stiffness.

Active Forces
We varied the input parameters of the active force model (described in The Numerical Solver ). For the control case, the scaling parameter of the active force, T max , was set to 100 kPa in both ventricles and 35 kPa in both atria. These values were chosen to obtain a control systolic LV pressure of 120 mmHg in the control geometry. Hoskins et al. [6] measured a 40 % decrease of the active force in HCM donor cells compared to controls. Therefore, we reduced the maximal active force to T max = 60 kPa in both ventricles and to 21 kPa in both atria.

Fiber Orientation
We defined two configurations of FO in the LV: one control case and one representing fiber disarray. The control FO was determined by a rule-based algorithm based on Bayer et al. [26] with angles changing transmurally from 60 • on the endocardium to -60 • on the epicardium [27]. The algorithm (Bayer et al. [26]) was adapted to eliminate a discontinuity of fibers in the free walls and to yield a fiber rotation that is approximately linear across the wall [1] . Ariga et al. [4] measured fiber disarray in HCM hearts. We modified the rulebased algorithm to yield disarrayed FO in the mid-wall ring of the LV (described in Creation of Fiber Disarray). On the epicardium and endocardium, the same angles were used as in the control case. We quantified the disarray by calculating the fractional anisotropy (FA) in the mid-wall ring, which was 0.95±0.11 for control and 0.81±0.25 for the disarrayed FO ( Figure 12). The minimum FA (0.3) for the control fiber was observed at the junction of LV and RV.
Overview Table 1 provides an overview of the cases covered in the sensitivity analysis. Each case is defined by a combination of the four model variants described above.

The Geometrical Model of the Heart
The control geometry was based on MRI data of the whole heart, acquired from a healthy volunteer at University Hospital Heidelberg with a 1.5 T MR tomograph (Philips Medical Systems). Voxel spacing was 0.7×0.7×1.8 mm. The volunteer gave informed consent and the study was approved by the institutional review board. Images were segmented to obtain the endocardial and epicardial surfaces of the four chambers which provide the boundaries for the volume mesh of the myocardium. Additionally, the convex hull of the four chambers was calculated to serve as an inner surface of the pericardium. The volume between the myocardium and the  pericardium was defined as fat tissue. The veins and arteries connected to the myocardium were represented in the model as trunks. The entire volume mesh consisted of 48780 nodes and 90801 cells. We used quadratic tetrahedral elements to discretize the volume and linear triangular elements for the surfaces. The nodes on the free ends of all trunks and the outside surface of the pericardial sac were fixated in all three directions to serve as a boundary condition for the model ( Figure 13A ,B).

The Numerical Solver
To calculate cardiac deformation, we used the mechanical solver CardioMechanics [28], which was previously verified [29]. The governing equation ensures that all forces are in balance at all times during the heart beat. External forces arise outside the myocardium, internal forces arise inside the myocardium. To calculate the external forces, we included a circulatory model and a pericardial model. The circulatory model provides a pressure-volume relation in the four chambers and delivers the pressure values, which are acting on the endocardial surfaces [30]. The pericardial model represents the pericardial sac, in which the heart is embedded, and the surrounding tissue [28]. The pericardial model limits the motion of the heart by reducing the myocardial radial contraction and increasing the atrioventricular plane displacement. It delivers the forces acting on the epicardial surface of the entire heart. The internal forces are calculated by the combination of passive and active force models. The passive force model delivers the force arising from the intrinsic material properties of the myocardial tissue and is described by a constitutive relation. In this study, we applied the model of Guccione et al. [31] describing a hyperelastic, transversely isotropic material by the following strain energy function: where C, b f , b t and b f t are the parameters of the Guccione model, E ij (i, j ∈ [1, 2, 3]) are elements of the Green strain tensor, det(F ) is the determinant of the deformation tensor and K scales the incompressibility term. For the contractile tissue, K = 10 6 Pa was chosen and the parameters for the pericardium were chosen as in [28]. The fat tissue had the same passive properties as the pericardium. For the trunks, we applied the hyperelastic model of Mooney-Rivlin [32] with C 1 = 14900 Pa and C 2 = 0 Pa. The active force model delivers the force acting along the fiber direction, which leads to fiber shortening and, therefore, to the contraction of the tissue. In this study, active force was described by a predefined curve as described by Stergiopulos et al. [33]. The normalized curve was scaled by a parameter T max which determines the maximal active force. The ventricles were simultaneously activated 150 ms after the atria were activated (also simultaneously at 0 ms).

Modeling of Fiber Disarray
The myocardial cells tend to align along their long axis to form bundles that are represented by fibers in the geometrical model. The FO determines the deformation of the tissue [34]. In our mode, the FO in the atria was determined by the rule-based algorithm of Wachter et al. [35] ( Figure 13C ,D). Fiber directions were assigned for each of the four quadrature points and for the centroid of each element. Ariga et al. [4] visualized the myocardial microstructure of HCM hearts with diffusion tensor MRI (DT-MRI). DT-MRI allows quantifying the direction diffusion of water molecules by measuring the fractional anisotropy (FA). A diffusion-weighted signal intensity is measured to construct the diffusion tensor (DT) [36]. The DT is a 3×3 matrix obtained for each voxel and can be transformed to a diagonal matrix with its eigenvalues λ 1 , λ 2 and λ 3 as diagonal elements [36]. The eigenvector belonging to λ 1 , indicates the orientation of the long axis of the myocytes and λ 1 , the magnitude of the diffusion in this direction. The other two eigenvectors are orthogonal to the primary one and define a transverse orthogonal plane. FA is calculated from the eigenvalues of the DT as follows [36]: where D av is the mean diffusivity, D av = (λ 1 + λ 2 + λ 3 )/3. A FA value close to 0 corresponds to isotropic diffusion and therefore indicates tissue with variable FO. A FA value close to 1 corresponds to anisotropic diffusion and therefore indicates coherently aligned tissue [4]. Ariga et al. [4] measured reduced FA in the mid-wall ring (circumferentially aligned fibers) in the hearts of HCM patients compared to controls. We constructed a virtual DT to measure the FA in our computational heart model. In the geometrical model, the FO is known for each element. Therefore, we estimated the diffusivity of the fibers λ 1 in finite element regions to construct the virtual DT.
The diffusivity in the other two directions was set to λ 2,3 (v i ) = 0.5(1 − λ 1 (v i )). Finally, the values obtained for λ 1 , λ 2 and λ 3 were used in Equation 2 to obtain FA for the provided fiber configuration. Here again, for a coherent fiber arrangement in a specific region, we obtain λ 1 = 1, λ 2,3 = 0 and therefore FA = 1. In a region of strongly disarrayed FO, we obtain λ 1,2,3 = 1/3 and therefore, FA = 0. We adapted the fiber assignment algorithm to generate disarrayed FO in the midwall ring of the LV (Figure 12). The mid-wall ring was defined to enclose all elements with transmural coordinates between 0.34 and 0.66. The transmural coordinate ranged from 1 on the endocardial surface to 0 on the epicardial surface and was obtained by solving the Laplace's equation in the volume. In this ring, the gradient value in each element was multiplied by a random number from a uniform distribution on the interval [0, 1] to obtain the distorted FO in this element. The sheet and sheet-normal directions were calculated to yield an orthonormal system together with the distorted fiber. Outside the mid-wall ring, the FO was assigned as in the control case. The FO were generated on a fine mesh (around one million elements). A nearest neighbor interpolation transferred the FO to the coarse geometry used for the simulations (Figure 11, right).

Evaluation Metrics
We introduced metrics to evaluate the deformation of the LV and one metric for the left atrium (LA) based on common imaging-derived features [37,38].
The following measures were evaluated globally (one value for the entire ventricle per time point) and regionally (one value per one of the 17 AHA segments [39] per time point, Figure 14): strain, strain rate, velocity and wall thickening. Figure 14 The division of the LV in 17 AHA segments. Each segment is numbered in both the anterior and posterior views. On the left of the gray separation line is the control geometry and on the right, the hypertrophic geometry (HCM 2). The apex segment (17) includes the endocardial apex.
The strain, strain rate and velocities were calculated in a local heart coordinate system (R-Lo-C), spanned by radial, longitudinal and circumferential directions [37] ( Figure 15). For every finite element in each geometry, these axes were calculated at the initial time point of the simulation and preserved over the heart beat. Regional and global measures were derived as the mean over all elements in the respective regions. All measures were calculated during the systolic and during the dias- tolic period. The regional measures at end-systole and end-diastole are visualized in bull's-eye displays [39]. Time of end-systole and end-diastole was determined based on the pressure-volume relation. For each global measure, we calculated the RMSD (root mean squared deviation) between each pair-wise combination of cases (Table 1) during the systolic and diastolic period.

Strain
Strain ε (%) describes the change in length relative to the initial length (in one dimension). Positive strain values correspond to lengthening and negative to shortening [37].
In the numerical simulation, the deformation of the heart in each element of the mesh is characterized by a deformation gradient tensor F calculated for each time step of the heart beat. The FO f d in the deformed element is F f init , where f init is the initial FO. The stretch of the deformed fiber is λ(f d ) = f T init F T F f init and the strain is ε(f d ) = λ(f d ) − 1. The strain in the deformed sheet ε(s d ) and sheet-normal ε(sn d ) directions is calculated likewise. The strain in the R-Lo-C system is then obtain by a coordinate transformation: The matrix T transforms the Cartesian coordinates (with respect to the standard basis of R 3 ) into the local R-Lo-C coordinates. The rows of T are the radial, longitudinal and circumferential vectors in the current element. The strain in radial direction is the first entry of the transformed strain: ε RLoC (1), the strain in longitudinal direction is ε RLoC (2) and the strain in circumferential direction is ε RLoC (3).

Strain Rate
Strain rate (%/s) is the speed at which the strain changes [37]. We calculated the strain rateε t at time t with ∆t = 0.01 s: where ε t is the strain at time t.

Velocity
Velocity (m/s) is the temporal change of the displacement. By solving the governing equation of the heart mechanics, we obtain the displacement and thus, the velocity for every node and each time step [40].
To obtain the velocity of each finite element, the mean of the velocity over its four nodes was calculated for each of the three directions in the Cartesian coordinate system: v = (v x , v y , v z ). To convert the Cartesian velocity into the local R-Lo-C system, the coordinate transformation with the matrix T was conducted analogous to the transformation of the strain: v RLoC = T v T . Then, for the observed finite element at the current time point, the absolute value of the velocity in radial direction is |v RLoC (1)|, in longitudinal direction |v RLoC (3)| and in circumferential direction |v RLoC (3)|. In the following, the absolute value of the velocity will be referred to as velocity.

Wall Thickening
The wall thickness (WT, ω in mm) of the LV was calculated as proposed by Yezzi et al. [41]. WT was obtained for every time step of the heart beat and each surface node. The wall thickening (in percent) for time t is then: (ω t − ω 0 )/ω 0 , where the upper index corresponds to the time with 0 denoting the initial WT.

Mechanics Of Left Atrium
We defined a main longitudinal axis of the LA by calculating the mean of the local longitudinal directions over all elements in the LV. For every time point t, all nodes of the LA were orthogonally projected on the main axis. With the maximal euclidean distance between any two projected points l t LA , the longitudinal strain for time t is then: (l t LA − l 0 LA )/l 0 LA , where the upper index 0 correspond to the initial geometry configuration at time 0.

Declarations Ethics approval and consent to participate
The volunteer gave informed consent and the study was approved by the institutional review board.

Consent for publication
Not applicable