A Simplified 3D Model of Whole Heart Electrical Activity and 12-Lead ECG Generation

We present a computationally efficient three-dimensional bidomain model of torso-embedded whole heart electrical activity, with spontaneous initiation of activation in the sinoatrial node, incorporating a specialized conduction system with heterogeneous action potential morphologies throughout the heart. The simplified geometry incorporates the whole heart as a volume source, with heart cavities, lungs, and torso as passive volume conductors. We placed four surface electrodes at the limbs of the torso: V R, V L, V F and V GND and six electrodes on the chest to simulate the Einthoven, Goldberger-augmented and precordial leads of a standard 12-lead system. By placing additional seven electrodes at the appropriate torso positions, we were also able to calculate the vectorcardiogram of the Frank lead system. Themodel was able to simulate realistic electrocardiogram (ECG) morphologies for the 12 standard leads, orthogonal X, Y, and Z leads, as well as the vectorcardiogram under normal and pathological heart states. Thus, simplified and easy replicable 3D cardiac bidomain model offers a compromise between computational load and model complexity and can be used as an investigative tool to adjust cell, tissue, and whole heart properties, such as setting ischemic lesions or regions of myocardial infarction, to readily investigate their effects on whole ECG morphology.


Introduction
Biophysically based mathematical models of whole-heart electrical activity are becoming increasingly detailed and complex, with high-resolution anatomically accurate models requiring extensive computation times, dedicated software, and even the use of supercomputers [1][2][3]. We have developed a simplified, computationally highly efficient threedimensional (3D) torso-embedded whole heart model, capable of reproducing realistic 12-lead surface electrocardiograms (ECGs) on the torso. The model incorporates spontaneous activation in the natural pacemaker of the heart, the sinoatrial node (SAN), as well as the His-Purkinje specialized conduction system. Also included are heterogeneous action potential (AP) morphologies in various regions of the heart, whose parameters such as duration and amplitude can be readily adjusted to simulate their influence on ECG waveforms.
This simplified whole heart model can be utilized as a starting point for the rapid a priori formulation, testing, and refinement of hypothesis on the forward and inverse relations between the surface ECG and cardiac cell or tissue properties, including action potential duration, refractory period, upstroke velocity, amplitude, conduction velocity, and tissue conductivity in various subregions of the heart. This will allow the ready visualization of the relation between specific myocardial electrical properties and the resulting body surface potentials, expressed as a standard 12-lead ECG system, a vectorcardiogram (VCG), an electrogram (EGM), a body surface potential map (BSPM), or any other lead system. Such a model could readily be utilized to simulate regions of the heart with reduced excitability and/or conduction disturbance, for the development of new diagnostic tools able to inversely assess key cardiac properties or localization of ischemic or infarcted regions by processing features of the ECG.
It has been suggested that it may not be necessary to use overly complex cardiac models to gain insights into cardiac electrical activity and ECG characteristics; simplified models can be developed first and progressively expanded to more complex and detailed models [4].
To date, there has been a lack of a general purpose, easy replicable, and simplified 2D or 3D model of the torsoembedded heart [5,6], perhaps because simplified models are considered unable to generate realistic ECG morphologies, despite the fact that such simplified models have been mainly used for cardiac tissue modeling [1,2]. However, simplified models are easy to replicate, offering a good compromise between computational load and the electrophysiological complexity. If possible, they should be considered as an important stepping stone towards anatomically and functionally more detailed computer models of cardiac electrical activity [4,7].

Methods and Materials
The cardiac torso bidomain model developed in this study comprised a simplified 3D description of the torso, lungs, and the whole heart including atria, ventricles, and blood chambers, as shown in Figure 1. The torso and lungs were extruded from ellipses, with a cavity inserted into the lung domains for positioning the heart. Dimensions of all shapes were approximations to real human anatomy from the Visible Human Project gallery. All model domains have an assigned extracellular conductivity to account for volume conductor effects (Table 1), with values assigned from the literature [8,9].
Modelling of the standard Einthoven and Goldbergeraugmented leads was achieved using four ECG electrodes at the vertices of the torso corresponding to positions of the left arm ( ), right arm ( ), left foot ( ), and right leg ( GND ), as shown in Figure 2(a). Additional six electrodes ( 1 , 2 , 3 , 4 , 5 , and 6 ) were positioned at the chest to model the precordial leads, completing all channels of a 12lead ECG system. Furthermore, seven additional electrodes ( , , , , , , and ) were placed on the surface of the chest (Figure 2(b)) to determine the VCG using a Frank lead system. The right leg electrode ( GND ) was set as the ground reference, so the potential at that node was fixed at 0 V. The governing equation for the extracellular voltage in the passive volume conductor regions (excluding the myocardium) was given by the Laplace formulation where is the electrical conductivity of respective outside heart subdomain: torso, lungs, and cardiac blood chambers, with values given in Table 1. All exterior boundaries of the torso were set to be electrically insulating (zero normal component of current density), and all the interior boundaries in contact with the heart were set to = where is the extracellular voltage in the myocardial walls.
The three Einthoven leads I , II , and III from the standard 12-lead ECG were determined from The other leads were obtained directly from electrodes on the torso from the six precordial leads 1 , 2 , 3 , 4 , 5 , and  transmembrane potential (TMP) at the heart surface (a) and in a frontal plane cross section midway through the heart (b), when the depolarization wave front first excites the left and right ventricles and when atrial repolarization has just begun. 6 (Figure 2(a)) or by implementing additional calculations: the three Goldberger-augmented leads were calculated from while the precordial leads were calculated as differences of each precordial electrode from the Wilson central terminal = ( + + )/3. Furthermore, the orthogonal leads , , and for the VCG were determined from the Frank lead system by placing additional six electrodes , , , , , and at the anterior surface of the torso and additional electrode at the posterior (Figure 2(b)). The orthogonal leads were computed as in accordance with existing methods [10,11]. Moreover, with an additional calculation, the root mean square curve can be derived as RMS = √( 2 I + 2 II + 2 III )/3, and a similar equation can be used for the signal averaged ECG (SAECG), The heart itself was divided into seven subdomains or regions with heterogeneous cardiac cell properties and tissue conductivities representing specialized cells of the conduction system and the myocardium (Figure 3). An electrical isolation gap exists between the atria and ventricles except at a junction in the septum that links the atrioventricular node (AVN) with the His bundle.
The bidomain model of cellular activation of the heart, including the SAN, was defined at the cellular level by three dependent variables: : the extracellular potential, : the intracellular potential, and : a recovery variable governing cellular refractoriness. The bidomain equations were based on modified FitzHugh-Nagumo equations [12,13]. For each region of the heart, they were defined according to with , denoting the extracellular and intracellular conductivities, respectively, = − , and , , 1 , 2 , , , , , and are region-specific parameters, while ion is defined according to within the SAN and within the walls of the atria, ventricles, AVN, His bundle, bundle branches, and Purkinje fibers. Parameters of the model region-specific values are listed in Table 2. The parameter mainly regulates action potential duration (APD), whilst conductivity parameters , control conduction velocities in the tissue which are well known from the literature [9]. For example, a lower conductivity in the AVN performs an appropriate delay in impulse conduction to the ventricles. The initial values of all model variables are given in Table 3.
Boundary conditions on all interior boundaries in contact with the torso, lungs, and cardiac cavities are zero flux for ; therefore, −n ⋅ Γ = 0 where n is the unit outward normal vector on the boundary and Γ is the flux vector through  that boundary for the intracellular voltage, equal to Γ= − ⋅ / n. For the variable , the inward flux on these boundaries is equal to the outward current density J from the torso/chamber volume conductor; therefore, − ⋅ / n = n ⋅ J.
The simplified 3D cardiac model was simulated using COMSOL Multiphysics (COMSOL AB, Switzerland, v4.2a) finite element software. The resulting finite element mesh consisted of 21106 tetrahedral elements with 51680 degrees of freedom to be solved at each time step. Simulations took approximately 1 hour to solve one second of cardiac activity with 1 ms output time resolution. The simulations were performed on an Intel Core i7-970 PC workstation, with processing power of about 100 Gflops.
In order to validate the model against typical cardiomyopathies encountered in clinical practice, we simulated two myocardial infarcts (MIs) and their resulting effect on the ECG: (a) an anterior MI in the apical section of the heart and (b) an inferior MI. We also examined the possibility of localizing and identifying these two common infarction types with our model. For simulating the infarcts, the intracellular conductivity and parameter were set to zero in the infarcted regions, with initial values for and set to −60 mV and −20 mV, respectively. Both types of infarct were positioned apically, spanning the epicardium, midmyocardium, and endocardium (Figure 7).

Results and Discussion
In all simulations, spontaneous and periodic rhythmic activation occurred in the SAN pacemaker region, located in the right atrial wall of the heart. The electrical activation impulse then spread throughout the atria and through the atrial septum before reaching the AVN, where the excitation wavefront was delayed until the atria were entirely activated. Subsequently, the AVN activated the His bundle from where the activation spread to the bundle branches, the Purkinje fibers, and the whole ventricles ( Figure 4). Figure 5 illustrates the relationship between the Einthoven lead II ECG and the state of excitation in the heart in a frontal plane cross section, in which characteristic timepoints are labeled on the ECG with the matching activation    there is an overlap between ventricular depolarization and atrial repolarization. The peak of the wave (sequence 2) is seen to occur when most of the atria have been depolarized, whilst the QRS complex (sequences 4, 5, and 6) corresponds to the activation of the His and Purkinje regions and the entire depolarization of the ventricles. The peak of the wave (sequence 9) occurs when both ventricles are fully in their repolarization phase. Figure 6 illustrates the time course of the simulated ECG and the transmembrane action potentials (TMPs) at various myocardial locations over an interval of 2 seconds. As noted, the excitation in the SAN is spontaneous, stable, and periodic. The ECG morphology, action potential durations, and activation times throughout the heart agree well with known behavior. Figure 7 illustrates the transmembrane potential at the epicardial surface when the ventricles are completely depolarized for the three cases: (1) normal heart, (2) heart with anterior myocardial infarction (anterior MI), and (3) heart with inferior myocardial infarction (inferior MI). The simulations indicate the absence of activation in both types of infarct. Figure 8 shows the simulated ECG morphology in both lead II and precordial lead 1 for the three cases shown in Figure 7. The lead II waveform reveals ST segment depression (which starts at the point and ends at the beginning of the wave) for both infarct types, with the anterior MI exhibiting deeper ST segment depression than the inferior type (Figure 1, left). The difference between anterior and inferior MI is even more pronounced in precordial lead 1 (Figure 8(b)), since the former leads to ST segment elevation and the latter exhibits depression. These findings of ST segment depression or elevation match their clinical usage as a diagnostic criterion for detection of myocardial infarction [14], whilst MI localization from this criterion is reported to be not very accurate [15]. It is important to emphasize that the exact heart angle as well as the position and orientation of the septum within the torso substantially determines the level of ST segment change. Simplified 3D cardiac models, such as the one described in this study, could be useful as a basis for the refinement and enhancement of ECG-based methods for detection and localization of myocardial infarctions. Figure 9 illustrates the simulated ECG for all 12 leads (Einthoven, Goldberger-augmented, and precordial leads) for the three cases considered above: normal heart, anterior MI, and inferior MI. For most of the ECG leads, the presence of myocardial infarction is disclosed as either ST segment elevation or depression relative to the normal ECG, in agreement with standard clinical findings [14]. Figure 10 illustrates the three simulated orthogonal leads , , and along with reconstructed VCG for the same three cases: normal heart, anterior MI, and inferior MI. Although the ST segment depressions and elevations for the infarction cases are observable in all three orthogonal leads, it is not as easy to clearly identify these slight ST segment differences in the VCG.

Conclusion
We have developed a simplified 3D model of the torsoembedded whole heart capable of generating spontaneous and morphologically realistic 12-lead ECG signals at the torso surface, offering a good compromise between computational efficiency, model complexity, and simulated ECG signal quality.
Two cases were simulated in which myocardial infarctions were imposed at different heart regions, and the impact on ECG morphology was observed. The simulation results showed a change in the ST segment level, confirming wellestablished current clinical diagnostic criteria for identification of myocardial infarction based on the ECG. Our cardiac 3D model offers wide possibilities for modeling various cardiac myopathies and electrical abnormalities, in order to explore and develop new ECG-based diagnostic methods. By providing a good compromise between computational load and model complexity, the model can be used as a stepping stone towards anatomically more detailed models with realistic myocardial architecture, including anisotropic fiber conductivity.
In the future, the simplified 3D cardiac model could even be useful for inverse estimation, whereby cardiac electrical function could be assessed from body surface potentials by estimating model parameters from forward model solutions to obtain meaningful and unique reconstructions of electrical activity at the level of the heart.