## Abstract

Electrical heterogeneities play a role in the initiation of cardiac arrhythmias. In certain pathological conditions such as ischaemia, current sinks can develop in the diseased cardiac tissue. In this study, we investigate the effects of changing the amount of heterogeneity and intercellular coupling on wavefront stability in a cardiac cell culture system and a mathematical model of excitable media. In both systems, we observe three types of behaviour: plane wave propagation without breakup, plane wave breakup into spiral waves and plane wave block. In the theoretical model, we observe a linear decrease in propagation velocity as the number of heterogeneities is increased, followed by a rapid, nonlinear decrease to zero. The linear decrease results from the heterogeneities acting independently on the wavefront. A general scaling argument that considers the degree of system heterogeneity and the properties of the excitable medium is used to derive a dimensionless parameter that describes the interaction of the wavefront with the heterogeneities.

## 1. Introduction

Each heart beat is associated with the propagation of an electrical wave through the cardiac tissue in a coordinated manner. The wave of activity moves through a heterogeneous tissue at both the macroscopic and microscopic scale. Nevertheless, in healthy tissue, the electrical conduction occurs as if through a homogeneous medium. Certain pathological states, such as ischaemia, fibrosis and cardiac sarcoidosis, can change the anatomical properties of the tissue and the dynamical properties of the cardiac impulse, heightening the impact of heterogeneities on wave propagation. Conditions in which the inception or propagation of the cardiac impulse is abnormal can lead to cardiac arrhythmias (de Bakker *et al*. 1988; Kawara *et al*. 2001; Hsia & Marchlinski 2002).

A common cause of cardiac arrhythmias are reentrant waves (Grant & Whalley 1998). A reentrant circuit involves a pathway that bifurcates into two branches. One pathway is blocked to anterograde conduction, but is excited in a retrograde fashion by the electrical impulse that travels through the unblocked path. Spiral waves have also been seen in cardiac tissue and preparations with no anatomical obstacle (Pertsov *et al*. 1993; Witkowski *et al*. 1998; Bub *et al*. 2002). The rotating waves of excitation have similar properties to the spiral waves seen in models of excitable media. In particular, wavebreaks in models of excitable media can generate spatiotemporal patterns similar to rotating waves of excitation associated with certain cardiac arrhythmias (Allessie *et al*. 1973; Winfree & Strogatz 1984; Pertsov *et al*. 1993). Accordingly, wavefront stability has been studied extensively in models of excitable media.

Wavebreak can occur both in homogeneous (Karma 1993) and heterogeneous (Bub & Shrier 2002) systems. Of particular interest to the present investigation is the effect of heterogeneities on the stability of propagating waves (Starobin *et al*. 1996; Bub *et al*. 1998, 2002; Xie *et al*. 2001; Bub & Shrier 2002). In a model of atrial fibrillation, wavelets of excitation propagate around areas of conduction block resulting from spatially fixed differences in refractory time (Moe 1964). Alternatively, anatomical obstacles with a space scale of comparable size relative to the wavefront may also result in a broken wave (Starobin *et al*. 1996; Bub & Shrier 2002). Such large obstacles are fortunately not seen in healthy tissue. However, in ischaemic tissue where blood flow has been compromised, cell death and/or a reduction in cell-to-cell connectivity can occur. As a result, wave stability is diminished by local asymmetries in cell coupling and by the heterogeneity associated with cellular death (de Bakker *et al*. 1988). In extreme cases, a considerable reduction in blood flow can cause a myocardial infarction where large areas of tissue death force electrical waves to navigate through tissue with considerable anatomical heterogeneity (de Bakker *et al*. 1988, 1993).

Previous studies from our group and others (Bardou *et al*. 1996; Bub & Shrier 2002; Panfilov 2002; Arutunyan *et al*. 2003; Bub *et al*. 2003; ten Tusscher & Panfilov 2003) have investigated the effects of varying connectivity, cell density and obstacle size on wavefront stability in experimental and theoretical systems. In the present work, we investigate the effects of changing the number of heterogeneities on planar wavefront stability in a cardiac cell culture system and a mathematical model of excitable media. In both the experimental and theoretical systems, we observe three types of behaviour as parameters are varied: plane wave propagation without breakup, plane wave breakup into spiral waves and propagation block. In the theoretical model, we observe a linear decrease in propagation velocity as the number of heterogeneities is increased, followed by a rapid, nonlinear decrease in propagation velocity to zero.

## 2. Experimental results

Heart cell monolayers are thin layers of tissue grown in culture dishes from embryonic or neonatal cardiac cells. Cardiac cells from very young animals have the capacity to easily form gap junctional connections with neighbouring cells in culture. After a few days in culture, embryonic cardiac cells are capable of supporting propagating waves of excitation over long distances. Cardiac monolayers allow controlled environments for studying conduction on microscopic and macroscopic scales.

In brief (Bub *et al*. 2003), cardiac monolayers are created by enzymatically isolating chick embryonic ventricular heart cells and plating them on an appropriate substrate for 72 h in maintenance medium 818a at 36 °C. In the present study, cells were plated on glass cover-slips treated with varying amounts of collagen (rat tail collagen Type 1, BD Biosciences, 8–1 μg cm^{−2}) and were loaded with a calcium sensitive dye to monitor activity. The calcium sensitive dye Calcium Green (5 μM, loaded for 30 min) was used to allow long recording times with high signal to noise ratios while minimizing phototoxic damage.

The activity was monitored using the macroscope-based optical mapping system discussed in Bub *et al*. (2003). The macroscope was used to perform low light level measurements at low-magnification scales (objective: Nikon 80 mm, imaging lens Chromicar Zoom 130, excitation filter 460 nm, dichroic beamsplitter at 510 nm, imaging filter 540 nm, Omega optical). Images were collected using a cooled charge-coupled device (CCD) camera (Princeton Instruments, Model TE CCD 576) at 20 frames per second. Adjacent pixels in the CCD were binned (2×2) and consecutive images were transferred to a computer for storage and analysis. Image data from each binned pixel was scaled based on its maximal range over 20 frames and the image was spatially averaged, background subtracted and viewed using Custom-Written software.

Figure 1 presents fluorescent images captured at 100 m intervals at each concentration of collagen. At high levels of collagen, monolayers display planar waves with no wavebreaks (figure 1*a*,*b*). For intermediate levels of collagen (figure 1*c*) we observe wavebreaks that result in reentrant activity. For low levels of collagen, monolayers display numerous holes and heterogeneities. Waves are fractured and form irregular shaped activation fronts (figure 1*d*). Activation fronts frequently break and form partial reentrant circuits, but fail to propagate throughout the whole medium.

The experimental data provides a challenge for theory. The collagen can affect the wave propagation in several different ways. Low collagen results in cells having a relatively higher affinity for cells than the cover-slip, which produces local variations in cell density. This may produce high-density regions that require more current to be excited than low-density regions, as well as low regions that are poorly connected to neighbouring cells. In addition, wavefront propagation may be affected by several other factors including local changes in ionic currents and indirect effects that might arise as a consequence of altered propagation velocity through the medium. In this study, we focus on the effect of high-density regions on the stability of the advancing wavefront by adding current sinks to a simple model of excitable media.

## 3. Mathematical model

The dynamics of wave propagation in inhomogeneous systems has been studied by analysing the movement of current between regions of current sinks and sources (Shaw & Rudy 1997; Bub & Shrier 2002). Any region of tissue that drains current from the wavefront can be defined as a current sink, while any region that supplies current to its neighbours is defined as a current source. As the propagating wave sources current to its neighbours, its neighbours become activated and the wavefront is propagated. The effect of a current sink on the wavefront depends in part on the connectivity between cells. In the case of low connectivity, the sink drains less current from the propagating wave causing the wavefront to be stabilized (Wang & Rudy 2000).

To simulate cardiac tissue, we make use of the nonlinear partial differential FitzHugh–Nagumo (FHN) equation. This prototypical model of excitable media, coupled to the diffusion equation, is given by(3.1)(3.2)where *D* is the diffusion coefficient and *I*(*r*) is an applied perturbation at the position *r*. The time course of the fast variable *v* in response to a supra-threshold stimulation has a similar form to the action potential of certain biological tissues. Owing to this property, in addition to the model's simplicity which allows for efficient computational and analytical explorations, many theoretical investigations of wavefront stability and propagation have used the FHN equations or its derivates as their model system (Glass & Josephson 1995; Aliev & Panfilov 1996; Rabinovitch *et al*. 1999; Berenfeld *et al*. 2001; Arutunyan *et al*. 2003). While there is no direct relation between variables of the FHN system and the biophysical properties of a cardiac cell, the variable *v* can be loosely associated with a cell's electrical potential and the slow variable *w* with the permeability of the cell's ion channels. The parameter values used in this investigation are *α*=0.02, *β*=0.5, *γ*=1.0, *δ*=0.0 and *ϵ*=0.01. Simulations were run using a diffusion coefficient of 0.0007 cm^{2} s^{−1} unless otherwise stated.

A two-dimensional sheet of cardiac tissue was simulated with a no-flux boundary condition. The computation used a forward time centred space finite-difference approximation, carried out on an *N*×*N* array, where *N*=200. The tissue was taken to be isotropic, with a constant space step of 0.02 cm; the time integration step was 0.1 ms. To elicit a planar wavefront to move across the lattice, a current (a perturbation) was applied along the left boundary with a magnitude of *I*=1.0 for 100 time units. Simulations were run until all cells had returned to their resting state following the initial perturbation or until the 30 000th time-step.

In this study, the effects of electrical heterogeneity on wave propagation were explored. Since the above represents an electrically homogeneous system, excitable cells were removed from the lattice and replaced with current sinks in order to model a form of electrical heterogeneity. A current sink was defined as a cell clamped at the globally stable steady state *v*=0 and *w*=0. Each current sink remained electrically coupled to its four neighbours. Through the coupling, the current sink could draw electrical potential from its active neighbours. The current sinks were randomly chosen across the lattice using a uniform probability distribution. Each node had an equally likely chance of being set as a current sink. Nodes were selected at random until the desired sink density on the lattice was obtained, with the density of current sinks varied across simulations.

### (a) Plane wave breakup in a heterogeneous medium

Computer simulations of the FHN system were first run with varying densities *ρ* of current sinks for *D*=0.0007 cm^{2} s^{−1}. Note that *ρ* is equal to the number of current sinks divided by the area of the system.

Three regimes were observed as the number of current sinks was increased. At low heterogeneity (approx. *ρ*<42.5 cm^{−2}) the plane wave propagates across the lattice without the plane wave breaking up into spiral formations or being blocked. While the current sinks produced small irregularities in the wavefront, the wave of excitation maintained an approximately planar geometry. For example, as the wavefront moved over a region containing current sinks, the excitable cells in close proximity to the sinks would not be excited because of the large amount of current being drained in that region. This would result in a transient break in the wavefront. If the group of current sinks was sufficiently small, the resulting wavebreak would gradually be filled by the diffusion of excitation between excitable nodes as the front continued towards the right boundary. The formation of small transient breaks in the wavefront occurred across the lattice wherever current sinks were located. Since the diffusion allowed the holes produced by the sinks to be filled as the wavefront moved forward, the stability of the wave was not compromised and plane wave breakup did not occur.

At a current sink density of *ρ*=42.5 cm^{−2}, plane wave breakup into spiral waves was first observed. For spiral formation to occur, a localized group of current sinks must cause a sufficiently large wavebreak. The excitation must then be able to propagate around the region containing the current sinks and back towards a previously excited region. Effectively, the propagation of the planar wave causes a directional asymmetry in the current sinks, which allows for a unidirectional conduction block. As described above, this is a defining characteristic of a reentrant circuit. A single heterogeneity is not sufficient to generate reentry in our model. Rather, several nearby heterogeneities are required to generate the substrate for reentry. Although we do not at the present time have a quantitative description of the geometry of the cluster of heterogeneities necessary to generate reentry, the effect of the distance between heterogeneities is discussed in §3*c*. Figure 2*a* shows the dynamics observed at 100 ms intervals of a simulation run with *ρ*=42.5 cm^{−2}. A spiral formed in the lower left corner of the lattice. The spiral formation persisted until the 30 000th time-step at which point the simulation was stopped. Due to the relatively high degree of heterogeneity in this system, the spiral arm contained several breaks (see second panel from the left) which in this case did not generate additional spirals.

Not all simulations run with current sink fractions of *ρ*>42.5 cm^{−2} resulted in spiral formation. For spiral patterns to be generated, the randomly distributed current sinks had to be arranged such that they produced a sufficiently large wavebreak about which a spiral could form. If the region near the initial plane wavebreak contained numerous current sinks, then the wave of excitation was unable to generate a spiral. In this situation, the broken plane wave would continue to propagate towards the right boundary with an approximately planar geometry. In other cases, current sinks would cause portions of the plane wave to break into fragments that would then propagate through the system. These fragments did not necessarily develop into stable spirals, but would either move towards the system's boundary or would be dissipated by the system's current sinks. The latter dynamic is depicted in figure 2*b* where a value of *ρ*=51.25 cm^{−2} was used.

Plane wave block was first observed with *ρ*=60 cm^{−2}. In plane wave block, the current sinks drain enough current to prevent the wavefront from propagating forward to excite upcoming cells. Figure 2*c* shows the dynamics observed with *ρ*=60 cm^{−2}. Extensive fragmentation and disruption of the planar wavefront can be seen. The irregularity in wavefront geometry becomes greater as *ρ* is increased.

### (b) Effects of heterogeneity on wavefront propagation velocity: the linear regime

The electrical heterogeneities clearly impact on the propagation of the wavefront. A decreased conduction is necessary to allow cells proximal to the unidirectional block to recover their excitability and ultimately be reexcited by the reentrant wave. Accordingly, any factor that decreases the plane wave's propagating velocity may increase the stability of a reentrant circuit and stabilize spiral waveforms. Similarly, the intercellular coupling impacts on the plane wave's propagation velocity. To help characterize their influence on wave propagation and stability, we analysed how the current sinks and intercellular coupling affect the wavefront's propagation velocity, *v*_{p}.

*v*_{p} was measured by recording the time course of the excitation variable *v* along the columns *x*=50 and *x*=200. This gave 200 pairs of equally spaced points perpendicular to the wavefront. The times at which the 400 cells were activated, defined as reaching a value of *v*=0.3 on the upstroke of the excitation profile, were recorded. The *v*_{p} for each pair of cells was computed as the distance between the cells divided by the time between the activation of the two cells. When the current sinks were positioned in a way that a recorded cell could not become activated, the *v*_{p} for the given pair was defined to be zero. The injection protocol described above causes a transient acceleration of the plane wave as it moves to the right. As a result, to compute *v*_{p}, *v* was recorded at cells with *x*=50 to allow for the wave velocity to stabilize.

Figure 3 shows the measured mean *v*_{p} of simulations run at a given *ρ* for *D*=0.0005, 0.0007, 0.0010 and 0.0012 cm^{2} s^{−1}. For low levels of heterogeneity, there is an approximately linear decrease in *v*_{p} for increasing *ρ*. As *ρ* is further increased, the data deviate from the linear trend, ultimately resulting in an abrupt transition to wave block. The increased heterogeneity disrupts the organization of the plane wave. For example, in the case of *D*=0.0007 cm^{2} s^{−1}, the deviation becomes more apparent around *ρ*=42.5 cm^{−2}, corresponding to where spiral formations were first observed, as described above. Interestingly, the transition occurs at lower values of *ρ* as *D* is increased, suggesting decreased plane wave stability at increased *D*. This result is further discussed below.

The linear decrease in *v*_{p} for increasing *ρ* suggests that for low numbers of current sinks, each sink is acting independently on the wavefront. To test this hypothesis we proceeded as follows. We constructed a 200×200 lattice with no-flux boundaries on the left and right borders and periodic boundary conditions for the top and bottom. The periodic boundary conditions connect the bottom nodes (*y*=1) with the top nodes (*y*=200) by diffusion, making the lattice effectively a cylinder with its axis parallel to the *x*-axis. As a result, the *y*-coordinate becomes arbitrary. With the above periodic boundary conditions a single current sink was placed on the lattice at *x*=1. Using *D*=0.0007 cm^{2} s^{−1}, a plane wave was initiated using the same protocol as above and the average propagation velocity between columns *x*=50 and *x*=200 was computed. The horizontal position of the current sink was incremented by one, the simulation run and *v*_{p} computed as before. This procedure was continued until the end of the lattice was reached. Combining the results of the simulations gives the average *v*_{p} that results from a single sink located at a given *x*-coordinate. The computed profile of the decrement in *v*_{p} (Δ*v*_{p}) relative to the *v*_{p} in a homogeneous system (*ρ*=0) is shown in figure 4.

If we assume that the current sinks are uniformly distributed across the lattice and that each sink acts independently on the wavefront, we can calculate the expected propagation velocity as(3.3)where is *v*_{p} in the homogeneous system (*ρ*=0), *ρ* is the current sinks density, *f*(*n*) is the computed profile given in figure 4, *X* and *Y* are the number of horizontal and vertical lattice nodes, respectively, and 0.02^{2} is the area of an individual lattice node. Using this equation we obtain as a function of *ρ*. The predicted propagation velocity along with the measured propagation velocity is shown in figure 5. There is an extremely close fit of the theoretical expression to the data over a limited range of densities of sinks. Thus, we believe that for low sink densities, the assumptions underlying the derivation of equation (3.3) are valid and that each sink is exerting a slowing of the velocity independent of the other sinks. However, as the sink density increases, the velocity decreases more rapidly. This regime is further discussed below.

A possible interpretation of our use of current sinks is that the current sinks act locally to decrease the excitability of the surrounding region. This suggests that the current sinks may be equivalent to uniformly decreasing the excitability in all cells. To address this possibility, simulations were carried out on a homogeneous system (*ρ*=0) while applying a uniform, global sink to all nodes in the lattice. This was done by setting *I*(*r*) of equation (3.1) to a non-zero value at each lattice position for the duration of the simulation. As the size of global sink was increased, a linear decrease in propagation velocity was observed. At *I*(*r*)≈−0.0012, the plane wave could no longer propagate and wave block occurred. As expected, spiral wave formation was not observed for any value of the global sink. The linear decrease in propagation velocity produced by increasing the magnitude of the global sink is shown in figure 5.

### (c) A dimensionless number characterizing wave breakup in heterogeneous excitable media

In order to gain insight into the underlying mechanism of the breakup, we identify a dimensionless number that reflects both the heterogeneity as well as the properties of the excitable medium. In the present simulations, a single isolated heterogeneity generates a transient wavebreak that seals once the wave propagates some distance from the heterogeneity. If we assume that wavebreaks only persist due to the interaction of the wave with more than one heterogeneity, we can define two relevant time-scales. One relates to the time it takes a wave to propagate the average distance between heterogeneities. In excitable media, the velocity of propagation is proportional to where *D* is the diffusion coefficient and *τ* is the rise time of the excited phase (see p. 236 in Winfree 1993). Consequently, the transit time of the excitation between heterogeneities is proportional to , where *ρ* is the current sink density as defined above. The other relevant time-scale is the refractory period *R*, which is approximately equal to the duration of the active phase of the excitation. Dividing the two time constants, we obtain the dimensionless parameter(3.4)

We believe that *σ* captures the features governing the interaction of the propagating wave with the heterogeneities, and that consequently, there will be similar qualitative properties in the dynamics for the same values of *σ*. To support this, keeping the kinetics of the FHN model constant gives *R*≈90 ms and *τ*≈24 ms, but varying *D* and *ρ* we find from figure 3 that the value of *σ* where breakup occurs is *σ*≈0.10, 0.10, 0.11 and 0.11 for *D*=0.0005, 0.0007, 0.0010 and 0.0012 cm^{2} s^{−1}, respectively.

Thus, as diffusion and heterogeneity vary, the breakup of the wave appears to occur at approximately constant values of *σ*. Although our current interest is in wave propagation in heterogeneous cardiac tissue, we believe that the current characterization of the instability should have implications for wave propagation in other heterogeneous excitable media. However, the constant value of *σ* in the present study assumes that more than one heterogeneity is required to generate reentry, that different heterogeneities have equivalent effects on wavefront stability and that the diffusion coefficient and refractory period do not modify wavefront stability in the absence of heterogeneities. Thus, further studies are required to test this relationship in other types of excitable media and with other types of heterogeneities.

## 4. Discussion

In this study, we investigated two systems with randomly distributed heterogeneities: an experimental cardiac monolayer system, where heterogeneities were added by altering the adhesive substrate and a theoretical model where heterogeneities were modelled by randomly distributed sink cells. In both the experimental and model systems, we observed transitions from stable waves to spiral waves to block as the degree of heterogeneity was increased.

In the theoretical model, we examined the impact of the current sinks on *v*_{p} and intercellular coupling on propagation. At low values of *ρ*, we observe a linear decrease in *v*_{p} as *ρ* was increased. In this linear regime, the current sinks act independently on the wavefront to decrease *v*_{p}. By summing the individual effects of current sinks at a given *x*-coordinate, we were able to predict the decrease in *v*_{p} using equation (3.3), suggesting that the current sinks act independent of each other at low *ρ*. Our prediction slightly overestimated the effect of the current sinks on *v*_{p}. The overestimation may be a result of using periodic boundary conditions to construct the profiles given in figure 4 compared to the no-flux boundaries of the initial simulations.

Following the linear regime, there is an abrupt transition towards wave block. This transition was observed in simulations for all values of *D* tested. Interestingly, this transition occurred at different sink densities, with wavebreak occurring more readily at the higher *D*. We account for the paradoxical decrease of wavefront stability by deriving a dimensionless scaling parameter *σ* (equation (3.4)). Since the variables used to derive equation (3.4) are common to all excitable media with randomly distributed heterogeneities, we expect to see this scaling relationship in a wide variety of excitable systems. Paradoxical increases in spiral wave stability were reported by Panfilov (2002) in an excitable media with local variations in connectivity and ten Tusscher & Panfilov (2003) in excitable media with randomly distributed heterogeneities. In both cases, increased stability was attributed to increased spiral period. Decreased connectivity has also stabilized waves in the presence of spatial heterogeneity (Rohr *et al*. 1997) and decreased excitability (Shaw & Rudy 1997); however, in these cases wavefront stabilization was accounted for by an analysis of current sources and sinks.

The present results on wave speed as a function of sink density are qualitatively similar to simulation results obtained on a one-dimensional strand of cardiac cells where excitability was varied (Shaw & Rudy 1997). The addition of current sinks in our model may be interpreted as lowering the local excitability of the lattice, which suggests that the results shown here may be equivalent to decreasing excitability equally in all cells. Simulations where a global sink current was applied to all cells in the simulation display a similar reduction in propagation velocity and eventual block without spiral wave formation.

In other experimental and theoretical systems, the scaling relationship defined in equation (3.4) does not appear to hold. Notably, the cardiac monolayer preparation displays increased wavebreaks as connectivity is decreased (Bub *et al*. 2002). A possible explanation for this discrepancy is that a monolayer with low-connectivity and a large number of heterogeneities is more appropriately modelled by a discrete system. Simulations with discrete cellular automata models display similar behaviour to the results found in monolayer experiments (Bub *et al*. 1998, 2002; Bub & Shrier 2002).

The present study has several limitations. First, we simulate local heterogeneity by adding current sinks, which greatly oversimplifies the effects of varying the collagen substrate on propagation. Alternative approaches, such as modelling heterogeneities as holes, or modelling the effects of local changes in connectivity, should be assessed in the context of the present work. Second, the FHN model does not directly simulate ionic currents that are known to affect wave speed and wave stability. The simulation results in this paper should be confirmed in ionic models of cardiac conduction.

In conclusion, the effect of randomly distributed heterogeneities on wave stability was characterized as a function intercellular coupling and heterogeneity density. We propose a new dimensionless parameter to characterize the interaction of the propagating wave with the heterogeneities. While this work made use of a simple model of cardiac tissue, the mechanism should be applicable to the study of wave propagation in other heterogeneous excitable media.

## Editors' note

Please see also related communications in this focussed issue by Fink *et al*. (2006) and Marée *et al*. (2006).

## Acknowledgments

The authors thank Katrin Rohlf for helpful comments. This research has been supported by funding from the Canadian Heart and Stroke Foundation, CIHR, MITACS and the NIH National Center for Research Resources.

## Footnotes

One contribution of 13 to a Theme Issue ‘Biomathematical modelling I’.

- © 2006 The Royal Society