Impact of mechanical deformation on pseudo-ECG: a simulation study

Introduction

The spatio-temporal morphology of the electrocardiogram (ECG) is the result of the motion of the cardiac extracellular action potentials with respect to the position of the electrodes. In the computation of the ECG, often a fixed geometry of the heart in the torso is assumed. Based on this approximation, patient-tailored ECGs can be simulated in a few minutes on high performance computing (HPC) architectures with anatomically detailed cardiac electrophysiology models.^1,2 However, the heart is mechanically active and several model studies showed that the T-wave in the ECGs is influenced by the volume of blood in the cavity³ and hence by the current configuration of the heart.^4–6 This indicates that the commonly used pure electrophysiology (PE) scenario may be too simple.

Actually, mechanics plays a double role in possible modifications of the ECG signals in contrast to a PE scenario: one related to the change of distance from the electrodes during the cardiac cycle and one due to the effects of mechanical contraction on the electrical activation front.⁵

The aim of this study is to investigate the impact of two above-mentioned mechanical effects on simulated pseudo-ECGs. We employ a mono-directional (MD) coupling scenario to isolate the role of the distance from the electrodes and a bi-directional (BD) coupling scenario to understand the combined role of the two effects. Mono-directional allows to remove all mechano-electric couplings (MECs) which may affect the action potential propagation. Additionally, this study results also in a comparison of the three computational scenarios.

Methods

Fully coupled bi-directional model

Mono-domain equation is derived from macroscopic balance laws of intra- and extra-cellular charge distributions.⁷ As in similar multiphase theories,⁸ this approach allows to extend the static mono-domain¹ system to moving cardiac domains. Cellular electrophysiology is described by adopting the properties of endocardial cells from the Bueno-Orovio model.⁹ The conductivity tensor $G$ is assumed to be transversely isotropic. Electrophysiology material properties are homogeneous in space and taken from the ‘N-benchmark’ paper.¹⁰

MECs are not considered in this work, except for the evolution in time of the cardiac configuration.⁷ To account for this, a convenient approach, routinely adopted in computational mechanics, is to map the mono-domain equation on a fixed reference geometry. The mapping yields to an equation equivalent to the mono-domain but with conductivity $F-1GF-T$⁠, where $F$ is the deformation gradient map, thanks to the Piola transformation.¹¹ As a matter of fact, the electric conductivity in the reference geometry is decreased (respectively, increased) if the tissue is locally stretched (respectively, contracted). Moreover, mapping activation patterns arising from different models into the same reference geometry allows for easy comparison of them.

The mono-domain system is coupled with the hyperelastic mechanical model that we presented and validated in.¹² The strain energy density is transversely isotropic and, differently from several other models,^4,13,14 we assume full incompressibility of the cardiac tissue. The passive stress is described by the Guccione–Costa strain-energy density function.¹⁴ Since the material is assumed to be fully incompressible, the latter depends on the deviatoric component of the strain.¹¹

Coupling of electrical activation to mechanics is modelled by assuming a transient of active stress development,^13–15 exerted only along the fibre direction. The temporal evolution of the force generated by myocytes is described by the model introduced in.¹⁵ In this simplified model, the active force depends only on a scaled trans-membrane potential which is multiplied by a prescribed maximum active force what we set to $Tamax=30 N·m-2$⁠. The time relaxation constant has been set to $ε=10 ms-1$⁠.

Experimental setup

We aim to systematically study the effect of the mechanical deformation on ECGs. Thus we reduce the complex effects due to geometry and fibre orientation, considering the following three test cases:

• A. a cube with edge of $7$ mm and fibres aligned along x-axis,

• B. a slab of size 20 $×$7 $×$3 mm and fibres aligned along x-axis, and

• C. a slab of size 20 $×$7 $×$3 mm and fibres aligned along y-axis.

The test case (B) is the electro-mechanical enrichment of the purely electrophysiological test studied in the ‘N-benchmark’ paper.¹⁰ In order to avoid rigid-body motion, the three faces along the planes x = 0, y = 0, and z = 0 are fixed in the normal direction.

We computed the ECG signals at six unipolar electrodes $Ei, ±$⁠, two for each Cartesian axis $i=x,y,z$⁠. They are arranged symmetrically with respect to the centre of the domain and are 10 cm apart from the boundary faces. In addition, the unipolar signals from electrodes aligned in the same direction are combined into a single bipolar ECG, with the electrode in the positive direction being the cathode and the opposite one being the anode, defining $Ei=Ei, + -Ei, -$⁠.

Reduced models: mono-directional and pure electrophysiology

An MD is employed in order to isolate the influence of the relative position of tissue to the electrodes from the change of shape of the tissue. In order to achieve that, the MD model does not consider the change in conductivity due to the deformation, i.e. the mono-domain equation is solved in its original form. This scenario is particularly appealing from a computational point of view, since electric activation and mechanical deformation can be computed in sequence. This fact motivates the ‘MD’ nomenclature.

The MD approach is equivalent to the ‘static-dynamic’ approach,⁴ in which the evolution of the potential is computed on a reference geometry (static step) and then the deformation due to the active force is calculated to obtain the current configuration (dynamic step). Potential from the static step is mapped on this deformed configuration where ECGs are actually computed.

Finally, a PE^1,2 is employed in order to remove also the effects due to the change of configuration in space. In this scenario, propagation of the action potential and the ECGs are both computed on a reference domain.

It is worth to remark that PE and MD provides the same activation pattern. In the latter one, mechanics is used to deform the geometry in order to compute ECGs. Instead BD considers also the changes of geometry in the monodomain system. Since, the potentials provided by MD and BD are different, also their current configuration will be different.

Solution methods

A novel software HART has been developed for the efficient computation of electrophysiological and electro-mechanical problems. Developed inside the finite element framework MOOSE,¹⁷ HART allows for an HPC solution of the three presented scenarios. A massively parallel variational transfer,^18,19 that enables seamless communication between two arbitrary grids, allows to exchange potential and deformation between the two solvers. This allows to employ coarser meshes to reduce the computational cost of mechanical problem and high resolution meshes to capture the upstroke in the action potential.²⁰

Three different numerical algorithms have been developed for the solution of the different scenarios (see Figure 1). Time integration of the mono-domain system is performed by means of an implicit–explicit (IMEX) scheme.²¹ First, gating variables are explicitly integrated employing a Rush–Larsen scheme. Then, the resulting currents are used as a source term for the diffusion problem, which instead is implicitly solved. The non-linear mechanical model has been solved with Newton’s method.

Figure 1

Schematic representation of the three simulation strategies. (Top) Pure electrophysiological approach. (Centre) Mono-directional (static-dynamic) approach: in order to remove the effects of the deformation, electrophysiology is solved on a fixed reference configuration. At each time step, the obtained action potential is transferred to the mechanical solver to update the active force and, hence, to compute the new configuration. Then the deformed configuration is used to computed the ECG taking into account the relative motion of the tissue with respect to the. (Bottom) The BD approach: same as the MD coupling but displacements are also transferred to the electrophysiological problem in order to solve the mono-domain equation in the deformed configuration. This allows to consider the changes of conductivity due to current configuration.

Open in new tabDownload slide

The computational meshes of the electrophysiology have a resolution of 0.125 mm. Since no steep gradients in the deformation are present, the mechanical solver employs a coarser mesh of step-size 0.5 mm. Both solvers employ a temporal resolution of 0.05 ms. The linear system arising from the IMEX linearization of the non-linear mono-domain equation has been solved with an algebraic multigrid solver (BoomerAMG) and the linear system arising from non-linear mechanics has been solved with a parallel direct solver (SuperLU).

Results

The action potential and the pseudo-ECGs calculated with the three computational scenarios are shown in Figures 2–4 for test cases A, B, and C, respectively. In order to compare zero potential isosurfaces from MD and BD approaches, the simulated action potentials have been mapped to the reference geometry. The isosurfaces of PE are equal to those of the MD since the electrophysiology problems are solved on the same domain for both the scenarios.

Figure 2

Test Case A. (Top Left) Snapshots of the spatial distribution of the action potential from the MD and BD for a cube of tissue with fibres aligned along the x-axis (marked in red). Configurations are reported 25 ms after applying an electrical stimulus at the location of the asterisk. The grid reflects the mechanical mesh of the static configuration. (Top right) Comparison of zero isopotential surfaces obtained from the MD scenario (red surface) and BD scenario (blue surface) for different time steps (10, 25, 30 ms). In order to compare them, the different activations have been projected on the reference configuration. (Bottom) Obtained pseudo-ECGs. The three bipolar ECGs along the Cartesian axes are shown here. Green, red, and blue lines represent results for the pure-electrophysiology approach, the MD approach, and the fully coupled approach, respectively.

Open in new tabDownload slide

The two leftmost upper panels of Figure 2 depict the deformed states (coloured) of the two electro-mechanical approaches and they are compared with the undeformed state (grid). The activation starts from the rightmost corner (asterisk) closed to the three clamped faces x = 0, y = 0, and z = 0. On the other faces, a stress-free boundary condition is imposed. The action potential is depicted as the membrane potential after 25 ms: the green region has been depolarized and the red region not yet. It can be seen that in the MD simulation the wavefront has progressed less, especially in the XZ- and XY-planes (parallel to the fibre direction, marked in red). Front of the BD method is also flatter.

The influence of deformation on the activation wavefronts at different times is reported in the right upper panel. In general, due to the contraction the BD gave a faster activation front in the interior of the tissue and close to the stress-free faces along fibre direction. Hence, with the fully coupled approach the same point activated earlier. Instead, in the cross fibre direction a slower activation has been noticed. An opposite behaviour has been observed close to the clamped edge along fibre direction. Specifically, the point located at (0,0,7) activated 0.2 ms earlier with PE and MD.

The influence of mechanical deformations is then reflected in the pseudo-ECGs. As a consequence of the simple axis aligned geometry, ECG potential was basically dominated by the activated areas orthogonal to positive-to-negative electrode direction, weighted by the distance from the electrodes. The differences in the propagation are mostly visible in the X-direction where larger deformations occur. Initially, the potential increased since the activated area on the $x=0$ face increased. Since no changes in configuration are present at this point, ECG are almost equal. Then, after reaching its peak value, when the opposite face started to activate the potential abruptly decreased. Here differences are reflected in the notching present in the QRS complex of X-direction around the S-wave. Mono-directional was responsible for a sharp notch, compared with the PE and MD. This was a consequence of the varying distance of the electrodes from the tissue. The activated areas between the two approaches were comparable, but they were located at different distances from the electrode $Ex,+.$ On the other hand, the difference between BD and PE approaches is less appreciable, due to the change of conductivity, that compensated for the relative variation of the electrode location with slower areal velocity. For the other two electrodes, the QRS complex duration and shape were similar among PE, MD, and BD approaches.

In general, PE and BD produced T-waves of different amplitude. The discrepancy varied among the electrodes and test cases, with either higher or lower amplitudes. Mono-directional always provided a good approximation of BD.

In order to show more clearly the effect of conduction and contraction along the fibres, simulations were performed in two flatter and longer tissue slabs. For test case B (Figure 3), the effect of the boundary condition on the propagation was more apparent due to the orientation of the fibres along the major axis. Mono-directional and bi-directional provided similar activation patterns until the signal reached the top face (around 7 ms after initiation). Then, since the signal has reached the stress-free top face, the tissue slab started to contract and thus the interplay between mechanics and electrophysiology became manifest. The propagation for BD approach resulted faster on the top face and slower close to the bottom face of the preparation. In particular, the zero-potential isosurfaces from MD and BD intersected on the reference geometry, certifying a face-to-face smooth variation in conduction velocity (see right panel of Figure 3).

Figure 3

Activation after 25ms, wavefronts, and ECG for test case B. See caption of Figure 2 for explanation. In the central panel a zoom of the three fronts.

Open in new tabDownload slide

Concerning the ECG, the PE simulation gave a slightly smaller potential at the peak of the QRS complex compared with the other two approaches while the Y-lead of the pseudo-ECG showed an opposite effect (see Figure 3). No large differences were observed between MD and BD. Actually, the tissue is long enough along the fibre direction to decouple the activation of the two opposite faces, and thus easing the ECG interpretation. Differently from the previous case, accounting for deformation in the BD resulted in a larger and slightly earlier T-wave along the x-axis than in the MD and PE approaches.

Finally, the test case C (reported in Figure 4) confirmed the results of the previous two cases, with a faster activation front of BD along fibre directions. In the ECG, differences were noticeable in the S-wave between the different approaches in the fibre direction (see bottom panel of Figure 4). As for test case A, larger potential in the QRS complex were given by the BD and even larger by the MD. A delayed and amplified T-wave was observed for the y-axis. Again, the MD approximation is not satisfactory for the computation of ECGs. Even if the overall shape of the QRS complex, overestimated values are observed (for example, peak of S-wave for y-axis was 0.09270 mV) and surprisingly a change of the sign of the T-wave was also present.

Figure 4

Activation after 25ms, wavefronts, and ECG for test case C. See caption of Figure 2 for explanation. In the central panel a zoom of the three fronts.

Open in new tabDownload slide

Discussion

Electrocardiogram is usually computed as static process, without considering the mechanical deformation due to the beating heart and breathing. Our results confirm the importance to account for mechanical deformation while simulating the ECG. Moreover, they significantly extend previous knowledge indicating that different ways of coupling electrophysiology and cardiac mechanics may result in different propagation patterns and ECGs.

Mechanical changes have been shown to have major effects on potential propagation both at the macroscopic level,²² depending on the position of the heart in the torso, and at cellular level.^5,23,24 In this study, we focused on the changes in the ECG signals due to the different propagation pattern given by the mechanical deformations and due to the changes in the relative distance from the electrodes. The three different approaches had the role to the separate these two effects.

At the T-wave, the shortening along the fibre direction has its maximal value, and hence the major differences between the electro-mechanical models and the PE one are noticeable. On the other hand, smaller differences are usually observed between the two electro-mechanical scenarios. This demonstrates that the change in conductivity is mostly a dynamical effect: deformations are slower during repolarization. Also Cluitmans et al.³ confirmed that the differences in the repolarization phase are mostly due to the difference location of the heart, performing a static comparison of electrograms in systolic and diastolic cardiac geometries.

Differences between the MD and BD scenarios, when present, become more relevant at the QRS complex. During the depolarization phase, faster and larger changes occur in the active force and hence in the deformation. The position of the tissue in the ideal torso appears to be counteracted in the BD approach by the effect of deforming myocardium on its conduction velocity. This can be seen observing that the calculated ECGs in the BD simulations were closer to those from the PE simulations.

Changes in both the QRS complex and the T-wave were observed also by Smith et al.⁵ In a two-dimensional comparison between BD and PE they motivated the early T-wave partially by the geometrical changes and partially by phenomena at the cellular level. In our simulations, an altered T-wave both in time and amplitude had been observed but, differently from the previous studies, time shift depends on the geometry and fibre orientation. Actually for the cubic case, reported in Figure 2, the T-wave with mechanics has a larger amplitude while for the tissue slab, reported in Figure 3, an opposite effect was noticeable.

Our results confirmed the findings of Smith et al. also concerning the QRS complex. The BD approach is responsible for a larger potential in QRS complex. The MD approach in this cases is responsible also for an early S-wave and for a notch that are in particular visible for the leads located along the propagation direction (see Figures 2 and 4).

The expected larger conduction velocities are usually visible on the non-clamped sides of the preparations. On the other sides, PE and MD had a faster velocity. This means that mechanical boundary condition may have a significant impact on the propagation of the activation potential. This effect is mostly challenging in its understanding but, in any case, it suggests that mechanical boundary conditions have to be accurately chosen in fully coupled simulations in order not to introduce disturbance on the propagation, and hence invalidate the results.

Notching and impact on local conduction velocity

The use of the MD coupling, and to a lesser extent the BD coupling, resulted in some cases in QRS complexes characterized by a notch that was less visible in the PE simulations. This observation suggests that the occurrence of notches in the QRS complex may in part be related to the effect of deformation.

A well-known ECG condition, where notching appears, is left bundle branch block (LBBB). In LBBB, due to a delay (or disruption) of conduction of the left bundle, the right ventricle is activated first and then the activation front spreads from right-to-left through the septum. The typical notching observed in the ECG of patients with LBBB, may be due to a slow trans-septal conduction. Our simulations, however, provide a possible additional explanation for notching in LBBB. The septum undergoes a complicated paradoxical deformation, often referred to as ‘septal flash’²⁵ while the action potential is travelling through it LV. This is similar to ECG along fibre direction occurring, for example, reported in Figures 2 and 4. In this simulation, propagation was orthogonal to the fibres, and thus the cube or the slab started to contract before the entire block had been activated. Because of the incompressibility of the material, a contraction in the fibre direction translated into an expansion in the cross-fibre plane, thus introducing a relative motion of the action potential with respect to the electrodes.

Limitations

The model adopted in this paper comes with some simplifications: the active force generation is phenomenological and not integrated into the cellular model. Moreover, used fibre orientations were uniform and aligned along one axis of the ECG electrodes. Also, the rectangular arrangement of the six ECG electrodes and the space they surround are a simplification of the actual shape of a human torso. These simplifications were used in order to more clearly demonstrate the proof of principle of the impact of deformation on the ECG. Further investigations are, hence, necessary in order to understand the role of fibres distribution, e.g. employing slabs with non-uniform fibre directions, and of the geometry, e.g. employing a truncated ellipsoids or realistic hearts.

Conclusions

In this study, we developed a BD fully coupled model of cardiac electro-mechanics. Depending on the tissue thickness and fibre orientation, also the QRS complex, and not only T-wave, was shown to be influences by the mechanical deformation. Introducing the MD approximation, we showed that the V-shape present in some ECG signals may be due a wrong simplification of the BD model or, reversely, may suggest a possible mechanism for the notching observed in LBBB patients.

Funding

The authors gratefully acknowledge financial support by the Theo Rossi di Montelera Foundation, the Metis Foundation Sergio Mantegazza, the Fidinam Foundation, and the Horten Foundation. This work was also supported by the Swiss National Science Foundation, project 205321_149828 ‘A Flexible High Performance Approach to Cardiac Electromechanics’.

Conflict of interest: none declared.

References

Appendix 1: Note on choice of the active force

The choice of the active stress parameters, in some cases, led to an unstable behaviour of the mechanical solver. This instability generally occurred during the repolarization phase and it was due to a singularity in the linearized mechanical problem. A possible explanation comes from a buckling phenomenon given by active stress model under compression:^26,27 fibres can exert an active force only under extension.