An Inverse Finite Element Method for Determining the Tissue Compressibility of Human Left Ventricular Wall during the Cardiac Cycle

The determination of the myocardium’s tissue properties is important in constructing functional finite element (FE) models of the human heart. To obtain accurate properties especially for functional modeling of a heart, tissue properties have to be determined in vivo. At present, there are only few in vivo methods that can be applied to characterize the internal myocardium tissue mechanics. This work introduced and evaluated an FE inverse method to determine the myocardial tissue compressibility. Specifically, it combined an inverse FE method with the experimentally-measured left ventricular (LV) internal cavity pressure and volume versus time curves. Results indicated that the FE inverse method showed good correlation between LV repolarization and the variations in the myocardium tissue bulk modulus K (K = 1/compressibility), as well as provided an ability to describe in vivo human myocardium material behavior. The myocardium bulk modulus can be effectively used as a diagnostic tool of the heart ejection fraction. The model developed is proved to be robust and efficient. It offers a new perspective and means to the study of living-myocardium tissue properties, as it shows the variation of the bulk modulus throughout the cardiac cycle.


Introduction
One of the tools that can potentially lead to the early detection of human heart failures is the understanding of the mechanical behavior of the human left ventricle (LV) in normal and diseased states. This leads to continuous interest in the determination of the material properties of myocardium through the mechanical testing of excised strips of myocardium. These strips under prescribed homogeneous loading conditions give stress-strain relationships. Originally, these tests were done uniaxially but more recently, biaxial tests were also performed. Uniaxial test is used to define passive stress-strain relationships in the fiber direction [1]. It is very useful test in determining the general characteristics of the behavior of cardiac tissue in both healthy and diseased states, but is not sufficient to provide a unique description of threedimensional (3D) constitutive behavior of the myocardium. Due to the incompressibility of cardiac tissue, biaxial tests can be used to determine certain multidimensional stress-strain relationships for the fiber and cross-fiber [2]. Despite the frequent assumption that human myocardial tissue is incompressible, the fact remains: all myocardial tissues have some degree of compressibility. Furthermore, the issue of the compressibility of myocardium tissue is raised due to the systolic intra-and extra-vascular blood displacements [3]. As a result of the previously mentioned fact, it is clear that the myocardium tissue compressibility changes during the cardiac cycle.
The information embodied in the myocardial tissue bulk modulus adds further insight to the mechanical nature of the soft tissue. Bulk modulus is very important as a standalone parameter, and as additional information to the shear/Young's modulus. Precise values of the myocardium bulk modulus are especially required to improve the accuracy of finite element (FE) simulation for the modeling of the human heart. In the last several decades, publications related to cardiac modeling have handled the myocardial bulk modulus in many different approaches. The published values of myocardial bulk modulus recorded by researchers who were interested in simulating the performance of the LV during the diastolic phase are quite small, for instance 28 kPa [4] and 160 kPa [5]. This may be attributed to small changes in ventricular wall volume. Some studies evaluated the bulk modulus under the assumption that during the rapid filling phase, the volume change of the ventricular wall should be less than 10%. Additionally, relatively high values of bulk modulus were used by researchers who analyzed the systolic phase, for example 380 kPa [6], 600 kPa [7], and 25 MPa [8]. Using high values for bulk modulus during systolic phase are due to the systolic intra-and extra-vascular blood displacements that give rise to the compressibility of the tissue. Mean constant value of myocardial bulk modulus was assumed during each cardiac phase.
Despite the widespread use of uniaxial and biaxial tests for determining the myocardium characteristics, there are four major problems that arise from these kinds of studies [9,10]: 1. Tests were carried out on non-human tissue, and as a result may not be directly applicable to humans. 2. The properties may not be directly applicable to FE models of the heart due to heterogeneous behavior of the myocardium. 3. The mechanical properties of myocardium changed drastically immediately after death. 4. The variation in the values of the mechanical properties according to the experimental loading conditions.
All the above limitations lead researchers to find different ways to run their experiments without having to excise samples of myocardium. Hence, a group of researchers moved to use Magnetic Resonance Imaging (MRI) and FE or mathematical methods to determine the mechanical properties in vivo [11,12], while others applied the scanning acoustic microscope with high frequency ultrasound to measure the bulk modulus and describe the mechanical properties of the myocardium [13]. However, most bulk modulus experimental values obtained from these studies are very high (<3 GPa), and cannot be used directly to FE modeling. To overcome these shortcomings, the FE approach is suggested to determine, in vivo, the myocardial tissue bulk modulus during the cardiac cycle.
Usually, FE analysis can be done directly once the input parameters such as LV geometry, LV internal cavity pressure and the myocardium tissue properties are known. Once the model is constructed, the desired outputs such as deformation behavior, LV cavity volume, LV wall stresses and strains can be predicted from the model. However, it is not uncommon in reality that some or all of the output values are known from experiments beforehand, while some of the input parameters still need to be determined. This requires the FE analysis to be done in an inverse way, where iteration of FE simulation is performed to find the material properties that give the best fit between the computed and experimentally measured LV internal cavity volumes. An inverse FE approach is a complex engineering process that can determine the unknown causes of known consequences. This approach has the advantage that the determination of the dynamic properties is measured non-invasively [14][15][16][17].
A novel method combining the experimentally measured LV pressure-volume curves and an inverse FE method is proposed to determine myocardial bulk modulus. The main purpose of this research is to develop an inverse FE procedure with ANSYSH computer code for the determination of the bulk modulus of human LV during cardiac cycle. The proposed inverse technique is based on published experimental measured LV pressure-volume curves. By using these outputs of published LV experimental data, the bulk modulus versus time curve is traced through inverse technique. Based on the obtained results, the repeated changes of myocardium tissue bulk modulus in the LV wall during cardiac cycle result in a highly efficient global function of the normal heart. Therefore, the myocardium bulk modulus can be effectively used as a diagnostic tool of the heart ejection fraction.

Left Ventricle Geometry and FE Model
An inverse FE model was adopted to evaluate the LV tissue compressibility during cardiac cycle. 3D-FE model was built to simulate the deformation mechanics of the LV using ANSYSH commercial software. To simplify the analysis, FE simulation model was represented by an ellipsoid, truncated at two-thirds of the major axis including two sets of fibers (myocardial fibers bound by a mesh of collagen fibers), which were attached to each other to form a spatial network. The geometric parameters and dimensions of the LV model in the initial undeformed configuration (at hypothetical zero pressure applied inside the LV cavity) are shown in Figure 1. The wall thickness of the LV model, in the reference unstressed state, was divided into seven equal thickness layers. Figure 2A shows the initial shape of a typical FE mesh used for the present computations, while Figure 2B shows the end-diastolic deformed shape of the FE mesh. The wall of LV model was discretized with 20-node hexahedral element; with the exception of the apical region which was meshed using 10-node tetrahedral element. The current FE mesh consisted of 22,080 total number of elements and 29,777 nodes. This discretization of the current LV model was sufficient and any further mesh refinement showed very little improvement [18]. The LV blood cavity was modeled by the hydrostatic fluid 3D solid element; this element is well suited for calculating fluid (blood) volume and pressure for coupled problems involving fluid-solid interaction. Hydrostatic fluid elements were overlaid on the faces of 3D solid element enclosing the fluid volume. Figure 2C shows the section view in the FE mesh in order to clarify the shape of the elements used for modeling the LV internal cavity. The hydrostatic fluid element was defined by nine nodes; eight nodes on the internal surface of LV cavity (endocardium) and the remaining ninth node at the base center, which is also called ''the pressure node''. This pressure node was used to define the LV pressure which was assumed to be uniform through the LV cavity; the predefined value of pressure was automatically moved to the centroid of the fluid volume. In all FE computations the LV cavity and LV wall volumes were kept constant at 50 ml and 73.6 ml respectively. The circumference of the LV internal cavity was divided into 48 equally spaced divisions, i.e. 48 elements along the circumferential direction.
Two separate parallel sets of 3D fiber network; contractile muscle fibers bundles (myofibers) bound by a mesh of collagen fibers were embedded within continuum 3D solid element to reproduce the globally anisotropic behavior of cardiac tissue. Computationally, these fibers were modeled as layers of uniformly spaced reinforcement bars (rebar) within the continuum 3D elements; each layer was set to be parallel to two of the isoparametric directions in the element's local coordinate system.
The collagen fibers were arranged in the radial direction, while the myofibers orientation changed with position within the LV wall. The 3D reinforcing element was used to simulate both myofiber and collagen fiber. The continuum 3D element is suitable for simulating reinforcing fibers with arbitrary orientations and used to model the myofiber force. The force was restricted in the direction of the fiber only (uniaxial fiber tension). The reinforcing element was firmly attached to its base element, i.e. no relative movement between the reinforcing element and the base was allowed. FE computations were conducted with myofibers and collagen fibers volume fractions of 0.7 and 0.015, respectively [19,20].

LV Myofiber Architecture
The human muscle fibers oriented at different angles throughout the ventricle wall in the form of sheets that are separated by a complex structure of cleavage surfaces (see Figure 3) [21][22][23]. The myofibers could be fully described by two inclination projection angles; the helix angle (b) and the transverse angle (g) in the two perpendicular planes. The transmural distribution of helix angle (b) varies in a linear manner through wall thickness from the lowest negative value at the epicardium to the highest positive value at the endocardium, while the transverse angle (g) varies in a linear manner through along the longitudinal axis of the LV from the lowest negative value at the apex to the highest positive value at the base [24,25,21,26,27]. The corresponding variations of the fiber helix angles (b) through the eight regions are as follows (see Figure 4); septum-basal region (260u : +40u), anterior-basal region (240u : +60u), lateral-basal region (220u : +50u), posterior-basal region (220u : +60u), septum-apical region (250u : +40u), anteriorapical region (220u : +60u), lateral-apical region (220u : +50u), and posterior-apical region (220u : +60u) vary smoothly across the LV wall thickness from a negative angle at epicardium to positive angle at endocardium respectively [28].
The distribution of transverse angle (g) was taken as a variable value in a linear manner from +15u at the base to the circumferential direction (g = 0u) in the equatorial region to 215u at the apex [24]. Meanwhile, the collagen fibers were arranged in the radial directions.

Loading and Boundary Conditions
To demonstrate the performance of the proposed FE model, LV pressure versus time curve for a healthy human heart was used (see Figure 5). This case was adapted from the experimental measurements carried out by Hall [29]. The LV pressure versus time curve given in Figure 5 represents the FE model with applied loads, while the accompanying LV internal cavity volume were adopted as the FE model target values.
The external surface of the heart is affected by the surrounding organs (lungs, ribcage and diaphragm). In order to simulate the boundary conditions imposed by these surrounding organs and tissues, an elastic foundation 3D structural surface effect element with a stiffness K f = 0.02 kPa was used [4,30]. Due to the lack of information about the influence of the surrounding organs and tissue on the deformation of the heart, a uniform elastic foundation was assumed.
To prevent rigid body motion of the model, degrees of freedom for all nodes at the base were suppressed in the longitudinal direction (UY = 0). To avoid possible excessive deformations of FE mesh elements, the pressure node was fixed laterally (UX = UZ = 0) (see Figure 2).

Myocardium Active and Passive Material Properties
The LV pressure response to its changing volume during ejection phase relies on the active elastance property activated by the muscle action. During systolic phase, the muscle generate adequate contractile force (muscle active force) to provide sufficient LV pressure to open the aortic valve, and pump an appropriate volume of blood [31]. At the beginning of the isovolumic contraction, the LV internal cavity pressure increases rapidly until the peak pressure and consequently the muscle contraction force continually increase, to sustain increasing LV pressures. Both the LV pressure and muscle contraction force increased simultaneously until its peak. During cardiac cycle the LV wall is subjected to dual forces; the active force generated by the myocardium muscles and the force generated by the blood pressure in the LV cavity. The active muscle force is not only governed by the myocardium passive properties but also depend mainly on myocardium active elastance operating throughout the cardiac cycle.
In this study we suggest a simple method to calculate the active myofiber elastance properties by taking into consideration that the maximum value of myofiber Young's modulus does not exceed 0.5 MPa [32,33]. The active myofiber elastic properties during cardiac cycle can be simply calculated by multiplying the LV pressure value (given in Figure 5) with constant value that equals to 29.5, i.e. the value of active myofiber Young's modulus depending linearly on the intracavital pressure. This constant value is equal to 0.5 (maximum myofiber Young's modulus) divided by the maximum value of LV pressure. The computation of active myofiber Young's modulus based on the above calculated constant represents the hypothesis of the present study. Figure 6 shows the calculated active myofiber Young's modulus versus time during one cardiac cycle. The time dependent calculated values of active myofiber Young's modulus is applied on the FE model in order to calculate the contraction force in the myocardial wall.
The myocardium tissue (matrix) was represented as an isotropic, slightly compressible hyperelastic material with relatively soft properties. The Ogden model with two parameters was used in the present study and having the following parameters; m 1 = 0.22 MPa, m 2 = 0.11 MPa, a 1 = 11.77 and a 2 = 14.34 [34]. The tissue bulk modulus K of the LV model was tuned so that the simulation behaves in accordance with the measured pressure vs. volume curve (patient-specific datasets).
The collagen fiber behavior was represented by isotropic linear elastic with large displacements, to simulate the large strains occurring in the collagen fiber during LV filling. Collagen properties are as follows; the Young's modulus (E) = 50 kPa, Poisson's ratio (n) = 0.49, and density (r) = 1000 kg/m 3 [35]. Figure 7 shows the inverse FE computation sequence procedures for the evaluation of the myocardial bulk modulus. The myocardium tissue bulk modulus satisfying the required output was inversely identified. The LV pressure versus time curve (see   Figure 2C). The other passive material properties of myocardium tissue (see item 2.4) were fixed as material constants. The calculated active myofiber Young's modulus (see Figure 6) was used to simulate the muscle active contraction force generated through the LV wall during the cardiac cycle.

Solution Procedures and Inverse Identification of Myocardium Tissue Compressibility
An initial guess value of tissue compressibility for the myocardium was then applied and by successive computations this was refined until the calculated LV cavity volume matched the measured volume. Iterative FE computation for the LV cavity volume during the cardiac cycle was carried out. At the end of each computation step, the predicted LV cavity volume of the FE model was compared to the measured value. The computation ended if the relative error between the computed and measured values #1%. For this calculation, the total computation time for each run was taken about 1000 s using PC Intel Core i7 (2.93 GHz with RAM 2.00 GB). At the start of all FE computations the LV wall initially was assumed to be stress free (LV cavity pressure equals zero).

FE Model Merits and Limitations
A new FE inverse model is presented that features the ability to determine the myocardial tissue bulk modulus during cardiac cycle. Other advantages of the proposed method include: 1. Slightly compressible hyperelastic tissue properties. 2. Realistic boundary conditions based on MRI observations. 3. Myofibers orientation simulated based on the obtained data from MRI. 4. The effect of interaction between blood and the internal cavity of LV wall. 5. The effect of surrounding organs and tissues on the deformation of the heart. However, this model is limited in a number of important ways, including: 1. Simplified geometry for the LV.
2. An incomplete understanding of some heart diseases. 3. The model's inability to study effects of electrical activation, blood flow, porous medium, and cardiac metabolism. 4. The used measured data taken for a healthy human heart ''ideal proband''. 5. More realistic tissue mechanical properties of LV are still needed.
These limitations and weaknesses can be used as the bases in future model improvements. Figure 8 shows the comparison between the predicted FE and experimentally measured LV cavity volumes. The LV cavity volumes increased rapidly from 110 ml to 130 ml (end diastolic volume EDV = 130 ml), shortly after the beginning of atrial systole phase through time = 0.1 sec. Then, LV cavity volume remained constant during the isovolumic contraction phase until time = 0.21 sec. After that, a sudden decrease in LV cavity volume can be seen at the onset of the rapid ejection phase until time = 0.3 sec, followed by a slight decrease during the reduced ejection phase until time = 0.43 sec, at the end systolic phase. The size of LV cavity volume at the moment was equal to that at the end systolic volume (ESV = 50 ml). The LV cavity volume remained constant during the isovolumic relaxation phase until time = 0.5 sec, followed by a rapid increase during the rapid filling phase until time of 0.65 sec. Finally, the LV volume slightly increased during the reduced filling phase up to the end of cardiac cycle.

Volume Variations of LV Cavity During One Cardiac Cycle
The LV cavity volume was locally adapted in response to several stimuli, such as systole phase (isovolumic contraction and ejection) and diastole phase (isovolumic relaxation, rapid filling, reduced filling and atrial contraction). Comparing the present simulation results of LV cavity volumes with corresponding values reported in groups of normal healthy subjects in previous studies [36][37][38][39], it was found that the present model can accurately predict the change in the LV cavity volume during the cardiac cycle. It can also be found that the various events (cardiac cycle phases) can be sharply defined. It should be noted that most of those studies ignored the effect of atrial systole on LV volume while others failed to accurately represent all cardiac cycle phases precisely isovolumetric contraction and isovolumetric relaxation.

Variations of Tissue Compressibility during One Cardiac Cycle
During cardiac cycle, the myocardium wall tissue is exposed to successive active contraction and relaxation in consequence of depolarization and repolarization, respectively. Due to heart beating and LV pressure dynamics response, a significant amount of energy can be expended to compress the heart wall, essentially squeezing the myocardium cells closer together. Figure 9 shows the FE results for the myocardium tissue compressibility variations during one cardiac cycle. It can be noticed that the myocardium tissue is nearly incompressible during a short period of time about 0.16 sec throughout the cardiac cycle (approximately not exceeding 20% of total cardiac cycle time). The tissue showed incompressible behavior during the reduced ejection and isovolumic relaxation phases. The myocardium tissue incompressibility is a common notion used in numerical simulations and analytical studies [40][41][42][43]. This is mainly for the purpose of simplifying the analytical formulations and the interpretation of experimental data.
From Figure 9, it can be noticed that myocardium tissue compressibility decreased rapidly at the beginning of the rapid ejection phase. After that, myocardium tissue compressibility increased gradually at the rapid filling phase and remained nearly constant during reduced filling, followed by further increase during the atrial systole (increased about 30% of the maximum of tissue compressibility). Increasing of myocardium tissue compressibility leads to LV wall growth and remodeling. The ability of the present FE model to simulate the observations of growth and remodeling consistent with the previous theoretical study obtained by Kroon et al. [40] and also with the empirical observations using diffusion tensor MRI carried out by Teske et al. [44]. The difference between the maximum myocardium tissue compressibility is 3 kPa 21 and its minimum value (< 0 kPa 21 ) can be used as a good index for the normal ventricular function. The decrease in this value means abnormal cardiac function occurrences, maybe due to myocardial infarction or heart dysfunction. The myocardium tissue compressibility increased with decreasing in the contractile force of the ventricles. This lead to an increase in LV cavity size which means that the heart cannot pump blood efficiently and structural alterations of the myocardium (i.e. heart is enlarged and the heart's pumping ability is impaired). Hence, the myocardial tissue compressibility should be considered if myocardial performance, myocardial deformation and heart wall stresses, response time, is critical.
Like many biological tissues, myocardium tissue is able to adapt to changes in mechanical load through growth (change in mass) and remodeling (change in tissue properties) [40]. This hypothesis has confirmed by the results obtained using the present FE model. The LV wall volume was locally adapted (change in tissue compressibility) in response to several stimuli, such as early systolic fiber stretch, fiber shortening during ejection, and contractility. Figure 10B shows the FE results for the variations of myocardium bulk modulus during one cardiac cycle. Great variations in the values of the myocardium bulk modulus occurred during the rapid ejection and isovolumic relaxation phases. The myocardium bulk modulus reached to its maximum value at end of ejection and began to decline at the beginning of isovolumic relaxation, i.e. the tissue of LV wall was stiffened by contraction and softened by relaxation. It can be noticed that the myocardium bulk modulus increased exponentially during ejection phase till its peak value, then followed by a linear decrease during isovolumic relaxation phase (i.e. the myofibrils return to their original length). The peak value of bulk modulus occurred at 0.43 sec and the duration time for bulk modulus changes from 0.3 sec to 0.5 sec.

Variations of Bulk Modulus during One Cardiac Cycle
With regard to the timing, the computed duration times for bulk modulus changes were compared to the electrocardiogram (ECG) during one cardiac cycle (see Figure 10C). It can be noticed that the duration times for the onset and ending of LV repolarization marked by T wave (on surface ECG) agree very well with the onset and ending duration times of bulk modulus changes (see Figure 10C). Also, it can be noticed that the onset and ending times of LV repolarization from 0.3 sec to 0.5 sec are sharply defined. Actually, the electrical signals were not included in the present FE analysis, but their impacts are included in the present FE model by introducing the myofiber active elastance. From Figure 10A-C, it can be noticed that there is a good correlation (synchronization) between the instantaneous variation of myocar- dium tissue bulk modulus and the onset and ending of LV repolarization (T wave).
The bulk modulus for myocardium tissue predicted by the current FE model is considerably less than the experimental values measured by Masugata et al. [45]. This is potentially a result of the bulk modulus for myocardium tissue changed drastically immediately after death. Also, there are challenges in the proper acquisition of human myocardium tissue samples and protocols for appropriate experimentation. Diversity of results can also be explained by the differences in the FE models simulation and experimental conditions or the choice of model parameters.
In conclusion, it was concluded that cyclic variation of bulk modulus as predicted by FE model exists during the cardiac cycle with the myocardium tissue being stiffer in systole than it is in diastole. Such behaviors was observed in the experimentally measured shear stiffness for normal myocardium throughout the cardiac cycle using Magnetic Resonance Elastography (MRE) by Kolipaka et al. [46].

Comparison between Predicted FE Bulk Modulus and Ejection Fraction
The predicted myocardium bulk modulus (K) and ejection fraction (E f ) were compared. Four different LV pressure-volume diagrams with different sets of physiological conditions were used, see Figure 11A, B. The initial LV cavity volumes of the used models were 50 ml, 50 ml, 65 ml, and 70 ml for Model_1, Model_2, Model_3, and Model_4 respectively, while the LV wall thickness was kept constant as described in Figure 1. Figure 11C shows the variations of the predicted myocardium bulk modulus during one cardiac cycle. It can noticed that a discrepancy among the peak values of myocardium bulk moduli exists; 1500 kPa, 2855 kPa, 4760 kPa, and 8000 kPa, which correspond ejection fractions of 61.5% [29], 58.3% [47], 56.3% [48] and 53.6% [49] respectively. Figure 12 shows the variations of the maximum bulk modulus versus ejection fraction. It is clear that the ejection fraction increased with decreasing peak values of myocardium bulk modulus. Such decrease (i.e. increases of myocardial tissue compressibility) caused an increase of myocardial contraction, which led to an increase of heart ejection fraction. Further study is still required to verify the correlation among the myocardium tissue bulk modulus as a marker for heart function, strain and  strain rate. Finally, the above results confirm the hypothesis, in the computing the active myofiber Young's modulus and the usage of inverse FE analysis for the determination of myocardium tissue bulk modulus.

Conclusions
A novel method combining FE inverse model and LV pressurevolume diagrams was developed to determine the myocardial tissue bulk modulus. The human LV wall (for simplicity) modeled as a thick-walled ellipsoid truncated at two thirds of major axis with spatial myofiber angle distribution was used. The ellipsoidal geometry has been chosen to model the human LV, as it is close to the real anatomical shape, and yet quite simple. The myocardial bulk modulus of different LV pressure-volume diagrams (with four different sets of physiological conditions; stroke volume, LV maximum pressure and ejection fraction), was determined. Based on the results and discussion presented in the preceding section, the following conclusions can be drawn: 1. The myocardium bulk modulus can be used as a diagnostic tool (clinical indicator) of the heart ejection fraction. 2. Our results show that the present FE model is sensitive to the overall cardiac function parameters expressed in terms of LV pressure-volume variations during cardiac cycle and ejection fraction. 3. The calculations of active myofiber Young's modulus (myocardium active properties) based on the LV pressure proved their correctness.
With this research, the authors would like to recommend further investigations on the subject of compressibility of myocardium tissue, which is still debated and remains a challenging experimental topic.