New developments in the auckland heart model

P.J. Hunter and A.J. Pullan
Bioengineering Institute, University of Auckland, New Zealand

ABSTRACT

An anatomically accurate mathematical model of the right and left ventricular myocardium is described, based on measurements of the geometry and fibrous-sheet structure of the pig heart. The model is used for mechanical analysis of ventricular filling with Galerkin finite element techniques and an orthotropic constitutive law based on the fibrous-sheet structure. Bidomain electrical activation of the myocardium is also modeled, with ionic current based electrophysiological equations and reaction-diffusion equations based on orthotropic conductivity tensors referred to the fibrous-sheet material axes. The heart model is also solved with the eikonal wavefront equation and coupled into a torso model to predict the potential distributions on the body surface. 

INTRODUCTION

Computational modeling of the heart has the potential to link molecular processes at the genomic and proteomic level, through cellular and tissue level models to integrated physiological processes at the whole organ and whole body levels. To achieve this ambitious goal the models must be anatomically accurate, so that the detailed tissue structure can be incorporated, be biophysically detailed, so that tissue material properties can be incorporated, and be linked down to cell electrophysiology, signal transduction pathways and metabolic pathways. In this paper we describe the Auckland heart and torso models which are being developed as part of a world wide “cardiome project” effort [1,2] to use computational modeling of the heart to link cell and tissue behaviour to intact organ and whole body physiological function.

VENTRICULAR ANATOMY

Our initial finite element model of the geometry and fibrous-sheet structure of the left and right ventricles of the heart was based on measurements from dog hearts and used prolate spheroidal coordinates [3]. Fig.1 shows the geometry and orientation of fibres around the canine heart at epicaridial, midwall and endocardial locations.

        (a)                                 (b)                                      (c)

Figure 1. Muscle fibre orientations at three transmural depths: (a) epicardial, (b) midwall, and (c) endocardial. The fibre angle at the LV free wall (on the right in these figures) varies continuously from about –60° to the horizontal plane at the epicardium to +80° at the endocardium. The RV outflow tract is at top left in these figures.

Myocardial sheets are formed by layers of tightly coupled myocytes about 4 cells thick [4], which facilitate wall shear and also influence electrical wave propagation (see below).

More recently a new Auckland heart model has been developed in rectangular cartesian coordinates with geometry and structural measurements from pig hearts. The fibre orientations are very similar to the previous dog heart measurements but the sheets are more radially oriented. More accurate anatomical measurement in the basal and apical regions of these hearts has resulted in: (i) a better representation of the fibrous skeleton supporting the heart valves (see Fig.2a); (ii) a more accurate parameterisation of the very thin wall at the apex; and (iii) the non-zero imbrication angle at the base and apex which allows the fibres to spiral inward from epicardium to endocardium. 88 elements with tricubic Hermite basis functions are used to capture the myocardial anatomy to within the geometric measurement error of about 0.5mm.

 

                                                        (a)                                                  (b)

Figure 2. (a) The epicardial surface of the 88 element finite element model of the right and left ventricles. The fibre orientations are shown as vectors on the surface.  (b) Principal strain vectors shown at end-diastole at midwall points of the finite element mesh. The transmurally directed vectors are compressive.

MECHANICS

An orthotropic “pole-zero” constitutive law has been developed for modeling the behaviour of the passive myocardium in which strain energy terms are based on deformation relative to the fibre, sheet and sheet-normal material axes [5]. The HMT model for active tension generation [6], based on data from cardiac trabeculae, is driven by Ca-binding to troponin-C (see below). The stress equilibrium equations of mechanics are approximated with Galerkin finite element equations, using a hydrostatic pressure variable as Lagrange multiplier for the incompressible constraint, and solved with Newton-Raphson and sparse matrix linear solvers [7]. Fig.2b shows the principal strains at end-diastole (1kPa left ventricular pressure). This solution took 1 hour on a single processor of an IBM EServer p690 using 4 load steps with 5 Newton-Raphson iterations to achieve convergence at each load step.      

CELLULAR MODELS

Cardiac cell electrical activity is based on one of a number of published electrophysiological models, such as the Luo-Rudy II [8] and Noble98 [9] models. These models and many more are encapsulated in the CellML markup language at www.cellml.org, a public domain repository of cellular models being maintained by the Auckland Bioengineering Institute and Physiome Sciences Inc. Models of coupled cellular electromechanics are described in [10] and include stretch-activated ion channels [11]. 

ELECTRICAL ACTIVITY OF THE MYOCARDIUM

The model for the spread of electrical activity in the specialised conduction pathways and myocardium currently includes a (fairly crude) geometric representation of the Purkinje fibres. The ends of the Purkinje network are the activation points for the myocardial model as shown in Fig.3.

Figure 3. Activation of the septum from the Purkinje fibres. The second and third figures are 5ms and 10ms, respectively, after the first. The darker shade indicates the active tissue.

The bidomain equations underlying myocardial activation and extracellular current flow are formulated with a diffusion tensor based on the initially orthotropic fibrous-sheet material axes [12]. Ion channel density (reflected in the open channel conductivity) for I_to, I_K1, I_Kr and I_Ks can be varied transmurally to reflect the longer action potential duration (APD) for subepicardial cells and the very long action potential duration seen in M-cells under slow pacing [13]. The computational grid point mesh for solving the activation equations is generated at material points defined by the tri-cubic Hermite finite element mesh described above. Typically there are 60x60x60 grid points per element, giving a spatial resolution of 0.2mm and about 4.10⁷ computational points for a 0.3kg heart.

We are currently developing multigrid schemes for solving the bidomain equations which only use the fine spatial resolution when the spatial gradients of potential are large – in the vicinity of the moving wavefronts which typically occupy only 0.1% of the myocardial volume. Without this grid adaptivity our estimate of the time taken to solve a complete ventricular sequence on 8 processors of an IBM EServer p690 machine with 32GB of RAM is at least 30 hours.

EIKONAL EQUATION MODELING

Incorporating cellular models into a whole-heart framework generates a computational problem that is large, even by today’s standards. However, if details of the action potential upstroke are simplified to that of a step function, then it is possible to simulate the electrical activity of the ventricles with rather modest computational resource.  Such an eikonal approach reduces the problem to finding the location of the activating wave, at the expense of ignoring such features as repolarisation and the variation of the electrical properties of the different cell types within the myocardium. Despite its limitations, the eikonal-based solution can still offer insights into the electrical activity of the heart. Eikonal-based solutions on the Auckland heart model are illustrated (for the canine heart) in Fig. 4 [14].

Figure 4. Wavefront locations using an eikonal equation to simulate propagation from the distal ends of the Purkinje tree. For each sample time anterior (top) and posterior (bottom) views are given.

COUPLING TO THE TORSO

An anatomically detailed model of the torso has been developed in order to couple cardiac electrical activity to the body surface [15]. Physically the cardiac electrical activity is connected to current flow in the torso via the extracellular potential field. Treating the combined current flow in the heart and torso as a single problem requires an estimated 100GB of storage at a converged resolution [16]. To work around these computational demands the problem is split into two parts: the full bidomain simulation on the heart is used to derive a lower resolution multiple dipole source model to produce body surface signals. The resultant signals are quantitatively different from those obtained in a fully combined approach, due to the limited resolution offered by the dipoles, but qualitatively the pattern is the same. Results from such an approach are given in Fig. 5 [17].

Figure 5. Anterior images of torso potentials generated from a  bidomain solution using the Noble 98 electrophysiological cell model. Images are drawn at 15, 35, 55, 75, 95, 405, 420 and 435 ms after the initial endocardial stimulus. The first 5 images correspond to the QRS interval and the remaining 3 are for the T wave. A full bidomain model (with unequal anisotropy) is used and contains over 250,000 grid points which equates to about 1.2mm spacing. The time step is 0.01ms and at each 10 time steps, 40 equivalent moving dipoles are calculated (one for each of the 40 heart elements). These dipoles are then placed into the torso model containing 1120 cardiac elements, 2 lungs each with 320 elements and the main torso volume is described using 2816 elements. The reference is on the anterior surface and the lower right of the torso. The heart is activated at 3 sites at different times to try to get a reasonable activation sequence. The lightly shaded patches are negative with a maximum negative potential of -0.456mV and the darker patches indicate a positive potential with a maximum potential of 0.252mV.

MODELING OTHER CARDIAC PROCESSES

An anatomically realistic model of coronary blood flow has been developed in which the Navier-Stokes equations governing arterial and venous flow have been reduced to one dimension [19]. These flow equations are then combined with a pressure-area relationship which accounts for the compliant vessel wall and time varying compressive force produced by myocardial contraction on these embedded vessels [20]. Fig.6 demonstrates the change in blood flow velocities through the heart cycle, a central result from this modelling study. The next step is to use the predicted distribution of coronary blood flow with a recently developed cellular model of cardiac metabolism [2] such that the supply of oxygen determines the nature of energy production and the resulting effects on myocyte production.

Figure 6. Coronary arteries used for computation of blood flow, pressure and oxygen transport. Calculated coronary arterial blood flow velocities are shown at the end of the four phases of the cardiac cycle. The darker shade indicates the coronary flow with a maximum value of 250mm/s.

APPLICATIONS

The anatomically and biophysically based models of the heart described here are being used for a number of clinical applications. For example, the time-varying distributions of stress and strain throughout the myocardium can be determined for an individual patient, using tagged MRI and measurement of ventricular pressure, and used for diagnosis of aberrant wall motion.

Model-based diagnosis of cardiac arrhythmias from body surface measurements and applications to drug discovery will be developed once the models include sufficient detail on signal transduction pathways and metabolic pathways in addition to the current level of electrophysiological detail. For example, class III anti-arrhythmic drugs such as amiodarone, which prolong the action potential plateau by inhibiting I_Kr or I_Ks, can give rise to ‘Early After Depolarisations’ (EADs), especially under conditions of hypokalaemia, through reactivation of the I_CaL channels [13]. Due to the longer time available for inactivation processes to recover, these channels can reactivate during the plateau. This situation is exacerbated in exercise when b-receptor activation enhances I_CaL.  The question of whether these EADs can cause Torsade de Pointe arrhythmia is still debated since the electrotonic load in 3D tissue may prevent propagation of the EADs. This is a prime example of the need for anatomically and biophysically detailed heart models.

Another example is the use of patient-specific models to aid in the understanding of multiple ECG signals - the so-called inverse problem of electrocardiography.  Experiments are currently underway [18] to validate the model-based predications.

The heart modeling described above is contributing to the international cardiome project being conducted in collaboration with many of the leading cardiac groups around the world. Ultimately an accurate beating heart is envisioned in which patient-specific profiling (genetic as well as anatomical data) may be incorporated, opening the way for personalised medicine based around model-based predictions.

ACKNOWLEDGEMENTS

We would like to acknowledge the contributions from many members of the Bioengineering Institute to the work presented here.  In particular, much of the work reported here is based on the recent doctoral research of Carey Stevens, Martin Buist, Karl Tomlinson and Nic Smith. The graphical environment for displaying and interacting with the 3D models is the work of Dr David Bullivant, Dr Richard Christie and Shane Blackett.

REFERENCES

[1] P. Kohl, D. Noble, R.L. Winslow and P.J. Hunter. “Computaitonal modelling of biological systems: tools and visions.” Phil. Trans. Roy. Soc. Lond. A358, pp.579-610, 2000.

[2] N.P. Smith, P.J. Mulquiney, M.P. Nash, C.P. Bradley, D.P. Nickerson and P.J. Hunter. “Mathematical modelling of the heart: cell to organ.” Chaos Sol. Frac. Vol.13, pp.1613-1621, 2002

[3] I.J. LeGrice, P.J. Hunter and B.H. Smaill. “Laminar structure of the heart: a mathematical model.” J. Physiol. Vol.272, pp.H2466-H2476, 1997.

[4] I.J. LeGrice, P.J. Hunter, A.A. Young and B.H. Smaill. “The architecture of the heart: a data-based model.” Phil. Trans. Roy. Soc. Vol. 359, pp.1217-1232, 2001.

[5] M.P. Nash and P.J. Hunter. “Computational Mechanics Of The Heart.” J. Elasticity Vol.61(1-3), pp.113-141, 2001.

[6] P.J. Hunter, A.D. McCulloch and H.E.D.J. ter Keurs. “Modeling the mechanical properties of cardiac muscle.” Prog. Biophys. Molec. Biol.  Vol.69, pp.289-331, 1998.

[7] K.D Costa, P.J. Hunter, J.S. Wayne, L.K. Waldman, J.M. Guccione and A.D. McCulloch. “A three-dimensional finite element method for large elastic deformations of ventricular myocardium: Part II - Prolate spherical coordinates.” J. Biomech. Eng. Vol.118, pp.464-472, 1996.

[8] C. Luo and Y. Rudy. “A Dynamic model of the cardiac ventricular action potential - simulations of ionic currents and concentration changes.” Circ. Res., Vol.74, pp.1071-1097, 1994.

[9] D. Noble, A. Varghese, P. Kohl and P. Noble. “Improved guinea-pig ventricular cell model incorporating a diadic space, IKr and IKs, and length- and tension-dependent processes.” Can. J. Cardiol., Vol.14, pp.123-134, 1998.

[10] D.P. Nickerson, N.P. Smith and P.J. Hunter. “A model of cardiac cellular electromechanics.” Phil. Trans. Roy. Soc. Vol. 359, pp.1159-1172, 2001.

[11] P. Kohl, P.J. Hunter and D. Noble. “Stretch-induced changes in heart rate and rhythm: clinical observations, experiments and mathematical models.” Prog. Biophys. Molec. Biol. Vol.71(1), pp.91-138, 1999.

[12] P.J. Hunter and B.H. Smaill. “Electromechanics of the heart based on anatomical models.” In Cardiac Electrophysiology: from cell to bedside. 3rd Edn.. Eds Zipes and Jalife. Saunders. Vol.32, pp.277-283, 2000.

[13] C. Antzelovitch, V.V. Nesterenko, A.L. Muzikant, J.J. Rice, G. Chien and T. Colatsky. “Influence of transmural gradients on the electrophysiology and pharmacology of ventricular myocardium. Cellular basis for the Brugada and long-QT syndromes.” Phil. Trans. Roy. Soc. Lond. A359, pp.1201-1216, 2001.

[14] K.A. Tomlinson, P.J. Hunter and A.J. Pullan. “A finite element method for an eikonal equation model of myocardial excitation wavefront propagation.” SIAM J. Appl. Math. (In press), 2002.  

[15] C.P. Bradley, A.J. Pullan and P.J. Hunter. “Geometric modeling of the human torso using cubic Hermite elements.” Ann. Biomed. Eng.  Vol.25, pp.96-111, 1997.

[16] M.L. Buist and A.J. Pullan. “From cell to body surface: a fully coupled approach.” J. Electrocardiology, Vol.34, pp.191-195, 2001.

[17] M.L. Buist. Modelling cardiac activity from cell to body surface, PhD Thesis, University of Auckland, 2001.

[18] M.P. Nash, C.P. Bradley, A. Kardos, A.J. Pullan, D.J. Paterson. Chaos Sol. Frac. Vol.13(8), pp.1735-1742, 2002.

[19] N.P. Smith, A.J. Pullan and P.J. Hunter. “Generation of an anatomically based geometric coronary model.” Ann. Biomed. Eng. Vol. 28(1), pp.14-25, 2000.

[20] N.P. Smith, A.J. Pullan and P.J. Hunter. “An anatomically based model of transient coronary blood flow in the heart.” SIAM J. Appl. Math., Vol.62, No.3, pp.990 –1018, 2001.