Investigating the mechanisms underlying the genesis and conduction of electrical excitation in the atria at physiological and pathological states is of great importance. To provide knowledge concerning the mechanisms of excitation, we constructed a biophysical detailed and anatomically accurate computer model of human atria that incorporates both structural and electrophysiological heterogeneities.
The three-dimensional geometry was extracted from the visible female dataset. The sinoatrial node (SAN) and atrium, including crista terminalis (CT), pectinate muscles (PM), appendages (APG) and Bachmann's bundle (BB) were segmented in this work. Fibre orientation in CT, PM and BB was set to local longitudinal direction.
Descriptions for all used cell types were based on modifications of the Courtemanche et al. model of a human atrial cell. Maximum conductances of , and were modified for PM, CT, APG and atrioventricular ring to reproduce measured action potentials (AP). Pacemaker activity in the human SAN was reproduced by removing , but including , , and gradients of channel conductances as described in previous studies for heterogeneous rabbit SAN.
Anisotropic conduction was computed with a monodomain model using the finite element method. The transversal to longitudinal ratio of conductivity for PM, CT and BB was 1 : 9. Atrial working myocardium (AWM) was set to be isotropic.
Simulation of atrial electrophysiology showed initiation of APs in the SAN centre. The excitation spread afterwards to the periphery near to the region of the CT and preferentially towards the atrioventricular region. The excitation extends over the right atrium along PM. Both CT and PM activated the right AWM. Earliest activation of the left atrium was through BB and excitation spread over to the APG. The conduction velocities were 0.6 m s−1 for AWM, 1.2 m s−1 for CT, 1.6 m s−1 for PM and 1.1 m s−1 for BB at a rate of 63 bpm.
The simulations revealed that bundles form dominant pathways for atrial conduction. The preferential conduction towards CT and along PM is comparable with clinical mapping. Repolarization is more homogeneous than excitation due to the heterogeneous distribution of electrophysiological properties and hence the action potential duration.
Atrial dysfunction is the most common pathology in the heart. Atrial fibrillation (AF), for example, contributes significantly to cardiac mortality and is the single most important cause of ischaemic stroke in people older than 75 years (Nattel 2002). Atrial dysfunction is also accompanied by other severe cardiac diseases, e.g. congestive heart failure. A better understanding of fundamental mechanisms underlying atrial function and dysfunction is of great benefit to therapeutic approaches. Anatomical structure, interatrial coupling, fast conduction bundles and electrophysiological heterogeneity seem to play an important role for atrial excitation conduction during physiological and pathological conditions. To support these hypotheses, we present a biophysical detailed and anatomically accurate computer model of human atria in this work that incorporates both structural and electrophysiological heterogeneities.
The most important anatomical structures for initiation and conduction of atrial excitation are the sinoatrial node (SAN), the crista terminalis (CT), pectinate muscles (PM) and the interatrial connections, e.g. Bachmann's bundle (BB). The left and right atrial working myocardium (AWM) is also important to conduct electrical excitation, but their main function is mechanical contraction. The atrial walls have a thickness of 2–3 mm with bundles being the most dominant structures. Compared to the ventricles, the AWM has less organized fibre orientation. Only the fast conduction pathways in right (RA) and left atrium (LA), the regions around the ostia of the vessels, the limbus of the fossa ovalis, and some regions of the AWM, have some organized fibres (Ho et al. 2002).
Cardiac excitation is first initiated in the SAN, the primary pacemaker. Then, the fast right atrial pathway, i.e. CT and PM conducts the excitation preferentially towards the atrioventricular region. During that phase the RA becomes depolarized. At the same time, BB conducts the excitation towards LA. Further connections between both atria are the limbus of the fossa ovalis and an area at the coronary sinus. The fossa ovalis itself is usually fibrous, being populated by few myocytes and thus mainly non-conducting. Finally the whole atrium is excited after 120 ms (Lemery et al. 2004). All conduction pathways are characterized by a large electrical conductivity along the fibre direction leading to a faster conduction velocity along the main axis compared to AWM. The electrical coupling in SAN is stronger in the periphery.
Lemery et al. (2004) investigated recently the interatrial coupling with simultaneous bi-atrial measurement from humans achieved with non-contact mapping systems. They showed that earliest activation of the LA is during sinus rhythm via BB. Furthermore, the time to activation after SAN activation is 19 ms for BB and 24 ms for the septum in the RA. In the LA the first activation is measured after 41 ms near BB and after 43 ms in septal regions. The endocardial activation times were: 23 ms for BB, 81 ms for RA, 80 ms for LA and 120 ms for the whole atrium.
Understanding atrial electromechanical function requires knowledge of anatomy, electrophysiology and excitation conduction, as well as active and passive mechanics. Models of these functions and properties have been improved in recent decades. In this work, we present a computer model of human atrium including SAN with detailed anatomical and electrophysiological information. The model provides a platform to study how excitation is conducted in the atria during normal physiological conditions. It can also be used to investigate the genesis of abnormal excitation patterns, such as AF. Anti-arrhythmic strategies like specific drug administration or ablation interventions can be investigated in such a model considering electromechanical effects of the treatment.
2. Material and methods
The mathematical description of the excitation conduction in human atrium was based on an anatomical model of the Visible Female dataset, cellular models for describing the electrophysiological properties and a model to calculate the electrical interaction of cells in the whole atrium.
(a) Anatomical model
A foundation of this work was the computer based anatomical model of the whole human heart (Sachse et al. 2000). This model was derived from the Visible Female dataset. The images were taken from the frozen corpse of a 59 year old female and are transversal cryosections with a resolution of 0.33 mm. The distance between the images was 0.33 mm. The images were pre-processed to obtain a three-dimensional dataset. Afterwards, it was segmented and classified using different techniques of digital image processing, e.g. region-growing and interactively deformable contours. Figure 1 depicts the model consisting of cubic volume elements in a surface based visualization with 19 tissue classes.
(b) Electrophysiological model
The employed electrophysiological model of Courtemanche et al. (1998) consists of a set of coupled nonlinear differential equations to characterize the electrophysiological behaviour of human atrial cells. This Courtemanche–Ramirez–Nattel (CRN) model describes the changes of concentration and flow of ions, the kinetics of ion channels, the properties of intracellular proteins and the transmembrane voltage . Ionic channels, exchangers and pumps determine the ion flow through the cell membrane. The sum of these membrane currents defines the temporal change of the transmembrane voltage with the membrane capacity and the intercellular stimulus current as:
(i) Heterogeneous atrial electrophysiology
Experimental measurements of atrial tissue from different anatomical regions of canine (appendage; APG, CT, PM and atrioventricular ring; AVR) show electrophysiological differences (Feng et al. 1998). Myocytes from CT have the most prominent notch and the longest action potential duration (APD) of 270 ms at 1 Hz. Myocytes from PM with 190 ms and those from the APG with 180 ms have intermediate APDs compared to those from the CT. Myocytes from the AVR have the shortest APD with 160 ms and no notch.
These regional differences were considered by adjusting the parameters of the CRN model to the measured data (Seemann et al. 2004a). The maximum conductances of the transient outward potassium current , the L-type calcium current and the rapid delayed rectifier potassium current were adjusted (see table 1) in the different electrophysiological descriptions in the same way as reported by Feng et al. (1998). Figure 2 depicts the varying action potentials (AP) in the measured and simulated case and illustrates the qualitative good agreement between measured and simulated AP heterogeneity.
(ii) Heterogeneous human sinoatrial node electrophysiology
In a previous work, we developed an electrophysiological model of human SAN behaviour (Höper 2003). The CRN model was modified to incorporate several rate constants and the maximum conductances of ion channels of the SAN properties. Some additional channels were included based on the Zhang et al. (2000) rabbit SAN (ZRS) model, which were not reported in human AWM. These changes and inclusions are discussed in short in the following.
is responsible for the stable resting voltage in the AWM. If this current is reduced or blocked, the cells start to depolarize spontaneously. is not present in SAN cells (Irisawa et al. 1993). Thus, was set to 0% in the human SAN model blocking .
The presence of has been shown in cardiac SAN cells (Irisawa et al. 1993). If this current is blocked, a negative chronotropic effect is present in SAN. is not implemented in the CRN model. The equations of the ZRS model were used to describe in the human SAN. Also, was adopted from the ZRS model, since this current was not present in the CRN equations, but in the ZRS model.
The sustained outward K+ current was included in the ZRS model additionally to the transient outward current . The activation variable of and are the same in the ZRS model. The inactivation variable is only used for since has no inactivation. Consequently, we included in the human SAN model with the corresponding characteristics of the CRN model.
The pacemaker current , having a Na+ and a K+ component was adopted from the equations of the ZRS model. The transmembrane voltage , where 50% of the channels are in the open state, is −77 mV. A between −65 and −90 mV was measured in human atrial cardiomyocytes (Ludwig et al. 1999).
The maximum conductance of some ion currents was adjusted for the human SAN model. was reduced to take the prominent reduction in the upstroke velocity of the AP of SAN cells compared to AWM into account (Kodama et al. 1997). The maximum diastolic potential (MDP) of pacemaker cells is determined mainly by the amplitude of . A MDP of mV was measured in human SAN cells (Drouin 1997). was increased to reach this value. The pacemaker potential, which is the phase where the cells activate slowly towards depolarization, the MDP and thus the frequency is influenced by the current , for which the maximum conductance was adjusted for the human SAN model (table 2).
The heterogeneity of the SAN was implemented in the new model in the same way as in the ZRS model. The variation of the maximum conductances in dependence of the position of a cell in the SAN is as follows (Zhang et al. 2000):with the scaling factor and the location of cells (=0 central; =1 peripheral). The maximum conductances of channel type x for central and peripheral cells and , respectively, were set to the same ratios as in the ZRS model. The parameters of central and peripheral cells are listed in table 2.
The CRN model considers intracellular Ca2+ handling including the sarcoplasmic reticulum (SR). The network SR (NSR) and the junctional SR (JSR) have the so-called uptake () and release volume (), respectively. The volumes of NSR and JSR were changed in the human SAN model compared to the CRN model, because the SR in SAN cells is smaller than in AWM (Boyett et al. 2001):with the volume fraction of uptake and release part and , respectively, and the cell volume , and were adopted from data of Boyett et al. (2001) with =0.075, =0.41, =0.06 and =0.3. The varying volume fraction of cells in between was considered by linear interpolation.
A schematic description of the new human SAN model, which was derived from the CRN model and the ZRS model, is depicted in figure 3. The total membrane current is given by:
The simulated transmembrane voltages for central and peripheral cells of human SAN are illustrated in figure 4a. The isolated central SAN cell had a longer cycle length (796 ms) than the peripheral (697 ms). The MDP of a central cell was −51.2 mV and the overshoot 0.3 mV. The MDP of a peripheral cell was −63.3 mV and the overshoot was 9.7 mV. The measured values for peripheral cells are mV and mV (Drouin 1997). Measured values of central cells were not available and were estimated by comparing the results with the ZRS model.
The total Na+, K+ and Ca2+ transmembrane currents for central and peripheral cells are shown in figure 4b,c. The ionic currents of central cells are smaller than of peripheral cells.
(c) Excitation conduction model
The myocardium is a syncytium of excitable cells. To reconstruct this property mathematically a model describing the intercellular current flow must be defined. For this, the bidomain model can be used to represent the anisotropic current flow through gap junctions and the intra- and extracellular space (Henriquez 1993). By assumption of equal anisotropy ratios for the intra- and extracellular space the bidomain model can be translated into the monodomain model, which we used in this work. The intercellular current is calculated with aid of Poisson's equation for stationary electrical fields with the finite element method:with the combined conductivity tensor σ, the membrane surface to cell volume ratio β and an externally applied current source .
The combined conductivity σ consists of conductivities for intra- and extracellular components and gap junctions. σ was adapted to the resolution of 0.33 mm to get a sound velocity. The transversal conductivity was set to 0.15 S m−1 for AWM and CT and to 0.11 S m−1 for PM and BB. The ratio for cross-axis to long-axis conductivity was set to 1 : 9 for anisotropic tissue classes, i.e. CT, PM and BB.
The density of gap junctions decreases from periphery towards centre of the SAN (Honjo et al. 2002). This characteristic was included in the model by setting the conductivity in peripheral cells to 0.11 S m−1 and in central cells to 0.05 S m−1. The conductivity was linearly interpolated for cells inbetween.
(d) Simulation environment
The presented models were combined to an electrically coupled atrial model considering heterogeneous electrophysiological properties in the SAN and the atrium. Different anisotropic features in the bundles were included. The simulation was performed on 10 Apple XServe G5 dual 2 GHz processor cluster nodes using MPI as multiprocessing technique. The model consisted of ca 18.8 million elements. The simulation required 3.8 GB of main memory and needed 37 h of calculation time for a 1 s interval with a temporal increment for the integration of the ionic models and the calculation of the transmembrane voltage distribution of 20 μs. The ordinary differential equations of the ionic model were solved with the forward Euler method. The Poisson's equation is discretized with the finite element method to achieve a system of linear equations.
(a) Anatomical model
The three-dimensional anatomical model (figure 5) includes—besides the already segmented AWM—the SAN, the CT, 10 PMs, the fossa ovalis in the septum, the BB and the right and left atrial APG. They were segmented using interactive deformable meshes and classified manually with aid of anatomical knowledge. This strategy was necessary since the structures are too fine for automatic segmentation and even in the cryosection images not easy to detect.
The SAN lies near the ostium of the superior vena cava and has a heterogeneous distribution of tissue classes. A linear gradient describes this SAN heterogeneity from central to peripheral cells. The periphery of the SAN was everywhere connected directly to adjacent tissue. The CT starts next to the SAN and extends towards inferior vena cava. The PMs are arranged comb like on the right atrial endocardium starting at the side of CT. BB originates next to SAN and extends into the LA towards the left atrial APG.
The anatomical model of the human atrium consisted of 18.8 million cubic voxel including 1.58 million of excitable tissue. The remainder was blood and surrounding tissue. Ten tissue classes were included additionally in this model.
To set fibre orientation of the conduction pathways the principal component analysis was used to find the main axis of each segment of the pathways. Fibre orientation was set for CT, PMs and BB. The resulting orientation was always set parallel to the local main axis of the chosen segment. AWM was set to be isotropic because of conflicting data concerning fibre orientation in human AWM (Boineau et al. 1988). Furthermore, the resolution of the Visible Female data was too low for automatic detection of fibre orientation in AWM.
(b) Sinoatrial node properties
Isolated SAN cells of the periphery showed a higher frequency of autorhythmicity (figure 4). This feature was the same in the model of whole human SAN, but without surrounding atrial tissue (figure 6). This model included heterogeneous tissue distribution from central to peripheral cells. The excitation was initiated primary in peripheral cells (figure 6a). Afterwards central cells got activated (figure 6b). The distribution was relatively homogeneous during the plateau phase (figure 6c–e). The repolarization started in peripheral cells (figure 6f) and propagated afterwards towards central cells (figure 6g), where final repolarization was recognized (figure 6h). The autorhythmic frequency of the isolated SAN model was 79 bpm.
The sequence of activation changed when the cells were surrounded by AWM (figure 7). The of peripheral cells was shifted towards more negative values due to electrotonic effects of the surrounding atrial cells, which had a resting voltage of ca −80 mV. This reduced the autorhythmic frequency of peripheral SAN cells and the final excitation of the SAN started in central cells with 63 bpm (figure 7a). The central activation was driving the peripheral cells to get excited (figure 7b). The repolarization sequence of the SAN model coupled to atrium was comparable to the model without coupling to AWM, especially during plateau phase (figure 7c–h). The repolarization began in peripheral cells due to the longer APD of central cells (figure 7f,g). Furthermore, the surrounding AWM was already repolarized and shifted peripheral cells towards negative voltage. The conduction velocity in central SAN cells was 0.15 m s−1 and in peripheral 0.5 m s−1.
(c) Excitation conduction and repolarization in the atrium
The patterns of the simulation of the whole atrium including SAN are presented in figures 8 and 9. Initiation of excitation started in the centre of the SAN (figure 7a–c). The excitation spread afterwards to the periphery near to the region where CT is coupled to the SAN. Once CT was activated, the excitation propagated preferentially towards the atrioventricular region along the CT (figures 8a–c and 9). The excitation spread on the endocardium of RA along PMs and both CT and PM activated right AWM (figures 8b–d and 9d–f). Earliest activation of LA was found in the region of BB and afterwards excitation spread over to the APG (figure 8d–f). The whole left AWM was activated from BB towards atrioventricular region. During the plateau phase the is varying only little within the atrium (figure 8g). Repolarization started in the area of first activation (figure 8h) and spread similar to the depolarization, but in a more homogeneous manner (figure 8i).
The atrial model was excited periodically due to the pacemaker with a frequency of 63 bpm. The excitation conduction velocities were 0.6 m s−1 for AWM, 1.2 m s−1 for CT, 1.6 m s−1 for PM and 1.1 m s−1 for BB. The time to activation of the RA was 5 ms for BB and 21 ms for the septum. In LA the first activation was found after 24 ms near BB and after 38 ms in septal regions. The endocardial activation times were: 26 ms for BB, 83 ms for RA, 79 ms for LA and 103 ms for the whole atrium. The whole atrium was repolarized after 250 ms.
4. Discussion and conclusion
The simulation of the whole atrium showed the importance of including realistic anatomical structures and fibre orientation combined with realistic heterogeneous human electrophysiology and anisotropic electrical properties. The importance of BB as a dominant interatrial connection during sinus rhythm was verified.
The conduction velocities of the tissue types were in the range of measured data. Hansson et al. (1998) reported 0.68–1.03 m s−1 for AWM. Boineau et al. (1988) measured a CT velocity of 0.7–1.3 m s−1. PMs have been measured with 1.17–1.54 m s−1 (Hayashi et al. 1982). Dolber & Spach (1989) reported a velocity in BB of 0.92–1.67 m s−1. The simulated time to activation and total endocardial activation time were equivalent to the measurement data (Lemery et al.2004) mentioned in §1.
The first activation of BB was earlier in our model compared to the findings of Lemery et al. (2004). This effect can be explained partly by the fact that we were analysing transmembrane voltages and Lemery et al. (2004) determined activation from extracellular potentials. Furthermore, geometrical differences might influence the results significantly. Total endocardial activation time of the whole atrium was shorter in our model with 103 ms. This activation time was longer in the patient study of Lemery et al. (2004) with 120 ms, possibly due to effects of reduced atrial gap junction density of patients with AF (van der Velden & Jongsma 2002). This can additionally explain the shorter times of first activation in our model.
The two breakthrough points in figure 8b have not been reported by Boineau et al. (1988). This might be due to an uncoupled region between PMs and AWM in the model, the neglect of fibre orientation in AWM or that this kind of detail in activation is not visible in the extracellular recordings. The differences between repolarization in atrial endocardium and epicardium reported by Anyukhovsky & Rosenshtraukh (1999) was not integrated in the model because of the lack of channel and cellular data related to this phenomenon. Furthermore, a widely distributed atrial pacemaker complex (Boineau et al. 1988) was not included. This complex is of interest when investigating the atria in unphysiological states.
No measured repolarization data were yet reported. In the model, repolarization was more homogeneous compared to the depolarization because of the heterogeneous APD distribution having longest APD near to the SAN and shortest in distal areas. Repolarization is an important topic to investigate. Its disturbance can lead to arrhythmia due to pathologies reducing the APD or disturbing the dispersion of repolarization.
Our work showed that detailed anatomical structures combined with fast conducting bundles have to be considered in order to reproduce the normal excitation pattern in the atrium seen experimentally. CT, PM and BB are preferential pathways for atrial conduction. Other groups have also worked on modelling whole human atrial activity (Harrild & Henriquez 2000; Blanc et al. 2001; Zemlin et al. 2001), which have not included the anatomical and electrophysiological details presented in this work. Furthermore, they stimulated in the SAN region to initiate atrial excitation. The work of Blanc et al. (2001) and Zemlin et al. (2001) was mainly focused on efficiency for long-time simulations, for which some details were consequently neglected. Blanc et al. (2001) investigated arrhythmias and their treatment by ablation in a surface model, which had no fast-conducting bundles. Zemlin et al. (2001) included in their surface model fibre orientation of AWM but neglected BB. Harrild & Henriquez (2000) have included anatomical structures without heterogeneous electrophysiology. Furthermore, they set slowed conduction in some regions.
There were several limitations of our model. First, the anatomical model of the SAN is inaccurate due to lack of experimental details of structural properties of the human SAN. Second, we neglected the interatrial connection of the coronary sinus, which is not important during sinus rhythm but may contribute during AF. Also, the structure and fibre orientation of the ostia have to be included, e.g. of pulmonary veins. For direct comparison of extracellular measurement data with the model, a bidomain approach has to be implemented to determine the extracellular electrical field. Another limitation is that the heterogeneity is based on canine data for the atrium (Feng et al. 1998) and on rabbit data for the SAN (Kodama et al. 1997, Honjo et al. 2002). We assumed that the human SAN is also electrophysiological heterogeneous since a shift in the site of primary pacemaker in diseased human hearts has been reported similarly to those reported in the rabbit heart (Boyett et al. 2000). Also James (2002, 2003) has reported structural and functional heterogeneity on a microscopic scale in the human SAN. As this model is a macroscopic approach, the detailed interaction between primary pacemaker cells, transitional cells and the adjacent tissue might be considered in future model development as described in the rabbit model of Dobrzynski et al. (2005).
The SAN model was validated partly by single cell and two-dimensional simulations (Höper 2003). For further studies, the cellular models of SAN have to be validated in details against experimental data, to confirm the role of each individual ion channel. Especially the properties of the central SAN cell have to be readjusted if further human measurement data are available. Also the differences between mosaic and gradient approaches to specify the heterogeneity of SAN tissue have to be based on future experimental evidence on human. The heart rate of 63 bpm was found to be in the lower range of expectation, which might be addressed with respect to the neglect of neuronal and other influences on the model.
The detailed model of human atria provides a useful and complementary tool to investigate the dynamical behaviour of atria in such a detailed way that currently cannot be achieved by in vivo and in vitro studies. All electrophysiological parameters during normal and abnormal excitation can be traced and analysed. As such, the model can be used to examine electrophysiological remodelling effects of different cell types induced by chronic AF on genesis and control of AF as described in, for example, the single cell simulation study of Zhang et al. (2003, 2005). The effects of mutations on channel properties, e.g. S140G leading to familial AF and the resulting AP as described in the model of Seemann et al. (2004b) can be included in the whole atrial model to investigate the arrhythmogenic effects of these defects.
In future, a detailed description of anisotropic conductivity of all tissue types can be included and examined. Internodal bundles can be segmented and the complex fibre orientation of AWM can be assigned. The interaction between SAN and adjacent tissue can be examined in more detail. To validate the model, experimentally recorded ECG maps of atrial activation can be compared to simulated ECGs in physiological, pathological, and pacing cases. Future improvements of the model can include electromechanical coupling in the atria or even the whole heart. This model can also incorporate individual anatomy to get more information about the influence of varying geometry on conduction and its role in initiation of arrhythmia.
This work was funded by the German Research Foundation (DFG) under the project H4 of the Collaborative Research Centre SFB414, by the Engineering and Physical Sciences Research Council (EPSRC) with the grant GR/S03027/01 and by the Richard A. and Nora Eccles Fund for Cardiovascular Research and awards from the Nora Eccles Treadwell Foundation. The authors thank the scientific committee of the Oxford meeting in October 2004 to select this work for the award of the Physiological Society.
One contribution of 15 to a Theme Issue ‘Biomathematical modelling II’.
- © 2006 The Royal Society