The ES pressure values for LBBB patients measured through aortic valves.
Detection of regional ventricular dysfunction is a challenging problem. This study presents an efficient method based on ultrasound (US) imaging and finite element (FE) analysis, for detecting akinetic and dyskinetic regions in the left ventricle (LV). The underlying hypothesis is that the contraction of a healthy LV is approximately homogeneous. Therefore, any deviations between the image-based measured deformation and a homogeneous contraction FE model should correspond to a pathological region. The method was first successfully applied to synthetic data simulating an acute ischemia; it demonstrated that the pathological areas were revealed with a higher contrast than those observed directly in the deformation maps. The technique was then applied to a cohort of eight left bundle branch block (LBBB) patients. For this group, the heterogeneities were significantly less pronounced than those revealed for the synthetic cases but the method was still able to identify the abnormal regions of the LV. This study indicated the potential clinical utility of the method by its simplicity in a patient-specific context and its ability to quickly identify various heterogeneities in LV function. Further studies are required to determine the model accuracy in other pathologies and to investigate its robustness to noise and image artifacts.
- acute ischemia
- left ventricle
- cardiac mechanics
- GE Healthcare US
Realistic finite element (FE) simulation of regional defects caused by pathologies such as acute ischemia remains a significant problem for physiologically based FE analysis of the heart. From the moment that the cell death starts to the thinning and the expansion of the scar surface of myocardium tissue, it takes approximately 21 days before the remodeling phase starts . Therefore, an accessible method to detect and characterize the influenced region can provide a better plan for diagnosis of these patients with medications and myocardial injection to avoid tissue progressive remodeling [2, 3].
Tracking the deformation of the LV has been actively pursued for the last decades, but important information is lacking about in-vivo wall stress reconstruction in medical imaging modalities for computational purpose . Accurate FE models based on high-resolution medical images can be applied for stress reconstruction, but these models are still limited to medical research and not generally available in the clinic.
Generally, the primary tool for ischemia detection is the analysis of standard electrocardiogram (ECG) signals. The majority of research focus regarding mathematical and computational tools of ischemic zone detection were addressed by electro-physiology models [5, 6, 7, 8]. An alternative method for ischemia detection is to assess cardiac function with magnetic resonance imaging (MRI) modalities. However, the abnormal zone detection in MRI modalities is limited and remains a challenging problem [9, 10, 11]. In this context, the speckle tracking technique which employs the non-invasive ultrasound (US)  has become a commonplace technique in clinical daily practice  due to the accessibility of this system.
In this chapter, we study the possibilities of computing a region of abnormal function from 4D strain acquired with the EchoPac® system. The strain estimates are based on speckle tracking, and capture 3D strain through time for 17 LV regions based on the standardized segments of the American Heart Association (AHA) .
Coupling US data with computations based on the FE method may provide more detailed information about the mechanical state of the patient’s heart. The purpose of the present chapter is to apply FE simulations to detect regions of abnormal function in the LV, based on the hypothesis that a homogeneous contraction is a regional assumption for a healthy LV, and large deviations between measured deformations and FE simulations based on homogeneous contraction will correspond to an abnormal region. We aim to investigate the validity of this hypothesis using the 4D strain data acquired with EchoPac®.
The chapter is organized as follows. First, we present an overview of our pipeline and of the equations for the systolic phase of the heart cycle. Next, we test the accuracy of the method using synthetic data, simulating zones of acute ischemia in the LV. Finally, the methodology is applied to a cohort of patients with left bundle branch block (LBBB), and we discuss the results and the strengths and limitations of the proposed method.
2. Material and method
Akinetic or dyskinetic regions such as ischemic regions or infarcts may be detected in medical images as regions of abnormal deformation. However, healthy tissue surrounding a pathological region may also show abnormal patterns of deformation, which makes it challenging to locate the pathological region. The purpose of the present paper is to investigate an alternative method for detecting abnormal regions, by comparing measured strains with strains resulting from a homogeneous contraction. As noted earlier, the underlying hypothesis is that the healthy contraction is approximately homogeneous [15, 16], and the largest deviation will therefore be associated with the pathological region. To test our hypothesis, we propose to:
Build a patient-specific mesh from the images (using an in-house mesh morphing method ),
Apply a uniform contraction and tune this to give the closest match to the patient’s strain data,
Evaluate the spatial deviations, and compare the region of largest deviation to the (known) regions of abnormalities.
These main steps of the methodology are illustrated in Figure 1. When applying the algorithm to the patient data, it is necessary to generate a patient-specific FE mesh for each individual case. Besides, the accuracy of the geometry associated with the generic LV model is greater than the patient-specific geometries, hence, the need to adapt the generic model to each patient. This task is accomplished by a previously developed mesh morphing method, which takes a generic LV mesh and adapts it to the patients’ image data . The rest of the methodology consists of a minimization loop on the active tension af. For each patient-specific mesh, the algorithm assumes a constant value for the actively generated stress, denoted as average active stress (AAS), and performs a systolic FE simulation. Then, the resulting regional strains are compared with the patient’s data, and a minimization procedure is performed to minimize the deviation between measured and simulated regional strains. The final, minimized, deviation is then visualized and analyzed in order to locate regions of dysfunctional LV segments.
The main reason why we focus on the systolic phase is that in acute ischemia, the passive behavior of the abnormal zone during diastolic phase stays intact, whereas the first appearance of the defected region is during the active contraction . As a result, remodeling of abnormal tissue induces a different deformation due to the average stress in regional LV motion.
2.1. Reference FE simulation model for systolic phase
To examine the abnormal tissue size effects on global and regional deformation, we need to explain the model which was used for FE systolic phase. We employed a healthy volunteer LV extracted from 4D US  at end diastole (ED). A 3D mesh with eight-noded linear hexahedral elements was created using Abaqus® mesh generation software. The reference mesh had 155,172 nodes and 141,405 elements.
To incorporate the known material anisotropy of the LV during the systolic phase, we defined a local curvilinear coordinate system aligned with the local fibre orientation, to model myocyte contraction [19, 20]. We assumed a helical fibre distribution, with a helical angle varying from −70° at the epicardial surface to 60° at the endocardial surface (0° is the circumferential direction). At each point, we defined the local orthonormal basis as (), where () is aligned with the local fibre .
We were interested in computing the mechanical state of the tissue at ES. The tissue deformation in the systolic phase results from a combination of active and passive mechanical properties. We modeled this interaction using a so-called active stress approach, where the total tissue stress is decomposed into an active () and a passive () part.
Introducing Rf as the matrix of transformation from global to local coordinate system, F as the deformation gradient, and (J = 1 for incompressible tissue), the active stress in global coordinates can be expressed as
where af is the active tension developed by the contracting fibres.
Following our assumption of homogeneous contraction, the value of af is assumed constant in space, and is tuned to provide the best possible match with the measured strain data for each patient. The chosen maximum value for af was set to 135 kPa [22, 23]. The passive component in Eq. (1) was modeled as hyperelastic, with stress derived from a strain energy function W;
Here E is the Green–Lagrange deformation tensor.
In the literature, different strain energy functions have been proposed for the cardiac tissue [24, 25, 26]. To reduce the complexity and number of parameters in the model, and based on the assumption that at ES the tissue stress is dominated by the active component, we used a simple Mooney-Rivlin strain energy model, as previously used by Marchesseau et al. :
Here c1, c2 are material parameters (MPs), and are the invariants of the right Cauchy-Green tensor and is the bulk modulus.
We performed a reference FE simulation using the following MPs: c1 = 0.0187 MPa, c2 = 0.0198 MPa, d1 = 0.1697 MPa, and the contraction model outlined earlier. The boundary conditions were as follows: a pressure of 11.24 kPa was applied to the endocardial surface , the epicardial surface was unloaded, and the basal nodes were assigned to remain coplanar during deformation . This reference FE simulation was compared with results reported by Dorri et al.  and showed good agreement in terms of LV ejection fraction (EF) (33.83%), wall thickness change (18.7%) and Von Mises stresses, which are in the range of 100–150 kPa [23, 28, 29, 30].
We used these boundary conditions, MPs and contraction model for the patient-specific FE simulations of the provided cohort with the Abaqus® software solver.
2.2. Application on synthetic acute ischemia
Two reference meshes have been used for application of our developed pipeline.
2.2.1. Ellipsoid reference mesh
A truncated ellipsoid mesh has been generated using Abaqus® software. The height of the mesh was 9 cm, and the wall thickness was 1 cm. The fibre directions were assigned as outlined earlier. We defined a small zone at the equatorial level of the septal wall (see Figure 2) as acutely ischemic, i.e. non-contractile. The active stress in this region was set to zero, while the rest of the reference FE model, af = 160 kPa, was attributed. The same condition as explained in Section 2.1 in terms of boundary conditions, ES pressure, contraction model and MPs was used to perform this model. In order to relate the results to the echo images, we defined a midwall surface mesh that was split into 17 segments as shown in Figure 3 . We will refer to this mesh as the midwall mesh throughout this chapter. Regional deformations were calculated as average strains in each segments (longitudinal, circumferential and area strain), similar to the output of the EchoPac strain analysis.
To provide a first test of the hypothesis underlying this research, we performed an FE simulation with an AAS, and tuned the active stress parameter to match the strains of the synthetic reference model as closely as possible. The cost function was defined as the sum of the squared differences (LSQ) between the AAS strains and the reference strains, and a simple minimization was performed by gradually increasing af from 90 to 180 kPa in 5 kPa steps. The smallest cost function obtained with this procedure was used as the closest match.
2.2.2. Healthy reference mesh
The healthy mesh obtained from US imaging in Section 2.1 was used in this study case. The abnormal zone for this mesh was defined on the mid-cavity anterior wall (Figure 2). The boundary condition, systolic blood pressure, contraction model and the MPs were attributed as in the ellipsoid study case. The midwall regions were also generated in order to follow the regional deformations through the systolic phase. Reorientation of the fibre directions in abnormal LV tissue is an open question [31, 32, 33, 34]. Therefore, we leave the fibres unchanged in the abnormal zone as defined for a healthy LV case in Section 2.1.
2.3. Prediction of abnormal tissue zone for LBBB patients
A cohort of patients suffering from LBBB was available for studying our proposed method for abnormal zone detection (with chronic infarcts). Left bundle branch block (LBBB) is a cardiac conduction abnormality which causes the LV to contract later than the right ventricle. We presented, previously, our FE method for a healthy LV contraction. In order to introduce the patient data to the FE software, we developed an in-house code to morph the reference FE mesh to the patient geometries. We demonstrated the steps for estimating homogeneous active contraction value data of acute ischemia in Section 2.2.
Eight LBBB patients were recruited after informed consent at Rikshospitalet, Oslo, Norway. LV triangulated mesh geometries were acquired by a standard GE Healthcare echocardiography examination . The patient’s LV cavity pressures were measured via an arterial catheter inside the left ventricle’s cavity. The 17 regional strains were measured using 4D strain tracking EchoPac® US system manufactured by GE Healthcare . EchoPac® algorithm calculates the classical longitudinal, circumferential and area strains for each segment at the midwall between endocardial and epicardial surfaces from the segment’s dimensions  (Figure 4).
From area strain (AS), we can calculate the radial strain (RS) as detailed by Heimdal  from the incompressibility assumption (R = Volume/Area) of the cardiac tissue:
We morphed the reference mesh (refer to Section 2.1) onto the geometry of each patient at ED time step for which US geometries were acquired. For this, endocardial and epicardial nodes of the deformable mesh were projected onto the patient’s triangulated surfaces from a defined LV centerline  with rigid and non-rigid transformation methods employing an in-house developed Matlab® code. In order to deform the reference bulk model with intermediate nodes, we used FE elastic rigid body deformation employing the displacement vectors obtained from boundaries projection trajectories . An in-house program coded in Matlab® found the closest nodes to the midwall mesh nodal positions obtained from the EchoPac® system for each patient at ED. Therefore, we can follow the evolution of LV during its systolic phase by calculating the deformations of each segment’s mesh in area, radial circumferential and longitudinal strains.
For each patient, after application of the mesh morphing method, we employed the pressure values at ES from pressure curves obtained during the medical intervention (Table 1). This cavity pressure was applied as the boundary condition on the endocardial surface as previously mentioned in Section 2.1. For the sake of simulation costs, for each patient, we considered a set of af values from the literature starting at 60–280 kPa with 5 kPa of increments. Then, the LSQ simply takes the af which returns the minimum cost value.
|ES pressure (kPa)||16.8||14||14.27||15.47||12.67||16.53||14.53||8666|
In this section, we describe the results of applying the proposed method to synthetic reference data and patient data from a group of LBBB patients.
3.1. Application to synthetic data
For the ellipsoid mesh, the closest match was achieved for af=155 kPa, which gave a cost function value of 0.0096. The EF for the reference ischemic case was 41.32% and the closest match was 43%. Figure 5 shows the strains for the reference case (left) and the closest matching case with homogenous contraction (middle). The difference between the two is shown in the right panel, and demonstrates that the method identifies the correct location of the ischemic zone (starred).
Figure 6 shows the results of applying the method on the healthy LV geometry, with synthetically generated acute ischemic data. In this case, a minimum cost function value of 0.0387 was obtained for af=130 kPa. Strain values are shown for the reference ischemic case (left) and the closest matching case (middle), while the difference between the two (error) is shown on the right. We observe that the ischemic zone is easily detectable from the error plots.
3.2. Application on LBBB patients’ data
We applied the mesh morphing method explained in Section 2.3 to construct patient-specific FE meshes for 8 LBBBB patients, shown in the figure by Behdadfar et al. . Table 2 shows the minimum cost value and the AAS values for each patient. In addition, the measured and simulated ED and ES cavity volumes are shown, as well as the EFs for both cases.
|Case||CostValue||af (kPa)||EDV-S (ml)||ESV-S (ml)||EF-S%||EDV-P (ml)||ESV-P (ml)||EF-P%|
The deformations obtained from US are in circumferential, longitudinal and area strain of the midwall mesh. Then, the RS was obtained by incompressibility hypothesis of the LV segments as explained by Heimdal . The radial, circumferential and longitudinal strains are shown in Figure 7 for patients 4 and 7 where the pipeline detected the abnormal zone. The FE strain results were subtracted from the patient’s respective data and are shown in this figure as well.
Figure 8 shows the transversal and longitudinal cuts of LBBB patients’ geometries superimposed by the FE simulation results for the optimal AAS. These cuts are tuned to show the maximum error zone (from 17 segments), obtained from LSQ application of regional deformation. Patient 5 had convergence issues and was not considered for further analysis.
In this study, a new approach for detecting potentially infarcted areas was first validated on synthetic data and then applied to a cohort of 8 LBBB patients before resynchronization. The results revealed zones of moderate heterogeneities which were often of smaller thickness but the heterogeneities were significantly less pronounced than that revealed for a synthetic case of acute ischemia. The method also revealed that the tissue tends to be stiffer in the lateral wall of LBBB patients. It is very interesting to notice that the obtained results confirm what was already observed by other researchers using standard medical imaging. More specifically, Veress et al.  showed the typical septum-related abnormal wall motion and impairment of wall thickening during systole, which caused LV remodeling. Often, technologies such as MRI or PET (positron emission tomography) were used to characterize the infarct heterogeneities as well as histological images [36, 37, 38]. We examined if such ability to detect abnormal tissue can stand for novel 4D strain EchoPac® system.
In order to detect the abnormal tissue, we identified an AAS for all FE meshes which was the closest match to the patient’s midwall mesh deformation (Figure 1). This AAS represents a homogeneous contraction for each patient [15, 16]. In this condition, it is expected that for a healthy and homogeneous contraction, the FE deformation responses will be different in abnormal tissue than the rest of the LV model. In addition, this value represents the average movement of remote (uninfarcted) and injured tissue.
Often, patients suffering from abnormal contractility are evaluated by an optimized gradient of activation stress map. Wenk et al.  have developed an animal-specific FE model for such reconstruction by minimizing the regional deformation error between experimental data and the FE model. In this study, a border zone was considered for early stage of tissue damaging process (due to calcium concentrations) which relates the injured to the remote tissue. They confirm that the infarct zone has no contractile action (zero active value).
4.1. Impact of MPs on the FE LV deformation
In Table 2, the AASs are mentioned for optimal cost values. These identified values are strongly coupled to the wall thickness, the blood pressure and the cavity volume. The FE results are based on the midwall strain deformation which is not representative data for the complete wall motion and might be affected by artifacts. Some meshes are inflated such as in cases 4 and 7 in Table 2, even if, these cases are the most successful results of the abnormal tissue detection. The reason for this inflation might be the compliance MPs for the studied cases. To study this hypothesis, we increased the elastic parameters (c1 and c2) by factors [2, 4, 5, 7, 9, 11] to analyze the behavior of the patient FE simulation to this increase in tissue rigidity.
For this test, patient 7 has been selected as the posterior and neighboring walls were deformed dramatically in FE simulation while these walls had not been moved during ED to ES movement. The cost value for this test from 0.7901 with 175 kPa of AAS was reduced to [0.1893, 0.0775, 0.0626, 0.0592, 0.0591, 0.0564], respectively, and it is shown in Figure 4 for factor 12 (c1 and c2) in 7*. Novak et al.  showed that the aneurysmal wall is stiffer than the remote tissue as well. In this test, we observed that the results (7*) were improved by increasing the rigidity of the tissue MPs. It has been previously investigated, with a mathematical model and experimental tests, that the infarct tissue stiffen (despite dilated infarctions) . However, the infarct tissue properties are still under investigation for various accompanying pathologies such as LBBB.
Figure 8 shows several successful zone detections wherein the maximum error happens to be at the abnormal region. In case 1, the maximum deformation was observed at the Inf-Sept-AntSept regions where this neighborhood wall slightly inflates (blue lines at Max-error in longitudinal cut) and the method detected this region successfully. However, in case 6, the inhomogeneous strain was observed at the septal wall wherein the method detected the posterior wall as the maximum error. These results show that the strain data in some cases can be difficult to be judged on the detected zone and is subjective to the visual detection of the patient’s strain and geometry data. However, a FE simulation of pathology cases permit us to narrow down the potentially defeated zone to less regions, especially, in case 4 and 7 in Figure 7 and to study the mechanical state of the cardiac tissue by analyzing the wall stress as a post-processing step.
The results showed a good agreement with several patient deformations in RS, which has been shown to be a better indicator than circumferential or longitudinal deformations due to the homogeneous distribution of these deformation values, for discrimination of possible heterogeneities in the active stress and abnormal tissue movement. We observed that the identified active stress is highly coupled to the patient geometry and wall thickness . This is the advantage of the new generations of US system in extracting the RS from measuring the surface of a given segment that is not possible for modalities such as tagging MRI .
4.2. Remodeling in LBBB patients and its impact on the pipeline results
The LBBB, itself, causes tissue remodeling due to the redistribution of LV workload in the long run, especially, in circumferential shortening, cardiac mass, septal hypoperfusion and myocardial blood flow [44, 45, 46, 47]. However, this pathology is, often, the result of other cardiac diseases. Statistically, 40% of patients with cardiomyopathy (CMP), weakening of the heart muscle and congestive heart failure have unsynchronized ventricular contraction. The resynchronization procedure helps in reverse remodeling and hospitalization-free survival rate of 70–90% in these patients [48, 49, 50].
The observed regional differences did not exist in the case of ellipsoid and healthy meshes illustrated previously, where the wall thickness is nearly homogeneous. Thus, the change in MPs or wall thickness in LBBB patients with previous infarct history should reduce the efficiency of this method for some patients with high heterogeneity in the wall deformation due to tissue remodeling. In Figure 8 and Table 2, we can observe that the FE simulation responses in the maximal discrepancies with the patient’s data are mostly over or under estimates of the tissue deformation. These maximal discrepancies have been observed to be at the high heterogeneities in wall thickness.
Patient 5 had convergence issues due to large cavity volume (319 ml at ED) and thin wall thickness (remodeled tissue). The large cavity volume is related to an aneurysmal bulging which made contraction simulation fail to reach a solution.
Several studies showed that LBBB patients respond poorly (30–50%) to resynchronization procedure for multiple reasons and cases with developed scarred tissue tend to not respond at all [51, 52, 53]. This study is a tool for patient cases which are not yet affected by this remodeling phase and are in the early stages of LBBB pathology. Information on areas with pathology within the left ventricular wall may be helpful in planning therapy with a resynchronization pacemaker .
One major limitation in this study is the choice of material and contraction models to simplify the complexity and computational costs. In addition, we defined an af value which is also a brutal simplification for AAS identification procedure. However, with all this simplifications, the time of each simulation on a 12-core system was 72 h in average. One issue might be the small elements generated by developed mesh morphing procedure that reduced the time steps in the FE solver. One natural step can be improving the material model to anisotropic from the state-of-the-art . The mesh morphing method should also be improved to reduce the computational costs in further work.
Another limitation is the lack of patient’s history information. In this study, the patients were known to have LBBB pathology with infarct history but there was no other information on co-morbidities that may cause depressed myocardial contraction. The large cavity volume, low EF and low wall thickness can be related to various pathologies such as idiopathic dilated CMP, and ischemic cardiomyopathy. So in this work, it is not possible to categorize the patients with their pathologies (such as patient 5) and the identified AAS. We considered a stress-free configuration to simulate the systolic phase for our cohort so this issue should be further addressed.
In this study, a new approach for detecting potentially infarcted areas was introduced, validated and then applied onto a cohort of eight LBBB patients before their cardiac resynchronization therapy. The results revealed zones of moderate heterogeneities which were often of regions of thinner wall. However, the heterogeneities were significantly less pronounced than what was revealed for a synthetic case of acute ischemia, which we interpreted as an effect of remodeling induced by the therapy. The method also revealed that the tissue tends to be stiffer in the lateral wall of LBBB patients, as the deformations are not as pronounced as they are in the simulation. This study is promising for the assessment of LBBB and its quantification using FE simulations. Again, this illustrates the importance of patient-specific FE simulations in the domain of cardiac biomechanics.
Further work is required to transfer the promising synthetic results to real acute ischemic patient cases. A possible improvement is to reconstruct the patient fiber orientations from MRI imaging modalities such as diffusion tensor and tagging MRI to improve the deformation trajectories during contraction . The possibility to compare this study results with cardiac MRI for identifying the abnormal tissue is also an interesting future step. Finally, it would be interesting to compare the reconstructions made with our methodology on LBBB patients at different stages of remodeling after the cardiac resynchronization treatment.
The authors gratefully acknowledge the support and data provided by the Center for Cardiological Innovation, funded by the Research Council of Norway. SA is grateful to the European Research Council for grant ERC-2014-CoG BIOLOCHANICS.