Effects of dispersed fibres in myocardial mechanics, Part II: active response

1.   Introduction

In this part, we will study how myofibre rotation and its dispersion affect left ventricle (LV) active contraction. As demonstrated in many studies ^[1,2,3,4] and also the Part I ^[5], fibre dispersion can have significant effect on passive behaviours of the myocardium. Here we will mainly review studies on myocardial active contraction with myofibre dispersion. For a general introduction of myocardial mechanics and fibre dispersion in passive constitutive modelling, readers can refer to Part I ^[5] and recent reviews ^[6,7].

For the first time, Lin and Yin ^[8] reported that there was 40% cross-myofibre active stress in their biaxial contraction tests in rabbit myocardium. Since then, modelling studies began to account for cross-fibre active contraction using the active stress approach, in which total stress tensor consists of passive and active parts. For example, Guccione and co-workers included cross-myofibre active stress in the total stress to better match in vivo cardiac pump function ^[9,10]. Later, Sack et al. ^[11] added cross-myofibre active tension with cross-myofibre contraction ratios that were inversely determined using a healthy porcine heart and a failing heart. Recently, Guan et al. ^[1] investigated the effects of cross-myofibre active stress on cardiac pump function in which a generalised structural tensor (GST) based dispersed active tension model was developed by assuming myofibres follow a non-fully dispersed distribution, which was informed by experimental data, such as the diffusion-tensor magnetic resonance imaging ^[12]. They found that the dispersed active stress along the sheet direction, i.e., the transmural direction from endocardium to epicardium, would reduce the cardiac pump function, whilst the active stress along the sheet-normal direction enhanced the pump function. The GST-based active stress approach was also used by Eriksson et al. ^[13], which was based on the $\kappa$–model ^[3] by assuming myofibres following a fully dispersed distribution. While experimentally measured data demonstrates that the in-plane and out-of-plane myofibre dispersion is different ^[14], suggesting an anisotropic myofibre dispersion will be more appropriate to be considered for modelling dispersed active stress, such as the non-fully dispersed distribution.

There are a few studies ^[15] on how myofibre rotation angle variation affects LV passive behaviour, but very limited on active contraction. In our previous study ^[1], we derived the dispersed active tension from experimentally measured data and studied the effects of a few extremely dispersed distributions on myocardial contraction. While there is no existing study about how myofibre rotation angle variation affects LV contraction, and how the dispersed active tension in the sheet-normal direction affects longitudinal contraction, we will address the two questions mentioned above in this study, firstly by varying myofibre rotation angles from $-40^\circ \sim 40^\circ$ to $-80^\circ \sim 80^\circ$, and secondly by increasing the sheet-normal active tension proportion from 0 to its analytical maximum using the GST-based active tension model ^[1].

2.   Methods

Similar to Part I ^[5], dispersed myofibre (${\bf{f}}_n$) in a local fibre system (${\bf{f}}_0, {\bf{s}}_0, {\bf{n}}_0$) can be defined by the two spherical polar angles $\Theta$ and $\Phi$,

and the domain of dispersed myofibres at each location is an unit hemisphere with $\mathbb{S} = \left\{(\Theta, \Phi) | \Theta \, \in \left[0, \pi/2\right], \Phi\, \in\left[ 0, 2\pi \right] \right\}$. By assuming the myofibre dispersion distribution following a probability density function, then we have

in which $\rho_{\rm{in}}(\Theta, b_1)$ describes the in-plane dispersion, $\rho_{\rm{op}}(\Phi, b_2)$ describes the out-of-plane dispersion, and $G$ is a constant. We choose the $\pi$-periodic Von Mises distribution for $\rho_{\rm{in}}$ and $\rho_{\rm{out}}$, and please refer to Part I ^[5] for further details.

2.1.   Passive constitutive modelling

Following the results from Part I ^[5], the passive strain energy function for the myocardium is

where $a_{({\rm{g, f, s, fs}})}$, $b_{({\rm{g, f, s, fs}})}$ are material parameters, $I_1 = {\rm{trace}}({\bf{F}}^{\rm{T}} {\bf{F}})$, ${\bf{F}}$ is the deformation gradient tensor, and $I_{\rm{8fs}} = {\bf{f}}_0 \cdot ({\bf{F}}^{\rm{T}}{\bf{F}}\, {\bf{s}}_0)$, in which ${\bf{f}}_0$ is the mean myofibre direction and ${\bf{s}}_0$ is the sheet direction, i.e., the transmural direction from endocardium to epicardium. $\Psi^{I_{\rm{4f}}*}_{\rm{aniso}}$ and $\Psi^{I_{\rm{4s}}*}_{\rm{aniso}}$ are the sums of strain energy of all dispersed myofibres (${\bf{f}}_n$) and sheet fibres (${\bf{s}}_m$), respectively, by using the discrete fibre dispersion method ^[2], in which $I_{\rm{4f}}^n = {\bf{f}}_n \cdot ({\bf{F}}^{\rm{T}} {\bf{F}}\, {\bf{f}}_n)$, $I_{\rm{4s}}^m = {\bf{s}}_m \cdot ({\bf{F}}^{\rm{T}} {\bf{F}}\, {\bf{s}}_m)$, $\mathcal{H}(\cdot)$ is the Heaviside function to ensure the only stretched fibres can bear the load, and discrete probability $\rho_n$ of ${\bf{f}}_n$ and $\rho_m$ of ${\bf{s}}_n$ are determined by their respective probability density functions. Further details can be found in Part I ^[5].

2.2.   Active stress modelling

If without myofibre dispersion, then the generated active stress is only along ${\bf{f}}$, that is

where $T_{\rm{a}}$ is the active tension and $I_{\rm{4f}} = {\bf{f}} \cdot {\bf{f}}$ is the square of myofibre stretch with ${\bf{f}} = {\bf{F}}{\bf{f}}_0$. Cross-myofibre active contraction has been observed in the biaxial active contraction experiment conducted by Lin and Yin ^[8], which has inspired modelling studies ^[9,16] to include cross-myofibre active contraction along with the myofibre contraction. Here, the dispersed active stress approach ^[1,13] is employed to take into account myofibre dispersion, that is

in which $\hat{{\bf{f}}} = {\bf{F}}{\bf{f}}_0/|{\bf{F}}{\bf{f}}_0|$, $\hat{{\bf{s}}} = {\bf{F}}{\bf{s}}_0/|{\bf{F}}{\bf{s}}_0|$ and $\hat{{\bf{n}}} = {\bf{F}}{\bf{n}}_0/|{\bf{F}}{\bf{n}}_0|$ are unit vectors in the current configuration, and ${\bf{n}}_0$ is the sheet-normal direction and defined by ${\bf{n}}_0 = {\bf{s}}_0 \times {\bf{f}}_0$, $n_{\rm{f}}$, $n_{\rm{s}}$ and $n_{\rm{n}}$ can be derived from a generalised structure tensor ${\bf{H}}$ that is defined over an unit hemisphere using the predefined probability density function Eq (2.2) ^[3,13,17],

where $H_{11}, H_{22}, H_{33}$ are diagonal components with respect to the orthogonal material axes (${\bf{f}}_0$, ${\bf{s}}_0$ and ${\bf{n}}_0$). Thus we have

and $n_{\rm{f}}+n_{\rm{s}}+n_{\rm{n}} = 1$. The tensor ${\bf{H}}$ is diagonal because of the $\pi-$period Von Mises probability density functions for the fibre dispersion and the independence assumption between the in-plane and out-of-plane dispersion. In other words, the myofibre dispersion distribution is symmetric with respect to (${\bf{f}}_0$, ${\bf{s}}_0$) plane, (${\bf{f}}_0$, ${\bf{n}}_0$) plane and (${\bf{s}}_0$, ${\bf{n}}_0$) plane. We would like to mention that the tensor ${\bf{H}}$ may not be diagonal for general fibre dispersion.

The active tension ($T_{\rm{a}}$) is determined by a well-established time-varying elastance model ^[11,18],

in which $T_{\rm{max}}$ is the maximum isometric active tension, ${\rm{Ca}}_0$ is the peak intra-cellular calcium, and the length-dependent calcium sensitivity (${\rm{ECa}}_{50}$) is given by

where $B$ and ${\rm{Ca}}_{\rm{0max}}$ are constants, $l_0$ is the minimum sarcomere length to produce active stress, and $l$ is the deformed sarcomere length

in which $l_r$ is the stress-free sarcomere length, and $E_{\rm{ff}}$ is the Lagrange strain in the myofibre direction. The time function after onset of contraction in Eq (2.8) is

where $t_0$ is the time to peak tension, and $t_r$ is the duration of relaxation

where $m$ and $b$ are constants. Finally, the total myocardial Cauchy stress is formulated using the active stress approach,

where $p$ is the Lagrange multiplier to enforce the incompressibility and ${\bf{I}}$ is the identity tensor.

2.3.   In vivo left ventricle modelling

The same human LV model ^[19] as in the Part I is used here to study how different myofibre dispersion and rotation can affect the LV pump function, as shown in Figure 1(a). The mean myofibre direction at each location is generated using a rule-based method ^[15], in which the myofibre rotation angle varies linearly from the epicardium ($\bar{\Theta}^{\rm{epi}}$) to the endocardium ($\bar{\Theta}^{\rm{endo}}$). A simplified close-loop lumped-parameter circulation system including the aorta and the left atrium is coupled to the LV model to provide physiologically accurate pressure boundary conditions, which is implemented using the surface-based fluid cavities and the fluid exchanges ^[20] as shown in Figure 1(a). This type of boundary condition is equivalent to a two-element Windkessel model ^[21], and more details can be found in Appendix.

LV simulations are performed by using the nonlinear finite-element software ABAQUS EXPLICIT (Dassault Systemes, Johnston RI, USA). In brief, the dynamic system of the LV model is

where $\mathbb{M}$ is the mass matrix, $\mathbb{C}$ is the damping matrix, $\mathbb{K}({\bf{x}})$ is the nonlinear stiffness matrix depending on the current configuration, $\mathbb{F}_{\rm{ext}}$ is the external forces, and ${\bf{x}}$ is the displacement field. To eliminate unrealistic transient behavior, the damping matrix is introduced via Rayleigh damping ^[20], i.e.,

where $\alpha$ is the damping factor. We find that $\alpha = 160 {\rm{s}}^{-1}$ can simulate LV dynamics without noticeable oscillation, which is also the value used in the LivingHeart Project ^[16].

The material parameters of the passive strain energy function are re-scaled from the fitted parameters in Part I ^[5] following the approach proposed in ^[22,23], that is

where $a_j^0$ and $b_j^0$ are the parameters of Case 4 in Part I ^[5]. We find that the filling volume in diastole is around 60 ml when $\alpha = 0.4$ and $\beta = 0.5$, which lies in the normal physiological range of human cardiac function. Table 1 summarizes the passive constitutive parameters of the myocardium in this study, along with active contraction parameters. Values of $t_0$, $m$, $b$, $l_0$, $B$, $l_r$, ${\rm{Ca}}_0$ and ${\rm{Ca}}_{0{\rm{max}}}$ are referred from ^[11,18], and $T_{\rm{max}}$ is inversely determined by achieving an ejection fraction within $50$–$75\%$. Note in this study, we do not intend to simulate personalized LV dynamics, which requires extra measurements for presonalization ^[19].

To simulate LV dynamics, we first preload the LV to 8\, mmHg within $0.5\, s$ to reach the end of diastole, then the iso-volumetric contraction begins as the whole LV is contracting simultaneously, then followed by the systolic ejection at $t\approx 0.6\, s$ when the LV cavity pressure exceeds the aortic pressure (about 80 mmHg), and finally the ejection ends when the LV pressure is lower than the aortic pressure.

To describe LV contraction, we further introduce the long-axis shortening, radial shortening and apex twist angle with respect to the end of diastole (beginning of systole), as shown in Figure 1(b). The long-axis length is the distance between the centre at the base plane and the apex point at the inner surface, then the long-axis shortening ($L_r$) is

in which $L_0$ and $L_t$ are the long-axis length at the beginning of systole and the current time $t$, respectively, and the radial-axis shortening ($R_r$) is the mean radius change of the inner wall at the middle LV,

where $R_i^0$ is the mean radius at the end of diastole, and $R_i^t$ is the corresponding value at time $t$. Apex twist angle ($\theta$) is the average angle change of the apex in related to the end of diastole.

3.   Results

In this section, we first investigate the LV pump function with different myofibre rotation angles in Section 3.1 with the same myofibre dispersion ($b_1 = 4.5$ and $b_2 = 3.9$). The values of the corresponding scaling factors $n_{\rm{f}}$, $n_{\rm{s}}$ and $n_{\rm{n}}$ in Eq (2.5) are 0.879, 0.009 and 0.112 obtained from ^[1]. Then, the effect of myofibre dispersion on the cardiac pump function is studied in Section 3.2 with the same myofibre rotation angle whilst values of $b_1$ and $b_2$ vary according to the predefined values of $n_{\rm{f}}$, $n_{\rm{s}}$ and $n_{\rm{n}}$. Please note the sheet dispersion remains the constant fully dispersed distribution with $b_3 = 4.5$ and $b_4 = 0.0$.

3.1.   Transmural myofibre rotation

Seven cases are considered here with different myofibre rotation angles from the epicardium to the endocardium, which are

● rot-80: $-40^\circ \sim 40^\circ$, and the rotation angle: 80$^\circ$;

● rot-100: $-50^\circ \sim 50^\circ$, and the rotation angle: 100$^\circ$;

● rot-120: $-60^\circ \sim 60^\circ$, and the rotation angle: 120$^\circ$;

● rot-140: $-70^\circ \sim 70^\circ$, and the rotation angle: 140$^\circ$;

● rot-160: $-80^\circ \sim 80^\circ$, and the rotation angle: 160$^\circ$;

● rot-unsym1: $-30^\circ \sim 90^\circ$, and the rotation angle: 120$^\circ$;

● rot-unsym2: $-90^\circ \sim 30^\circ$, and the rotation angle: 120$^\circ$,

and the corresponding myofibre structures are shown in Figure 2. Cases rot-80 to rot-160 have the same angle to the circumferential direction with opposite signs at the epicardium and the endocardium, denoted as symmetrical cases, and Cases rot-unsym1/2 have different angles to the circumferential direction at the epicardium and the endocardium, denoted as unsymmetrical cases.

Figure 3(a) shows the Pressure-volume (PV) loops of the seven cases, and the corresponding values of end diastolic volume (EDV), end systolic volume (ESV) and ejection fraction (EF) are summarised in Table 2. For the symmetric rotation cases (rot-80 to rot-160), myocardial contraction is significantly improved when increasing myofibre rotation angle from $80^\circ$ (rot-80) to $100^\circ$ (rot-100), with a further decrease of 5.9 ml in ESV from 55.8 ml of the Case rot-80. Less improvement in myocardial contraction with further increased myofibre rotation angle, the smallest ESV (46.1 ml) happens with $140^\circ$ rotation angle. Interestingly, EDV increases 1.0 ml per $20^\circ$ increment of the myofibre rotation angle. Except for EF in Case rot-80 that is 46%, Cases rot-100 to rot-160 achieve an EF of greater than 50% with the biggest one (56.5%) occurring in Case rot-140.

Long-axis shortening is shown in Figure 3(b) where the greater negative value suggests larger longitudinal contraction. Again, for the symmetric rotation cases, the long-axis shortening ratio at the end of systole ($t = 0.2\, s$) increases with increased myofibre rotation angle from Case rot-80 (5.0%) to Case rot-160 (-16.0%). Interestingly, for Cases rot-80 and rot-100, due to the small rotation angle, the long-axis shortening ratio is positively increasing in the early-systole rather than shortening compared to other cases. This long-axial elongation during systole suggests that the LV does not contract longitudinally compared to the end-diastolic shape, which means the LV contraction may not physiologically accurate. The radial shortening is similar for all cases as shown in Figure 3(c). The rotational angle affects the radial shortening magnitude, the less the rotational angle, the more the radial shortening, which can be explained by more myofibres aligned along the circumferential direction, leading to deeper contraction circumferentially, thus more radial shortening. Figure 3(d) describes the LV apex twist, in which the small rotation angles can cause large twists, such as from 34.2$^\circ$ in Case rot-80 to 20.0$^\circ$ in Case rot-160. Summary of LV pump function can be found in Table 2, including EDV, ESV, EF, etc.

The unsymmetrical rotation angle can also significantly affect cardiac pump function as demonstrated by Cases rot-unsym1 and rot-unsym2 that have the same rotation angle as Case rot-120. Compared to Case rot-120 with 55.6% of EF, both unsymmetric cases have reduced contractility with smaller EFs (50.0 % in Case rot-unsym1 and 46.7% in Case rot-unsym2) due to much larger ESVs. The large myofibre angle at the endocardium in case rot-unsym1 leads to a slight increase in EDV (108.1 ml) compared to Case rot-120 (105 ml), whilst a slightly reduced EDV occurs in Case rot-unsym2 (104.4 ml). Small myofibre angle at the endocardium in Case rot-unsym2 achieves both higher radial (Figure 3(c)) and long-axis shortening (Figure 3(b)) compared to Case rot-unsym1 at the end of systole, whilst they both experience less circumferential contraction and larger longitudinal contraction compared to Case rot-120, which may be explained by more longitudinally aligned myofibres in those two cases. Furthermore, the unsymmetrical myofibre rotation significantly reduces the apex twist as shown in Figure 3(d), in which the apex twist angle at the end of systole is 3.0$^\circ$ for Case rot-unsym1 and 16.3$^\circ$ for Case rot-unsym2, much less than that of Case rot-120 (25.8$^\circ$).

Distributions of end-systolic myofibre stress (${\boldsymbol{\sigma }}_{\rm{ff}}$) with deformed geometries are shown in Figure 4, and mean values are listed in Table 2. The mean myofibre stress ${\boldsymbol{\sigma }}_{\rm{ff}}$ slightly increases from Case rot-80 (33.8 $\pm$ 31.1 kPa) to Case rot-160 (37.0 $\pm$ 32.4 kPa) in Figure 5. Compared to Cases rot-120/140/160, higher myofibre stress can be found at the basal region of LV in Cases rot-80, rot-100 and rot-unsym2 with more heterogeneous stress distributions. Corresponding statistic density distributions of ${\boldsymbol{\sigma }}_{\rm{ff}}$ are shown in Figure 5. It can be found that Case rot-unsym1 and Case rot-unsym2 both have more heterogeneous stress distributions compared to other cases, while Case rot-120 has the least standard deviation suggesting the most homogeneous distribution, see Figure 5.

In sum, myofibre rotation angles can significantly affect LV pump function. For example, a myofibre structure with a large rotation angle can lead to greater longitudinal contraction than with a small rotation angle. For the human LV model, it shows the myofibre rotation angle in the range of $120^\circ \sim 140^\circ$ can achieve a more physiologically accurate pump function, i.e., EF $>$ 50%. In fact, the myofibre rotation angle from 60$^\circ$ to -60$^\circ$ is the most widely used range in the literature ^[15,19,24], and also within the range reported by experimental studies ^[25].

3.2.   Sheet-normal active tension

The reduced myocardial contraction due to the dispersed active tension along the sheet direction has been studied in ^[12], in which the LV ejection fraction was reduced by 12.8% when increasing $n_{\rm{s}}$ from 0 to 0.2 with a constant $n_{\rm{f}} = 1$. One reason is that since the myocardium is incompressible, the active tensions along the myofibre and the sheet-normal directions need to overcome the non-zero active tension along the sheet direction before the myocardium can contract. Our previous study also has shown that $n_{\rm{s}} = 0.008$ based on Sommer et al.'s measurements from human myocardium ^[25] and $n_{\rm{s}} = 0.0009$ based on Ahmad et al.'s measurements from neonatal porcine myocardium ^[14]. Both values are very close to zero, thus we assumed the active tension along sheet direction is negligible by setting $n_{\rm{s}} = 0$. Based on Eq (2.6), value of $n_{\rm{s}} = 0$ drops to less than 0.01 for any $b_1$ when $b_2\geq16$, thus we set $b_2 = 16$, whilst values of $n_{\rm{f}}$ and $n_{\rm{n}}$ are determined by $b_1$. Finally, five cases of myofibre dispersion are compared as following

● Case f10n0: $n_{\rm{f}} = 1.0$, and $n_{\rm{n}} = 0.0$, corresponding to no-dispersion;

● Case f8n2: $n_{\rm{f}} = 0.8$, and $n_{\rm{n}} = 0.2$, ($b_1 = 2.9$, and $b_2 = 16$);

● Case f6n4: $n_{\rm{f}} = 0.6$, and $n_{\rm{n}} = 0.4$, ($b_1 = 1.35$, and $b_2 = 16$);

● Case f5n5: $n_{\rm{f}} = 0.5$, and $n_{\rm{n}} = 0.5$, ($b_1 = 0.85$, and $b_2 = 16$);

● Case f4n6: $n_{\rm{f}} = 0.4$, and $n_{\rm{n}} = 0.6$, ($b_1 = 0.35$, and $b_2 = 16$).

Note the minimum value of $n_{\rm{f}}$ is 0.33 when $b_1 = 0.0$ and $b_2 = 16$, thus we only consider $n_{\rm{f}} \geq 0.4$ in this study. The corresponding myofibre dispersion distributions ($b_1$ and $b_2$ in Eq (2.2)) are determined according to the values of $n_{\rm{f}}$, $n_{\rm{s}}$ and $n_{\rm{n}}$. Myofibre structures of those above cases are the same as that of Case rot-120 ($-60^\circ \sim 60^\circ$) in Section 3.1 but with different in-plane and out-of-plane dispersion.

Figure 6(a) shows the PV loops of the LV model with different dispersed active tension, and the related EDV, ESV and EF are summarised in Table 3. In general, the case with large sheet-normal active tension has small ESV and large EDV. For example, ESV in Case f4n6 is 31.6 ml, which is much smaller than that in Case f10n0 with an ESV of 45.9 ml. Similarly, the long-axis shortening is greater with increased $n_{\rm{n}}$ as shown in Figure 6(b). Long-axis shortening ratios at the end of systole are negative for all cases, from -8.1% in Case f10n0 to -39.2% in Case f4n6. Conversely, the radial-shortening is smaller with more active tension along the sheet-normal direction, from -34.3% in Case f10n0 to -27.0% in Case f4n6 as shown in Figure 6(c). It also can be found that the sheet-normal contraction can prevent long-axis elongation during systole in Case f10n0.

Apex twist is significantly influenced by the sheet-normal active tension as shown in Figure 6(d). When $n_{\rm{f}}$ is greater than $n_{\rm{n}}$, the apex twist increases to the peak value quickly at early-systole and then decreases slightly at late-systole, while the larger $n_{\rm{n}}$ in Cases f5n5 and f4n6 leads to continuously increased apex twist in systole.

Although very strong pump function can be found in Case f5n5/f4n6 according to the PV loops in Figure 6(a), excessive longitudinal contraction with compressed myofibres around the endocardium can be found in Figure 7, which is not physiologically correct. Compared to Case f10n0, Case f8n2 has the more homogeneous stress distribution at the base with greater long-axis shortening. Figure 8 shows the statistics of stress distributions for the five cases, in which Case f8n2 has the most homogeneous myofibre stress distribution. Case f6n4 achieves the better pump function with smaller myofibre stress than Case f10n0 and Case f8n2 while maximum ${\boldsymbol{\sigma }}_{\rm{ff}}$ locating at the middle LV wall and preventing the myofibre compression in Cases f4n6 and f5n5. Mean values of ${\boldsymbol{\sigma }}_{\rm{ff}}$ of the whole heart at end-systole are summarised in Table 3. The myofibre stress ${\boldsymbol{\sigma }}_{\rm{ff}}$ reduces dramatically from 38.5 kPa in case f10n0 to 4.5 kPa in Case f4n6 (Figure 8), which is presumably because the active tension aligns more along the sheet-normal direction, less along the myofibre direction.

In sum, we find that $n_{\rm{n}} > 0$ can lead to greater long-axis contraction, smaller ESV, and much-reduced myofibre stress at systole compared to the case without dispersion. However, using larger $n_{\rm{n}}$, i.e., $n_{\rm{n}}\ge 0.4$, in the LV model, excessive long-axis contraction can occur which would not be physiological correct.

4.   Discussion

In this study, we have investigated the effects of transmural myofibre rotation angles and their dispersion distributions on cardiac pump functions, and we demonstrate that they both can significantly affect the LV pump function. Firstly, five symmetric transmural myofibre rotation angles, from $-40^\circ \sim40^\circ$ to $-80^\circ \sim80^\circ$, are compared, and the results show that the larger myofibre rotation angle can enhance myocardial contractility with greater longitudinal contraction compared to cases with smaller myofibre rotation angle, while the two unsymmetric myofibre rotation cases not only have reduced pump function but also reduced apex twist, especially with a small rotation angle at the epicardium. Secondly, the results from different combinations of $n_{\rm{f}}$ and $n_{\rm{n}}$ show that dispersed active tension along the sheet-normal direction can significantly increase longitudinal contraction, even beyond the physiological range when its portion is greater than 0.4. Thus, myofibre rotation and dispersion both can moderate cardiac pump function by adjusting the longitudinal and circumferential active tension, which have not been quantified in other studies. This also would suggest that a personalized cardiac model will need to take into account fibre dispersion and rotation angle variation rather than assuming a general fibre structure without dispersion, such as the widely used $-60^\circ \sim 60^\circ$ myofibre rotation angle ^[10,15,22].

In the cases of myofibre rotation, the dispersion is constant, which means the cross-fibre active tension is the same in those considered cases. Different myofibre rotation angles will change the active tension along the circumferential and longitudinal directions, and further affect the contraction along the two directions accordingly. For example, a larger myofibre rotation angle leads to more myofibres aligning along the longitudinal direction and a larger fraction of active tension longitudinally, thus greater long-axis shortening. It can be found the maximum contraction can be achieved as shown in Figure 3, when the myofibre rotation angle is in the range of $120^\circ \sim140^\circ$ with the same angles at the endocardium and the epicardium.

The myofibre rotation angle $-60^\circ \sim60^\circ$ is widely used in rule-based myofibre structure (Case rot-120 in this study) in many studies ^[1,15,19]. Our results show that such a myofibre structure can achieve a physiologically accurate LV pump function, which may suggest that if experimentally measured fibre structure is not available, $-60^\circ \sim60^\circ$ could be a good approximation. For the cases with unsymmetric myofibre rotation angles, both cases have reduced pump function, i.e., larger ESV compared to Case rot-120 that has the same and symmetric myofibre rotation angle. It may suggest that a remodelled myofibre structure leads to inefficient myocardial contraction under pathological conditions, for example, more circumstantially aligned myofibres ^[26].

For myofibre dispersion cases, $n_{\rm{n}}$ determines the longitudinal proportion of active tension and can significantly affect the long-axis contraction. Large $n_{\rm{n}}$, in other words, small $n_{\rm{f}}$ leads to reduced myofibre stress and increased long-axis shortening, and thus reduced ESV. This may suggest that the long-axial shortening is highly sensitive to $n_{\rm{n}}$. Case f8n2 has the similar myofibre dispersion as the measured data from the human heart ^[25] and achieves the similar PV loop as Case f10n0 without myofibre dispersion, which is consistent with the results in ^[1], but with a more homogeneous stress distribution and greater long-axial shortening. We further find that a much-enhanced contraction can also be achieved when $n_{\rm{n}}$ is in the range of $0.2\sim0.4$. In fact, $0.4$ is the cross-fibre active tension portion reported by Lin and Ying ^[8]. In our LV model, a 6% increase of $T_{\rm{max}}$ is required for the case f10n0 in order to match the pump function of the case f6n4. Non-fully dispersed myofibre distribution is used in this study, because the experimental data measured by Sommer et al. ^[25] and Ahmad et al. ^[14] has demonstrated myofibres disperse differently in-plane and out-of-plane. Our previous study ^[12] showed that active tension along the sheet direction could reduce myocardial contraction, which should be minimal in a healthy heart. Thus, we mainly focus on the variation of $n_{\rm{f}}$ and $n_{\rm{n}}$ with $n_{\rm{s}} = 0$ by setting a large value of $b_2 = 16$.

Together with Part I ^[5], we have comprehensively studied how fibre dispersion affects LV dynamics in diastole and systole, respectively. Our results suggest it is necessary to incorporate fibre dispersion in the passive constitutive law of the myocardium, both the myofibre and sheet fibre dispersion. Comparable results of LV passive dynamics can be achieved with only myofibre dispersion given that the sheet fibre dispersion can be challenging to measure experimentally. In this study, we focus on the effects of fibre dispersion on myocardial contraction. Our results further suggest that it is essential to include realistic fibre structure with both realistic transmural myofibre rotation angle and its dispersion.

The active stress approach (Eq (2.5)) was employed to simulate active stress response in the myocardium, which included the effects of myofibre dispersion using the generalised structure tensor approach, from which three scalar factors ($n_{\rm{f}}$, $n_{\rm{s}}$ and $n_{\rm{n}}$) can be analytically derived to characterize dispersed active tension along each material axis. To take into account length-dependency, the active tension $T_a$ was computed using the myofibre stretch along the mean fibre direction (Eqs (2.8)–(2.12)). We further tested one simulation when the active tension for each dispersed myofibre was calculated using its own stretch rather than the mean myofibre stretch, the dispersed active stress tensor was then evaluated using the discrete fibre dispersion approach as done in the passive response, and the results showed that the LV pump function was almost identical to the simulation using the GST approach. One advantage of using the GST-based active tension model is that it allows analytical integration of dispersed active tension with a given distribution density function, thus much higher computational efficiency can be achieved.

In the ABAQUS implementation, we treated the myocardium as nearly incompressible through a Lagrange multiplier which is derived from an additional term ($\frac{1}{D}(\frac{J^2-1}{2}-\ln (J))$) added to the strain energy function with a multiple of bulk modulus $D$. In the implementation, the deformation gradient tensor (${\bf{F}}$) was further decomposed to a deviatoric part ($\bar{{\bf{F}}} = J^{-1/3} {\bf{F}}$) and a volumetric part, and the volumetric part was used to compute the Lagrange mutliplier to ensure the incompressibility constraint. In our simulations, $J = 1\pm0.008$ at the end of diastole and $J = 1\pm0.009$ at the end of systole, suggesting the incompressibility was ensured. We did not meet the volumetric locking problem in all simulations.

We would like to mention the limitations of this study. Firstly, a simple circulation system is attached to the LV model to achieve physiologically-realistic cardiac boundary conditions. The hemodynamics characteristics in the aorta and left atrium are ignored, which can be further included to investigate the interaction between ventricles and blood flow ^[21,27]. Furthermore, electrophysiology is not modelled here because it will dramatically increase the complexity of the model ^[28]. Considering that the propagation of action potential is much faster than the mechanical contraction, the assumption of simultaneous contraction in the whole LV is feasible as done in many other studies ^[11,12,19].

5.   Conclusions

In this study, we have systemically investigated the effects of myofibre rotation and its dispersion on cardiac pump functions with a focus on myocardial contraction. Our results suggest that they both can significantly affect LV pump function, especially myofibre dispersion. Large transmural myofibre rotation angle can lead to enhanced long-axis shortening, while a small rotation angle can lead to enhanced circumferential contraction due to more myofibres aligned circumferentially. Myofibre dispersion redistributes active tension along the myofibre and cross-fibre directions, which determines how the LV contacts, especially when the dispersed active tension along the sheet-normal direction is greater than 0.4 for which excessive long-axis shortening can happen. Our two studies have highlighted that it is necessary to include realistic fibre structures in personalized cardiac modelling.

Acknowledgements

We are grateful for the funding provided by the UK EPSRC (EP/S030875, EP/S020950/1, EP/S014284/1, EP/R511705/1) and H. G. further acknowledges the EPSRC ECR Capital Award (308011). L. C. acknowledges the National Natural Science Foundation of China (11871399, 11471261, 11571275). D. G. and Y. W. also acknowledge funding from the Chinese Scholarship Council and the fee waiver from the University of Glasgow. Many thanks to Mr. Yuzhang Ge for proofreading.

Conflict of interest

The authors declare that there is no conflict of interest.

Data accessibility

The datasets supporting this article have been uploaded to GitHub as part of the electronic supplementary material, https://github.com/HaoGao/FibreDispersionMyocardialMechanics.git.

Appendix

A.   Analogy to a Windkessl model

A simplified close-loop circulatory model is implemented here by using the fluid-cavity and exchange modules in ABAQUS ^[20], which is similar to a Windkessel model. Specifically, by neglecting the inertia effects in the circulatory model, the flow rate between the two neighboring cavities ($Q_{c\rightarrow c+1}$) is driven by ^[20],

in which $\Delta p$ is the pressure difference between the two connected cavities, $A$ is the effective area of the fluid exchange, $C_v$ is the viscous resistance coefficient, and $C_H$ is the hydrodynamic resistance coefficient that is 0 in this study. Let $\xi$ be the blood density, then $\dot{m} = \xi \, Q_{c\rightarrow c+1}$ and

that is

Thus, the fluid exchange corresponds to a resistor in the Windkessel model, and its resistance is $\xi C_v/A$.

For a fluid cavity with one face fixed and the opposite one attached to a linear spring. The stiffness of the linear spring is $k$, and the area of the attached face is $\tilde{A}$. Again by ignoring the inertia effect and applying the force balance condition to the face attached to the spring at current time ${\rm{tc}}$, we have

where $P$ is the fluid cavity pressure, $\Delta x_{\rm{tc}}$ is the distance the spring is compressed, $V$ is the fluid cavity volume. Similarly for the next time step ${\rm{tn}}$, we have

Let $Q = (V_{\rm{tn}} - V_{\rm{tc}}) / ({\rm{tn}} - {\rm{tc}})$, it leads to

By substituting $\Delta x_{\rm{tn}} - \Delta x_{\rm{tc}} = (P_{\rm{tn}} - P_{\rm{tc}})\, \tilde{A} /k$,

Let ${\rm{tn}} - {\rm{tc}} \rightarrow 0$, the above equation can be written as

with the compliance $\tilde{C} = {\tilde{A}^2}/{k}$. The equivalent Windkessel model in this study is shown in Figure A1.