1 Introduction

Computational models have been developed over decades aiming to better understand the mechanism of the heart function under physiological and pathological conditions, ranging from lumped parameter models [1, 2] to single ventricles [3], to bi-ventricles and whole heart [4], and towards multiphysics personalized models for treatment planning and risk stratification [5], etc. One of the essential parts in computational cardiac modelling is the strain energy function, accurate and reliable stress/strain prediction essentially relies on the constitutive modelling of myocardium. The current common practice is to use fibre-reinforced hyperelastic material models by taking into account layered myofibre structures, such as the widely used strain-invariant-based function proposed by Holzapfel and Ogden [6], the so-called HO model.

Micro-structurally informed constitutive modelling in soft tissue has attracted tremendous interest in this area since its introduction in the 1970s [7]. With advanced imaging techniques, such as diffusion tensor magnetic resonance imaging, detailed three-dimensional (3D) fibre distribution for the whole organ, such as the heart, can be acquired in ex/in vivo [8, 9]. Existed measurements have shown that fibres in myocardium are dispersed in space with a predominant mean fibre direction [8, 10]. Constitutive modelling of soft tissue with dispersed fibres have found that dispersed collagen fibres can have a significant effect on the overall mechanical response of the soft tissue [11, 12].

To incorporate fibre dispersion into its constitutive law, one way is to assume a probability distribution functions (PDF) with respect to the mean fibre direction, such as using the \(\pi \)-periodic von Mises distribution [11, 13]. Then the strain energy function for the tissue can be the summarization of each fibre contribution along with other constituents. Broadly speaking, there are two approaches for counting collagen fibre contributions: (1) the angular integration (AI) approach [7], in which the stress from each fibre is added together; and (2) the other one is the generalized structure tensor (GST) approach, which was first proposed by Gasser et al. [11]. By assuming a rotational symmetry for the fibre distribution around the mean direction, Gasser et al. derived a \(\kappa \)-model based on a modified squared mean fibre stretch. Later, Holzapfel and co-workers employed this GST approach to depict mechanical responses of various soft tissues (arterial walls, myocardium, etc.) [12, 14, 15]. Recently, Melnik [16] further extended the GST model to a coupled invariant from two fibre families in the HO model [6].

In soft tissue mechanics, it is often considered that collagen fibres will not bear load when compressed, thus excluding compressed fibres is necessary, or the so-called tension–compression switch [17]. Such exclusion is simple if fibres are not dispersed by simply zeroing out fibre stress if compressed, but can be challenging when fibres are dispersed because a stretched domain needs to be determined at each loading step according to its PDF. With a dispersed fibre distribution, exclusion of compressed fibres in the AI approach is relatively simple by only adding stress contribution from each stretched fibre [18,19,20] compared to the GST approach. For example, Federico et al. [19] excluded compressed fibres from a planar von Mises distribution using a Heaviside function with value to 1 when the fibre is stretched, otherwise zero. In general, the AI approach will require significant computational resource to take into account each fibre’s contribution at each loading step for each location. On the contrary, the GST approach can be much more computational efficient because of one evaluation of the fibre potential derivative with the precomputed structural tensor, one angular integration compared to hundreds of integration in the AI approach [21]. In recent years, a few studies have tried to address the compressed fibre exclusion in GST. When firstly developing the GST model, Gasser et al. [11] suggested to only include fibre contribution when the mean squared fibre stretch is greater than 1. However, as being discussed extensively in the literature [17, 22], even though the mean fibre stretch is less than 1, there is a portion of fibres being stretched depending on the deformation state, thus, the original GST treatment will redistribute the fibre tensile stress over all the fibres. Significant different consequences may happen if the formulation of the squared mean fibre stretched is based on the volumetric/isochoric split as discussed in [23]. To exclude compressed fibres in GST, Melnik et al. [22] introduced a Heaviside function into the integration of the structural tensor which only includes the stretched fibres. In a recent study, Holzapfel and Ogden [13] modified the original \(\kappa \) model [11] to have \(\kappa \) only depending on the stretched fibre domain, rather than the whole PDF. They further compared this modified \(\kappa \) model with a AI model, and concluded that both the GST and the AI models have equivalent predictive power for characterizing various fibre-reinforced soft tissues. In a similar way, Li et al. [24] proposed a general fibre invariant by integrating \((I_4 -1)^2\) only over each stretched fibre.

In both the AI and GST approaches, to exclude compressed fibres under complex dynamics will generally require a two-dimensional integration over a unit sphere at each computational location at each loading step, except for some special cases where analytical solutions may exist. The numerical realization of this two-dimensional integration over a unit sphere may require hundreds of integrations, thus the computational demand can be very high [25]. To improve computational efficiency of this 2-D integration of stretched fibres in finite element simulation, Li et al. [21] developed a discrete fibre dispersion model (DFD). Instead of integrating fibre contributions from thousands of fibres over a unit sphere, the DFD method will firstly divide this unit sphere into finite triangles, each triangle in the sphere surface will associate one uniform fibre bundle which will contribute to the total stress weighted by its corresponding density distribution determined from the corresponding PDF. Li et al. [21] further demonstrated that the computational demand was significantly reduced for excluding compressed fibres in their DFD approach, a speed-up of 224 times was observed in their study than using a traditional AI approach.

It is also necessary to take into account myofibre dispersion in myocardial active stress, which is evidenced in a recent experimental study by Ahmad et al. [8], who reported that myofibres disperse within the ventricular wall. This also agrees with the historical experimental findings from Lin and Yin [26], who measured around 40% cross-fibre active stress in rabbit myocardium. There are three different approaches for modelling active stress in myocardium: the active stress approach [14, 27], the active strain approach [27, 28], and the hybrid approach, or the Hill model, as summarized in [29]. In the active stress approach, the total stress tensor is consisting of passive and active parts, which is convenient when calibrating parameters because of separately measured passive and active stresses [30]. However, coupling two stress tensors from different concepts may have issue on mathematical convexity [27]. In the active strain approach, multiplicative decomposition of total deformation gradient tensor into an elastic part and an activation part is carried out in order to fulfil thermodynamic constraints, which is more inherent to the ‘sliding filament theory’  [28]. In this study, we choose the active stress approach for modelling myocardial stress similar as in our previous work [31] because it is conceptually straightforward to implement by simply counting each individual myocyte’s contribution and less complex for parameter calibration.

There are limited numerical studies which have considered fibre dispersion in active contraction models of myocardium using the active stress approach. For example, Guccione and co-workers explicitly added cross-myofibre active tension into the total active stress tensor, up to 40% in order to match in vivo cardiac function [32,33,34]. Eriksson et al. [14] first incorporated myofibre dispersion in both the passive and active stresses using the GST-based \(\kappa \)-model when modelling a rabbit left ventricle, in which a rotationally symmetric fibre distribution was assumed. They found that the large dispersion in diseased heart can greatly affect ventricular pump function. In a recent study [31], we have studied dispersed active tension in a bi-ventricular model using a non-symmetric fibre distribution derived from an ex vivo DT-MRI dataset, the dispersed active tension was derived by following the GST approach, the results showed that with rule-based myofibre structure, there is a need to employ a dispersed fibre model when simulating cardiac contraction, but fibre dispersion in the passive response of myocardium was not considered.

Published studies on fibre compression-exclusion have been mainly applied to passive responses under simplified loading conditions, i.e. simple extension or simple shear with a rotationally symmetrical fibre distributions or an in-plane dispersion [21, 35]. There is lack of study of how compressed fibre exclusion could affect the ventricular passive filling and further affecting to its pump function overall. Thus there is a need to systematically quantify to which extent the dispersion can affect ventricular dynamics both at diastole and systole under physiological loadings. In this study, we firstly extend the DFD approach [21] into the myocardium with a non-rotationally symmetric fibre distribution, compressed fibres are excluded in computing myocardial passive response, and then the dispersed active tension is modelled as in the previous study based on the GST approach. Finally, a comprehensive study of in-plane and out-of-plane dispersions on ventricular pump functions is carried out, first on a myocardial strip under uniaxial stretch, followed by a rabbit bi-ventricle model, and then a human left ventricle (LV) derived from in vivo clinical data.

2 Constitutive law

2.1 Passive stress

It is a common practice to model myocardium as a hyperelastic fibre-reinforced incompressible material as evidenced in various experimental studies [8, 10] and modelling studies [27, 34, 36]. Here, a reduced form of invariant-based strain energy function for incompressible soft tissue is used by only including the contributions from the ground matrix and the fibres (both collagen and myofibres), that is

$$\begin{aligned} \varPsi =\varPsi _\text {g}+\varPsi _\text {f}=\frac{a}{2b}\{\exp [b( {I}_1-3)]-1\} +\frac{a_\text {f}}{2b_\text {f}} \{\exp [b_\text {f}(\text {max}({I}_\text {4f},1)-1)^2]-1\}, \end{aligned}$$
(1)

where \(a, b, a_\text {f}, b_\text {f}\) are material constants, \( {I}_1 = \text {trace}( \mathbf{C })\) and \( {I}_\text {4f} = \mathbf {f}_0 \cdot \mathbf{C } \mathbf {f}_0\) are strain invariants with \( \mathbf{C }= {\mathbf {F}}^\text {T} {\mathbf {F}}\) and \(\mathbf {F}\) the deformation gradient, \( \mathbf {f}_0 \) is the mean fibre direction in the reference state. The \(\max ( )\) in Eq. (1) is to ensure that only the stretched fibres can bear the loads. Eq. (1) has been widely used for modelling collagenous tissue [20, 37] and myocardium [38, 39].

In general, fibres do not align perfectly along the mean fibre direction, but disperse [8]. To describe such dispersed fibres, we first introduce a spherical polar coordinate system based on the layered fibre structure, the so-called \(\mathbf {f}_0-\mathbf{s} _0-\mathbf{n} _0\) system, shown in Fig. 1a, thus a single fibre can be described in terms of the two spherical polar angles \(\varTheta \) and \(\varPhi \) in the reference configuration as

$$\begin{aligned} \mathbf {M}(\varTheta , \varPhi ) = \cos \varTheta \, \mathbf {f}_0 +\sin \varTheta \, \cos \varPhi \, \mathbf{n} _0 + \sin \varTheta \sin \varPhi \, \mathbf{s} _0, \end{aligned}$$
(2)

within the domains of \(\varTheta \) and \(\varPhi \) defined over a unit hemisphere

\(\mathbb {S}=\left\{ (\varTheta , \varPhi ) | \varTheta \,\in \left[ 0,\pi /2\right] , \varPhi \,\in \left[ 0, 2\pi \right] \right\} \) as shown in Fig. 1b. Since the two fibres lying in one line are mechanically identical, thus \(\mathbb {S}\) only needs to be defined over a unit hemisphere.

Fig. 1
figure 1

Dispersed fibre field. a A unit vector \(\mathbf {M}\) (red) representing the fibre direction defined by \(\varTheta \) and \(\varPhi \) with respect to the fibre system \(\mathbf {f}_0\), \(\mathbf{n} _0\), and \(\mathbf{s} _0\). b Illustration of a discrete triangular discretization of the unit hemisphere domain centralized with the mean fibre direction \(\mathbf {f}_0\) (the red arrow) with N representative fibre directions \(\mathbf {M}_q\) (blue arrows) at the centroid of each triangular surface. (Color figure online)

Full size image

Fibres at any location can be described by a density distribution function \(\rho (\mathbf {M})\), further assuming this distribution is essentially separately [40] with respect to \(\varTheta \) and \(\varPhi \), we have

$$\begin{aligned} \rho \big (\mathbf {M}(\varTheta , \varPhi )\big )= \rho _\text {in}(\varTheta )\,\rho _\text {op}(\varPhi ) \end{aligned}$$
(3)

in which \(\rho _\text {op}(\varPhi ) \) describes the out-of-plane dispersion, and \(\rho _\text {in}(\varTheta )\) describes the in-plane dispersion. Note that \(\rho \big (\mathbf {M}(\varTheta , \varPhi )\big ) \) is defined in the reference configuration, its density at current configuration can be derived using \(\mathbf {F}\). Following [15], the \(\pi \)-periodic von Mises distribution is adopted here for \(\rho (\varTheta )\) and \(\rho (\varPhi )\), respectively,

$$\begin{aligned} \rho (\gamma ,b) = \frac{\exp (b \cos (2\gamma ))}{2G\int _0^{\pi } \exp (b \cos (x)) \mathrm{{d}}x} \end{aligned}$$
(4)

in which \(\gamma \) is the fibre angle, \(b>0\) is the concentration parameter, \(({1}/{\pi }) \int _0^{\pi } \exp (b \cos (x)) \mathrm{{d}}x\) is the modified Bessel function of the first kind of order zero, and G is a constant to ensure

$$\begin{aligned} \int _{\mathbb {S}} \rho (\varTheta ,\,\varPhi )\sin \,\varTheta \, \mathrm{{d}}\varTheta \,\mathrm{{d}}\varPhi =1. \end{aligned}$$
(5)

Note large concentration parameter b suggests less dispersion, and vice versa. We write \(\varPsi _\text {f}\) in Eq. (1) as

$$\begin{aligned} \varPsi ^*_\text {f}=\int _\mathbb {S} \rho (\mathbf {M}(\varTheta , \varPhi )) \varPsi _\text {f}({I}_\text {4M}(\varTheta , \varPhi )) \sin \varTheta \,\mathrm{{d}}\varTheta \,\mathrm{{d}}\varPhi , \end{aligned}$$
(6)

where \(I_{4\mathrm{{M}}}(\varTheta , \varPhi )=\mathbf {M} \cdot \left( \mathbf {C} \mathbf {M} \right) \), and \(\varPsi _\text {f}(I_\text {4f}) = {a_\text{ f }}/({2b_\text{ f }})\{ \exp [b_\text{ f } (\max (I_\text{4f },1) -1)^2 ] -1\}\).

Integrating Eq. (6) analytically can be very challenging because of the \(\max ()\) function, in other words, the exclusion of non-stretched fibres. In a recent study, Li et al. [21] divided the hemisphere space domain \(\mathbb {S}\) into N spherical triangle elements using representative fibre bundles for excluding non-stretched fibres, the so-called discrete fibre dispersion (DFD) approach, shown in Fig. 1b. For detailed description of the DFD approach, please refer to [21]. In brief, the representative fibre direction of a triangle (q) in Fig. 1b is denoted as \(\mathbf {M}_q(\varTheta _q,\,\varPhi _q)\) with \(\varTheta _q\) and \(\varPhi _q\) the spherical coordinates of the centroid of the qth triangle, \(\varDelta \mathbb {S}_q\) is the triangular area, and the fibre density distribution at this triangle is approximated as

$$\begin{aligned} \rho _q=\int _{\varDelta \, \mathbb {S}_q}\,\rho (\varTheta ,\,\varPhi )\,\sin \,\varTheta \, \mathrm{{d}}\varTheta \,\mathrm{{d}}\varPhi , \quad \,q=1,\ldots ,N, \quad \text {subject to} \quad \sum _{q=1}^N \, \rho _q\,=1, \end{aligned}$$
(7)

in which N is the number of discretized triangles for the unit hemisphere. Thus \(\varPsi ^*_\text {f}\) can be further approximated as

$$\begin{aligned} \begin{aligned} \varPsi ^*_\text {f} \approxeq \sum _{q=1}^N\,\rho _q\,\varPsi _\text {f}(I_{4\mathrm{{M}}}^q) \quad \text {with}\quad \varPsi _\text {f}(I_{4\mathrm{{M}}}^q)=\left\{ \begin{array}{ccc}{\varPsi _\text {f}(I_{4\mathrm{{M}}}^q)} &{} { \text{ if } } &{} {I_{4\mathrm{{M}}}^q > 1}, \\ {0} &{} { \text{ if } } &{} {I_{4 M}^q\le 1},\end{array}\right. \end{aligned} \end{aligned}$$
(8)

where \(I_{4\mathrm{{M}}}^q = \mathbf {M}_q \cdot (\mathbf{C} \mathbf {M}_q)\).

Finally, the passive stress of myocardium with dispersed fibres is

$$\begin{aligned} \begin{aligned} \varvec{\sigma }^\text {p}&= J^{-1} a\exp [b( {I}_1-3)] \, {\mathbf {b}} \\&\quad \,+\, \underbrace{ 2 J^{-1}\,\sum _{q=1}^N\,\rho _q\, a_\text {f}\,({I}_\text {4M}^{q*}-1)\,\exp [b_\text {f}({I}_\text {4M}^{q*}-1)^2] \,( {\mathbf {m}}_{q} \otimes {\mathbf {m}}_{q}\,) }_{\varvec{\sigma }_\text {f}} -p\mathbf {I}, \end{aligned} \end{aligned}$$
(9)

where \(J=\text {det}(\mathbf {F})\), \(\mathbf {b}=\mathbf {F}\mathbf {F}^\mathrm{{T}}\), \({I}_\text {4M}^{q*} = \max ({I}_\text {4M}^{q}, 1)\), \(\mathbf {m}_{q}=\mathbf {F}\mathbf {M}_q\), p is the hydrostatic-like pressure to ensure incompressibility, and \(\mathbf {I}\) is the identity matrix. Algorithm 1 lists the detailed steps of numerically calculating \(\varvec{\sigma }_\text {f}\) in Eq. (9).

figure a

2.2 Active stress

To take into account dispersed active stress due to dispersed myofibres, we followed the GST approach as in our previous study [31] by introducing a structural tensor \(\mathbf{H} \), which is also defined over the same unit hemisphere [11, 14, 15]

$$\begin{aligned} \mathbf{H} = \int _{\mathbb {S}} \rho (\varTheta , \varPhi ) \sin \varTheta \, \mathbf {M} \otimes \mathbf {M} \, \mathrm{{d}}\varTheta \,\mathrm{{d}}\varPhi \end{aligned}$$
(10)

in which we assume the same \(\rho (\varTheta , \varPhi )\) as defined in Eq. (3). Because of this symmetry assumption of fibre distributions, \(\mathbf{H} \) is a diagonal second order tensor as \(\mathbf{H} =\text {diag}\{H_{11},H_{22},H_{33}\}\) and \(H_{11}+H_{22}+H_{33}=1\) in the local \(\mathbf {f}_0-\mathbf{s} _0-\mathbf{n} _0\) system. We further define the myocardial active stress as

$$\begin{aligned} \varvec{\sigma }^\text {a}= T_\text {a} H_\text {11} \hat{\mathbf {f}} \otimes \hat{\mathbf {f}} + T_\text {a} H_\text {22} \hat{\mathbf{n }} \otimes \hat{\mathbf{n }} + T_\text {a} H_\text {33} \hat{\mathbf{s }} \otimes \hat{\mathbf{s }} \end{aligned}$$
(11)

in which \(T_\text {a}\) is the active tension along the myofibre direction, and described by a time-varying elastance model that has been illustrated extensively in the literatures [31, 32, 41], \(\hat{\mathbf {f}} = \mathbf {F}\mathbf {f}_0/|\mathbf {F}\mathbf {f}_0|\), \(\hat{\mathbf{s }} = \mathbf {F}\mathbf{s} _0/|\mathbf {F}\mathbf{s} _0|\), and \(\hat{\mathbf{n }} = \mathbf {F}\mathbf{n} _0/|\mathbf {F}\mathbf{n} _0|\). None zero \(H_\text {22}\) and \(H_\text {33}\) represent the cross-fibre contraction.

Accordingly, the total myocardial Cauchy stress is

$$\begin{aligned} \varvec{\sigma }=\varvec{\sigma }^\text {p}+\varvec{\sigma }^\text {a}. \end{aligned}$$
(12)

Note, the additive framework is employed here to model myocardial active contraction, alternative approaches exist, such as the active strain framework, please refer to [27, 42, 43] for more details on the active strain framework.

2.3 Estimating \(I_{4\mathrm{{M}}}\) in the eigenvector space of \(\mathbf{C} \): the AI approach

In Eq. (8), \(I_{4\mathrm{{M}}}^{q}\) is calculated directly using \(\mathbf{C} \) and the representative fibre direction \(\mathbf {M}_q\), an alternative approach to calculate \(I_{4\mathrm{{M}}}\) is to project \(\mathbf {M}\) into the eigenvector space of \(\mathbf{C} \), and then determine whether the fibre bundle \(\mathbf {M}_q\) should be excluded or not. The eigenvector space of the right Cauchy–Green tensor is

$$\begin{aligned} \mathbf {C}=\mathbf {F}^\mathrm{{T}} \mathbf {F}=\mathbf {U}^{2}=\lambda _{1} \mathbf {v}_{1} \otimes \mathbf {v}_{1}+\lambda _{2} \mathbf {v}_{2} \otimes \mathbf {v}_{2}+\lambda _{3} \mathbf {v}_{3} \otimes \mathbf {v}_{3}, \end{aligned}$$
(13)

where \(\lambda _i\) is the ith eigenvalue and \(\mathbf {v}_i\) is the corresponding ith eigenvector. Since there is no guarantee to ensure that \(\mathbf {v}_1\), \(\mathbf {v}_2\), and \(\mathbf {v}_3\) will form a right orthogonal system, we redefine \(\mathbf {v}_3\) as

$$\begin{aligned} \mathbf {v}_3 = \mathbf {v}_1 \times \mathbf {v}_2. \end{aligned}$$

Then a fibre vector \(\mathbf {M}\) can be rewritten with respect to this eigenvector space (\(\mathbf {v}_1\), \(\mathbf {v}_2\), \(\mathbf {v}_3\)), that is

$$\begin{aligned} \mathbf {M}(\theta , \phi )= \cos \theta \, \mathbf {v}_1 + \sin \theta \, \cos \phi \, \mathbf {v}_2 + \sin \theta \sin \phi \, \mathbf {v}_3, \end{aligned}$$
(14)

where \(\theta ,\,\phi \) are polar angles in \(\{\mathbf {v}_1,\,\mathbf {v}_2,\,\mathbf {v}_3\}\) shown in Fig. 2a, and the squared fibre stretch is

$$\begin{aligned} I_{4\mathrm{{M}}}(\theta , \phi )= \mathbf {M}(\theta , \phi ) \cdot (\mathbf{C} \mathbf {M}(\theta , \phi )) =\sin ^{2} \theta \left( \lambda _{2} \cos ^{2} \phi +\lambda _{3} \sin ^{2} \phi \right) +\lambda _{1} \cos ^{2} \theta . \end{aligned}$$
(15)

For incompressible myocardium, when \(\mathbf {F}\ne \mathbf{I} \), we have

$$\begin{aligned} {\text {det}} \mathbf {C}=\lambda _{1} \lambda _{2} \lambda _{3} \equiv 1\quad \Rightarrow \lambda _{1}>1>\lambda _{3}. \end{aligned}$$
(16)

Replacing \(a=\lambda _1-1\), \(b=\lambda _2-1\), \(c=\lambda _3-1\), and \(x=\sin ^2 \theta \), \(y=\sin ^2 \phi \) in Eq. (15), we have

$$\begin{aligned} \begin{aligned} I_{4\mathrm{{M}}}(\theta , \phi )-1&=\sin ^{2} \theta [\cos ^{2} \phi \left( \lambda _{2}-1\right) +\sin ^{2} \phi \left( \lambda _{3}-1\right) ]+\left( 1-\sin ^{2} \theta \right) \left( \lambda _{1}-1\right) \\&=x[(1-y)b+cy ]+a(1-x)=(b-a)x+a-(b-c)xy>0, \end{aligned} \end{aligned}$$
(17)

and it can be further simplified when \(b-c \ne 0 \), we arrive

$$\begin{aligned} {\left\{ \begin{array}{ll} \frac{b-a}{b-c}+\frac{a}{b-c}\,\frac{1}{x}>y, \quad &{}\text{ if } \,\, x> 0 \\ {a > 0 } &{}\text{ if } \,\, x = 0. \end{array}\right. } \end{aligned}$$
(18)

Note \(0 \le x=\sin ^2\,\theta \le 1\) and \(0 \le y=\sin ^2\,\phi \le 1\), \(b-a=\lambda _2 - \lambda _1 \le 0 \) and \(b-c=\lambda _2 - \lambda _3\ge 0\). Thus, from Eq. (18), we now consider the following four scenarios:

  1. (1)

    \(x=0\), \(a=( \lambda _1-1)>0\) is always satisfied, see Eq. (16)

  2. (2)

    \(b\ge 0 > c\), that is \(\lambda _2 \ge 1 > \lambda _3 \), the \(x-y\) curve defined in Eq. (18) is illustrated in Fig. 2b. The \(x-y\) will stay above the x axis except when \(b=0\) in which the \(x-y\) curve crosses the x axis at (1,0), and no crossing point with the y axis. Therefore, the valid domain for stretched fibres \(\varOmega ^e\) is the shaded region defined in Fig. 2b with two sub-regions (\(\varOmega ^e_1\) and \(\varOmega ^e_2\)), that is

    $$\begin{aligned} \varOmega ^e_1=\left\{ \begin{array}{l} {y \in \left[ \frac{b}{b-c}, 1\right] } \\ {x \in \left[ 0,\,f(y)\right] } \end{array}\right. , \quad \varOmega ^e_2=\left\{ \begin{array}{l} {y \in \left[ 0,\,\frac{b}{b-c} \right] } \\ {x \in \left[ 0,\,1\right] } \end{array}\right. ,\quad \text {and} \quad \varOmega ^e=\varOmega ^e_1\,\cup \, \varOmega ^e_2, \end{aligned}$$
    (19)

    where \(f(y)=a/[(b-c)y+a-b]\). The corresponding polar angle domains \(\varOmega \) are

    $$\begin{aligned}&\varOmega _1=\left\{ \begin{array}{l} {{\phi } \in \left[ \eta , \,\pi -\eta \right] \,\cup \,\left[ \pi +\eta , \,2\pi -\eta \right] } \quad \text {with}\quad \eta =\arcsin \sqrt{\frac{b}{b-c}}, \\ {\theta \in \left[ 0,\, \arcsin \sqrt{\frac{a}{(b-c){\sin ^2 \phi }+a-b}} \,\right] }, \end{array}\right. \nonumber \\&\varOmega _2=\left\{ \begin{array}{l} {{\phi } \in \left[ 0, \,\eta \right] \cup \left[ \pi -\eta ,\,\pi +\eta \right] \cup \left[ 2\pi -\eta ,\,2\pi \right] }, \\ {{\theta } \in \left[ 0, \frac{\pi }{2} \,\right] }, \end{array}\right. \quad \text {and} \quad \varOmega =\varOmega _1\,\cup \,\varOmega _2; \end{aligned}$$
    (20)
  3. (3)

    \(0> b > c\), the \(x-y\) curve crosses the x axis at \(({a}/(a-b)\), 0) as shown in Fig. 2c, the valid domain for stretched fibres is

    $$\begin{aligned} \varOmega ^e=\left\{ \begin{array}{l}{y \in \left[ 0, 1 \right] }, \\ {x \in \left[ 0,\, f(y)\right] }, \\ \end{array}\right. \end{aligned}$$
    (21)

    and the corresponding domain in terms of the polar angles is

    $$\begin{aligned} \varOmega =\left\{ \begin{array}{l} {{\phi } \in \left[ 0, 2\pi \,\right] }, \\ {{\theta } \in \left[ 0,\, \arcsin \sqrt{\frac{a}{(b-c){\ \sin ^2 \phi }+a-b}} \,\right] }, \end{array}\right. \end{aligned}$$
    (22)
  4. (4)

    \(b=c<0\), a special case of the scenario 3 by setting \(f(y)= (a/(a-c))\), and the domain is

    $$\begin{aligned} \varOmega =\left\{ \begin{array}{l} {{\phi } \in \left[ 0, 2\pi \,\right] } \\ {{\theta } \in \left[ 0,\, \arcsin \sqrt{\frac{a}{a-b}} \,\right] } \end{array}\right. ; \end{aligned}$$
    (23)
Fig. 2
figure 2

a Fibre directions (\(\mathbf {M}\) ) defined using the eigenvectors (\(\mathbf {v}_1\), \(\mathbf {v}_2\), and \(\mathbf {v}_3\) ) of the right Cauchy–Green tenor \(\mathbf{C} \) with the two fibre angles \(\theta \) and \(\phi \). b the stretched fibres \(\varOmega ^e\) of (x, y) represented by the shaded area for case 2 and c case 3

Full size image

Alternatively, the fibre strain energy function per unit volume (\(\varPsi _\text {f}\)) with only stretched fibres can be re-defined with respect to \(\theta \) and \(\phi \), that is

$$\begin{aligned} \varPsi _\text {f}=\int _{\varOmega } \rho ^*(\theta , \phi ) \varPsi _{\mathrm {f}}\left( {I}_\text {4M}(\theta , \phi )\right) \sin \theta \; \mathrm{{d}}\theta \,\mathrm{{d}}\phi , \end{aligned}$$
(24)

where \(\rho ^*(\theta , \phi )=\rho (\mathcal {M}(\varTheta ), \mathcal {M}(\varPhi )) \), and \(\mathcal {M}\) is a mapping between \((\varTheta , \varPhi )\) and \((\theta , \phi )\) using the following identity;

$$\begin{aligned}{}[\mathbf {v}_1 \,\mathbf {v}_2 \,\mathbf {v}_3]\,[\cos \theta \,\, \sin \theta \cos \phi \,\, \sin \theta \sin \phi ]^\mathrm{{T}} = [\mathbf {f}_0 \,\mathbf{n} _0 \,\mathbf{s} _0]\,[\cos \varTheta \,\, \sin \varTheta \cos \varPhi \,\, \sin \varTheta \sin \varPhi ]^\mathrm{{T}}. \end{aligned}$$
(25)

Finally, the total passive stress in terms of \(\theta \) and \(\phi \) is

$$\begin{aligned} \begin{aligned} \varvec{\sigma }^\text {p}&= J^{-1} a\exp [b( {I}_1-3)] \, {\mathbf {b}} \\&\quad +\, 2 J^{-1}\int _{\varOmega } \rho ^*(\theta , \phi ) a_\text {f}\,({I}_{4\mathrm{{M}}}(\theta , \phi )-1)\exp [b_\text {f}({I}_{4\mathrm{{M}}}(\theta , \phi )-1)^2] \\&\quad \times \,\big ( {\mathbf {m}(\theta , \phi )} \otimes {\mathbf {m}(\theta , \phi )}\big )\sin \theta \;\mathrm {d}\theta \,\mathrm {d}\phi - p\mathbf {I}, \end{aligned} \end{aligned}$$
(26)

where \(\mathbf {m}(\theta , \phi )=\mathbf {F}\mathbf {M}(\theta , \phi )\). Algorithm 2 illustrates the numerical evaluation of the dual integration in Eq. (26), identical to the AI approach. An example of theoretical solution for a special case is explained in Appendix C.

figure b

As has been reported by Li et al. [21] and others, the dual integration in (Eq. (26)) can be very computationally demanding. Since the strain invariant \(I_{4\mathrm{{M}}}\) will be same either evaluated at the \(\mathbf {f}_0-\mathbf{s} _0-\mathbf{n} _0\) system or the eigen-space of \(\mathbf{C} \), and the expression of \(I_{4\mathrm{{M}}}\) according to (Eq. (15)) can be readily obtained, we could replace \(I_{4\mathrm{{M}}}^q (\varTheta , \varPhi )\) in Algorithm 1 with \({I_{4\mathrm{{M}}}^{q}}(\theta , \phi )\) which is evaluated in the eigen-space of \(\mathbf{C} \), thus an updated approach based on Algorithm 1 can be realized (Algorithm 3), see Appendix D for details.

3 Results

In this section, we first studied the influences on stress distributions resulted from different fibre dispersions in a multi-element myocardial strip (Sect. 3.1) under uniaxial stretching, and compared the DFD approach (Algorithm 1) with other two algorithms (2 and 3) in Sect. 3.2. We then studied pump function in a dynamic bi-ventricular rabbit heart model (Sect. 3.3) using the DFD approach for passive response and the GST approach for active response, finally in a dynamic human LV model (Sect. 3.4).

3.1 Uniaxial test on multi-element strip

Fig. 3a schematically illustrates different combinations of \(b_1\) and \(b_2\) as in [44]. Different values of N (40, 80, 160, 640) are chosen for integrating Eq. (7), the relative differences of the numerical integrations compared to the analytical value (1.0) are shown in Fig. 4a. The differences with \(N=640\) are almost negligible for all chosen \(b_1\) and \(b_2\). While for \(N=160\), high accuracy of \(\sum _{q=1}^N \, \rho _q\) can be achieved even for the combination of \(b_1=8\) and \(b_2=8\) in which fibres are highly aligned both in-plane and out-of-plane. When \(N=40\), \(\sum _{q=1}^N \, \rho _q\) is less accurate whenever fibres are highly aligned either in-plane or out-of-plane, while good accuracy can be achieved when \(b_1\le 2\) and \(b_2\le 2\). Therefore, in the following studies, unless explicitly stated, \(N=160\) is used when either \(b_1=8\) or \(b_2=8\), otherwise \(N=40\) is chosen for the sake of computational efficiency.

Fig. 3
figure 3

Surface plots of \(\rho (\varTheta , \varPhi )\mathbf {M}(\varTheta , \varPhi )\) for different combinations of \(b_1\) and \(b_2\) \( \in \{0,1,2,4,6,8\}\) (a), and b a uniaxially stretched strip with linearly rotated fibres represented by red arrows

Full size image

Figure 3b is the uniaxial test sample with a dimension of \(15\times 5\times 3\) mm\(^3\) as used in ex vivo experiments [10]. The fibres lie in the \(\mathbf{e} _1 - \mathbf{e} _2\) plane and linearly rotate from \(-60^{\circ }\) to \(60^{\circ }\) with respect to \(\mathbf{e} _1\). The strip was virtually stretched along \(\mathbf{e} _1\) by 1.2 by fixing one end. Two sets of passive parameters were used in this study, namely the rabbit and the human myocardium. Parameter values for the rabbit myocardium were determined by fitting Eq. (1) to the equal-biaxial experiments from [26] using the same optimization procedure in [36]. Detailed results can be found in Appendix A. Parameter values for the human myocardium were adopted from [27] inferred by matching the measured in vivo heart dynamics. All parameter values are listed in Table 1. Note that fibre dispersion was not taken into account when estimating passive parameters for both the rabbit and human myocardium because of lack of experimental data.

Table 1 Parameters of rabbit myocardium by fitting to the experimental data in [26] and parameters for human LV model adopted from [27]
Full size table

Figure 4b shows \({\sigma }_\text {11}\) using \(N=640\) under uniaxial stretching with four different fibre dispersions, they are

  1. Case 1:

    \(b_1=b_2=0\), the isotropic fibre distribution;

  2. Case 2:

    \(b_1=0,b_2=2\), the in-plane isotropic fibre distribution;

  3. Case 3:

    \(b_1=2, b_2=0\), the transversely isotropic fibre distribution, in other words, the rotationally symmetric distribution [11];

  4. Case 4:

    \(b_1=b_2=2\), the general fibre distribution with both in-plane and out-of-plane dispersion;

Significant differences can be found when using the human passive material parameters, with the highest stress in the transversely isotropic fibre distribution, and the lowest in the planar-isotropic distribution. Similar trend can be found for the rabbit myocardium, but with much less differences, which could be explained by much smaller \(a_\text {f}\) for the rabbit myocardium (rabbit: 0.097 kPa vs human: 2.4 kPa). Figure 4c and d are the contours of \({\sigma }_\text {11}\) at the maximum stretch with different \(b_1\) and \(b_2\) for the rabbit and human myocardium, respectively. With more aligned in-plane fibres (\(b_1 \rightarrow 8\)), \({\sigma }_{11}\) becomes much higher suggesting stiffening along \(\mathbf{e} _1\), the mean fibre direction. For example, \({\sigma }_{11}\) is increased by \(10.6\%\) from 3.85 Pa (\(b_1=b_2=0\)) to 4.26 Pa (\(b_1=2,b_2=0\)), followed by a \(7.5\%\) increment from \(b_1=2,b_2=0\) to \(b_1=8,b_2=0\). Similar trends can be found for the human myocardium, a \(62.9\%\) increase from \(b_1=0,b_2=0\) to \(b_1=2,b_2=0\), and a further \(21.5\%\) increase from \(b_1=2,b_2=0\) to \(b_1=8,b_2=0\). On the contrary, when \(b_1=0\), increasing \(b_2\) seems having little effect on \(\sigma _{11}\) as shown in Fig. 4c and d. Thus, a transversely isotropic fibre distribution will reinforce the material stiffness along the mean fibre direction overall, but not a planar-isotropic distribution.

Fig. 4
figure 4

Relative errors of evaluating Eq. (7) using \(N=40\), \(N=160\), and \(N=640\) with different \(b_1\) and \(b_2\) combinations (a), the relative error is calculated in relate to the analytical value (1.0). b Stress-stretch responses of the uniaxial tension in four dispersion cases by using human (blue lines) and rabbit (red lines) material parameters, respectively. c Contours of stresses at the maximum stretch when using rabbit material parameters and (d) the human material parameters. (Color figure online)

Full size image

A compression experiment using human material parameters was further carried out with the same strip up to 20% shortening along the mean fibre direction. Two dispersions were considered here, one with \(b_1=b_2=1\) and the other one with \(b_1=b_2=2\). Figure 5 shows the resultant \(\sigma _\text {11}\). It can be found that the more dispersed fibre distribution leads to higher resultant stress in magnitude compared to the case with less dispersion. This can be explained by the more stretched dispersed fibres along the cross-fibre direction in the case with \(b_1=b_2=1\) compared to the case with \(b_1=b_2=2\), although the mean fibre direction is under compression. Therefore, dispersed fibres could lead to increased stiffness in compression.

Fig. 5
figure 5

Compression tests with \(b_1=b_2=1\) and \(b_1=b_2=2\), respectively. A maximum of 20% compression is applied along the mean fibre direction, the negative values mean the resultant compressive stress

Full size image

3.2 Comparison between Algorithms 1, 2 and 3

We first compared the computational efficiency between the DFD method (Algorithm 1) and the AI method with the dual integration in the eigenvector system (Algorithm 2), and the hybrid DFD method (Algorithm 3) in a single hexahedron element as shown in Fig. 6a. The mean fibre direction is along the \(\mathbf{e} _1\) axis with a general fibre dispersion (\(b_1=b_2=2\)). Under a maximum stretch of 1.3, all algorithms converge to one identical stretch-stress response, shown in Fig. 6b. Algorithm 1 with \(N=160\) took 21.1 seconds, followed by Algorithm 3 (\(N=160\), 22.3 s), however, Algorithm 2 took about 75.7 s with \(\varDelta \theta =\varDelta \phi =0.0982\) in order to achieve the same accuracy as Algorithms 1 and 2. The much longer computational time used by Algorithm 2 agrees with the findings by Li et al. [21]. All computations were performed using ABAQUS Explicit (Dassault Systemes, Johnston RI, USA) with user-coded subroutine for calculating myocardial stress in a Windows workstation (CPU E5-2680 v3@2.50 GHz and 64.0 GB memory), and only one CPU was used. Therefore, Algorithm 1 has been used in the following studies because of its high computational efficiency.

Fig. 6
figure 6

A single element model with non-rotationally symmetric fibre distribution under uniaxial stretching (a), and b stretch-stress responses using Algorithm 1 (the black line), Algorithm 2 (the red cycles), and Algorithm 3 (the blue squares). (Color figure online)

Full size image

3.3 The dynamic rabbit bi-ventricle model

A three-dimensional bi-ventricle rabbit heart model was constructed from an ex vivo dataset of a healthy rabbit heart, which was acquired using a 7 Tesla Bruker Pharmascan magnetic resonance imaging (MRI) system. The acquired volumetric image had a size of \(160\times 110\times 110\) with voxel dimensions of 0.282 \(\times \) 0.282 \(\times \) 0.282 mm\(^3\). Figure 7a shows one processed diffusion tensor (DT) MRI image. Details of MRI acquisition protocols can be found in the supplementary material. The bi-ventricular wall was first segmented from the ex vivo MRI data using ITK-SNAP Footnote 1 (Fig. 7a), then the boundary contours were exported into SolidWorks (Dassault Systemes, MA USA) for geometry reconstruction, and finally meshed by ICEM (ANSYS, Inc. PA USA). To take into account varied fibre distribution at different regions, the rabbit geometry was further divided into four regions: the left ventricular free-wall (LVFW), the right ventricular free-wall (RVFW), the septum, and the apex, as shown in Fig. 7b. The layered fibre architecture (Fig. 7c) was generated using a rule-based method (RBM) [31, 45], and the fibre rotation angles varied linearly from the epicardium to endocardium according to the measurements from rabbit hearts in [46], and further summarized in Table 2.

Fig. 7
figure 7

The reconstructed bi-ventricle rabbit heart geometry from a DT-MRI data (296,785 linear tetrahedral elements and 55,957 nodes). a An example DT-MRI rabbit heart with delineated ventricular wall enclosed by the red lines, b four regions are defined for the rabbit heart with different colours, \(\mathbf {f}_\mathbf{0 },\,\mathbf{s} _\mathbf{0 },\,\mathbf{n} _\mathbf{0 }\) are the fibre–sheet–normal system, in which \(\mathbf {f}_0\) is the mean fibre direction, \(\mathbf{s} _0\) is the sheet direction, in general along the transmural direction from endocardium to epicardium, and \(\mathbf{n} _0\) is the sheet–normal direction, c rule-based fibre architecture with fibre rotation angles defined in Table 2, d a schematic illustration of the bi-ventricular rabbit model coupled with a circulatory system. MV mitral valve, AV aortic valve, RA right atrium, TV tricuspid valve, PV pulmonary valve, LA left atrium, RA right atrium, AO aorta, Sys systemic circulation, Pul pulmonary circulation, and PA pulmonary artery. Grounded spring with a stiffness of k is tuned to provide the appropriate pressure–volume response (i.e. compliance) for the corresponding cavity. \(C_V\) is the viscous resistance coefficient for describing resistance between cavities. Flow through valves is realized by only allowing uni-directional fluid exchanging between two cavities. (Color figure online)

Full size image
Table 2 Fibre rotation angles at different regions of the rabbit heart adopt from [46]
Full size table

The dynamic bi-ventricular rabbit model was implemented in ABAQUS Explict. In order to simulate diastolic filling and systolic ejection, a lumped model of the pulmonary and systemic circulation systems was attached to this bi-ventricle, which was realized through a combination of surface-based fluid cavities and fluid exchanges [47] as shown in Fig. 7d. After preloading both ventricles with the prescribed end-diastolic pressures (LV: 8 mmHg and RV: 3 mmHg) within 0.1 s, the iso-volumetric contraction starts and the ventricular pressures increase sharply, followed by the systolic ejection at around \(t = 0.15\) s after the valves open when the LV and RV pressures exceed the pressure in the aorta (90 mmHg) and the pulmonary artery (10 mmHg), respectively. The systolic ejection ends at around \(t = 0.3\) s. We further assumed myocardium contracts simultaneously across the whole wall immediately after end-diastole. Further details of this bi-ventricle rabbit model can be found in Appendix B. Cardiac dynamics were simulated with the four representative cases as mentioned in the Sect. 3.1, and an additional case with perfectly aligned fibres (case 5, without dispersion).

Figure 8 depicts the end-diastolic and end-systolic fibre stress distributions for the five representative cases. Both the end-diastolic stress distributions and the deformed shapes are very similar for all cases, but large differences can be found at end-systole. For example, nearly no contraction for case 1 with isotropic fibre distribution (\(b_1=b_2=0\)). Excessive longitudinal shortening (\(\approx -20\%\)) and wall thickening (\(\approx 120\%\)) in case 2 with \(b_1=0,b_2=2\), which has in-plane isotropic fibre distribution. The heart contracts less in case 3 with \(b_1=2, b_2=0\) compared to the general fibre structure in case 4 with \(b_1=2,b_2=2\), which behaves similarly to case 5 (no dispersion) at end-diastole and end-systole. Figure  9a illustrates the pressure–voulme loops (p–v) for all cases, again, different combinations of in-plane and out-of-plane dispersions can significantly affect the pump function of both the LV and RV. Because of no contraction in case 1, its pressure loop degenerates to a point as indicated in Fig. 9a. Case 2 has the highest cardiac work, that is the area enclosed by the p–v loop, in particular in the RV. Cases 4 and 5 have similar cardiac outputs with better performance than case 3. Figure 9b, c and d shows the long-axis shortening, average myofibre stress at mid-ventricle (\(\sigma _\text {ff}\)) and the apical twist for the five cases during systole. For case 1, due to the isotropic fibre distribution, there is little contraction in systole with the smallest long-axis shortening and the apex twist, and nearly zero myofibre stress due to the counteraction of symmetrical fibre distributions. While case 2 is on the other extreme, with the highest long-axis shortening (\(\approx -20\%\)), the lowest myofibre stress due to the isotropic in-plane distribution, and the highest apex twist. For cases 3, 4, and 5, the long-axis shortening, the myofibre stress, and the apex twist follow similar trends as shown in Fig. 9b, c, and d. Considering the largest long-axis shortening and the apex twist, and excessive wall thickening in case 2 compared to experimentally reported values (around 46% in [48], 40% in [28]), it may suggest that the dynamics of case 2 can be unphysiological, especially in systole.

Fig. 8
figure 8

Myofibre stress (\(\varvec{\sigma _\text {ff}}\)) distributions with deformed shapes at end of diastole (top) and end of systole (bottom) in the rabbit bi-ventricle model for case 1 (\(b_1=b_2=0\)), case 2 (\(b_1=0, b_2=2\)), case 3 (\(b_1=2,b_2=0\)), case 4 (\(b_1=b_2=2\)), and case 5 without fibre dispersion

Full size image
Fig. 9
figure 9

Systolic function for the five cases with different dispersed fibres in the rabbit heart model. a Pressure–volume loops, b long-axis shortening, and the long axis is defined as the line segment between the LV basal centre and the apex (the black line in Fig. 8), c average myofibre stress \(\varvec{\sigma _\text {ff}}\) in the middle ventricle enclosed by the black rectangle in Fig. 8, and d the apex twist angles

Full size image

Figures 10a and b show the differences of end-diastolic volume (EDV) of the LV and RV for different combinations of \(b_1\) and \(b_2\) by comparing with case 5 without fibre dispersion. It can be found that the differences in LV EDV and RV EDV are not significant compared to case 5, with a maximum difference of 2% and concentrated close to the isotropic distribution (\(b_1\rightarrow 0\) and \(b_2\rightarrow 0\)). However, significant differences exist for end-systolic volume (ESV) of LV and RV as shown in Fig. 10c and d, in particular when \(b_1\rightarrow 0\) and \(b_2\rightarrow 0\). The solid and dashed lines in Fig. 10c and d indicate the differences within the 5% range compared to case 5. The \(b_1\)–\(b_2\) space within 5% range is much narrower in the LV than RV, indicating that including fibre dispersion in the LV is necessary for this rabbit model. Our results further suggest that fibre dispersion can significantly affect both the LV and RV systolic function for this rabbit heart model, but less in the diastolic filling. With high in-plane dispersion (small \(b_1\)) and less out-of-plane dispersions (\(b_2\rightarrow 8\)), both the LV and RV can contract more, while with high in-plane and out-of-plane dispersions, both the LV and RV pump functions are reduced significantly especially when \(b_1=b_2=0\), despite slightly larger end-diastolic volumes for both the LV and RV chambers.

Fig. 10
figure 10

Relative differences of ESV and EDV values with different dispersion parameters compared to case 5. a EDV differences of the rabbit LV and b RV, c ESV differences of the rabbit LV and d RV, respectively

Full size image

3.4 The human LV model

A human LV model from our previous work is used here [27] to study how different fibre dispersions can affect the pump function in the human LV. Figure 11a shows the LV model with a rule-based fibre structure as illustrated in (b) with linearly varied fibre rotation angle from the epicardium (\(-60^{\circ }\)) to the endocardium (60\(^{\circ }\)). A similar simplified circulation system was attached to the LV model as in the rabbit model, whereas only the aorta and the left atrium (Fig. 11c) were included. To simulate LV dynamics, we first preload the LV to 8 mmHg within 0.3 s, then the iso-volumetric contraction begins, followed by the systolic ejection around at \(t=0.4\) s when the LV pressure exceeds the aortic pressure (75 mmHg), and finally the ejection ends when the LV pressure is lower than the aortic pressure. We also assumed simultaneous contraction in this LV model. The same five cases with varied fibre dispersions are analysed in this section.

Fig. 11
figure 11

A human LV model with 133,042 linear tetrahedral elements and 26,010 nodes (a) and its fibre architecture (b) generated by a rule-based method with fibre rotation angles from \(-60^{\circ }\) at epicardium to \(60^{\circ }\) at endocardium. c A schematic illustration of the human LV model with a circulation system, and the definitions of various symbols are same as those in Fig. 7d

Full size image

The end-diastolic and end-systolic fibre stress distributions for the five cases are shown in Fig. 12. Myofibre stresses at end-diastole are similar across the five cases with slight difference in the deformed shapes. The end-systolic shapes are largely different among the five cases. Similar as in the rabbit model, no contraction happens in case 1, and excessive contraction experienced by case 2. Although cases 3, 4, and 5 have similar end-systolic shapes, the stress distribution is different in case 3 with stress concentration near the base, which may be caused by the less long-axis shortening compared to cases 4 and 5. Figure 13a further summarizes the p-v loops for the five cases. For the cases either with \(b_1=0\) or \(b_2=0\), a larger end-diastolic volume can be achieved, suggesting an increased compliance in diastole because of dispersed fibres. But those dispersed fibres seem compromising active contraction because of the counteracting effects as in case 1 in which all fibres are equally distributed over the unit hemisphere, resulting no contraction at all under incompressible assumption (Fig. 13a). Interestingly, the cardiac output from case 4 with a general fibre dispersion is slightly larger than the non-dispersed case 5, that may be because of the larger end-diastolic volume in case 4, which leads to higher active contraction according to the ‘Frank–Starling’ law. Fig. 13b, c, and d shows the long-axial shortening, the fibre stress, and the apex twist in systole, respectively. Again, an excessive long-axis shortening can be found for \(b_1=0\) and \(b_2=2\) with a peak value of -36.8% (case 2), nearly no shortening in case 1, and the long-axis shortening in case 4 is around 20%, higher than cases 3 and 5. The myofibre stress at mid-ventricle is extremely low in case 1 due to the counteraction of dispersed myofibres, and followed by case 2 with in-plane isotropic fibres. Similar trends can be found in cases 3, 4 and 5 but with different levels. Note because of dispersed myofibres, there are cross-active tension along the sheet and sheet–normal direction, which also affect myocardial contraction. For example, even though the myofibre stresses in cases 2 and 3 are low, because of the in-plane myofibre distribution, large contractile stress can be found along the sheet–normal direction, which will further contribute to contraction as suggested in [31]. The maximum apex twist angles are similar for cases 2, 3, and 5 (around 25\(^{\circ }\)), 15\(^{\circ }\) for case 4, and lowest in case 1 as expected. Thus for this human LV model, a fully dispersed myofibre structure (case 1) is non-physiological because of nearly no contraction, while an in-plane isotropic fibre structure (case 2) can lead to non-physiological dynamics in systole because of the excessive long-axis shortening.

Fig. 12
figure 12

Myofibre stress distributions (\(\varvec{\sigma }_\text { ff }\)) at the end of diastole (top) and at the end of systole (bottom) for the human LV model with different combinations of \(b_1\) and \(b_2\)

Full size image
Fig. 13
figure 13

Systolic function for the five cases in the human LV model. a Pressure–volume loops, b long-axis shortening, the long axis is defined as the link between the LV basal centre and the apex (the black line in Fig. 12), c average myofibre stress \(\varvec{\sigma _\text {ff}}\) in the middle ventricle indicated by the black rectangle in Fig. 12, and d the apex twist angles

Full size image

Figure 14 is the contour plot of the relative EDV and ESV differences with varied \(b_1\) and \(b_2\) compared to the case with non-dispersion fibres (case 5), the superimposed lines indicate the parameter space with differences in the range of ± 5%. The parameter space with \(>5\%\) difference for EDV mainly locates near the \(b_1\) axis and \(b_2\) axis, with a maximum difference of 18% when \(b_1=0\) and \(b_2=0\). The differences in ESV are much more significant as shown in Fig. 14b. Still, the regions near the \(b_1\) and \(b_2\) axis have reduced myocardial contraction with much larger ESV compared to case 5. In Fig. 14b, we can also find a region enclosed by the dashed line which has much less ESV compared to case 5, this can be explained by the large in-plane dispersion which can enhance pump function as found in our previous study [31]. Thus for this human LV model, when (\(b_1, b_2 \in [0,\, 2]\)), there is a need to incorporate fibre dispersion, beyond that, the differences of LV EDV and ESV are in general within 5% range compared to the case without considering fibre dispersion. In fact, the same parameter space of \(b_1, b_2 \in [0,\, 2]\) can be found for the rabbit model, where the large difference (i.e. \(>5\%\) ) exists as shown in Fig. 10. The measured fibre dispersion in human myocardial samples from [10] is represented by the white dot in Fig. 14, which lies in the \(<5\%\) region of the EDV and ESV differences compared to case 5.

Fig. 14
figure 14

Relative differences of EDV (a) and ESV (b) with different dispersion parameters compared to case 5 for the human LV model. The contour lines indicate \(\pm 5\)% difference, and the white dot is the measured dispersion (\(b_1=4.5, b_2=3.9\)) from Sommer et al. [10]

Full size image

4 Discussion

In this study, we have focused on how different fibre dispersions will affect myocardial mechanics both passively and actively. Fibre dispersion is described by a non-rotationally symmetric distribution based on recent experimental studies [8, 10]. In order to exclude compressed fibres for passive response, we adopted the discrete fibre dispersion model developed in [21], and then the general structural tensor [11] was employed for describing dispersed active tension as in our previous study [31]. We first studied the numerical accuracy of the integration of fibre contributions using the DFD approach, then studied the different mechanical response in a uniaxially stretched myocardial sample with varied fibre dispersions. Two heart models were further employed in this study, the rabbit bi-ventricle model and the human LV model. Cardiac dynamics from diastole to systole were simulated with different fibre dispersions. Our results show that fibre dispersion can have significant effects on myocardial mechanics and the pump function, which highlights the necessity of using appropriate dispersion models when modelling myocardial mechanics, especially when fibres are largely dispersed.

The strain energy function for myocardium (Eq. (1)) only includes two strain invariants \(I_1\) and \(I_4\) . We further fitted this strain energy function to the biaxial stretching experimental data of human tissues reported by Sommer et al. [10], see Appendix F. Good agreements were achieved along the fibre direction and the cross-fibre direction. Myocardium was further assumed to be incompressible, a widely adopted assumption in literature [6, 27, 32,33,34]. In our human LV model, \(J=1\pm 0.008\) at end-diastole and \(J=1\pm 0.009\) at end-systole, which suggests myocardial incompressibility was achieved in our simulations, and the very small standard deviation is due to the numerical realization of the incompressible penalty. Interested readers can refer to [49] for a recent contemporary review of constitutive modelling of myocardium.

To model fibre dispersion, the AI approach and the GST approach have been widely used in various soft tissue mechanics, such as arterial wall, myocardium, tendon, valvular tissue, and skin. One striking feature of collagen fibres is that they cannot bear the loading when compressed, thus it is necessary to exclude compressed fibres. In this study, we have adopted the DFD approach developed in [21], which is a numerical approximation of the AI approach. By discretizing a unit hemisphere into N spherical triangles, the DFD approach can then integrate stress contributions from stretched fibre bundles in a much faster way as demonstrated by Li et al. [21], and also in this study. The advantage of the DFD approach is that by replacing the double integration in the AI approach with finite representative fibre bundles, it preserves the straightforward compressed fibre exclusion from the AI approach.

In our previous study [31], we have determined the dispersed active tension from a DT-MRI acquired fibre structure using the GST approach along with a rule-based fibre structure. We found that this GST-based dispersed active stress model can well approximate the systolic function in a bi-ventricle model compared to the DT-MRI derived fibre structure, the most realistic fibre structure. Since there is no need to exclude any myofibre during active contraction, the GST approach is naturally convenient to take into account dispersed active tension for the sake of computational efficiency. In this study, we further tested the DFD approach for the active stress with one hexahedra element, and the total active tension is \(\varvec{\sigma }^\text {a}= T_a\,\sum _{q=1}^N\,\rho _q\,\hat{\mathbf {m}}_q \otimes \hat{\mathbf {m}}_q\) where \(\hat{\mathbf {m}}_q=\mathbf {m}_q/|\mathbf {m}_q|\). Because the active tension is much higher compared to the passive stress part at systole, a much larger N was needed compared to the value used for the passive response during the numerical integration of the fibre contribution. We then compared the DFD-based and GST-based active tension in the human LV model with \(N=160\), the end-systolic volume difference was negligible for the two approaches, while the computational time was much longer for the DFD-based active tension, which took 6 times longer for one cardiac cycle than the GST-based active tension.

According to the studies from Sommer et al. [10] and Ahmad et al. [8], the in-plane and out-of-plane dispersions are largely different, thus a non-rotational symmetry assumption is more appropriate than the rotationally symmetry distribution as assumed in the \(\kappa \) model [14]. Furthermore, the non-rotational symmetric fibre distribution allows us to study how different degrees of in-plane and out-of-plane dispersions can affect ventricular dynamics, for example from the fully isotropic to the in-plane isotropic fibres, the transversely isotropic fibres, the general dispersion, and the highly aligned fibre structures as depicted in Fig. 3a.

In the DFD approach, the dispersed fibre distribution within an hemisphere is first divided into N triangles before excluding compressed fibres, see Fig. 4. In Li et al.’s study [21], they used \(N=640\) for the hemisphere discretization when modelling the simple shear and the uniaxial tests. Because the computational demanding for the rabbit and human heart models can be very high even without fibre dispersion, usually hours for one cardiac cycle, thus we first determined the possible minimum N for integrating fibre contributions with different \(b_1\) and \(b_2\). We found that in general \(N=40\) can achieve a fairly good agreement with the theoretical value if fibres are not highly aligned (Fig. 4a), but except for \(b_1=8\) or \(b_2=8\) in which a larger N will be needed. A further test on the human LV model with \(N=160\) and \(N=40\) showed that the computation time with \(N=160\) was almost 3.2 times longer than with \(N=40\), while the pump function was nearly identical.

Fibre dispersion seems having little influence on the passive filling of the rabbit bi-ventricle, but not for the human heart. This agrees with the our previous work using a neonatal porcine bi-ventricle model [31], in which the LV and RV end-diastolic volume differences were around \(1.4\%\) between a rule-based fibre structure without dispersion and a DT-MRI fibre structure which is naturally dispersed. The size of the neonatal porcine heart in [31] is similar to the rabbit heart in this study. The possible reasons for the different impacts from fibre dispersion between the rabbit and LV models are (1) the much thinner wall in the rabbit heart (4 mm) compared to the human heart (10 mm), thus the changing of mean fibre angle is more rapid in the rabbit heart than the human heart, which may indicate that the diastolic filling in the rabbit heart could be more sensitive to the mean fibre rotation angle; (2) different fibre structures, i.e. different fibre rotation angles; (3) different passive material properties; (4) the much smaller size of the rabbit heart compared to the human heart.

Because of the wavy structure of the collagen network in the soft tissue, collagen fibres are initially crimped, and gradually recruited to bear the loading with increased stretch [50, 51]. Only recently, Cheng et al. [51] assessed collagen fibre recruitment in bladder tissue using advanced bioimaging, and experimentally demonstrated that the low resistance in the toe regime, corresponding to the low stretch regime, can be explained by the no-discernible recruitment of collagen fibres. This will support the fundamental hypothesis in this study, also among others  [21, 52], that is a straight fibre under compression will buckle and not support load because of its crimped configuration. This assumption is also necessary for reasons of stability as discussed in [52]. Including recruitment into the strain energy function would be more physiologically relevant compared to the simple tension–compression switch [53]. Another way to take into account the crimped wavy collagen fibre network is to adopt a multiscale approach from nanoscale up to the macro-scale using homogenization techniques as in [54]. In this study, the tension–compression switch is used because of its simple numerical implementation.

The importance of convexity of a strain energy function has been studied in [52] for ensuing material stability and meaningful mechanical behaviour. Here we will briefly discuss the convexity of the strain energy function in the DFD approach when only considering myocardial passive response, see Eqs. (1) and (7). The convexity of the isotropic part in Eq. (1) has been demonstrated in [6], thus we only discuss the convexity of the anisotropic part (Eq. (7)). For each fibre bundle, \(\rho _q\) is constant, thus for the local deformation tensor \(\mathbf{C} \), we have

$$\begin{aligned} \frac{\partial \varPsi _\text {f}(I_{4\mathrm{{M}}}^q)}{\partial \mathbf{C} }= \varPsi ^{'}_\text {f}(I_{4\mathrm{{M}}}^q) \mathbf {M}_q \otimes \mathbf {M}_q, \qquad \frac{\partial ^2 \varPsi _\text {f}(I_{4\mathrm{{M}}}^q)}{\partial \mathbf{C} \partial \mathbf{C} }= \varPsi ^{''}_\text {f}(I_{4\mathrm{{M}}}^q) \mathbf {M}_q \otimes \mathbf {M}_q \otimes \mathbf {M}_q \otimes \mathbf {M}_q, \end{aligned}$$
(27)

and

$$\begin{aligned} \begin{aligned}&\varPsi ^{'}_\text {f}(I_{4\mathrm{{M}}}^q) = a_\text {f} \,(I_{4\mathrm{{M}}}^q-1)\,\exp [b_\text {f}(I_{4\mathrm{{M}}}^q-1)^2],\\&\varPsi ^{''}_\text {f}(I_{4\mathrm{{M}}}^q)= a_\text {f} \,\exp [b_\text {f}(I_{4\mathrm{{M}}}^q-1)^2]\,[1+2b_\text {f}(I_{4\mathrm{{M}}}^q-1)^2]. \end{aligned} \end{aligned}$$

Because \(a_\text {f}\) and \(b_\text {f}\) are positive parameters, when the fibre bundle (\(\mathbf {M}_q\)) is under stretch, \(I_{4\mathrm{{M}}}^q > 1\) ensures both \(\varPsi ^{'}_\text {f}(I_{4\mathrm{{M}}}^q)>0\) and \(\varPsi ^{''}_\text {f}(I_{4\mathrm{{M}}}^q)>0\); when the fibre bundle is under compression, by setting \(I_{4\mathrm{{M}}}^q = 1\) or equivalently \(\varPsi _\text {f}(I_{4\mathrm{{M}}}^q)=0\), then \(\varPsi ^{'}_\text {f}(I_{4\mathrm{{M}}}^q)=\varPsi ^{''}_\text {f}(I_{4\mathrm{{M}}}^q)=0\). Therefore, \(\sum _{q=1}^N \varPsi ^{'}_\text {f}(I_{4\mathrm{{M}}}^q) > 0\) and \( \sum _{q=1}^N\varPsi ^{''}_\text {f}(I_{4\mathrm{{M}}}^q) > 0\) ensure the convexity of Eq. (7) . During active contraction, the adopted active stress approach may violate the thermodynamic constraints [27], and might further lead to non-convexity and instability issue. As discussed in other studies [27, 28], the active strain approach would be an alternative approach if thermodynamic constraints need to be enforced. Nevertheless, we have not met instability issue in this study.

We further found that active contraction is sensitive to fibre dispersions for both the human and rabbit models. For example, the isotropic fibre distribution leads to almost no contraction, while the in-plane isotropic fibre distribution results purely in-plane active stresses with a high proportion of active tension along the sheet–normal direction, which leads to excessive longitudinal contraction in the human model (−36.8%) than the reported normal range in human (\(-16.7\pm 2.2\%\) [55], and \(-17.75\pm 5.44\%\) [56]). For both the LV and rabbit models, a general dispersed fibre structure, i.e. \(b_1=2\) and \(b_2=2\) can achieve a slightly larger cardiac output than the non-dispersed fibre structure, that is because the cross-fibre active stress can enhance the pump function as demonstrated in our previous study [31] and other studies [26, 34].

We would like to mention limitations. First, the passive strain energy function only incorporates the matrix and fibre contributions, but not including the terms associated with the sheet direction and the sheet–normal as in other studies [6, 36], that is because of lack of experimental data for identifying all parameters and dispersed fibre measurements in those two directions. Secondly, fibre dispersion is only considered along the fibre direction \(\mathbf {f}\), it is possible to include the sheet dispersion and the sheet–normal direction as in [16], but it would make the computation very challenging. Thirdly, a lumped circulation model is used for the ventricular models to provide physiologically accurate pressure boundary conditions. This lumped circulation model is similar to the Windkessel model and realized by a combination of surface-based fluid cavities and uni-directional fluid exchanges [31, 34]. Using a more realistic circulation model, such as one-dimensional models, will allow us to systematically investigate the interactions between ventricles and blood flow in vessels [57, 58]. A further limitation is that the electrophysiology is not modelled in this study, assuming all myocytes contract simultaneously following our previous studies [31, 45, 57] and other studies [34]. Reasons are that in healthy hearts (1) the action propagation in the left ventricle is much faster than the mechanical contraction, (2) an electromechanical model will significantly increase the modelling complexity, such as the Purkinje fibre network [59], and further difficulties in parameter calibration. We refer interested readers to the review [60] for multiphysics modelling. Last but not the least, the experimental data for the rabbits were from different studies, and rule-based fibre structures were assumed for both heart models. A combined experimental measurements (biaxial and shear mechanical tests, DT-MRI fibre acquisition, ventricular pressure measurements, in vivo dynamics imaging, etc.) with the computational modelling from the same heart will be desirable to gain a deeper understanding of how different fibre structures affect ventricular pump function.

Despite those above limitations, the present computational framework can be readily extended to subject-specific multiphysics simulations [60] by including other heart chambers, electrophysiology, ventricular blood flow, and perfusion within myocardium. Future studies shall also explore state-of-the-art assimilation methods in computational cardiology for clinical translation [5], and fast computation using cutting-edge machine learning approaches [49, 61].

5 Conclusion

This study systematically investigates the impact of fibre dispersions on myocardial mechanics both passively and actively, first on a myocardial strip, then a rabbit bi-ventricle model, and finally a human LV model. The fibre dispersion in myocardium is characterized by a non-rotationally symmetric distribution using a \(\pi \)-periodic von Mises distribution. To exclude compressed fibres, two different approaches are compared, including the discrete fibre model and the angular integration-based approach within the eigen-space of the right Cauchy-Green tensor. The dispersed active tension is derived from the general structural tensor approach. Our results show that the discrete fibre model is preferred for excluding compressible fibres because of high computational efficiency as already demonstrated in the literature. Our results further suggest that both diastolic filling and systolic contraction can be largely affected by dispersed fibres depending on the in-plane and out-of-plane dispersion degrees, especially for systolic contraction. The in-plane dispersion seems to affect myocardial mechanics more than the out-of-plane dispersion, an inappropriate dispersed fibre structure will result in a non-physiological dynamics (i.e. in-plane isotropic fibres). Despite different effects in the rabbit and human models caused by the fibre dispersion, large differences in pump function exist when \(b_1, b_2 \in [0, 2]\), suggesting the necessary including fibre dispersion in cardiac models when the fibre dispersion is high, especially for pathological myocardium, i.e. fibrosis.