Abstract

Aims

Evidences of asynchrony between epicardial and endocardial activation in the atrial wall have been reported. We used a computer model of the atria and torso to investigate the consequences of such activation delay on P wave morphology, while controlling for P wave duration.

Methods and results

We created 390 models of the atria based on the same geometry. These models differed by atrial wall thickness (from 2 to 3 mm), transmural coupling, and tissue conductivity in the endocardial and epicardial layers. Among them, 18 were in baseline, 186 had slower conduction in the epicardium layer and 186 in the endocardial layer. Conduction properties were adjusted in such a way that total activation time was the same in all models. P waves on a 16-lead system were simulated during sinus rhythm. Activation maps were similar in all cases. Endo-epicardial delay varied between −5.5 and 5.5 ms vs. 0 ± 0.5 ms in baseline. All P waves had the same duration but variability in their morphology was observed. With slower epicardial conduction, P wave amplitude was reduced by an average of 20% on leads V3–V5 and P wave area decreased by 50% on leads V1–V2 and by 40% on lead V3. Reversed, lower magnitude effects were observed with slower endocardial conduction.

Conclusion

An endo-epicardial delay of a few milliseconds is sufficient to significantly alter P wave morphology, even if the activation map remains the same.

What's new?

  • A computer model of the atria was designed to investigate the changes in P wave morphology resulting from an endo-epicardial delay.

  • A tool was developed to compute the separate endo- and epicardial contributions to the P wave.

  • Slower conduction in the epicardial layer decreased P wave amplitude and area, notably on leads V2–V4, even when P wave duration and activation map remained the same.

  • These changes in P wave morphology were associated with the orientation of the dipole current sources that rotated when the wave front was slanted.

Introduction

P wave descriptors such as P wave duration, amplitude, area and terminal force have been proposed to quantify atrial electrical function.1 These indices have been associated with atrial fibrillation1 and have notably been used for localizing ectopic foci,2–4 predicting atrial fibrillation recurrence after catheter ablation,5,6 and assessing the susceptibility to paroxysmal atrial fibrillation.7 Understanding the factors affecting P wave indices is therefore relevant to the development of predictors of atrial arrhythmias.

Atrial endocardial–epicardial (endo-epi) activation delays have been reported in animal models of atrial fibrillation,8–11 in patients,12,13 and in computer models of the atria.14,15 This delay is hypothesized to be caused by fibrosis, decoupling and dissociation between tissue layers. Histological analysis in a goat model suggested that endomysial fibrosis was stronger in the epicardium than in the endocardium.11 These structural changes are arrhythmogenic16 and alter the current sources that constitute the ionic basis of the electrocardiogram (ECG).17

This arises the question of whether that delay alone can influence P wave morphology. Endo-epi activation delay has been shown to affect the asymmetry of atrial electrograms even if the conduction velocity remains the same.18,19 We hypothesized that P waves would be affected as well. Our aim was to estimate the relation between endo-epi delay and P wave morphology in a computer model in which confounding factors have been as far as possible eliminated.

Methods

Model design

Our objective was to create a mathematical model of the heart and torso in which the atrial endo-epi delay can be controlled and P waves can be computed. When studying the link between endo-epi delay and P wave morphology, the three main confounding factors are (i) the spatial variation in tissue thickness, (ii) the spatial heterogeneity in conduction properties, and (iii) inter-subject variability in P wave duration. To avoid these issues, we created a model with constant thickness throughout the atrial tissue (except in the fast conduction system) and uniform conduction properties in the epicardial and the endocardial layers. In addition, conduction properties were optimized to ensure that the total activation time (and thus the P wave duration) was exactly the same in all simulated conditions.

Atrial geometry

The geometry was derived from an imaging dataset of a patient with atrial fibrillation and included the atria, atrial and ventricular blood cavities, the lungs and the torso, as presented in our previous works.20,21 From the triangulated surface describing the atria, three geometrical models were created with a wall thickness of 2 mm, 2.5 mm, and 3 mm, respectively. These selected values of thickness covered the range of thickness variability observed in most regions of the human atria.22–24

These models were composed of an epicardial layer and an endocardial layer, both with the same thickness. The resulting three-dimensional models were represented by a cubic mesh with a resolution of 0.33 mm (which means 3–5 computational nodes in each of the two layers) and incorporated a simplified fast conducting system (Bachmann’s bundle, crista terminalis, and pectinate muscles, from our previous works20,21) that remained the same in the three geometrical models.

Tissue properties

Propagation of the electrical impulse was simulated in the monodomain framework. To reproduce the conditions of a tissue electrically remodelled by episodes of atrial fibrillation, the Courtemanche et al.25 model was used with ICaL current inhibited by 70%, the Ito current by 50%, and the IKur current by 50%. The numerical methods were the same as previously published.20

The conductivity tensor was determined by three parameters: the longitudinal conductivity in the epicardial layer (σepi) and in the endocardial layer (σendo), and the transmural conductivity (σtm). In the working myocardium, rule-based fibre orientation was introduced as in Krueger et al.26 The anisotropy ratio was 4:1. The conductivity in the fast conducting system, taken from the baseline model in Jacquemet,20 was kept the same in all models.

Conduction properties were adjusted such that the total activation time was 150 ms (see section ‘Activation maps and endo-epi delay’). Because the first and the last few activated cells were difficult to detect on the ECG, this translated to a P wave duration of about 140–150 ms depending on the measurement method.21 This choice corresponds to the category ‘very long P waves’ in Nielsen et al.27 This category is associated with increased risk of atrial fibrillation and is coherent with the dilated atrial geometry used in our study (left atrial volume was about 90 mL) and with our aim to simulate reduced coupling.

For the three thicknesses and six values of σtm between 0.75 and 5 mS/cm, 21 sets of (σepi, σendo) conductivity values (23 sets for 2 mm thickness) were sought that resulted in a total activation time of 150 ms. For that purpose, a regula falsi root finding algorithm was applied to iteratively determine σepi as a function of σendo (or σendo as a function of σepi), generalizing our method designed for two-dimensional tissues.19

Baseline models were characterized by tissues with no endo-epi difference in conduction properties (σepi = σendo). A baseline model was created for each value of σtm by simultaneously optimizing the parameters σepi and σendo under the constraint that they must be equal.

Simulation protocol

In total, there were 390 tissue substrates, including 18 baseline cases: three thicknesses, six transmural conductivities (σtm), and 21 or 23 longitudinal conductivities (σepi). Normal propagation was initiated by injecting intracellular current at the location of the sinoatrial node. The stimulus was 10-ms long with low amplitude to reduce stimulus artefact on the ECG.

Activation maps and endo-epi delay

Local activation time (LAT) was defined as the time the transmembrane potential crosses the −40 mV threshold (computed by linear interpolation with 0.1 ms sampling interval). The total activation time was defined as the time interval between the 0.5% percentile of LAT (near the sinoatrial node) and the 99.5% percentile of LAT (in the posterior left atrium). Exclusion of extreme points improved robustness against the stimulus characteristics and the mechanistic details of the final wave front merging.

The activation map constructed from the LAT was characterized by two components: macroscopic propagation and millimetre-scale delay within wall thickness. The macroscopic component was extracted by Gaussian spatial filtering of the activation map with a characteristic length scale equal to wall thickness.28 This map will be referred to as the thickness-averaged activation map. It enabled the comparison of activation patterns between cases.

The delay map was computed as the difference between the activation map and the thickness-averaged activation map. For each node of the computational grid, the location (z) within the thickness was known from the mesh building process; z varied from −T/2 to T/2, where T is the wall thickness (z > 0 is the epicardium). After exclusion of data points from the septum, the rims of the veins and valves, and the fast conducting system, the local delays were aggregated by z values using bins of width 0.1 mm to reconstruct an average profile representing wavefront shape. The endo-epi delay was defined as the difference between the delays at z = T/2 (epi) and z = −T/2 (endo). The endo-epi delay was positive when the epicardium was activated later than the endocardium. Owing to constant tissue thickness, this averaging approach improved the accuracy of endo-epi delay estimation.

Dipole current sources and body surface potentials

Volume conduction theory was applied to the inhomogeneous torso model to compute the transfer matrix between current dipole sources located within the atria and body surface potentials, as in Jacquemet.20 This approach was based on the boundary element method and assumed isotropic conductivity within each compartment such as blood cavities and lungs. A total of 11 598 dipoles were spread over the atria: 5622 in the epicardium, 5622 in the endocardium, and 354 in the fast conducting system. At each time instant and at each of the 11 598 sites, the dipole source vector was computed as the integral of the local current density over the neighbourhood of the dipole. The local current density was derived from the membrane potential gradient using monodomain theory.20 As a result, the local electrical activity was represented by an equivalent dipole.

The dipole current vector is generally oriented in the direction of wavefront propagation, as illustrated by past studies using oblique dipole layer models.29 The angle between the dipole current vector and the tangent plane to the atrial tissue surface is therefore a measure of wavefront slantedness. This angle would be zero for plane wave propagation in a uniform tissue and 90° for a pure breakthrough wave. The angle of the local current dipole (at the time when dipole magnitude was maximal) was averaged over the 2 × 5622 epi and endo sites and was documented for each simulation.

Body surface potentials can be expressed as the sum of the potentials generated on the torso by each dipole. The advantage of this formulation is that the individual contribution of different sub-parts of the atria to the ECG can be quantified. Here, the contributions of the epicardium, the endocardium, and the fast conducting system were separated. In addition, each dipole current vector was decomposed into a tangent and a normal component (with respect to the atrial surface). This made it possible to determine to what extend the epicardial contribution to the ECG was due to the normal or the tangent component of the dipole sources.

P wave components

P waves, including their epicardial, endocardial, and fast conduction system contributions, as well as their normal and tangent components were computed using a 16-lead system extending the standard 12-lead system by adding leads V3R, V4R, V8, and V9.21 The Wilson central terminal was used as reference.

P wave amplitude was measured in all leads from baseline to the top of the positive deflection. Area under the P wave was summed up algebraically from the onset to the end of the P wave, which implies that it was equal to P wave duration × mean value.

To quantify the relative contribution of each P wave component, the energy of a component was defined as the mean square of the P wave component signals averaged over the 10 precordial leads (squared Frobenius norm). Because the P waves had the same duration in all simulations, this energy can be compared between simulations. The relative importance of the contributions was measured as ratio of energy, for example, energy of the epicardial contribution divided by the energy of the endocardial contribution.

Results

Conductivity parameter optimization

In the 390 tissue substrates, conductivity parameters were optimized such that they all had the same total activation time (an accuracy of 0.1 ms was achieved). The resulting values of conductivity are displayed in Figure 1.

Sets of conductivity parameters that give a total activation time of 150 ms in an atrial model with a wall thickness of (A) 2 mm, (B) 2.5 mm, and (C) 3 mm. Each curve representing endocardial (σendo) vs. epicardial longitudinal conductivity (σepi) corresponds to a value of transmural conductivity (σtm) of 0.75, 2, 3, 4, 5 mS/cm as indicated by the upward arrow. Each data point (circles) is associated with a simulation and the resulting average endo-epi delay is colour-coded. The diagonal dotted line (σendo = σendo) shows conditions with uniform properties along the thickness. The downward diagonal dotted line illustrates the theoretical limit when σtm→ ∞. σendo, endocardial longitudinal conductivity; σepi, epicardial longitudinal conductivity; σtm, transmural conductivity.
Figure 1

Sets of conductivity parameters that give a total activation time of 150 ms in an atrial model with a wall thickness of (A) 2 mm, (B) 2.5 mm, and (C) 3 mm. Each curve representing endocardial (σendo) vs. epicardial longitudinal conductivity (σepi) corresponds to a value of transmural conductivity (σtm) of 0.75, 2, 3, 4, 5 mS/cm as indicated by the upward arrow. Each data point (circles) is associated with a simulation and the resulting average endo-epi delay is colour-coded. The diagonal dotted line (σendo = σendo) shows conditions with uniform properties along the thickness. The downward diagonal dotted line illustrates the theoretical limit when σtm→ ∞. σendo, endocardial longitudinal conductivity; σepi, epicardial longitudinal conductivity; σtm, transmural conductivity.

In baseline models (σepi = σendo), the longitudinal conductivity had to be increased by up to 7% when σtm was reduced from 5 to 0.75 mS/cm. Similarly, an 8% increase in longitudinal conductivity was needed when thickness was reduced from 3 to 2 mm. Because of the curvature of the atrial surface, transmural coupling helped propagation even in the absence of endo-epi differences. The shortest path between two points in the atria may indeed lie partly on the endo- and partly on the epicardial surface.30 Transmural propagation from the pectinate muscles in the right atrium and the high curvature in the septal region also contributed to that effect. Because of these differences, comparison to baseline will always be performed with respect to the baseline model with the same wall thickness and the same transmural conductivity.

There was an inverse relationship between σepi and σendo (Figure 1). A reduced coupling in one layer can be compensated by stronger coupling in the other layer. In the limit, where the tissue is very thin or when layers are strongly coupled, only the mean of σepi and σendo matters,31 so this relationship would become linear (downward dotted line in Figure 1). When transmural coupling is weak, propagation is driven by the faster layer, so this relationship would become almost independent from the slower conductivity19; geometrically, the curve converges to a square when σtm tends to zero.

The range of σepi and σendo obtained was smaller in the thicker models. In thin models, larger endo-epi differences in conductivity were needed to generate similar conduction delays. This was confirmed by the calculated endo-epi activation delays (colour coded in Figure 1). While the maximum delay was 3 ms in the 2-mm thick model, delays over 5 ms were obtained in the 3-mm thick model.

Thickness-averaged activation maps

Although total activation time was the same in all cases, activation patterns may differ, with potential impact on P wave morphology. Thickness-averaged activation maps were compared with quantify these differences. Figure 2A and B show the baseline case and a case with an endo-epi delay of 3 ms (epi was activated later than endo). The activation maps seem identical visually; Figure 2C shows the difference between the two. Because the septum route to the left atrium relied more on the epicardial layer (due to the curvature), activation occurred later in the right pulmonary veins when compared with baseline. In contrast, the Bachmann’s bundle route was faster than baseline in such a way that the total activation ended at the same time. As demonstrated on the histogram of Figure 2D, the difference between thickness-averaged activation maps was of the order of a few milliseconds.

Activation maps averaged across the thickness. (A) Baseline model with a wall thickness of 2.5 mm and a transmural conductivity (σtm) of 4 mS/cm. LAT is colour-coded. Posterior view is shown. (B) Thickness-averaged activation map with σepi = 1 mS/cm, σendo = 8.28 mS/cm, and σtm = 4 mS/cm. (C) Difference between activation maps (A) and (B), displayed using the colour code of the histogram (D) of the differences in LAT. Blue means that activation occurs later in the non-uniform model when compared with baseline. The asterisk identifies the latest activated point where the final wave front fusion occurs. LAT, local activation time; σepi, epicardial longitudinal conductivity; σendo, endocardial longitudinal conductivity; σtm, transmural conductivity.
Figure 2

Activation maps averaged across the thickness. (A) Baseline model with a wall thickness of 2.5 mm and a transmural conductivity (σtm) of 4 mS/cm. LAT is colour-coded. Posterior view is shown. (B) Thickness-averaged activation map with σepi = 1 mS/cm, σendo = 8.28 mS/cm, and σtm = 4 mS/cm. (C) Difference between activation maps (A) and (B), displayed using the colour code of the histogram (D) of the differences in LAT. Blue means that activation occurs later in the non-uniform model when compared with baseline. The asterisk identifies the latest activated point where the final wave front fusion occurs. LAT, local activation time; σepi, epicardial longitudinal conductivity; σendo, endocardial longitudinal conductivity; σtm, transmural conductivity.

The median difference in LAT between each simulated activation map and its corresponding baseline case ranged from −0.2 to 1 ms and its interquartile range was <1.5 ms for endo-epi delays <1 ms and up to 2.8 ms for the largest endo-epi delays. Overall, the differences were always at most of the order of a few milliseconds. Considering that at the clinical sampling rate of 1 kHz, the uncertainty on LAT estimation would be at least 1 ms, we will consider these thickness-averaged activation maps to be essentially identical for practical purpose.

Activation delay and wave front shape

The activation map minus the thickness-averaged activation map gave the local delay. Figure 3A shows an example corresponding to the case from Figure 2B. Owing to the uniform tissue properties, the local delay was rather uniform, both in the epicardium and the endocardium. Notable exceptions were at the sites of wave front collisions and in the fast conduction system that was excluded from the analysis.

Average wave front shape and determination of the endo-epi delay. (A) Difference between the activation map and the thickness-averaged activation map, displaying the local endo-epi delay (blue means that the epicardium is activated later). (B) Wave front shape averaged over the entire tissue. The local delay (x axis) is shown as a function of the position within the wall (y axis). The dark grey area shows the standard deviation of the local delay. The scaling between x and y axes is chosen such that at the average propagation velocity (66 cm/s), the curve would correspond to the shape of an isopotential line. Among the nine examples shown, the wall thickness is the same in each row and the transmural conductivity is the same in each column (σtm = 0.75, 2, and 4 mS/cm from left to right). σtm, transmural conductivity.
Figure 3

Average wave front shape and determination of the endo-epi delay. (A) Difference between the activation map and the thickness-averaged activation map, displaying the local endo-epi delay (blue means that the epicardium is activated later). (B) Wave front shape averaged over the entire tissue. The local delay (x axis) is shown as a function of the position within the wall (y axis). The dark grey area shows the standard deviation of the local delay. The scaling between x and y axes is chosen such that at the average propagation velocity (66 cm/s), the curve would correspond to the shape of an isopotential line. Among the nine examples shown, the wall thickness is the same in each row and the transmural conductivity is the same in each column (σtm = 0.75, 2, and 4 mS/cm from left to right). σtm, transmural conductivity.

Even in baseline, endo-epi delays were observed due to atrial surface curvature and three-dimensional structures. The effect was stronger in a thicker tissue and when transmural coupling was low. The mean delay was almost zero and its standard deviation was about 0.5, 0.8, and 1 ms for a wall thickness of 2, 2.5, and 3 mm, respectively.

Since the local delay in baseline was location-dependent, which component was subtracted from the local delay map of the non-baseline cases to reduce spatial variability. This enabled averaging of the transmural profile of the delay over space, thus increasing resolution and accuracy. Examples of the resulting activation profiles are shown in Figure 3B. The wave fronts appear to be more slanted in the slower layer than in the faster layer that drives propagation.

Decomposition of the P wave in baseline

Figure 4 shows the P wave on lead V3 for all baseline models (with different transmural conductivities and a thickness of 2.5 mm). The P wave was decomposed into three components: the contribution of the epicardium, the endocardium, and the fast conducting system. The epi and endo contributions were similar in shape and amplitude. Their differences, mostly seen during activation of the right atrium, were associated with the wavefront coming from three-dimensional structures such as the pectinate muscles. The nature of propagation was highlighted by further decomposition into normal and tangent components. The normal components, associated with transmural propagation, were very small, reflecting that propagation in the epi and endo layers was essentially synchronous in baseline.

The P wave on lead V3 at baseline and its different components. The curves for all six values of transmural conductivity (wall thickness: 2.5 mm) are superimposed. The P wave is the sum of the contribution of the epicardium (epi), the endocardium (endo), and the fast conduction system (bundles). Each epi and endo component is further decomposed into the contribution of the dipole current sources that are normal and tangent to the tissue surface.
Figure 4

The P wave on lead V3 at baseline and its different components. The curves for all six values of transmural conductivity (wall thickness: 2.5 mm) are superimposed. The P wave is the sum of the contribution of the epicardium (epi), the endocardium (endo), and the fast conduction system (bundles). Each epi and endo component is further decomposed into the contribution of the dipole current sources that are normal and tangent to the tissue surface.

P wave morphology in the presence of endo-epi delay

A uniform endo-epi delay affected P wave morphology but not its duration, as illustrated by the changes on V3 in Figure 5. In simulations with slower epi, the amplitude of the epi contribution to the P wave was reduced (Figure 5A). This was largely due to the decrease of the tangent component. The wavefront shape was more slanted in the epicardium, leading to a predominantly transmural propagation (Figure 3B). This translated into a much larger normal component. In contrast, the amplitude of the endo contribution (both the tangent and normal components) increased uniformly. In simulations with slower endo, opposite effects were observed (Figure 5B).

Variations in P wave morphology on lead V3 in the presence of endo-epi delays. The same display format as Figure 4 is used. (A) Slower conduction in the epicardium; (B) slower conduction in the endocardium. The P waves from all 126 simulations with a wall thickness of 2.5 mm are included in the plot. The dark lines display the baseline P waves from Figure 4. The interquartile interval and the min–max range of all P waves (computed separately for each time instant) are shown.
Figure 5

Variations in P wave morphology on lead V3 in the presence of endo-epi delays. The same display format as Figure 4 is used. (A) Slower conduction in the epicardium; (B) slower conduction in the endocardium. The P waves from all 126 simulations with a wall thickness of 2.5 mm are included in the plot. The dark lines display the baseline P waves from Figure 4. The interquartile interval and the min–max range of all P waves (computed separately for each time instant) are shown.

The changes in P wave morphology can be quantified using P wave amplitude and area. Figure 6 presents the statistics of these parameters for all leads in baseline (n = 18) and in models with slower epi (n = 186) or slower endo (n = 186). Models with slower epi had P wave amplitude reduced by approximately 20% vs. baseline on leads V3–V5 and by 12% on leads V1–V2. The most dramatic effect was a decrease in P wave area of over 50% on leads V1–V2 and 40% on lead V3. Overall, P wave area was smaller than baseline (in absolute value) in all leads except V9. In models with slower endo, P wave amplitude was significantly smaller on leads V4R, V3R, and V8–V9. In contrast with the slower epi cases, P wave area was increased on lead V2 and V3 by 29% and 18% respectively, while it was reduced by 37% in the back on lead V9.

Amplitude (A) and area (B) of the P waves in the 10 precordial leads. Means and standard deviations are computed over all baseline cases (n = 18), all models with slower epi (n = 186), and all models with slower endo (n = 186). Asterisk indicate statistical significance (t-test vs. baseline, *P < 0.05, **P < 0.001). AUC, area under the curve; ECG, electrocardiogram.
Figure 6

Amplitude (A) and area (B) of the P waves in the 10 precordial leads. Means and standard deviations are computed over all baseline cases (n = 18), all models with slower epi (n = 186), and all models with slower endo (n = 186). Asterisk indicate statistical significance (t-test vs. baseline, *P < 0.05, **P < 0.001). AUC, area under the curve; ECG, electrocardiogram.

Relative contributions of the P wave components

The weights of the epi and endo contributions were quantified using their energy computed from all leads. Figure 7A shows epi/endo energy ratio as a function of epi/endo conductivity ratio. These two ratios were of the same order of magnitude. Lower conductivity in a layer resulted in a smaller contribution. Tissue conductivity was indeed a multiplicative factor in the formula for computing current sources.20 Near baseline (σepi/σendo ratio of 1), the relation was a power law (energy ratio ∝ (σepi/σendo)1.7).

Relative contribution of the different P wave components, measured using the energy based on all precordial leads. (A) Ratio of the energy of the epi and endo components as a function of the ratio of their conductivity. Wall thickness is colour-coded. (B) Energy of the normal component (epi and endo combined) divided by the total energy of the P wave, displayed as a function of the angle between the equivalent dipole current source and the tissue surface.
Figure 7

Relative contribution of the different P wave components, measured using the energy based on all precordial leads. (A) Ratio of the energy of the epi and endo components as a function of the ratio of their conductivity. Wall thickness is colour-coded. (B) Energy of the normal component (epi and endo combined) divided by the total energy of the P wave, displayed as a function of the angle between the equivalent dipole current source and the tissue surface.

Slanted wavefronts resulted in an angle between the local equivalent dipole and the atrial surface, which in turns generated a normal component on the P wave (Figure 5). This principle was tested in Figure 7B, where the energy of the normal component is displayed as a function of the equivalent dipole angle. A parabolic relationship was found. The dipole angle varied in the range −30° to 40° (positive means outward-oriented dipole). Notice that the parabola is slightly shifted to the left since there is a small angle in baseline. Therefore, the normal component is associated with the angle of the equivalent dipole.

This arises the question about the factors that affect the angle of the equivalent dipole. From a previous two-dimensional study,19 plausible candidates were the endo-epi delay, tissue thickness, and also transmural coupling. Multivariate analysis (Figure 8) showed that the dipole angle was correlated to the product of σtm and endo-epi delay divided by wall thickness (correlation coefficient > 0.98; mean square error: 2.2°).

Factors affecting the angle of the equivalent dipole current source. The angle of the equivalent dipole is displayed as a function of the product between the transmural conductivity (mS/cm) and the endo-epi delay (ms), divided by the wall thickness (mm). Positive values correspond to slower conduction in the epicardium.
Figure 8

Factors affecting the angle of the equivalent dipole current source. The angle of the equivalent dipole is displayed as a function of the product between the transmural conductivity (mS/cm) and the endo-epi delay (ms), divided by the wall thickness (mm). Positive values correspond to slower conduction in the epicardium.

Discussion

Main findings and mechanisms

The central result of this article is that a few milliseconds of delay between epicardial and endocardial atrial activation is sufficient to generate significant changes in P wave morphology, notably in terms of amplitude and area under the curve (Figure 6). This hypothesis was tested in a model specifically designed for investigating the effect of endo-epi delay, wall thickness and conductivity. The endo-epi delay was rather uniform by design (Figure 3A). This might be seen as a worst-case scenario, for instance, when compared with a situation where a delay is found only in the left atrium. On the other hand, there may be conditions (possibly during atrial fibrillation) in which the combination of positive and negative delays in different parts of the atria creates an even larger effect on some leads.

Thanks to the ability of computer models to decompose the P wave into epi and endo components (Figure 5), the mechanism underlying the alteration of P wave morphology caused by endo-epi delay was found to be the change in the orientation of dipole current sources in the myocardium. Because the atrial wall was thin, a delay of a few milliseconds was sufficient to significantly slant the wavefront (Figure 3B), thereby rotating the dipole source vectors by up to 40° (Figure 7B). The angle of the dipole vector appeared to be proportional to the endo-epi delay divided by tissue thickness (Figure 8), which, at constant conduction velocity, corresponds to the angle of wavefront shape. Another factor was transmural coupling (σtm), a parameter involved explicitly in the formula for computing the normal component of the dipole source.

Thickness-averaged activation maps were almost the same in all cases (Figure 2A) up to a few milliseconds. In clinical signals, a local difference in activation time of a few samples would not be considered relevant because of the limited accuracy of activation time estimation in the presence of fibrosis or electrogram fractionation. Still, the differences vs. baseline were of the same order of magnitude as the endo-epi delay. This is where the decomposition into tangent and normal components was insightful. Alteration of the thickness-averaged activation would only generate a tangent contribution, so the increase of the normal component can be definitely attributed to the endo-epi delay.

Comparison with experimental data

Schuessler et al.9 measured an endo-epi delay of 0.8 ± 5.7 ms in the isolated canine right atrium during sinus rhythm. The standard deviation increased to 6–8 ms during pacing and up to 10–28 ms during atrial fibrillation. The same range of endo-epi delay was observed in the human right atrium during atrial fibrillation, with typically up to 26% of the activations having >15 ms endo-epi delay.12 Similar values were obtained in a bilayer computer model of endo-epi dissociation.32 However, during atrial fibrillation, decoupling of layers had also significant impact on activation patterns (e.g. breakthroughs), in contrast with our simulations.

Simulated endo-epi delays were in the range of 0–5.5 ms during sinus rhythm for a thickness up to 3 mm, which is within the range of Schuessler et al.9 Endo-epi delays were derived from the activation time computed based on the membrane potential time course. Estimation from epicardial and endocardial electrograms as performed in experiments may be less accurate.

Implications for atrial modelling

The results are relevant to non-invasive determination of atrial activation map from body surface potentials.33,34 They suggest that a priori information or hypotheses about possible endo-epi delay (e.g. its absence in most of the atrial tissue) is needed to solve the inverse problem in the atria. Clinical late gadolinium-enhanced magnetic resonance imaging currently does not provide yet sufficient resolution to estimate fibrosis density separately in the epicardial and the endocardial layers.35 This information may have to be inferred from or combined with other measurements.

Along the same lines, patient-tailoring of conduction properties in computer models of the atria are typically based on reproducing clinically-measured endocardial activation maps and body surface potentials.36–38 Unknown local endo-epi delays would increase the number of degrees of freedom and complicates the optimization procedure.

Limitations

Slower conduction in a layer was simulated by uniform reduced conductivity. This reproduced macroscopic activation patterns but ignored the microstructural changes involved in fibrosis development. The arbitrary distinction of two layers of equal thickness was also a simplifying hypothesis that came from the limited resolution of the computational mesh and from our objective to restrain the number of adjustable parameters. Nevertheless, due to the smoothing effect of the volume conductor, the local current dipole sources aggregated at the scale of a few millimetres are expected to be the main determinants of P wave morphology.

Model limitations also include simplified three-dimensional structures and constant thickness, both introduced by design to improve reliability of statistical analysis. Local variations in wall thickness, notably in the left atrium, may play a role in the genesis of the P wave as well as in the wave dynamics of atrial fibrillation.23 Future works would benefit from more realistic spatial variations in left atrial wall thickness.

The volume conduction model was composed of uniform isotropic regions. Extracellular myocardial anisotropy was neglected in ECG computations. Our previous works suggested that this effect was small in thin-walled (<3 mm) models.39 This conclusion may have to be revised in models with spatially-varying thickness. Monodomain conductivity was used in the simulations and the figures without specifying the relative contributions of the intra- and extracellular conductivities. The situation was complicated by the necessarily unequal anisotropy ratios created by modifying transmural coupling. It may be hypothesized that extracellular conductivities (σe) remain the same in all models and that slower propagation is attributed to reduced intracellular coupling (σi). In the transmural direction (the component that drives the rotation of the current dipole), σi may be obtained as 1/(σtm−1σe−1) > σtm. Because current density is proportional to σi, this is equivalent to scaling up the normal components in Figure 5, which provides another application of P wave decomposition. More accurate volume conduction models incorporate both intra- and extracellular conductivity tensors.40,41

Transmural coupling was not decreased to very low values in order to avoid numerical artefacts caused by cubic mesh discretization. Because of coarse spatial discretization, the conductivity curves in Figure 1 should be seen as qualitative relations. In the slow layer, conductivity is overestimated. It is therefore easier (in terms of epi/endo conductivity ratio) to create an endo-epi delay. A finer mesh would be needed to reach convergence. However, the fundamental mechanisms leading to changes in P wave morphology only involved the orientation of dipole current sources and would remain valid in more complex models.

Because P wave duration was controlled, this study does not directly apply to the evolution of P wave morphology when fibrosis progressively structurally remodels the epicardium. In such case, P wave duration is expected to be prolonged with potential impact on P wave morphology caused by a new activation pattern. Accordingly, clinical validation of endo-epi delay-induced changes in the P wave would be difficult because of the inability to specifically modify the delay in the same subject while keeping all other parameters intact. This was the argument that prompted this computational study.

Conclusion

The clinical relevance of this study to ECG interpretation is that in addition to the changes in P wave morphology caused by activation patterns slowed down by reduced coupling, endo-epi delays alone can further affect the P waves, particularly on leads V2 to V4. This may constitute a confounding factor in the analysis of P wave morphology.

Funding

This work was supported by the Natural Sciences and Engineering Research Council of Canada [NSERC RGPIN-2015–05658].

Conflict of interest: none declared.

References

1

Magnani
JW
,
Zhu
L
,
Lopez
F
,
Pencina
MJ
,
Agarwal
SK
,
Soliman
EZ
et al.
P-wave indices and atrial fibrillation: cross-cohort assessments from the Framingham Heart Study (FHS) and Atherosclerosis Risk in Communities (ARIC) study
.
Am Heart J
2015
;
169
:
53
61.e51
.

2

Fitzgerald
DM
,
Hawthorne
HR
,
Crossley
GH
,
Simmons
TW
,
Haisty
WK
Jr.
P wave morphology during atrial pacing along the atrioventricular ring. ECG localization of the site of origin of retrograde atrial activation
.
J Electrocardiol
1996
;
29
:
1
10
.

3

Kistler
PM
,
Roberts-Thomson
KC
,
Haqqani
HM
,
Fynn
SP
,
Singarayar
S
,
Vohra
JK
et al.
P-wave morphology in focal atrial tachycardia: development of an algorithm to predict the anatomic site of origin
.
J Am Coll Cardiol
2006
;
48
:
1010
7
.

4

Tang
CW
,
Scheinman
MM
,
Van Hare
GF
,
Epstein
LM
,
Fitzpatrick
AP
,
Lee
RJ
et al.
Use of P wave configuration during atrial tachycardia to predict site of origin
.
J Am Coll Cardiol
1995
;
26
:
1315
24
.

5

Lim
TW
,
Wu
G
,
Ross
DL
,
Thomas
SP.
P-wave measurements and electrical reconnection of the posterior left atrium after catheter ablation for atrial fibrillation
.
Pacing Clin Electrophysiol
2010
;
33
:
1324
34
.

6

Okumura
Y
,
Watanabe
I
,
Ohkubo
K
,
Ashino
S
,
Kofune
M
,
Hashimoto
K
et al.
Prediction of the efficacy of pulmonary vein isolation for the treatment of atrial fibrillation by the signal-averaged P-wave duration
.
Pacing Clin Electrophysiol
2007
;
30
:
304
13
.

7

Conte
G
,
Luca
A
,
Yazdani
S
,
Caputo
ML
,
Regoli
F
,
Moccetti
T
et al.
Usefulness of P-wave duration and morphologic variability to identify patients prone to paroxysmal atrial fibrillation
.
The Am J Cardiol
2017
;
119
:
275
9
.

8

Eckstein
J
,
Verheule
S
,
de Groot
NM
,
de Groot
N
,
Allessie
M
,
Schotten
U.
Mechanisms of perpetuation of atrial fibrillation in chronically dilated atria
.
Prog Biophys Mol Biol
2008
;
97
:
435
51
.

9

Schuessler
RB
,
Kawamoto
T
,
Hand
DE
,
Mitsuno
M
,
Bromberg
BI
,
Cox
JL
et al.
Simultaneous epicardial and endocardial activation sequence mapping in the isolated canine right atrium
.
Circulation
1993
;
88
:
250
63
.

10

Verheule
S
,
Eckstein
J
,
Linz
D
,
Maesen
B
,
Bidar
E
,
Gharaviri
A
et al.
Role of endo-epicardial dissociation of electrical activity and transmural conduction in the development of persistent atrial fibrillation
.
Prog Biophys Mol Biol
2014
;
115
:
173
85
.

11

Verheule
S
,
Tuyls
E
,
Gharaviri
A
,
Hulsmans
S
,
van Hunnik
A
,
Kuiper
M
et al.
Loss of continuity in the thin epicardial layer because of endomysial fibrosis increases the complexity of atrial fibrillatory conduction
.
Circ Arrhythm Electrophysiol
2013
;
6
:
202
11
.

12

de Groot
N
,
van der Does
L
,
Yaksh
A
,
Lanters
E
,
Teuwen
C
,
Knops
P
et al.
Direct proof of endo-epicardial asynchrony of the atrial wall during atrial fibrillation in humans
.
Circ Arrhythm Electrophysiol
2016
;
9
:
e003648
.

13

Lemery
R
,
Birnie
D
,
Tang
AS
,
Green
M
,
Gollob
M
,
Hendry
M
et al.
Normal atrial activation and voltage during sinus rhythm in the human heart: an endocardial and epicardial mapping study in patients with a history of atrial fibrillation
.
J Cardiovasc Electrophysiol
2007
;
18
:
402
8
.

14

Gharaviri
A
,
Verheule
S
,
Eckstein
J
,
Potse
M
,
Kuijpers
NH
,
Schotten
U.
A computer model of endo-epicardial electrical dissociation and transmural conduction during atrial fibrillation
.
Europace
2012
;
14
:
v10
6
.

15

Labarthe
S
,
Bayer
J
,
Coudiere
Y
,
Henry
J
,
Cochet
H
,
Jais
P
et al.
A bilayer model of human atria: mathematical background, construction, and assessment
.
Europace
2014
;
16
:
iv21
9
.

16

Schotten
U
,
Verheule
S
,
Kirchhof
P
,
Goette
A.
Pathophysiological mechanisms of atrial fibrillation: a translational appraisal
.
Physiol Rev
2011
;
91
:
265
325
.

17

Jacquemet
V
,
van Oosterom
A
,
Vesin
JM
,
Kappenberger
L.
Analysis of electrocardiograms during atrial fibrillation. A biophysical model approach
.
IEEE Eng Med Biol Mag
2006
;
25
:
79
88
.

18

Houben
RP
,
de Groot
NM
,
Smeets
JL
,
Becker
AE
,
Lindemans
FW
,
Allessie
MA.
S-wave predominance of epicardial electrograms during atrial fibrillation in humans: indirect evidence for a role of the thin subepicardial layer
.
Heart Rhythm
2004
;
1
:
639
47
.

19

Irakoze
E
,
Gowda
CHR
,
Jacquemet
V.
Asymmetry of unipolar electrograms in a thin tissue with epicardial-endocardial activation delay
.
Comput Cardiol
2017
;
44
:
019
5
.

20

Jacquemet
V.
Modeling left and right atrial contributions to the ECG: a dipole-current source approach
.
Comput Biol Med
2015
;
65
:
192
9
.

21

Saha
M
,
Conte
G
,
Caputo
ML
,
Regoli
F
,
Krause
R
,
Auricchio
A
et al.
Changes in P-wave morphology after pulmonary vein isolation: insights from computer simulations
.
Europace
2016
;
18
:
iv23
34
.

22

Platonov
PG
,
Ivanov
V
,
Ho
SY
,
Mitrofanova
L.
Left atrial posterior wall thickness in patients with and without atrial fibrillation: data from 298 consecutive autopsies
.
J Cardiovasc Electrophysiol
2008
;
19
:
689
92
.

23

Song
JS
,
Wi
J
,
Lee
HJ
,
Hwang
M
,
Lim
B
,
Kim
TH
et al.
Role of atrial wall thickness in wave-dynamics of atrial fibrillation
.
PloS One
2017
;
12
:
e0182174.

24

Whitaker
J
,
Rajani
R
,
Chubb
H
,
Gabrawi
M
,
Varela
M
,
Wright
M
et al.
The role of myocardial wall thickness in atrial arrhythmogenesis
.
Europace
2016
;
18
:
1758
72
.

25

Courtemanche
M
,
Ramirez
RJ
,
Nattel
S.
Ionic targets for drug therapy and atrial fibrillation-induced electrical remodeling: insights from a mathematical model
.
Cardiovasc Res
1999
;
42
:
477
89
.

26

Krueger
M
,
Schmidt
V
,
Tobón
C
,
Weber
F
,
Lorenz
C
,
Keller
D
et al. Modeling atrial fiber orientation in patient-specific geometries: a semi-automatic rule-based approach. In
Metaxas
DN
,
Axel
L
(eds)
Functional Imaging and Modeling of the Heart
.
Berlin, Heidelberg
:
Springer
;
2011
. pp.
223
32
.

27

Nielsen
JB
,
Kuhl
JT
,
Pietersen
A
,
Graff
C
,
Lind
B
,
Struijk
JJ
et al.
P-wave duration and the risk of atrial fibrillation: results from the Copenhagen ECG Study
.
Heart Rhythm
2015
;
12
:
1887
95
.

28

Aswath Kumar
AK
,
Drahi
A
,
Jacquemet
V.
Fitting local repolarization parameters in cardiac reaction-diffusion models in the presence of electrotonic coupling
.
Comput Biol Med
2017
;
81
:
55
63
.

29

Colli-Franzone
P
,
Guerri
L
,
Viganotti
C
,
Macchi
E
,
Baruffi
S
,
Spaggiari
S
et al.
Potential fields generated by oblique dipole layers modeling excitation wavefronts in the anisotropic myocardium. Comparison with potential fields elicited by paced dog hearts in a volume conductor
.
Circ Res
1982
;
51
:
330
46
.

30

van Dam
PM
,
Oostendorp
TF
,
van Oosterom
A.
Application of the fastest route algorithm in the interactive simulation of the effect of local ischemia on the ECG
.
Med Biol Eng Comput
2009
;
47
:
11
20
.

31

Coudiere
Y
,
Henry
J
,
Labarthe
S.
A two layers monodomain model of cardiac electrophysiology of the atria
.
J Math Biol
2015
;
71
:
1607
41
.

32

Gharaviri
A
,
Verheule
S
,
Eckstein
J
,
Potse
M
,
Kuklik
P
,
Kuijpers
NH
et al.
How disruption of endo-epicardial electrical connections enhances endo-epicardial conduction during atrial fibrillation
.
Europace
2017
;
19
:
308
18
.

33

Figuera
C
,
Suarez-Gutierrez
V
,
Hernandez-Romero
I
,
Rodrigo
M
,
Liberos
A
,
Atienza
F
et al.
Regularization techniques for ECG imaging during atrial fibrillation: a computational study
.
Front Physiol
2016
;
7
:
466
.

34

Zhou
Z
,
Jin
Q
,
Yu
L
,
Wu
L
,
He
B.
Noninvasive imaging of human atrial activation during atrial flutter and normal rhythm from body surface potential maps
.
PloS One
2016
;
11
:
e0163445
.

35

Marrouche
NF
,
Wilber
D
,
Hindricks
G
,
Jais
P
,
Akoum
N
,
Marchlinski
F
et al.
Association of atrial tissue fibrosis identified by delayed enhancement MRI and atrial fibrillation catheter ablation: the DECAAF study
.
JAMA
2014
;
311
:
498
506
.

36

Boyle
PM
,
Zahid
S
,
Trayanova
NA.
Towards personalized computational modelling of the fibrotic substrate for atrial arrhythmia
.
Europace
2016
;
18
:
iv136
45
.

37

Krueger
MW
,
Schulze
WH
,
Rhode
KS
,
Razavi
R
,
Seemann
G
,
Dossel
O.
Towards personalized clinical in-silico modeling of atrial anatomy and electrophysiology
.
Med Biol Eng Comput
2013
;
51
:
1251
60
.

38

Krueger
MW
,
Seemann
G
,
Rhode
K
,
Keller
DU
,
Schilling
C
,
Arujuna
A
et al.
Personalization of atrial anatomy and electrophysiology as a basis for clinical modeling of radio-frequency ablation of atrial fibrillation
.
IEEE Trans Med Imaging
2013
;
32
:
73
84
.

39

Jacquemet
V.
Equivalent dipole sources to estimate the influence of extracellular myocardial anisotropy in thin-walled cardiac forward models
.
Math Biosci
2017
;
286
:
31
8
.

40

Potse
M
,
Dube
B
,
Vinet
A.
Cardiac anisotropy in boundary-element models for the electrocardiogram
.
Med Biol Eng Comput
2009
;
47
:
719
29
.

41

Potse
M
,
Lankveld
TA
,
Zeemering
S
,
Dagnelie
PC
,
Stehouwer
CD
,
Henry
RM
et al.
P-wave complexity in normal subjects and computer models
.
J Electrocardiol
2016
;
49
:
545
53
.

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://dbpia.nl.go.kr/journals/pages/open_access/funder_policies/chorus/standard_publication_model)