Modelling of fibre dispersion and its effects on cardiac mechanics from diastole to systole

Detailed fibre architecture plays a crucial role in myocardial mechanics both passively and actively. Strong interest has been attracted over decades in mathematical modelling of fibrous tissue (arterial wall, myocardium, etc.) by taking into account realistic fibre structures, i.e. from perfectly aligned one family of fibres, to two families of fibres, and to dispersed fibres described by probability distribution functions. It is widely accepted that the fibres, i.e. collage, cannot bear the load when compressed, thus it is necessary to exclude compressed fibres when computing the stress in fibrous tissue. In this study, we have focused on mathematical modelling of fibre dispersion in myocardial mechanics, and studied how different fibre dispersions affect cardiac pump function. The fibre dispersion in myocardium is characterized by a non-rotationally symmetric distribution using a π\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi $$\end{document}-periodic Von Mises distribution based on recent experimental studies. In order to exclude compressed fibres for passive response, we adopted the discrete fibre dispersion model for approximating a continuous fibre distribution with finite fibre bundles, and then the general structural tensor was employed for describing dispersed active tension. We first studied the numerical accuracy of the integration of fibre contributions using the discrete fibre dispersion approach, then compared different mechanical responses in a uniaxially stretched myocardial sample with varied fibre dispersions. We finally studied the cardiac pump functions from diastole to systole in two heart models, a rabbit bi-ventricle model and a human left ventricle model. Our results show that the discrete fibre model is preferred for excluding compressed fibres because of its high computational efficiency. Both the diastolic filling and the systolic contraction will be affected by dispersed fibres depending on the in-plane and out-of-plane dispersion degrees, especially in systolic contraction. The in-plane dispersion seems affecting myocardial mechanics more than the out-of-plane dispersion. Despite different effects in the rabbit and human models caused by the fibre dispersion, large differences in pump function exist when fibres are highly dispersed at in-plane and out-of-plane. Our results highlight the necessity of using dispersed fibre models when modelling myocardial mechanics, especially when fibres are largely dispersed under pathological conditions, such as fibrosis.

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 κ model [11] to have κ only depending on the stretched fibre domain, rather than the whole PDF. They further compared this modified κ 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 crossmyofibre 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 κ-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 Fig. 1 Dispersed fibre field. a A unit vector M (red) representing the fibre direction defined by Θ and Φ with respect to the fibre system f 0 , n 0 , and s 0 . b Illustration of a discrete triangular discretization of the unit hemisphere domain centralized with the mean fibre direction f 0 (the red arrow) with N representative fibre directions M q (blue arrows) at the centroid of each triangular surface. (Color figure online) is modelled as in the previous study based on the GST approach. Finally, a comprehensive study of in-plane and outof-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.

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 where a, b, a f , b f are material constants, I 1 = trace(C) and I 4f = f 0 · Cf 0 are strain invariants with C = F T F and F the deformation gradient, 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 f 0 − s 0 − n 0 system, shown in Fig. 1a, thus a single fibre can be described in terms of the two spherical polar angles Θ and Φ in the reference configuration as within the domains of Θ and Φ defined over a unit hemisphere S = {(Θ, Φ)|Θ ∈ [0, π/2] , Φ ∈ [0, 2π ]} as shown in Fig. 1b. Since the two fibres lying in one line are mechanically identical, thus S only needs to be defined over a unit hemisphere. Fibres at any location can be described by a density distribution function ρ(M), further assuming this distribution is essentially separately [40] with respect to Θ and Φ, we have in which ρ op (Φ) describes the out-of-plane dispersion, and ρ in (Θ) describes the in-plane dispersion. Note that ρ M(Θ, Φ) is defined in the reference configuration, its density at current configuration can be derived using F. Following [15], the π -periodic von Mises distribution is adopted here for ρ(Θ) and ρ(Φ), respectively, (4) in which γ is the fibre angle, b > 0 is the concentration parameter, (1/π ) π 0 exp(b cos(x))dx is the modified Bessel function of the first kind of order zero, and G is a constant to ensure Note large concentration parameter b suggests less dispersion, and vice versa. We write Ψ f in Eq. (1) as where 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 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 M q (Θ q , Φ q ) with Θ q and Φ q the spherical coordinates of the centroid of the qth triangle, ΔS q is the triangular area, and the fibre density distribution at this triangle is approximated as in which N is the number of discretized triangles for the unit hemisphere. Thus Ψ * f can be further approximated as where I q 4M = M q · (CM q ). Finally, the passive stress of myocardium with dispersed fibres is where J = det(F), b = FF T , I q * , 1), m q = FM q , p is the hydrostatic-like pressure to ensure incompressibility, and I is the identity matrix. Algorithm 1 lists the detailed steps of numerically calculating σ f in Eq. (9). Algorithm 1 Excluding compressed fibres using the DFD framework for computing passive stress from stretched fibres

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 H, which is also defined over the same unit hemisphere [11,14,15] (10) in which we assume the same ρ(Θ, Φ) as defined in Eq. (3). Because of this symmetry assumption of fibre distributions, H is a diagonal second order tensor as H = diag{H 11 , H 22 , H 33 } and H 11 + H 22 + H 33 = 1 in the local f 0 − s 0 − n 0 system. We further define the myocardial active stress as in which T 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],f = Ff 0 /|Ff 0 |,ŝ = Fs 0 /|Fs 0 |, andn = Fn 0 /|Fn 0 |. None zero H 22 and H 33 represent the cross-fibre contraction. Accordingly, the total myocardial Cauchy stress is 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 4M in the eigenvector space of C: the AI approach In Eq. (8), I q 4M is calculated directly using C and the representative fibre direction M q , an alternative approach to calculate I 4M is to project M into the eigenvector space of C, and then determine whether the fibre bundle M q should be excluded or not. The eigenvector space of the right Cauchy-Green tensor is where λ i is the ith eigenvalue and v i is the corresponding ith eigenvector. Since there is no guarantee to ensure that v 1 , v 2 , and v 3 will form a right orthogonal system, we redefine v 3 as Then a fibre vector M can be rewritten with respect to this eigenvector space where θ, φ are polar angles in {v 1 , v 2 , v 3 } shown in Fig. 2a, and the squared fibre stretch is For incompressible myocardium, when F = I, we have Replacing a = λ 1 − 1, b = λ 2 − 1, c = λ 3 − 1, and x = sin 2 θ , y = sin 2 φ in Eq. (15), we have and it can be further simplified Eq. (18), we now consider the following four scenarios: 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 Ω e is the shaded region defined in Fig. 2b with two sub-regions (Ω e 1 and Ω e 2 ), that is , The corresponding polar angle domains Ω are Fig. 2c, the valid domain for stretched fibres is and the corresponding domain in terms of the polar angles is (4) b = c < 0, a special case of the scenario 3 by setting f (y) = (a/(a − c)), and the domain is Alternatively, the fibre strain energy function per unit volume (Ψ f ) with only stretched fibres can be re-defined with respect to θ and φ, that is where ρ * (θ, φ) = ρ(M(Θ), M(Φ)), and M is a mapping between (Θ, Φ) and (θ, φ) using the following identity; Finally, the total passive stress in terms of θ and φ is where m(θ, φ) = FM(θ, φ). 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.

Algorithm 2 Calculation of Cauchy fibre stress in the eigen-space of C by only counting stretched fibres
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 4M will be same either evaluated at the f 0 − s 0 − n 0 system or the eigenspace of C, and the expression of I 4M according to (Eq. (15)) can be readily obtained, we could replace I q 4M (Θ, Φ) in Algorithm 1 with I q 4M (θ, φ) which is evaluated in the eigen-space of C, thus an updated approach based on Algorithm 1 can be realized (Algorithm 3), see Appendix D for details.

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).   Table 1 Parameters of rabbit myocardium by fitting to the experimental data in [26] and parameters for human LV model adopted from [27] Parameters  While for N = 160, high accuracy of N q=1 ρ 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, N q=1 ρ 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 ≤ 2 and b 2 ≤ 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. Figure 3b is the uniaxial test sample with a dimension of 15 × 5 × 3 mm 3 as used in ex vivo experiments [10]. The fibres lie in the e 1 − e 2 plane and linearly rotate from −60 • to 60 • with respect to e 1 . The strip was virtually stretched along 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. Case 2: b 1 = 0, b 2 = 2, the in-plane isotropic fibre distribution; Case 3: b 1 = 2, b 2 = 0, the transversely isotropic fibre distribution, in other words, the rotationally symmetric distribution [11]; 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 f for the rabbit myocardium (rabbit: 0.097 kPa vs human: 2.4 kPa). Figure 4c and d are the contours of σ 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 → 8), σ 11 becomes much higher suggesting stiffening along e 1 , the mean fibre direction. For example, σ 11 is increased by 10.6% from 3.85 Pa (b 1 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 σ 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.
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   Figure 5 shows the resultant σ 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.

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 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 Δθ = Δφ = 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.

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 Fig. 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, f 0 , s 0 , n 0 are the fibre-sheet-normal system, in which f 0 is the mean fibre direction, s 0 is the sheet direction, in general along the transmural direction from endocardium to epicardium, and 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)  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 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.
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 (≈ −20%) and wall thickening (≈ 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 (σ 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 (≈ −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. 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 → 0 and b 2 → 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 → 0 and b 2 → 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 Fig. 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 σ ff in the middle ventricle enclosed by the black rectangle in Fig. 8, and d the apex twist angles out-of-plane dispersions (b 2 → 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.

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 • ) to the endocardium (60 • ). 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.
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 endsystolic shapes are largely different among the five cases. Similar as in the rabbit model, no contraction happens in  and its fibre architecture (b) generated by a rule-based method with fibre rotation angles from −60 • at epicardium to 60 • 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 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 • ), 15 • 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. 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 ∈ [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 ∈ [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.

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 Fig. 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 σ ff in the middle ventricle indicated by the black rectangle in Fig. 12, and d the apex twist angles  [10] 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 ± 0.008 at end-diastole and J = 1 ± 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 σ a = T a N q=1 ρ qmq ⊗m q wherem q = m q /|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 κ 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, ρ q is constant, thus for the local deformation tensor C, we have and ) > 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 ± 2.2% [55], and −17.75 ± 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 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].

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 π -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 ∈ [0, 2], suggesting the necessary including fibre dispersion in cardiac models when the fibre dispersion is high, especially for pathological myocardium, i.e. fibrosis.