Introduction

First idealized models describing the normal activation sequence in the human heart were reported over two decades ago [1, 2]. Mainly due to computational limitations, the models did not include myocardial anisotropy nor physiological propagation. On the basis of anisotropic bidomain theory [3] and cellular automata theory [4], development of more realistic whole-heart models have become feasible [5].

A ventricular model that produces a correct normal activation sequence is a prerequisite for simulating pathological conditions, such as ischemia, infarction or ventricular arrhythmias. A comprehensive model should feature an anatomically accurate geometry, intramural fibrous structure, and a conduction system. Implementation of the ventricular conduction system is a challenge, since it should be flexible enough to allow "rewiring" of the conduction system for each individual case where endocardial mapping data are available. Advances in computer technology have only recently made this approach possible.

Detailed studies on the anatomy of the human conduction system are scarce, and often report a high degree of interindividual variability. Therefore, we followed the criterion that the conduction system should rather accurately reproduce correct the pattern of ventricular activation sequence than to follow any predetermined anatomical design. The conduction system model was validated, first by comparing the simulated activation isochrones to those reported in the literature, and second by comparing the resulting body surface potential maps and magnetocardiographic distributions to our recordings obtained in normal subjects. 

Methods

Anisotropic propagation model

Our hybrid model of the anisotropic myocardium describes the subthreshold behaviour of the elements according to the bidomain theory [3], while in the suprathreshold region the elements behave as cellular automata [4]. Each element is assigned a specific type and a vector of local fiber direction [5, 6]. Each element can be in one of four macrostates, corresponding to different phases of an action potential: resting, excitatory, absolute refractory, or relative refractory. During simulation, elements undergo a series of macrostate transitions, which have been described in detail in Ref. [7]. The subthreshold transmembrane potential, v_m is calculated from the generalized cable equation, derived from the anisotropic bidomain theory [3]:

Here we assume that the assumption of equal anisotropy ratios is valid [5, 7]. In Eq. 1, c_m is the specific membrane capacitance, D_i is the intracellular conductivity tensor, and i_ion and i_app are, respectively, ionic and applied currents. Above a threshold value, the cell enters the absolute refractory state and v_m is estimated from a predefined logarithmic function [7]. Physiological parameters for the model are adopted from literature [8, 9].

Extracardiac field calculation

In the forward computations, the ventricles are assumed to be embedded in a torso-shaped piecewise homogeneous volume conductor, including the lungs and the intracardiac blood masses. The infinite medium potential, f_o, in the extracardiac region is determined from the discretized equation

where the integrals are evaluated over the ventricular volume, v_m is as in Eq. (1), a is the local direction of the fiber axis, s is the extracellular conductivity (s = 2.0 mS/cm), s₁ is the transverse conductivity (s₁ = 0.5 mS/cm) and s₂ is the difference between the axial and transverse conductivity (s₂ = 1.5 mS/cm) [7]. The first term of the potential in Eq. (2) represents the contribution of the isotropic component and the second term accounts for the anisotropic properties of the cardiac tissue.

The extracardiac magnetic field, H_o, is evaluated from a corresponding equation,

To compute the body surface potential and magnetic field in the inhomogeneous torso, a "fast forward solution" is used [7]:
 

In Eq. (4), W is a geometry matrix consisting of the solid angles subtended by a surface element at each node, and A is another (vector-value) geometry matrix taking into account the influence of the conductivity change surfaces on the magnetic field (for details, see Ref. [10]).

Fiber arrangement

The model consists of 2.000.000 excitable elements on a cubic lattice with 0.5 mm spacing [11]. The geometry was reconstructed from 1 mm sections of the human heart [1]. The assignment of the principal fiber direction was performed separately for left-ventricular, right-ventricular and papillary-muscle cells, with the fiber direction rotating from endocardium to epicardium [12]. Figures 2 and 3 illustrate the amount of anatomical detail that was included in the model.

Figure 1. Fiber arrangement on an yz-cross section approximately in the middle of the ventricles. The principal fiber directions are marked with arrows. The angles, measured from the local tangent vector at each element, rotate from about -50 degrees on the endocardium to about +45 degrees on the epicardium.

Figure 2 .Three-dimensional presentation of the ventricles as seen from the base of the heart.

Results

Simulations in slab geometry

In order to test the propagation algorithm and difference approximations to calculate Eq. (1), we first run simulations in a three-dimensional slab (100 x 100 x 100 cells, h = 0.05 cm). The fiber direction was kept constant in each simulation. Propagated isochrones for three cases are shown in Figure 3. They agree very well with theoretical predictions, i.e., the ratio of the longitudinal and transversal axes of an isochrone are proportional to the square-root of the conductivities in the corresponding directions. The propagation velocities defined from the isochrones yield 0.69 m/s in the axial direction, and 0.31 m/s in the transversal direction.

Figure 3. Propagation in slab geometry. The fiber direction in each case was constant, the angle from the x-axis being 0, 45 and 65 degrees from left to right, correspondingly. Each shaded zone corresponds to 5 ms. The ellipsoid plotted on top of the isochrones denotes the theoretically defined 25-ms isocrone.

Conduction system

We have created a computer interface for the interactive editing of the conduction system, utilizing the OpenGL capabilities of a SGI Visual PC. The organization of the conduction system was created by using information available from textbooks on human anatomy and other literature on the topic [13, 14]. The geometry of the conduction system is shown in Figure 4. The His bundle consists of a single branch in the right ventricle, whereas it resembles a fan-like sheet of fibers in the left ventricle. The His bundle continues on both sides as the Purkinje network that contains the Purkinje-myocardial junction (PMJ) sites. On the right, the most prominent feature of the conduction system is the single bundle that carries the activation from the septum to the free wall and the papillary muscles along the moderator band. On the left, there are three major areas of activation: the septum and the inferior and superior free wall. There is no conduction network in the posterior free wall.

Figure 4. The geometry of the conduction system observed from different angles. The right branch of the His bundle (not shown) connects the topmost nodes of the left and right conduction system. In the RV apical view (a), notice how the moderator band contains the conduction system activating the RV free wall and the papillary muscles (arrow). The RV free wall (b) and septal (c) views show the right conduction system, consisting of one main bundle with several smaller branches.
The views of LV septum (d) and LV free wall (e) show the more complicated structure of the left conduction system. Notice the fan-like structure of the septal conduction system in (d), and the vertical spiral gap (dashed area) in the conduction system in (e). The PMJs are shown as light spheres.

Propagation in the 3D ventricles

The simulated activation sequence resulting caused by the activation of the conduction system is shown in Figure 5. The result agrees well with isochrones obtained from an isolated human heart [15]. The ventricular activation starts in the left ventricular septum (layer 110), matched by a right ventricular septal activation 20 ms later (layers 70-90). Almost simultaneously with the RV septal activation, the inferior (in body coordinates) and anterior LV activation appear (layers 90-110). These are followed by the activation of the RV free wall (layers 90-130). The RV and LV breakthroughs take place at 30 ms and 45 ms, respectively. In the final stages of the QRS, activation propagates through the posterior LV free wall and the pulmonary conus.

Figure 5. Simulated activation isochrones. The layers are with a spacing of 10 mm, and correspond to those presented by Durrer et al. [15]. The isochrone spacing is 5 ms.

Simulated BSPM and MCG patterns

The geometry in the extracardiac field computations is depicted in Fig. 6. Figs. 7 and 8 show, respectively, simulated BSPM and MCG distributions during the QRS complex. In Fig. 7, the initial maximum resulting from left septal activation is anterior. Then the minimum on the back moves upward and travels over the right shoulder onto the right anterosuperior region, indicating apical activation and masking of the left septal activation by the corresponding right-septal activation. The area of positive potentials then drifts to the back. In the end, the activation propagates toward the pulmonary conus, reflecting as a superior positive potential in BSPM. The sequence of MCG patterns in Fig. 8 is also in good agreement with the behaviour in measured MCG maps.

Figure 6. The geometry in extracardiac BSPM and MCG computations. BSPMs were evaluated on the nodes of a triangulated torso model of an adult male, while the normal component of the cardiomagnetic field (MCG) was computed on a plane array in front of the chest.

Figure 7. Simulated BSPM distributions at 10 ms steps during the QRS complex. The read and blue areas denote, respectively, positive and negative values.

Figure 8. Simulated MCG distributions at 10 ms steps during the QRS complex. The red areas denote positive values, where the magnetic flux enters the chest, while in blue areas the flux comes out of the chest.

Conclusions

We created an interactive method for implementing the ventricular conduction system in anatomically accurate model of the human ventricles. The conduction system model was validated by comparing simulated isochrones, body surface potential maps and magnetocardiograms to existing data. Our validated conduction network represents an advancement of the forward solution of electrocardiography and magnetocardiography. It will be indispensable for simulating ventricular arrhythmias involving the conduction system, such as Purkinje-muscle re-entrant arrhythmias. The flexibility of our method of conduction system generation will allow modifications required for adapting our model to match measured data obtained during endocardial mapping. 

References