Effect of myofibre architecture on ventricular pump function by using a neonatal porcine heart model: from DT-MRI to rule-based methods

1. Introduction

Cardiac disease remains the leading cause of mortality and morbidity worldwide, as a result, extensive research has been carried out to develop computational cardiac models to understand mechanical behaviours of the heart [1–3]. For instance, finite-element (FE) method has been widely used to model heart function physiologically or pathologically, and to develop novel therapies [2,4,5]. The remaining challenges are to deal with the complex geometry, myofibre structure and material characterization of the myocardium [6,7]. Recent reviews on heart modelling can be found in [3,8,9].

The spatial architecture of myofibres plays a central role in electrical propagation, myocardial expansion and contraction [10]. Early studies relied on fibre dissections and histological slices [11] to determine local fibre structure. Currently, the cardiac fibres can be imaged via diffusion tensor magnetic resonance imaging (DT-MRI) [12] that allows a direct description of the three-dimensional myofibre architecture. To reconstruct myofibres in computational models, two different approaches have been developed. One is directly mapping myofibres from ex/in vivo datasets to the models, i.e. reconstructing models directly from DT-MRI [6], or using atlas-based methods to warp DT-MRI data into different models [13]. The other approach is the rule-based method (RBM), in which myofibres rotate from endocardium to epicardium with prescribed angles concerning the circumferential direction, varied linearly across the wall in most of the studies [14–16]. One key step in RBM is to parametrize the normalized wall thickness ( $\bar{e}$) in order to assign local fibre angles, from $\bar{e}=0$ at endocardial surface to $\bar{e}=1$ at epicardial surface. With measured fibre angles at endocardium θ^endo and epicardium θ^epi, the local fibre angle can then be assigned by varying linearly or nonlinearly with $\bar{e}$. Bayer et al. [17] proposed a Laplace–Dirichlet RBM, in which the circumferential-radial (transmural)-longitudinal directions and normalized wall thickness are determined by solving a series of Laplace equations. They demonstrated that the Laplace-Dirichlet rule-based fibre could achieve almost identical electrical activation patterns in a whole heart model as a DT-MRI-based model. Additionally, regionally varied RBM has been developed to take into account spatial variations of myofibre rotation angles [18].

Three-dimensional FE mechanics models of the heart have been used extensively to investigate the role of myofibre architecture in cardiac function under normal and abnormal function, including ischaemia, ventricular pacing, myofibre disarray and heart failure. For example, By using a rule-based approach for myofibre reconstruction in an left ventricle (LV) model, Wang et al. [15] found that changes in myofibre rotation angle can dramatically affect the stress and strain distributions during diastole. Using a bi-ventricular model, Palit et al. [19] also demonstrated that changes in myofibre angle can significantly affect myofibre stress–strain distribution within the LV wall in diastole. Pluijmert et al. [20] found that a change of 8° in myofibre orientation along transmural direction can cause a considerable increase in cardiac pump work (17%). In a recent study, Gil et al. [21] compared three different myofibre architectures in an electromechanics bi-ventricular model, one is from a DT-MRI dataset [22], the other two are reconstructed using a rule-based approach [14] with histologically measured myofibre angles [23]. Their results showed that the model with realistic myofibre structure from DT-MRI produces functional scores much closer to healthy ranges than rule-based approaches. By using the polynomial chaos expansion method, Rodríguez-Cantano et al. [24] studied the uncertainty in myofibre orientation and demonstrated that a realistic myofibre structure is necessary for a personalized cardiac model, such as DT-MRI-acquired myofibres.

Furthermore, myofibrils do not align perfectly along one direction at any location within a ventricular wall, but dispersed as reported by Ahmad et al. [25], who measured in-plane and out-of-plane myofibre and collagen fibre dispersion using two-photon-excited fluorescence and second harmonic generation microscopy on neonatal heart samples. To incorporate fibre dispersion in material constitutive law, Gasser et al. [26] introduced a structural tensor to account for collagen fibre dispersion in arterial tissue, they assumed a rotational symmetry for fibre distribution and a compact form of characterizing fibre dispersion was then given with one dispersion parameter, the so-called κ-model. In a series of studies [27,28], Holzapfel and co-workers used this generalized structural tensor to characterize the passive response of fibre-reinforced soft tissues. Later on, Pandolfi and co-workers [29–31] extended Gasser’s general structural tensor approach by including the second-order term of the Taylor expansion on the mean invariant along the fibre direction, to improve the accuracy of a structural tensor with large dispersions. In a recent study, Melnik [32] further extended the generalized structural tensor to include fibre dispersion in a coupled strain invariant.

Although there are several studies on passive constitutive responses of soft tissue [26,27], very few studies included fibre dispersion in active contraction models for the myocardium. There are two commonly used approaches for modelling active contraction in biological tissue: the active stress formulation [5,6,18,33] and the active strain formulation [9,34,35]. In the active stress formulation, the total stress tensor is decomposed into passive and active parts [36]. This approach has been widely used in personalized cardiac modelling because of its easy implementation, and the fact that there are abundant experimental data for the parameter calibration [5,6,18,33]. In the active strain approach, the total deformation gradient F is multiplicatively decomposed into an elastic part (F^pass) for passive response and an activation part (F^act), which could be more inherent to the ‘sliding filament theory’ [34]. The same structural tensor for the passive response could be linked to F^act to account for the active response [35]. This seems to be an elegant approach, though fitting personalized parameters to experimental data remains a challenge [34]. It is for this reason, that the active stress approach is still adopted here. To take into account active contraction caused by dispersed myofibres, Guccione and co-workers introduced cross-fibre active contraction in cardiac models [33,37] based on experiments by Lin & Yin [38]. Recently, Sack et al. [6] inversely determined cross-fibre contraction ratio in a healthy porcine heart and a failing heart. It has been argued that cross-fibre active contraction may be related to myofibre dispersion. However, no detailed studies have reported this connection. Eriksson et al. [18] incorporated myofibre dispersion in both the passive and active mechanics in an electromechanically coupled idealized left ventricular model. Their model, based on the κ-model [26], showed that large dispersion in the diseased heart could greatly affect the ventricular pump function. On the other hand, Ahmad’s study [25] demonstrated that in-plane dispersion is different from out-of-plane dispersion, which suggests the rotational symmetry assumption used in the κ-model may not be appropriate. Therefore, for the active stress formulation, a better approach would be to use the non-symmetric dispersion model developed in [27]. In this study, myocardial contraction is modelled following the active stress approach similar to that in [6].

Overall, there is a lack of studies on how different myofibre generation approaches, DT-MRI derived or RBM, affecting ventricular pump functions. One particular question is whether the difference between DT-MRI- and RBM-based models can be rectified using a proper consideration of fibre dispersion. We hypothesize that incorporating a non-symmetrical dispersed active tension model in an RBM-generated myofibre architecture can approximate the DT-MRI-based approach when simulating the heart pump function.

2. Material and methods

2.1. Geometry and fibre construction

A three-dimensional FE bi-ventricular model from [39] is used in this study (figure 1a), which is reconstructed from a computed tomography (CT) data of a neonatal porcine heart. Details of the data acquisition can be found in [40]. The three-dimensional CT data are first segmented using Seg3D,¹ then the boundary contours are exported into SolidWorks (Dassault Systemes, MA, USA) for geometry reconstruction, and meshed (figure 1a) using ICEM (ANSYS, Inc., PA, USA).

• Download figure
• Open in new tab
• Download PowerPoint

Because the myofibre structure of the neonatal porcine heart is not available, it is interpolated from a canine heart obtained from the public dataset of Cardiovascular Research Grid² [22]. We first reconstruct a bi-ventricular geometry for the canine heart with myofibres extracted from the primary eigenvector of the DT-MRI tensors, as shown in figure 1b. Clearly, the neonatal bi-ventricle geometry is different from the canine geometry, as shown in figure 1. Therefore, we cannot directly interpolate the measured canine myofibre structure for the neonatal bi-ventricle model. Instead, Deformetrica³ is employed to register the two bi-ventricular geometries by warping a template ( $C_{\alpha }$: the canine bi-ventricle) to a target ( $C_{\beta }$: the neonatal porcine heart) by minimizing a loss function that measures the distance between the template and target. Deformetrica is an open-source package based on a large deformation diffeomorphic metric mapping (LDDMM) framework [41,42]; further details about Deformetrica are given in the electronic supplementary material.

After warping $C_{\alpha }$ into $C_{\beta }$, the displacement fields u for all nodes on the external surface of $C_{\alpha }$ are obtained, denoting ${\mathbf{u}}_{\mathrm{Ex}}^{\mathrm{LDDMM}}$ as shown in figure 1c. The displacement vectors on the nodes lying within the ventricular wall are then interpolated by solving a Laplace system with Dirichlet boundary conditions (equation (2.1)) in Fenics,⁴

Following the finite deformation theory, the deformation gradient of warping the canine bi-ventricle model into the porcine model is

in which I is the identity matrix. Note that F and u are associated with the canine bi-ventricle model. Myofibre orientation in the warped canine model is where ${\mathbf{f}}_{\mathrm{template}}^{\mathrm{canine}}$ is the unit myofibre direction from the DT-MRI canine dataset. Finally myofibres in the porcine model f₀ are assigned according to the nearest neighbours between the warped canine and porcine geometries, such that in which x^procine is a position vector in the porcine model, and ${\mathbf{x}}_{\mathrm{warp}}^{\mathrm{canine}}$ is the position vector in the warped canine model. The sheet direction s₀ is defined transmurally across the wall, and the sheet-normal is n₀ = f₀ × s₀.

We further generate two different myofibre structures in the left side of the bi-ventricle using a rule-based approach [15], septum included. By projecting f₀ into the c₀ − l₀ plane to have ${\mathbf{f}}_0^{\parallel }$, we define the myofibre angle as the angle between ${\mathbf{f}}_0^{\parallel }$ and c₀, as shown in figure 2a. The average myofibre angles in the porcine model are then summarized at endocardium ( ${\theta }_{\mathrm{endo}}^{\hspace{0.17em}\mathrm{ave}}$) and epicardium ( ${\theta }_{\mathrm{epi}}^{\mathrm{ave}}$) in two ways: (i) across the whole LV, and (ii) at each ventricular segment according to the AHA17 (American Heart Association) definition [43] as shown in figure 2b,c based on right ventricular insertion points. A rule-based approach is used to generate two different myofibre structures: (i) one set of myofibre rotation angles varies linearly from endocardium to epicardium for the whole LV; (ii) for each AHA17 segment, myofibre rotates linearly based on the average rotation angles from that segment, which means myofibre angles are different at different segments. Note that the myofibre structure in the right ventricle (RV) of the bi-ventricle model, excluding the septum, is generated by the same rule-based approach but using one set of rotation angles due to the lack of DT-MRI data for the right side. We further assume the myofibre rotation angles at RV are the same as the angles when averaged across the whole LV.

• Download figure
• Open in new tab
• Download PowerPoint

With these three myofibre structures generated (figure 2d–f), we consider the following cases:

—	case LDDMM: the LV with the mapped ex vivo DT-MRI acquired myofibre architecture (equation (2.4));
—	case RBM17: myofibre rotates linearly from endocardium to epicardium for each LV segment according to the average rotation angles at 17 segments, derived from case LDDMM; table 1 lists the myofibre rotation angles at each segment, including the angles for the RV;
—	case RBMuni: myofibre uniformly rotates between endocardium and epicardium in the whole LV with one set of the average rotation angles (endocardium: 40°, epicardium: −30°), which are also derived from case LDDMM.

Note that case RBM¹⁷ has a heterogeneous myofibre structure in the whole LV but homogeneous within each segment. We also have not smoothed rotation angles between segments since those variations are within the range of local angle variations in case LDDMM as suggested in figure 4a. Thus, case RBM¹⁷ is a simplification of case LDDMM. Case RBM^uni has the same myofibre structure across the whole LV, a further simplification compared with case RBM¹⁷.

2.2. Constitutive model

2.2.1. Passive stress response

The passive behaviour of myocardium is described by a strain-invariant-based function [39], which is reduced from the model proposed by Holzapfel & Ogden [7] by fitting to an experimental study of neonatal porcine myocardium [40]. The strain energy function consists of a deviatoric (Ψ_dev) and a volumetric (Ψ_vol) parts,

where a, b, a_i, b_i, a_ij, b_ij are material constants and D is a multiple of the bulk modulus K, i.e. D = 2/K. J = det(F), $\mathbf{F}=J^{1/3}\bar{\mathbf{F}}$ and $\bar{\mathbf{C}}={\bar{\mathbf{F}}}^{\mathrm{T}}\bar{\mathbf{F}}$. The isochoric invariants are defined as ${\bar{I}}_1=\text{trace}(\bar{\mathbf{C}})$, ${\bar{I}}_{4\mathrm{f}}={\mathbf{f}}_0\cdot \bar{\mathbf{C}}{\mathbf{f}}_0$, ${\bar{I}}_{4\mathrm{n}}={\mathbf{n}}_0\cdot \bar{\mathbf{C}}{\mathbf{n}}_0$, ${\bar{I}}_{8\mathrm{fs}}={\mathbf{f}}_0\cdot \bar{\mathbf{C}}{\mathbf{s}}_0$ and ${\bar{I}}_{8\mathrm{fn}}={\mathbf{f}}_0\cdot \bar{\mathbf{C}}{\mathbf{n}}_0$, in which f₀, s₀, n₀ are the myofibre, sheet and sheet-normal directions in the reference state. In this study, we assume the collagen fibres follow the layered myocyte structure. Thus, myofibres represent both myocyte and collagen fibres. The max () in equation (2.5) will ensure the collagen fibres can only bear load when in tension. The passive Cauchy stress tensor is given by in which ${\bar{\psi }}_i=\mathrm{\partial }{\mathit{\Psi }}_{\mathrm{dev}}/\mathrm{\partial }{\bar{I}}_i,\hspace{0.17em}i\in \{1,\text{4f},\text{4n},\text{8fs},\text{8fn}\}$, $\bar{\mathbf{f}}=\bar{\mathbf{F}}\hspace{0.17em}{\mathbf{f}}_0$, $\bar{\mathbf{s}}=\bar{\mathbf{F}}\hspace{0.17em}{\mathbf{s}}_0$, $\bar{\mathbf{n}}=\bar{\mathbf{F}}\hspace{0.17em}{\mathbf{n}}_0$, $\bar{\mathbf{b}}=\bar{\mathbf{F}}{\bar{\mathbf{F}}}^{\mathrm{T}}$, p_vol = ∂Ψ_vol/∂J, and $\text{dev}(\bullet )=(\bullet )-(1/3)[(\bullet ):\mathbf{I}]\hspace{0.17em}\mathbf{I}$ denotes the deviatoric operator.

2.2.2. Active stress

Biaxial investigations on actively contracting rabbit myocardium [38] suggest that a large portion of active stress exists in cross-fibre direction. This has motivated computational efforts to include a proportion of the active stress to the cross-fibre direction when RBM generated myofibres are used [4,37]. In this study, we employ the active stress approach for myocardial active stress along the myofibre, sheet and sheet-normal directions

in which $\hat{\mathbf{f}}=\mathbf{f}/|\mathbf{f}|$, $\hat{\mathbf{s}}=\mathbf{s}/|\mathbf{s}|$ and $\hat{\mathbf{n}}=\mathbf{n}/|\mathbf{n}|$, n_f, n_s and n_n (all positive and sum up to 1) are the proportions of the active tension in their respective directions. T_a is the active tension generated along the myofibre direction, which is described by a time-varying elastance model that has been described extensively in the literature [6,36,37] where T_max is the maximum allowable active tension, Ca₀ is the peak intracellular calcium concentration, ECa₅₀ represents length-dependent calcium sensitivity, t is time and l is myofibre stretch. Further details are provided in the electronic supplementary material.

We assume the cross-fibre contraction in the RV is zero, i.e. n_f = 1, n_s = 0, and n_n = 0. This is because RV has a much thinner wall thickness, and Ahmad et al. [25] reported the fibre dispersion in the RV is much less than in the LV (9.3° versus 19.2°). We also performed simulations for the RBM^uni case, using the LV’s non-zero cross-fibre contraction for the RV. Our results show the differences of ejection fraction are 0.7% and 4.1% for the LV and RV, respectively. Thus assuming no cross-fibre contraction for the RV seems to be reasonable.

As for the LV (septum is included), since DT-MRI-derived myofibres are naturally dispersed in case LDDMM (figure 2), we set n_f = 1, n_s = 0 and n_n = 0. But for the RBM cases, it is necessary to include cross-fibre active tension, and n_s and n_n are calculated based on a dispersed fibre structure tensor. This will be explained in the following section.

2.2.3. Determination of n_f, n_n and n_s using DT-MRI derived myofibres for case RBM^uni

We first introduce {e₁, e₂, e₃} to denote the axes of a Cartesian coordinate system as shown in figure 3, and then define myofibre direction of the reference configuration to be M with a density distribution ρ(M). M can be further characterized by two angles Θ ∈ [0, π] and Φ ∈ [0, 2 π] (figure 3), that is

Θ is the angle between e₃ and M, and Φ is the angle between e₁ and the projected vector of M in the e₁–e₂ plane.

• Download figure
• Open in new tab
• Download PowerPoint

We assume the dispersions in different planes are essentially independent [44], i.e.

in which ρ_op(Φ) describes the out-of-plane dispersion, and ρ_in(Θ) describes the in-plane dispersion. Note in the ventricular model, in-plane is the plane defined by c₀ − l₀, and out-of-plane is the plane defined by l₀ − r₀. This is consistent with experimental studies when measuring in-/out-of-plane fibre angles [25,45]. The normalization of ρ(Θ, Φ) over a unit sphere requires in which N is a normalization factor.

When there is no dispersion, the structure tensor M ⊗ M can be directly used for constructing I_4f = C : M ⊗ M and active stress tensor T_a M ⊗ M/I_4f. With dispersion, a generalized structure tensor H can be defined over an unit sphere [18,26,27],

π-periodic von Mises distribution is then used for ρ(Θ) and ρ(Φ) [27],

in which θ is a variable representing Θ or Φ, b > 0 is the concentration parameter, $\frac{1}{\pi }{\int }_0^{\pi }\exp (b\cos (x))dx$ is the modified Bessel function of the first kind of order zero.

From figure 4a,b, we can find that in-plane angle (Θ) varies linearly from endocardium to epicardium for both RBM cases, but the fibres are much dispersed for case LDDMM, especially near the endocardium and epicardium, where myofibres align more longitudinally (l₀). The out-of-plane angle (Φ) is zero for both RBM cases since RBM generated myofibres only lie in the c₀ − l₀ plane. However, out-of-plane dispersion can be seen in case LDDMM shown in figure 4b. We now determine the in/out-of-plane dispersions from the angle differences between case LDDMM and RBM^uni. Figure 4c,d shows the histograms of in/out-of-plane dispersion in the LV, both Θ and Φ centre around 0°. The maximum-likelihood method mle() from Matlab is used to fit ρ_ip and ρ_op to the histograms of the in/out-of-plane dispersions, with b₁ = 1.6153 for the in-plane dispersion, and b₂ = 1.2144 for the out-of-plane dispersion.

• Download figure
• Open in new tab
• Download PowerPoint

Without loss of generality, we consider the mean fibre direction along e₃, the sheet direction along e₁ and the sheet-normal direction along e₂. Then the in-plane distribution is ρ_ip(Θ − 0, b₁), the out-of-plane distribution is ρ_op(Φ − (π/2), b₂), and

Similar to [18], we assume the active Cauchy stress with dispersed myofibres is Thus we have n_s = H₁₁ = 0.086, n_n = H₂₂ = 0.268 and n_f = H₃₃ = 0.646 for case RBM^uni.

2.3. Boundary conditions and implementations

The bi-ventricular model is implemented using the nonlinear FE software Abaqus (Dassault Systemes, Johnston, RI, USA). In order to simulate diastolic filling and systolic ejection, a lumped model for the pulmonary and systemic circulation systems is attached to this bi-ventricular model, which is realized through a combination of surface-based fluid cavities and fluid exchanges [46] as shown in figure 5. We define the mass flow rate between two different cavities as

where ρ is the blood density, A is the effective area between the two connected cavities and $\stackrel{\dot{}}{\bar{V}}$ is the fluid flux. $\stackrel{\dot{}}{m}$ is further related to the pressure difference where Δp is the pressure difference between two connected cavities, C_V is viscous resistance coefficient, and C_H is hydrodynamic resistance coefficient, and C_H = 0 in this study. This type of boundary conditions is equivalent to a simplified two-element Windkessel model. Parameters for the lumped circulation system are listed in table 2 and scaled from [6] by taking into account the dimensions of the neonatal porcine heart. For example, the total blood volume is around 80 ml for a newborn piglet [47], much less than in an adult porcine (67.2 ± 4.12 ml kg⁻¹) [48], the valvular area in a newborn heart is about one-tenth of the area in an adult heart [49–51], and the diameter of blood vessel is also much smaller in the newborn piglet compared with an adult porcine [52], which suggests that under similar pressure loadings, the vessel compliance, calculated as ΔV/ΔP will be much less in a newborn porcine because of much smaller ΔV in a newborn piglet.

• Download figure
• Open in new tab
• Download PowerPoint

Parameters for passive strain energy function and the maximum active tension from myocytes (T_max) are listed in table 3. Initial values for passive response are from [39], a_i (i ∈ {1, 4f, 4n, 8fs, 8fn}) are further reduced by half together with chosen T_max to ensure both LV and RV can achieve ejection fraction (EF) within the physiological range ( $\text{EF}>50\text{\%}$). Values of b_i (i ∈ {1, 4f, 4n, 8fs, 8fn}) are kept the same as in [39]. Note that because of missing measured data (wall motion, ventricular pressure) for the porcine heart, rather than constructing a personalized model [5,6], we only aim to obtain a set of parameters with which the bi-ventricle behaves physiologically.

The FE nodes on the top basal plane are constrained along the longitudinal axis but free to move within the basal plane. The longitudinal axis is defined as the line passing the LV basal centre and perpendicular to the basal plane. To start the simulation, linearly increased blood pressures from 0 to end-diastolic values are first applied to the inner surfaces of the bi-ventricular model, 8 mmHg in the LV and 4 mmHg in the RV. Typical diastolic pressures inside the pulmonary, left atrium, aorta and right atrium are also applied to those four cavities (10, 8, 67.5 and 4 mmHg [53]). Then the bi-ventricular model starts iso-volumetric contraction (t = 0 s), followed by systolic ejection when the ventricular pressure is higher than that of the aorta (around t = 0.045 s), and then the iso-volumetric relaxation. Systolic ejection ends at 0.12 s; 1 s is chosen for a whole cardiac cycle for computational convenience. In order to ensure the end-diastolic pressures in both LV and RV are same at next cardiac cycles, end-diastolic pressures in both atria are maintained constant.

3. Results

We first compare the heart pump function for cases LDDMM, RBM¹⁷ and RBM^uni without cross-fibre active tension. We then analyse the effect of cross-fibre active tension in case RBM^uni. Finally, we include dispersed active tension derived from DT-MRI myofibres in case RBM^uni and compared with case LDDMM.

3.1. No cross-fibre active tension

Figure 6a shows the pressure–volume loops from the three cases with no cross-fibre active tension. Although they all have the same end-diastolic pressure, the LV end-diastolic volume from case LDDMM (2.87 ml) is slightly larger than the other two rule-based cases (2.83 ml), the relative difference is around 1.4%. The LV end-systolic volume in case LDDMM is also the smallest (1.38 ml). Interestingly, though myofibre structures in the RV for the three cases are same, however, due to the difference in LV dynamics, the RV end-systolic volume from case LDDMM is also the smallest (0.87 ml).

• Download figure
• Open in new tab
• Download PowerPoint

Figure 6b shows ejection fractions for the three cases. Again, case LDDMM achieves higher ejection fraction both at LV (51.92%) and RV (55.47%) than the two rule-based cases. Furthermore, the LV ejection fractions for cases RBM¹⁷ and RBM^uni are less than 50%, which are below literature reported normal range (50%–75%), indicating the LV pump function is suboptimal in those two cases.

Figure 6c shows the average end-systolic stress for the entire LV along the circumferential, radial and longitudinal directions, respectively. Although the circumferential stress from case LDDMM is lower near endocardium and epicardium than RBM cases, it is much higher in the midwall, with the lowest value from case RBM^uni. Contrary to the circumferential stress, the longitudinal stress is higher in case LDDMM at endocardium and epicardium, while lowest at part of the midwall. The opposite trends of the circumferential and longitudinal stress levels in case LDDMM may compensate each other to achieve a deeper systolic contraction than cases RBM¹⁷ and RBM^uni. The radial stress is negative for all three cases with the lowest in case LDDMM.

Figure 6d is the apex twist angle within one cardiac cycle. The twist angle is defined as the rotation of the apex with respect to the basal plane at end-diastole. The apex from case LDDMM twists more compared with cases RBM¹⁷ and RBM^uni, with a peak value of 11°, which is well within the reported ranges in healthy hearts (10.2 ± 7.6°) [54]. Therefore, a more efficient pump function is achieved in case LDDMM compared with the RBM cases. Difference between the two rule-based cases are subtle, only slightly improved pump function can be found in case RBM¹⁷, compared with case RBM^uni, but it has a reduced apex twist.

Figure 7a–c shows the end-systolic myofibre stress distributions for the three cases. In case LDDMM, higher myofibre stress ( $\hat{\mathbf{f}}\cdot (\sigma \widehat{\mathbf{f})}$) can be found at both the endocardial and epicardial surfaces, especially in the LV side, while its distribution is less uniform compared with the two RBM cases. Figure 7d–f shows the strains along myofibre at end-systole. Strain distributions are similar in the two RBM cases, but the great difference is seen from the LDDMM case. The less uniform distributions of stress and strain in case LDDMM may be partially explained by much dispersed myofibre structures. The angle between the long-axis and the longitudinal axis at end-systole, defined in figure 7d–f, is largest in the LDDMM case (8.7°) and lowest in RMB^uni (4.2°), also suggesting different deformed end-systolic shapes.

• Download figure
• Open in new tab
• Download PowerPoint

3.2. RBM^uni with cross-fibre active tension

Based on case RBM^uni, five different sets of n_s and n_n are chosen to investigate how they affect ventricular dynamics. These are: (i) n_s = 0, n_n = 0, (ii) n_s = 0.2, n_n = 0, (iii) n_s = 0.4, n_n = 0, (iv) n_s = 0.0, n_n = 0.2 and (v) n_s = 0.0, n_n = 0.4. For all simulations n_f = 1.0. Figure 8 shows the pump functions with varied n_s or n_n. If we only consider cross-fibre active tension along the sheet direction, then the pressure–volume loop enclosed area is reduced as shown in figure 8a, suggesting that the active tension along the sheet direction will counteract the myofibre contraction. On the other hand, non-zero n_n increases the area enclosed by the pressure–volume loop and enhances the cardiac work. For example, with n_s = 0.4, the LV EF is around 29.97%, which is much less than the case with n_s = 0 (46.08%), while with n_n = 0.4, LV EF is increased by 10% as shown in figure 8b. Therefore, active tension along the sheet-normal direction is beneficial to the pump function, but contraction along the sheet direction has the opposite effect.

• Download figure
• Open in new tab
• Download PowerPoint

Figure 9 shows results from case RBM^uni with dispersed active contraction, modelled by the structural tensor from equation (2.14). In this case, case RBM^uni has nearly the same LV P–V loop as case LDDMM, and the apical twist is also similar to case LDDMM (figure 9a,b). Only a small difference in end-diastolic volume ( $\approx 1.4\text{\%}$) is observed between the two models. On the other hand, figure 9c,d shows that the end-systolic circumferential stress is much lower compared with case LDDMM, particularly in the midwall. The longitudinal and radial stresses are also slightly higher in the midwall because of non-zero n_n and n_s.

• Download figure
• Open in new tab
• Download PowerPoint

In summary, compared with case LDDMM, case RBM^uni shows a lower and more homogeneous stress level but achieves a similar pump function if using a suitable general structural tensor approach for the cross-fibre contraction.

It is interesting to see if similar results could be obtained without any knowledge of the patient-specific fibre field. To this end, we run extra simulations based on RBM^uni using literature-based values for n_f, n_s and n_n. Specifically, we consider (i) no dispersion n_f = 1, n_s = n_n = 0, (ii) n_f = 0.879, n_s = 0.009, n_n = 0.112 [45] and (iii) n_f = 0.646, n_s = 0.086, n_n = 0.268, derived from DT-MRI in this study. The fibre rotation angles are also chosen from 30° to −30° (exRBM¹), 45° to −45° (exRBM²), or 60° to −60° (exRBM³) [15]. The results are summarized in figure 10 in terms of the LV and RV ejection fractions. Clearly, EFs increase with fibre rotation angles, as more myofibres align longitudinally which enhance the active contraction. Different dispersion parameters also affect the pump function. Compared with case LDDMM, the EFs are lower in exRBM¹ (39.37% (LV), 45.89% (RV)), and still lower in exRBM² (47.86% (LV), 51.72% (RV)). Only exRBM³ with DT-MRI-derived dispersion parameters can achieve the similar pump functions as in case LDDMM, though the myofibre rotation angles (60° to −60°) are much greater than case LDDMM (mean angles 40° to −30°). This would suggest that subject-specific myofibre structure is necessary for cardiac mechanic modelling, as using literature-based myofibre structures seems to underestimate the pump function.

• Download figure
• Open in new tab
• Download PowerPoint

4. Discussion

In this study, LDDMM-based Deformetrica [41] is used to warp a canine bi-ventricle to a neonatal porcine heart, and DT-MRI-measured myofibre structure is then mapped to a porcine heart by solving a Laplace system. Base on the mapped DT-MRI measured myofibre structure, two simplified fibres are further generated using a rule-based approach. Our results show that under the same pre-/after-loading conditions, both LV and RV have a higher pump function in the case with LDDMM-mapped fibres compared with the rule-based cases, while case LDDMM experiences higher myofibre stress and more heterogeneous stress pattern than rule-based cases. Large differences can be expected when using literature-based fibre structures and dispersion parameters compared with case LDDMM. Those different results highlight the necessity of use realistic myofibre structure for personalized cardiac modelling as demonstrated in other studies [15,20,21,24].

In case LDDMM, the high active fibre stresses at both epicardial and endocardial surfaces (figure 7a) can potentially enhance the long-axis shortening and also apical twist (figure 6d). In fact, long-axis shortening in systole with respect to end-diastole is slightly higher in case LDDMM (−7.3%) than other two cases (−6.8% for RBM¹⁷ and −7% for RBM^uni). Our results show (figure 4a,b) that DT-MRI-derived myofibres do not lie in the c₀ − l₀ plane but are dispersed. Thus the active tension in case LDDMM is generated along fibres dispersed with both in-plane and out-of-plane components. In §2.2.3, we firstly quantify myofibre dispersion with in-plane and out-of-plane distributions, and then introduce a structural tensor H [26,55] by fitting to the measured in/out-of-plane dispersions. The π-period von Mises distribution is used to describe myofibre dispersion, good agreement can be achieved as shown in figure 4c,d. While it may not be guaranteed that the von Mises distribution can be applied to pathological tissues, such as myocardial infarction [56].

To take into account myofibre dispersion in the active stress formulation, we further assume the generated active force T_a is dispersed along the fibre, the sheet and the sheet-normal directions with proportions determined by H. We find that cross-fibre contraction is highest along the sheet-normal direction compared with that of the sheet direction, but much lower than along mean fibre direction. Furthermore, active contraction in the sheet-normal direction can facilitate contraction, but not in sheet direction. This is because myofibres dominantly lie in c₀ − l₀ plane, in which f and n are defined, and contraction along f and n causes circumferential and long-axial shortening [57], so the wall thickens to maintain the constant wall volume if the material is incompressible. While transmural contraction along s causes wall thinning, which counteracts myofibre contraction. Note that in this study, the sheet direction is defined transmurally across the wall, which is consistent with studies from [7,11,18], though some studies define it as the sheet-normal direction [6]. Unlike the myofibres which rotate from endocardium to epicardium, here the sheet direction is assumed to align the radial direction in all cases. In other words, the sheet rotation angle is chosen to be zero. To evaluate this assumption, we have tested three sets of sheet rotation angles as in [15]: 30° to −30°, 45° to −45° and 60° to −60°, based on case RBM^uni with dispersed active tension. The results show that the sheet rotation angle has little effect on ventricular pump function, and the differences in ejection fraction between different sheet rotation angles are within 1%. This agrees with observations from other groups. For example, Wang et al. [15] found that the sheet rotation angle nearly has no influence on passive mechanics in an LV model.

We now compare our values of cross-fibre proportions (n_s = 0.086, n_n = 0.268 and n_f = 0.646) with previous studies. Based on the experimental study by Lin & Yin [38], Guccione and co-workers introduced cross-fibre active contraction with n_s = 0.0, n_n = 0.4 and n_f = 1.0 [37]. In a recent study, Sack et al. [6] inversely determined cross-fibre contraction ratios⁵ in a healthy porcine heart (n_n = 0.07) and a failure heart (n_n = 0.14) with n_f = 1.0 and n_s = 0. In our study, n_n (0.268) is higher than that of Sack et al.'s study [6]. This could be due to (i) subject variation; (ii) higher n_f = 1.0 used in their study (our n_f = 0.646), leading to a higher contraction along the averaged myofibre direction so a lower n_n could match the measured pump function; (iii) they inversely determined n_n and T_a, which are not from measurements. In this study, proportions of cross-fibre contraction are derived directly from intrinsic fibre structures, which have a clear biological explanation. When normalized by n_f, the ratio between the sheet-normal and myofibre direction is 41%, which agrees with the ratio reported by Lin & Yin (40%) [38]. We further calculate the dispersion parameters from a recent study on neonatal porcine heart by Ahmad et al. [25], n_f = 0.68 and n_n = 0.32 with nearly negligible n_s ≈ 0.0009, again very close to our values in this study. We are not aware of any available experimental measurements for estimating n_n and n_s in the myocardium.

Rodríguez-Cantano et al. [24] argued that RBM tends to exaggerate myofibre-layered architecture and the passive stiffness of the ventricle, while DT-MRI-measured fibres may underestimate ventricular stiffness due to measurement noise and uncertainties. We find that when taking into account the cross-fibre contraction in the case RBM^uni, we can achieve similar systolic contraction as case LDDMM (figure 9) with less heterogeneous stress patterns. Because of challenging of in vivo DT-MRI acquisition, rule-based myofibre structures will continue to be used when modelling cardiac mechanics, even in personalized models. Our results suggest by incorporating fibre dispersion using a structural tensor, RBM-based model can be a good approximation of the most realistic myofibre structure as derived from DT-MRI, and the structural tensor may be determined either from limited in/ex vivo DT-MRI data [58] or inversely estimated, while caution needs to be paid when myofibre structures are from different subjects or species. There is a small difference (around 1.4%) in end-diastolic volume in figure 9a, presumably because the dispersion is not included in the passive constitutive law. Given that exclusion of compressed fibres using structural tensor approach is non-trivial in the passive modelling [59,60], we will leave the work of including dispersion in the passive model in future.

Using material parameters estimated from ex vivo measurements to describe in vivo material behaviours is a standing challenge. Published studies have suggested passive parameters estimated from ex vivo experiments can overestimate the stiffness in vivo [15,33,61]. Hence, most of the studies, ours included, scaled the parameters from ex vivo data to match the in vivo dynamics [5,6,61]. Here, the initial passive parameters are adopted from our previous study [39] which were inferred from ex vivo neonatal myocardial stretching experiments [40]; then a, a_f, a_n, a_fs and a_fn are scaled to achieve the targeted end-diastolic volumes. The myocardial contractility T_max is determined by matching the targeted ejection fractions (greater than 50%) for both the LV and RV. We further assume the passive scaling factor is the same for the LV and RV. Thus only three parameters need to be determined: the passive scaling factor, T_max for the LV and T_max for the RV. The sensitivity study on the passive parameters and T_max, and an illustration of their inferences are provided in the electronic supplementary material.

The convexity of the HO type strain energy function requires all parameters greater than zero as suggested in [7], which is satisfied in our approach. However, as pointed out by Giantesio et al. [62], the polyconvexity of the total energy function (passive and active) may not be ensured even though the passive strain energy function is convex. Although we have not experienced stability issues using the active stress approach, we must point out this approach may not be thermodynamically consistent. For generalized thermodynamically consistent approaches, the reader is referred to [9,34,62].

Due to lack of DT-MRI data for the RV from the canine experiment, a rule-based approach is used for generating fibre structure in the RV, and zero cross-fibre contraction is assumed. This can be readily improved if measured RV fibre structure becomes available. We note there is a difference in the RV systolic function even though the RV model is identical in all three cases. In particular, the RV contracts more in case LDDMM than in the two RBM cases. We think this is due to the different LV contraction in the three cases. For instance, the end-systole angle between the long-axis and longitudinal axis is different in each case. Palit et al. [19] also found that there are strong interactions between the LV and RV dynamics in diastole. This highlights the importance of LV-RV interaction on cardiac pump function, which is why the bi-ventricle model is used. In addition, the LDDMM framework [41] relies on geometrical features for warping the two different geometries, a bi-ventricular model has much richer information compared with a stand-alone LV model, in particular in the RV-LV insertion regions.

It is expected that there are differences in myofibre structure between the porcine heart and the canine heart, but this is difficult to assess as we do not have measured DT-MRI fibre structure for the porcine heart. However, despite the species difference, we find that the mapped canine myofibre structure agrees well with other studies in terms of mean values [6,23,25] (table 1). For instance, Ahmad et al. [25] measured myofibre rotation angles in LV free wall of neonatal hearts (anterior 51.1 ± 3.8° to −51.1 ± 3.8° and posterior 40.2 ± 2.9° to −40.2 ± 2.9°). Sack et al. [6] reported fibre rotation angles for a normal adult porcine heart based on DT-MRI measurements (endocardium: 66.5 ± 16.6°, epicardium: −37.4 ± 22.4°). Myofibre rotation angles from published experimental and numerical studies are also summarized in the electronic supplementary material.

The spatial variations of the material properties have not been considered in this study, and the same averaged dispersed active contraction model is applied across the whole LV for case RBM^uni. This approximation may be reasonable for healthy hearts, but questionable for pathological cases. For example, the myocardium is known to be more heterogeneous post-myocardial infarction [56].

Finally, we would like to mention other limitations of our study. In the boundary conditions we used, the basal plane of the models is constrained along the longitudinal direction, and the rest nodes in the basal plane are free to move. This type of boundary conditions does not represent in vivo conditions due to the lack of the pericardium and great vessels. Under in vivo situation, with the constraints imposed by the pericardium, the apex does not move much. Instead, the basal plane moves downward towards the apex in systole and moves upward in diastole. In a recent study, Pfaller et al. [63] demonstrated that simulated cardiac mechanics could be much closer to the measured heart motion by including the pericardium influences, which highlights the necessity of pericardial–myocardial interaction. A simplified lumped circulation model is used to provide pressure boundary conditions, which is a simplification of pulmonary and systemic circulations. Coupling to a more realistic circulation model, such as one-dimensional systemic models [64,65], will allow us to simulate more detailed cardiovascular function in pathological situations [66]. Furthermore, we have not coupled the blood flow inside ventricle, only applied a spatially homogeneous pressure to the endocardial surface, nor have we considered contraction delay due to the action potential propagation [9]. Tremendous efforts will be needed to address all those limitations, which is beyond the scope of this study.

5. Conclusion

In this study, we have developed a bi-ventricular porcine heart computational model from a neonatal dataset, with mapped myofibre architecture from an ex vivo canine DT-MRI dataset using an LDDMM framework. Different approximations of myofibre architecture based on widely used rule-based approaches are analysed in terms of cardiac pump function. Our results show that using DT-MRI derived myofibre architecture can enhance cardiac work, achieve higher ejection fraction and larger apical twist compared with rule-based myofibre models, even though they are all derived from the same DT-MRI dataset. Our work shows that the major difference between the LDDMM and RMB approaches is due to the fibre dispersion, which enables cross-fibre active tensions. These are not captured by standard RBM-based models. Introducing regional dependent fibre structure in RBM is not sufficient to improve the model. However, when the myofibre dispersion is taken into consideration, a simplified RBM-based cardiac model can achieve similar pump function as the LDDMM-based model. We further note that in RBM-based cardiac models, the cross-fibre active tension along the sheet-normal direction can enhance active contraction, but the opposite is true along the sheet direction.

Data accessibility

Electronic supplementary material is provided. The datasets supporting this article have been uploaded to github as part of the electronic supplementary material, https://github.com/HaoGao/FibreGeneration-LDDMM, and the DOI is 10.5281/zenodo.3458254.

Authors' contributions

D.G. performed the modelling, analysed results and drafted the manuscript. J.Y. contributed to the model development. H.G. and X.L. conceived the study, analysed the results and supervised the project. H.G. coordinated the study. All authors contributed to the writing of the manuscript and gave final approval for publication.

Competing interests

Prof. Xiaoyu Luo is a member of the board of Royal Society Open Science at the time of submission.

Funding

We are grateful for the funding provided by the UK EPSRC (grant no. EP/N014642/1). D.G. also acknowledges funding from the Chinese Scholarship Council and the fee waiver from the University of Glasgow. H.G. also thanks for the funding from the National Natural Science Foundation of China (grant no. 11871399). Thanks also extend to all members of the Living Heart Project.

Footnotes

Electronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.4915563.