I
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. VOL. 39. NO. 12. DECEMBER 1992
1232
A Model Study of Electric Field Interactions Between Cardiac Myocytes HeWC Hogues, L . Joshua Leon, and Fernand A . Roberge, Senior Member, IEEE
Abstract-The transmission of excitation via electric field coupling was studied in a model comprising two myocytes abutted end-to-end and placed in an unbounded volume conductor. Each myocyte was modeled as a small cylinder of membrane (10 pm in diameter and 100 pm in length) capped at both ends. A Beeler-Reuter model modified for the Na’ current dynamics served to simulate the membrane ionic current. There was no resistive coupling between the myocytes and the intercellular junction consisted of closely apposed pre- and post-junctional membranes, separated by a uniform cleft distance. The membrane current crossing the prejunctional membrane during the action potential upstroke tends to flow out of the cleft, but it is partly prevented from doing so by the shunt resistance constituted by the cleft volume conductor. The prejunctional upstroke gives rise to a pulse of positive potential within the cleft which induces a small capacitive current across the post-junctional membrane to yield a small positive change in the intracellular potential in the post-junctional cell. The net result is an hyperpolarization of the post-junctional cleft membrane and a slight depolarization of the rest of the cell membrane since the extracellular potential outside of the cell is zero. The magnitude of this depolarization is quite small for a flat junctional membrane and it can be increased by membrane folding and interdigitation, so as to increase the junctional membrane area by a factor of 10 or more. Even then the post-junctional depolarization does not reach threshold when the extracellular potential around the post-junctional cell is effectively zero. Threshold depolarization occurs in the presence of a large decrease of post-junctional load, by increasing the junctional membrane capacitance andlor decreasing the volume of the post-junctional cell. Assuming that the normal resistive coupling between two cardiac myocytes is 1-4 MQ, our model study indicates that electric field coupling would then be about two orders of magnitude smaller. However, substantial enhancement of the efficacy of electric field transmission was observed in the case of cells with substantial junctional membrane folding.
I. INTRODUCTION HE spread of electrical activity in cardiac tissue involves a length constant which is much longer than a single myocyte, implying a sufficiently low resistance between neighboring cells to allow the flow of local circuit current. Anatomically cardiac tissue is made up of tightly packed rod-shaped cells appraximately 100 X 10 pm [ 11, [2]. The cells appear to be separated at their ends by in-
T
Manuscript received November 12, 1991; revised March 30, 1992. This was work supported by the Medical Research Council and The Natural Sciences and Engineering Research Council of Canada, and the Ministere de I’Enseignement Superieur et de la Science du Quebec. Dr. Leon was also supported by the Whittaker Foundation. The authors are with the Institut de Genie Biomedical, Ecole Polytechnique and Universite de Montreal, Montreal, P.Q. H3C 357, Canada. IEEE Log Number 9204046.
tercalated disks which define an intercellular gap which is on average somewhere between 8 and 20 nm [3]. These gap junctions contain tiny channels which permit the flow of ions between cells and provide the primary resistive electrical coupling [4].Studies of pairs of isolated ventricular myocytes have shown that the junctional resistance, Rj, linking the cytoplasm of the two cells is of the order of 1-4 M Q and that it is not changed by the transmembrane potential or the trans-junctional potential during the conduction of an action potential [5], [ 6 ] . Capacitive coupling between cardiac cells has also been postulated [7], based on observations that the junctional complex displays a time constant that is about 10% of the surface membrane time constant [8]. It has been proposed that the anatomical basis for the existence of a junctional capacitance could be the very close apposition of the preand post-junctional membranes in the area of the gap, combined with an extensive folding and a high degree of interdigitation of these membranes. The possibility that the connective tissue surrounding the cells may have capacitive properties has also been considered [7]. No junctional capacitive effects have been reported in studies of coupled isolated myocytes [6]. Previous model studies of the putative electric field interactions between cardiac myocytes have assumed that the pre- and post-junctional membranes were ordinary excitable membranes, and that the junctional resistance Rj was infinite [9].We maintain these assumptions here in order to examine the possibility of electrical field coupling in the absence of ionic local circuit current flow through R,. Instead of deriving an equivalent circuit to model the junction we use a model developed from electric field theoretic principles. Our model includes a threedimensional construct of the cell system, the conductivities of the intra- and extracellular milieu, and a HodgkinHuxley type representation of the membrane dynamics. Our purpose is to examine how the field coupling is affected by the physical characteristics of the intercellular gap, namely its size and the extent of the membrane folding, and the possible development of excitation in the post-junctional cell when an action potential reaches the junction. 11. METHODS Our study is based on a two-cell assembly lying in a large volume of conducting fluid. Each cell is modeled as a closed surface which separates a homogeneous internal
0018-9294/92$03,00 0 1992 IEEE
I
1233
HOGUES et a l . : FIELD INTERACTIONS BETWEEN CARDIAC MYOCYTES
current. The transmembrane potential is the difference between the intracellular potential 4j and extracellular potential + e , i.e., V = c$i The divergence of the current density is zero everywhere in the internal and external media. For the case of a single cell, we can write [lo]:
........._. >
n ’
,
JUNCTION
++
(a)
subject to boundary conditions
d .................. ~................. ~
% = -1 1, 87
(4)
0,
where 9 is the outward normal to the surface element dS. Green’s theorem [ 111 gives an expression for the potential at a point p on the inside or outside of the membrane surface S as follows: SIDE VIEW
(b)
END VIEW
4,(p)
=
--1
27r s a7 r
1 + 27r -
-+
s
@ ,rr ., d7 S
(6)
where 7 is the vector joining the source point to the field point, r is its length, and dS refers to the surface of integration S. This approach can be extended to consider n disjoint , Si,* , S”, and to include cells with surfaces S’, * a point source of strength Zstim in the intra or extracellular medium. We then have
-
--
(c) Fig. 1. Diagram of a two-cell system without resistive coupling. (a) Cell 1 of length LI and cell 2 of length L2 are aligned axially and have the same radius r . The junctional membranes are uniformly separated by a distance d. Stimulation is applied on the proximal cap of cell 1. (b) Pattern of junctional membrane folding, in the form of regular concentric peaks and troughs. The side view shows a longitudinal cross-section along the diameter of the cells where s is the radial distance between two adjacent peaks or troughs, and a is the amplitude of folding. We define the pattern in terms of the number of full cycles, which is 4 in this case with s = 2r/4. (c) Perspective view of the folded membrane at the junction.
medium of conductivity ai from an unbounded homogeneous external medium of conductivity a,. The two cells are abutted end-to-end to form a narrow cleft of volume conductor separating the prejunctional and post-junctional end-caps of excitable membrane (Fig. 1). We first derive a general formulation of the problem in relation to an arbitrary number of cells, using a boundary element method. Next we specify the parameters of the model and the simulation approach.
A. Boundary Element Method The transmembrane current i, (pA/cm2) is the sum of the ionic and capacitive currents Cm dV/dt + iion= i, (1) where Cm (pF/cm2) is the membrane capacitance, V(mV) the transmembrane potential, and iion(pA/cm2) the ionic
where 4: and 4; denote the potential distribution on the inside and outside surface of cell j . We discretize each surface Si into mi patches {S”, . . . Si” , Si,’}, and assume that &, 4,, and iion are constant, on any given patch S’”. Denoting these quantities by r#$’, 45, and i$’,,, then the system can be rewritten as 9
9
- -
(9)
I234
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 39. NO. 12. DECEMBER 1992
where R is the distance between the observation point and the source and a is the conductivity of the milieu. The average rise in potential caused by such a source on surface element A' is
JA, R
dS
(10) where d Q = V( 1 / r ) VdS is the solid angle subtended by dS at the center of patch S"'. Because the terms inside the integrals depend only on geometry they can be integrated to give a system of 2N equations for the potentials on either side of the N membrane patches, at each instant in time [ 101. The system can be written succinctly in matrix form. Ordering the potentials and currents in column vectors:
4 4.7
=
2na
1
(13)
dS
AI
Consider now a small resistive channel linking the intracellular media of two neighboring cells. The current crossing the channel is proportional to the potential difference between the inner surface of the two cells at the extremities of the channel. The current flowing across a channel of resistance R, joining cells k and i , at patches SkJ and Si'is
This current constitutes intracellular point sources of strength - I , on patch Sk*' and +Z, on patch Si'. The effect of an intracellular connection is included by introducing point sources of opposite signs at the two interconnection sites.
Equations (9) and (10) can be written as a system of first order non-linear ordinary differential equations of the form:
a+
M , - + M29 at
+ M31ion
= 0
(11)
where MI, M2,and M3 are matrices calculated from (9) and (10) (see Appendix A). Replacing the time derivative by a first order implicit finite difference with t, = ti-, At, the matrix system yields:
B. Simulations The above method allows to model several cells of virtually any geometry. For simplicity we represent each cell as a cylinder of membrane of radius r , capped at each end either by a disk of membrane or a pattern of folded membrane. The folding pattern chosen transforms the flat-disk gap junction into a set of two interlocking rippled surfaces of rotation. It is obtained by rotating the curve a 2
+
+(ti) = (MI
+ AtM2)-I(MI@(ti-l)
- AtM3Zio,).
(12)
Equation (12) allows an iterative solution for the potential distribution from an initial condition 9(0). Each element of the current vector Zion is obtained from the membrane at time model using the transmembrane potential (ai ti - I . The matrix (MI + AM,) is nonsingular and can be decomposed in LU form. Solution vectors are then obtained from successive forward and back substitutions, The resulting system of 2N coupled ordinary differential equations can be integrated in time to yield the potential on each side of the membrane. Point Sources and Intercellular Resistive Coupling: An important advantage of this approach is the possibility of locating current sources in either media, or on the boundary [ l o ] , [ 1 2 ] , [ 1 3 ] . The potential created by a point source of current of intensity I, is
4
4s = 2naR
I
+ a-2 cos nnp R
where a / 2 is the amplitude of the ripples, R the cleft radius, p the distance from the cell axis, and n the number of ripples. Hence the amplitude of the ripples is zero at the edge of the cleft and varies between + a / 2 and - a / 2 along the diameter. We define s as the length of a full cycle, that is the distance between two successive peaks or troughs along a diameter, Fig. 1 (b). The number of full cycles is also equal to n , so that s = 2 R / n . The surface area of the gap is S(a, n) =
sd
nRwJ4R2
+ a2n2n2sin (nnw) du
(16)
where R is the radius of the cell and w = p / R . Since it can be shown that S ( a / k , kn) = S ( a , n ) , there is flexibility in choosing the parameters a and n to obtain a given surface area. A perspective view of the folded membrane at the junction is shown in Fig. l(c). For example, with a = 1 pm and n = 13 (or alternatively with a = 0 . 5 pm and n = 26) the membrane area of one end-cap is in-
HOGUES et a l . : FIELD INTERACTIONS BETWEEN CARDIAC MYOCYTES
creased by a factor of about 9 in comparison with a flat disk of the same radius. Such a figure is in agreement with estimates of membrane folding at ventricular myocardial gap junctions [ 141. The cylinder surface was discretized into cylindrical patches of 2.5 pm in length. The flat disk end-cap with r = 5 pm was divided into 12 patches of roughly the same area by concentric circles. When the membrane area was increased by folding or by increasing the radius, the number of end-cap patches was increased proportionally. There is a constant separation of thickness d between the pre- and post-junctional membranes, Fig. 1 , irrespective of the shape of the end-caps. The cleft thickness in ventricular myocardium has been estimated to be on average between 8 and 20 nm [3]. The cleft might very well be narrower on the cell axis, however, where maximum confinement of the prejunctional membrane current would occur. For this reason the nominal value of d = 1 nm was chosen here to insure that maximum conditions of confinement were considered. The Beeler-Reuter model of the cardiac action potential [lS], as modified by Drouhard and Roberge [16], was used to describe the transmembrane current. An action potential was elicited in cell 1 by applying a current to the membrane of the proximal end-cap. It was carefully verified that the spatial and temporal discretization was fine enough to ensure that no numerical distortions occurred during propagation when nominal parameter values (Table I) were used. The time step used in all simulations was constant at 2 ps. Current and potential values were written to a file every SO ps. The method described above was programmed in Fortran on a Silicon Graphics 4D/25 workstation (Silicon Graphics, Mountain View, CA). All simulation were carried out without a resistive interconnection between the two cells, i.e., with Rj = 03, in order to emphasize the role of electric field interaction. Variables q5i, q5e, I,,,, and I,,, were monitored on all patches of the two junctional end-caps, and also on the cylindrical membrane at the periphery of the cleft. The mean values of potential and current in the gap (8e,j,, ,Im2) were computed from the value on each patch, multiplied by the number of patches, and divided by the total membrane area of the end-cap. 111. RESULTS As the action potential elicited on cell 1 reaches the junction, the transmembrane current crossing the prejunctional membrane generates changes in the extracellular potential within the cleft. There results post-junctional transmembrane potential changes which may or may not be sufficient to trigger an action potential in cell 2.
Extracellular Potential in the Cle@ on the cell axis, We consider first the cleft potential midway (d/2) between the two junctional membranes (Fig. 1). Current flows across the prejunctional membrane during the action potential upstroke and produces a smooth
1
1235
TABLE I NOMINAL PARAMETER VALUES Membrane capacitance, C,,,(pF/cm’) Cell radius, r (pm) Length of cell I , L , (pm) Length of cell 2 , L, ( p m ) Intracellular conductivity, U, (mu/cm) Extracellular conductivity, U , (mvlcm) Cleft thickness, d (nm) Amplitude of folding, a (pm)
1 .o 5 1000 IO0 5
66.66 1 1
monophasic waveform, as depicted in Fig. 2(a) for cleft thicknesses between 1 and SO nm. Since the current tends to flow radially out of the cleft, the cleft space can be represented by a shunt resistor (Rsh)whose value is directly proportional to the length of the radial pathway for current flow, and inversely proportional to the cleft thickness. With flat junctional membranes (no folding), the pathway length is equal to the radius of the cleft. Fig. 2(a) shows that a large Rsh, corresponding to d = 1 nm, allows a peak potential of nearly 4 mV. But q5e falls rapidly as d increases and becomes negligible at d = SO nm in this case. This behavior of peak q5c with increasing cleft thickness is summarized in Fig. 2(b). Still without junctional membrane folding, and with a fixed uniform spacing of d = 10 nm, we see in Fig. 2(c) is largest on the cell axis ( r = 0) where the radial that pathway length and Rsh are maximum. Closer to the edge of the cleft, R,, decreases and the amplitude of is lower. The changes in peak with radial distance for the conditions of panel (c) are summarized in Fig. 2(d). Folding of the junctional membranes, while keeping d uniform throughout the cleft, increases both the length of the radial pathway for current flow and the membrane surface area. There is then more current available within the cleft, and the current is relatively more confined because of the increased R,,. For simplicity, we consider only a radially symmetric folding pattern with the amplitude of folding and the number of full cycles as parameters (see “Simulations”). Simulation results with a folding of five cycles [see Fig. l(b)] are shown in Fig. 3. Except for an increase in amplitude by a factor of about 3 due to the larger Rsh, there is no other difference due to folding as the thickness d is varied (compare (a) and (b) of Figs. 2 and 3). In other words, since i s h varies linearly withd or radial distance, changes in these parameters affect only the amplitude of q5e and not its general shape. Folding affects the local current density in the cleft and at various radial locations, as thus the amplitude of illustrated in Fig. 3(c) and (d). The actual pattern of five cycles (S = 2 pm) gives a folding peak on the cell axis ( r = 0) on the side of cell 2. There are also two other folding peaks on the side of cell 2 at r = 2 and r = 4 pm. Paradoxically, the amplitude of at the two latter folding peaks is appreciably larger than that at a trough located closer to the axis where R,, is larger. Thus the waveform at r = 2 is larger than that at r = 1, and the waveform at r = 4 is larger than that at r = 3 [Fig. 3(c)]. Variations
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 39, NO. 12. DECEMBER 1992
1236
4 1
1
0
(')
2
TIME (ms)
0.4
0.3
-i
0.2
E
0.1
.0.1
I 1
0
2
0
1
TIME (ms)
2
3
4
5
RADIAL DISTANCE (pm)
Fig. 2. Changes in q4
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. VOL. 39. NO. 12. DECEMBER 1992
1232
A Model Study of Electric Field Interactions Between Cardiac Myocytes HeWC Hogues, L . Joshua Leon, and Fernand A . Roberge, Senior Member, IEEE
Abstract-The transmission of excitation via electric field coupling was studied in a model comprising two myocytes abutted end-to-end and placed in an unbounded volume conductor. Each myocyte was modeled as a small cylinder of membrane (10 pm in diameter and 100 pm in length) capped at both ends. A Beeler-Reuter model modified for the Na’ current dynamics served to simulate the membrane ionic current. There was no resistive coupling between the myocytes and the intercellular junction consisted of closely apposed pre- and post-junctional membranes, separated by a uniform cleft distance. The membrane current crossing the prejunctional membrane during the action potential upstroke tends to flow out of the cleft, but it is partly prevented from doing so by the shunt resistance constituted by the cleft volume conductor. The prejunctional upstroke gives rise to a pulse of positive potential within the cleft which induces a small capacitive current across the post-junctional membrane to yield a small positive change in the intracellular potential in the post-junctional cell. The net result is an hyperpolarization of the post-junctional cleft membrane and a slight depolarization of the rest of the cell membrane since the extracellular potential outside of the cell is zero. The magnitude of this depolarization is quite small for a flat junctional membrane and it can be increased by membrane folding and interdigitation, so as to increase the junctional membrane area by a factor of 10 or more. Even then the post-junctional depolarization does not reach threshold when the extracellular potential around the post-junctional cell is effectively zero. Threshold depolarization occurs in the presence of a large decrease of post-junctional load, by increasing the junctional membrane capacitance andlor decreasing the volume of the post-junctional cell. Assuming that the normal resistive coupling between two cardiac myocytes is 1-4 MQ, our model study indicates that electric field coupling would then be about two orders of magnitude smaller. However, substantial enhancement of the efficacy of electric field transmission was observed in the case of cells with substantial junctional membrane folding.
I. INTRODUCTION HE spread of electrical activity in cardiac tissue involves a length constant which is much longer than a single myocyte, implying a sufficiently low resistance between neighboring cells to allow the flow of local circuit current. Anatomically cardiac tissue is made up of tightly packed rod-shaped cells appraximately 100 X 10 pm [ 11, [2]. The cells appear to be separated at their ends by in-
T
Manuscript received November 12, 1991; revised March 30, 1992. This was work supported by the Medical Research Council and The Natural Sciences and Engineering Research Council of Canada, and the Ministere de I’Enseignement Superieur et de la Science du Quebec. Dr. Leon was also supported by the Whittaker Foundation. The authors are with the Institut de Genie Biomedical, Ecole Polytechnique and Universite de Montreal, Montreal, P.Q. H3C 357, Canada. IEEE Log Number 9204046.
tercalated disks which define an intercellular gap which is on average somewhere between 8 and 20 nm [3]. These gap junctions contain tiny channels which permit the flow of ions between cells and provide the primary resistive electrical coupling [4].Studies of pairs of isolated ventricular myocytes have shown that the junctional resistance, Rj, linking the cytoplasm of the two cells is of the order of 1-4 M Q and that it is not changed by the transmembrane potential or the trans-junctional potential during the conduction of an action potential [5], [ 6 ] . Capacitive coupling between cardiac cells has also been postulated [7], based on observations that the junctional complex displays a time constant that is about 10% of the surface membrane time constant [8]. It has been proposed that the anatomical basis for the existence of a junctional capacitance could be the very close apposition of the preand post-junctional membranes in the area of the gap, combined with an extensive folding and a high degree of interdigitation of these membranes. The possibility that the connective tissue surrounding the cells may have capacitive properties has also been considered [7]. No junctional capacitive effects have been reported in studies of coupled isolated myocytes [6]. Previous model studies of the putative electric field interactions between cardiac myocytes have assumed that the pre- and post-junctional membranes were ordinary excitable membranes, and that the junctional resistance Rj was infinite [9].We maintain these assumptions here in order to examine the possibility of electrical field coupling in the absence of ionic local circuit current flow through R,. Instead of deriving an equivalent circuit to model the junction we use a model developed from electric field theoretic principles. Our model includes a threedimensional construct of the cell system, the conductivities of the intra- and extracellular milieu, and a HodgkinHuxley type representation of the membrane dynamics. Our purpose is to examine how the field coupling is affected by the physical characteristics of the intercellular gap, namely its size and the extent of the membrane folding, and the possible development of excitation in the post-junctional cell when an action potential reaches the junction. 11. METHODS Our study is based on a two-cell assembly lying in a large volume of conducting fluid. Each cell is modeled as a closed surface which separates a homogeneous internal
0018-9294/92$03,00 0 1992 IEEE
I
1233
HOGUES et a l . : FIELD INTERACTIONS BETWEEN CARDIAC MYOCYTES
current. The transmembrane potential is the difference between the intracellular potential 4j and extracellular potential + e , i.e., V = c$i The divergence of the current density is zero everywhere in the internal and external media. For the case of a single cell, we can write [lo]:
........._. >
n ’
,
JUNCTION
++
(a)
subject to boundary conditions
d .................. ~................. ~
% = -1 1, 87
(4)
0,
where 9 is the outward normal to the surface element dS. Green’s theorem [ 111 gives an expression for the potential at a point p on the inside or outside of the membrane surface S as follows: SIDE VIEW
(b)
END VIEW
4,(p)
=
--1
27r s a7 r
1 + 27r -
-+
s
@ ,rr ., d7 S
(6)
where 7 is the vector joining the source point to the field point, r is its length, and dS refers to the surface of integration S. This approach can be extended to consider n disjoint , Si,* , S”, and to include cells with surfaces S’, * a point source of strength Zstim in the intra or extracellular medium. We then have
-
--
(c) Fig. 1. Diagram of a two-cell system without resistive coupling. (a) Cell 1 of length LI and cell 2 of length L2 are aligned axially and have the same radius r . The junctional membranes are uniformly separated by a distance d. Stimulation is applied on the proximal cap of cell 1. (b) Pattern of junctional membrane folding, in the form of regular concentric peaks and troughs. The side view shows a longitudinal cross-section along the diameter of the cells where s is the radial distance between two adjacent peaks or troughs, and a is the amplitude of folding. We define the pattern in terms of the number of full cycles, which is 4 in this case with s = 2r/4. (c) Perspective view of the folded membrane at the junction.
medium of conductivity ai from an unbounded homogeneous external medium of conductivity a,. The two cells are abutted end-to-end to form a narrow cleft of volume conductor separating the prejunctional and post-junctional end-caps of excitable membrane (Fig. 1). We first derive a general formulation of the problem in relation to an arbitrary number of cells, using a boundary element method. Next we specify the parameters of the model and the simulation approach.
A. Boundary Element Method The transmembrane current i, (pA/cm2) is the sum of the ionic and capacitive currents Cm dV/dt + iion= i, (1) where Cm (pF/cm2) is the membrane capacitance, V(mV) the transmembrane potential, and iion(pA/cm2) the ionic
where 4: and 4; denote the potential distribution on the inside and outside surface of cell j . We discretize each surface Si into mi patches {S”, . . . Si” , Si,’}, and assume that &, 4,, and iion are constant, on any given patch S’”. Denoting these quantities by r#$’, 45, and i$’,,, then the system can be rewritten as 9
9
- -
(9)
I234
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 39. NO. 12. DECEMBER 1992
where R is the distance between the observation point and the source and a is the conductivity of the milieu. The average rise in potential caused by such a source on surface element A' is
JA, R
dS
(10) where d Q = V( 1 / r ) VdS is the solid angle subtended by dS at the center of patch S"'. Because the terms inside the integrals depend only on geometry they can be integrated to give a system of 2N equations for the potentials on either side of the N membrane patches, at each instant in time [ 101. The system can be written succinctly in matrix form. Ordering the potentials and currents in column vectors:
4 4.7
=
2na
1
(13)
dS
AI
Consider now a small resistive channel linking the intracellular media of two neighboring cells. The current crossing the channel is proportional to the potential difference between the inner surface of the two cells at the extremities of the channel. The current flowing across a channel of resistance R, joining cells k and i , at patches SkJ and Si'is
This current constitutes intracellular point sources of strength - I , on patch Sk*' and +Z, on patch Si'. The effect of an intracellular connection is included by introducing point sources of opposite signs at the two interconnection sites.
Equations (9) and (10) can be written as a system of first order non-linear ordinary differential equations of the form:
a+
M , - + M29 at
+ M31ion
= 0
(11)
where MI, M2,and M3 are matrices calculated from (9) and (10) (see Appendix A). Replacing the time derivative by a first order implicit finite difference with t, = ti-, At, the matrix system yields:
B. Simulations The above method allows to model several cells of virtually any geometry. For simplicity we represent each cell as a cylinder of membrane of radius r , capped at each end either by a disk of membrane or a pattern of folded membrane. The folding pattern chosen transforms the flat-disk gap junction into a set of two interlocking rippled surfaces of rotation. It is obtained by rotating the curve a 2
+
+(ti) = (MI
+ AtM2)-I(MI@(ti-l)
- AtM3Zio,).
(12)
Equation (12) allows an iterative solution for the potential distribution from an initial condition 9(0). Each element of the current vector Zion is obtained from the membrane at time model using the transmembrane potential (ai ti - I . The matrix (MI + AM,) is nonsingular and can be decomposed in LU form. Solution vectors are then obtained from successive forward and back substitutions, The resulting system of 2N coupled ordinary differential equations can be integrated in time to yield the potential on each side of the membrane. Point Sources and Intercellular Resistive Coupling: An important advantage of this approach is the possibility of locating current sources in either media, or on the boundary [ l o ] , [ 1 2 ] , [ 1 3 ] . The potential created by a point source of current of intensity I, is
4
4s = 2naR
I
+ a-2 cos nnp R
where a / 2 is the amplitude of the ripples, R the cleft radius, p the distance from the cell axis, and n the number of ripples. Hence the amplitude of the ripples is zero at the edge of the cleft and varies between + a / 2 and - a / 2 along the diameter. We define s as the length of a full cycle, that is the distance between two successive peaks or troughs along a diameter, Fig. 1 (b). The number of full cycles is also equal to n , so that s = 2 R / n . The surface area of the gap is S(a, n) =
sd
nRwJ4R2
+ a2n2n2sin (nnw) du
(16)
where R is the radius of the cell and w = p / R . Since it can be shown that S ( a / k , kn) = S ( a , n ) , there is flexibility in choosing the parameters a and n to obtain a given surface area. A perspective view of the folded membrane at the junction is shown in Fig. l(c). For example, with a = 1 pm and n = 13 (or alternatively with a = 0 . 5 pm and n = 26) the membrane area of one end-cap is in-
HOGUES et a l . : FIELD INTERACTIONS BETWEEN CARDIAC MYOCYTES
creased by a factor of about 9 in comparison with a flat disk of the same radius. Such a figure is in agreement with estimates of membrane folding at ventricular myocardial gap junctions [ 141. The cylinder surface was discretized into cylindrical patches of 2.5 pm in length. The flat disk end-cap with r = 5 pm was divided into 12 patches of roughly the same area by concentric circles. When the membrane area was increased by folding or by increasing the radius, the number of end-cap patches was increased proportionally. There is a constant separation of thickness d between the pre- and post-junctional membranes, Fig. 1 , irrespective of the shape of the end-caps. The cleft thickness in ventricular myocardium has been estimated to be on average between 8 and 20 nm [3]. The cleft might very well be narrower on the cell axis, however, where maximum confinement of the prejunctional membrane current would occur. For this reason the nominal value of d = 1 nm was chosen here to insure that maximum conditions of confinement were considered. The Beeler-Reuter model of the cardiac action potential [lS], as modified by Drouhard and Roberge [16], was used to describe the transmembrane current. An action potential was elicited in cell 1 by applying a current to the membrane of the proximal end-cap. It was carefully verified that the spatial and temporal discretization was fine enough to ensure that no numerical distortions occurred during propagation when nominal parameter values (Table I) were used. The time step used in all simulations was constant at 2 ps. Current and potential values were written to a file every SO ps. The method described above was programmed in Fortran on a Silicon Graphics 4D/25 workstation (Silicon Graphics, Mountain View, CA). All simulation were carried out without a resistive interconnection between the two cells, i.e., with Rj = 03, in order to emphasize the role of electric field interaction. Variables q5i, q5e, I,,,, and I,,, were monitored on all patches of the two junctional end-caps, and also on the cylindrical membrane at the periphery of the cleft. The mean values of potential and current in the gap (8e,j,, ,Im2) were computed from the value on each patch, multiplied by the number of patches, and divided by the total membrane area of the end-cap. 111. RESULTS As the action potential elicited on cell 1 reaches the junction, the transmembrane current crossing the prejunctional membrane generates changes in the extracellular potential within the cleft. There results post-junctional transmembrane potential changes which may or may not be sufficient to trigger an action potential in cell 2.
Extracellular Potential in the Cle@ on the cell axis, We consider first the cleft potential midway (d/2) between the two junctional membranes (Fig. 1). Current flows across the prejunctional membrane during the action potential upstroke and produces a smooth
1
1235
TABLE I NOMINAL PARAMETER VALUES Membrane capacitance, C,,,(pF/cm’) Cell radius, r (pm) Length of cell I , L , (pm) Length of cell 2 , L, ( p m ) Intracellular conductivity, U, (mu/cm) Extracellular conductivity, U , (mvlcm) Cleft thickness, d (nm) Amplitude of folding, a (pm)
1 .o 5 1000 IO0 5
66.66 1 1
monophasic waveform, as depicted in Fig. 2(a) for cleft thicknesses between 1 and SO nm. Since the current tends to flow radially out of the cleft, the cleft space can be represented by a shunt resistor (Rsh)whose value is directly proportional to the length of the radial pathway for current flow, and inversely proportional to the cleft thickness. With flat junctional membranes (no folding), the pathway length is equal to the radius of the cleft. Fig. 2(a) shows that a large Rsh, corresponding to d = 1 nm, allows a peak potential of nearly 4 mV. But q5e falls rapidly as d increases and becomes negligible at d = SO nm in this case. This behavior of peak q5c with increasing cleft thickness is summarized in Fig. 2(b). Still without junctional membrane folding, and with a fixed uniform spacing of d = 10 nm, we see in Fig. 2(c) is largest on the cell axis ( r = 0) where the radial that pathway length and Rsh are maximum. Closer to the edge of the cleft, R,, decreases and the amplitude of is lower. The changes in peak with radial distance for the conditions of panel (c) are summarized in Fig. 2(d). Folding of the junctional membranes, while keeping d uniform throughout the cleft, increases both the length of the radial pathway for current flow and the membrane surface area. There is then more current available within the cleft, and the current is relatively more confined because of the increased R,,. For simplicity, we consider only a radially symmetric folding pattern with the amplitude of folding and the number of full cycles as parameters (see “Simulations”). Simulation results with a folding of five cycles [see Fig. l(b)] are shown in Fig. 3. Except for an increase in amplitude by a factor of about 3 due to the larger Rsh, there is no other difference due to folding as the thickness d is varied (compare (a) and (b) of Figs. 2 and 3). In other words, since i s h varies linearly withd or radial distance, changes in these parameters affect only the amplitude of q5e and not its general shape. Folding affects the local current density in the cleft and at various radial locations, as thus the amplitude of illustrated in Fig. 3(c) and (d). The actual pattern of five cycles (S = 2 pm) gives a folding peak on the cell axis ( r = 0) on the side of cell 2. There are also two other folding peaks on the side of cell 2 at r = 2 and r = 4 pm. Paradoxically, the amplitude of at the two latter folding peaks is appreciably larger than that at a trough located closer to the axis where R,, is larger. Thus the waveform at r = 2 is larger than that at r = 1, and the waveform at r = 4 is larger than that at r = 3 [Fig. 3(c)]. Variations
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 39, NO. 12. DECEMBER 1992
1236
4 1
1
0
(')
2
TIME (ms)
0.4
0.3
-i
0.2
E
0.1
.0.1
I 1
0
2
0
1
TIME (ms)
2
3
4
5
RADIAL DISTANCE (pm)
Fig. 2. Changes in q4
