In the following, we describe the heart geometry and the numerical solver used to obtain the deformation of the entire heart for three heart beats. Then, we describe how we modeled fiber disarray as measured in HCM patients and introduce the physiological and pathological mechanisms included in the sensitivity analysis. Finally, we present the metrics which measure alteration of left ventricular mechanics.
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 48,780 nodes and 90,801 cells. In the LV, there were two elements transmurally, which suffices to obtain correct deformation according to the convergence analysis conducted by Gerach et al. [34]. 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 (Fig. 13A ,B).
The numerical solver
To calculate cardiac deformation, we used the mechanical solver CardioMechanics [35], which was previously verified [36]. To describe the cardiac mechanics, the equation of balance of the linear momentum is solved by the Finite-Element Method. 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 closed-loop 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 [34]. The closed-loop model ensures that the total blood volume in the circulatory system is preserved over several heart beats. The model is strongly coupled to the finite-element model as described by Gerach et al. [34]. The input parameters and the initial conditions are provided in Additional file 1: Tables S1 and S2). The pericardial model represents the pericardial sac, in which the heart is embedded, and the surrounding tissue [35]. A sliding boundary condition is imposed between the inner surface of the pericardial model and the outer surface of the heart model. 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. [37] describing a hyperelastic, transversely isotropic material by the following strain energy function:
$$\begin{aligned} W&= \frac{C}{2} \left( e^Q - 1\right) + \frac{K}{2} \left( \text {det}(F) - 1\right) ^2, \nonumber \\ Q&= b_f E^2_{11} + b_t \left( E^2_{22} + E^2_{33} + E^2_{23} + E^2_{32}\right) + b_{ft} \left( E^2_{12} + E^2_{21} + E^2_{13} + E^2_{31}\right) , \end{aligned}$$
(1)
where C, \(b_f\), \(b_t\), and \(b_{ft}\) are the parameters of the Guccione model, \(E_{ij}\) \((i,j \in \left[ 1,2,3\right] )\) are elements of the Green strain tensor, \(\text {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 [35]. The fat tissue had the same passive properties as the pericardium. For the trunks, we applied the hyperelastic model of Mooney–Rivlin [38] 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. [39]. The normalized curve was scaled by a parameter, 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 [40]. In our model, the FO in the atria was determined by the rule-based algorithm of Wachter et al. [41] (Fig. 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 DT-MRI. It allows quantifying the direction diffusion of water molecules by measuring the FA. A diffusion-weighted signal intensity is measured to construct the diffusion tensor (DT) [42]. The DT is a 3\(\times\)3 matrix obtained for each voxel and can be transformed to a diagonal matrix with its eigenvalues \(\lambda _1\), \(\lambda _2\), and \(\lambda _3\) as diagonal elements [42]. The eigenvector belonging to \(\lambda _1\) indicates the orientation of the long axis of the myocytes and \(\lambda _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 [42]:
$$\begin{aligned} FA = \sqrt{\frac{3}{2}}\sqrt{ \frac{(\lambda _1 - D_{av})^2 + (\lambda _2 - D_{av})^2 + (\lambda _3 - D_{av})^2}{\lambda _1 ^2 + \lambda _2 ^2 + \lambda _3 ^2} }, \end{aligned}$$
(2)
where \(D_{av}\) is the mean diffusivity; \(D_{av} = (\lambda _1 + \lambda _2 + \lambda _3) / 3\). An FA value close to 0 corresponds to isotropic diffusion and therefore indicates tissue with variable FO. An 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 \(\lambda _1\) in finite-element regions to construct the virtual DT. We subdivided the LV in N regions (\(v_i\), \(i = 1, \dots , N\), N = 1500) of similar size (around \(6 \times 20 \times 20\) mm) and calculated in each region the mean FO \(f^{mean}_i\). For each element in the region (\(e^k_i\), \(k = 1, \dots , M\) with M the number of elements in the current region), we calculated the length of the projection of the fiber on the mean FO (\(l^k_i\)). Then, we set \(\lambda _1\) of the region \(v_i\) to the mean of these lengths across all elements in the region
$$\begin{aligned} \lambda _1 (v_i) = \frac{1}{M}\sum _{k=1} ^M l^k_i. \end{aligned}$$
(3)
The diffusivity in the other two directions was set to \(\lambda _{2,3} (v_i) = 0.5( 1 - \lambda _1 (v_i))\). Finally, the values obtained for \(\lambda _1\), \(\lambda _2\), and \(\lambda _3\) were used in Eq. 2 to obtain FA for the provided fiber configuration. Here again, for a coherent fiber arrangement in a specific region, we obtain \(\lambda _1 = 1\), \(\lambda _{2,3} = 0\) and therefore FA \(= 1\). In a region of strongly disarrayed FO, we obtain \(\lambda _{1,2,3} = 1/3\), and therefore, FA \(=0\).
We adapted the fiber assignment algorithm to generate disarrayed FO in the mid-wall ring of the LV (Fig. 14). 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 1 million elements). A nearest-neighbor interpolation transferred the FO to the coarse geometry used for the simulations (Fig. 15, right).
Sensitivity analysis
We created geometries with increased WT of the LV, varied internal forces (passive and active), and distorted 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\,\pm \,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\,\pm \,3.3\) mm and 2) \(17\,\pm \,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 15 shows the three geometries with distinct WT. WT of the right ventricle and both atria were not modified compared to the initial geometrical model.
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 [43] and obtained the following parameters for the Guccione model: \(C =\) 309 Pa, \(b_f =\) 17.8, \(b_t =\) 7.1, and \(b_{ft} =\) 12.4. In HCM myocardium, Villemain et al. [5] measured a fivefold 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^{\text {V}}_{\text {max}}\), was set to 100 kPa in both ventricles and \(T^{\text {A}}_{\text {max}}\) = 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^{\text {V}}_{\text {max}} = 60\) kPa in both ventricles and \(T^{\text {A}}_{\text {max}} = 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. [44] with angles changing transmurally from 60\(^\circ\) on the endocardium to \(-60^\circ\) on the epicardium [45]. The algorithm (Bayer et al. [44]) 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.Footnote 1
Ariga et al. [4] measured fiber disarray in HCM hearts. We modified the rule-based algorithm to yield disarrayed FO in the mid-wall ring of the LV (described in “Modeling 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 FA in the mid-wall ring, which was 0.95 ± 0.11 for control and 0.81 ± 0.25 for the disarrayed FO (Fig. 14). The minimum FA (0.3) for the control fiber was observed at the junction of LV and RV.
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 [46, 47].
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 [48] per time point, Fig. 16): strain, strain rate, velocity, and wall thickening.
The strain, strain rate, and velocities were calculated in a local heart coordinate system (R-Lo-C), spanned by radial, longitudinal, and circumferential directions [46] (Fig. 17). 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 diastolic period. The regional measures at end-systole and end-diastole are visualized in bull’s-eye displays [48]. 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 \(\varepsilon\) (%) describes the change in length relative to the initial length (in one dimension). Positive strain values correspond to lengthening and negative to shortening [46].
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 \(\lambda (f_{d}) = \sqrt{f^T_{init} F^T F f_{init}}\) and the strain is \(\varepsilon (f_{d}) = \lambda (f_{d}) -1\). The strain in the deformed sheet \(\varepsilon (s_{d})\) and sheet-normal \(\varepsilon (sn_{d})\) directions is calculated likewise.
The strain in the R-Lo-C system is then obtain by a coordinate transformation with matrix T. The matrix T transforms a vector with Cartesian coordinates (with respect to the standard basis of \(\mathbb {R}^3\)) into the local R-Lo-C coordinates. The rows of T are the radial, longitudinal, and circumferential vectors in the current element. By multiplication with the matrix T from the left, the normalized vectors pointing in the fiber (\(f_{d}\)), sheet (\(s_{d}\)) and sheet-normal (\(sn_{d}\)) directions in a deformed element are projected on each local direction vector, pointing in longitudinal, circumferential, and radial directions. Then, the projections are scaled by the strain obtained in the deformed element \(\lambda (f_{d})\), \(\lambda (s_{d})\), and \(\lambda (sn_{d})\). The following equation provides the strain in the R-Lo-C system:
$$\begin{aligned} \varepsilon _{RLoC} = \left| T(f_{d}/\Vert f_{d}\Vert _2) \right| \varepsilon (f_{d}) + \left| T(s_{d}/\Vert s_{d}\Vert _2) \right| \varepsilon (s_{d}) + \left| T(sn_{d}/\Vert sn_{d}\Vert _2) \right| \varepsilon (sn_{d}). \end{aligned}$$
(4)
The strain in radial direction is the first entry of the transformed strain: \(\varepsilon _{RLoC}(1)\), the strain in longitudinal direction is \(\varepsilon _{RLoC}(2)\), and the strain in circumferential direction is \(\varepsilon _{RLoC}(3)\).
Strain rate
Strain rate (%/s) is the speed at which the strain changes [46]. We calculated the strain rate \(\dot{\varepsilon }_t\) at time t with \(\Delta t =\) 0.01 s
$$\begin{aligned} \dot{\varepsilon }_t = (\varepsilon _t - \varepsilon _{t - \Delta t}) / \Delta t, \end{aligned}$$
(5)
where \(\varepsilon _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 [49].
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} = Tv^T\). Then, for the observed finite element at the current time point, the absolute value of the velocity in radial direction is \(\left| v_{RLoC}(1)\right|\), in longitudinal direction \(\left| v_{RLoC}(3)\right|\), and in circumferential direction \(\left| v_{RLoC}(3)\right|\). In the following, the absolute value of the velocity will be referred to as velocity.
Wall thickening
The wall thickness (WT, \(\omega\) in mm) of the LV was calculated as proposed by Yezzi et al. [50]. 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: \((\omega ^t - \omega ^{0})/\omega ^{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 corresponds to the initial geometry configuration at time 0.
