Improving access to greenspaces in the Mediterranean city of Faro

1.   Introduction

Cardiac disease remains one of the tops kills worldwide, in particular the adverse remodelling of cardiac function. Many factors have been acknowledged which are responsible for the deterioration in heart function, of which is stress. A few studies have demonstrated that an imbalanced biomechanical environment can have significant effects on triggering myocardial pathology ^[1,2]. However, it is nearly impossible to measure stress in vivo. On the contrary, patient-specific biomechanical models can predict detailed stress fields, in which the myocardial constitutive law is a crucial component. Many myocardial constitutive models have been proposed to capture myocardial mechanical behaviours, ranging from linear elastic to hyperelastic, from isotropic to anisotropic, and from phenomenological to microstructurally informed constitutive laws ^[3]. In these mathematical models, myofibre structure plays the important role to determine the spatial passive and active stress responses of myocardium.

Nowadays, the prevalent practice is to treat myocardium as an anisotropic and hyper-elastic material. To this end, the invariant-based Holzapfel and Ogden (H-O) model ^[3] has been widely used in the cardiac modelling community for personalized modelling ^[4,5,6,7], which incorporates strain invariants from two families of fibres, one for the myofibre and the other one for the fibre along the sheet direction or the transmural direction. For example, the LivingHeart Project ^[6] used the H-O model in a four-chamber heart model. In a series of studies, Gao et al. ^[5,8] studied myocardial biomechanics both in passive diastole and active systole by implementing this H-O model into an immersed-boundary based finite-element LV model, and later in a poroelastic heart model ^[9]. Wang et al. ^[10] explored the effects of myofibre orientation on the diastolic filling process of the left ventricle (LV) by using the H-O model. Recently, Guan et al. ^[11] studied how accurate the general H-O model is when fitting to various ex-vivo experimental data, such as biaxial tests or simple shear tests. They found that the H-O model has good descriptive and predictive capability for characterizing myocardial mechanical behaviours. Alternatively, Gao et al. ^[12] firstly explored the inverse estimation of material parameters of the H-O model using in vivo data, with recent extension using machine-learning based statistical emulators for faster parameter inference ^[13,14].

Often existing studies assumed myofibres or collagen fibres align perfectly along one unique direction at a specific location. Experimental data has clearly suggested that both myofibres and collagen fibres are dispersed in the myocardium ^[15,16]. With the fast development of imaging technologies, detailed data of collagen network can be measured and quantified, which has led to micro-structurally informed constitutive modelling by taking into account dispersed fibres ^{[17,18,19,20]}. In specific, cardiac modelling studies also ^{[19,21,22,23,24]} have begun to consider dispersed myofibres rather than assuming all fibres aligning perfectly along one direction. To model fibre dispersion, a probability density function is usually used, such as the $\pi$-periodic von Mises distribution ^[18,25,26]. Then, the total strain energy function is the sum of each dispersed fibre's mechanical contribution, such as the angular integration method ^[17]. In general, accounting for each fibre contribution can be very computationally expensive. Another approach to incorporate dispersed fibre contributions is the generalised structure tensor method that analytically determines the proportion of fibre dispersion along each material axis ^[18], while to exclude compressed fibres can be very difficult. To overcome the computational cost in the angular integration approach and the difficulty in excluding compressed fibres in the generalised structure tensor model, Li et al. ^[26] proposed a discrete fibre dispersion model (DFD) that is an approximation of the angular integration model using pseudo-fibre bundles. In this study, we will follow the DFD approach for taking into account fibre dispersion at the same time for the convenience of excluding compressed fibres.

Although a few studies have shown that myofibre dispersion can have significant effects on myocardial mechanics ^{[19,22,24,27]}, a few studies have investigated the effects of sheet dispersion, i.e., the collagen fibre dispersion around the sheet direction, except Eriksson et al. ^[19] who included myofibre and sheet dispersion for the two anisotropic terms in the H-O model. However, they only considered the fully dispersed fibre dispersion by using the $\kappa$-model based on the generalised structure tensor method, and assumed the coupling term was not affected by the fibre dispersion. Full dispersion indicates the dispersed fibres rotationally symmetrically around the mean fibre axis ^[28,29], in other words, fully dispersed in the cross-section plane. Moreover, the data of myofibre and sheet dispersion were cited from different species, and compressed fibre exclusion was not considered in their study.

In this study, we will focus on how fibre dispersion affects myocardial passive behaviours, in particular the sheet dispersion. We firstly extend the DFD-based dispersed model to the sheet dispersion for the H-O model. We then calibrate the H-O model with/without considering the sheet dispersion using the simple shear data from Sommer et al. ^[15] to estimate material parameters. A human LV model in diastole is further simulated to quantify how sheet dispersion affects the LV passive filling. Finally, the effects of myofibre and sheet rotation angle are studied together with dispersed myofibres and sheet fibres.

2.   Methods

2.1.   Invariant-based HO models

To describe the mechanical properties of the myocardium, the invariant-based strain energy function proposed by Holzapfel and Ogden ^[3] is used in this study,

in which $\Psi_{\rm{iso}}$ accounts for the isotropic ground matrix, $\Psi_{\rm{aniso}}$ describes the anisotropic behaviours associated with the two families of fibres and their interactions,

where $a_{({\rm{g, f, s, fs}})}, b_{({\rm{g, f, s, fs}})}$ are material parameters, $I_1 = {\rm{trace}}({\bf{C}})$ is the first invariant of the Cauchy-Green deformation tensor ${\bf{C}} = {\bf{F}}^{\rm{T}} {\bf{F}}$ with ${\bf{F}}$ the deformation gradient tensor, $I_{\rm{4f}} = {\bf{f}}_0 \cdot ({\bf{C}}\, {\bf{f}}_0)$ and $I_{\rm{4s}} = {\bf{s}}_0 \cdot ({\bf{C}}\, {\bf{s}}_0)$ are the invariants representing squared stretches along each fibre direction, in which ${\bf{f}}_0$ is the mean myofibre direction at the reference configuration and ${\bf{s}}_0$ is the sheet direction, $I_{\rm{8fs}} = {\bf{f}}_0 \cdot ({\bf{C}}\, {\bf{s}}_0)$ is the coupling effect between the two families of fibres, and $\mathcal{H}(\cdot)$ is the Heaviside function to ensure the only stretched fibres can bear the load. For example, if ${\bf{f}}_0$ is stretched, then $\mathcal{H}(I_{\rm{4f}}-1) = 1$ with $I_{\rm{4f}}-1 > 0$, and the contribution of ${\bf{f}}_0$ is included into the total strain energy, otherwise $\mathcal{H}(I_{\rm{4f}}-1) = 0$.

2.2.   Discrete fibre dispersion

Studies have found that fibres are spatially dispersed around the mean directions ^[15,16], and the mean myofibre and sheet directions usually form a local material coordinate system (${\bf{f}}_0, \, {\bf{s}}_0, \, {\bf{n}}_0$) at each material point with ${\bf{n}}_0 = {\bf{s}}_0\times {\bf{f}}_0$. Thus a corresponding local spherical polar coordinate system can be defined as Figure 1(a), and a single myofibre, ${\bf{f}}_n(\Theta, \Phi)$, can be defined by the two spherical polar angles $\Theta$ and $\Phi$ according to ${\bf{f}}_0, \, {\bf{s}}_0, \, {\bf{n}}_0$, that is

Given that two fibres lying in one straight line have identical mechanical response, thus the domain of all myofibres in related to the mean direction ${\bf{f}}_0$ can be reduced to be a unit hemisphere with $\mathbb{S}_{\rm{f}} = \left\{(\Theta, \Phi) | \Theta \, \in \left[0, \pi/2\right], \Phi\, \in\left[0, 2\pi \right] \right\}$. In a similar way, a dispersed sheet fibre with respect to ${\bf{s}}_0$ can be defined as

in which $\theta$ and $\phi$ are the two polar angles as shown in Figure 1(a), and its domain is $\mathbb{S}_{\rm{s}} = \left\{(\theta, \phi) | \theta \, \in \left[0, \pi/2\right], \phi\, \in\left[0, 2\pi \right] \right\}$ with respect to the ${\bf{s}}_0$. Please note the myofibre dispersion and sheet dispersion are independent, and we do not consider crosslinks between myofibres and sheet-fibres.

We assume that dispersed fibres at any location can be described by a probability density function $\varrho(\Theta, \Phi)$, and further assume $\varrho (\Theta, \Phi)$ to be composed by two independent functions in terms of $\Theta$ and $\Phi$ following ^[24,27,30]. Specifically, for dispersed myofibres around ${\bf{f}}_0$, we have

in which $\rho_{\rm{in}}(\Theta, b_1)$ denotes the in-plane myofibre dispersion and $\rho_{\rm{op}}(\Phi, b_2)$ describes the out-of-plane myofibre dispersion, $b_1$ and $b_2$ are the concentration parameters, and $G$ is a constant to ensure

The $\pi$-periodic von Mises distribution ^[30] is then used for $\rho_{\rm{in}}$ and $\rho_{\rm{op}}$,

where $\tau$ is the dispersed myofibre angle, $\eta > 0$ is the concentration parameter which can be estimated from measured in-plane and out-of-plane myofibre distributions ^[15,31], and $I_0(\eta) = \frac{1}{\pi} \int_0^{\pi} \exp(\eta \cos(x)) d x$ is the modified Bessel function of the first kind of order zero. Note larger $\eta$ value suggests less dispersion. Figure 2 shows the density function $\rho$ defined in Eq (2.7) with $\eta = 4.5$, $\eta = 3.9$ and $\eta = 0.0$.

Using Eq (2.5), the fibres related strain energy functions by taking into account their dispersion are

where $I_{\rm{4f}}(\Theta, \Phi) = {\bf{f}}_n \cdot ({\bf{C}}\, {\bf{f}}_n)$ and $I_{\rm{4s}}(\theta, \phi) = {\bf{s}}_m \cdot ({\bf{C}}\, {\bf{s}}_m)$.

To exclude non-stretched fibres efficiently in Eq (2.8), Li et al. ^[26] proposed the DFD method by dividing the surface of a hemisphere space domain $\mathbb{S}$ into $N$ spherical triangle elements with representative fibre bundles at each triangle element, see Figure 1(b). Note the DFD approach can also be considered as a quadrature formula for evaluating the integrals in Eq (2.8) but based on spherical triangular elements. In brief, the fibre in the centroid of the $n^{\rm{th}}$ spherical triangle element is defined by a representative fibre ${\bf{f}}_n(\Theta_n, \, \Phi_n)$ for the triangular area $\Delta \mathbb{S}_n$. Then, the fibre distribution probability at this triangle is

where $N$ is the number of spherical triangle elements for the unit hemisphere after discretization. Then $\Psi^{I_{\rm{4f}}\dagger}_{\rm{aniso}}$ can be accordingly approximated as

Similarly, for the sheet dispersion, the total strain energy is

Some studies have tried to consider fibre dispersion in the coupling term $I_{\rm{8fs}}$. For example, Melnik et al. ^[21] used the generalised structure tensor approach ^[18] to account for fibre dispersion in the myocardium with two fully dispersed families of fibres. However, the physical meaning of dispersed $I_{\rm{8fs}}$ between two dispersed fibre bundles has not been studied well. Some studies considered it through cross-link fibres, whilst some studies completely ignored the dispersion in $I_{\rm{8fs}}$. Given the lack of detailed experimental data for $I_{\rm{8fs}}$, we also do not consider fibre dispersion on $I_{\rm{8fs}}$, which shall be studied in the future in particular when modelling fibrosis. Finally, the approximated total strain energy function with two dispersed families of fibres is

The importance of convexity of a strain energy function has been studied in ^[32]. Here we will briefly analyse the convexity of the proposed strain energy function in Eq (2.12). Because the convexity of the two terms $\Psi_{\rm{iso}}$ and $\Psi^{I_{\rm{8fs}}}_{\rm{aniso}}$ has been demonstrated in ^[3], we only discuss the convexity of $\Psi^{I_{\rm{4f}}*}_{\rm{aniso}}$ and $\Psi^{I_{\rm{4s}}*}_{\rm{aniso}}$. For each myofibre bundle, $\rho_n$ is a positive constant, thus for the local Cauchy-Green tensor ${\bf{C}}$, we have the following derivatives

with

Because $a_{\rm{f}}$ and $b_{\rm{f}}$ are positive material parameters, when the myofibre bundle (${\bf{f}}_n$) is under stretch, $I_{\rm{4f}}^n > 1$ ensures both $\Psi^{'}_{\rm{f}}(I_{\rm{4f}}^n) > 0$ and $\Psi^{''}_{\rm{f}}(I_{\rm{4f}}^n) > 0$; when the myofibre bundle is under compression, $\mathcal{H}(I_{\rm{4f}}^n-1) = 0$, then $\Psi^{'}_{\rm{f}}(I_{\rm{4f}}^n) = \Psi^{''}_{\rm{f}}(I_{\rm{4f}}^n) = 0$. Therefore, $\sum_{n = 1}^N \Psi^{'}_{\rm{f}}(I_{\rm{4f}}^n) \ge 0$ and $\sum_{n = 1}^N\Psi^{''}_{\rm{f}}(I_{\rm{4f}}^n) \ge 0$. Similarly, $\sum_{m = 1}^M \Psi^{'}_{\rm{s}}(I_{\rm{4s}}^m) \ge 0$ and $\sum_{m = 1}^M\Psi^{''}_{\rm{s}}(I_{\rm{4s}}^m) \ge 0$ for all positive $a_{\rm{s}}$ and $b_{\rm{s}}$. Finally, the convexity of the strain energy function (Eq (2.12)) can be ensured.

Our previous studies of myocardium only considered the dispersion along myofibre ^[24,27] but not in the sheet direction, and the corresponding strain energy function is

In this study, we further investigate how dispersion in the sheet direction affects passive myocardial mechanic behaviours by using Eqs (2.12) and (2.15) in terms of the fitting to the experimental data and the heart dynamics in diastole.

2.3.   Parameter estimation using shearing experimental data

Shearing experimental data is obtained from Sommer's study ^[15], in which six types of shear tests were performed on human myocardial samples. As shown in Figure 3, the specimen was cut from the LV free wall and followed by six different shear modes using this sample. From Eq (2.1), the passive Cauchy stress is

where $p$ is the Lagrange multiplier to enforce incompressibility of the myocardium, and ${\bf{I}}$ is the identity tensor. Thus, the derived total stress using the dispersed strain energy function Eq (2.12) is

where ${\bf{B}} = {\bf{F}}{\bf{F}}^{\rm{T}}$, ${\bf{f}}^*_{n} = {\bf{F}}{\bf{f}}_n$, ${\bf{s}}^*_{m} = {\bf{F}}{\bf{s}}_m$, ${\bf{f}} = {\bf{F}}{\bf{f}}_0$, ${\bf{s}} = {\bf{F}}{\bf{f}}_0$, and

Similarly, the total stress derived from the strain energy function Eq (2.15) which only consider myofibre dispersion is

Similar to our previous study ^[33], we first estimate material parameters in Eqs (2.12) and (2.15) using a non-linear least square minimization function (fmincon from MatLab, MathWorks 2021) with the following loss function

where $K$ is the total number of data points, ${\bf{\Lambda}}$ denotes the set of unknown parameters, the scalar $\sigma_k$ is the model-predicted stress component according to the corresponding experiment, and $\sigma_k^{\rm{exp}}$ is the measured value. In specific, $\sigma$ is ${\sigma}_{21}$ for the shear mode (fs), ${\sigma}_{31}$ for the shear mode (fn), ${\sigma}_{32}$ for the shear mode (sn), ${\sigma}_{12}$ for the shear mode (sf), ${\sigma}_{13}$ for the shear mode (nf), and ${\sigma}_{23}$ for the shear mode (ns). In this study, the range for each parameter is set to be 0.001–60 ^[33]. To further quantify the fitting goodness, the relative and absolute errors (${\rm{err}}^{\rm{Relative}}$ and ${\rm{err}}^{\rm{Absolute}}$) between the experimental and model-predicted stress-shear curves are introduced,

where $\Delta \gamma_k = |\gamma_k - \gamma_{k-1}|$ is the step size of measured shear amount. In other words, $\sum_{k = 1}^{K-1} \Delta \gamma_k \sigma_k^{\rm{exp}}$ approximates the area under the experimental stress-strain curves, and $\sum_{k = 1}^{K-1} \Delta \gamma_k |\sigma_k({\bf{\Lambda}}) - \sigma_k^{\rm{exp}}|$ is the area enclosed by the measured and model-predicted stress-strain curves. The closer the error value to zero, the more accurate the fitting to the experimental data.

Due to the lack of experimental data for sheet dispersion, we first assume that the dispersion distributions along ${\bf{f}}_0$ and ${\bf{s}}_0$ are same, i.e., $b_1 = b_3$ and $b_2 = b_4$, denoted as Case 1, in which the values of $b_1 = 4.5$ and $b_2 = 3.9$ are adopted from the study by Sommer et al. ^[15]. Based on Case 1, two special cases are considered, which are Case 2 with the same in/out-of-plane dispersion along both the myofibre and sheet direction, and Case 3 with fully-dispersed fibre in the out-of-plane, the so-called fully dispersed case. Since most existing studies only consider myofibre dispersion, including ours, thus we further include three cases based on Cases 1–3 by only considering myofibre dispersion. All simulated cases are summarized in Table 1, and corresponding illustrations of fibre dispersion distributions are shown in Figure 4. Case 1 is considered to be the most realistic one with experimentally measured in/out-of-plane dispersion; Case 4 is the up-to-date model with measured myofibre dispersion ^[15,24]; while Case 6 is the simplest one yet prevalent dispersion model for myofiber and other soft tissue with a fully dispersed distribution ^[19,21,25] which can be considered to be generalised from Case 3.

2.4.   In vivo left ventricle modelling

A subject-specific human LV model from our previous study ^[5] is used here to study passive diastolic filling process ^[10] as shown in Figure 5(a). The LV model has a rule-based myofibre structure with linearly varied myofibre rotation angle from the epicardium ($\bar{\Theta}^{\rm{epi}}$) to the endocardium ($\bar{\Theta}^{\rm{endo}}$). Considering the average wall thickness of the LV model is 8.7 mm and the mean myofibre rotation angle is $14.8^\circ/{\rm{mm}}$ as measured in the human myocardium ^[15], we thus set $\bar{\Theta}^{\rm{endo}} = -\bar{\Theta}^{\rm{epi}} = 60^\circ$ as shown in Figure 5(b). The same myofibre rotation angle has been widely used in the literature ^[4,5,10]. The sheet fibre is along the transmural direction (${\bf{s}}_0$) from endocardium to epicardium at each material point, in other words, the sheet angle is zero.

Following ^[10], the myofibre rotation angle at ventricular wall thickness $\bar{e}$ is

where $\bar{e} = 0$ at endocardial surface and $\bar{e} = 1$ at epicardial surface. The circumferential direction $č_0$ and the orthogonal direction $\tilde{{\bf{n}}}_0$ at each material point can be determined by sheet ${\bf{s}}_0$ and longitudinal direction ${\bf{l}}_0$, which are

Then, the myofibre direction ${\bf{f}}_0$ is defined by

Corresponding sheet-normal direction is ${\bf{n}}_0 = {\bf{s}}_0 \times {\bf{f}}_0$. Similarly, if we would like to include sheet rotation, then the rotated sheet direction ${\bf{s}}^{'}_0$ is

where $\bar{\theta}$ is the sheet rotation angle with respect to ${\bf{s}}_0$ in the ${\bf{s}}_0-{\bf{n}}_0$ plane, and the sheet rotation angle at the thickness of $\bar{e}$ is

in which $\bar{\theta}^{\rm{endo}}$ and $\bar{\theta}^{\rm{epi}}$ are sheet angle at the endocardial and epicardial surfaces, respectively. Note when $\bar{\theta}^{\rm{endo}} = \bar{\theta}^{\rm{epi}} = 0$, all sheets along transmural direction with ${\bf{s}}^{'}_0 = {\bf{s}}_0$.

The LV passive diastolic filling is simulated using ABAQUS 2019 (Dassault Systemes, Johnston RI, USA), and the LV model is loaded with 8 mmHg within $0.5\, s$ with a zero-displacement constraint on the longitudinal movement of nodes on the top base surface. In this study, we consider LV passive mechanics to be quasi-static, and the system of equations to be solved are

where $\Omega$ is the computational domain occupied by the LV geometry, ${P}$ is the LV cavity pressure, ${\bf{n}}$ is the unit normal direction on the endocardial surface, and $u_z$ is the zero-displacement Dirichlet boundary condition on the basal plane.

To measure the LV filling process, we introduce two ratios, they are 1) the radial expansion ratio ($R_r$) that is estimated using the internal diameter ($d$) measured by the two points at the base plane as indicated in Figure 5(a), and 2) the longitudinal elongation ratio ($L_r$) that is calculated using the distance ($l$) between the central point at the base plane and the endocardial apex point, see Figure 5(a). The definitions of $R_r$ and $L_r$ are

where $d_0$ and $l_0$ are the initial values at $t = 0$, and $d_t$ and $l_t$ are the values at time $t$.

3.   Results

3.1.   Fitting to the experimental data

We fit all six cases to the simple shear experiments from Sommer et al. ^[15]. Note Cases 1–3 have dispersion along both ${\bf{f}}_0$ and ${\bf{s}}_0$ using the strain energy function Eq (2.12), while Cases 4–6 only have dispersion along ${\bf{f}}_0$ using the strain energy function Eq (2.15). From our previous study ^[24], we found that when $N \ge 640$, the integration of Eq (2.10) converged to the analytical solution with errors below 5e-5, thus we set $N = 640$ when accounting for dispersed fibre bundle contributions. Figure 6 shows the final fitting results for all six cases. In general, all cases can well describe the mechanical behaviours of the six shear modes. The relative and absolute errors are summarized in Table 2. Case 3 has the least relative and absolute errors with fully dispersed fibre dispersion along both ${\bf{f}}_0$ and ${\bf{s}}_0$, then followed by Case 6 with full dispersion along ${\bf{f}}_0$ only. The errors for other cases are much higher than Cases 3 and 6. This comparison would suggest that different myofibre dispersion can potentially have large effects on myocardial passive response, and the full dispersion could be a good approximation of fibre dispersion if out-of-plane measurement is unavailable. Moreover, incorporating fibre dispersion along ${\bf{s}}_0$ can marginally improve the descriptive capability of a constitutive law to the experiential data, i.e., the H-O model studied here, which may further indicate the necessary to include sheet dispersion if a high-fidelity myocardial model is needed. Inferred material parameters are listed in Table 3. It can be found parameters from one case can vary from other cases, which suggests that shear experimental data alone may not be sufficient to uniquely determine the 8 parameters in the H-O type strain energy function due to parameter correlation, a common and not-resolved issue in personalized cardiac modelling ^[12,33].

3.2.   The human LV model

The LV model in diastole has been simulated with the six cases using the estimated material parameters from Table 3. In order to reduce the computational time, we further compare $N = 40$ and $N = 640$ for discretizing the unit hemisphere of the myofibre dispersion model, nearly identical results are obtained, while the computational time is much reduced for $N = 40$ in a Windows workstation (CPU E5-2680 v3@2.50 GHz and 64.0 GB memory). Thus we set $N = 40$ for all six cases when simulating the LV model in diastole. Note $N = 40$ has also been used in ^[24,26].

Distributions of stress components along myofibre (${\bf{ \pmb{\mathsf{ σ}}}}_{\rm{ff}}$) and sheet (${\bf{ \pmb{\mathsf{ σ}}}}_{\rm{ss}}$) directions at end of diastole in Case 1 are shown in Figure 7. Peak ${\bf{ \pmb{\mathsf{ σ}}}}_{\rm{ff}}$ mainly occurs at the endocardium surface near the LV base, and value of ${\bf{ \pmb{\mathsf{ σ}}}}_{\rm{ff}}$ gradually decreases from the endocardium to the epicardium. Most of sheet fibres are in compressed state with negative ${\bf{ \pmb{\mathsf{ σ}}}}_{\rm{ss}}$ values. Differences of ${\bf{ \pmb{\mathsf{ σ}}}}_{\rm{ff}}$ at each material point between Case 1 and other cases are calculated by $\delta{\bf{ \pmb{\mathsf{ σ}}}}_{\rm{ff}} = {\bf{ \pmb{\mathsf{ σ}}}}^{{\rm{Case}}\, i}_{\rm{ff}}-{\bf{ \pmb{\mathsf{ σ}}}}^{{\rm{Case}}\, 1}_{\rm{ff}}$ with $i\in\{2, 3, \cdots, 6\}$, and their distributions are also shown in Figure 7. These absolute differences are minor with peak value of 0.5 kPa, which could be explained by the optimized material parameters from the same set of experimental data, while the peak relative error can reach 20% with respect to the mean ${\bf{ \pmb{\mathsf{ σ}}}}_{\rm{ff}}$ in Case 1. Cases 2, 4, and 5 with non-full out-of-plane dispersion show almost same stress responses as Case 1, while the fully dispersed Cases 3 and 6 have smaller ${\bf{ \pmb{\mathsf{ σ}}}}_{\rm{ff}}$ and greater ${\bf{ \pmb{\mathsf{ σ}}}}_{\rm{ss}}$ compared to Case 1. Their mean values and standard deviations are summarized in Table 4. Compared to other cases, full dispersion in Cases 3 and 6 leads to slightly larger end-diastolic volume and with higher radial expansion but reduced longitudinal elongation.

Reducing sheet dispersion in the constitutive model contributes to improving computing efficiency, such as saving 42.3% time in Case 4 compared to Case 1. The six cases are further validated by comparing the end-diastolic pressure-volume relationship generated by our models to the measurements from the human heart ^[34]. Normalised end diastolic volume is computed by $({\rm{EDV}}-{\rm{EDV}}_{0})/({\rm{EDV}}_{30}-{\rm{EDV}}_{0})$, in which ${\rm{EDV}}_{0} = 50.3$ ml is the unloading volume, and ${\rm{EDV}}_{30}$ is the LV cavity volume with the diastolic pressure 30 mmHg. Figure 8 compares the model predictions and the human experimental data. It can be found that our simulated pressure-volume curves are overlapped with good agreements with experimental data, which suggests the H-O model with any one of the six fibre dispersion cases can well predict LV cavity volume in diastole.

3.3.   Variation of the mean myofibre and sheet rotation angles

In this section, we further study the effects of mean myofibre and sheet angle variation from the endocardium to the epicardium using Case 1 (with dispersion along both myofibre and the sheet direction). In total, four fibre rotation tests are performed as following

● Test 1: the myofibre angle $\in [-80, 80]$, the sheet angle $\in [0, 0]$.

● Test 2: the myofibre angle $\in [-40, 40]$, the sheet angle $\in [0, 0]$.

● Test 3: the myofibre angle $\in [-60, 60]$, the sheet angle $\in [-30, 30]$.

● Test 4: the myofibre angle $\in [-60, 60]$, the sheet angle $\in [-60, 60]$.

Figure 9 shows the differences of myofibre and sheet stress distributions for the four tests compared to Case 1. Myofibre rotation causes a significant difference in myofibre stress, while sheet stress almost remains constant as in Case 1 when the sheet rotation angle is the same. The mean ${\bf{ \pmb{\mathsf{ σ}}}}_{\rm{ff}}$ of Tests 1 and 2 in Table 5 is smaller than that of Case 1. While compared to Case 1, Test 1 with small myofibre rotation angle has higher ${\bf{ \pmb{\mathsf{ σ}}}}_{\rm{ff}}$ at the endocardium, but Test 2 with large myofibre rotation angle has lower ${\bf{ \pmb{\mathsf{ σ}}}}_{\rm{ff}}$ at the endocardium. The sheet rotation variation has little influence on myofibre stress, while greater sheet stress can be found across the LV wall with increased sheet rotation angle as shown in Figure 9. Table 5 further summarizes EDV, ${\bf{ \pmb{\mathsf{ σ}}}}_{\rm{ss}}$, $R_r$ and $L_r$ for the four test cases. It can be found that the fibre rotation variations can affect the overall LV passive mechanics as reported by Wang et al. ^[10], while the sheet rotation angle mainly affects ${\bf{ \pmb{\mathsf{ σ}}}}_{\rm{ss}}$ with little influence on LV passive filling, i.e., EDV and ${\bf{ \pmb{\mathsf{ σ}}}}_{\rm{ff}}$.

4.   Discussion

In this study, we first implement a DFD-based dispersion model using the H-O model for both the myofibre and sheet directions, and then study the effects of fibre dispersion when fitting to the experimental data and later in a realistic human LV model in diastole. The focus of this study is on the myocardial passive behaviour, and the active contraction will be studied in the accompanied paper. In general, the dispersed H-O models considered in this study can match experimental data ^[15], agree well with Klotz's study ^[34], and produce very similar LV passive dynamics. While when including both myofibre and sheet dispersion with fully dispersed distributions, the best fitting results to the simple shear experimental data can be achieved as shown in Table 2, followed by the case only consider myofibre dispersion with a fully dispersed distribution. However, when the out-of-plane dispersion is not fully dispersed either for myofibres or the sheet fibres, both the relative and absolute errors are increased by twice around. The simulated LV dynamics also have some differences in terms of out-of-plane dispersion. Both Cases 3 and 6 are different from other cases which have non-full out-of-plane dispersion. Therefore, our simulation results demonstrate that myofibre and sheet dispersion can have large effects on myocardial passive response, and the full out-of-plane dispersion could be a good approximation given the very sparse measurements and not-improved fitting to the experimental data.

Angular integration approach and the generalised structure tensor approach have been often used to study fibre dispersion in soft tissue mechanics ^{[18,19,25,26,35]}. It is widely accepted that compressed fibres cannot bear the loading, thus it is necessary to apply such criteria to each dispersed fibre for both myofibres and collagen fibres in the myocardium at every loading time step. Considering the high computing cost of the angular integration approach ^[26] and the extreme complexity of excluding compressed fibres in the generalised structural tensor approach ^[36], we adopt the DFD approach ^[26] in this study, which is a numerical approximation of the angular integration approach with much fewer fibre bundles than the full collagen network, thus it can achieve high computing efficiency ^[26], also observed in our previous study ^[24].

Existed measurements have found that myofibres have both in-plane and out-of-plane dispersion ^[15,31]. As suggested in ^[30], a non-rotational symmetric myofibre dispersion can well describe different in-plane and out-of-plane dispersion. Because of the lack of measured data, a few assumptions have been made in this study, for example, the same dispersion for myofibres and the collagen fibres along the sheet direction in Case 1, and further assumptions of the full distribution for the out-of-plane dispersion. By comparing all six cases, we have quantified the effects of fibre dispersion on myocardial passive behaviours. It can be found that including fibre dispersion can improve the fitting to experimental data, and full out-of-plane dispersion seems a good approximation if there is no measured data of the out-of-plane fibre distribution. Because the material parameters are re-calibrated using the same shear experiments for all cases, the overall behaviours of the LV model are all similar, while differences are still evident as shown in Figure 7. In general, Cases 1, 2, 4, 5 with non-full out-of-plane dispersion have similar stress distribution, but very different from the Cases 3, 6 with full out-of-plane dispersion, the variations being as high as 20% by comparing the peak $\delta \sigma_{\rm{ff}}$ with respect to the mean ${\bf{ \pmb{\mathsf{ σ}}}}_{\rm{ff}}$ in Case 1. Cases 3 and 6 have the highest goodness-of-fit to the experimental data, which would suggest the stress prediction could be more accurate than other cases. Even though Case 1 has measured in-plane and out-of-plane dispersion, presumably the poor fitting result may be subject to the measurement noises in out-of-plane dispersion. Future studies shall include both accurately measured in-plane and out-of-plane dispersion to quantify the effects of out-of-plane dispersion on myocardial mechanics.

Due to the lack of experimental data on the sheet collagen fibre structure, reduced strain energy functions based on the H-O model have been proposed by only including two invariants $I_1$ and $I_{\rm{4f}}$ ^[24,28,37]. We find that such reduced strain energy function can fit uni-axial or bi-axial data but has difficulty in fitting the six different simple shear responses ^[15]. Therefore, the four invariant-based H-O model is used in this study, which also allows to incorporating the sheet dispersion.

According to the microstructural measurements reported by Sommer et al. ^[15], only a few fibres disperse along the sheet direction (${\bf{s}}_0$). Similar results are also shown in Ahmad et al.'s study ^[16], in which the out-of-plane dispersion is much smaller than the in-plane dispersion, suggesting most fibre dispersed in the $({\bf{f}}_0, {\bf{n}}_0)$ plane. For a materiel with preferred fibre direction, i.e., the myocardium, the mean fibre direction is the direction along which the majority of fibres will align that direction. To describe myofibre dispersion, given that more dispersed myofibers in the $({\bf{f}}_0, {\bf{n}}_0)$ plane than those in the $({\bf{f}}_0, {\bf{s}}_0)$ plane, thus the axis ${\bf{f}}_0$ is defined as the first primary direction while ${\bf{n}}_0$ is the second primary direction, see Eq (2.3), and the probability density functions of $\Theta$ and $\Phi$ are both described by Eq (2.7). As mentioned before, more myofibres are dispersed in the $({\bf{f}}_0, {\bf{n}}_0)$ plane, thus it would also suggest that the probability for $\Phi = 0$ is higher than $\Phi = \pi/2$, and they are equal only with full out-of-plane dispersion ($b_2 = 0$).

The myofibre dispersion and sheet dispersion are two independent fibre families in the myocardium and thus are analysed separately. To the authors' best knowledge, no experiments have observed that one dispersed myofibre will have one corresponding orthogonal dispersed sheet fibre. The interaction between the two families of dispersed myofibres is ignored in this study, which means the $I_{\rm{8fs}}$ term remains a phenomenological form. Melnik et al. ^[21] used the generalised structure tensor method to study the effects of myofibre and sheet dispersion, and they found that including fibre dispersion in $I_{\rm{8fs}}$ caused softer material responses than the model without $I_{\rm{8fs}}$ dispersion when using the same set of material parameters. Furthermore, they approved the softening effects caused by fibre dispersion was greater in $I_{\rm{4f}}$ and $I_{\rm{4s}}$ than $I_{\rm{8fs}}$. It is possible to use the discrete fibre dispersion approach to take into account $I_{\rm{8fs}}$ dispersion if assuming that one dispersed myofibre interacts with one corresponding sheet fibres. While the biological explanation of $I_{\rm{8fs}}$ is still unclear, an alternative approach to including dispersion in $I_{\rm{8fs}}$ is to model the cross-link between two families of fibres ^[38], which is beyond the scope of this study.

Finally, we would like to mention limitations. Firstly, only the simple shear data is used to estimate material parameters. The combination of the bi-axial and simple shear data should provide extra information for more accurate parameter inference. However, published data is usually average values of many different samples, leading to the difficulty in matching all data at the same time. Secondly, experimental data of sheet dispersion is lacking, for which we have assumed the sheet dispersion is similar to the myofibre dispersion, and future experiments shall measure that dispersion separately. Thirdly, the mean fibre structure is constructed using a rule-based method without considering spatial heterogeneity. Including more realistic fibre rotation from different regions would be necessary to further improve our understanding of how fibre dispersion affects passive myocardial response. The zero displacement boundary condition on the ventricular base is a simplified implementation. In in vivo, the apex does not move much, instead the basal plane moves up-down. In our LV model, the basal plane can not move along longitudinal direction, while the other regions including the apex are free to move. This is equivalent to fixing the apex and allowing the basal plane free movement, depending on the observer's position either in the apex or in the basal plane. Various studies ^[5,7] have used the fixed basal plane along the longitudinal axis, and also have shown that main features of heart dynamics can be reproduced. A more realistic basal boundary condition may need the measurements of myocardial motion in the basal plane, and also the pericardium needs to be included in order to keep the apex in place ^[39,40].

5.   Conclusions

This study has investigated myofibre and sheet fibre dispersion in passive myocardial mechanics using a widely-used strain energy function, the so-called H-O model. The discrete fibre bundle dispersion model is used to exclude compressed fibres. Our results demonstrate that the H-O model can match ex vivo experimental data very well by including fibre dispersion, in particular when assuming fully dispersed dispersion for the out-of-plane fibre distributions of myofibres and the sheet fibres. Noticeable differences can be found in LV diastolic mechanics when comparing the cases between full and non-full out-of-plane dispersion. Our results seem to suggest that the full out-of-plane dispersion could be a good approximation considering the difficulty in measuring out-of-plane dispersion, and it is necessary to include both dispersion for myofibres and the sheet fibres for the improved descriptive capability to the experimental data and potentially more accurate stress prediction.

Acknowledgements

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

Conflict of interest

The authors declare that there is no conflict of interest.

Data accessibility

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