IJBEM, Vol. 2, No. 2, pp. 321-340, 2000.

Myocardial Anisotropy and Multiphasic Electrograms -
Simulation Study in a Monoventricular Mode

* Department of Mathematics, University of Piemonte Orientale, Alessandria, Italy
** Department of Mathematics, University of Pavia, Pavia, Italy
*** Institute of Numerical Analysis - CNR, Pavia, Italy
**** CVRTI - University of Utah, Salt Lake City, UT, USA

Abstract. Unipolar epicardial electrograms are affected by the fiber structure and anisotropy of the myocardial tissue. The electrogram shape can vary from monophasic to pentaphasic. This problem is investigated by means of a mathematical model of the left ventricle with the myocardium represented by an anisotropic bidomain. Simulated and recorded electrograms are compared. The analysis of the electrogram shape is facilitated by a splitting of the sources, the potential distribution and the electrograms into axial and conormal components.

Keywords: Electrograms, bidomain model, excitation wavefront, cardiac sources

Introduction

The excitation process in the myocardium, the associate potential distribution and the electrograms are influenced by the fiber structure and the anisotropy of the myocardial tissue. In fact propagation is faster along than across fiber, and positive potentials are more frequently observed in those areas toward which excitation spreads along fibers. As a consequence, the shape of the unipolar electrograms (EGs) depends on whether the recording site is reached by an excitation wave moving along or across or obliquely to the fiber direction [13]. Only few simulations of EGs based on the full anisotropic bidomain model have appeared in literature [5,14,1,10,21] A preliminary study of the effect of the reference potential on unipolar electrograms (EGs) in a large slab of myocardium was carried out in [21]. In this paper we use a mathematical model of the left ventricle, based on the anisotropic bidomain, taking into account the intramural rotation and obliqueness of the fibers and a simplified Purkinje network. The computation of EGs requires a refined numerical technique to avoid spurious oscillations and peaks. This adaptive procedure, which has some affinity with that proposed in [5], was first used in [2] and was adapted to the model of this study. The analysis of the EG shape is facilitated by the use of the splitting technique applied to sources, potential fields and EGs decomposing them into axial and conormal components. Also the effect of the drift of the reference potential must be considered.

Methods

Mathematical Model

The numerical simulations were carried out on a monoventricular model previously used to simulate the spread of excitation and the associated potential distributions. The model, representing the left ventricle as an ellipsoid symmetric with respect to the z-axis and truncated at the basis and the apex, takes into account the fiber rotation counterclockwise (CCW) from epicardium (-45^o) to endocardium (75^o). and the epi- endocardial obliqueness see [18]. Thus, the fibers do not lie on the nested ellipsoidal surfaces but intersect them at a small angle. For more details on the geometry and on the fiber architecture of our model of the left ventricle see [3].

The ventricular model, used in this simulation study, incorporates the main electrical and structural features of the myocardium such as :

a) the anisotropic electrical conductivities of the intra- and extra-cellular media with unequal anisotropic ratio;

b) the curvature and variable thickness of the ventricular wall;

c) the intramural fiber rotation and obliqueness, whose combined effect results in fiber pathways developing on nested toroidal surfaces inside the ventricular wall (see [3,18]);

d) a simplified subendocardial Purkinje network;

e) two fluid layers, one having conductivity similar to the blood and lining the endocardium, the other in contact with the epicardial surface, with conductivity similar to the average conductivity of chest tissues (i.e., 2 $\times 10^{-3}\Omega^{-1} \hbox{cm}^{-1}$ ).

In our model the ventricular cavity is not completely filled by blood and the fluid layer in contact with the epicardium has a thickness varying from 2.5 to 5 mm, not comparable with that of the conducting fluid filling the experimental tank. Nevertheless, the model is a good approximation of the experimental setup related to isolated dog hearts, especially for epicardial and intramural observation points.

Therefore, most of the essential factors affecting the EGs are present in both the experimental and the modeling framework. This allowed us to perform a qualitative comparison between measured and computed EGs.

We briefly recall that the mathematical model used to compute EGs is based on the bidomain approach where two superposed and connected media, the intracellular (i) and the extracellular (e), are considered [7,9].

In the following we shall denote by H the heart tissue, by Ω_b and σ_b the cavitary blood and its conductivity, by Ω_f and σ_f the fluid volume adjacent to the epicardium and its conductivity, by $\Omega_0=\Omega_b\cup\Omega_f$ the whole extracardiac medium in contact with H and by $\bar \Omega=\bar H \cup \bar \Omega_0$ . Let $\Gamma=\partial \Omega$ , $\Sigma=\partial H \cap \partial \Omega_0$ , i.e. $\Sigma$ represents the parts of the epi and/or endocardium in contact with $\Omega_0$ , and $\Gamma_0=\partial \Omega_0-\Sigma$ , $\Gamma_H=\partial H-\Sigma$ , hence $\Gamma=\Gamma_H \cup \Gamma_0$ and $\partial H=\Sigma \cup \Gamma_H$ . In particular the epicardial surface will be denoted by $\Sigma_{epi}$ .

We define M=M_i+M_e the ``bulk'' conductivity tensor and

$\begin{displaymath}\widehat{M} = \widehat{M}({\bf x})=\left \{\begin{array}{ll......u_0({\bf x},t), & {\bf x}\in \Omega_0\end{array} \right .\end{displaymath}$

with M₀ conductivity tensor in $\Omega_0$ given by $M_0({\bf x})=\left \{\begin{array}{ll}\sigma_f~ I, & {\bf x}\in \Omega_f \\\sigma_b~ I, & {\bf x}\in \Omega_b\end{array} \right.$ where I denote the identity matrix. Let $\Gamma=\partial \Omega$ and $\Sigma=\partial H \cap \partial \Omega_0$ , i.e. $\Sigma$ represents the parts of the epi and/or endocardium in contact with $\Omega_0$ . In particular the epicardial surface will be denoted by $\Sigma_{epi}$ .

The bidomain is characterized by:

a) the local fiber direction ${\bf a}={\bf a}({\bf x})$ and the global fiber structure

b) the conductivity coefficients $\sigma^{i,e}_{l,t}$ along and across fiber in the (i), (e) media. $\sigma^{i,e}_t$ is the same for all directions orthogonal to ${\bf a}({\bf x})$ . The conductivity tensors $M_i({\bf x}),~M_e({\bf x})$ are defined by:

$\begin{displaymath}M_{i,e} ({\bf x}) = \sigma^{i,e}_t I +\left( \sigma^{i,e}_l - \sigma^{i,e}_t \right) \, {\bf a}{\bf a}^T\end{displaymath}$

(1)

where ${\bf a}$ is a column vector and ${\bf a}^T$ its transpose. In the bidomain model the following potentials are considered: $u_i({\bf x},t),~u_e({\bf x},t)$ , $v({\bf x},t)= u_i({\bf x},t)-u_e({\bf x},t)$ and $u_0({\bf x},t)$ representing respectively the intra-, extracellular transmembrane and extracardiac potentials and the current vector densities associated to the intra, extracellular and extracardiac potentials : $\begin{displaymath}{\bf J}_i=-M_i\nabla u_i,~{\bf J}_e=-M_e\nabla u_e ~~ and ~~{\bf J}_0=-M_0\nablau_0\end{displaymath}$

Moreover setting

$\begin{displaymath}{\bf J}_v({\bf x},t)=-M_i \nabla v({\bf x},t)\end{displaymath}$

(2)

and applying the current conservation to the bidomain model we obtain :

$\begin{displaymath}\left \{\begin{array}{l}\hbox{div}~( {\bf J}_i + {\bf J}_...... n}^T {\bf J}_0=0 ~~ \hbox{ on }\Gamma_0\end{array}\right .\end{displaymath}$

(3)

Assuming known the transmembrane potential $v({\bf x},t)$ then in terms of the extracellular or extracardiac potential $u({\bf x},t)$ the bidomain current conservation equations (3) are equivalent to :

$\begin{displaymath}\left \{\begin{array}{l}\hbox{div}~\widehat{M}\nabla u({\......_v({\bf x},t)\qquad \hbox{ on }\Gamma_H\end{array}\right .\end{displaymath}$

(4)

where $\left\lbrack\! \left\lbrack \right. \right.\varphi\left.\left.\right\rbrack \! \right\rbrack_\Sigma$ denotes the jump of $\varphi$ through $\Sigma$ , i.e. $\left\lbrack\! \left\lbrack \right. \right.\varphi\left.\left.\right\rbrack \! \right\rbrack_\Sigma = \varphi_{_{\Sigma^+}} - \varphi_{\vert _{\Sigma^-}}$ with $\varphi_{\vert _{\Sigma^\pm}}$ the traces taken on the positive and negative side of $\Sigma$ with respect to the oriented normal. Therefore problem (4) provides the extracellular or extracardiac potential $u({\bf x},t)$ from the knowledge of $v({\bf x},t)$ , see ([2],[4]). As a consequence of the bidomain approach and assuming known the transmembrane potential distribution v, the term $\hbox{div}~{\bf J}_v({\bf x},t)$ plays the role of the ``impressed'' or ``source'' current density generating the extracellular potential. Different mathematical approaches can be used to simulate EGs (see e.g. [6,8,15,21]). Because EGs are usually recorded from a limited number of sites, an integral formulation (see [23,7]) is more convenient than differential or variational representations, for large scale simulations as required for an accurate computation of extracardiac, epicardial or intramural EGs [15,2]. Each EG represents the time course of the potential difference between one point of the domain, called observation point, and a ``reference potential''. For unipolar EGs, the reference potential is that at a remote site or, more frequently, is obtained by averaging the potential values over a set of three or more points or over a surface, e.g. the entire insulated surface of the conducting volume as suggested in [17]. Recent experimental observations with isolated canine hearts in a torso-shaped tank showed that the potential of Wilson's central terminal is close to the average potential [19] .

In this work EGs are computed by using the latter reference potential, i.e. the average potential on the epicardial surface $\Sigma_{epi}$ .

Then we consider

$\begin{displaymath}w({\bf x},t)= u({\bf x},t)-\frac{1}{\vert\Sigma_{epi}\vert} \......Sigma_{epi}}u({\mbox{\boldmath$\xi$ }},t)~d\sigma_{{\!\xi}}.\end{displaymath}$

(5)

For ${\bf x}\in\Sigma$ , we introduce the solution $\psi$ of the following elliptic problem :

$\displaystyle \left \{ \renewcommand{\arraystretch}{1.7}\begin{array}{ll}-\hbo...... n}^T \widehat{M} \nabla_{\!\xi}\psi= 0& \hbox{ on }\Gamma\end{array}\right .$

(6)

with ${\mbox{\Large$\chi$ }}_{\Sigma_{epi}}({\mbox{\boldmath$\xi$ }})$ characteristic function of $\Sigma_{epi}$ , i.e. ${\mbox{\Large$\chi$ }}_{\Sigma_{epi}}({\mbox{\boldmath$\xi$ }})=\left \{\begi......dmath$\xi$ }}\in \Sigma_{epi} \\0, & \hbox{ elsewhere}\end{array} \right.$ and $\left\lbrack\! \left\lbrack \right. \right.\phi \left.\left.\right\rbrack \! \right\rbrack_\Sigma$ denotes the jump of $\phi$ through $\Sigma$ , i.e. $\left\lbrack\! \left\lbrack \right. \right.\phi \left.\left.\right\rbrack \! \right\rbrack_\Sigma= \phi_{\Sigma^+} -\phi_{\Sigma^-}$ with $\phi_{\Sigma^{\pm}}$ the traces taken on the positive and negative side of $\Sigma$ with respect to the oriented normal.

Applying the second Green formula to the couple $(u,\psi)$ in $\Omega_0$ and in H with $\psi$ solution of (6), adding these two relations we obtain the following integral representation:

$\begin{displaymath}w({\bf x},t) = \int_H {\bf J}_v^T ~\nabla_{{\!\xi}} \psi({\...... \nabla_{{\!\xi}} \psi({\mbox{\boldmath$\xi$ }},{\bf x}) ~d\xi\end{displaymath}$

(7)

The boundary condition imposed for the ``lead field'' actually reflects the property that $u({\bf x},t)$ , defined by (5), has zero average on the epicardium.

The same integral representation can be extended for observation points ${\bf x}\in \Omega_0\cup H$ , i.e. within the wall or in the extracardiac medium. In these cases the lead field $\psi$ is defined by:

$\displaystyle \left \{ \renewcommand{\arraystretch}{1.7}\begin{array}{ll}-\hbo...... n}^T \widehat{M} \nabla_{\!\xi}\psi= 0& \hbox{ on }\Gamma\end{array}\right .$

(8)

For more details on the integral representation of w see [2] .

The integral approach easily allows to change the reference potential of the EGs without recomputing them. For instance let us consider, as suggested by Spach et al. in [17], as reference potential the average potential on the whole insulated surface $\Gamma$ . If we denote by

$\begin{displaymath}w_\Gamma({\bf x},t)=\int_H {\bf J}_v^T \nabla \psi_\Gamma ~d\xi\end{displaymath}$

(9)

with $\psi_\Gamma$ solution of:

$\displaystyle \left \{ \renewcommand{\arraystretch}{1.7}\begin{array}{ll}-\hbo......splaystyle -\frac{1}{\vert\Gamma\vert}& \hbox{ on }\Gamma\end{array}\right .$

(10)

then it is easy to verify that $w_\Gamma({\bf x},t)$ has zero average on the whole surface $\Gamma$ . Thus, from (7) and (9) we have

$\begin{displaymath}w_\Gamma({\bf x},t)-w({\bf x},t)= \int_H {\bf J}_v^T \nabla(\psi_\Gamma -\psi)\end{displaymath}$

and, defining $\psi_d=\psi_\Gamma- \psi$ , we obtain that $\psi_d$ is solution of

$\displaystyle \left \{ \renewcommand{\arraystretch}{1.7}\begin{array}{ll}-\hbo......splaystyle -\frac{1}{\vert\Gamma\vert}& \hbox{ on }\Gamma\end{array}\right .$

(11)

In conclusion, defining $w_d(t)= \int_H {\bf J}_v^T \nabla\psi_d ~d\xi$ , we have

$\begin{displaymath}w_\Gamma({\bf x},t)=w({\bf x},t)+ w_d(t)\end{displaymath}$

(12)

hence to change the reference potential fixed for the EGs, i.e. in this case considering the average potential on the whole $\Gamma$ instead of the average potential on $\Sigma_{epi}$ , we have to compute $\psi_d$ and $\int_H {\bf J}_v^T \nabla \psi_d ~d\xi=w_d(t)$ and to add w_d(t) to $w({\bf x},t)$ . Note that $\psi_d$ , solution of problem (11), can be computed only once since it does not depend on the observation point ${\bf x}$ ; moreover $\psi_d$ is a smooth function without singularity, i.e. with no delta function as the lead fields (6) and (10); hence it can be computed using a uniform mesh without refining it near the observation point ${\bf x}$ (see [2] for more details). This procedure can easily be applied to obtain other reference potentials, e.g. the average on a subset of $\Gamma$ or $\Sigma_{epi}$ .

In a previous paper [2], we developed and validated an efficient numerical technique for the computation of EGs, at any point inside or outside the myocardium, free from spurious oscillations and peaks.

We now want to use this numerical procedure to elucidate the mechanisms that produce the frequently observed multiphasic EG waveform; more generally, we want to analyze the factors that influence the shape, amplitude and polarity of unipolar EGs.

Source splitting

We briefly recall the main formulae related to the split form which will be applied to the analysis of the EG wave shape. The ``bulk'' conductivity tensor is defined by M=M_i+M_e. Setting $\sigma_{l,t}= \sigma_{l,t}^i + \sigma_{l,t}^e$ it is easy to verify that:

$\begin{displaymath}M_i = \alpha\displaystyle \left (~M +\beta ~{\bf a}~{\bf a}^T......{\sigma_t}{\sigma_l}\frac{\sigma_l^i}{\sigma_t^i}-1\right )\end{displaymath}$

(13)

and that the case of equal anisotropy ratio $\rho_i=\sigma_l^i/\sigma_t^i=\rho_e=\sigma_l^e/\sigma_t^e$ is characterized by $\beta=0$ . Using (13), since ${\bf J}_v = -M_i \nabla v$ we obtain:

$\begin{displaymath}{\bf J}_v = {\bf J}_a + {\bf J}_c\end{displaymath}$

(14)

where

$\begin{displaymath}{\bf J}_a=-\alpha~\beta\left({\bf a}^T\nabla v\right){\bf a}\quad,\quad {\bf J}_c=-\alpha~ M \nabla v\end{displaymath}$

(15)

We remark that ${\bf J}_a$ , ${\bf J}_c$ are distributed current densities with dipole axes parallel respectively to the fiber direction ${\bf a}$ and to the conormal vector $M \nabla v$ . For this reason ${\bf J}_a$ and ${\bf J}_c$ are called axial and conormal current densities. From (7) it follows that :

$\begin{displaymath}w({\bf x},t)= w_a({\bf x},t)+ w_c({\bf x},t)\end{displaymath}$

(16)

where

$\begin{displaymath}w_a({\bf x},t)=\int_H {\bf J}_a ({\mbox{\boldmath$\xi$ }},t)^T ~ \nabla_{\!\xi}\psi({\mbox{\boldmath$\xi$ }},{\bf x})~d\xi\end{displaymath}$

(17)

$\begin{displaymath}w_c({\bf x},t)=\int_H {\bf J}_c({\mbox{\boldmath$\xi$ }},t)^T ~\nabla_{\!\xi}\psi({\mbox{\boldmath$\xi$ }},{\bf x})~d\xi\end{displaymath}$

(18)

For t fixed and ${\bf x}$ variable (16) defines a split of the potential distribution into the axial and conormal component of the field, while for ${\bf x}$ fixed and t variable (16) defines a split of the EG into the axial and conormal EG waveforms.

We recall that for a bidomain with equal anisotropy ratio we have $\beta=0$ and consequently in our decomposition (16) only the conormal component is present, i.e. the full EG $w({\bf x},t)$ reduces to the conormal EG component $w_c({\bf x},t)$ .

Assuming $\alpha= \sigma^i_t/\sigma_t$ constant in H, i.e. $\nabla\alpha=0$ , which holds if $\sigma^{i,e}_t$ are constant, using the Green formula $\int_H \nablav^T~M\nabla\psi~d\xi= -\int_H v ~\hbox{div}~_{\!\xi}M \nabla_{\!......h$\xi$ }}+ \int_{\partial H} v ~{\bf n}^T M \nabla_{\!\xi}\psi~d\sigma_{\!\xi}$ and taking into account the properties of $\psi$ in (6) we obtain the following further decomposition of the conormal field $w_c({\bf x},t)$ :

$\begin{displaymath}w_c({\bf x},t) = u_J({\bf x},t) + u_d(t) + u_{hs}({\bf x},t)\end{displaymath}$

(19)

where

$\begin{displaymath}\begin{array}{c}u_J({\bf x},t) =\left \{ \begin{array}{l......a v ~{\bf n}^T M \nabla_{\!\xi}\psi~d\sigma_{\!\xi}\end{array}\end{displaymath}$

(20)

where ${\vert\Sigma_{epi}\vert}$ denotes the epicardial surface area.

The component u_J is usually called ``jump component'' because its effect on the conormal potential w_c reduces to a shift proportional to the jump of the transmembrane potential from the resting to the plateau value through the excitation front for t fixed or at the activation time for ${\bf x}$ fixed.

We note that $u_d(t)=\displaystyle \frac{\alpha}{\vert\Sigma_{epi}\vert}\int_{\Sigma_{epi}} v ~d\sigma$ is independent of the observation point ${\bf x}$ where the EGs are simulated, positive and increasing with time. Therefore, it needs to be computed only once and can be considered as a drift component related to the chosen reference potential; we call this term the drift of the conormal component of the EG.

We remark that $u_{hs}({\bf x},t)$ contributes to the conormal potential $w_c({\bf x},t)$ only when the heart is in contact with an extracardiac conducting medium (i.e. $\Sigma\neq \emptyset$ ) and depends on the trace on $\Sigma$ of the transmembrane potential $v({\bf x},t)$ . Since the component (20) represents the potential generated by a conormal dipole layer lying on the non-insulated heart surface $\Sigma$ we refer to this field component as usual (see e.g. [6,7,15]) as heart surface source model, i.e. this component is the potential generated by a double layer on $\Sigma$ with dipoles parallel to $M {\bf n}$ and moment proportional to the transmembrane potential v.

In conclusion we get:

$\begin{displaymath}w({\bf x},t) =w_a({\bf x},t) + u_J({\bf x},t) + u_d(t) + u_{hs}({\bf x},t)\end{displaymath}$

(21)

For a derivation of the split form under more general conditions and for more details see [] and [2]. In the special case of H fully insulated (i.e. $\Gamma_H=\partial H$ and $\Sigma=\emptyset$ ) and $\alpha$ constant the conormal component reduces to:

$\begin{displaymath}w_c({\bf x},t)= u_J({\bf x},t) + u_d(t) = -\alpha v ({\bf x},......_{\Sigma_{epi}} v({\mbox{\boldmath$\xi$ }},t) ~d\sigma_{\!\xi}\end{displaymath}$

(22)

Computation of $v({\mbox{\boldmath$\xi$ }},t)$

The integral formulation (7) requires the knowledge of the time-space distribution of the transmembrane potential $v({\mbox{\boldmath$\xi$ }},t)$ . Our numerical procedure can be used to compute EGs over the whole cardiac cycle from ventricular excitation to recovery; however in this paper we limit our simulations to waveforms associated with the depolarization phase, i.e. to the so called QRS complex of EGs.

We use the eikonal model (see [1,3]) to build an approximation of the transmembrane potential $v({\mbox{\boldmath$\xi$ }},t)$ ; in this approximation, assuming uniform membrane excitability properties, $v({\mbox{\boldmath$\xi$ }},t)$ is given by:

$\begin{displaymath}v({\mbox{\boldmath$\xi$ }},t)=\hbox{\Large\it a}(t-\varphi({\mbox{\boldmath$\xi$ }})).\end{displaymath}$

(23)

where

ii): the activation time $\varphi({\mbox{\boldmath$\xi$ }})$ is obtained by solving the eikonal equation (see []) describing the motion of the excitation wave front;
iii): the action potential upstroke is given by:

(24)

with C_m membrane capacitance per unit area, v_r, v_p, v_thresting, plateau and threshold value of the transmembrane potential v, Gmaximum membrane conductance per unit area and s = (v_p-v_r)/(2(v_th-v_r))-1 (see []) .

Numerical approximation

The numerical computation of the integral (7) is carried out by FEM and requires special care because of the singularity of $\psi$ at the observation point ${\bf x}$ and the presence of the moving steep wave front characterizing $v({\mbox{\boldmath$\xi$ }},t)$ during the depolarization phase.

We sketch here the main steps of the procedure used to improve the accuracy of $w({\bf x},t)$ :

b): refinement of the Finite Element Method mesh in computing the lead field $\psi({\mbox{\boldmath$\xi$ }};{\bf x})$ near its singular point ${\bf x}$ . Since we know that $\psi({\mbox{\boldmath$\xi$ }},{\bf x})$ is a smooth function far from ${\bf x}~$ we build a grid having small elements in proximity of the singularity of $\psi$ and elements of increasing dimensions further away. A structured grid refinement around the singularity was performed maintaining approximately the same number of elements of the mesh associated to the computation of the activation time $\varphi({\mbox{\boldmath$\xi$ }})$ so that the computer time and memory remain unchanged.
c): at each time instant, identification of the elements crossed by the wave front and successive subdivision of these elements;
d): nonlinear interpolation of $v({\mbox{\boldmath$\xi$ }},t)=\hbox{\Large\it a}(t-\varphi({\mbox{\boldmath$\xi$ }}))$ ; The activation time $\varphi({\mbox{\boldmath$\xi$ }})$ , a smooth function without sharp variations, is well approximated on each sub-element using a trilinear polynomial $\varphi_h({\mbox{\boldmath$\xi$ }})$ defined by the values of $\varphi$ at the nodes of the macro hexahedral element. Hence using this trilinear interpolation $\varphi_h$ by means of (23) we are led to consider, for the function $v({\mbox{\boldmath$\xi$ }},t)$ exhibiting high gradient, the following nonlinear, but more accurate, interpolation formula:
e): use of the split form (16) of $w({\bf x},t)$ ; The axial and conormal component w_a, w_c can be computed by exploiting the procedure outlined in steps b), c) for the evaluation of the volume integral (17) and (18) respectively. Moreover, the accuracy of the computation of w can be increased in the case of H fully insulated by combining the split form (22) with the numerical technique outlined above. In fact, from (22), the conormal component $w_c({\bf x},t)$ can be evaluated by computing the surface integral $u_d(t)=\displaystyle \frac{\alpha}{\vert\Sigma_{epi}\vert}\int_{\Sigma_{epi}} v ~d\sigma$ and the jump $u_J({\bf x},t)=-\alpha v({\bf x},t)$ by using (23). For the evaluation of the surface integral we used steps b), c) adapted to surface elements. In conclusion in the fully insulated case we apply the sub-element technique to evaluate the volume integral (17) yielding $w_a({\bf x},t)$ and we compute $w_c({\bf x},t)$ using (19).

In this way it is possible, by using the outlined numerical procedure, to obtain EGs free from numerical artifacts like spurious oscillations or peaks, with a limited computational effort, as established in [2].

Results

Experimental EGs showed that within 10-15 mm from the pacing site the shape of unipolar EGs markedly depends on whether excitation reaches the recording site by traveling mainly along or across fibers [16,20,21]. In the following, these recording sites will be referred to as along-fiber or cross-fiber sites. Fig. 1 shows the morphological differences in epicardial EGs recorded from 4 sites located just beyond the 18 msec isochrone, along longitudinal (along-fiber) or transverse (cross- fiber) pathways. Due to higher longitudinal conduction velocity, the 18 msec isochrone is elongated, quasi-elliptical, centered about the pacing site, with major axis oriented along the local direction of epicardial fibers.

Figure 1. Isolated dog heart immersed in an electrolytic tank shaped as a 10-year boy's torso, filled with a 500 $\Omega$ cm solution. CX = circumflex coronary branch. LAD = left anterior descending coronary branch. Ventricle was paced from star-marked site. Elliptical contour is the 20 ms isochrone. Right side : unipolar EGs were recorded from sites denoted by large dots near the epicardial elliptical isochrone at 18 msec after pacing . Left side: computed EGs in monoventricular model with the epi- and endocardium in contact with conducting media. Epicardial EGs No 1 through No 4 were computed for sites marked by dots around the simulated epicardial isochrone at 20 msec after stimulus.

We define the EGs as monophasic, biphasic etc. according to the number of relative maxima and minima that occur before and immediately after the intrinsic deflection (the major downstroke).

The EGs displayed in Fig. 1 show that:

the cross-fiber EGs are monophasic. The portion preceding the intrinsic deflection is monotonically negative-going, as previously described in [19,20,21];
the along-fiber EGs are biphasic, with an initial positive R-wave followed by the negative-going intrinsic deflection [22].

At later stages of ventricular excitation, when approximately 50 $\%$ of the ventricular surface has been excited, the experimental EGs, recorded from sites close to the 60 ms isochrone, (Fig. 2 ), are biphasic, with the positive and negative phases having almost equal magnitude.

Figure 2. Measured and computed epicardial electrograms. Measured and simulated epicardial isochrones at 60 and 80 msec after the stimulus are traced in the right and left Panel respectively. Same layout and explanations as in Figure 1.

For a positive and negative R and S wave respectively we introduce the ratio between the R wave amplitude and the pick to pick drop given by: $\displaystyle \frac{R}{R+\vert S\vert}$ . The progressive positive trend of the $\displaystyle \frac{R}{R+\vert S\vert}$ wave amplitude ratio, which reaches its maximum value in Fig.2, is an effect of the drift of the reference potential, which moves from near the positive to near the negative extreme of the potential range during ventricular excitation, as recently recognized and discussed in our previous work [19,21]. These studies showed that all unipolar EGs recorded directly from the heart, when referenced to Wilson's Central Terminal or to the average epicardial potential, contain a positive-trending component, associated with the drift of the reference potential. This component causes the $\displaystyle \frac{R}{R+\vert S\vert}$ ratio to increase progressively from zero to 1 as a function of the amount of excited epicardial area.

As a result, EGs recorded from areas that are excited at the end of the QRS interval, near the extinction region, are entirely positive (Fig. 3 ).

Figure 3. Measured and computed epicardial electrograms. Measured and simulated epicardial isochrones at 85 and 110 msec after the stimulus are traced in the right and left Panel respectively. Same layout and explanations as in Figure 1.

Proceeding from the pacing site along cross-fiber pathways, the EGs shift from monophasic to tetraphasic. These multiphasic waveforms occur mostly at the epicardial and subepicardial levels in experimentally recorded EGs, in regions where excitation spreads across fibers. In regions located outside the 40 msec isochrone (see Fig. 4):

the along-fiber EGs still have a biphasic shape with an R wave preceding the intrinsic deflection.
the cross-fiber EGs show a more complex time course; in fact we have:

EGs that, after an initially negative- going phase, show a minimum followed by a maximum before the intrinsic deflection (Fig. 4 A, EG No 1).
EGs that, after a flat initial phase, show a positive- negative wave followed by a hump. The hump remains below the baseline (Fig. 4 A, EG No 2, arrow) or, at greater distances from the pacing site, becomes a positive spike (Fig. 4 A, EG No 3, arrow). Along-fiber EGs still have a biphasic shape with an R wave preceding the intrinsic deflection.

Figures 1,2,3,4 display the computed EGs, which were simulated by applying a local epicardial stimulus to the monoventricular model incorporating the Purkinje network; these EGs were obtained from sites which are similarly located, relative to the relevant isochrones, to the experimental recording sites shown in Figures 1-4. The entire variety of morphologic features observed in measured EGs was also found in the computed EGs and the agreement was maintained over the entire ventricular surface. The main similarities relate to:

a) the simple and almost invariant biphasic shape of along-fiber EGs;

b) the evolving complexity of EG waveforms recorded at increasing across-fiber distance from the pacing site;

c) the biphasic shape in areas excited at approximately 60 ms and the totally positive QRS complexes in the extinction regions.

Some discrepancies between experimental and simulated EGs have been observed. For instance the QRS duration of the experimental and simulated EGs was not expected to be the same since the model geometry used concerns only the left ventricle and does not match the actual dog heart dimensions. Moreover the model incorporates a simplified Purkinje network that covers only 1/3 of the left ventricular endocardium. The septum has no Purkinje network. This can explain the longer simulated QRS duration, compared with the experimental ones. Other discrepancies may be due to simplifying assumptions in the model, e.g. the flat plateau of the action potential, which does not reproduce the initial fast repolarization typical of epicardial fibers. In all previous simulations we used the average epicardial potential as the reference potential. If we choose a different reference potential, e.g. the average of the potential on the epi and endocardium or on all the insulated boundary, we produce only a slight modification in the EG wave forms, which does not alter the multiphasic structure of EGs (not shown).

Figure 4. Panel A : isolated dog heart immersed in an electrolytic tank shaped as a 10-year boy's torso, filled with a 500 $\Omega$ cm solution. CX = circumflex coronary branch. LAD = left anterior descending coronary branch. Ventricle was paced from star-marked site. Elliptical contour is the 40 ms isochrone. Unipolar EGs No 1 through No 4 were recorded from sites 1 - 4 (large dots near the epicardial elliptical isochrone at 20 msec after pacing ). Panel B: computed EGs in monoventricular model with the epi- and endocardium in contact with conducting media. Epicardial EGs No 1 through No 4 were computed for sites marked by dots around the simulated epicardial isochrone at 50 msec after stimulus.

Interpretation of EG waves With a view to interpreting the multiphasic shape of the QRS complexes in terms of intracardiac excitation events, we simulated EGs elicited by epicardial pacing of a fully insulated ventricle, in order to eliminate the effect of the conducting media on potential fields and EGs. One significant feature of the fully insulated model is that, in the split procedure that decomposes the electric sources into an axial and a conormal component, only the axial component of the electrical sources associated with the wave front contributes to the extracellular potential field, while the conormal component is silent (except for the potential jump associated with the wave front), as predicted by the classical uniform dipole layer model . A selection of epicardial EGs is displayed in Fig. 5; the two sets of EGs reported in Panels A and B were obtained from 6 cross-fiber and 6 along-fiber sites, located at increasing distance from the pacing site along two diagonal directions.

Figure 5. Monoventricular model fully insulated. Central epicardial stimulation. Panels A, B : superimposed, computed epicardial electrograms from points marked by "A" and "B" dots in Panel D. Dots are located on semidiagonal pathways along and across fiber respectively, at increasing distance from the stimulation site. Panel C : four superimposed epicardial electrograms related to observation points marked by dots on a vertical line (labeled C) in Panel D ; in this region, isochrones exhibit a dimple-like inflection. Panel D : simulated excitation time map depicting the spread of excitation on the epicardial surface. Time interval between successive isochrones is 5 msec.

The main features of the EGs are:

1.: a biphasic wave form for along-fiber sites. More precisely, the EGs display a positive R wave followed by the intrinsic deflection, see Fig. 4 Panel A .
2.: proceeding along the cross-fiber diagonal (see Fig. 5 Panel B from left to right), we observe at first, near the pacing site, monophasic EGs (the first and second EGs) which become progressively triphasic (the third EG) and then pentaphasic (the fourth and fifth EGs), because the initial portion exhibits a quite weak negative minimum (-0.1 mV) (phase 1), see the EG enlargement in Fig. 4 Panel E, and subsequently shows 4 additional phases: a positive hump (phase 2) followed by a negative-going segment, a minimum (phase 3) and a spike (phase 4) before the intrinsic deflection that terminates with a negative spike (phase 5). Finally the last (sixth) EGs is tetraphasic because the initial weak minimum is missing.

Notice that pentaphasic EGs are found, Fig. 5 Panels B, C, in those cross-fiber areas where the epicardial isochrones show a dimple-like inflection i.e. a change of curvature. At increasing distances from the pacing site, these EGs show, before the intrinsic deflection, a spike of increasing magnitude as the dimple-like inflection becomes more pronounced.

These results show that the main features of the measured and computed EGs, previously described for hearts surrounded by conducting media (see Figs. 1-4), are also present in the case of a fully insulated ventricle; therefore, most of the polyphasic behavior of the EGs cannot be attributed to the influence of the interfaces with the extracardiac conducting media. The media, however, do modify the EGs to some extent, as shown later in this Section.

The initial positive R wave, which is characteristic of along-fiber EGs, is the electrographic counterpart of the positive areas facing those portions of the wave front that propagate mainly along fibers (see Fig. 8 Panel A, site 2). A site reached by the wave front as it moves along fibers exhibits increasingly positive potentials and records a positive R wave (Fig. 8 Panel B); when the wave front crosses the site, the intrinsic deflection is recorded. The subsequent rising phase, where the tracing returns to the baseline, is due to the drift of the reference potential, as previously discussed in [19,21].

The negative-going initial deflection in EGs from cross-fiber points near the pacing site is consistent with the epicardial potential distribution in the relevant area, which shows an increasing negativity for an observer that moves from the resting region toward the wave front. More details of the potential patterns in the early stages of propagation can be found in our previous articles [20,4].

If we record EGs at increasing transverse distances from the pacing site, we see that the initial negative-going deflection (Q wave) is progressively replaced by a flat segment followed by a positive- going hump, whose amplitude increases with increasing distance from the pacing site (see Fig. 6 Panel B).

The evolving hump is in part an epicardial reflection of the rotating deep positivity generated by those intramural portions of the wave front which propagate along deep fibers, whose direction rotates CCW from epi- to endocardium. In fact at all depths, two maxima are observed near those portions of the deep wave front that propagate along the direction of deep fibers; this produces two helical positive potential volumes that develop CCW from epicardium to endocardium and are associated with the twisted left-handed helical shape of the wave front (see [3,4]). The epicardial reflection of this deep positivity produces the first, smooth positive hump (phase 2) in cross-fiber EGs. The rule that R waves in unipolar EGs are associated with an excitation wave approaching the observation point (see [22]) invariably applies to locations reached by an excitation wave front that moves mainly along fibers, but does not apply for cross-fiber points near the pacing site . However, at several cm from the pacing site, in some cross-fiber locations, a spike was observed following a first, smooth hump (see Fig. 6 Panel C). Both the hump and the spike became more and more positive with increasing distance from the pacing site.

To ascertain whether this spike was the epicardial manifestation of a localized axial current source, arising where the wave front is not parallel to the local fiber direction and therefore the current sources associated with the wave front have an axial component we analyzed the morphology of the split form of the EGs.

Figure 6. Components of the split form of the electrograms in the fully insulated monoventricular model. The epicardial potential is used as the reference potential. The axial and conormal electrograms are displayed in the left and right column respectively In first, second and third row are reported the electrograms from sites A, B and C in Figure 6, Panel D, rspectively.

Columns in Fig. 6 show the axial and conormal EG components for the wave trains of the full EGs displayed in Fig. 5 A, B and C. For all cross-fiber and along-fiber observation sites, the conormal EGs have a biphasic shape with a positive R-wave followed by an identical downstroke, see the right column of Fig. 6. The positive trend (R-wave) preceding the downstroke increases in magnitude as a function of increasing excitation times, that is, for sites at increasing distances from the stimulus location.

From the theoretical considerations about the further split of the conormal EG developed in Methods, it follows that, for a fully insulated ventricle, the conormal waveform is the superposition of the jump and of the drift component of the conormal EG, see equation (22). The drift is given by ${u}_d(t)=\frac{\alpha}{\vert\Sigma_{epi}\vert}~\int_{\Sigma_{epi}} v d\sigma$ for our reference potential; because the upstroke of v approximates a step function, the conormal drift is nearly proportional to the area of the excited epicardial surface at time t, therefore it is the same for all sites, positive and monotonically increasing as a function of time. The addition of this drift to the jump EG component results in the biphasic structure of the conormal waveforms; in fact, from the right column of Fig. 5, it is evident that the conormal EG consists of a moving downstroke superimposed to a positive trend which is the same for all recording sites.

It follows from the analysis of the conormal EGs that an insulated bidomain structure, with equal anisotropy ratio between the intra and extracellular media, can yield only biphasic EGs, since the axial component vanishes, and we are left with the conormal component alone. This contention is confirmed by EG waveforms published by Simms and Geselowitz (Fig. 4 in [15]). Another interesting consideration is that in such a structure there is no current flowing except within the wave front itself since from equation (22) $\nabla w_c$ results proportional to $\nabla v$ . Therefore, in a fully insulated heart, unipolar EGs that deviate from a biphasic waveshape indicate that the underlying bidomain has unequal anisotropic ratio.

The axial components of the EGs (Fig. 6, left column), have widely different shapes, depending on the geometrical interrelationships between the propagating wave front and the local fiber direction at the various observation points. The first row of Fig. 6 left column, relates to along- fiber sites and shows biphasic EGs. This configuration is maintained and emphasized when the biphasic conormal EG is added in order to obtain the full EG. The axial waveforms for cross-fiber sites are displayed in the left column of Fig. 6 second and third row. These EGs demonstrate that the multiphasic shape of the full EGs originates from the axial waveforms because the full and axial EGs share the same polymorphic configuration from mono- to pentaphasic EGs with a positive hump followed by a spike. Adding the conormal to the axial component to produce the full cross-fiber EGs increases the amplitude of the intrinsic deflection, damps the negative- going trend after the hump and enhances the monotonically increasing trend during and after the spikes.

We recall that the axial current density ${\bf J}_a=-\alpha~\beta~\left( {\bfa}^T \nabla v \right) {\bf a}$ strongly depends on $\vert {\bfa}^T{\bf n}\vert$ since $\nabla v$ results parallel to the unit vector $\bf n(\bf x)$ normal to the excitation wave front surface at point x; as a result of this geometric interaction between the direction of propagation and the local fiber direction , especially for cross-fiber sites, the axial EGs exhibit multiphasic features.

In the more realistic conditions, in which the epi- and endocardial surfaces are in contact with fluids, the multiphasic EGs, with humps and spikes, are still present, although the magnitude of the epicardial EGs is greatly attenuated relative to the insulated case.When the heart is in contact with a conducting medium, the heart surface component u_hs(t) arises and produces its own contingent of currents. Extension of the previous analysis of the individual EG waveshape can be performed based on the split form and shows that the multiple phases dispayed by the cross-fiber EGs derive mostly from the axial EGs because the conormal EGs, a part from an initial weakly negative-going phase, displayes a biphasic behaviour. In the model with fluids, too, the presence of the drift associated with the reference potential prevents us from interpreting the spike as a characteristic marker of an underlying axial current source. To check whether the spike is the expression of a local current source we must first remove the drift from the EGs.

Finally the splitting analysis could be applied to the EGs wave shape after removing the reference or drift component that reflects the motion, during QRS, of the reference potential from near the positive to near the negative extreme of the potential range.

Conclusions

The aim of this study was:

1.: to assess whether our monoventricular model, based on the anisotropic bidomain representation of the heart muscle, with unequal anisotropy ratio, was able to reproduce the variety of mono- to multiphasic unipolar electrograms recorded from real hearts.
2.: to investigate the electrophysiology of the multiphasic unipolar EGs.

With regard to point 1) the simulated results showed that the bidomain model, with unequal anisotropy ratio, accurately reproduces the great variety of multiphasic shapes experimentally recorded both near the pacing site and at distance from the pacing site.

Concerning point 2), it was shown in [21], that some important aspects of the EGs can be explained by considering them as the sum of two components: i) a field component that reflects the changes of the local potential produced by an approaching or receding wave front. ii) a reference or drift component that reflects the motion, during QRS, of the reference potential from near the positive to near the negative extreme of the potential range. Interpretation of the mechanisms that generate the different wave forms, including the humps and spikes in the multiphasic EGs, was greatly facilitated by using the split form of the model. This enabled us to distinguish the EG features that reflect the perturbations of the potential field brought about by an approaching or receding wave front from the features that reflect the variable position of the reference potential within the potential range (drift component). Removing the drift should be possible to re-establish consistency between the EG wave forms and the potential distributions and suppressed artificial spikes produced by the drift. Further work will extend the analysis to EGs recorded from intramural and subendocardial locations, and also to beats produced by intramural stimulation and spontaneous sinus rhythm over the entire heart beat including also the repolarization phase.

References