## New developments in a strongly coupled cardiac electromechanical model

## Abstract

**Aim** The aim of this study is to develop a coupled three-dimensional computational model of cardiac electromechanics to investigate fibre length transients and the role of electrical heterogeneity in determining left ventricular function.

**Methods** A mathematical model of cellular electromechanics was embedded in a simple geometric model of the cardiac left ventricle. Electrical and mechanical boundary conditions were applied based on Purkinje fibre activation times and ventricular volumes through the heart cycle. The mono-domain reaction diffusion equations and finite deformation elasticity equations were solved simultaneously through the full pump cycle. Simulations were run to assess the importance of cellular electrical heterogeneity on myocardial mechanics.

**Results** Following electrical activation, mechanical contraction moves out through the wall to the circumferentially oriented mid-wall fibres, producing a progressively longitudinal and twisting deformation. This is followed by a more spherical deformation as the inclined epicardial fibres are activated. Mid-way between base and apex peak tensions and fibre shortening of 40 kPa and 5%, respectively, are generated at the endocardial surface with values of 18 kPa and 12% at the epicardial surface. Embedding an electrically homogeneous cell model for the same simulations produced equivalent values of 36.5 kPa, 4% at the endocardium and 14 kPa, 13.5% at the epicardium.

**Conclusion** The substantial redistribution of fibre lengths during the early pre-ejection phase of systole may play a significant role in preparing the mid-wall fibres to contract. The inclusion of transmural heterogeneity of action potential duration has a marked effect on reducing sarcomere length transmural dispersion during repolarization.

- cardiac electromechanics
- computational model
- finite element

## Introduction

The integration of cardiac electrical and mechanical experimental data and hypothesizes across multiple spatial scales and functions via mathematical modelling is arguably the most advanced example of a ‘Physiome’ [1] style framework.

Detailed electrophysiological models of myocyte cellular [2,3] and subcellular processes incorporating genetic mutations [4] are now being linked to a large body of existing research for modelling the spread of electrical excitation throughout cardiac tissue and the whole heart (for reviews see [5–7]). Numerous methods have been developed to model the spread of electrical excitation varying from discrete models in which individual cells are electrically coupled to their neighbours via variable resistive elements (gap junctions) through to continuum based models. Similarly, detailed finite element continuum models incorporating anatomy, structure and passive mechanics have been coupled to cellular active tension (see [5,7] for reviews).

The development of models of these two physical processes has occurred largely independently. However, coupled models of cardiac electromechanics which trigger mechanical contraction based on the calculated activation sequence have recently been developed to investigate the relationship between activation and contraction in the normal [8,9] and paced heart [10,11]. It is indicative of both the computational load and modelling complexities that in each of these studies the activation sequence has been calculated independently and then later used as the excitation wavefront to initiate contraction. Such an approach negates the possibility of incorporating the effects of electro-mechanical coupling due to, for example, length-dependent calcium binding or stretch-dependent conductance changes that are characterized in some of the more recently developed myocyte models [2,12] and have also been shown to be important in tissue models [13,14].

While we can expect computational performance to continue to improve, in order to incorporate these coupling processes, or indeed new advances such as signal transduction models [15], gene regulation models [4] or metabolic models [16], a method is required which strongly couples tissue stretch and activation to cellular tension and membrane potential.

The method we have developed uses cellular models of electromechanics to drive the dynamic active material properties of a continuum organ model while the cellular models themselves get feedback from both the electrical excitation and the deforming mechanical model. The mechanical deformation of the geometric models also feeds into the excitation model affecting the propagation of electrical activation and cellular electrophysiology.

In this study we use the developed framework to examine explicitly the effect of adding transmural electrical heterogeneity in the embedded cell models on global myocardial mechanics.

## Methods

A computationally efficient cellular model was developed by coupling the Fenton and Karma [17] simplified activation model to the cardiac mechanics of Hunter et al. [18] which characterizes the length and calcium dependence of tension generation. This coupled model is termed below as FK-HMT. For the purposes of coupling, a calcium transient (not explicitly included in the Fenton and Karma formulation) was calculated from the slow inward current based on the original intracellular calcium formulation in the Luo and Rudy [19] model on which the Fenton and Karma is based. Following the method described by Nickerson et al. [12], the activation model was then used to provide the calcium transient input required to drive the active and dynamic development of cellular tension.

Three different parameter sets were used in this formulation to characterize the heterogeneous behaviour of myocardium produced by midmyocardial or M cells, located in the middle of the ventricular wall [20]. Parameter values of τ_{r}=150 ms for midmyocardial cells, 144 ms for epicardial and 142 ms for endocardial cells resulted in calculated single cell action potential durations (APD) of 437 ms, 352 ms and 329 ms, respectively, matching values found in human ventricular tissue [21].

This cellular model was then embedded in a simplified geometric model of the left ventricle (shown in Fig. 1a) based on a prolate geometry with 50 mm outer radius (at the base) and 60 mm longitudinal base apex dimensions represented using 128 tricubic Hermite finite elements in rectangular Cartesian coordinates. Wall thickness was set at 10 mm at the base tapering at the apex. The simple distribution of the epicardial, midmyocardial, and endocardial cell types prescribed transmurally in the model is illustrated in Fig. 1b. The microstructure of the models was defined such that the fibre angle varied sigmoidally over 120° between the endocardium and the epicardium.

The material properties for the electromechanical stimulations are referred to in this local embedded microstructure. Recent studies have shown that conclusive evidence for the use of fully orthotropic conductivities in tissue level models is not yet available [22]. Thus, the myocardial model was treated as transversely isotropic using a ratio of 10:1 [23], longitudinal:transverse to the fibre direction; parameters for the mono-domain activation model are shown in Table 1.

Material axis | Value (mS mm^{−1}) |
---|---|

Fibre conductivity | 0.263 |

Sheet conductivity | 0.0263 |

Normal conductivity | 0.0263 |

The mechanical constitutive properties were also referenced to the microstructure using the non-linear pole-zero law [18] with the parameter values listed in Table 2 taken from Remme et al. [24]. This parameter set was applied to the regions of the model shown in Fig. 2 to characterize increased stiffness at the apex and base due to higher collagen and valve ring effects, respectively.

Type | Parameter | Apex | Sub-apex | Normal | Base |
---|---|---|---|---|---|

Coefficient | k_{11}, k_{22}, k_{33} | 2.22 kPa | 2.22 kPa | 2.22 kPa | 2.22 kPa |

k_{12}, k_{21}, k_{13}, k_{31}, k_{32}, k_{23} | 2.0 kPa | 2.0 kPa | 1.0 kPa | 2.0 kPa | |

Pole | a_{11} | 0.136 | 0.227 | 0.475 | 0.423 |

a_{22} | 0.136 | 0.227 | 0.619 | 0.555 | |

a_{33} | 0.136 | 0.227 | 0.943 | 0.845 | |

a_{12}, a_{21}, a_{13}, a_{31}, a_{23}, a_{32} | 0.3 | 0.4 | 0.8 | 0.4 | |

Curvature | b_{11}, b_{22} | 2.22 | 1.67 | 1.5 | 1.5 |

b_{33} | 2.22 | 1.67 | 0.442 | 0.442 | |

b_{12}, b_{21}, b_{13}, b_{31}, b_{23}, b_{32} | 1.5 | 1.0 | 1.2 | 1.0 |

Spatial location | Homogeneous model (ms) | Heterogeneous model (ms) |

A | 362 | 437 |

B | 360 | 434 |

C | 351 | 418 |

D | 350 | 417 |

E | 344 | 367 |

F | 339 | 367 |

The spatial locations of points A–F are illustrated in Fig. 7a.

Following the work of Tomlinson et al. [25] the initial activation times at points on the left ventricular endocardial surface are specified to simulate Purkinje fibre activation. A stimulus current is then applied at these points following the specified activation profile. The applied stimulus distribution is shown in Fig. 3 where the stimulus times are derived from the initial activation times from the canine ventricular model of Tomlinson et al. [25]. As shown in Fig. 3, the earliest activation occurs in the lower mid-septal wall and the latest activation is toward the base of the free wall.

The boundary conditions for simulating the four phases of the cardiac mechanical cycle were applied based on the temporal relationship between global contraction and electrical transients as described by a classical Wigger's diagram [26]. At all time points throughout the cycle spatially constant pressures were applied over the entire endocardial surface. Prior to the initiation of activation the model was inflated to a cavity of pressure of 1.0 kPa. The isovolumic phase was implemented by adding a finite element mesh which represented the volume of blood in the left ventricle and was constrained such that these elements had a constant total volume. This cavity region was then coupled to the ventricular wall to act as feedback mechanisms to constrain weakly the deformation of the ventricular wall. Once pressure exceeded 10 kPa the ejection phase was initiated by allowing a central node in the cavity mesh to displace in the base–apex axis of the model, thereby reducing cavity volume. Displacement of the cavity node was controlled to match the left ventricular cavity volume transient described by the Wigger's diagram, scaled to give a peak ejection fraction of 50%.

The separate solution algorithms of the grid based finite element method (GBFE) for the mono-domain equations (activation) and finite element method (FEM) of the finite deformations equations (mechanics), have previously been described elsewhere [27,28]. The solution process for coupling these two solution techniques together is summarized in Fig. 4. From the activation solutions solved using a time step of 0.01 ms and an average grid point spacing 200 μm the cellular fields were locally interpolated at the lower resolution gauss points in the mechanics mesh. Using these interpolated values the cell model was solved at the gauss point to determine regional active tension, which was then incorporated into the finite deformation mechanics model. The local strain was calculated at each activation grid point from the deformation and used to update the model cellular lengths. The relatively slower temporal dynamics of myocardial deformation means that this coupling update process could be implemented at a significantly larger time step, 1 ms, than the activation solution process.

## Results

The membrane potential and myocardial deformation of the finite element mesh are shown in Fig. 5a–h at specific time points marked on Fig. 6 through the cardiac cycle to illustrate stages of activation and contraction.

At the eight locations in the finite element mesh shown in Fig. 7a the calculated strain, active tension and membrane potentials are illustrated through the cardiac cycle in Fig. 7b–d.

Using the Purkinje-like stimulus the smooth transmural progression of electrical activation is clear, with full ventricular activation occurring after 70 ms. The repolarization wave begins at the apical epicardium approximately at 270 ms and progresses toward the endocardial surface and basal plane. The mid-wall at the base of the ventricle is the last tissue to repolarize approximately at 490 ms. With the simple transmural fibre variation illustrated in Fig. 1a and the sequence of contraction solutions given in Fig. 7, we show that the left ventricle initially deforms into a more spherical shape with just the endocardial region actively contracting. Then as the mechanical activation moves out through the wall to the circumferentially oriented mid-wall fibres, the deformation becomes more longitudinal and the twisting motion begins. This is followed once more by a move toward more spherical deformation as the inclined epicardial fibres are activated and the twisting increases.

An important feature of this heterogeneous model is the reversal of the direction in which the repolarization wave travels through the wall. During the repolarization phase of the cycle the model maintains a higher cavity pressure due to the prolonged action potentials of the midmyocardial cells, giving rise to increased active tension at the cellular level. The peak pressure achieved by the model during ejection is consistent with values reported in similar models [9].

To interpret the sequence of electrical and mechanical events taking place during the pre-ejection phase of systole, the results shown in Fig. 7 are reproduced on an expanded 200 ms time scale in Fig. 8. The electrical activation sequence is endocardium to epicardium with the free wall endocardial point shown in yellow (B) on the left of the ventricular model in Fig. 7a being activated 50 ms after the septal wall endocardial point shown in grey (A) on the right. The septal transmural points (endocardial to epicardial) shown as grey (A), orange (C) and black (E) are separated by 100 ms intervals (easiest to see in Fig. 8b) as are the (slightly later) transmural sequence of free wall points shown as green (B), light blue (D) and purple (F).

The earliest mechanical response is seen in Fig. 8a is fibre shortening at A, 15 ms after the electrical stimulus at A. The mechanical consequence of the active septal endocardial shortening is stretching of the still passive adjacent mid-wall fibres at C. The next point to activate is the free wall endocardial point B which begins contracting 15 ms later and similarly causes substantial lengthening of mid-wall fibres at point D. The third point to activate is the mid(septal)-wall point C at about 22 ms and 15 ms later its active response — which follows the ‘pre-stretch’ — can be seen in Fig. 8a as fibre shortening. The mid-wall fibres at C and D next begin to shorten (at 35 ms and 45 ms, respectively) and contribute to the very substantial lengthening of epicardial fibres at E and F which, following their passive pre-stretch, begin their contraction about 15 ms after excitation such that by 60 ms all points are activated. Note, however, that the fibre shortening at points A and B is temporarily reversed as the contracting mid-wall and then epicardial fibres exert their influence.

This substantial redistribution of fibre lengths during the early pre-ejection phase of systole may play a significant role in preparing the mid-wall fibres to contract. By 80 ms ventricular pressure has risen sufficiently to begin ejection.

To assess the potential influence of cellular heterogeneity on myocardial deformation, specifically during repolarization, the above simulation was repeated using a homogeneous distribution of cellular models (*τ*_{r}=144 ms and APD=352 ms). Figs. 9 and 10 contrast the model extension ratios and active tension transients in the heterogeneous and homogeneous models and Table 3 provides data on the transmural APD dispersion in the septum and free wall of the model. In the heterogeneous model the sequence of repolarization, as specified by the APD assigned to different regions through the left ventricular wall, is epicardial first (points E, F in Fig. 7a), then endocardial (points A, B), and then mid-wall (points C, D). The initial epicardial repolarization is reflected in the rapid early sarcomere length increase from 300 ms on in the sub-epicardium in Fig. 7b (and in expanded form in Fig. 9 for both homogeneous and heterogeneous models). The sarcomere lengths at the mid-wall points C, D begin increasing at about 360 ms for both the homogeneous and heterogeneous models. The major consequence of assuming homogeneous APDs (Fig. 9b) is a substantially more dispersed transmural sarcomere length variation. For example, at 400 ms fibre extension ratio distribution in the homogeneous model is 0.9 to 0.98, whereas in the heterogeneous model it is reduced to 0.91 to 0.96. Within the sub-epicardium and mid-wall region the range is only 0.91 to 0.92.

## Discussion

The simulation results presented above demonstrate that the modelling and simulation framework developed and implemented works well for coupled electromechanical models using a three-dimensional geometry representative of the cardiac left ventricle. Two specific insights have been provided by this study. The first is the important role that electrical heterogeneity may play in reducing the transmural sarcomere length variation during repolarization. The implications for this effect, including distribution of mechanical work and the metabolic efficiency of the whole heart, provide interesting avenues for future investigation via an electromechanical model. The second insight provided by the model is into the detailed relationship between activation times, shortening and spatial location within the myocardium. In particular, the possibility of endocardial stretch during late isovolumic contraction has not previously been reported. The results in Fig. 8a show a significant heterogeneity of strain at end isovolumic contraction which is consistent with the model results of Kerckhoffs et al. [8]. However, this result is inconsistent with experimental data that report a relatively homogeneous distribution of strain and work throughout the myocardium [29]. Thus, possible influences such as geometry, complex fibre angle variation and a spatial distribution of delay between excitation and activation may ultimately need to be included in the model. It is also important to note that the single beat simulated in this study is effectively modelling a very slow steady-state heart rate where the cells prior to stimulation are always at, or close to, the initial condition set in this model. This may explain the differences between model and experimental results. However, to characterize more effectively, these memory effects where ionic concentrations and gating variables at stimulation are dependent on heart rate, a more biophysical basis to the cell model is required.

The most significant drawback in the framework highlighted by these simulations is the sheer amount of time taken to obtain these results, the coupled electromechanics simulations shown in Figs. 6 and 7 took in the order of three weeks to achieve 1000 ms of simulation time using an IBM Regatta P690 high-performance computer implemented in parallel on eight processors. Furthermore, despite the large computational load imposed by high spatial resolution grids, numerical artefacts remain. A small error associated with grid discretization is shown in the brief hyperpolarization prior to upstroke. Also, the sensitivity of the length–active tension coupling produces rapid oscillations in calculated active tension during contraction induced by finite time stepping in the numerical method. However, on the time scale of myocardial mechanics these variations have little effect on ventricular deformation although ultimately physiological based viscous damping will need to be added completely to eliminate the oscillations.

It is this computational limitation which currently also precludes the inclusion of a more detailed cellular model. Such a model could potentially account for mechano-electric feedback and stretch dependence, which can now be explicitly linked to the tissue simulations via the strong coupling solution procedure outlined in Fig. 4. It is this coupling between cellular contraction and activation that Kerckoffs et al. [8] propose as a possible mechanism to produce a heterogeneous electromechanical delay that is required simultaneously to produce both physiological activation and contraction spatio-temporal sequences.

The framework is sufficiently general, such that new cell models can be accommodated without changing the existing software tools. Thus, as computational resources continue to increase rapidly, we anticipate that in the near future we will be using this framework with more detailed biophyisically based cell models. This generality has been implemented using the CellML language [30,31] to embed the cellular models in the continuum simulations. A general application programme interface (API) has been developed (freely available at http://cellml.sourceforge.net) to provide a layer of abstraction between an application and the implementation of a CellML processing library. Using the API a list of the cellular equations was generated for the model and a specific maths writer was developed to create a Fortran subroutine from the cellular equations which could then be compiled and linked with our CMISS simulation environment (freely available for academic use at http://www.cmiss.org). The CellML model repository (http://www.cellml.org/examples/repository/index.html) provides an extensive database of cellular models of cardiac electrophysiology and mechanics. Thus, it is largely only computational resources that prevent repeating the above study with a detailed biophysically based cellular model.

To assess the computational requirements of a more detailed cellular model we have run electromechanics simulations using a 4×4×4 mm cube of cardiac tissue using a 200 μm grid point spacing. With a reduced number of grid points it was possible to compare the computational performance with the simplified model FK-HMT model used above with a biophysically based coupled LRd-HMT model [4,18]. Profiling of these simulations showed the additional biophysical detail produced a threefold increase in memory usage and a tenfold increase in computational time.

Thus, applying the often quoted Moore's Law, a strongly coupled cardiac electromechanics simulation using a detailed cell model will be likely to be feasible within the next three years. Future developments within that timeframe also include solving electromechanical models using more realistic anatomical and microstructural models [32,33] and coupling three-dimensional ventricular fluid dynamics solutions to determine the endocardial boundary conditions.

Within the wider scope of the Physiome Project we aim to enhance further the modelling framework with the use of markup languages such as CellML (as discussed above) and other ontologies under development. For example, a markup language to represent the spatial variation of fields used in the models (FIELDML: http://www.physiome.org.nz) is also being developed. In doing so we aim to improve the accessibility of models across spatial scales and accelerate the development, validation and application of mathematical models of the heart and other organs.

## Acknowledgments

P.J.H. and D.P.N. acknowledge the support of the Centre for Molecular Biodiscovery. N.P.S. acknowledges support from the Marsden Fund and New Zealand Institute of Mathematics and its Applications. All authors gratefully acknowledge contributions by colleagues in the Auckland Bioengineering Institute.

- © 2005 The European Society of Cardiology. Published by Elsevier Ltd. All rights reserved.