Contents

Chapter 1
Introduction

The interaction of pacemakers is a broad area of study with applications in a variety of disciplines[97,241,265]. Coupled oscillators are investigated in pure and applied mathematics[147,188]. Autocatalytic reaction diffusion systems are studied in chemistry[84]. Dynamics of oscillating circuits are of interest to electrical engineers[73]. Biologists study entrainment in firefly populations[265], physiologists study coupled cardiac cells[65,255], and neurologists study the interaction of spontaneously firing neurons[124]. Pacemakers are of primary importance in the field of cardiology; hearts can be considered as a group of connected cells beating at a rhythm determined by few spontaneous pacemakers. Arrhythmias can result when normal pacemaker interactions break down. Nonentrained pacemakers give rise to relatively benign arrhythmias such as parasystole[180], or potentially life threatening conditions such as multifocal atrial tachycardia[161].

The dynamics of pacemakers under conditions of periodic forcing, and as globally and locally coupled populations, has been the subject of extensive theoretical investigation. However, pacemakers do not interact directly in cardiac preparations, but through a medium with its own dynamic properties. The present work investigates two cases in which pacemakers interact with a conductive excitable medium.

We study cardiac pacemaker/excitable media interactions in model systems using nonlinear mathematics and optical mapping techniques. This introductory chapter is a background on theoretical and experimental work on pacemakers and excitable media. The dynamics of two pacemakers coupled via a simplified excitable media is discussed in the context of a theoretical model for parasystole in Chapter 2. The assumptions introduced by the addition of an excitable medium dramatically alter the behaviour of the two oscillator system. Experimental investigation of coupled pacemaking systems and excitable media in cardiac preparations requires observing the behaviour of the system at many points in space simultaneously. Optical techniques for mapping activity in cardiac preparations are introduced in this chapter and discussed in detail in Chapter 3. Preliminary results and pilot studies using optical mapping techniques are presented in Chapter 4. The dynamics of a large number of interconnected oscillators is investigated in an experimental and theoretical model in Chapter 5. Here, the dynamics can be primarily driven not by the pacemakers themselves but by the properties of the underlying excitable media.

1.1  Nonlinear Dynamics.

Mathematical models can be formulated as differential equations or finite difference equations. Differential equations are equations of the form [dx/dt]=f(x), where f(x) defines the rate of change of variable x over time t. The above equation is one dimensional in the sense that x depends only on itself. Mathematical models of biological systems are often multidimensional, with several interacting variables varying continuously over time. Finite difference equations have the form x_t+1 = f(x_t). Here the value of variable x is determined some constant time in the future (t+1) by function f and the state of the variable at time t. A finite difference equation is iterated over time by repeatedly substituting x_t+1 for x_t. Despite their apparent simplicity, one dimensional difference equations (like the equation above) are capable of displaying a wide range of qualitatively different dynamics. For clarity, the remainder of the discussion in this section considers one dimensional difference equations, however the concepts can be applied to the study of differential equations as well.

An equation is considered linear if the right side of the equation describes a line a x_t + b, where a and b are constants. One dimensional linear equations have a limited range of behaviours. They either reach a steady state, a cycle or approach positive or negative infinity. A steady state is a solution where x_t+1=x_t A period P cycle occurs when x_t+P=x_t, where x_t+N ¹ x_t for N < P. Finite difference equations are non-linear if the right side of the equation describes a form other than a line. Non-linear finite difference equations can display chaos and quasiperiodicity in addition to the behaviours described above. Both chaos and quasiperiodicity are aperiodic (x_t+P ¹ x_t for all P.) In quasiperiodicity, two initial conditions initially close together stay close as time evolve. An example of quasiperiodic behaviour is two uncoupled clocks, running at slightly different speeds. The difference in phases (here the time read by the clock) change slowly over time, and a small change of phase of one clock does not change the overall behaviour of the system. Chaotic dynamics, on the other hand, are sensitive to small changes in initial conditions. In chaotic dynamics, two initial conditions that are close together will diverge after several iterations. Chaos is usually characterized by dynamics that appear random at first glance.

The trajectory of an equation refers to the sequence of values obtained by repeatedly iterating f(x_t). Nonlinear finite difference equations undergo large change in dynamics when their parameters are slowly varied under certain conditions. Such a change is referred to as a bifurcation. An example of a bifurcation is when a trajectory moves from a steady state solution to a cycle. Chapter 2 investigates the bifurcation structure of a finite difference equation formulated as a circle map. Circle maps are addressed later in this chapter and are discussed in detail in Chapter 2.

1.2  Pacemakers

Pacemaking activity is seen in many biological systems. Cells with pacemaking activity produce periodic signals. Some well known examples are the circadian rhythm, which controls the 24 hour cycle in animals[162], islet cells in the pancreas, which give periodic bursts of activity releasing insulin[49], and cardiac sinus[255] and atrioventricular (AV) node cells, which control the periodic rhythm of the heart[97].

Mechanisms of pacemaking activity has been extensively studied. Theoretical studies involving ionic mechanisms have been used to give detailed analysis of sinus node[65], AV node, b-pancreas cells[28] and other systems[196]. These models use differential equations to describe a slow current that causes membrane voltage to rise to threshold initiating an action potential. Such models have been used to study the effects of external stimulation on single or small groups of cells, and have been used to study the interaction of large numbers pacemaking cells[173]. As these models are complex, simple models of pacemakers are used as tools to develop a theoretical basis for pacemaker behaviour. These simple models use equations that do not represent ionic currents directly, but rather the topological characteristics of oscillators.

Any system that exhibits a periodic behaviour is called an oscillator. Mathematical models of oscillators are typically formulated as differential equations, in which the variables represent the key components needed to describe the oscillator. For example, for the ideal pendulum, the phase space consists of position and velocity and the trajectory is approximately elliptical for small amplitude oscillation.

Biological oscillators differ from the pendulum in that they tend to have a characteristic period and amplitude, which they return to following perturbations. Oscillators that return to a set trajectory are called limit cycle oscillators. A simple limit cycle oscillator, called the Poincaré oscillator can be given by:

dq/dt = 2p   mod (2p) dr/dt = Kr(1-r) (1.1)

where q gives the rotation angle and r gives its distance from the origin. K is a parameter that regulates the rate at which the value of r approaches 1. Equation (1.1) has an unstable steady state at r=0 and a stable limit cycle at r=1. Any perturbation results in a spiraling of the trajectory back to r=1. Despite its simplicity, this model has been used to characterize the behaviour paced cardiac oscillators[109].

1.2.1  Patterns of Behaviour of Coupled Oscillators.

1.2.2  Circle Maps.

Circle maps are reviewed in detail in Chapter 2. Briefly, circle maps can be characterized by their topological degree. The topological degree of F(q) counts the number of times F traverses the unit circle as q traverses the circle once. For example, F(q)=Mq+b, mod(1) has topological degree M. The dynamics of circle maps can be quantified by their rotation number, which gives the average change in F(q) over many iterations. For example, 2:1 phase locking results in a rotation number of 1/2. Rational rotation numbers are associated with periodic dynamics, while irrational rotation numbers are hallmarks of quasiperiodicity and chaos. The class of circle maps investigated in Chapter 2, however, display chaotic dynamics while having rational rotation numbers. This results in surprisingly regular chaotic trajectory called `banded chaos'. This unusual behaviour is a result of fractional topological degree caused by a jump discontinuity in the circle map.

1.2.3  Multiple Pacemaker Interaction.

The dynamics become more complex when more than one pacemaker is involved. Each oscillator can be coupled to its nearest neighbours or to all the oscillators in a large population. Global coupling interactions are observed in firefly populations[241], and synchronization of women's menstrual cycles[167]. Nearest neighbour interactions in multiple oscillator systems are studied in biological contexts with neuromuscular oscillators in the intestine[45], pacemaker neurons in the eye[124], spinal cord[33,256], and in culture[183].

Computer simulations and mathematical analysis demonstrated that globally coupled pacemakers necessarily synchronize or their phases remain evenly scattered if the oscillators are have identical intrinsic cycle lengths and the coupling is weak[241]. Under some coupling rules, the population of oscillators breaks down into subpopulations of synchronous oscillators called `clusters'. The number of clusters depends on the nature of coupling between oscillators[188]. Winfree[262] showed that coupled oscillators with different frequencies can still synchronize if their intrinsic frequencies are not far from a global mean value. Synchrony can persists even when some stochasticity is included in the form of random external forcing[188].

Locally coupled systems add a level of complexity by allowing the formation of spatial structures. Clustering still exists in locally coupled systems with identical frequencies [189], but the stability of the clusters depend on the conduction delay between oscillators. Simulations and analysis with one dimensional chains[147,256,33] of weakly coupled oscillators with different intrinsic frequencies have been conducted in the context of studying conduction in the lamprey spinal cord[33]. Frequency differences cause phase lags in conduction, and can cause a loss of entrainment if the differences in intrinsic frequencies are too great. Phase lags can be lengthened or shortened depending on the distribution of coupling strengths in the chain[33].

The above analyses on globally and locally coupled oscillators assume that coupling is extremely weak. This assumption ensures that oscillator timing maintain a relative order: an oscillator with phase lagging behind other oscillators can never overtake them. Under conditions of strong coupling, oscillators can change their order, resulting in mixing of their phases. Strongly coupled oscillator systems can display changes in the collective frequency and exhibit chaos[188].

1.3  Excitable Media.

Travelling waves in excitable media have been observed in many contexts. Travelling regenerative waves are observed in chemical systems, the most well studied being the Belousov-Zhabotinskii (BZ) reaction[74]. Other examples of waves in excitable media include autocatalysis reactions on metal surfaces[89], and propagation of forest fires[9]. In biological contexts propagating waves are used for communication within and across nerves, to generate contraction in a ordered fashion in the heart, and organize slime mold[201],in response to starvation. Despite the obvious differences in the above systems, propagation in all excitable media share many characteristics. The following briefly discusses features common to all excitable media systems.

A excitable media system can be considered as a group of individual elements coupled to each other. Each element can pass information to its neighbors. In physical system signals pass by diffusion; in the case of the BZ reaction chemical species (HBrO₂ and Fe³⁺) diffuse to neighbouring regions, in cardiac muscle current carried by sodium and potassium ions moves to neighboring cells. Each (coupled) point in the excitable media is characterized by a rest state that is it stable for small perturbations. A impulse with strength greater than a certain threshold can cause that point to undergo a large excursion from, and eventually return to, its rest state value. The length of time required to return close to the steady state value is a factor that determines the refractory time of that point; a refractory point cannot easily undergo another cycle until it recovers.

An impulse over a certain threshold initiates a wave of activity moving across the excitable media. As each element undergoes an excursion from steady state, it causes its neighbors to move over threshold at a rate determined by the diffusion coefficient (a `passive' property of the media), and the rate of rise of the diffused species of the excited element (a `active' property of the media.) A wave of activity propagates with a speed controlled by how fast elements ahead of the wave are induced to cross the threshold. For planar waves, wave speed is only a function of the active and passive properties of the media. However, wavefront curvature can affect its speed. The reason for this is that diffusion to points ahead of a curved wavefront distribute over an area that depends on the wavefront shape. For example, convex (curving outward) waves must drive a greater area of excitable media over threshold than concave waves do. Wave speed (c) for excitable media can be approximated by the following relation:

c = c0 + Dr,

(1.2)

where the curvature r is defined as the negative reciprocal of the local radius of curvature, D is the diffusion coefficient, and c₀ is the planar wave speed. Equation (1.2) is an approximation valid for r close to zero. Derivation of an exact relation depends on the system being studied, however Equation (1.2) has been verified experimentally for many physical and theoretical systems[84,85]. Equation (1.2) is a good first approximation over a large range of curvatures in these systems.

Since conduction velocity decreases as wavefront curvature becomes more negative, Equation (1.2) shows that velocity drops to zero for a critical radius of curvature (r_c):

rc = - c0 / D.

(1.3)

assuming that Equation (1.2) is valid at high r. It also follows that a circular nucleus of excitation does not result in outward propagation if r £ r_c.

A point seldom addressed in the study of excitable media is the effect of curvature on wavefront stability in heterogeneous media. Heterogeneities are more likely to prevent propagation at highly curved fronts than at wave fronts with shallow curves. Consider a case where small areas ahead of a target wave are refractory. The target wave with refractory holes is generated as the excited wave moves around refractory regions. This resulting patchy wave is less able to drive media ahead of it over threshold than a fully active wave, so a wave may not propagate even if r ³ r_c. Decreasing r lowers the relative contribution of refractory regions, allowing propagation. This mechanism is investigated in more detail in Chapter 7.

Another feature common to excitable media is dispersion effects[90]. Dispersion refers to the dependence of wave speed on recovery time of the medium. During periodic signaling at high rates, the medium does not have enough time to fully return to its rest state between waves. Consequently, the threshold increases as more excitation is required to trigger the next wave. If D is constant, it takes longer for a region ahead of the wavefront to reach threshold, resulting in slower wave speed. The dependence on conduction velocity c on recovery time T_rec can usually be fitted by an exponential function for diverse systems[90,123]:

c=1/aexp(-Trec/ b) + g,     for Trec ³    refractory time

(1.4)

where a, b and g are positive constants. Equation (1.4) is known as the dispersion relation.

1.3.1  Two Pacemaker Interaction via an Excitable Media: Parasystole.

In the healthy heart, the cardiac rhythm is set by a single pacemaking site called the sinoatrial node (see Section 1.5.1 for background information on cardiac conduction.) Occasionally, there can be pacemakers in abnormal locations that generate rhythms that compete with the normal sinus rhythm. The interaction between the ectopic pacemaking sites and the sinus pacemaker generate an arrhythmia called parasystole.

Parasystolic rhythms can be pure or modulated. In pure parasystole, the sinus and ectopic pacemakers do not phase reset each other. In modulated parasystole, the sinus rhythm acts to phase reset the ectopic focus. This introductory section discusses the dynamics of pure parasystole. Chapter 2 gives a detailed analysis of modulated parasystole.

The ectopic and sinus pacemakers interact via excitable cardiac tissue between them. An ectopic beat results if the ectopic focus depolarizes outside the the refractory period of the surrounding tissue, otherwise it is blocked. The sinus beat following an expressed ectopic beat is usually blocked due to refractoriness of the tissue surrounding the sinus node. This block is called a compensatory pause. In pure parasystole, the time intervals between ectopic beats are a multiple of a common divisor, and the time interval between sinus and ectopic beats is variable.

The above observations on the behaviour of parasystole allow the excitable media between the pacemakers to be modeled as a few simple rules. More complex models of excitable media are considered in the next section. A model for pure parasystole can be stated as a circle map[39,96] (see Section 1.2.2.) The sinus cycle is represented by a circle and the hearts refractory period is a segment of this circle. The ectopic pacemaker produces impulses that land within the sinus cycles phase. The ectopic beat is blocked if it occurs within the refractory period, and is expressed otherwise. Mathematically, this is formulated as a map of a circle into itself, f_t+1=f(f_t) mod 1, where f is the phase of the sinus cycle when the ectopic beat occurs. Glass and colleagues[96] looked at the number of sinus beats between ectopic beats obtained by iterating the circle map. Iterating the map results in a sequence of values that represent the number of interpolated sinus beats (NIBs) between each ectopic beat. Glass and colleagues[96] proved that the following four rules govern the NIB sequence: 1) there are three or less NIB values in a sequence of NIBs, 2) if here are three NIB values, the sum of the two smaller values is one less than the larger one, 3) only one of the values is odd, 4) only one of the allowed values can succeed itself in the sequence. The above rules hold only when the mathematical model of parasystole can be represented by a 1:1 invertible circle map. If the sinus node resets the ectopic focus, this condition no longer holds. Chapter 2 investigates the dynamics of discontinuous noninvertible circle maps, and applies these results to a mathematical model of modulated parasystole.

1.3.2  Mathematical Models of Excitable Media.

Excitable media are most naturally represented as partial differential equations. Here the excitable media is divided up into individual cells, where each cell is modeled by differential equations. The number of differential equations required to represent a cell in the system may be large. For example, the DiFrancesco-Noble model of Purkinje fibers[59] has 14 dimensions and over a hundred parameters. The rational for developing such complicated theoretical models is the assumption that inclusion of all known components is needed to achieve an accurate simulation of the behaviour. Although this modeling approach may be justified in some situations (where the effects of altering a specific variable are to be tested, for example), the high dimensionality of such models hinders intuitive understanding of the underlying dynamics.

A generic excitable media model can be represented in a simplified form by the interaction of two variables: a excitation variable (u) and recovery variable (v). The variables interact locally according to the ordinary differential equations du/dt = f(u,v), dv/dt = g(u,v), where the nullclines for functions f and g have a characteristic shape shown in Figure 1.1. In this model, there is a rest state where the nullclines intersect. A small perturbation is damped out, but a perturbation over a certain threshold triggers a long excursion as shown by the dotted line (Figure 1.1).

The excitation variable u increases rapidly until the trajectory approaches the rightmost arm of f(u,v)=0. At this point the trajectory moves slowly up the nullcline as the recovery variable increases. When the trajectory reaches the top of the left arm of f(u,v)=0, the trajectory moves to the leftmost arm as variable u rapidly decreases, which is followed by a slow decrease of v back to the rest state.

A simple mathematical model that behaves as described above is the FitzHugh-Nagumo[76] model of excitation in nerve and muscle tissue:

du/dt = f(u,v)=u-u3/3 -v dv/dt = g(u,v)=E (u+a-bg) (1.5)

where u is membrane potential and v approximates a slow current (when E << 1.) The excursion that u takes in response to super-threshold stimulus is similar in appearance to voltage changes called action potentials (APs) produced in nerve and muscle tissue. Although the model is unable to reproduce the detailed characteristics of APs in nerve and muscle fibers, it successfully encompasses many characteristics of excitable tissue including threshold behaviour, a refractory period, and dispersion effects.

Modeling wave propagation in two dimensions with a model similar to Equation (1.5) involves having several equations coupled to each other by diffusion:

du/dt = f(u,v)+ Dx d2(u) d x2 + Dy d2(u) d y2 , dv/dt = g(u,v)+ Dx d2(v) d x2 + Dy d2(v) d y2 (1.6)

where D_x and D_y are diffusion constants in the x and y directions. Equation (1.6) approximates continuous diffusion if the integration time step is small. More complex coupling relations can simulate preferential propagation in one axis[216,270] or take into account different diffusion rates for substances in different compartments (for example, the bidomain model of cardiac conduction uses different conductances for substances inside and outside cells[270].)

Much effort has been expended in creating different coupled differential equation models for different biological and physical systems. For example, the Hodgkin-Huxley model[114] and FitzHugh-Nagumo equations[77,187] model nerve tissue specifically, the Beeler-Reuter model[12] models heart tissue, the Oregonator model[74] models the BZ reaction, and the Martiel-Goldbeter equations[168] model slime mold aggregation. All differential equation models represent known microscopic quantities (such as the movement of a particular current) with equations. The different underlying mechanisms results in models with few obvious mathematical similarities.

Accurate simulation of travelling waves becomes computationally expensive as complexity of the underlying models increase. Investigating the dynamics of travelling waves where the wavefront is involved in subtle changes requires integration of `stiff' (small time step) equations which can result in thousands of iterations for each wave of excitation. If a large number of spatially discrete points is required, simulations may require a prohibitive amount of computer time[90]. In addition, changing one variable can change more than one characteristic of the model. For example, lowering the refractory period in Equation (1.5) decreases upstroke velocity as well (if one decreases the slope of g(u,v).) This problem is exasperated with higher dimensional models, and it becomes increasingly difficult to ask simple functional questions. Simulations in the Beeler-Reuter[12] model could easily determine what changes a 10 percent increase in extracellular sodium would cause, but cannot answer how a 10 percent increase in refractory period alone would change its dynamics. Complex models are therefore not appropriate for generating general results applicable to all excitable media.

Another approach for modeling excitable media that has recently been explored is the use of diffusively coupled difference equations. The coupled map lattice is discrete time (as each point in the lattice is a difference equation), continuous state model[11,131,235]. Each point is diffusively coupled to its neighbours and the lattice is updated using small time steps, making it continuous in space in the same sense a coupled differential equation model is. A coupled map lattice may be appropriate when each point in the physical system is well defined by map, and spatial averaging of the difference equation variable can be used to model diffusion.

1.3.2a   Cellular Automata.

Cellular automata models can be labelled as first generation or second generation for the purposes of this discussion. First generation cellular automata are characterized by cell-cell interaction limited to adjacent neighbours, a small number of update states, and cell spacing on a regular grid. First generation cellular automata have several shortcomings that are addressed by more complex second generation cellular automata models. First generation and second generation cellular automata are discussed below.

First generation cellular automata models usually have elements existing in one of three states, resting, excited or refractory. A resting cell updates its state based on the activity of its neighbours, while excited and refractory cells update their state based on the cells history. A resting cell remains at rest until a certain number of its neighbours are excited, in which case it becomes excited in the next time step. Excited cells become refractory and refractory cells return to rest.

A given cell has a number of neighbours which depend on the radius of interaction and the shape of the neighbourhood. Most first generation cellular automata use either a Neuman or Moore neighbourhood geometry with radius equal to one. A Neuman neighbourhood consists of the 4 cells on each of the cells edges. A Moore neighbourhood of radius N includes all cells encompassed by a square of side length 2N+1.

First generation cellular automata models of excitable media have two advantages over differential equation models; they are more intuitively transparent and computer simulations run far faster. Their simplicity come at a price, however. Radius of curvature and dispersion effects, which are common to all excitable media in nature, are absent in these models. Another problem encountered in first generation cellular automata models is that wave propagation direction is influenced by the underlying lattice. In contrast to the homogeneous system modeled by differential equations, waves in cellular automata models often have sharp edges not found in nature. The shape of the grid results in preferred conduction in the horizontal and vertical directions, while diagonal speed is reduced. This can result in pronounced artifact that can effect dynamics as parameters are changed. For example, Weimar and colleagues[252] investigated wave speed in a simple threshold cellular automata models with a extended Moore neighbourhood. They found that horizontal and vertical speed differs from diagonal speed in a way that depends strongly on threshold.

Second generation models are a response to the shortcomings of first generation models. Models with thresholds greater than 1 (more than one neighbour must be active to induce an resting cell to the excited state) generates curvature effects if neighbourhood size or the number of excited states is greater than one[90]. Fast and Etimov[72] add dispersion effects by adding a weighing function (where a cells relative contribution to excitability decreases as a function of distance from the wavefront) and another state called the relative refractory period (where cells have a higher threshold than resting cells but can still be excited.) Gerhardt and colleagues[90] add dispersion by making threshold a decreasing function of time. Several methods have been proposed to remove anisotropic effects. Marcus and Hess[166] randomize the grid by leaving many (random) grid points empty in circular neighbourhoods. Kurrer and Shulten[155] select a random subset of a square neighbourhood as the active grid. Several second generation models define diffusion coefficients[92,252] and use cellular automata models to simulate continuous differential equation models.

Cellular automata models offer a level of abstraction over differential equations as a number of underlying variables are encompassed in a single state. For example, the refractory period is not directly represented by differential equations directly but is a function of the interaction of several differential equations. In a cellular automata model, the refractoriness of a modeled cell is simply one of several states of that cell. This allows models to be constructed which encompass many different natural systems. Models by Gerhardt and colleagues[90] and Fast and Etimov[72] were designed as general models of excitable media, that can be applied to specific systems by fitting parameters.

Ito[122] introduced a different approach, called mesoscopic modeling, for formulating general models of excitable media. Coupled differential equation, coupled map lattice and cellular automata models divide up excitable media into a large number of cells, where each cell encompasses the microscopic characteristic of the system. Parameter fitting is difficult if the microscopic behaviour of the system is unknown. Dispersion and curvature relationships cannot be used directly to fit parameters in these models as they are macroscopic features caused by the interaction of many cells. Mesoscopic models [122,123] model the functional behaviour of tissue on a larger space scale. Here, each cell is modeled by the dispersion relation directly. The mesoscopic model is similar to differential equation models that its underlying dynamics are controlled by a continuous equation Equation (1.4), but it still models cell interaction using update rules as opposed to diffusion.

To summarize, several model types are available for modeling propagation in excitable media. Coupled differential equations can be considered continuous in time and space as the equations that drive local dynamics and cell interaction are continuous even if the mechanism for integrating them (discretely on a computer) is not. Coupled map lattices are similar to differential equations but the local dynamics are driven by discrete time (difference equations) processes. First generation cellular automata models are discontinuous in time and space, as a small number of states define cell activity and diffusion is approximated by simple update rules for a small local neighbourhood. Although these cellular automata models are computationally efficient and are general enough to encompass many systems, they do not display curvature or dispersion effects which are ubiquitous features of excitable media. Second generation models have managed to include these effects by increasing the complexity of the cellular automata in regard to local dynamics and update rules.

Several models are addressed in this thesis. Chapter 2 looks at the interaction of two pacemaking sites coupled by excitable media. Here, dispersion and geometric effects are assumed to be unimportant, allowing the excitable media to be approximated by simple rules. Chapter 7 looks at cellular automata and differential equation models coupling excitable media and multiple pacemaking sites. Limitations of the specific modeling approaches used are discussed.

1.3.3  Reentry in excitable media: Spiral waves.

As discussed above, a superthreshold stimulus on an excitable medium can elicit a wave of excitation travelling from the initiation site. In homogeneous media, the wave of excitation moves out as a expanding circle, called a target pattern. Initial conditions in the media can give rise to a different pattern of excitation called a spiral wave. Spiral waves form naturally when a wavefront travels around a pivot point repeatedly reexciting itself. A spiral is created as points further from the pivot have to travel farther to make a circuit, and therefor lag behind the wave at the center.

Spirals waves of excitation are observed in a wide variety of natural systems. Spirals have been observed in the BZ reaction[262], intact[46] and cultured (Chapter 7) cardiac tissue, retinal[101] and cortical[231] neural preparations, and aggregating slime mold [93]. The existence of spiral waves of excitation (also called vortices and rotors) is a general property of excitable media of different nature. The mechanisms of onset and stability of spirals is the subject of extensive ongoing investigation. This section reviews some basic results, however a detailed summary is beyond the scope of this work.

Initiation.

One mechanism was proposed by Krinsky in 1968[149], who studied propagation in cellular automata models. He assumed a region D with a prolonged refractory period in the path of two closely spaced excitation fronts. After the first wave passes, region D remains refractory while surrounding tissue becomes excitable. Passage of the second wave results in a broken wave at the boundary of D. The wave front moves around D until it becomes excitable again, at which point it reexcites the region to form a spiral wave.

Winfree[266] proposed a mechanism of spiral formation in initially homogeneous media. Here, two waves excite an excitable media within a short interval, with the second wave travelling perpendicular to the first. The first wave sets up a gradient of refractoriness in the medium. The second wave breaks on the refractory wake of the first, but passes through excitable media farther behind the first wave. The resultant broken wave forms a spiral as the refractory media becomes excitable again. This mechanism has been used to generate spirals in excised cardiac tissue[46].

A third mechanism is thought to play a significant role in systems displaying spontaneous oscillations leading to frequent target excitation patterns. Here, a region in the wake of travelling wave becomes excited, generating a wave moving outward from the excitation point. The excitation is blocked in the direction of first wave, which results in a semicircular wavefront centered at the excitation point with two free ends. The free ends curve inward to generate a pair of mirror-image spiral waves. This mechanism has been observed in spontaneously oscillating BZ reactions and slime mold preparations[201].

A different mechanism of initiation in heterogeneous media with spontaneous oscillations is discussed in Chapter 7. Here, initiation sites have a high chance of generating a broken wavefront due to heterogeneities in the media. Radius of curvature effects stabilize the wave front as it propagates outward, allowing it to pass through locally heterogeneous media and form a stable spiral wave.

Existence and Geometry.

Tyson and Keener analytically determined the relationship between core size and period for stationary spirals in[246]. Rotating spiral waves generate (approximately) planar waves far from the core. As these waves follow each other periodically, their dynamics are governed only by the dispersion relation, Equation (1.4). Close to the core, where wave front curvature is high, Tyson and Keener used the curvature relation to derive a second constraint between wave speed and period based on the size of the core. The intersection of this new constraint and the dispersion relation gives unique values of planar wave speed and period in terms of core size. The spiral solution derived in[246] assumes that excitation occurs much (infinitely) faster than recovery, which is not valid in real systems[90].

The stability and shape of the spiral depend on the relative size of the core and refractory period of the tissue. If the refractory period is small relative to the period of the spiral wave, then the media ahead of the wave is always excitable, and the path of the wave tip is not affected by refractoriness of tissue. If the inward rotation angle of the tip is large then the core collapses, and the tips trajectory is controlled by refractory media directly ahead of it. In this case, the period of the spiral is simply equal to the refractory period of the media. Collapsed cores are always seen in first generation cellular automata models as radius of curvature effects are not present. There is some evidence that reentry in intact cardiac tissue follows this regime[87,72]. Finally, intermediate states exists, where the core follows a trajectory primarily limited by radius of curvature effects but interacts with partially refractory media ahead of it. This results in a `flower-petal' shaped trajectory called meandering. Computer simulations have shown transitions between circular cores, meandering, and collapsed linear cores as a function of threshold in second generation cellular automata models[72] and coupled differential equation models[276].

Spiral Instability.

Spiral wave break up does not result in annihilation of activity, but multiplies the number of cores resulting in disorganized activity[10]. Two mechanisms of spiral wave annihilation have been studied: meandering of the wave tip in homogeneous media can cause the wave to collide with a boundary[72], and spirals can drift to the boundary in heterogeneous media[206].

A different form of spiral wave annihilation is discussed in Chapter 7. Here, wave propagation is inhibited by a gradual decrease in the excitability of the cells. A cellular automata model with curvature effects is used to simulate the effects of a gradual change in parameters of the system (in this case threshold). Spiral wave annihilation is preceded by a gradual increase in period and core size, which is also observed experimentally.

1.4  Bursting in Other Systems

1.4.1  Conceptual Models of Bursting Behaviour

There are two ways that bursts can be generated in groups of cells: one or more cells in a population have the property of endogenous bursting, or else the bursting arises from interactions between cells.

A cell can be considered a true endogenous burster if it produces bursting pacemaker potentials when completely isolated from its neighbours as well as from bloodborne substances[228]. The most important feature of these cells is that they contain intrinsic conductances that produce slow potential waves. Voltage-sensitive conductances necessary to generate spikes are activated by the slow potential waves to generate bursts.

Bursting activity can be generated in networks of cells without endogenous bursters as a function of the connections between cells. Most network bursters depend on the activity of inhibitory connections on spontaneously firing cells: low frequency patterns are generated by alternate phases of rapid spontaneous firing and suppression of activity by activated neighbours. Other variants use networks to average ionic currents or rely on a fixed delay excitation mechanism to generate bursting. Examples of each class is given in the sections below.

1.4.2  Endogenous Bursters - Mechanisms

There are a large number of theoretical studies on endogenous bursting[26,28,176,195,196,215,230]. In general, these models have a spike generating mechanism of the Hodgkin-Huxley[114] type, and the slow process include slow inward or outward currents which either activate currents which induce spiking or deactivate currents resulting in cessation of spiking[215,230].

An illustrative example of the general form of these models is taken from[230]:

¶v ¶t = (-Ifast(v) - IK(V)(v,n) - Islow(v,Sl)) / cm

(1.7)

where c_m is membrane capacitance, I_fast is a voltage dependent inward current (usually sodium or calcium), I_K(V) is a delayed rectifier potassium current, and I_slow is a slowly varying current. Variables for v (voltage) and n (the gating variable) are fast relative to the burst period, while S_l varies slowly with the burst period. The variables (v,n) generate the spikes within the bursts while changes in S_l cause the model to move from active to silent dynamics by activating an inhibitory outward current. The concentration of S_l gradually increases during the active phase, hyperpolarizing the cell. The downward shift in membrane potential raises the threshold for spiking, and firing stops. S_l returns to baseline during the silent phase and cell firing resumes.

Theoretical and experimental studies have generated several variations of the above equation with different dynamic properties. A classification scheme for the different mechanisms was put forth[215] that groups models into three classes. Classification is based on theoretical analysis of how the slow and fast subsystems interact. In each case, fast variables can be considered at steady state (between bursts) or on a stable limit cycle (during bursts). The system of equations undergoes a bifurcation to switch between the two stable states as parameters controlled by slow changing variables are varied.

Class I models are bistable for a range of S_l, with a stable fixed point and stable periodic cycle. S_l increases during a burst until the trajectory undergoes a homoclinic bifurcation and moves to a stable fixed point. S_l then decreases until the trajectory jumps via a saddle node bifurcation back to the stable cycle. Class I models demonstrate an increasing period through each burst as the trajectory moves toward the homoclinic bifurcation. Class II models do not have bistability, but generate bursting via two slow variables with their own oscillation. A homoclinic bifurcation is crossed at the beginning and end of each burst, causing the period of oscillation to start high, decrease, and then increase with each burst. This characteristic change in period is called parabolic bursting. Class III systems are bistable like class I systems, but here S_l increases during a burst, causing the trajectory to loose stability in a Hopf bifurcation and move to a stable cycle. S_l now decreases until the trajectory undergoes another bifurcation (a saddle node bifurcation of periodic orbits, where one stable and one unstable limit cycle are created) and moves to the stable cycle. The period of oscillation during a burst is not determined by movement to this bifurcation point as it is with homoclinic (Class I and II) bifurcations.

1.4.3  Aplysia Models

A commonly studied endogenous bursting pacemaker is the R-15 cell in the abdominal ganglion of mollusk Aplysia. The general scheme for bursts in Aplysia is as follows: Calcium ions enter the cell with each action potential. During a burst, internal calcium accumulates at a rate faster than the cell can remove it. Eventually, internal calcium concentrations reach levels which both activate a calcium activated potassium outward current as well as reduce calcium current by Ca-induced inactivation of calcium channels. These two effects hyperpolarize the cell, shutting off spiking activity. Inactivity halts calcium entry allowing the cell to slowly lower internal calcium levels. The potassium channels shut when internal calcium concentration reaches baseline levels, and bursting resumes.

A Aplysia model was proposed by Plant[196]. Plants scheme falls under the type II class of parabolic bursting models. The fast variables are contained in sodium, potassium and leak currents, and the two slow variables are calcium current (Ca) and a factor x that controls calcium activation. Slow oscillation in Ca and x moves Ca backward and forward through a homoclinic bifurcation. The rate of this oscillation is controlled in a time constant for activation of Ca; Ca in the model activates and deactivates slowly in response to changes in voltage.

1.4.4  Single b-Cell Models

The pancreatic b-cell is also studied using the endogenous bursting framework. b-cells secrete insulin, which regulates use and uptake of glucose in target tissues. The character of electrical bursting in b-cells is in turn partially regulated by glucose concentrations surrounding the cells. The period of a burst monotonically increases during the burst. Models for b-cells follow a Class I scheme.

Early models for b-cell bursting were developed by Chay and colleagues[28,29]. The proposed mechanisms are similar to those for Aplysia, but here Class I dynamics involving hysteresis and bistability are responsible for switching between silent and active states. Fast variables are controlled by voltage and calcium activated potassium currents. The early models have the internal calcium concentration [Ca²⁺]_i as S_l and I_slow as either a calcium/potassium current[28] or a calcium inactivated calcium current[29].

The models proposed by Chay et. al[28,29] have been challenged by experiments. The early models predict that I_slow varies with a slow sawtooth pattern (increasing slowly during a burst, and decreasing slowly between bursts). Direct measurements of [Ca²⁺]_i showed that it varied as a square wave. This indicates that [Ca²⁺]_i is fast, not slow, relative to the burst period, eliminating calcium as the slow variable. More recent models have set S_l to variables that could respond slowly to [Ca²⁺]_i. Satin and Cook[225] postulate that I_Ca is slowly inactivated by voltage, however Smolen and Kaizer[233] show this hypothesis is incompatible with voltage clamp data. Recently, Smolen and colleagues proposed K(ATP) for I_slow and ATP for S_l[233], but these candidates have been ruled out by experiment[217,232].

1.4.5  Coupled b-cell models

In addition to the uncertainty associated with the mechanisms of bursting, there is debate as to whether isolated b-cells burst at all. b-cells comprise the Islets of Langerhans, a organelle of 10³ to 10⁴ cells. Although bursting is a robust behaviour in whole islets, isolated b-cells spike erratically[230]. It has been proposed that electrical coupling between b-cells is necessary for bursting behaviour[234,229].

One coupled cell model proposes that channel noise in individual b-cells is high, masking intrinsic bursting dynamics. Coupling cells acts to remove noise by averaging, allowing b-cells to burst[29,229]. An alternate hypothesis is based on cell heterogeneity[234]. Bursting dynamics in most models occurs over a narrow range of parameters, with silence at one end of the range and continuous spiking at the other[230]. The heterogeneity hypothesis assumes that currents from different cells with non bursting modes of behaviour sum in such a way that the population as a whole bursts.

The coupled b-cell models depend on cell interactions to burst but they do so in a way different from other network models. The underlying assumption in both coupled b-cell models is that each individual cell bursts under the right (low noise or parametric) conditions.

It is interesting to speculate that bursting in b-cells may be caused by a mechanism similar to that observed in the monolayer. The possibility that bursting in b-cells is driven by patterns of excitation is suggested by the following: First, b-cells only burst in groups and show irregular spontaneous activity in isolation[232]. It is possible that individual b-cells do not have the capacity to endogenously burst even under low noise conditions. Irregular spiking could trigger spiral waves as demonstrated in the cellular automata and differential equation models. Eddlestone and colleagues[61] recorded potential from two sites within an islet. They discovered that bursting within each islet is synchronized, though the spikes are not. More recently waves of electrical activity have been observed in Islets of Langerhans[192] using optical recording techniques. These observations indicate he islets are not isopotential, and raise the possibility that spiral or scroll waves exist in the islet. Finally, gap junction conductance decreases significantly during bursts in b-cell clusters, which may produce fatigue like effects, and heterogeneity in gap junctional conductance exists in the islet[169], which can give rise to unidirectional propagation.

1.4.6  Network Bursters

Bursting can be generated as a function of the connections between cells in the absence of endogenous bursting mechanisms. Network driven bursters are often studied in neuronal systems, where they are commonly referred to as central pattern generators (CPGs). Network bursters in these systems are characterized by a known set of neurons with excitatory or inhibitory connections which have well defined roles in generating bursting behaviour. Network bursters are well studied in invertebrate systems such as the Lampry pattern generator for locomotion[33], the Tritonia swim oscillator[257], and the Locust flight motor pattern generator[258] but are also studied in vertebrates in context of respiratory rhythm generation[200].

A simple example of a network burster is the sea slug Tritonia swim oscillator. Here three groups of premotor interneurons are responsible for bursting. The principle synaptic mechanisms underlying the swim oscillator is reciprocal inhibition with delayed excitation. Three broad groups of cells are involved in this system; dorsal swim interneurons (DSI), ventral swim interneurons (VSI) and a cerebral cell (C2). The swim cycle begins when the DSI depolarizes and begin to fire. These cells inhibit the VSI and excite C2. C2 fires synchronously with DSI, and stimulate VSI, which eventually start firing despite inhibition from DSI. The VSI acts to inhibit DSI and C2, terminating their bursts. VSI stops firing as stimulation via C2 is lacking, releasing DSI from inhibitory input, allowing the cycle to start again.

Inhibition of one group of neurons by activity of another is observed in complex networks. A mechanism similar to Tritonia is seen in the vertebrate respiratory rhythm generator, which involves six interconnected neuronal groups[200]. The lamprey swim pattern generator may involve chains of dozens of interconnected oscillators[33,34]. The lamprey pattern generator consists of groups of endogenously bursting elements that use inhibition to synchronize dorsal and lateral elements 180 deg out of phase.

1.4.7  Reentrant bursting models

Bursting dynamics are sometimes symptoms of disease states. An example is repetitive paroxysmal tachycardia, which is characterized by the onset and offset of rapid cardiac activity. The mechanism of bursting in this system is not well understood. The tachycardia could be set by the rhythm of a spontaneously beating focus. Alternatively, the tachycardia could be caused by activation along a reentrant pathway (see Chapter 1.) In this case, both termination and initiation could be caused by physiological properties of the conductive tissue itself.

A simple experimental model for paroxysmal tachycardias is presented in Kunysz and associates[154]. The experimental model consists of spontaneously beating aggregates of chick embryonic atrial heart cells. A reentry loop is modeled by causing the aggregate to be excited a fixed time after it fires. The fixed delay is set so that the aggregate is forced to beat faster than its control cycle length. Prolonged rapid activity results in the build-up of ion species (thought to be sodium) which act to increase the cycle length of the aggregate (the increase in cycle length is called overdrive suppression[153].) The reexcitation pulse falls at progressively earlier times in the cycle of the aggregate until excitation fails. The aggregate beats again after a few seconds (its reset intrinsic cycle length) and the bursting pattern resumes.

The above model assumes the reentrant loop has properties that are invariant over time. Termination of reentry can occur due to changes in the conductive tissue itself. Dispersion effects as well as time dependent changes in the refractory period and action potential duration can cause excitable tissue to block under rapid stimulation conditions. There is a large body of work discussing the dynamics of reentry and block on rings of cardiac tissue[210,43,249,99]. Frame and colleagues[86,18] conducted experiments on rings of cardiac tissue of fixed and variable (using a fixed delay protocol[86]) diameters. Reentry becomes unstable and eventually blocks for small rings. Similar dynamics are observed in conductive strips made from cultured embryonic chick cardiac cells stimulated using a fixed delay protocol (unpublished results.) Gradual decrease in fixed delay results in oscillation and finally block. This mechanism for termination of reentry coupled with spontaneous initiation of reentry (via unidirectional block in the ring) could generate bursting patterns in a ring of excitable tissue with external pacing.

1.5  The Heart as an Excitable Medium.

The propagation of electrical activity in cardiac muscle involves the interaction of different ion species across a combination of active and passive ion channels and diffusion of charge through a heterogeneous substrate with dynamically changing conductances[71,140]. Despite complexity inherent in conduction at microscopic scales, the heart has been approximated as a continuous excitable media in many studies. A variety of cardiac tachycardias have been attributed to formation of large scale patterns of excitation such as the interaction of several oscillators or the formation and break up of spiral waves. The following sections outline some of the differences between cardiac tissue and simple excitable media, and discuss reentry and initiation of ectopic pacemaking activity in diseased heart.

1.5.1  Conduction in Healthy Hearts.

The chambers of the heart contract in a precise sequence controlled by the oscillatory and conductive properties of a number of different cells types. All cardiac cells are excitable and some can sustain spontaneous oscillatory activity. In healthy heart, the cardiac rhythm originates from a small area (1-3 mm²) called the sinoatrial (SA) node located in the upper right atria. The spontaneous activity in the SA node acts to suppress the activity of other groups of slowly beating cells[248]. These latent pacemakers can drive the activity of the heart if the the sinus node fails to beat.

The impulse originating from the SA node is conducted through the atria, which causes them to contract and move blood to the ventricles. The impulse then passes through the atrioventricular node, which is the sole conductive pathway connecting the atria and ventricles. The low electrical coupling and slow rise time of AV node cells imposes a delay, which allows efficient filling of the ventricles. Once the impulse passes the AV node, it flows through rapidly conducting Purkinje fibers to the ventricles. The ventricles forcefully contract, completing the cycle.

There is strong spatial specialization in cells that constitute the myocardium. Action potential characteristics, such as duration, upstroke velocity, and threshold, change with anatomical location. Action potential duration is a factor controlling the refractory period and duration of contraction of cardiac cells[182]. The action potential duration is different for every cell type in the heart. The upstroke velocity (which is a factor controlling conduction velocity) is fast in ventricular and atrial cells (300-400 V/s) and slow in AV and SA nodal cells (4 - 30 V/s)[184]. In addition, cells in different regions have different excitation thresholds and can display graded or partial excitation in response to superthreshold stimulus[184]

There are several additional differences which distinguish the heart from the simple homogeneous excitable media examples discussed in Section 1.3.2. Cardiac tissue has appreciable thickness, allowing for propagation in three dimensions, and cells are aligned in a preferred orientation. Conduction in cardiac tissue can also be affected by heterogeneities at the cellular level. Finally, cardiac tissue can display long term effects based on its stimulation history. These issues are discussed below.

1.5.2  Anisotropy.

Simulations with models with different diffusion coefficients gives elliptical targets from point source excitation, and spirals `stretched' in one axis. Patterns created by anisotropic simulations can be rescaled to have (approximately - see reference [270]) radial symmetry. Several interesting points relating to changes in conduction velocity along a elliptical target due to change in front shape are discussed in reference [270]. Anisotropy can give rise to subtle changes in regards to the direction of conduction, interaction of the wave with boundaries, and in the critical radius of curvature. The effects of anisotropy on conduction is an area of ongoing research [270,158].

1.5.3  Propagation in 3D.

Spirals in three dimensions manifest as scroll waves. A simple scroll wave can be visualized as a loosely wound sheet with its ends connected. The spiral core becomes an elongated filament connected at both ends. A scroll wave has a circular filament, however complex geometries with highly twisted filaments are possible. Simulations show that the filament follows a drift that depends on both the curvature and the twist of the filament. In addition to moving through the medium, the filament can contract and expand in length. Scroll waves extinguish themselves if the filament collapses to a point or expand past the mediums boundaries[205,138].

1.5.4  Local Heterogeneity.

Discontinuities can be due to the distribution of gap junctions, and the presence of small clefts and anatomical boundaries such as blood vessels. The effects of heterogeneity on conduction has been investigated in models[219,237] and in experimental systems[71]. Experimental studies have shown that intracellular clefts (on the order of one cell in length) cause increases in the upstroke velocity at either side of the cleft and reductions of the upstroke velocity as the waves of excitation passes through narrow channels between clefts. Changes in the upstroke velocity are implicated as sites of reduced conduction velocity and unidirectional block.

The effect of small discontinuities is considered in greater detail in Chapter 7. An experimental system is presented which is highly heterogeneous at microscopic scales but appears to behave as a continuous (isotropic) medium at macroscopic scales. The relationship between discontinuous conduction and continuous theoretical models is unclear. Winfree[270] gives several `rules of thumb' that test whether a continuous representation is valid at macroscopic scales (see Section 7.4.3). The continuity conditions are met in normal intact cardiac tissue and may hold for the model system under investigation as well. The distinction between continuous and discontinuous theoretical models is partially avoided by using a discontinuous cellular automata model.

1.5.5  Overdrive Suppression and Fatigue.

Overdrive suppression has been studied in a variety of spontaneous cardiac tissues preparations. Overdrive suppression results in a gradual increase of the intrinsic cycle period in response to prolonged rapid pacing. Overdrive suppression has been modeled in aggregates[153] by assuming that each AP produces an outward (hyperpolarizing) current which decays exponentially with a long time constant.

Fatigue has been extensively studied in AV nodal preparations, where prolonged rapid stimulation produces pronounced effects. The dominant effect of fatigue on function is to shift the AV nodal recovery curve¹ upward, increasing the conduction time. Fatigue has been modeled by assuming that a fatigue term increases with each AP and decays away exponentially with a long (30 seconds) time constant. The fatigue term is summed with the recovery curve, increasing conduction time. Fatigue is discussed in greater detail in Chapter 7.

1.5.6  Conduction in Diseased Heart: Reentry.

In healthy myocardium, rapid conduction velocity and long refractory periods protects the medium from reexitation from the propagating wave front. For example, assuming a conduction velocity of 50 cm/s[143] the conduction time across the ventricle is 25-30 ms, which is 8 times shorter than the ventricles refractory period.

In diseased heart, changes in conduction velocity and refractory period can result in conditions that allow reentry. The most common form, proposed by Mines[175] involves reentrant conduction around a non excitable obstacle. The conceptual model involves conduction along a one dimensional ring of tissue. Reentrant excitation is possible because the tissue ahead of the exciting wavefront regains excitability. This can occur if the length of the refractory zone (sometimes called the wavelength of excitation) is shorter than the circumference of the ring. Therefore, shorter refractory times, slower conduction velocities, and longer path lengths increase the ability of the circuit to support reentrant excitation. Mines's ring model may explain forms of reentry with distinct pathways such as Wolff-Parkinson-White syndrome[142], but is not easily applied to intact myocardium.

Allesie[3,4] showed that reentry is possible without an anatomical obstacle. Here, the mechanism of reexcitation is roughly analogous to that of spiral waves as discussed in Section 1.3.3. The mechanism proposed by Allesie is called `leading circle reentry' where a core region displays a slightly elongated refractory period than the surrounding region. Partial excitation of the core result in increased refractoriness in these cells. This mechanism has the impulse circulating around the smallest possible pathway at the highest speed supported by the tissue. It is unclear whether the length of the reentrant pathway in this conceptual model is a function of the circumference of the functionally refractory core or is limited by refractory tissue directly in front of the excitation front. Another related mechanism does not require heterogeneity of refractory times; the progress of a wave of excitation is only limited by refractory tissue directly ahead of it. The core (in this case a collapsed line) circumference is the minimum path length allowed by the refractory period. The period of rotation is therefore equal to or slightly greater than to the refractory time of the cells. As discussed above, the period of spiral wave rotation can also be a function of the excitability of (fully recovered) cells ahead of the wavefront. Here, the diameter of the core is a determined by the curvature angle of the spiral tip, and is not limited by the refractoriness of the myocardium. Both leading circle and excitation limited reentry can in theory result in an excitable gap²[270]: a region of tissue ahead of the wavefront that is excitable. Excitable gaps have been observed experimentally in cardiac preparations[119].

1.5.7  Changes in Refractory and Excitability Properties in Diseased Hearts.

Biological models for reentrant arrhythmias in intact hearts involve inducing ischemia (a reduction of blood flow.) Ischemic conditions can be produced by ligation of a coronary artery, reperfusion of acute ischemic myocardium, or production of chronic infarction[142]. These models correspond to frequent clinical situations where tachyarrhythmias are encountered. Ischemic conditions alter ionic properties of cardiac cells and the functional properties of the tissue, allowing reentry.

Conduction velocity can be decreased by decreasing cell-cell coupling or by decreasing cell excitability. Decreased cell coupling is known to occur in ischemic tissue, and is probably due to an increase in intracellular Ca⁺⁺[142,165]. Increases in Ca⁺⁺ occur at later phases of ischemia (10-20 minutes after arrest of coronary flow[144].) Other changes occur earlier in ischemic models and are more likely to play important roles in induction of reentry. One of the first changes is a increase in excitability, brought on by a depolarization of the membrane and a decrease in the amplitude and duration of the action potential. Action potential shortening and membrane depolarization decreases the refractory period of cells, allowing reentry in ischemic tissue. On the other hand, infarction results in decreased cell coupling, while maintaining normal action potential shape[142]. This allows reentry to occur by decreasing conduction velocity.

1.5.8  Spontaneous Activity in Diseased Heart: Focal Mechanisms.

Injury currents are due to local potential differences between ischemic and non ischemic cells (recall ischemia depolarizes cells and decreases action potential duration.) If ischemic and healthy cells are electrically coupled the difference in membrane current produces local circuits that source excitation current to cells, inducing premature APs.

DADs are related to Ca⁺⁺ overload which occurs late in ischemia. The cellular Ca++ overload causes a cyclic release of calcium from intracellular stores[142,67]. Cyclic changes in calcium levels cause cyclic changes in membrane potential, as increased calcium levels cause inward movement of positive charge (electrogenic inward movement of sodium via the sodium/calcium exchanger and inward current through calcium activated cation channel are thought to play a role.)

1.5.9  Cultured Monolayers as Models for Cardiac Conduction.

Israel and colleagues measured intercellular time delays in monolayers of 10-day chick embryo myocytes grown on a microelectrode array. Their aim was to test the theoretical model of Diaz and coworkers[55,219] which indicated that conduction within a cell is an order of magnitude higher than the average conduction velocity over many cells, and predicted delays of 0.45 ms between adjacent cells. Israel and coworkers found delays as long as 0.41 ms between adjacent cells under conditions of reduced cell-cell coupling.

Fast and Kleber[69] developed an experimental system that allows them to control the pattern of cell growth on a collagen substrate. The collagen substrate is brushed so that cells plated on the collagen align themselves in a parallel fashion. Impulses from the monolayers are recorded with a photodiode array. Their experimental system allows the study of anisotropy and heterogeneity in cell culture (see Sections 1.5.2, 1.5.4 and 7.4.3 for further discussion.)

Fast and Kleber[70] calculated the minimum diameter of a strand of tissue necessary for conduction into a large area. Block of conduction during conduction from a small to large tissue volume occurs due to critical curvature beyond the point of expansion (see Equation (1.2).) Rhor and Salzberg[213] confirmed this theory with optical mapping studies of conduction in cultured neonatal rat myocytes. Rhor and colleagues developed a cell culture system that allows for patterns of cell adhesive collagen to be laid out using photolithographic techniques.

Rhor, Kleber and colleagues[214] altered cell-cell coupling and cell excitability on cultured strands of neonatal rat myocytes. The motivation for their study was to investigate slowed conduction, which is a contributing factor for the creation of reentrant arrhythmias. By lowering conductance the experimenters were able to slow conduction from 47 cm/s to less than one cm/s, a speed which would allow reentry in tissue areas of less than 1 mm².

The experimental work detailed in Chapter 7 investigates the effects of low coupling and heterogeneity on macroscopic conduction. Poorly coupled and spontaneously beating monolayers generate broken wave fronts which evolve to spiral waves. Heterogeneities in local cell-cell coupling may play an important role in generating unidirectional block at the site of initiation. Alternative theories are discussed in Chapter 7.

1.6  Mapping Techniques.

The electrical activity of interacting cells has been of interest to researchers for many years. The primary methods for detecting electrical activity has been extracellular electrodes that measure potential differences between different areas and intracellular microelectrodes that penetrate the cell and measure potential across the cell membrane. Arrays of extracellular electrodes remain popular for detecting patterns of propagation in whole[64] and cultured[120] cardiac preparations, and slices[100] and cultured neural[183] preparations.

Recently, optical methods have been developed that give direct information on membrane potential or the concentration of particular ion species without requiring cell penetration. A sample is stained with a voltage or ion sensitive dye whose fluorescence, absorption, or bifringence changes linearly with changes in the measured quantity in the cell. Over 1800 substances have been tested for optical responses, yielding 200 dyes with significant optical changes with acceptable levels of phototoxic damage and pharmacological side effects. The mechanism producing optical changes controls the response time of the dyes. Dyes either act slowly, by diffusing in response to changes in the local environment, or act rapidly by changing orientation or conformation.

A number of fast response dyes have been used to map potential changes in a variety of preparations. Early mapping studies were in giant squid axon[48,243] and invertebrate central neurons[223,222]. These techniques have more recently been applied whole brain[106] and sections of the visual cortex[105] and olfactory system[199]. Intact[95,156], sectioned[220,46] and cultured[213,71] cardiac preparations have also been studied. Reentry has been initiated and observed optically in epicardial slices and whole heart[46,95]. The effects of structural anisotropy and heterogeneity have also been studied using optical mapping techniques[71]. Propagation abnormalities due to extracellular stimulation have been investigated in intact heart[94]. The effect of geometry on block[213] has been investigated in cultured preparations, and measurements of the effects of wavefront curvature on conduction velocity and block[141,24] has been investigated in intact and sectioned preparations. Optical mapping techniques have been used to track changes in individual ion species such as Ca⁺⁺[207,66,160] and metabolites (i.e. NADH[221]) in a variety of substrates.

All optical mapping studies involve specialized illumination and detection apparatus. In most studies, stained samples are illuminated over a relatively large area. Appropriate sources of illumination vary with the characteristics of the dye used, but is usually a tungsten halogen or mercury arc lamp fitted with appropriate filters to select optimal excitation wavelengths. Fluorescent light from the preparation is filtered to remove background reflected light. A detector with multiple elements, such as a photodiode array or CCD, is placed in the objective image plane and detects the signal of interest. The image is then signal processed and stored. Variations on this general scheme have used a single detector with a scanning excitation light source[181]. The specifics of the above schemes are discussed in detail in Chapter 3.

1.7  Oscillators and Excitable Media, Experiments and Theory.

Oscillators responses to periodic pacing have been studied extensively. A extension of these studies to intact cardiac preparations is limited to direct stimulation of pacemaking sites by extracellular electrodes. Simple interaction of two pacemaking sites located a distance from each other requires that the intervening excitable media be factored in. In some cases, theoretical understanding of the rhythms produced by a two oscillator system can be obtained by using gross approximations of the effect of the excitable media on the interaction of the oscillators.

Chapter 2 investigates the interaction of two pacemakers in an excitable media in the context of modulated parasystole. Parasystole is a benign arrhythmia caused by intermittent beats originating from an ectopic pacemaking site located some distance from the sinus node. In pure parasystole, the pacemakers do not phase reset each other, but the effects of either pacemaker can be masked if the surrounding excitable tissue is refractory. Modulated parasystole assumes the ectopic pacemaker can be reset by the sinus rhythm.

The modulated parasystole model is presented as a circle map which incorporates the Poincaré oscillator as the ectopic pacemaker. The Poincaré oscillator has well understood dynamics under conditions of direct forcing. Addition of assumptions associated with an excitable media introduces a discontinuity in the circle map which dramatically alters the dynamics. Investigations into the behaviour of the discontinuous parasystole map led to a general analysis of a broad class of discontinuous maps which have not previously been characterized. The class of maps investigated are chaotic but display a type of periodicity (called banded chaos) on a macroscopic scale. The analytic results from the chaotic maps are extended by carrying out simulations on the more complicated modulated parasystole map, which also displays quasiperiodic and periodic dynamics for certain parameters.

Experiments on oscillators coupled by an excitable substrate are performed using monolayer cultures of embryonic chick cardiac cells. The pacemaking behaviour of embryonic chick cells have been extensively documented. Aggregates of embryonic chick cardiac cells spontaneously beat ³ and their dynamics under periodic forcing are relatively well understood. Aggregates are considered to be isopotential systems and therefore do not incorporate effects of conduction into their dynamics. However, conduction can be observed in spatially extended cultured monolayers of cardiac cells.

Here monolayers are used to investigate the interaction of many oscillators dispersed in space. Monolayers of coupled pacemaking cells can potentially provide insight into fundamental issues related to initiation and conduction in spontaneously active media. The behaviour of the model system may also help understand diseased cardiac tissue, which displays spontaneous oscillations when ischemic.

Cultured cells couple to each other by forming gap junctions after they are plated. Investigations of cultured cells `by eye' (under a microscope) Suggests that functional connectivity is increases over time. Older cultures (2-3 days) appear to fire in relative synchrony. A simple explanation is that a periodic target wave originating from one point in the culture acts to entrain the entire population. Younger cultures (1 day old) display complex dynamics. Cells occasionally synchronize, but often beat in a disorganized fashion. The dynamics of the entire system are impossible to characterize by eye as the field of view (where movement artifacts are visible) is only a few tenths of a millimeter in diameter. The dynamics of loosely coupled monolayer cultures is investigated using macroscopic mapping techniques and mathematical models in Chapter 7.

position/velocity phase space.

Chapter 2
Pacemaker Interactions: A Mathematical Model for Parasystole

2.1  Introduction.

A circle map is a function that maps the circumference of a circle onto itself. Such functions have been intensively studied in recent years, both for their mathematical interest, as well as their relevance to physical and biological problems. Our interest in circle maps has been stimulated by evidence from both experimental and clinical studies that circle maps are appropriate approximations for the periodic stimulation of spontaneous cardiac oscillations.

Circle maps are often classified by their topological degree. Say that f is a point on the circumference of a circle and f(f) is a circle map. As f increases from 0 to 1, f(f) traverses the circumference of the circle M times where M is the topological degree. It is usual to assume that M is an integer. In this case the circle map is continuous, as f(0) = f(1). A large body of work has analysed the bifurcations of continuous circle maps.

Circle maps can also be discontinuous, but the bifurcations of these maps are not as well studied. One type of discontinuous circle map that has been investigated is a map that is monotonically increasing except at the discontinuity. Keener[135] has given results relating the possible dynamics of these maps to the nature of the map at the discontinuity. However, Keeners analysis for maps that are noninvertible at the discontinuity applies only to maps with at least one fixed point. The dynamics of noninvertible maps without fixed points have been investigated for piecewise linear maps in connection to the study of neuronal dynamics[21] and analog to digital converters[73]. Our present interest is in the application to modulated parasystole, which can be represented as a discontinuous circle map with no fixed points.

Modulated parasystole is a cardiac arrhythmia caused by an abnormal (ectopic) pacemaking site located in the ventricles of the heart that competes with the normal rhythm generated in the sinus node. The two pacemakers do not fire independently of each other, but rather interact in a complicated fashion depending on the electrophysiological properties of the cardiac tissue. Despite the apparent complexity of the problem, the reduction of the problem to mathematical questions involving circle maps, has enabled one to make theoretical predictions, that in some cases can be confirmed in clinical studies. Previous theoretical studies of modulated parasystole employed discontinuous circle maps that for some parameter values displayed deterministic chaos[39].

In this paper we consider the possible dynamics of a general class of discontinuous maps with no fixed points and apply our findings to a mathematical model for modulated parasystole. In Section 2.2, we review relevant terminology and definitions concerning circle maps. In Section 2.2.1, we consider the dynamics of a general class of discontinuous, noninvertible circle maps with slope greater than one and no fixed points. In section 2.2.2 we apply these results to understanding the bifurcations in parameter space of a piecewise linear map. During the course of this investigation we became aware that some of the results in this section have already been presented using a different approach employing symbolic dynamics[73]. In section 2.3 we describe a theoretical model for modulated parasystole using discontinuous circle maps without fixed points but with slope less than one for part of its domain. Results of simulation of the theoretical model are presented in Section 2.3.1. The results are discussed in Section 2.4.

2.2  Circle Maps.

We are interested in the dynamics of discontinuous circle maps. Such maps can be defined in the following way. Consider the functions F :[0,1)® R, and G =F (mod 1), where F is a single valued function. Then the map G is a circle map, G:S¹ ® S¹.

Given an initial condition f_i, f_i+1=F (f_i). Consider the sequence of f₀,f₁, ¼, found from iteration of a circle map. There are three important concepts that are necessary to characterize the dynamics: periodicity, rotation number, and Lyapunov number.

There is a periodic orbit of period n if f_n=f₀ but f_j ¹  f₀ for j=1,2, ¼, n-1.

The rotation number r(f₀) is defined[14]

r(f0) = limsup n ® ¥  1 n n å i=0  F (fi) - fi.

In applications in which circle maps are models of periodically forced nonlinear oscillators, the rotation number reflects the relative frequencies of the spontaneous oscillator and the periodic forcing. The rotation number is independent of the initial conditions for any monotonic circle map, but may depend on initial conditions, for a nonmonotonic map[88]. If r is sensitive to initial conditions the set of r(f₀) will span a non zero interval, called the rotation interval.

Another concept associated with rotation number is the resonance interval, or mode-locking interval, I_p/q[113,6,7,129]. I_p/q can be defined as an interval on an arbitrary parameter a such that r[f _a(f)] = p/q[129]. For a large class of circle maps, resonance intervals will have a Farey order in respect to the rotation number on the parameter a[8]; for two resonance intervals with rotation numbers [N/M] and [(N¢)/(M¢)] the most prominant resonance interval between them will have r = [(N+N¢)/(M+M¢)].

The Lyapunov number, l is defined

l = lim N ® ¥  1 N N å i=l  ln|G ¢(fi)|

(2.1)

where negative values of the Lyapunov number indicate periodicity, and zero represents quasiperiodicity. Although it is usual to define deterministic chaos based on a positive Lyapunov number, several studies define deterministic chaos based on a positive rotation interval[135,27,88]. These 2 definitions are not equivalent. We will consider both the Lyapunov number and the rotation interval in the following.

2.2.1  The Dynamics of Discontinuous maps.

1.   F(0) < F(1). 2.   F¢(f) > 1, f Î [0,1]. 3.   F(f) ¹ f,   f Î [0,1].

(2.2)

A sketch of F(f) is shown in Figure 2.1. Keener discussed the dynamics of discontinuous circle maps F that are noninvertible at the discontinuity, and the slope everywhere else is greater than one. If F has at least one fixed point then the rotation number r is sensitive to initial conditions resulting in a positive rotation interval. In this section we discuss possible solutions of F, a discontinuous noninvertible circle map with no fixed points. We use these results to analyse the bifurcations of a simple two parameter discontinuous map, which we will later use to help us understand the bifurcations of the more complex parasystole map.

F has F¢(f) > 1 for all f and will therefore be associated with a positive Lyapunov number and have chaotic dynamics. Although F will have no stable periodic orbits, the density of the trajectory can have a regular structure. Figure 2.2 shows a bifurcation diagram for the circle map,

ft+1 = k(f) = Mft+a+ b 2p sinft   (mod 1 )

(2.3)

where M > 1 and a,b are such that k(f) is of class F. For some values of a and b, the trajectory is restricted to discrete intervals, or bands, on [0,1]. The trajectory moves in a periodic fashion between the bands, but the trajectroy is chaotic within the bands. This behaviour has been defined as `noisy periodicity' or `banded chaos' by Lorenz, who observed it in the logistic map [163]. Here we 1) describe conditions for banded chaotic behavior in discontinuous circle maps, and 2) argue that if the map displays banded chaos, then the rotation interval goes to a point and has a rational value.

Consider a sequence S generated from an unstable period N orbit of F, S = f¹, f², ... f^N-1, f^N . Between any two points f^a and f^b, where f^b immediately follows f^a in S, there is a third point f^c that maps within i iterates, i < N, to the discontinuity at 0(1), as shown in Figure 2.3A. The trajectory will return every N iterates to an interval [a,b] on [0,1] where a > f^a and b < f^b if

i) FN(ft ­ fc) < fa ii) FN(ft ¯ fc) > fb

(2.4)

where F^N(f₀) is the N^th iterate of F starting with initial f₀. If conditions (i) and (ii) are met, then a trajectory will be restricted to regions between points on an unstable period N orbit, and will move in a periodic fashion between these regions.

We can further characterize the dynamics of F by constructing maps G_a and G_b defined as

Ga(x) = infx £ y F(y),   Gb(x) = supx ³ y F(y)

which bound F as shown in Figure 2.3C. If G_a and G_b are continuous monotonic circle maps they each have a single r for all initial f₀[110]. Since G_a and G_b have a region where the function is a constant (a `flat') , and are increasing functions with slope greater than one everywhere else, the rotation number assumes rational value[19]. The rotation interval I of a nonmonotonic map f(x) is bounded by the rotation numbers of two monotonic maps G_a,G_b[88] where

Ga (x) ³ f(x) ³ Gb (x),   I[f(x)] = [r(Ga(x)),r(Gb(x))]

(2.5)

Figure 2.3D shows G_a and G_b, where f_t+N is plotted against f_t. G_a (or G_b) will have a stable period N orbit if the flat of the N^th iterate of G_a(f) (or G_b(f)) intersects with the f_t = f_t+N line, as shown. It is easy to see that conditions (i) and (ii) above are equivalent to the conditions for G_a and G_b to have a stable period N orbits. In addition, conditions (i) and (ii) ensure that G_a and G_b share unstable period N points f^a and f^b, as shown. Continuous monotonic increasing maps have a constant rotation number for all periodic points, stable or unstable. Consequently, both G_a and G_b will have the same rotation number as the orbit including f^a and f^b. If conditions (i) and (ii) are met for F, we can conclude that the rotation interval I[F(f)] goes to a point, where r(F) = r(G_a) = r(G_b). As r(G_a) and r(G_b) are rationals[19], r(F) will assume a rational value.

The above results allow us to determine whether a given map of class F will have a regular structure in the density of the trajectory over many iterates. If G_a and G_b both have N:M orbits, conditions (i) and (ii) will then hold at N locations on [0,1]. The trajectory will then be restricted to N discrete regions, or bands, on [0,1]. The existence and the number of bands as well as r(F) depends on parameters controlling the structure of F. In the following we use G_a and G_b to study the bifurcations of a simple two parameter discontinuous map of class F.

2.2.2  Bifurcation Structure of a Discontinuous Map

The equation

ft+1=Mft+a    (mod 1),(1 < M < 2),(0 < a < 1), (M < 2-a)

(2.6)

is a noninvertible discontinuous map with slope M always greater than one. The discontinuity at f = 0 (=1) is such that points f_t+1(f_t ® 1)  >  f_t+1(f_t ®0). If M < 2-a Equation (2.6) will have no fixed points on [0,1]. The bifurcations of Equation (2.6) have been investigated by Feely and Choa[73] using a symbolic dynamics approach. We use different techniques to obtain a more general understanding of the bifurcations in the parameter space for Equation (2.6), which will allow application to more complex maps.

Figure 2.4 shows the results of iterating Equation (2.6) for different values of M and a. The trajectory can either fill the entire phase space, or be bounded into narrow bands. As the slope of the map is everywhere greater than 1 the trajectory within the bands is necissarily chaotic, however the orbits move in a periodic fashion between bands. The location and number of bands depends on parameters M and a.

missing...

For Equation (2.6), G_a and G_b can now be given by

Ga (x) = infy ³ x F(x) = ì ï ï ï ï ï í ï ï ï ï ï î ft+1 = M ft + a       if ft[0,     1 M ] ft+1 = a       if ft[ 1 M ,    1]   Gb (x) = supy £ x F(x) = ì ï ï ï ï ï í ï ï ï ï ï î ft+1 = M + a - 1 if ft[0,   M-1 M ] ft+1 = M ft + a if ft[ M-1 M ,  1]

(2.7)

Figure 2.5A and 2.5B show entrainment zones (equivalent to Arnold tongues) calculated for monotonic maps G_a and G_b. Figure 2.5C shows the entrainment zones of G_a superimposed on that of G_b. The regions where tongues of the same rotation number overlap correspond to regions in the parameter space where conditions (i) and (ii) are met for Equation (2.6) and banding of the trajectory occurs.

We can immediately state the following results for Equation (2.6):

1. Parameters that result in banding are organized into discrete zones, which we will refer to as banding zones, in the parameter space.
2. The rotation number within each zone is a constant rational value M/N, and the trajectory will be restricted to N bands on [0,1].
3. The banding zones have a Farey ordering in the parameter space with respect to the rotation number.
4. The banding zones cannot overlap.

The above results are a consequence of the structure of the entrainment zones for the bounding maps G_a and G_b. Banding occurs in regions in the parameter space M,a where r(G_a) = r(G_b). For constant M, this will occur only when equivalent resonance intervals of G_a and G_b overlap. Therefore, a (M/N) banding zone will occur only on a subset of a that gives rise to a M/N resonance interval of G_a. As G_a has a Farey order with respect to r this will also be true for the banding zones for Equation (2.6). Also as G_a is a continuous invertible map for all M and a, the resonance intervals of G_a do not overlap each other. We conclude that the banding zones for Equation (2.6) do not overlap in the parameter space.

The above results are not dependent on the local structure of Equation (2.6). These results should hold for a class of circle maps F(f) = f(f_i) + t if the entrainment zones for the bounding maps (G_a and G_b) of F(f) exist and have a Farey ordering in the parameter space.

We now discuss the conditions for the existence of a given N/M banding zone for Equation (2.6). Consider the f^N vs f map in Figure 2.6. Conditions (i) and (ii) will hold if the slope of the line segments (M^N) is less than two. Therefore, all N banded zones, with r = p/N, p < N, will first appear as M decrease below A[N]2, and disappear as M approaches 1. Also, for a given M, the highest order banding zone (N¢) is given by the integer part of log(2)/log(M). As M ® 1, N¢® ¥, and as there are at least N¢ different banding zones, the number of different banding zones also goes to ¥ as M decreases to 1.

Figure 2.7 is a plot created by determining the differences in the period (N) of the trajectory for G_a and G_b for many values of M and a. Each point in Figure 2.7 is a colour representation for the difference in the periods N for G_a and G_b (D[N(G_a),N(G_b)]). The intensity of the colour used is proportional to the differences in the period . (Here D[N(G_a),N(G_b)]=0 is black). The banding zones are clearly visible. Additional black regions are caused by the the overlap of N:p and N:q zones, where p ¹ q.

2.3  A Theoretical Model for Parasystole

Parasystole is a cardiac arrhythmia in which a normal heartbeat generated by spontaneous electrical activity coming from the normal pacemaker of the heart located in the sinus node in the atria, interacts with an abnormal ectopic pacemaker located in the ventricles. The normal sinus beat acts to phase reset the ectopic pacemaker.

A mathematical model for parasystole was proposed by Jalife and coworkers, and several subsequent studies have interpreted these studies in terms of circle maps. In this paper we follow the development of Courtemanche et al.[39] The mathematical model of modulated parasystole used in this study is based on the following rules:

• The sinus and ectopic pacemakers each have their own intrinsic periods, designated t_S, t_E, respectively.
• A beat initiated by the sinus node may modulate the time when the ectopic pacemaker next fires;
• There is a refractory period (q) after each sinus beat.
• After each ectopic beat the next sinus beat is blocked.

To express this model in mathematical terms, we call f_i the phase of the sinus beat in the ectopic cycle. The ith sinus beat sinus beat at phase f_i leads to a phase resetting of the ectopic pacemaker so that immediately after the sinus beat the phase of the ectopic pacemaker is g(f_i), where g is a function (sometimes called the phase transition curve) that is either measured experimentally or modelled in an approximate fashion. The model for modulated parasystole can be expressed in terms of the following difference equations:

ì ï ï ï ï ï ï ï í ï ï ï ï ï ï ï î ft+1   = a(ft) a(ft) = ft + t if ft [0, t - q TE ] a(ft) = G(ft,B) + t if ft [t - q TE , 1]

(2.8)

where G(x,B) gives the degree of modulation. The general form of this model was first used by Courtemanche et al.[39], where phase resetting curves G(x,B) were elucidated from clinical data. In this study G(x,B) is modelled by the Poincaré oscillator, which has been studied as a model for biological oscillators[274]. The Poincaré oscillator can be given by:

cos(2pG(f,B)) = B +cos2pf (1 +2Bcos2pf + B2)[1/2]

(2.9)

The Poincaré oscillator may not be a realistic model for phase resetting in parasystole, however it does allow constant changes in stimulation amplitude. The periodically stimulated Poincaré oscillator can display periodicity, quasiperiodicity, as well as chaos under different conditions[274]. Chaotic dynamics are observed at high stimulation strengths (B > 1) in the periodically stimulated Poincaré oscillator.

Equation (2.8) is shown in Figure 2.8. It has a discontinuity at d = t - [(q)/(T_E)]. Equation (2.8) is monotonically increasing at points not on the discontinuity if the stimulation strength B is less than one. The nature of Equation (2.8) at the discontinuity depends on t. If t < [1/(2(1-q))] the map will be nonmonotonic at the discontinuity as a(f_t ­ d) > a(f_t ¯ d). If t > [1/(2(1-q))] the map will be monotonic at the discontinuity as a(f_t ­ d) < a(f_t ¯ d). In addition, the slope of the map is less than one for a large part of its domain.

2.3.1  The Dynamics of Parasystole

In order to test for banded chaotic behaviour such as that seen in equation 2.6, we collected symbolic sequences as well as sequences of phases from continued iteration of the parasystole map. A symbolic sequence S consists of a list of symbols (s₁,s₂,s₃...) where s represents a discrete region of the map. A symbolic sequence for maps like Equation (2.8) would typically involve three or more symbols to represent the individual branches of the map. However, we collect a symbolic sequence M = m₁,m₂,m₃... where m equals the number of iterates before the trajectory passes to the leftmost branch (Figure 2.8) of the parasystole map. This is equivalent to counting the number of successive sinus beats between each ectopic beat, a method that is frequently used to characterize parasystolic rhythms[39], however here we count both expressed and suppressed beats. The symbolic sequence M also gives information on the rotation number; if M repeats every N symbols, then r = N/(m₁ + m₂ + ... m_N).

The N banded chaotic trajectory would appear the same as a period N trajectory when counted as symbolic sequence M, as long as the trajectory is not dense around the origin. Values of B and q that give rise to banding should then have periodic symbolic sequences with aperiodic sequences of phases. Regions where no banded chaos or periodicity exist should have aperiodic sequences in the phase and symbolic sequences. Regions that are periodic will have both periodic phase and symbolic sequences.

Equation (2.8) was iterated over 20000 times for a large range of values of B and t. Periodicity was determined for both symbolic sequences and phase sequences. The results are summarized in Figure 2.9. In Figure 2.9, black dots were plotted at points in the parameter space where periodicity was detected in the phase (maximum period 1000 for 5000 collected values). Regions shaded in grey correspond to where there is periodicity in M but not in the phase. The period of M found in each of these zones was constant throughout the zone. The relative placement of these zones in respect to the rotation number in the B,t parameter space is the same as that seen for Equation (2.6).

In contrast to relatively simple situation found in equation (2.6), the parasystole map has a slope that is less than one over a range of f_t, which allows for periodic and chaotic dynamics[135]. Here we investigate the relationship between rotation interval, Lyapunov number and periodicity for equation (2.8). Figure 2.10 summarizes the results obtained by determining the Lyapunov number and the rotation number for t from 0 to 1 (400 steps using multiples of an irrational number) for a constant weak stimulation strength (B= 0.5), for a range of initial conditions. Figure 2.10A shows the bifurcation diagram , 2.10B shows the Lyapunov number vs t, 2.10C shows the rotation number vs t, and 2.10D shows the range of rotation numbers obtained for ten initial conditions vs t.

The dynamics can be seen to be periodic from the bifurcation diagram (Figure 2.10A) and the Lyapunov number (Figure 2.10B) for [1/(2(1-q)))] < t < 1 (t > 0.71). Equation 2.8 is monotonic increasing at the discontinuity for t in this range (Sec. 4), and the dynamics are known to be periodic (r is irrational on a set of measure zero) [135] for these maps. Another region of periodicity can be seen between t = 0.45 and 0.55. Bifurcation diagrams identical to this periodic region have been published separately by Courtemanche (from a model of parasystole[41]) and Bressloff and Stark (from a discontinuous circle map model of neural networks[21]). Regions in this diagram corresponding to the chaotic banding zones seen in Figure 2.5 are recognizable by the formation of distinct bands in this figure (ie the 3:1 region can be located at t = 0.3).

The Lyapunov number is positive for a large range of values, indicating chaotic dynamics. For most values of t where the Lyapunov number is positive, the rotation number is sensitive to initial conditions (Figure 2.10D), resulting in a range of rotation numbers for different initial conditions. The rotation number is seen to plateau at certain regions corresponding to a local decrease in the Lyapunov exponent (Figure 2.10C). Local decreases in the Lyapunov exponent also correspond with a reduction in the rotation interval to zero even if the Lyapunov exponent is still positive. When the Lyapunov number is negative, there is periodicity in the phase of Equation (2.8). If the dip in the Lyapunov number remains above 0 there is still a plateau in the rotation number, along with a constant rotation number, however there is no periodicity in the phases of the system.

The relative stability of the rotation number in the regions where there is a dip in the Lyapunov number is due to a narrowing of the range of rotation numbers. The rotation numbers in these regions are constant (one hundred initial conditions ranging from phase = [0.01:0.99] were checked running the model for 800000 iterations). In addition, it is apparent when comparing Figures 2.10B and 2.10D that there is a positive relationship in the size of the rotation interval and magnitude of the Lyapunov number above zero.

The shrinking of the rotation interval can be predicted using the same technique as used for Equation (2.6) The two maps bounding Equation (2.8) were generated and the rotation numbers computed for 1000 points along the t axis. The resulting diagram (Figure 2.10E) completely bounds the rotation number spread seen in the numerical simulations. In addition, local increases and decreases in rotation number spread are closely tracked using this technique. It should be noted that there is an order of magnitude difference between predicted and observed rotation number spread. This difference does not significantly change if more initial values are chosen (data not shown).

As the analytic results shown in Figure 2.10E give the upper bound for the magnitude of the rotation interval, we confirm the numeric results from our estimate of the rotation interval (Figure 2.10D.) It is important to note that there are regions where the rotation interval is zero, but there is still a positive Lyapunov exponent. The downward spikes in the Lyapunov exponent correspond to parameters where the chaotic trajectory becomes bounded (Figure 2.10B), and where the rotation interval drops to zero. The dynamics in these regions, as well as the organization in the parameter space is analogous to that seen in Equation (2.6).

The parasystole map can generate a sequence that correspods to the number of expressed sinus beats between successive expressed ectopic beats. The number of intervening sinus beats (NIBs) have been used to classify parasystolic rhythms both in theoretical models and in clinical situation[39,96,250]. In Figure 2.11 we plot the different NIB values found from iterating Equation (2.8) against t for the same stimulation strength and refractory period used in Figure 2.10. Ectopic beats are concealed for t < 0.08 and t > 0.92. Values of t that result in zero rotation interval width but with positive Lyapunov exponent (the `banding zones') give rise a large number of possible NIB values (91 different NIB values at t = 0.337857) . NIB values of greater than 1000 have been found in these regions (NIB=1313 for t = 0.09946). In between these regions, we find that the number of different NIB values for a given t decreases. The generation of NIB sequences for maps with parameters outside of the `banding zones' appear to follow the rules governing NIB sequences for pure parasystole, where the coupling strength (B) is zero, as characterized by Glass et al.[96]: 1) there are at most 3 NIB values, 2) only one of the three values is odd; 3) if there are three NIB values, the largest value is one less than the sum of the smaller values; and 4) in any sequence with three possible NIB values, only one can follow itself. A generalization that is valid for Equation (2.8) for parameters within the banding zones is that if there are more than two possible NIB values, the any larger NIB values found are equal to the sum of two smaller NIB values plus one. These observations concerning patterns in NIB sequences may be valid only for the parameter range tested.

2.4  Discussion

We have analysed the bifurcations of a simple noninvertible discontinuous circle map, and applied these results to understanding the bifurcations of a general model for modulated parasystole. In Section 2.2.1 of this paper, we found that noninvertible discontinuous maps can display banded chaos, where a chaotic trajectory is restricted to narrow regions in the phase space. The bifurcations in parameter space associated with this behavior share characteristics associated with entrainment zones (Arnold tongues) characterized for the sine map (Equation 2.3)[6,8,129]; there are zones in the parameter space where all initial conditions result in the same rotation number, and these zones are organized such that between any two zones with rotation number N:M and N':M', there is a third zone with rotation number (N+N')/(M+M'). Unlike Arnold tongues, the trajectory for maps within the zones is not periodic. In addition, the zones cannot overlap, in contrast to Arnold tongues which overlap as parameters cross the critical line.

The above behaviour was demonstrated for a family of noninvertible discontinuous maps where the slope is always greater than one (Equation 2.6). The map associated with a model for modulated parasystole has slope less than one over a range of phases, and can be inverible or noninvertible, depending on the parameters used. Despite these differences, the modulated parasystole map displays banded chaotic behaviour, which was numericaly characterized by a positive Lyapunov exponent in the presence of a rotation interval of zero width, over a large range of parameters. The organization of parameters in the parameter space that result in banded chaos is analogous to that seen in equation 2.6.

The relevance of this work to understanding clinical data is not clear, however it is interesting to speculate that a banded chaotic trajectory could be mistaken as a periodic rhythm with noise. It should be noted, however, that the number of observed sinus beats between ectopic beats (NIBs) will not be periodic for a banded chaotic trajectory. A chaotic `band' will always cover regions of the map that correspond to missed vs blocked ectopic beats, as this is the site of the discontinuity. This results in aperiodicity in observed NIBs, although the total number of sinus beats (observed and unobserved) between ectopic beats can be periodic. A banded chaotic map can have observed NIB sequences characterized by long uninterrupted trains of sinus beats, and by a large range of possible NIB values. For parameters outside of the `banding zones' we find that the NIB values appear to follow patterns associated with pure parasytole[96]. A carefull analyses concerning the relation of NIB sequences to banded chaos would help in the comparison of our results to clinical data and the results of other theoretical models.

Chapter 3
Activation Mapping

The present chapter details the methods used for obtaining experimental results in the final chapter of this thesis. The primary experimental protocol used is optical recording of propagation in monolayers of embryonic cardiac cells using calcium sensitive dyes. Section 3.1 outlines the basic theory of optical recording in the context of designing a system for mapping activation patterns on cultured cardiac monolayers. Section 3.2 details the construction of the macroscope. Section 3.3 gives details on the camera used as well as general considerations involved in the choice of the CCD (charged coupled device). Section 3.4 outlines programs used to analyze camera data. Section 3.5 details photodiode array construction. Section 3.5.2 describes design and construction second stage (conditioning) circuitry. Results from optical mapping studies are presented in Chapter 6 and Chapter 7. Relevant programs associated with the above sections are included as appendices.

3.1  Introduction to Optical Recording.

In general, all mapping systems consist of a number of sensitive detecting elements which have their signals amplified, processed, displayed, and stored for later analysis. Monolayers are very thin (1-5 cells thick, depending on plating density) but can span a relatively large area (1.4 cm² in this study). This poses a challenge as signal strength is typically low for thin preparations and the light collecting power of conventional optical devices decreases with increasing field of view. In addition the type of experiments planned required long recording times (10-20 minutes) at high temporal and spatial resolution, putting additional demands on recording and storage apparatus.

3.1.1  Activation Mapping: Electrodes vs. Optical Techniques.

Activity can be mapped using extracellular electrodes or optical techniques. There are advantages and disadvantages in both systems. Extracellular electrodes have been used to map activation in whole heart[111,64] and in cultured tissue samples[120,71]. This system used on cultured preparations involves plating monolayers on electrodes which have been etched on glass using photolithographic techniques. Optical techniques use changes in fluoresced or transmitted light from samples to map activity.

A major advantage to using electrode arrays is that the preparation is not subjected to potentially phototoxic effects of photosensitive dyes. Also, electrode arrays are both sensitive and fast, allowing detection of upstrokes at 0.1 ms temporal resolution[120]. These techniques seem to be the only viable ones for looking at fast activity in thin cell preparations where phototoxicity would affect the results. One drawback is that electrode arrays measure change in membrane voltage as opposed to DC voltage levels. Optical mapping techniques can give direct information on action potential shape.

Extracellular electrode arrays are less flexible than optical techniques because the the spacing and number of detection sites is fixed once the array is constructed. Optical techniques allow the user to `zoom' on interesting areas without effecting the number of spatial points collected. In addition, a constructed electrode array would have to be replaced after a relatively few experiments as electrodes are sensitive to the saline cell bathing media (personal communication with Dr. Israel). Construction of the array involves specialized lithographic equipment, as well as sonic welding equipment not locally available Construction of the array as well as welding to external amplifiers would have to have been done by an external company. Finally, the number of extracellular electrodes is limited by space constraints on the substrate, as present fabrication techniques place electrodes and the wires leading to them on the same plane. The low number of detection sites make extracellular arrays a poor choice for mapping spatially complex activation patterns.

As the appropriate spatial resolution for mapping activation patterns in monolayer cultures was unknown design an extracellular array would have been impractical. We decided to use optical techniques as it provided a more general approach to mapping activity from an uncharacterised system.

3.1.2  Extrinsic vs. Intrinsic Signals.

Optical activity can be recorded by using intrinsic or extrinsic signals. Extrinsic signals are generated when a dye added to the preparation changes its optical characteristics in response to changes in the preparation. Intrinsic signals are activity dependent changes in the optical character of the tissue itself.

Extrinsic signals are generated by dyes which are loaded into the preparation prior to mapping. Hundreds of dyes have been developed which use a variety of mechanisms to respond to dynamic changes in the preparation. Their mechanisms of action will be discussed in more detail in sections below.

Intrinsic signals are generated by three mechanisms: light scattering at the cell level, changes in fluorescence or absorption of intrinsic chromophores, or changes in local blood flow at the anatomic level. The light scattering characteristics[60] of cells are affected by changes in cell volume, rapid changes in the membrane dipoles, or morphological changes associated by mechanisms of secretion. An example of measurement of scattering changes caused by cell shrinkage is demonstrated in studies by Mamose and colleagues[164] and in Kreisman and colleagues[148]. Intrinsic chromophores such as hemoglobin and cytochromes can change their fluorescent and absorptive characteristics during cellular activity. Intrinsic chromophores have been used to map activation in whole brain by Grinvald and associates[106]. Signals from Changes in blood flow and oxygen delivery in in vivo preparations can also change light scattering and reflectance. Examples of these techniques are seen in publications by Grinvald and colleagues[106] and Bonhoeffer and colleagues[20].

An advantage associated with intrinsic signal measurement is that it avoids potentially phototoxic effects of dyes. Dyes can produce phototoxic effects: Di-4-ANEPPS, a popular voltage sensitive dye, is not toxic at low concentrations in the absence of light, but degrades membranes in the presence of even low levels of light[259]. However, the time course of intrinsic signals are relatively slow (1 - 3 seconds depending on wavelength[106]), and often rely on metabolic effects secondary to those being studied. Intrinsic signal measurements are not appropriate for mapping activity in cultures of embryonic cells: The lack of vasculature prevents blood flow measurements, and movement artifacts are not large enough to be detected reliably at low magnification⁴

3.1.3  Fluorescence vs. Absorption Dyes.

Fluorescence and absorption dyes are available to measure extrinsic signals. Absorption dyes change the absolute level of light that reflects off or passes through a sample. Fluorescent dyes produce light when excited by light at a lower wavelength. In general, changes in the property being measured (voltage across the membrane, or concentration of a particular ion) cause a small change in light intensity (DI) that is accompanied by a large offset signal (I). The signal to noise ratio increases as the square root of total signal (see photodiode discussion following, and [105]), and is directly proportional to the fractional change in detected light caused by changes in the measured quantity [105]:

signal to noise ratio µ ÖI (DI) I

(3.1)

There are advantages and disadvantages associated with absorption and fluorescent measurements. Absorption dyes change their absolute level by a very small fraction (0.01 -0.1 %) and fluorescent dyes change by a much greater amount (5% for some voltage sensitive dyes to 100% for calcium sensitive dyes). However as the absolute light level of absorption signals is 3-4 orders of magnitude higher, the signal to noise ratio can be better for these dyes. In addition, absorption dyes differ from fluorescent dyes in respect to their sensitivity to nonspecific binding and sample thickness. Absorption dyes are less sensitive to the amount of nonspecific binding in the sample because background fluorescence (I) increases with nonspecific binding (thereby decreasing the signal to noise ratio), but transmitted light does not. Also, background fluorescence increases with increased sample thickness, but background transmitted light does not. Therefore the signal to noise ratio of a transmission experiment increases linearly with sample thickness (as DI increases), but increases only as the square root of thickness in a fluorescence experiment[105].

The advantages that absorption dyes have in regard to staining thick preparations with nonspecific binding do not apply to monolayer cultures, as they are relatively thin and are homogeneous⁵. In addition, the small percentage change in absorption signals requires a method for subtracting off a large light offset. This can easily be done with specialized circuitry that subtracts off a constant voltage offset, but it is not available in any commercially available camera technology. Absorption dyes therefore are better suited for photodiodes (to be apparatus discussed in detail later) but are not appropriate choices for charge coupled devices (CCDs) where the offset would be stored and digitized along with the signal⁶. Fluorescent dyes offer more flexibility in experimental design by allowing use of conventional camera technology. However, fluorescent signals can be very weak and are often at the noise floor of conventional camera equipment. The low light level necessitates design of specialized optical to capture the maximum amount of light as well as special care in selecting recording equipment (both to discussed in separate sections below.)

3.1.4  Slow vs. Fast Fluorescent Dyes.

Fluorescent dyes can be divided into slow and fast acting dyes. Slow dyes change their spatial distribution in cells with changes in cell potential or pH. As these dyes move through cell membranes in response to activity, the dyes time course is affected by the permeability of the dye, its aggregation characteristics, and its tendency to bind cell constituents. Studies that use such dyes are reviewed by Ebner and Chen[60]. Slow dyes have time constants on the order of 1 - 20 seconds and therefore are not useful for mapping fast activation changes. Fast dyes do not rely on dye movement, but change their optical characteristics by changing their conformation. Conformation changes occur when the dye binds to ion species released during metabolic activity or when the dye changes due to local changes in voltage. The time constant of such dyes are on the order of one millisecond. The mechanism responsible for the optical changes in these dyes have not been fully elucidated[60].

Fast voltage sensitive dyes are divided into four categories based on their chemistry: merocyanine-oxazolone, merocyanine-rhodanine, oxonol, and styryl. There are several mechanisms that account for a given dye's voltage sensitivity. The dye molecule can change its shape (molecular motion) or change its charge distribution (elecrochromism). Changes in the molecule's shape can change the depth of implantation, its orientation within, and its tendency to move on or off the membrane[79]. A variant on molecular motion is the rotation-dimer mechanism[82]. Here the dye switches from a monomeric form oriented perpendicular to the membrane to a dimer orientated parallel to the membrane with changes in voltage. Also, Fromherz and Lambacher[82] have suggested a molecular motion mechanism for styryl dyes where the `twist' of the molecule effects charge movement within the molecule. In contrast, electrochromism refers to the direct displacement of charges within a molecule by external electric fields. As there is no molecular movement involved, electrochromism effects have a nanosecond time scale. Dyes are designed to incorporate one or more of the above effects to optimize fluorescent and absorption characteristics[60,79].

3.1.5  Ion Sensitive Dyes.

Ion sensitive dyes change their conformation by reversibly binding with ion species. Dyes sensitive to sodium, pH, chloride, and calcium have been developed. Ion sensitive dyes can be broadly classified as ratiometric or non ratiometric. Nonratiometric dyes change their absolute level of fluorescence at a given wavelength. Examples of nonratiometric dyes are Fluo-3 used for Ca⁺⁺ and SPQ for Cl^-. Ratiometric dyes are either dual emission, which fluoresce at two different wavelengths, or dual excitation, which are excited by two different wavelengths. Examples of dual emission dyes are Indo-1 (for Ca⁺⁺ and Mg⁺⁺), DCH (for H⁺) and FCRYP-2 (for Na⁺⁺). Common dual excitation dyes are Fura-2 (for Ca⁺⁺, and Mg⁺⁺), SBFI (for Na⁺⁺), and BCECF (for pH).

Non-ratiometric dyes must be calibrated against known stock concentrations and controlled for light level variations between experiments if quantitative measurements are required. Ratiometric dyes give information on the absolute level of ion present without calibration. For example, Fura-2 increases fluorescence intensity in the presence of Ca⁺⁺ when excited at 350 nm, but the intensity decreases when excited at 380 nm. As the individual measurements are equally affected by variations in light levels and dye concentration, the ratio of fluorescence depends only on the Ca⁺⁺ concentration.

Ratiometric measurements use more equipment than non-ratiometric measures, and can be problematic when measuring rapid events. Dual excitation dyes must be stimulated at two wavelengths (at separate times), and sequential frames must be divided by each other for each measurement. Data collection rates are therefore half the acquisition rate of the imaging device. In addition, fast biological events can cause discrepancies between ratioed frames, resulting in blurring artifacts. Dual emission measurements can avoid these problems by splitting the image onto two detectors, however this greatly increases the cost of the measuring system. As we were interested more in measuring activation times, not absolute ion levels, and cost and complexity of the mapping system was a concern, nonratiometric dyes were considered.

3.1.6  Dye Choice.

The above considerations were used to choose two dye candidates for detection of signals in monolayers. The choice of dye also on tissue specific uptake: certain dyes are absorbed more easily by different cell types[60]. Di-4-ANEPPS, a fast voltage sensitive styryl dye and Fluo3, a non ratiometric calcium sensitive dye, were used. Di-4-ANEPPS was chosen as it had been successfully used in single cell and aggregate cultures of chick myocytes by other researchers[260]. Experiments using Fluo3 were initiated in Dr. Nelson Publicovers laboratory, where the correct cell densities and dye concentration were determined.

3.1.7  Detection Systems.

Fluorescent activity can be detected using a variety of different systems. In general, all transform light energy to an electrical signal which is displayed and stored in some fashion. Light detection systems considered are photomultiplier tubes, photodiodes, photocathode tubes, and CCDs. The choice of detector depends on application.

All detection systems are governed by fundamental limitations that govern spatial and temporal resolution in relation to signal to noise ratios. For example, given a constant level of system noise and illumination, increasing spatial resolution decreases signal to noise ratio by decreasing the amount of light hitting the detector. If spatial resolution is increased by dividing up the detector area into more elements, then signal strength decreases as the inverse square of element width [95]. Increasing temporal resolution can decrease signal levels by decreasing integration time (in the case of CCDs, see sections below) or by necessitating an increase in bandwidth which increases noise (as is the case with photodiodes, see below.)

3.1.8  Photomultiplier tubes.

Photomultiplier tubes (PMT's) are single element detectors that are well suited for very low light level applications. The photomultiplier works by converting light energy to electrons via a phosphor screen. The resulting charge is amplified within the tube as electrons are cascaded off charged plates. Through this amplification process the PMT is capable of detecting light levels as low as a single photon. Due to their sensitivity, PMTs are ideal choices when only a single spatial point need be detected. PMTs are too bulky to practically mount more than one or two on an imaging setup. A multiple site mapping method using a single PMT was devised by Dillon and Morad[56,181], where a laser is used to scan the surface of the sample, and the signal detected by PMT is later decoded into individual channels. This system is powerful and flexible, but we decided against it as building a laser scanning device was beyond our technological means. In addition, a primary advantage for laser systems is for scanning irregular or curved surfaces as no light collecting optics are involved: this allows scanning without consideration of the focal plane of the object or the depth of focus of the collecting optics. This is not a consideration when recording off flat monolayers of tissue.

3.1.9  Photodiodes.

Photodiodes are individual photosensing elements constructed from a layer of silicon mounted on a substrate that acts as an electron well. When light hits silicon, it causes the release of electrons, providing current. The current is small, and has to be carefully amplified for low light level applications. Arrays of photodiodes can be constructed (either from placing individual elements in a grid or other pattern, or by purchasing an array of elements with leads coming out the back and sides.) Photodiode arrays have been used for many years to detect activation patterns in a variety of preparations [33,220,222,46,105,107,94,71]. The primary advantage offered by photodiodes is that each element is amplified individually and in parallel so the signal degrades minimally after its been transduced from light to electricity (as opposed to charge coupled devices, discussed below.) Also, as each signal is treated individually, the experimenter can signal process each channel before storing the signal. The main difficulty in implementing a photodiode system is that each diode has amplification circuitry associated with it, necessitating significant space and construction time for even small arrays. This makes construction of large arrays impractical. The electronic back end for such arrays is not commercially available.⁷

3.1.10  Video tubes.

Video tubes are a common type of imaging device used for imaging applications[247]. Video tubes (also called Vidicon tubes) consist of an electron gun and a photoconductive faceplate housed in a cylindrical glass envelope. The photoconductive faceplate consists of a positively charged, transparent metal film. The image is focussed on the faceplate. The intensity of light on the metal film alters the charge carrying capacity of the film: dark areas on the faceplate retain a greater charge than light areas. A beam of electrons from the electron gun scans the faceplate at a certain resolution and frequency.⁸ A current proportional to the intensity of light on the faceplate is produced as electrons from the beam replace electrons in light areas of the faceplate. The resulting analog video signal can be displayed on a monitor or stored on tape. The video signal must be digitized using a frame-grabber board if quantitative measurements are required.

The characteristics of a video tube depend on the metal used as the photoconductor in the faceplate. The photoconductor effects dynamic range, spectral sensitivity and the speed and linearity of response to changes in the image. The dynamic range of video tubes is typically low, but can be increased if the gain (the strength of the electron beam) is varied using an automatic gain circuit. Varying gain allows for adjustment in signal strength between frames (allowing the camera to compensate for a gradual decrease in overall signal intensity, for example) but cannot increase the dynamic range within one frame[238]. The dynamic range within a frame for video cameras is at most 50 fold[245]. This value is low compared to individual photodiodes or CCDs (charge coupled devices, discussed below) which are resolvable to 12 bits (4096 levels) or greater.

Although video tubes have been designed with greater sensitivity and lower lag (by using different photoconductive substrates) the poor dynamic range within each frame is a concern. If large variations in signal levels exist over one frame (due to uneven lighting, patchiness in the monolayer, or heterogeneous dye loading) certain regions within a frame may be unmeasurable even if they are generating adequate signal. CCDs are a better choice if the appropriate dynamic range over a frame is not known.

3.1.11  Charged coupled devices.

Charged coupled devices (CCDs) are arrays of photosensitive elements. CCDs have underlying circuitry allowing the information from each element to be read off the array serially to be amplified and digitized. CCDs will be discussed in detail in Section 3.3.

3.1.12  Detector Choice.

As the proposed experiments demanded both the capacity to have high spatial and temporal resolution, it was decided to develop both photodiode array techniques and purchase a CCD imaging system. Imaging with a CCD offers high spatial resolution but low temporal resolution. A photodiode array offers the opposite. High speed and high temporal resolution CCDs and photodiode arrays exist but were too expensive. The system I will discuss involves using photodiodes and a CCD to provide high temporal and spatial resolution.

Light Collection and Data Acquisition

Fluorescent signals are relatively weak, necessitating the use of special optical apparatus to obtain enough light at low magnifications. Although microscopes are efficient light collectors at high magnification, they can be two orders of magnitude less efficient at low magnifications. Cultured monolayers of embryonic cardiac cells have previously been imaged by Rhor and colleagues[213] and Fast and colleagues[71], however these researchers were interested in small areas (250 mM squared) and could take advantage of high numerical aperture(NA) microscope objectives.

The macroscope used in this thesis work is based on designs outlined by Ratzlaff and associates[211] with some modifications. The macroscope is discussed in detail in Section 3.2.

Large amounts of data can be collected with imaging studies. Previously researchers relied on bulk storage media such as VCR tape or large scale RAM storage devices[260]. We used recordable CD ROM.

3.2  Optical Mapping System Design and Construction

3.2.1  Theory

Microscopes are optimized for collecting bright images from a small field of view by using high magnification lenses mounted close to the sample. Microscopes are not ideal for low magnification applications where long working distances are desirable. In general, the light collection efficiency of microscope objectives increases with increased magnification. For example, a typical 40× lens collects 100 times more light than a 1× objective. In addition, the distance between the lens and sample (working distance) is a few millimeters or less for microscope objectives. Commercial low magnification microscopes (called macroscopes) with long working distances are available, but are typically less efficient light collectors than microscopes [211]. Low magnification, long working distance imagers can be constructed from conventional camera equipment with relative ease. These macroscopes collect light as efficiently as high powered microscopes while allowing for a large (1 cm² ) field of view [211].

In general, image brightness in a conventional microscope imaging application is a function of the size of the field being imaged, the diameter of the lens, and distance between the lens and sample. The constructed macroscope has a similar light path as a microscope, but increases image brightness by increasing the diameter of the lens relative to sample size. In the following I will briefly review the microscope light path and go over some essential concepts related to macroscope construction. Background information on optics can be found in several basic texts[80].

The most familiar way of imaging an object is using a single lens, as shown in Figure 3.1. Here the object (O) sits past the focal point of the lens (L1) and a real image (I) is created on the opposite side. Magnification is given by:

M= Di Do

(3.2)

Where D_o is the distance from the object to L1 and D_i is the distance from the image to L1. Although simple, the above scheme is not used in microscopes as different magnifications require changing object and image distances.

Microscopes typically use a tandem lens scheme, as shown in Figure 3.2. Here two lenses are used. The object is placed at the focal plane of the objective lens (L1). Collected light is projected to infinity and is focussed by an imaging lens, or ocular, (L2) to form an image at its focal plane. The magnification is given by the ratio of focal lengths:

M = F(L1) F(L2)

(3.3)

where F(L1) is the focal length of the objective lens L1 and F(L2) is the focal length of the imaging lens L2. The tandem lens design allows changing magnifications by changing L1 without effecting the location of the magnified image.

In both the single and tandem lens imagers, brightness of the image is a function of the numerical aperture (NA) of L1. The NA is defined as:

NA=I sin(a/2)

(3.4)

where a is the maximum accepted light angle from the object at the optical axis, and I is the refractive index of the medium in which the lens is working (in our case I is always equal 1). NA can be approximated in simple lenses as:

NA= Diameter Focal Length

(3.5)

Light collected from L1 can be limited by the diameter of L2 and the distance between L1 and L2 in the tandem lens design. These points will be addressed below. Image brightness is also a function of the distance between L1 and the object in the single lens design (as that is a variable controlling magnification.) In general, image brightness decreases as the square of the distance to the object in the single lens scheme.

Fluorescence imaging using a single lens scheme would be implemented as shown in Figure 3.3. Monochromatic light is directed at the sample (by use of additional lenses or fiber optic light guides). Fluorescent and reflected light from the sample is collected by the lens and passes through a barrier filter(B) prior to forming an image.

Fluorescence imaging is implemented in the tandem lens scheme by adding a light source, excitation(E), and barrier(B) filters as shown in Figure 3.3. A dichroic mirror (D) acts as a wavelength sensitive reflector. Excitation light passes to the dichroic mirror and is reflected to the sample. Fluorescent light (at a longer wavelength) passes through the dichroic mirror to a barrier filter, and is imaged by L2. The use of illumination in this configuration is referred to as epi-illumination.

Both tandem and single lens imaging schemes have been used effectively for imaging fluorescence at low magnification. Single lens imaging is achieved by using a single macro camera lens, or by using a standard lens with macro adapters (which change the camera/image distance.) Single lenses have been used by Witkowski and colleagues[272] and Rosenbaum and colleagues[218]. Tandem lens designs have been implemented using camera lenses[211] or high NA aspheric lenses[46] to image cardiac and neural preparations. Single lens design is easier to implement as all components are commercially available without using specialized holders and mounts. The tandem lens design, however, has two advantages. Consider a 1X magnification imaging application. The single lens design has the object at 2× the focal length of the objective, while the tandem lens always has the object at the focal plane. The doubling in distance results in a 4× reduction in image brightness. More importantly, the barrier filter(B) may not function optimally as shown in Figure 3.3. Optical filters are designed to filter light when they are perpendicular to the light path. Some optical coatings block different wavelengths at even small (5%) angles. The choice of (B) in the single lens scheme is limited to coloured glass as light is not collimated. The above considerations led to the adoption of a tandem lens scheme.

3.2.2  Design.

The following design considerations are summarized from reference number [211]. Lenses L1, L2 and L3 in Figure 3.3 are standard compound 35mm camera lenses. Two lenses set to focus at infinity are placed face to face so the object and image conjugates are located at the back focal planes of the lenses. There are several advantages when Using camera lenses as opposed to single lenses has several advantages. All lens aberration concerns (spherical aberration, coma, astigmatism, distortion, chromatic aberration) can be ignored as these have been accounted for by the lens manufacturer. Camera lenses come premounted with standard adapters and have focusing sleeves. Camera lenses are high quality and have antireflection coatings on them as a matter of course.

As light is collimated after L1, the distance between L1 and L2 does not affect the location of the image or its magnification. This allows flexibility in mounting the dichroic mirror(D). As shown, the light path suggests that the distance between L1 and L2 is not a factor. However, this distance greatly affects image brightness due to vignetting.

Vignetting refers to light loss at image points off the central axis due to obstructions near the periphery of the light path[211]. Ratzlaff and colleagues explain that the collimated beam is ``somewhat divergent''[211], and that L2 should be have a slightly larger entrance pupil than the exit pupil of L1 to reduce vignetting. This is an oversimplification of the vignetting problem in this design. Part of the design strategy outlined by Ratzlaff and associates[211] suggests that an additional imaging stage for focusing should be added by putting a beam splitter in between L1 and L2. Ignoring losses due to reflection, this results in a reduction of image brightness at L2 for two reasons. First, a certain amount of light is diverted from L2 by the beam splitter. Second, adding a beam splitter increases the distance between L1 and L2 a minimum of a(2)× beam width, which results in light loss due to vignetting. Figure 3.4 demonstrates the effect of vignetting on the tandem lens imager. In order to prevent excessive vignetting (here defined as light loss greater than 50 percent for off-axis points), the image lens (L2) radius L2_r must increase proportionally to the distance between the lenses,

L2r µ d x f

(3.6)

where x is the diameter of the object, f is the focal length of the objective lens, and d is the distance between L1 and L2. If the radius of the image lens is not increased, the brightness of the image falls off as a function of d.

3.2.3  Macroscope Construction.

The macroscope was constructed by mounting lenses and filters on Spindler and Hoyer MacroBench (Spindler and Hoyer GmBH, 37070 Cottingen, Germany) optical components. The macroscope itself (Figure 3.5) was placed on horizontal rails, which allowed it to be slid over different preparations. The macroscope is mounted over an inverted microscope , allowing simultaneous monitoring of the sample from the top and bottom sides at different magnifications. Alternatively, it can be placed over a moving stage and heated bath for monitoring larger preparations. The illuminator can also be moved in the horizontal plane. The frame of the macroscope can be lowered by adjusting the height of specially machined legs. This allows for gross focusing of the macroscope. Fine focus is achieved by manipulating the focusing sleeve on the objective lens.

3.2.4  Optical Components.

The objective lens is a Nikon 80 mm F 2.0 lens. The imaging lens is a Sun Cosmicar videography lens with zoom. The variable focal length of the imaging lens allowed convenient change in magnification and field of view without changing the objective lens. Light from the epi-illuminator was collimated using an Oriel one Lens collimator and focusing sleeve (Oriel, Stratford,CT,USA), and focused using a Canon 50 mm lens. Filters and dichroic mirrors were purchased from Omega Optical (Omega Optical, Battleboro, VT, USA). Filters for monitoring Di-4-ANEPPS were a 1 inch emission filter centered at 530±20 nm, dichroic mirror at 560 nm and bandpass barrier filter set at 630±20 nM. Calcium sensitive Fluo-3 measurements were performed using 490±20 emission filter 515 nm dichroic and a 530 nm high pass filter. Dichroic mirrors and barrier filters were custom made to order, as the size required is greater than those used in conventional microscopes. Standard 1 inch filters were used for generating monochromatic light from the the light source. Additional filters for Di-4-ANEPPS and Fluo-3 were obtained as stock components for the Zeiss microscope (Axiovert 135m, Carl Ziess Inc, Thornwood, NY, USA).

3.2.5  Mechanical Components.

A list of standard Spindler and Hoyer macrobench components used is included as an appendix. Additional parts were machined at the McGill Physics mechanical engineering lab. Special parts were made to interface the camera lenses to the Spindler and Hoyer components. These were simple plates with a threaded hole in the center designed for use with different camera lenses. Nikon lenses were mounted using a bayonet mount taken from a camera body and attached directly to an adapter plate.

3.2.6  Illumination.

A 250 Watt lamp, 24 Volt housing was ordered from Dr Jalife's laboratory located in SUNY (State University of New York) research center in Syracuse NY. According to researchers at SUNY, conventional lamp housings have too much interior space, which allow convective air currents to be created as the housing heats up. We elected to follow their advice and purchased a quiet lamp housing of their design.

A tungsten halogen lamp is used for illumination. Different light sources (Xenon, Mercury arc, and laser) are relatively unstable and can add noise to measurements. Other labs have used different sources, however specific noise reduction circuitry was required[32] in some cases. Tungsten halogen is lower noise (in general, the output of the bulb fluctuates on the order of 10^-5[32], while arc lamps typically have fluctuations 20× as large.)

The majority of energy produced by tungsten lamps is in the UV range. This requires high power lamps for moderate illumination at the desired wavelengths. In addition, the lamp housing can become very hot, and stray UV light can be dangerous if care is not taken. The lamp housing was powered by two Eltech 12 V power supplies. Power supply stability is controlled by running a voltage monitor to the lamp leads within the housing.

Illumination intensity was measured to be roughly 100mW/cm² for the calcium and voltage imaging experiments. A light meter was not available in our laboritory. Estimation of illumination intensity was was carried out using a MRD500 photodiode coupled to a AD345 operational amplifier in the current to voltage configuration (see Section for details on this type of circuit.) The MRD500 photodiode is rated at 6.6 mA/W, or 3.3 mA/W/cm². The current output of the illuminated photodiode was determined by dividing the voltage of the circuit by the feedback resistance. Illumination intensity was determined by dividing the current by the quantum efficiency (see Section 3.3 for definition) of a typical silicon diode at the illumination wavelengths used[80], and converting the current to Watts based on the conversion factors above.

3.2.7  Performance

The light collecting capacity of the macroscope was estimated by using a scientific grade CCD camera (Princeton Instruments EEV 576×384 frame transfer CCD, see Section 3.3 for details.) The camera was focussed on an image using a Zeiss microscope (2.5× objective, NA 0.075). An identical image was obtained using the macroscope, where the macroscopes optics were adjusted so that the field of view was identical. The light collection capacity of the microscope and macroscope were compared by changing the integration time of the camera. The camera had to integrate the microscopes image 20 times longer than the macroscope image to achieve equivalent intensity. I therefore estimate that the macroscope collects 20 times more light than the microscope for similar fields of view.

3.3  Camera Theory.

Commercial video rate CCD cameras have been used to acquire fluorescence data from cardiac preparations. These CCDs are constructed to output analog video rate signals (charge is not digitized after it is shifted off the array.) Davidenko and colleagues[46,47] at SUNY used a Cohu camera to acquire images from a Di-4-ANEPPS stained epicardial sheet at 60 frames per second (fps). Such video rate cameras are usable for some applications, however the experiments performed were at the limit of the detection capabilities of this technology. Researchers at SUNY generously gave our group an opportunity to use their system on a embryonic chick monolayer preparation. Monolayers were stained with voltage sensitive dye Di-4-ANEPPS as described in Section 4.2. Although the camera is mounted on a high light throughput macroscope (see Section 3.2) no usable signal was detected despite repeated attempts. This negative result is not surprising considering the relative thickness of the monolayer in comparison with the epicardial sheet successfully imaged by Jalife's laboratory. Video out cameras in general have a dynamic range of only 8 bits (due to noise generated by the camera) and have low light collection efficiency (due to dead space between light collecting elements.) The majority of noise in the cameras output signal is generated by the process of converting data stored on large numbers of storage elements (on the order of 700×400) to an analog video signal at high rates. Data collected from conventional video can be averaged to reduce noise, however this results in a loss of temporal resolution.

Scientific grade CCDs improve on the design of conventional video chips in several ways. CCDs are programmable, allowing the user to use several processing techniques on chip to significantly lower noise. CCDs are designed to have good light collection characteristics and have a small dead space between light sensing elements. Individual light sensing elements are frequently relatively large, which improves dynamic range. A negative aspect of this technology is that users can no longer use conventional storage and viewing apparatus. Conventional video signals can be stored on VCR tape and viewed real-time on monitors. Storage and display of scientific CCD information is frequently done on a frame by frame basis, which complicates long term recording and analysis of signals. In addition CCD imaging systems differ greatly in acquisition speed, spatial resolution, noise levels, sensitivity, data storage and display characteristics. An understanding of how a CCD functions is essential for choosing a system for demanding applications. The following overview of scientific CCD characteristics is relevant to recording low light signals at high rates.

A CCD consists of an array of light gathering elements with support circuitry. Each element in the array is called a pixel. Light energy is transduced to electrical energy with a certain quantum efficiency, which varies for different CCDs. The well capacity is a measure of how much charge each pixel can hold before saturating. Electrical energy is transferred between pixels with a certain loss determined by the charge transfer efficiency of the CCD. Basic CCDs collect images in the following way: light hits a CCD, which collects charge in individual pixels. After a certain amount of time (the integration time) the charge is shifted, pixel by pixel, to horizontal shift registers. After the charge is collected in the shift registers, it is collected at an output node, preamplified, and eventually digitized. The output node is a capacitor where the charge is stored prior to amplification. This process is repeated for every frame collected.

The above sequence is associated with several noise sources. Shot noise is a random noise due to the particle nature of light. It is proportional to the square root of the number of electrons gathered. As the signal increases proportionately to charge, the signal to noise ratio increases as the square root of charge. Dark current is a source of noise that is caused by thermal interactions at each pixel, and is a function of the temperature of the chip, integration time and the CCD architecture. Frame read out noise is generated as accumulated charge in each pixel is shifted to the output node for digitization. It is proportional to the square root of the pixel readout rate. Electron shot noise and dark noise are generated at each pixel, and are generic properties of any light detector. Frame read out noise is created in the process of moving information from a large number of pixels to a single amplifier and digitiser, and is specific to CCDs.

The signal to noise ratio can be increased by functional and structural changes to the CCD. Structural changes can include cooling, use of multi-phased pinned (MPP) architecture, or back lighting the CCD. Cooling reduces read noise and dark current at the expense of charge transfer efficiency (the amount of charge transferred between pixels during readout.) MPP pixel architecture, reduces dark current at the expense of the light collecting capacity of each pixel in the CCD. Back lit CCDs decrease noise by increasing the quantum efficiency (the percentage of photons falling on the pixel surface that generate electrons.) Normally, some of the light hitting each pixel is absorbed by circuitry associated with each pixel. The back lit CCD architecture avoids this problem but at considerable financial expense. Functional changes involve on chip binning and sub-array read out (to be discussed below) which reduce noise by lowering read out rate at the expense of spatial information.

In general, faster readout rates contribute to noise in the following ways. Higher read out rates result in lower integration times, which reduce the signal levels in each pixel. This increases the relative contribution of shot noise (decreasing signal to noise ratio.) Dark currents contribution is actually lowered at high read out rates, as the amount of current collected is time dependent. High readout rates also increase read out noise, which originates at the processing end of the CCD. Importantly, both spatial and temporal resolution of the CCD affects read out noise levels. Higher spatial resolution require more pixels to be digitized per frame. As read noise is a function of the readout rate of individual pixels, increasing spatial resolution increases read noise at a fixed frame rate.

Another structural feature that must be considered is well capacity. The well size of an electronic component controls the number of electrons it can store before it saturates. Charge is accumulated at various points in the CCD prior to preamplification. The well size of these critical points in the CCD puts an upper limit on the dynamic range of the system. There are three areas of a CCD which are affected by well capacity. The pixel size is proportional to well capacity; large pixels can store more charge. The shift register capacity must be proportional to the well capacity of individual pixels. The output node of the CCD must be large enough to store the charge. Shift registers and output nodes typically hold more charge than a single pixel, but are often not much larger. As pixel binning increases the amount of charge moved to the shift registers and output node, intelligent CCD selection requires knowledge of these parameters.

Binning is a process in which adjacent pixels are summed. Binning can increase the dynamic range and reduce noise of a CCD at the expense of spatial resolution. Binned pixels are dumped to the shift register. The well capacity of the shift register is typically only 4 times that of a pixel. Therefore the dynamic range for certain CCDs is not increased for greater than 2×2 pixel binning. The greatest advantage of binning is that increased frame rates are possible (at the same read out rates) at the expense of spatial resolution. Also, the read noise for a binned vs. non binned image (at the same frame rate) is lower for the binned image (as there are less pixels to process.) Binning also collects more light per (binned) pixel, however this effect can be achieved by summing pixels offline. Offline summation is preferable to on chip binning as it does not permanently discard spatial information.

Sub-array processing is a process in which only a small set of pixels are used. This increases frame rate without increasing the read noise but it lowers the size of the field of view. It should be noted that because there is a certain amount of time and noise associated with moving pixels to the shift registers, a 64×64 pixel array will have lower noise and faster readout than a 64×64 sub-array on a 512×512 CCD.

Shifting information across the surface of the CCD to the output node takes a certain amount of time. Pixels are still light sensitive as charges are being shifted to the output node. Charge accumulated during the shifting process causes blurring, which becomes more noticeable as integration time is reduced as the charge collected during shifting increases relative to the measured signal. Three strategies can be used to reduce blurring. The most common method is using a shutter to cut light as the image is being shifted. This is not ideal as information is lost while the shutter is closed, and shutter movement is a potential source of noise. Another method is using line transfer CCDs. In these CCDs, every second pixel column is covered by opaque material. Charge can be transferred very quickly to these columns as charge moves only one pixel as opposed to being shifted across the whole array. The charge is then transferred without degradation to the output node. The main disadvantage with this system is that each CCD consists of 50% dead space, making low light measurements difficult. The final method, called frame transfer, shifts the entire image to a masked section of the CCD. The image is then read out as the next frame is collected. Information is shifted relatively quickly to the masked region as it is shifted row by row not pixel by pixel. The speed of the charge shift operation reduces blurring significantly. Shifting the image to the output node for amplification and digitization can occur over a relatively long period of time (the integration time of the CCD), reducing read noise. Even though half the array is used for storage, there is no dead space on the section of the CCD used for imaging, making this architecture ideal for low light, high frame rate applications.

Finally, commercially available CCD solutions come with appropriate control software that dictates display and storage of data while the information is being gathered. As a large amount of information can be collected by a CCD in a short period of time (512x512 CCD digitized at 12 bits comes to over 400Kb per frame) display and storage can be rate limiting steps in an experiment.

Scientific grade CCDs do not output video signals but transfer digitized frames to computers for storage and display. Camera manufacturers offer different storage options. The simplest, if least effective, storage technique is to move the information to a frame grabber board while the camera is collecting data. Information is later moved off the board to computer memory for offline display and storage. The main disadvantage of this system is that recording time is severely limited by the amount of RAM on the frame grabber board. A preferable method is moving the information directly to computer memory. This allows online viewing of the data and more rapid transfer to permanent storage. This is also a less expensive option, as frame grabber memory is more expensive than conventional computer RAM. Finally, the data can also be streamed to hard disk. Modern hard disk systems can accept 20-30 megabytes of data per second, which easily accommodates demanding camera applications. However, no camera manufacturer has implemented data streaming in this way. All cameras offered in a reasonable price range are designed to collect a few seconds worth of data at a time. Most are specifically designed for long integration time applications, where the CCD collects an image from a low light source for several seconds prior to digitization. Software as well as CCD choices reflect this bias: CCDs often have low dark current but poor well size or read noise characteristics.

3.3.1  Choosing an appropriate system:

A large number of CCD chips are available and many are ideal for high speed imaging applications. However, very few companies offer cameras which are optimized for low light, high speed recording. For example, Dalsa (Dalsa Inc. Waterloo, Ontario, Canada) offers a range of low resolution high speed camera systems. These cameras come with no software or support hardware, and are not cooled for low noise applications. Dalsa (personal communication) warned that these cameras are 8 bit resolution devices as sold. Very recently, Frank Witowski and colleagues[272] modified a Dalsa 64×64 array camera by adding a cooling stage and recorded potentials from intact heart preparations with millisecond resolution. The modifications involved were beyond the capabilities of our laboratory. Commercial options designed for optical mapping are available but at a prohibitively high price. Fuji Electric (Japan), makes a cooled 128×128 800 fps camera system that outputs differential measurements allowing offset subtraction prior to digitization. Kodak has intensified CCD systems that record more than 1000 fps at 128×128. Both systems cost in excess of a $100,000. Furthermore, these systems output digitized frames to special memory buffers, which limits continuous recording time to periods of a few seconds. Microphotonics Inc. offers a cooled 64×64 CCD array at low cost (about $40,000) however the camera operated with a frame grabber board, limiting recording time to 4000 frames. The company offered a S-bus model for continuous streaming. However this requires the use of a Sun workstation for downloading.

Additional options were not considered due to cost. Various thin film coatings increase quantum efficiency of CCDs. Uncoated CCDs have peak sensitivity at around 800nm, and are light sensitive into the UV range. Thin coatings with various substances can shift peak responsivity to lower wavelengths, improving detection of some dyes. Also, thinned, back lit CCDs show a 70% increase in quantum efficiency over front lit CCDs. Finally, CCDs with several readout nodes are available. These decrease read noise by decreasing the bandwidth at each node. However, cameras have to be designed with additional amplification circuitry, and this increases cost.

Most systems offered at a moderate price (lower than $40000) are general purpose devices that are optimized primarily for long integration time applications. Long integration time devices are useful for observing low light level phenomena that remains unchanged over long periods of time. These devices use CCDs with high spatial resolution (frames are read out very few seconds, reducing read noise and allowing higher resolution, see Section 3.5) and low dark current. Dark current optimized systems use an architecture that limits dynamic range (see Section 3.3), and are therefore not desirable. Finally, manufacturers sell general purpose software that is not optimized for high speed, continuous recording.

Several CCD systems were evaluated, and three camera systems were tested with monolayers loaded with Ca⁺⁺ sensitive dye Fluo3. None of the systems available at reasonable cost were optimized to continuously stream information to a storage device for long periods. In all cases, the rate limiting step in the storage and display process was the software as opposed to camera or computer hardware. In particular, camera control software is designed to be used with inappropriate operating systems; all cameras in a moderate price range operate under Microsoft Windows. Aside from placing unnecessary overhead on computer resources, the operating system caches memory prior to streaming to disk. As a result, the operating system tries to transfer cached RAM to disk in one operation. As the camera continues to download information to the computers direct memory access buffers, overload results, and the camera stops recording.

The configuration chosen was the Princeton Instruments (Trenton, NJ, USA) camera with the EEV 576×384 frame transfer CCD. The choice involved software considerations (it was one of two cameras tested that could stream data to disk as described above) and CCD options. Princeton Camera offered the EEV 576×384 frame transfer CCD which has several advantages over those offered by other companies. The CCD has larger than average pixels, and doesn't use AGP architecture (a common CCD configuration which limits dark current at the expense of well size.) Noise associated with frame transfer is comparable to more modern CCDs. The camera is cooled, lowering shot noise and dark current. The camera transfers information directly to computer memory via an interface board mounted on an EISA bus (an extended computer bus architecture designed for fast throughput). Each pixel is digitized at 12 bits, and the image is displayed in a window on the computer screen. The software (Winview, Princeton Instruments, Trenton, NJ, USA) has a mode of operation in which screen display is updated only once all essential storage and camera operation are complete. This allows continuous storage with no skipped frames while still permitting intermittent screen display for visual feedback during experimental runs. Images are stored directly on hard drive (Micropolis model 4AV) and then permanently stored on recordable CD ROM.

The camera can transfer a maximum of 1.4 megabytes per second to the computer due to limitations in the serial cable connecting the camera to the control board. This limits the camera to 8 frames per second at full resolution. By binning the image extensively, high temporal resolution was achieved (200 fps, using 5×5 binning on a 100×100 subarray.) Even at 8 frames per second, the software cannot store and display this information for more than a few frames without crashing due to direct memory access (DMA) overruns as described above. The operating systems interaction with the program limits the bandwidth of the camera system to 75 Kb/s. However, as experiments at low light required extensive binning, 75 Kb/s proved to be sufficient for some applications. Typical experiments used 3×3 binning from a sub array on the chip for a final spatial resolution of 50×50 pixels. Frame read out rate was set at 20 fps to limit throughput to 75 Kb/s. Although temporal resolution is low, continuous records of up to 30000 frames (the upper limit allowed by the software) were collected. Records at different temporal and spatial resolution were acquired within the described throughput limits.

3.4  Data Analysis.

The software sold with the Princeton Instruments camera was not capable of playing back the data files collected during experimental runs. The camera control software (Winview, Princeton Instruments) saves data from the camera in a proprietary format. Princeton Instruments provided details of the SPE format. Custom software was written that allowed playback, image processing, and data analysis.

The software (called gview.exe) runs under a 256 color DOS window in Windows 95. The interface allows mouse control and keyboard control of a variety of functions. In the following section summarizes image processing and data analysis methods in gview.

3.4.1  Pixel and Frame Transformations.

Each pixel undergoes the following transformation prior to display:

Final Value[x,y] = Offset + (Pixel[x,y]) Stretch

(3.7)

where Offset and Pixel Value are global variables (variables applied to each pixel regardless of spatial location) that can be adjusted by the user while the program runs. Offset functions as a brightness control, and Stretch controls contrast. The above transformation would be adequate if each pixel requires the same brightness and contrast adjustment, but this is rarely the case. For example one pixel can be changing its value between 1000 and 2000 in response to changes in the preparation, while another can be changing between 300 and 305. A single choice of Offset and Stretch would display appropriate information for only one of the pixels.

Pixel specific transformations allow automatic setting of offset and and stretch parameters for each pixel over the whole frame. A background frame and a pixel normalization frame can be calculated and applied to the raw data prior to display. The user can determine which frame or frames are used to calculate background or normalization frames and update them as required while viewing the data. Background frames are obtained by averaging each pixel in 20 successive frames and storing the result as an array of offset values. A normalization value for a pixel is determined by finding the maximum and minimum data values that the pixel must display over a certain time range (set at a default of 50 frames). The range (maximum - minimum value) for each pixel is stored in a array, which is then used to scale the data for each pixel in the display frame. Background subtract is a common feature found in most image processing packages, however pixel normalization is specific to this application.

In addition to the above transformations, adjacent pixel averaging was implemented. Data values from 8 nearest neighbours are averaged and used as instead of the raw pixel data. When this is done the background frame has to be spatially averaged in the same way to ensure that the spatial frequency of the background and displayed frame are the same. Spatial averaging is a common transformation used to reduce noise. Temporal averaging, where pixels are averaged over a number of subsequent frames, is implemented when calculating background frames. Temporal averaging is not applied on raw data continuously as this transformation tended to blur travelling wavefronts.

The final value at each pixel given the above transformations (background subtract, pixel normalization, and spatial averaging) is:

Final Value[x,y] = Offset-Background[x,y] + I å n=-I  J å n=-J  pixel[x+I,y+J] n (Stretch ×normalization[x+I,y+J]) .

(3.8)

where I=J is 1 if eight adjacent pixels are averaged.

Offset and Stretch parameters still operate as brightness and contrast that operate over the whole image.

Examples using the transformations in Equation (3.8) are shown in Figure 3.6. Data was obtained from a spontaneously beating monolayer preparation (see section 4.2). Figure 3.6A shows three frames of data with minimal processing, as only the Offset and Stretch are adjusted. Figure 3.6B shows the same frames with background subtraction. Figure 3.6C shows the effect of averaging, and Figure 3.6D shows the effect of pixel normalization. Details in the data become progressively more clear; pixel normalization brings out details hidden by uneven lighting in the lower half of the frame, as well as low dye loading and response in the top left corner.

3.4.2  Data Analysis within Gview.exe.

Relative timing of events is important when mapping activation patterns. Two techniques for determining the timing of events at different points in the frame were developed. The first involves selection of pixels of interest within the frame using the cursor or the mouse. Pixel intensity within the selected area is summed over the area, and the result plotted in a separate window and saved to disk as a separate file. The second technique involves detecting wavefront location automatically and generating isochronal maps (see [118]) from the transformed data.

Plotting pixel intensity over time gives a measure of the fluorescent signal at one point in the frame over time. As data is noisy, it is useful to average a number of adjacent pixels to reconstruct the underlying fluorescent signal. The user can select any number of rectangular areas with the mouse or cursor and the resulting intensity vs time plots are displayed. The location of the boxes can be adjusted, and additional regions can be selected and deselected as needed. The resulting intensity vs frame number plots are stored in separate files with a header giving frame number and the dimensions and location of the bounding rectangle for later analysis. The summed pixels values stored in the files are summed raw data and are not affected by the transformations in the above equations.

Frequently, useful timing information can be obtained from within gview by manipulating the intensity plots on the screen as they are generated. The user can move and scale traces within the display window. This allows traces to be lined up over each other in order to better resolve timing between events. This technique has been useful for rapid analysis between experimental runs enabling change of protocols if needed.

Figure 3.7 shows traces generated from the spiral shown in Figure 3.6. The location of the rectangles used is shown in Figure 3.7A. Spikes can be observed as the wave travels through individual boxes. This diagram indicates that the wave front travels in a circular path.

Data can also be exported from gview in the form of isochronal maps. Isochronal maps are contour maps where each contour gives the location of the wave front at constant time intervals. Isochronal maps have been used for displaying mapping data from optical[47] and electrode mapping[194] techniques. Gview can generate isochronal map files from image data directly ⁹.

Maps are generated as follows. Gview is operated in 2 colour mode, and the user manipulates parameters in Equation (3.8) to remove noise and enhance the signal. The user has direct access to the background and normalization frames in this process. Noisy areas can be selected with a mouse and blacked out, or especially faint regions can be enhanced. These manipulations are sometimes necessary as compression of data from 12 to 1 bit can result in artifacts. Also, regions of the frame that are not of interest can be removed from the final calculation using this technique. Parameters in Equation (3.8) cannot be changed once the isochron mapping routine is initiated, as this could generate false data.

Isochron generation starts once the user has a reasonably clear travelling wavefront on screen. The algorithm used removes high frequency temporal events using a heuristic approach without resorting to temporal averaging. Temporal averaging tended to blur the location of the wavefront. The user can change parameters that control how algorithm calculates the location of the wave front. For example, noise in the signal can create small transient `holes' in the recorded traveling wave. Holes in the wave can be filled depending on the size of the hole, the number of frames that the hole persists, and the intensity of the surrounding wave. Isolated large events (a region becomes active but persists for one or two frames only) can be also ignored. As the algorithm assumes the presence of smoothly travelling waves and tries to ignore isolated events, care must be taken to make sure that the resulting isochrons represents actual data.

The location of the wave is stored in an array which stacks the wave location over a user controlled number of frames prior to outputting the result to a file. Many sequential isochron files can be generated; for example a five second sequence consists of a 100 frames of raw data at 20 fps. Outputting an isochron file every second generates 5 isochron files, each containing the location of the wave over 20 frames.

A separate program (isoview.exe, see appendix) is used to view the isochron files. It reads in consecutive isochrone map files, and displays them as a movie by hi-lighting the location of the wave over time. The user can then select a frame from this movie and output it as a portable greymap file. This allows the user to choose the time window over which the final isochron map is generated. An example of an isochronal map generated from the same data set shown in Figure 3.6 is given in Figure 3.8.

3.4.3  Additional Output Modes.

Gview also can output files as a frame stack. This output mode is useful for compressing large amounts of data to a single file. Pixel values are summed over each column in a frame and the resulting row of summed values are outputted instead of the whole frame. An example of a spiral wave displayed as a frame stack is shown in Figure 3.9.

Gview outputs all graphics as portable greymap ¹⁰ files Portable greymap (PGM) is an open graphics standard which can be read and converted by several commercially available graphics programs.

3.5  Photodiode Array Construction.

This section discusses construction of a photodiode array with second stage amplification. Photodiodes were used in pilot studies of monolayer preparations stained with voltage and calcium sensitive dyes (see Chapter 6). Photodiode measurements were largely abandoned as the camera system (Section 3.3) was better suited for high spatial resolution mapping of the monolayer preparation. Photodiodes will be used in future studies where millisecond temporal resolution is required.

Equations governing noise and gain of photodiodes and operational amplifiers can be found in basic textbooks[23,117] and technical publications[22]. The photodiode array constructed in our lab is similar to those used by several researchers[32,60,69,107,220,222,259]. Some distinguishing features of the photodiode array used here are very high value feedback resistors and individual photodiode elements. The second stage unit introduces a novel circuit design to implement computer controlled gain and offset adjustments. Design and construction of the photodiode and second stage are discussed in Sections 3.5.1 and 3.5.2 below.

3.5.1  The Photodiode Array.

signal to noise ratio = RMS voltage of the signal RMS voltage of the noise

(3.9)

where RMS voltage is the root mean square (the square root of a mean of a set of squared values) of the voltage¹¹ The signal increases in direct proportion to R, so signal to noise ratio increases as the square root of R. For low light signals, R_gain is typically in the gigaohm range[22].

Noise from the operational amplifier that arises from the shot noise of the input bias current also influences the output, and is amplified the same amount as the signal current. The noise shows a characteristic frequency dependence known as gain peaking. Under DC conditions, the noise is amplified as 1+ R_gain/R_photodiode[22]. Although R_gain is in the gigaohm range, R_photodiode is larger, minimizing noise gain. However, the capacitance of the photodiode passes noise current at higher frequencies which effectively decreasing R_photodiode, adding significant noise. This noise can result in destabilizing oscillations as its amplified by high R_gain. The amplification of gain peaking noise is given by 1+ C_photodiode / C_gain[22], where C_photodiode and C_gain are capacitance of the photodiode and feedback resistor respectively. Gain peaking noise is attenuated by stray capacitances across the feedback resistor. As C_photodiode is usually in the hundreds of picofarads, and C_gain is low, noise is amplified by a factor of several hundred. To minimize this effect, low input bias operational amplifiers, and low capacitance, high resistance photodiodes are required[22]. Also, as gain peaking occurs at higher frequencies, it can be reduced by putting a capacitor across the feedback resistor. As picofarad capacitances can decrease bandwidth below 10 Hz if R_gain is large, this approach is limited by bandwidth requirements.

Another source of error is a DC offset current, created from high R_gain as it develops significant thermal DC voltage drift from the small current from the inverting input[22]. If the offset error is large, it can saturate the operational amplifier. To compensate for this error an equal resistance can be put across the operational amplifier noninverting input and ground. The remaining DC error is caused by mismatches between the two resistors and the amplifier input currents at the two inputs. However this approach causes a voltage drop across the photodiode as the noninverting input is no longer held at ground. This results in a leakage current from the photodiode, which can be as large as the original offset. In addition, matched gigaohm resistors are costly. The original problem is best avoided by using an operational amplifier with low input current specifications.

Finally, a major source of noise is interference from external noise sources. A current to voltage converter is sensitive to noise coupling from electrostatic, magnetic and radio frequency sources. Electrostatic coupling supplies noise signals through the mutual capacitances that exist between any two objects[22]. To minimize this error signal, shielding is used to intercept that current and shunt it to ground. A potential problem with shielding is that it can create parasitic capacitances which can add to gain peaking noise. Electrostatic noise can come from sources independent of the circuit, such as power lines, but can also come from neighbouring operational amplifier outputs such as operational amplifiers placed in an array. External noise sources are easily identified and shielding can be added until the noise disappears. However signal pick up from neighbouring amplifiers are more difficult to isolate. Amplifiers are in close proximity to each other as they must be mounted as close as possible to the photodiodes they monitor (long wires leading from the photodiodes would pick up external interference, and shielding these wires to decrease noise would add capacitance to the signal, increasing the gain peaking effects.) As a result, the output of one amplifier may be close enough to the input of others to be detected and amplified. Special attention must be made to avoid this possibility by keeping the outputs of the operational amplifiers as far as possible from the inputs of other amplifiers. One approach to minimizing this problem is using fiberoptic cables to carry the signal from different points in the image to the photodiodes. Photodiodes then can be placed farther apart. The sensitivity of the current to voltage circuit to noise makes careful wiring and soldering very important. It is necessary to keep the circuit as simple as possible around the amplifier inputs to avoid unnecessary interference and capacitance effects.

Design of a photodiode array involves defining the required performance of individual sensing elements, as well as taking into account noise sources outlined above. In general, increasing gain of a circuit decreases the frequency response, or bandwidth, of the circuit. Therefore, one has to consider the circuits expected input levels and required output levels as well as bandwidth requirements.

The circuit should be as sensitive to light as possible, so that excitation light levels can be kept low. If the circuit saturates for a given input level, excitation light can be lowered with the added benefit of lowered phototoxicity. Therefore sensitivity should be maximized given bandwidth constraints.

Resolving the upstroke of the action potential requires millisecond temporal resolution, which at first suggests a circuit with kilohertz bandwidth. Such a high bandwidth is not required, as action potentials have durations of hundreds of milliseconds. Action potential shape is distorted in a predictable way by low frequency filtering, and can be reconstructed if necessary. It is also possible to determine with millisecond resolution the action potential start time: the circuit responds at the start of a stimulus even though the amplitude of the response is attenuated. As cardiac action potentials are 100 ms in duration or more, 10 Hz bandwidth can be sufficient. As a decreased minimum bandwidth allows use of higher gain, circuit design should take advantage of this expected signal shape in order to optimize performance. It should be noted that neural action potentials have durations of only a few milliseconds. Detection of neural action potentials would require circuits with higher bandwidth.

The Motorola MRD500 photodiode was chosen as the detector element in the array because of its low capacitance[278]. This minimizes the effects of gain peaking. Its small photosensitive area (the capacitance of a photodiode is proportional to its size) is compensated for by a collecting lens built onto the photodiode. The resistance of the photodiode is large enough to minimize noise for low frequency signals. The wavelength is typical of most silicon photodiodes, and peaks at 800 nm. Although its possible to get diodes with sensitivities in other ranges, this frequency is close to the ideal range needed for Di-4-ANEPPS (630 nm) measurements. The feedback resistor chosen was a 10 gigaohm Eltec glass resistor. The resistor is very small, allowing it to be mounted directly across the inverting input and operational amplifier output, minimizing interference noise caused by longer leads. In addition it is expected to have small stray capacitance due to its size, preserving bandwidth.

The amplifier used is an Analog Devices AD645. It was chosen because it has very low offset voltage[278], which allows shunting the noninverting input directly to ground. Offset correction circuitry would have added an additional resistor and bypass capacitor to the inputs, which would increase the circuits tendency to pick up noise. The operational amplifier also has low input bias current, which minimizes the effects gain peaking.

Amplifiers are glued on their backs to bread-board. This allows attachment of the feedback resistors across the underside of the operational amplifiers. This arrangement minimizes the amount of wire attached to the input pin. The inverting input is attached to a photodiode, and the noninverting input is attached directly to a groundplane surrounding the breadboard. The photodiode is then attached to the groundplane in close proximity to the noninverting input. The output of each amplifier is as far away from the photodiode as its package allows. The output signal moves off the board and up through a shielded printer cable for further processing. This configuration allows 8 operational amplifiers to be mounted (on either side of the board) in close proximity to the photodiodes while minimizing electrostatic coupling effects from the operational amplifier inputs and outputs.

The response of one diode in the array is shown in Figure 3.11. A square pulse of light lasting 200 ms is generated from a LED powered by a signal generator (Wavetec80, San Diego, CA, USA). The resulting signal is shown. Distortion due to low bandwidth is visible near the top of the pulse. Frequency response can be estimated from this figure by obtaining a time constant t, where Photodiode Output µ 1- exp^[T/(t)]. If we assume that the dominant limit on frequency is due to stray capacitance across the photodiode, and we assume a one pole role off, we can estimate the stray capacitance to be 5 picofarads¹². The response of the diode to shorter signals is also shown. A plot of signal size vs signal duration is given in Figure 3.12. The range of the circuit is ± 10 V, approximately (rail to rail). The noise level is roughly 5 mV. There is negligible interference with shielding, however a large 60 cycle (100 mV peak to peak) signal is observed if the shielding is removed.

3.5.2  Second Stage Circuitry.

Useful optical signals often are small signals that are summed with a large offset. For example, voltage sensitive dye Di-4-ANEPPS gives has gives a large fluorescent signal which only changes 5-10% with membrane changes of 100mV[82]. Ideally the offset should be subtracted off and the useful signal amplified before the signal is recorded¹³.

Initial experiments with a small 3 element array demonstrated that this approach does not work. Offset signals vary greatly from point to point within the Di-4-ANEPPS stained monolayer: offsets could vary from 300 mV to more than 2 V depending on location, presumably due to differential staining and irregularities in the surface of the monolayer. There is also a significant signal drift associated with dye bleaching (which acts to reduce offset)[259]. More importantly, the system is very sensitive to light noise: it proved to be very difficult to isolate the array so that movement from the operator or output from computer screens did not affect the output. Also, the quality of the signal varied from point to point in the monolayer. Typical signals were 4-8% of total offset, however useful signals were sometimes as low as 2% of total offset. The above considerations resulted in signals that were not recordable. As the signals vary by as much as 10× and the range of the digitizing equipment has to be set to encompass the largest signal, digitizing equipment cannot record unconditioned signals reliably. For example, a 8% signal riding on a 300 mV offset would result in a useful signal of 24 mV, which would be digitized with 32 point resolution. Even though this may be sufficient for some purposes, interference noise make it impossible to move a 24 mV signal from detector to digitizer without loss.

The amount of noise inherent in the system makes AC coupling impossible. Preliminary experiments showed that action potentials are not easily resolvable with a AC coupled signal. This is not surprising considering the signal to noise ratio is sometimes as low as 3/1.

A simple three element second stage was constructed that hand tunable offset voltage controls and variable amplification (Figure 3.13). The circuit allowed the offset to be removed and the signal amplified by as much as 500× prior to digitization. However slowly changing offsets, as well as room light noise required frequent operator adjustment. Manual adjustments are impractical for larger arrays.

Automatic offset control circuits have been developed for use with photodiode arrays[222]. These circuits use a large resettable capacitor to store DC levels which is then subtracted off the incoming signal. The operator can hit a switch to reset the capacitors in the second stage. Although this technique is simple and effective, it requires a large capacitor and an analog switch for each channel. The number and size of components make laying out a large number of automatic offset circuits impractical. In addition, capacitors drift over time. Automatic gain circuits have also been developed[220,222] which allow the operator to switch between a few preset gain levels. This approach requires digital logic for each channel. Even though signals are eventually digitized, the above circuits do not take advantage of the recording computers ability to monitor and change incoming signal.

A circuit was developed that allowed computer controlled offset and gain control for each photodiode. The second stage circuit takes advantage of the fact that information on each channel is available to computer algorithms in real time. The computer outputs control signals to the second stage circuit, adjusting gain and offset automatically.

As every channel is eventually digitized by computer, this system was designed to allow automatic monitoring of levels with a computer program which outputs signals to second stage circuitry to optimize the voltage levels prior to digitization.

Offset and gain control is accomplished by the circuits shown in Figure 3.14. A digital signal is outputted from the computer and changed to an analog signal by multichannel 8 bit digital to analog converters(D-A's). D-A's can be arranged in parallel so that a single 12 bit word can control 32 D-A channels. As the D-A's generate analog signals from a preset reference signal, one D-A in the chain is used to provide reference signals to the other channels, increasing the D-As resolution to 10 bits[117,278]. The circuit uses one analog channel to control the offset, and one channel to control the gain. As there are 8 D-As per chip in the AD7064 package used, one chip is required to condition 4 channels. Automatic gain control is accomplished by the circuit shown in Figure 3.14A. The circuit uses an analog offset line from the D-A's to control gain, and sums another analog line with the output of the photodiode. The circuit uses an original design that allows use of very inexpensive audio circuitry and a small number of components to control gain over a large range with very little noise.

Voltage controlled amplifiers for low noise applications exist but are often expensive and come in large package sizes, making them impractical for multichanel applications. Lower cost VCAs in quad chips are available. A typical unit is Analog Devices SSM2024. These units function as current attenuators - an applied control voltage can decrease the current output of the circuit, not increase it. The resulting current has to be amplified by second amplifier. In addition, these circuits function optimally over a very narrow current input range requiring an input voltage less than 50 mV. Using these circuits as suggested by the manufacturer has three disadvantages: The weak incoming signal from the photodiode stage has to be scaled down before entering the VCA, which will increase the relative noise contribution from the second stage. The level of the signal is not known before entering the VCA, so optimal performance from the VCA is not possible. Finally the signal is attenuated by the VCA before being amplified, further increasing the contribution of noise.

The following circuit was designed to remove these disadvantages of the SSM2024¹⁴. The VCA is placed in the feedback loop of an amplifier in a non inverting gain configuration (see Figure 3.5.2A). Unlike conventional gain circuits, the level of the input signal is not known, but the desired level is - we require a peak to peak signal that fills the A-D range of the digitizing board. As the output peak to peak signal is known, its possible to scale the input to the VCA so that its input is optimal. The VCA output acts as a resistor in the feedback loop of the operational amplifier. The signal from the photodiode stage inputs directly into the operational amplifier, protecting the signal.

The circuit shown is laid out so that the VCA replaces the feedback resistor R_gain in a typical noninverting operational amplifier application (Figure 3.14). The gain of the circuit shown is inversely proportional to the current supplied to the VCA input. The VCA itself outputs a current that is proportional to V_in [278]:

IVCA = 8.17 Icontrol (200/(Rin + 200 ) Vin).

(3.10)

where I_control must be no more than 500 mA. The SSM2024 attenuates the signal it receives. The feedback resistor in a typical noninverting operational amplifier circuit essentially performs the same task. The VCA is used in the second stage circuit as a voltage controlled variable resistor.

The VCAs total harmonic distortion (THD) performance and signal to noise ratio are determined by the peak to peak input signal level. There is an inverse relationship between the signal to noise ratio and THD. The SSM2024 can give a signal to noise ratio of 80dB with a THD of 1% for a 50mV peek to peek signal[278]. Increasing the signal size degrades THD to as much as 10% which is unacceptable, while decreasing signal size decreases the signal to noise ratio to as low as 60 dB.

Resistor R₄'s value of 200 ohms is recommended by the manufacturer as offering the best offset and control rejection[278]. R₃ scales the input to the VCA. As the input to R₃ is the operational amplifier output, and we know that that equals the 10 V peek to peek, R₃ can be calculated to be

50mV=(200/R3+200)10 V, R3 » 40 KW

(3.11)

The output signal from the SSM2024's pin 4 is given by Equation (3.10). For a 500 mA control current, the output of the VCA I(VCA)is 0.0002A, for a 100 mA control current, the output is 0.00005A. The circuit is expected to deal with signals as low as 50 mV peak to peak, but as high as 1 V peak to peak. For the maximum input control current of 500 mA we desire a gain

Vout/Vin=10,   so  Vin=10V.

(3.12)

and for the minimum control current we require a gain of 200.

As the inverting pin accepts minimal current (input resistance is ideally infinite for operational amplifiers), the current leaving the VCA must pass through R₅. If the voltage drop across R₅ is 1V, and the output current is 0.0002A,

R5 = V/I = 1V/0.0002A =  5 KW

(3.13)

For example, if the input was 100 mV p-p, then a gain of 100 is required. For the values of resistances calculated above, the control current is

0.1V = R5 IVCA,  IVCA = 0.1V/5KW = 2e-5A

(3.14)

from Equation (3.10),

Icontrol = IVCA/(8.17(200/(Rin+200)Vin), Icontrol = 2e-5A/(8.17(200/40200)10, Icontrol = 50 uA.

(3.15)

If a 0-10V variable voltage is applied to the VCA, 10V gets scaled by R₆ to be 500 mA, so R₆ = 20KW . The VCA passes no current if less than 200 mV is applied to the control input (this information comes from the SSM2024 data sheet), which corresponds to an input current of 10 mA. This gives a maximum bound to the gain at 500×.

Performance of the circuit for different levels of V_control and V_in are shown in Figures 3.15 and 3.16 The test voltage being input is a sine wave generated from a Wavetec signal generator. The output is digitized across ± 5V range at 12 bits. The resultant square wave is not distorted at either the high or low input ranges (20mV -1V). In comparison, if operated outside of a feedback loop, the VCA would not tolerate such a large (50×) swing in input levels without having unacceptably large shifts in THD and signal to noise ratio.

Chapter 4
Additional Methods.

4.1  Methods.

4.2  Monolayer Cell Culture.

Cells were plated at just under confluent densities of 7.0-14.0 × 10³ cells/cm² within 12mm diameter glass retaining rings on lysine coated plastic cell culture dishes. Cells were kept in maintenance medium 818a at 5% CO₂ at 36^° for one to two days before experiments. Monolayers were loaded with calcium sensitive dye Fluo-3 (5 mM) (Molecular Probes, Eugene, OR) in Tyrode for 25 - 30 minutes, washed with fresh tyrode to remove excess dye. The Tyrode solution containted NaCl 124.0 mM, NaHCO₃ 15.0 mM, KCl 1.3 mM, CaCl₂ 2.2 mM, MgCl₂ 1.0 mM, NaH₂PO₄ 1.0 mM, glucose 5.5 mM, and pH was adjusted to 7.3 with 2M NaOH. The monolayers were then transfered to the imaging setup. All experiments were conducted at 36^°.

All solutions were filtered through a sterile 0.22 mm pore size filter prior to use.

4.3  Atrioventricular Node Preparation

The preparation was stained with Di-4-ANEPPS (4-[b-[2-(di-n-butylamino)-6-napthyl]vinyl]pyridinium, Molecular Probes, Eugene OR) at a concentration of 20 mm for 20 minutes. The Di-4-ANEPPS was prepared as a stock solution in dimethyl sulfoxide (DMSO) and Pluronic acid (1:1 ratio). After staining the preparation was rinsed with fresh oxygenated Tyrode to remove excess dye.

4.4  Electrophysiology

Chapter 5
Pilot Studies 1: Optical Mapping The Atrioventricular Node.

5.1  Introduction

My contribution to this study include design¹⁵ of a stage and perfusable, heated bath for large preparations, dye loading and optical mapping the AV node, and presenting and analyzing the data with the software described in Chapter 3. Disection of the AV node, calibration of the bath, design of stimulation protocols, and subsequent interpretation of the data were carried out by Dr. Andrew Munk and Kevin Petrecca. Additional immunohistochemical studies which related the sites of delay to the concentration of sodium channels were conducted on the preparations by Kevin Petrecca et al.[193]. As the immunohistochemical results do not relate to optical mapping they are not presented here.

5.1.1  Background.

Ablation techniques (where a section of the AV node is surgically destroyed by a pulse from a catheter) have been successful in treatment of atrioventricular nodal reentrant tachycardia (AVNRT)[197]. The theoretical basis underlying these techniques is the concept dual AV nodal pathways, hypothesized by Moe et al.[172,178] This early work showed that the recovery curve in dogs can show a distinct discontinuity, suggesting a switch in conduction between two pathways. Mapping studies[2,125,172] suggest that the slow pathway activates the node from the posterior direction while the fast pathway activates the node anteriorly. AVNRT supposedly occurs when an impulse travels down one pathway and up the other, forming a reentrant circuit[242]. Cardiologists target one of the pathways in order to break this circuit. Remarkably, despite extensive investigation over the last 50 years[62,128,253], the location and existence of the pathways is a subject of continued debate.

There are two competing views of conduction through the AV nodal region[31]. One hypothesis is that conduction is slowed decrimentally[115] due to high coupling resistance accross the nodal region. The competing view is that delay is localized in a region of the node[191]. The discontinuous gap in conduction was inferred from intracellular electrode studies[16,191]. The cause of the discontinuity in conduction is likely due to an anatomical inexcitable gap[31,50] which causes the impulse to conductuct electrotonically.

Optical mapping techniques are better suited than single electrodes for constructing activation maps of the AV node[63,170]. Conventional electrode mapping studies require replicating experimental protocols while changing electrode placement. The AV node has numerous cell types and a complex three dimensional morphology[16,170,184]. Extracellular recording cannot resolve activation times of different cell types at appropriate space scales. Accurate mapping of activation times using intracellular electrodes is difficult as electrodes impale at unkown depths, resulting in measurement of different activation times at similar electrode positions[31,16,17].

Choi and Salama[31] and Efimov et al.[62,63] have used optical mapping methods to record from the rabbit AV node. Both groups used photodiode arrays to record optical signals from Di-4-ANEPPS stained rabbit hearts. Choi and Salama[31] found that the location of the primary conduction delay within the AV nodal region is between atrium and AV node. Propagation is discontinuous in the sense that the atria, AV node, and ventricles activate with periods of electrical silence between them. Efimov et al.[62] mapped activation patterns to the node, and concluded that activation occurs through the posterior slow pathway during normal sinus rhythm. The discontinuous conduction observed in [31] were not obvious in Efimovs' records.

5.2  AV Node Preparation.

The AV node was prepared as described in 4. A high resolution image of the node is shown in Figure 5.1A. The image lines up with the labelled diagram shown in Figure 5.1B.

5.3  AV Node Results.

Optical signals were captured with the CCD camera and signal processed as described in previous sections. Images were collected every 10 ms. Figure 5.2 shows an application of the signal summation technique discussed in Section 3.4. Boxes used for signal reconstruction are shown in Figure 5.2A. As fluorescence in Di-4-ANEPPS stained preparations decreases as cells depolarize, downward spikes represent peak activity as the wave passes through the node. Box 0 and 1 are excited simultaneously under conditions of spontaneous activity. A 20 ms delay occurs between boxes 1 and 2. This corresponds to an area near the compact node. This site of delay corresponds well with those elucidated from microelectrode studies[16,17].

Panels showing the progression of a impulse triggered in the atria are shown in Figure 5.2C. The images presented are shown in reverse video (low fluorescence is shown as brighter pixels.) Images were collected every 10 ms (7x7 pixels binned.) Spontaneous impulses originating from the sinoatrial node approached the nodal region as a broad wavefront oriented parallel to the AV ring. The wave front became narrower and curved posteriorly towards a region below the coronary sinus. The wave front then paused in before the central region of the node. The final part of nodal activation was characterized by a rapid and synchronous activation of an anterior strip of cells along the AV ring.

Figure 5.3 shows a spontaneous escaped beat from the AV node. Impulses progress in a reverse direction from that shown in Figure 5.2. Here the impulse originates in a more posterior portion of the node instead of the anterior strip of cells observed in the previous figure. A delay analagous to that seen in Figure 5.2 exists between activation of the central node (Panels 2-5) and the region above the coronary sinus (Panels 6-8.)

5.4  Discussion

The data shows that the delay is concentrated in a narrow region of the AV node. This favours the inexcitable gap hypothesis of conduction; here the excitation has to pass over an electrically dead region in the AV node, which introduces delay. An alternate hypothesis, that delay is spread evenly over the whole node, is not supported by this study.

Figure 5.2 is similar in activation pattern and timing to Figure 6 in Choi and Salama[31]. Based on the data in Figure 6, Choi and Salama conclude that the delay occurs between atrial and AV node tissue. They also conclude that spread of activation within the AV node (the compact node) occurs over 10 ms. These conclusions are confirmed by our recordings. The records of Efimov et al. (Reference [62], Figures 5 and 6) have what appears to be a different activation pattern. The activation passes from the posterior side of the node, and does not display a distinct delay. This apparent diffence may be partially artifactual as the scale, location, and format of the record is dissimilar. The difference may also be due to differences in experimental protocol. Both groups used contraction blockers (diacetyl monoxime (DAM) or butanedione monoxime (BDM)) at low concentrations. However Choi and Salama washed the the contraction blocker out of the preparation after 10-12 minutes, while Efimov et al. maintained constant concentrations of contraction blocker in the bath. Prolonged (>45 minutes [31]) exposure to contraction blockers are known to cause changes in action potential shape and duration, We used DAM at constant concentrations for the duration of the experiments. The data presented in this Chapter was collected within 30 minutes of loading with contraction blocker. Efimov et al. performed experiments on tissue for one to three hours[62].

Interpretation of optical records is complicated by overlapping action potentials within one recorded optical signal. This effect is pronounced in optical records of the AV node because of its complex three dimensional structure. Both Efimov et al. and Choi et al. recorded distinct spikes from cells at different depths. This allowed construction of activation maps in the z, y and z axis. The CCD system used in our study does not have sufficient temporal resolution to resolve spikes from different layers in this fashion. This limitation could be addressed in future studies by using a small photodiode array in conjunction with the CCD.

The CCD based system successfully measured conduction through the AV node in several preparations. These experiments were not pursued due to time constraints. Optical mapping experiments of the AV node are presently underway, with the goal of resolving the presence of different conduction pathways at different pacing rates.

Chapter 6
Pilot Studies 2: Heart Cell Monolayer

6.1  Introduction.

Mapping activation by imaging fluorescent activity from living tissue is technically challenging. Commercially available imaging systems are expensive and require significant technical ability to operate. In addition, commercially available solutions are not optimized for different experimental preparations. Researchers often resort to using custom built apparatus. Differences in dye choice, image field size, preparation thickness, image resolution, and recording duration can make techniques appropriate for one preparation useless in another. In addition, design and construction of imaging solutions requires familiarity with optics, electronics, and computer programming. Due to these constraints, there are relatively few laboratories that perform high speed optical imaging.

Neither I, or anyone in the laboratory have the necessary expertise to design an optical mapping system from the ground up. In addition there are no laboratories (anywhere) that have imaging setups optimized for long term recording at fast rates over large areas from very thin cardiac preparations. The final experimental design is based on considerations outlined in Chapter 3 and on pilot studies using different technologies.

Several pilot studies using voltage sensitive dyes and photodiodes were carried out using the photodiode array described in Chapter 3. These experiments are outlined below. Although the photodiodes produced high quality signals, technical difficulties associated with building a large array prompted us to investigate other options. Several researchers were approached for advice. I visited the laboratories of Dr. Grinvald and colleagues, Dr. Windisch and colleagues, Dr. Jalife and colleagues and Dr. Nelson Publicover and colleagues. Ratslaff and Grinvald's macroscope design was used as a template for the macroscope built in our laboratory. Dr. Windisch's laboratory provided technical help on monolayer staining and photodiode amplifier construction. Complete pilot studies using the monolayer preparation were performed at Dr. Jalifes laboratory and Dr. Publicovers laboratory. Voltage sensitive dye measurements using the setup at Dr. Jalifes laboratory did not yield measurable signals. The lack of signal was not surprising: Dr. Jalife's system was optimized for measuring signals from preparations several millimeters thick. Voltage sensitive dye experiments at Dr. Publicover's laboratory also failed to produce results. Dr. Publicover suggested using calcium sensitive dye Fluo-3. The success of these measurements (described below) prompted us to change our experimental focus from photodiode measurement of voltage sensitive dyes to CCD based measurements of calcium sensitive dye. Data from calcium imaging experiments in our laboratory are presented in Chapter 7.

The CCD based system was subsequently modified to allow voltage mapping of intact cardiac preparations. Pilot studies were performed on rabbit atrioventricular (AV) node preparations. The CCD based system successfully measured conduction through the AV node in several preparations. These experiments (Section 5) were not pursued due to time constraints.

6.2  Calcium Mapping

Initial wide area mapping experiments were carried out with the help of Dr. Nelson Publicover in the University of Nevada School of Medicine. Young, sparse cultures were stained with calcium sensitive dye and mapped at low magnification using a very sensitive commercially available optical mapping system ¹⁶. The results from one experiment are shown in Figure 6.1. The optical system was not able to map propagation, but did capture periodic increases in activity that spanned the monolayer preparation. The results suggested that weekly coupled cardiac cells are able to periodically synchronize over an area many orders of magnitude larger than the range of local cell interaction. The behaviour observed was interesting on several levels. Periods of rapid activity followed by periods of silence have been observed in many biological contexts, and is referred to as bursting behaviour. Bursting is usually thought to be due to local dynamics of single and small networks of cells, not the group activity of many thousands of cells. Since cardiac cells do not burst in isolation ¹⁷ the behaviour of this system may provide a novel mechanism for bursting in large groups of coupled non bursting units. In addition, the monolayer is of interest in the context of spontaneous synchronization of locally coupled oscillators, which is the subject of extensive mathematical investigation. However, the low spatial resolution of the data provided little insight into the mechanisms of the synchronized activity.

Fluo 3 stained monolayer preparations monitored using the Princeton Instruments camera are the focus of Chapter 7. Preliminary work has also been done with voltage sensitive dye Di-4-ANEPPS on both the monolayer preparation and a whole heart AV node preparation. This data is presented here in the context of evaluating the optical mapping systems performance.

6.3  Potential Mapping in Monolayer Preparations.

Monolayers stained with Di-4-ANEPPS were monitored using photodiodes and camera as discussed above. Monolayers used for Di-4-ANEPPS studies were plated at high densities to maximize signal. The following photodiode records are obtained using a 10x lens. Signal levels from the Di-4-ANEPPS monolayers were too low to record using microscopic objectives less powerful than 10×.

Figure 6.2 shows the output of a single photodiode recording spontaneous activity from a monolayer preparation. As discussed before, the signal is distorted due to the low frequency response of the circuit, however signal quality allows clear estimation of key action potential characteristics. An estimate of action potential duration (200 ms) and the presence of phase 4 (the slow depolarization to threshold) can be resolved using this technique. Upstroke velocity should not be estimated from the raw signal however, for three reasons. Voltage signals do not give an absolute measure of voltage, so dV/dt would have to be made using an assumption of what V_min and V_max are. The circuit used to obtain the signal has a frequency response of 20 Hz which would have to be accounted for in a measurement. Finally, the photodiode is obtaining a measure of potentials from multiple cells. At 10 × magnification, the MRD500 photodiode detects light from an average of 20 cells. If these cells fired simultaneously, the resulting signal would give a average of the rise times of the individual cells. However, if we assume the cells fire in succession (a wave of excitation travels across the monolayer) then the rise time from individual cells is distorted. A estimate of this distortion (called spatial blurring) as a function of spatial resolution was carried out in Guinea pig hearts by S. Girouard et.al.[94]. Rise time from a single cell (RT_cell) can be estimated from the optical signal (RT_op) as follows[95],

RTcell=RTop - r/q.

(6.1)

where q is conduction velocity and r is spatial resolution. The rise time from Figure 6.2 can now be given as 1 ms. If we assume that the action potential covers 80 mV over this time, Dv/Dt is approximately 0.0125 mV/ms.

Figure 6.3 shows a simultaneous recording from a monolayer preparation impaled with a microelectrode. The signal obtained shows that the optical signal corresponds closely to the signal obtained from the electrode.

Figure 6.4 shows optical signals from 2 photodiodes simultaneously. The image was collected using a 10x lens. The photodiodes are measuring points approximately 0.5 mm apart. The difference gives an estimate of the conduction velocity of signals in this preparation.

Chapter 7
Bursting Behaviour in Cardiac Monolayer Preparations

7.1  Introduction

Bursting behaviour, where cells rapidly fire followed by periods of silence, is observed in single neurons[81,215] and networks of cells such as pancreatic b-cell[28,49,230], cortical slices[231], and cultured monolayers of neurons[183]. Similar types of behaviour are seen in disease states such as paroxysmal (reoccuring) reentrant tachycardias[208] and epilepsy[57].

Since early studies demonstrated that bursting rhythms could originate from single cells[81], theoretical analyses of bursting have focused on cellular ionic mechanisms[25,37,251,229], or small network models[34,228]. The current work demonstrates that bursting also occurs in dispersed monolayer cultures of embryonic heart cells.

Sparsely cultured cardiac cells form an increasing number of gap junctional connections with neighbouring cells as they mature. Cells beat synchronously after about three days in culture. The results presented here show that bursting activity occurs prior to the onset of regular beating in most cultures. This observation is unexpected as mechanisms responsible for bursting in single and network systems are not present in cardiac monolayers.

Optical mapping experiments demonstrate that bursting in this system is driven by spatial patterns of activation in the form of repeated initiation and termination of rotors. Rotors are observed in such diverse systems as chemical reactions[136], slime mold cultures[85], slices of brain tissue[231], and the developing retina[101]. Rotors are common to two dimensional excitable media; any sufficiently large continuous system of coupled excitable cells is thought to be able to generate rotors of excitation under the right conditions[268].

Previous studies of excised slices of cardiac tissue also demonstrate that sustained spirals can be initiated by electrical stimulation, and proposed that spirals might be a mechanism for abnormally rapid cardiac arrhythmias (tachycardias)[4,46]. Although tachycardias in people can be initiated by electrical stimulation, the more usual circumstance is that they start spontaneously. In most instances the rhythms also terminate spontaneously. On occasion there is a repeated spontaneous paroxysmal initiation and termination of tachycardias[38]. Although it is likely that the tachycardia are caused by reentrant activity and possibly by rotors of activation, the mechanism of initiation and termination of the tachycardia is unclear. The monolayer preparation may be a model system for this class of arrhythmia.

Rotors do not form spontaneously in homogeneous systems[265], but rather are induced by special initial conditions such as cross field stimulation[46], or physically induced wave breaks[179,205]. Spirals can also form when a well timed pulse is initiated in the wake of a passing wave[201,266]. More recently, attention was focussed on the role of heterogeneity[205] and geometrical boundaries[1,102] in the creation of rotors. Spiral initiation in the monolayer appears qualitatively different; unidirectional propagation evolves directly to a rotor without the interaction of other waves or obstacles.

The observed dynamics are simulated with modified two dimensional excitable media models. Excitable cells are modeled by mathematical constructs that can undergo an excursion from a resting state if stimulated, and become refractory to stimulation for a finite time after being activated. Excitable cells are coupled to their neighbours allowing excited cells to trigger unrefractory cells to activate. Existing excitable media models support rotors but do not display bursting behaviour observed in the monolayer preparation. Tissue culture of embryonic chick myocytes display intrinsic spontaneous oscillation[146,109]. In addition, in many cardiac preparations[190,17], including cultured cardiac myocytes[154] there is a decrease in excitability, or fatigue, due to rapid stimulation. Excitable media models are modified by adding spontaneous activity and fatigue seen in the monolayer.

Excitable media are modeled using cellular automata or coupled differential equation models. Cellular automata use discrete time and space and simple update rules to simulate an excitable medium. The simplicity of cellular automata make them popular as models for physical and biological systems[30,103,104,179,252,254,276]. but they fail to adequately model effects of wave front geometry (the curvature effect) or recent activation history (dispersion) on wave speed. Coupled systems of partial differential equations have also been used extensively to model excitable media[24,42,59,74,77,78,133,173,216,249,267]. Coupled differential equation models are generally better at modeling dispersion and curvature effects, but are less easily manipulated and are more computationally expensive[92].

Bursting is modeled using two cellular automata models and coupled partial differential equation model. The purpose of presenting different model types to simulate the same behaviour is twofold. First, the minimal set of assumptions used in the modified cellular automaton models generates behaviour that is easily understood; the differential equation model is then used to simulate the behaviour in a more realistic setting. Second, different models with similar behaviours suggest the mechanism of initiation and termination of spirals are not artifacts of the modeling environment. The minimal set of assumptions necessary to generate bursting in these simple models suggest that spatially driven bursting may be seen in other physical systems.

Results from optical mapping of cultured cardiac monolayers is given in Section 2. Section 3 gives results from computer simulations of bursting in spatially extended systems. The results are discussed in Section 4.

7.2  Experimental Results

7.2.1  Rotors Cardiac Monolayer Preparations.

7.2.2  Spontaneously Bursting Preparations

In order to determine the spatio-temporal patterns of activity underlying the burst, fluorescent activity was recorded from a 1 cm² area. These recordings demonstrated the presence of a spontaneously initiated rotor associated with the burst (Fig. 7.6). The top trace in Fig. 7.6 shows a 25 s burst of increased Ca²⁺ recorded from a 1 mm² region. A peak in activity is created every time the activation front passes through the 1 mm² recording region. In order to study the activation sequences leading to the initiation and the termination of the burst, contour plots of activation times were constructed. A wave of activity was spontaneously initiated (Fig. 7.6A) but did not propagate uniformly. The wave was blocked at a point near the site of initiation leading to unidirectional propagation to the left. After the wave propagated about 2 mm (Fig. 7.6B) it doubled back, invading the previously blocked regions leading to two mirror-image rotors, with a common pathway in between them. For the initial 14 rotations, the top rotor excites the common pathway first (Fig. 7.6C). Then the top and bottom rotors simultaneously excited the common pathway (Fig. 7.6D). After 44 rotations, the top rotor was blocked, after which there remains only a single rotor (Fig. 7.6E). This wave was eventually blocked leading to termination of the burst (Fig. 7.6F).

The sequence shown in Fig. 7.6 was observed in 72/120 bursts. In 24/120 cases there was a similar initiation to a pair of mirror-image rotors, but there was no evolution to a single rotor before termination. In 24/120 rotor bursts, unidirectional propagation evolved directly to a single rotor, and then terminated. An example of uniderectional propagation resulting in a single armed rotor is shown in 7.7.

Monolayers do not necessarily display bursting dynamics. Bursting rhythms occur in cultures incubated for relatively short times, while regular periodic beating is observed after about 3 days in culture. The effect of incubation time on dynamics was systematically explored in 3 culture preparations. Fluorescence data was collected at different incubation times from spontaneously beating and extracellularly paced monolayers. The results from representative experiments are shown in Figures 7.8 and 7.9.

The spontaneous rhythms from 6 monolayers plated at identical densities from the same culture with varying incubation times are shown in Figure 7.8. Infrequent bursting is observed after 24h in culture. Bursting frequency and duration increases from 36 to 48 hours. By 54h activity becomes regular, with infrequent short breaks. At 68h the dynamics are periodic. Rotors were clearly observed at every time point except the one measured at 68h. In this case the periodic rhythm was generated by a regularly beating focus. At 54h, the regular beating was generated by a stable rotor.

The effect of periodic stimulation on rhythm is shown in Figure 7.9. The preparations are paced every 800 ms at 2x the minimum stimulation current necessary to illicit a pulse. After 24 hours incubation time, the monolayer is able to follow every beat for 23 stimuli before skipping 3 beats, at which point it follows the pacing electrode for 14 beats before failing completely. After 42 hours, pacing is able to maintain a 1:1 rhythm for 40 beats prior to skipping one beat. At this point the dynamics follow a 3:2 rhythm for 6 beats, then a 2:1 rhythm for 5 beats, and then go to a regular 3:1 rhythm for the duration of the protocol. After 68 hours incubation time, the preparation beat with a 1:1 rhythm for the duration of the protocol.

Figure 7.10 gives calculations of curvature vs wave front velocity for several preparations. The dependence of wave speed on wave curvature is roughly linear for this preparation. Five different preparations were tested, with over 30 points taken from each. The diffusion coefficient averaged 0.07 cm²/s, and planar wave speed was estimated to be 20 mm/s. The above determination of diffusion coefficient and planar wave speed depend on certain continuum conditions being met. These are discussed in Section 7.4.3.

7.3  Modeling

A simple cellular automata model was developed which demonstrates that the addition of certain assumptions can result in bursting behaviour in groups of spontaneously active, excitable, but nonendogenously bursting cells. An intermediate cellular automata model (see Section 1.3.2), and a continuous differential equation model were developed to incorporate measurements of D(see figure 7.10) and to confirm that the behaviour is not an artifact of discrete cellular automata. The intermediate cellular automaton model and the continuous differential equation model simulate propagation in continuous excitable media. There are difficulties with treating the bursting monolayer as a continuum. These are discussed in Section 7.4.3. Due to these difficulties and to the complexity of the parameter space of these models, simulation results from the two continuum models are presented primarily as demonstrations of bursting in continuous excitable media models. A detailed investigation of the parameter space for these models is beyond the scope of this work.

Mathematical models of excitable media are coupled systems of equations in which each element represents the excitation properties of a local area (usually a single cell or a small group of cells.) The properties of a cell or group of cells is modeled by equations that mimic the properties of an excitable medium; the equations display excitability (where a small perturbation over a certain threshold leads to a large excursion from a steady state which eventually relaxes back to the steady state), and refractoriness (the models cannot be easily perturbed for a certain time interval after excitation.)

In cellular automata models, there is a lattice of sites with discrete states and a synchronous updating rule that is identical for all sites. Each grid point has a current state that is determined by its own state and by it's neighbors states at the preceding time. grid points undergo the following transitions in state; resting ® excited ® refractory ® resting. In general, a grid point spends a set number of time steps in the excited and refractory states before returning to a resting (excitable) state. Cellular automata have been used as simplified models for a large range of phenomena including neural networks[68], cardiac arrhythmias[179], oscillating chemical reactions[90], and forest fires[30].

Continuous differential equation models are systems of partial differential equations coupled by a diffusion term. The coupled differential equations either directly model individual ion species (as the Hodgkin-Huxley squid axon model[114] does) or model the phenomenological character of excitation (as the FitzHugh-Nagumo equations do[78].) The differential equations representing local dynamics are coupled via a diffusion term to simulate travelling waves of excitation. Coupled differential equations are used to simulate wave propagation in neural[77] and cardiac[133] as well as chemical reaction systems[136].

None of these excitable media models display the repeated onset and offset of rotors observed in the cardiac monolayer preparation. We modify two cellular automata models and a coupled differential equation model to simulate bursting dynamics. The Greenberg-Hastings (GH) model[103,75] is a simple cellular automata model. The Gehrardt-Shuster-Tyson (GST)[90,91] model is an intermediate cellular automata that simulates properties of continuous excitable media. The FitzHugh-Nagumo (FN)[78] model is a coupled differential equation model that is phenomenological model of excitable media. The models are modified to account for the physiological properties of spontaneous activity and fatigue.

Spontaneous activity is simulated either by giving each grid point a different intrinsic periodicity, or giving all grid points a probability of firing when excitable. A mechanism for fatigue in the monolayer preparation may involve the influence of Ca⁺⁺ on gap junctional conductance. For rapid firing, Ca⁺⁺ builds up to levels which decrease gap junctional conductance. The decrease in gap junctional conductance occurs over many seconds allowing for gradual buildup of fatigue[58,190]. Changes in gap junctional conductance affects the diffusion coefficient D. The diffusion coefficient in cardiac preparations can be given by D=[(S_v r)/((C_m(1 mF/cm²))][130,269], where S_v is cell surface to volume ratio, and r is resistivity of the interstitial fluid, cytoplasm and gap junctions. This is modeled in the modified FitzHugh-Nagumo equations by decreasing the diffusion coefficient under conditions of rapid cell firing or by lowering the excitability of each grid point. Lowered gap junctional conductance is simulated in the cellular automata models by decreasing a fatigued grid points ability to stimulate its neighbours.

Waves are initiated by the chance random activation of a sufficient number of sites in a given neighborhood. Once initiated, a wave can be blocked from propagating in all directions, it can propagate as a broken wave in one direction, or it can propagate symmetrically to generate a target. Waves of high curvature are blocked more easily by refractory sites than waves of low curvature. As high curvature activation fronts advance, they become less curved and less vulnerable to block. The relationship between curvature and stability allow stable rotors to be generated from broken activation fronts created at initiation sites.

7.3.1  Greenberg-Hasting Model

The update rule for the state of a site is as follows: if 1 £ u_ij(t) £ (E+R), then u_ij(t+1)=u_ij(t)+1. If u_ij(t) = 0, then u_ij(t+1)=0, unless one of two conditions hold: (i) the number of sites with an excited state at time (t) in the neighborhood of a site at (ij) is greater than q_ij(t)=q+f_ij(t), where K_exci is a positive number indicating the threshold for excitation in the absence of fatigue and f_ij(t) is a fatigue term defined below; or (ii) x_ij(t) < p, where x_ij(t) is a random number uniformly distributed on the interval [0,1], and p is a positive number representing the probability for spontaneous activation of a site.

The update rule for f_ij(t) is as follows. If u_ij(t)=0 and u_ij(t+1)=1, then f_ij(t+1)=f_ij(t)+d; otherwise f_ij(t+1)=gf_ij(t) where g is a positive constant, 0 \l g\l 1. Thus, there is an increment of the fatigue term associated with excitation of a site, and exponential decay of the fatigue at other times.

By setting p, f and d equal to 0 we recover the Greenberg-Hastings model[103,75] of excitable systems.

To carry out the simulations we need to set E, R, K_exci, d, g and p and the neighborhood size. This is a rich array of parameters, and we have not systematically investigated behaviors throughout parameter space. For simplicity, we assume the neighborhood consists of the eight nearest neighbors of a given site. For this neighborhood planar wave propagation requires that K_exci £ 2, and we set q = 2. By choosing K_exci at the maximum value that allows planar wave propagation, we increase the likelihood for block and the initiation of rotors. We require that broken wave fronts curve to generate rotors as observed in the experiments. With K_exci=2, rotors will form if E ³ 5 and we select E=5. R is set equal to E+1=6 to prevent a wave from exciting sites directly in its refractory wake. Simulations with the above parameters with no fatigue give a rotor with a period of about 36 iterations. If f_ij > 1, then site (ij) will not be excited by a planar wave. For these parameters, a rotor of period N will not persist if (1+d)g^N > 1. To obtain long burst we need a slow buildup of fatigue. We arbitrarily set d = .05 and g = 0.999 based on the above equation. The probability of random firing at a site is 7.5 ×10^-4 per iteration. The mean burst time and time between bursts increases as p decreases.

The model shows initiation and annihilation of bursting rotors similar to experimental observations. Impulses are initiated by the chance simultaneous firing of several sites in close proximity to each other, so that additional sites in the neighborhood fire. In the example shown (Fig. 7.11) local heterogeneity prevents spread of excitation in all directions, resulting in unidirectional propagation upwards (Fig. 9A). The wave advances and doubles back to excite the region previously blocked (Fig. 7.11B), resulting in reentrant excitation in the form of two mirror-image rotors (Fig. 7.11C). Heterogeneity, induced by random firing in the model, breaks the symmetry (Figs. 7.11D) leading to a single rotor (Fig. 7.11E). The buildup of fatigue eventually makes propagation impossible and the rotor is annihilated (Fig. 7.11F). These features are similar to experimental observations (Fig. 7.6).

In the cellular automata model each site corresponds to approximately 0.1 mm², or 7-14 cells in the biological preparation. During the burst the period of rotation of the rotor is approximately 36 iterations, so each iteration step corresponds to approximately 20 ms. With this scaling, the time interval between bursts in the cellular automata model is about 1000 iterations or 20 seconds, which is comparable experimental observations (Fig. 4).

Simulations were carried out for models with neighbourhood sizes greater than one. Table 1 summarizes the results. Parameter values are set to allow planar wave propagation from an initiation site proportional to neighbourhood size. K_exci is set to the maximum possible for propagation from an activated planar wave segment of size r grid points times (2r+1) grid points. The value in brackets is the largest value of K_exci that allows planar wave propagation. E is set to the minimum value that allows the broken wave front to curve inward, and R is set to 1+E. P is the measured rotor period for the parameter values used. V is wave front velocity measured at the initiation site. Wave speed is measured in the direction of planar wave propagation. For these simulations, planar wave propagation is measured parallel¹⁸ to the orientation of the grid. Propagation speeds in directions diagonal to the grid orientation is lower[252]. Planar wave speed for these K_exci values is 1 grid point per time step for unbroken wave segments.

r	k	E	R	P	dx/dt
1	3 (3)	5	6	36	1/5
2	8 (10)	6	7	56	2/5
3	17 (21)	9	10	76	3/7
4	29 (36)	18	19	47	4/13
5	43 (55)	16	17	61	5/12

The rotor periods for four different values of K_exci is given for neighbourhood sizes 1 through 5 in table 7.3.1. The same values of E and R for the corresponding neighbourhood sizes are used as in the previous table. For each neighbourhood size, K_exci is varied between the values shown in the second column. Values of K_exci that do not result in a stable rotor are denoted with an 'x'. The simulations indicate how rotor period is expected to change under the influence of fatigue for different neighbourhood sizes.

r	Kexci	Rotor Period (P)
1	0 ® 3	x	12	12	36
2	5 ® 8	14	14	19	56
3	14 ® 17	20	27	42	76
4	26 ® 29	x	x	43	47
5	40 ® 43	30	41	45	61

Rotor Initiation in Simple Cellular Automata Models.

Consider a GH model with a Moore neighbourhood of size 1, K_exci=3, E=5, and R=6. There are 56 ([(8*7*6)/6!]) different positions that three active grid points can have in a given 3 ×3 neighbourhood. Cells can have activation states ranging from 1 to 5 in each of these configurations. At least one grid point has a value of 1 at the moment of initiation. There are three possible locations for this grid point, giving 3 ×5 ×5-14=61 different permutations (where 14 is the number of identical repeats) for the 56 spatial configurations at the initiation site. Of this 3416 (56 ×61) different configurations, a certain portion give rise to unidirectional propagation and rotor waves, while the rest result in target patterns.

The probability that wave initiation results in a rotor can be determined by iterating the model for each of the 3416 possible configurations. As this is impractical, one spatial configuration (with 61 different configurations of activation state) is analyzed in detail. The representative spatial configuration is chosen so that the results can be generalized to the remaining 55.

The two examples presented below demonstrate that different spatial configurations generate different patterns of activation even when the state of each grid point at wave initiation is similar. The numbers give the activation state (1 to 5) of each grid point. The symbol `+' denotes the extent of the neighbourhood at t=0. Refractory grid points are represented as 0 and excitable grid points are left blank.

Example 1. t=0 t=5 t=10 1 12321 2345432 121 135000531 +++ 13431 240000042 141 25052 13500 00531 +++ 13431 240000042 121 135000531 2345432 12321 1

Example 2. t=0 t=5 t=10 +4+ 303 1350 0531 1+1 15451 240000042 +++ 232 135000531 1 2450542 1234321 121

The first configuration gives rise to a target. The second results in unidirectional propagation that eventually evolves to a double armed rotor.

Once initiated, the advancing wave is vulnerable to the effects of grid points becoming refractory only while the wave is close to (within one neighbourhood distance of) the initiating grid points. The configuration that results in the fastest propagating wave is more stable than others because the wave is in the neighbourhood of the initiating grid points for the least amount of time.

The configuration where three grid points are activated in a row (as is shown in the first example) is more stable than other configurations as it gives rise to more rapidly propagating wave fronts than other configurations. Other configurations activate grid points in the center of the neighbourhood before the wave is able to advance outward. An example of this behaviour is seen by comparing the two activation patterns above. By t=5, the propagating wave from the `three in a row' activation pattern has advanced farther (in all directions) than the second disjoint configuration.

Analysis of the `three in a row' configuration gives an upper bound on the number of targets for all possible spatial configurations as it is least susceptible to the effects of its initiating grid points becoming refractory. Of the 61 different state configurations, 47 gave rise to propagating waves, 24 of which were targets, and 23 of which gave rise to unidirectional propagation and rotors. Therefore, slightly less than half of all propagating waves seen in this system will give rise to rotors. The actual number is higher as other configurations generate proportionately more rotors than targets. For example, the configuration

 
                        ++x 
                        +++ 
                        xx+ 

7.3.2  Bursting in Other Models

The Gerhardt-Shuster-Tyson (GST) Model

The GST model[90,91] is formulated as a cellular automata, but addresses the issues raised above. Excitation spreads from to a grid point if more than a critical number (K_exci) of neighbors within Moore neighbourhood of radius R_n centered at the grid point are excited at the previous time.

The GST model displays important physiological properties commonly observed in excitable media: dispersion and curvature effects. Dispersion refers to the decrease in propagation velocity as the time elapsed since the passage of a previous wave decreases. Curvature effects refer to the slowing of propagation velocity for waves with a high curvature[90,265,270]. Although the GST model does support propagating rotor waves[90,92] it does not show the spontaneous bursting observed in our experiments.

We extend the GST model by incorporating spontaneous oscillation and fatigue. The spontaneous activity was generated by assigning a random probability for non refractory grid points to fire. Fatigue is modelled by increasing K_exci as a consequence of rapid activity.

The GST model has an excitation variable U, which has the values 0 (unexcited) and 1 (excited), and a recovery variable V. Transition between excited and unexcited states is controlled by V. V increases by G_up every time step when U is 1 to a V_max, at which point U is set to 0; V is decreased by G_down to 0. A grid point is excitable if U=0 and V is less than variable V_exci. A grid point will become excited in the next time step (U(t+1)=1) if U(t)=0, V(t) < V_exci and more than K_exci neighbors are excited. Similarly an excited grid point can become prematurely unexcited if V(t) is less than V_rec and more than K_rec neighbors are unexcited. Dispersion effects are included in the model by making K_exci an increasing function of V, K_exci(V)=K⁰_exci+[r(2r+1)-K⁰_exci](V/V_exci. Similarly, K_rec is made a decreasing function of V (K_rec(V)=K⁰_rec+[r(2r+1)-K⁰_rec](V-V_max)/(V_rec-V_max).

Spontaneous firing of cells was incorporated into the GST model by assigning a probability that each grid point will spontaneously fire once V drops to 0. Fatigue is modelled by adding a fatigue term F to K_exci. A large constant F_add is added to F every time the grid point fires and F decays exponentially by multiplying it by F_sub (F_sub < 1) every time step of the simulation. F_add and F_sub are set so the rapid firing resulting from rotor wave activity causes a gradual buildup of fatigue.

Space and time units for the GST model are set by scaling the model to known values of planar wave speed (c) and diffusion coefficient (D). Gerhardt and coworkers[92,91] show that the normal velocity (V) of the wave front in the GST model is a linear function of wave front curvature1.2. For the GST model, the slope and intercept of the curvature relation are functions of K_exci and r[90] only. Plane waves (c) propagate at a speed equal to the radius of the extended Moore neighbourhood for values of K_exci used in these simulations[92] (we use K_exci less than or equal to 2r+1, the value necessary for plane wave propagation):

c=r l time step

(7.1)

Gerhardt and colleagues[92] found that D depends primarily on r, and has the following approximate relationship:

D(r) = 0.032(2r+1)2 l2 time step

(7.2)

where l is the distance between adjacent grid points. Following the method employed in reference[91], we solve the equations for D and c simultaneously setting D and c to known values. D and c for the monolayer preparation were approximated in Figure 7.10 to be 0.07 cm²/s and 2 cm/s respectively.

t= 0.546875 r2 (2r+1)2

(7.3)

l= 1.0937 r (2r+1)2

(7.4)

G_up and G_down are chosen to give a maximum action potential duration of 400 ms and a maximum refractory period of 200 ms.

Two versions of the GST model were developed. The first version uses a Moore neighbourhood with cells positioned evenly on a grid. The second version assigns a random location to each cell by displacing each cell a random angle and distance from its original grid point location. The neighbours of each cell are then determined by checking which cells fall within a circular radius Rc of each other. Both versions gave nearly identical results for rotor initiation and annihilation. The second version gave rise to smoother wave fronts. The original model (with Moore neighbourhood and a rigid grid) is presented here, as the scaling relation determined above applies strictly only to cells on a grid as used in Ref.[91].

Parameters were set to reasonable estimates for embryonic heart cells. For the simulation shown in Fig 7.12 we assume that the average action we set V_max=1000, g_up = 250, g_down = 500, K_exci = 6, K_rec=8, Rc=3, V_exc=600, V_rec=900, S_pup=40. The probability of a grid point spontaneously firing is 0.005 per iteration. With Rc=3, the time units are 100 ms. per time step and the space unit is 0.67 mm for the Moore neighbourhood simulation. Results are shown in Figures 7.13. Figure 7.13 is analogous to Figure 7.8 shown for the simple cellular automata in the previous section.

The FitzHugh-Nagumo (FN) Model

The GST model answers some of the inadequacies apparent with the simple modified GH simulation. First, rotor initiation can occur in with large numbers of neighbours arranged in a non-grid formation. Second, the model can be appropriately scaled to known values of D and c. The final concern is whether initiation is an artifact caused by the discrete nature of grid point firing. It can be argued that individual grid points in a excitable media would not fire independently; rather they are continuously influenced by the state of their neighbours by diffusion. The simulations in this section attempt to broadly address this issue by using a continuous time, continuous space approximation to an excitable medium.

In general, excitable media can be described simply by two partial differential equations[90]. One equation describes an excitation variable (u), and the other describes a recovery variable (v),

¶u ¶t = f(u,v) E + D1 ¶2u ¶x2

(7.5)

¶v ¶t = E g(u,v) + D2 ¶2v ¶x2

(7.6)

In cardiac preparations, u represents electrical potential and v represents the state of ion channel openness. D₁ and D₂ are diffusion coefficients expressed in space units² over time units. Parameter E is a factor that quantifies the repetitive rates of excitation and recovery. In most excitable media, the excitation variable changes on a faster time scale than the recovery variable, so E is set to a (small) fraction less than one.

The FitzHugh-Nagumo equations[77,187] used as as an excitable medium,

f(u,v) = u -u3/3 -v

(7.7)

g(u,v) = u + b-gv

(7.8)

Here, D₁ quantifies the diffusion coefficient of variable u. D₂ quantifies the diffusion coefficient of the recovery variable v and is set to zero. (In cardiac simulations, v represents permeability of ion channels and does not diffuse. Computations with D₂ set to 1 display similar qualitative features but have some quantitative differences, see [267] and references within.) The local behaviour of the FitzHugh-Nagumo equations depend on parameters E, b and g, whereas the diffusion coefficient D is closely tied to time and space scales of the excitable media. Parameter b controls threshold, larger values of b cause the FitzHugh-Nagumo equations to require larger stimuli before they undergo an excursion from steady state (a simulated action potential.) Low values of b cause the FitzHugh-Nagumo equations to oscillate spontaneously (the exact ranges for these behaviours depend on the other two parameters.) Parameter E controls the ratio of excitation to recovery rates. Higher values give slower rise time kinetics and allow faster recovery. Lower values give sharper rise times (as are seen in physical excitable media such as heart and nerve tissue) and longer refractory times. Parameter g does not have a critical effect on the overall dynamics for small E[267]; increasing g increases the trajectories excursion speed for both recovery and excitation.

Computations of the FN equations are carried out by constructing a two dimensional grid of points which calculate the FN equations at small time steps (dt). Diffusion is simulated by summing the difference in u at each grid point with the value of its four nearest neighbours times a factor [D/(h²)]. The term [D/(h²)] is proportional to the number of samples per space unit. Here we use h=0.5, sampling 2 times every space unit. The first set of computations were carried out at b = 0.7, E=0.3, and g = 0.5. For this parameter range, rotor waves have a period of 11.4 time units, and a plane wave propagates at a maximum (fully recovered) speed of two grid points per time unit[267]. As rotors in the monolayer preparation have an initial period of approximately 400ms and a planar wave speed of 2 cm/s, the time and space units are scaled to be 40 ms per time unit, and 0.4 mm per space unit. Using these values, the diffusion coefficient is (space unit² / time unit) 0.04 cm²/s¹⁹, which is half the value determined in Figure 7.10.

Spontaneous activity was added to the FN model in two ways. For the first set of simulations, any grid point can spontaneously fire with a probability P_f if it is sufficiently recovered (recovery is defined as any grid point with u < -1 and v > 0.6.). This method is similar to the technique used in the simple cellular automata described in the previous sections. Although random activation is not realistic, it maintains the scaling calculations defined above.

For the second set of simulations, spontaneous activity was generated by assigning random values (within a range) to parameters b and E. The fatigue term was modified to decrease the spontaneous activity as well as increase E; simulations with fatigue only added E resulted in transient phase waves in the simulation, which are not evident in the monolayer preparation. E was varied evenly between 0.27 and 0.37, and b was varied evenly between 0.49 and 0.69. By randomly setting E and b, the model represents an assortment of coupled cells with a range of excitability and intrinsic periodicities, which we may expect to find in recently plated monolayer preparations. The main difficulty with this approach is that is harder to generalize these results to other homogeneous models of excitable media.

Both approaches resulted in the generation of rotors. The result from one representative simulation is shown in Fig 7.14. An isochronal map of a typical burst is shown in Figure 7.14A. Figure 7.14B shows a series of bursts, with a close view of one of them (Fig 7.14C).

7.4  Discussion

7.4.1  Rotor Initiation in Other Systems

Rotors are well known from experimental[3,46,157,227,262], and theoretical studies[151,179,90,265]. There are several known mechanisms for rotor formation in excitable media. Spirals can be created via external manipulation or arise spontaneously.

Rotor formation in homogeneous excitable media involves the creation of a heterogeneity in refractoriness in the medium with a stimulus followed by the creation of a wave break near the heterogeneity by a specially timed second stimulus. Winfree[261,267] first demonstrated that spirals can be created in homogeneous media with a procedure called the `pinwheel experiment.' It involves crossing a spatial gradient of phase with a spatial gradient of stimulus intensity. The gradient of phase is created by the recent passage of a travelling wave, and the gradient of stimulus intensity is generated by creating a temporary stimuli on the medium. This mechanism creates a pair of mirror image rotors. A related mechanism is cross field stimulation[265]. A wave of excitation passes perpendicularly though a gradient of refractoriness created by a previous wavefront. The two waves are timed such that the second wave blocks close to the previous wave but can propagate further back. This results in a wave break and rotor formation. This mechanism has been used to generate spirals in cardiac slices[47].

Rotors form spontaneously in a variety of different excitable media. Unlike spirals initiated by external manipulation, heterogeneity of some form is thought to be necessary for spontaneous spiral formation. The primary mechanism put forth for the creation of these waves was that it was due to heterogeneity of refractoriness[1,179]. Two waves (from a distant source) propagate toward a finite region of enhanced refractoriness. If the time between the waves is sufficiently short, the second wave cannot penetrate into the refractory region. A wave break is then formed at the boundary of the heterogeneity. The broken wave moves along the heterogeneity boundary until the tissue becomes excitable, at which point the break propagates into the previously refractory region forming a rotor. Another mechanism for wave generation involves chance initiation of a wave during the vulnerable period of an excitable media. Pacemakers at arbitrary sites in the excitable media results in initiation of multiple circular wave fronts at different locations. When wave is initiated by chance in the wake of a previous wave, a wave break occurs which can result in a rotor. This mechanism has been observed directly in Disciodium cultures[201]. A third mechanism of rotor formation in heterogeneous systems involves the presence of obstacles in the path of rapidly propagating waves[1,102,239]. When the period between travelling waves is long, the medium at the obstacle is well recovered and the wave splits at the obstacle and rejoin at the other side. When the wave period is shorter, the medium at the obstacle is not sufficiently recovered. Due to dispersion and curvature effects, wave speed is slowed at the broken tips so that they travel outward from the obstacle some distance before rejoining. If the obstacle is sufficiently large, the wave tips can curl inward before having the chance to rejoin, and form two rotors.

rotor formation in the above examples involve the breakup of established stable waves (as in the pinwheel and crossfield stimulation protocols[261,267] and the obstacle experiments[1]) or the creation of a broken wavefront in the wake of an established wave (as occurs premature pulse protocol[201]) In the experiments and models of the bursting monolayer, the focus is on initiation of a broken wave front in media with heterogeneities created by local refractory cells or patches of lower excitability. The media is heterogeneous on small space scale; unlike the mechanisms discussed above, a travelling plane wave propagates easily through these obstacles without breaking. However, an initiating pulse is going to be affected by the heterogeneity; if we assume that the distance between adjacent grid points is smaller than the critical radius of excitation for the medium then a wave will only be triggered if a certain number of those grid points fire in near synchrony. It is most likely that close to the minimum number of grid points will fire synchronously (unless these points are already entrained.) If this is the case then the initiation site will be equal to or only slightly greater than the critical radius of the media, and therefore maximally sensitive to heterogeneity. The degree of heterogeneity in the monolayer is unknown. In the simulations, the media is designed to be heterogeneous `enough': planar and slightly curved waves travel easily, and only waves of the high curvature break.

7.4.2  Rotor Annihilation in Other Systems.

Spirals in the monolayer preparation are not observed to break up as observed in the model systems discussed above. Spiral break-up in these models give rise in many small spirals; it is this condition that is thought to result in fibrillation in the heart. A typical annihilation sequence for the monolayer spiral is that its period slows, followed by wave front propagation failure. Present modeling assumes that a parameter in the biological system changes in a way that slows and eventually prevents conduction. However, it is possible that the core simply increases to a size greater than the diameter of the monolayer, and that propagation would have continued if the medium was larger. Alternatively, the spiral tip could have undergone meander and collided with the edge of the preparation. This latter possibility is unlikely, as in most cases the core was stationary or only meandered slightly before annihilating.

7.4.3  Continuum Considerations

An important issue that must be addressed is whether the monolayer preparation can be modeled as a anisotropic continuous medium. Most excitable media models of cardiac conduction assume that conduction is continuous. Continuous excitable media models of cardiac conduction are known as cable equations and have the following form[269]:

¶Vm ¶t = - Jm Cm + Dx ¶2 Vm ¶x2 + Dy ¶2 Vm ¶y2 + Dz ¶2 Vm ¶z2

(7.9)

V_m is local membrane potential, J_m is the sum of local ion-channel current densities, C_m is membrane capacitance, and D_x, D_y and D_z are diffusion coefficients that control the relative contribution of voltage from neighbouring regions (in the x,y and z axis) have on the local membrane potential. Models which employ the diffusion coefficient (such as the models presented in Section 7.3.2) assume continuum conditions. An underlying assumption in continuous cable theory is that the effect of discrete intercellular connections can be averaged out over the space scales important for conduction.

Although it is likely that propagation in healthy heart can be treated as continuous[269], propagation in certain cardiac preparations is saltatory[236]; that is waves of excitation propagate more rapidly through individual cells than between them. Isreal and colleagues[120,121] recorded electrical potentials from chick cardiac monolayers cultured on a microelectrode array. The elements in the array were spaced to allow simultaneous recording from points within individual cells and between adjacent cells. Conduction delays of up to 0.41 ms were observed between electrodes recording from adjacent cells under some experimental conditions. Similar results were obtained using a photodiode array recording from cultured monolayers of rat cardiac cells[71]. Rhor and colleagues[214] lowered conduction in strands of cultured neonatal rat ventricular myocytes using gap junction uncouplers. They were able to lower conduction velocity from 47 cm/s to under 1 cm/s. At high levels of cellular uncoupling, conduction jumped between adjacent patches of simultaneously activated cells instead of propagating smoothly. Rhor's measurements of discontinuous conduction is similar to visual observations of young chick monolayers loaded with calcium sensitive dye. Observations under high power (20x or 40x) reveal that conduction occasionally progresses as a rapid series of jumps, while at other times conduction progresses smoothly. It is interesting to speculate that cells within the chick monolayer after one day in culture are coupled to the same extent as Rhor's 3-8 day old rat cultures treated with gap junction uncoupling agents.

Art Winfree[267,268,269,270] collected several `rules of thumb' based on cable theory, continuum conditions and observations of rotor behaviour in different continuous media. Some of these rules are applied below to measured values from the monolayer preparation in order to test whether its behaviours are consistent with continuous cable theory.

One definition for D provided by Winfree is D = [(l²)/(t)]. Here l is the space constant defined as the distance where effects of potential drop off exponentially, and t is the rise time of space clamped excitation. If we take our measured value of D as the known quantity, and use t = 0.001s we get l = 0.00837cm. If we assume that t can be a factor of 5 larger (due to slow upstroke velocity) then we have l = 0.01878cm. From an estimation of planar wave speed c₀ being approximated by Ö([D/(t)]), we get t = 0.018s, which seems unreasonably high. Several electrophysiological measurements of monolayer preparations were conducted by myself and Dr. Arkady Kunysz (unpublished results) and upstroke velocities were on the order of tens of millivolts per millisecond. The electrophysiological trace in figure 6.3 is representative of action potential shapes observed over many preparations. However, I cannot rule out the possibility that the monolayers used in the optical mapping studies displayed slower upstroke velocities.

One criteria for treating an excitable medium as continuous is if l is greater than 10 cell lengths[216]. Clearly, this is not the case for the monolayer as l is on the order of one cell length. Another criteria is that observations should be made at a scale greater than wavefront thickness[236]. Wavefront thickness can be approximated by t×c₀[269]which in our case is 0.001-0.005s ×2cm/s= 0.002 - 0.01 cm. If we apply this criteria, the bursting monolayer preparation can be treated as continuous for the observations made using the optical setup, but may be discontinuous at scales important for rotor initiation.

Another criteria applied by Winfree to the monolayer preparation (personal communication) is that D > 1/(t×cell density). In our case, monolayer density is 7-14 ×10³ cells per cm², resulting in a value that ranges between 1/7 and 1/14 for t = 0.001ms and 5/7 and 5/14 for t = 0.005ms. The value of D calculated in Fig.7.10 only exceeds the smallest of these values.

Winfree presents several equations that relate rotor core diameter (d) and period (P) to planar wave speed, diffusion coefficient and length and space scales. The observed diameter of spiral cores in different monolayer preparations was approximately 1-2 mm, however a reliable estimate is difficult as some records are of poor quality. For argument sake, we consider the maximum diameter of the core to be d = [(c_r r)/(p)] = 3mm where c_r is wavespeed at curvature present at the spiral core (taken here to be roughly equal to planar wave speed.) Winfree states that the expressions Ö{(c₀ pd / D)} and Ö{(c₀² P/D )} are equivalent and approximately equal to 9, and reports that minimum and maximum values for different continuous excitable media are 4 and 11, respectively[269]. The value of the expression based on rotor diameter is 5.2 for d=3 mm and 3 for d=1mm for the bursting monolayer prepaparation. Using the expression based on rotor period, we get 4.8 for P=0.4 and 7.4 for P=1. Depending on the criteria used, the bursting monolayer is either slightly out of or comfortably within the limits of rotor behaviour in continuous media.

A final concern is the relationship between curvature and wave speed (see Equation 1.2). To the best of my knowledge, all derivations of the curvature relation are based on continuum assumptions. The curvature relation was used to calculate D for the monolayer preparation (see figure 7.10). The linear dependence of wave speed on curvature could be taken as an indication that the monolayer is best represented as a continuum. However, cellular automata display curvature dependent wave speed effects when their threshold (K_exci) is greater than one. Gerhardt and colleagues[90] calculated the relationship between curvature effects in his cellular automata model and the continuous curvature relation based on diffusion (see Section 7.3.2.) Even though the dependence of wave speed on curvature in Gerhardt's model is due to discontinuous threshold parameter K_exci, the relationship between wave speed and curvature is linear. If the monolayer is best described as a discontinuous system, the value of D calculated in figure 7.10 should be treated as an abstract measure of slope rather than a true diffusion coefficient.

Cable theory also has implications for rotor block under the influence of fatigue. The cable equations have wave speed decrease proportional to the square root of D, where D is proportional to the resistivity of the media. Therefore conduction velocity falls proportional to the square root of resistance, but only is zero for infinite resistance. This implies that rotors would slow but not block in response to increases in gap junctional resistance. Rotor annihilation would then depend on the core diameter increasing to a point where it exceeds the diameter of the monolayer, which is not seen in the present experiments. James Keener modified versions of the cable equations to account for the discontinuous nature of conduction between cells[137,139]. These modified cable equations predict that conduction blocks at finite gap junctional resistances, which is observed in experiments where resistance is increased using blocking agents[53,214,236].

In conclusion, use of continuous differential equation models based on cable theory to model the monolayer preparation requires more thought. The estimation of D in Figure 7.10 depends on Equation 1.2, which is strictly valid for waves of low curvature (although it has been shown to hold for highly curved waves in experiments[84,85].) Different estimates of D may confirm continuum conditions in the monolayer preparation. Use of discontinuous cellular automata models that do not depend on continuum conditions is the most reasonable modeling choice until these issues are resolved. Another safe approach would be to develop models similar to those used in reference [139].

7.4.4  Conclusion

These results should be applicable to cardiology where the paroxysmal onset and offset of reentrant cardiac arrhythmias presents significant clinical management problems[38]. Potential other applications involve bursting in neural systems[177,54] and in pancreatic beta cells[169], but more work is needed to determine the actual patterns of spatio-temporal activation in these systems. The current work underscores the delicate factors that lead to the stabilisation and destabilization of rotating waves of activity independent of external interventions. Moreover, since biological systems are spatially extended, patterns of spatial activation must necessarily be investigated before the mechanism of bursting behavior can be assessed.

Bibliography

[1]

Agladze K, Keener J P, Muller S C, Panfilov A. Rotating spiral waves created by geometry. Science, 264:1746-1748, 1994.

[2]

Anselme F, Hook b, Monahan K, et al. Heterogeneity of retrograde fast-pathway conduction pattern in patiens with atrioventricular nodal reentry tachycardia: Observations by simultaneous multisite catheter mapping of Koch's triangle. Circulation, 93:960-968, 1996.

[3]

Allessie M A, Bonke F I M, Schopman F J G. Circus movement in rabbit atrial muscle as a mechanism of tachycardia: I. Circ. Res., 33:54-62, 1973.

[4]

Allessie M A, Bonke F I M, Schopman F J G. Circus movement in rabbit atrial muscle as a mechanism of tachycardia: III The leading circle. Circ. Res., 41:9-18, 1977.

[5]

Anderson R H, Janse M J, van Capelle F J L, Billette J, Becker A E, Durrer D. A combined morphological and electrophysiological study of the atrioventricular node of the rabbit heart. Circ. Res., 35:909-922, (1974).

[6]

Anosov D V, Arnold V I. Dynamical Systems I: Ordinary Differential Equations and Smooth Dynamical Systems. Springer-Verlag, Berlin 1980

[7]

Arnold V I. Small denominators I: mapping the circumference onto itself. AMS. Trans. Ser. 2, 46:213-284, 1965.

[8]

Bak P. The devils staircase. Physics Today, Dec, 38-45, 1986.

[9]

Bak P, Chen K, Tang C. Phys. Lett.A, 147:297-300, 1990.

[10]

Bar M, Eiswirth M. Turbulence due to spiral breakup in continuous excitable media Phys. Rev. E. 48 R1635-R1637, 1993.

[11]

Barkley D, Kress M, Tucherman S Spiral-wave dynamics in a simple model of excitable media: transition from simple to compound rotation. Phys. Rev. A., 42:2489-2491, 1990.

[12]

Beeler G W, Reuter H. Reconstruction of the action potential of ventricular myocardial fibres. J. Physiol. (Lond.), 268:177-210, 1977.

[13]

Bélair J, Glass L. Self-similarity in periodically forced oscillators. Phys. Let., 96A(3):113-116, 1983.

[14]

J. Bélair, L. Glass Universality and self similarity in the bifurcations of circle maps Physica 16D, 143-154, 1985.

[15]

Billette J, Janse M J, van Capelle F J L, Anderson R H, Touboul P, Durrer D. Cycle-length-dependent properties of AV nodal activation in rabbit hearts. Am. J. Physiol., 231:1129-1139, 1976.

[16]

Billette J. Atrioventricular nodal activation during premature stimulation of the atrium. Am. J. Physiol. (Heart Circ. Physiol. 21), 252:H163-H177, 1987.

[17]

Billette J, Metayer R, St-Vincent M. Selective functional characteristics of rate-induced fatigue in rabbit atrioventriular node. Circ. Res., 68: 790-799, 1988.

[18]

Bernstien R C, Frame L H. Ventricular reentry around a fixed barrier: resetting with advancement in an in vitro model. Circulation, 81:267-280 1990.

[19]

C. Boyd On the structure of the family of cherry fields on the torus. Ergod. Th. Dynam. Sys., 5: 27-36,1985.

[20]

Bonhoeffer T, Grinvald A. Iso-orientation domains in cat visual cortex are arranged in pinwheel-like patterns. Nature, 353(6343): 429-431, 1991.

[21]

P. C. Bressloff,J. Stark Neuronal dynamics based on discontinuous circle maps. Phys. Lett. A, 150(3,4) 1990.

[22]

Photodiode monitoring with op amps. Burr-Brown Application Bulletin AB075, Burr-Brown Corporation, Tuscon, AZ, USA. 1994.

[23]

Brown P, Maxfield B, Moraff H. Electronics for Neurobiologists. MIT Press, Boston, USA, 1973.

[24]

Cabo C, Pertsov A M, Baxter W T, Davidenko J M, Gray R A, Jalife J. Wave-front curvature as a cause of slow conduction and block in isolated cardiac muscle. Circ. Res., 75(6):1014-1028, 1994.

[25]

Canavier C C, Clark J W, Byrne J H. Simulation of the bursting activity of neuron R15 in Aplysia: role of ionic currents, calcium balance, and modulatory transmitters. J. Neurophys. 66(6):2107-2124, 1991.

[26]

Carpenter A. Bursting phenomena in excitable membranes. SIAM J. Appl. Math., 36:334-372, 1979.

[27]

Casdagli M, Rotational chaos in dissipative systems. Physica D, 29:365-386, 1988.

[28]

Chay T R, Keizer J. A minimal model for membrane oscillations in the panreatic B-cell. Biophys. J., 42: 181-190, 1983.

[29]

Chay T R, Kang H S. Role of single channel sochastic noise on bursting clusters of pancreatic b-cells. Biophys. J., 54:427-435, 1988.

[30]

Chen K, Bak P, Jensen M H. A deterministic critical forest fire model. Phys. Lett. A., 149 no 4:207-210, 1990.

[31]

Choi B R, Salama G. Optical mapping of atrioventricular node reveals a conduction barrier between atrial and nodal cells. Am. J. Physiol., 274(3 pt 2):H829-H845, 1998.

[32]

Cohen L B, Lesher S. Optical monitoring of membrane potential: methods of multisite optical measurement. in Optical Methods in Cell Physiology, P. de Weer and B. Salzberg eds. Wiley-Interscience, New York, (1986).

[33]

Cohen A H. Effects of oscillator frequency on phase-locking in the lamprey central pattern generator. J. Neurosci. Methods, 21: 113-125, 1987.

[34]

Cohen A H, Dobrov T A, Li G, Kiemel T, Baker M T. The development of the Lamprey pattern generator for locomotion. J. Neurobiol., 21(7):958-969, 1990.

[35]

Cohen J, Fozzard H, Shue S S. Increase in intracellular sodium ion activity during stimulation in mammalian cardiac muscle. Circ. Res., 50:651-662, 1982.

[36]

Colizza D, Guevara M R, Shrier A. A comparitive study of collagenease- and trypsin- dissociated embryonic heart cells: reaggregation, electrophysiology, and pharmacology. Can. J. Physiol. Pharmacol., 61:408-419, 1983.

[37]

Cook D L, Satin L S, Hopkins W F. Pancreatic B cells are bursting, but how? TINS 14:411-414, 1991.

[38]

Coumel, P., Leclercq, J.-F., & Slama, R. Chapter 50, in Cardiac Electrophysiology and Arrhythmias eds. Zipes, D. P, Jalife, J. (Grune & Stratton, Orlando), 1985.

[39]

Courtemanche M, Glass L, Rosengarten M D, et. al. Beyond pure parasystole: promises and problems in modelling complex arrhythmias. Am. J. Physiol., 257:693-706, 1989.

[40]

Courtemanche M, Glass L, Belair J, Scagliotti D, Gordon D. A circle map in a human heart. Physica D., #40, 299-310, (1989).

[41]

Courtemanche M, Glass L, Rosengarten M D. Modelling ventricular parasystole. Mathematical approaches to cardiac arrhythmias Annals of the New York Academy of Sciences, vol 591 1990.

[42]

Courtemanche M, Winfree A T. Re-entrant rotating waves in a Beeler-Reuter based model of two dimensional cardiac conduction. Int. J. Bifurcation and Chaos, 1:431-444, 1991.

[43]

Courtemanche M, Glass L, Keener J P. Instabilities of a propagating pulse in a ring of excitable media. Phys. Rev. Lett., 70:2182-2185, 1993.

[44]

Courtemanche M. Complex spiral wave dynamics in a spatially distributed ionic model of cardiac electrical activity. Chaos, 6(4):579-600, 1996.

[45]

Daniel E E, Bardakjian B L, Huizinga J D, Diamant N E. Relaxation oscillator and core conductor models are needed for understanding of GI electrical activities. Am. J. Physiol., 266(3 pt 1):G339-G349, 1994.

[46]

Davidenko J M, Kent P F, Chialvo D R, Michales D C, Jalife J. Sustained vortex-like waves in normal isolated ventricular muscle. Proc. Natl. Acad. Sci. USA, 87:8785-8789, 1990.

[47]

Davidenko J M. Spiral wave activity: a possible common mechanism for polymorphic and monomorphic ventricular tachycardias. J. Cardiovasc. Electrophysiol., 4(6):730-746, 1993.

[48]

Davila H V, Salzberg B M, Cohen L B, Waggoner A S. A large change in axon fluorescence that provides a promising method for measuring membrane potential. Nature New Biol., 241:159-160, 1973.

[49]

Dean P M, Matthews E K. Electrical activity in pancreatic islet cells: effect of ions. J. Physiol., 210:265-275, 1970.

[50]

DeFelice L J, Challice C E. Anatomical and ultrastructural study of the electrophysiological atrioventricular region of the rabbit. Circ. Res., 24:457-474, 1969.

[51]

DeHaan R L. Regulation of spontaneous activity and growth of embryonic chick heart cells in tissue culture. Developmental Biol., 16:216-249, 1967.

[52]

DeHaan R L, Fozzard H A. Membrane response to current pulses in spheroidal aggregates of embryonic heart cells. J. Gen. Physiol., 65(2):207-222, 1975.

[53]

Delmar M, Michaels D C, Johnson T, Jalife J. Effects of increasing intercellular resistance on transverse and longitudinal propagation in sheep epicardial muscle. Circ. Res., 60:780-785, 1987.

[54]

Destexhe A, Contreras D, Sejnowski T J, Steriade M. A model of spindle rhythmicity in the isolated thalamic reticular nucleus. J. Neurophys., 72:803-818, 1994.

[55]

Diaz P J, Rudy Y, Plonsey R. Intercalated discs as a cause for discontinuous propagation in cardiac muscle: A theoretical simulation. Ann. Biomed. Eng., 11:177-189, 1983.

[56]

Dillon S, Morad M. A new laser scanning system for measuring action potential propagation in the heart. Science, 214(4519):453-456, 1981.

[57]

Dudec F E, Snow R W, Taylor C. Role of electrical interactions in synchronization of epileptiform bursts. Advances in Neurology, 44:593-617, 1986.

[58]

De Mello W C. Effect of intracellular injection of calcium and strontium on cell communication in heart. J. Physiol., 250:231-245, 1975.

[59]

DiFranscesco D, Noble D. A model of cardiac electrical activity incorperating ionic pumps and concentration changes. Phil. Trans. R. Soc. Lond., B307:353-376, 1985.

[60]

Ebner T J, Chen G. Use of voltage-sensitive dyes and optical recordings in the central nervous system. Prog Neurobiol., 46(5):463-506 1995.

[61]

Eddlestone G T, Goncalves A, Bangham J A, Rojas E. Electrical coupling between cells in Islets of Langerhans in mouse. J. Membr. Biol., 77:1-14, 1984.

[62]

Efimov I R, Fahy G J, Yuanna C, van Wagoner D R, Tchou P J, Mazgalev T N. High-resolution fluorescent imaging does not reveal a distinct atrioventricular nodal anterior input channel (fast pathway) in the rabbit heart during sinus rhythm. J. Cardiovasc. Electrophysiol., 8:295-306, 1997.

[63]

Efimov I R, Mazgalev T N. High-resolution, three dimensional fluorescent imaging reveals multilayer conduction pattern in the atrioventricular node. Circulation, 98:54-57, 1998.

[64]

El-Sherif N, Smith A, Evans K. Canine ventricular arrhythmias in the late myocardial infarction period. 8. Epicardial mapping of reentrant circuits. Circ. Res., 49:255-265, 1981.

[65]

Endresen B. Chaos in weakly-coupled pacemaker cells. J. Theor. Biol., 184(1):41-50, 1997.

[66]

Fabiato A. Fluorescence and differential light absorption recordings with calcium probes and potential-sensitive dyes in skinned cardiac cells. Can. J. Physiol. Pharmac., 60(4):556-567, 1982.

[67]

Fabiato A. Calcium-induced release of calcium from the sarcoplasmic reticulum. J. Gen. Physiol., 85:247-320, 1985.

[68]

Farley B G. Computers in Biomed. Res.,1:265-294, 1965.

[69]

Fast V G, Kleber A G. Anisotropic conduction in monolayers of neonatal rat heart cells cultured on collagen substrate. Circ. Res., 75:591-595, 1994.

[70]

Fast V G, Kleber A G. Block of impulse propagation at an abrupt tissue expansion: Evaluation of the critical strand diameter in 2- and 3- dimensional computer models. Cardiocasc. Res., 30(3):449-459, 1995.

[71]

Fast V G, Darrow B J, Saffitz J E, Kleber A G. Anisotropic activation spread in heart cell monolayers assessed by high-resolution optical mapping. Circ. Res., 79:115-127, 1996.

[72]

Fast V G, Etimov I R. Stability of vortex rotation in an excitable cellular medium. Physica D, 49:75-81, 1991.

[73]

Feely O, Chua L. Nonlinear dynamics of a class of analog to digital converters. Int. J. Bifurcation and Chaos, 2(2):325-340, 1992.

[74]

Field R J, Koros E, Noyes R M. Oscillations in chemical systems II, Thorough analysis of temporal oscillations in the bromate-cerium-malonic acid system. J. Am. Chem. Soc., 94:8649-8664, 1972.

[75]

Fisch R, Gravner J, Griffeath D. Statistics and Computing, 1:23-29, 1991.

[76]

FitzHugh R. Mathematical models of threshold phenomena in the nerve membrane. Bull. Math. Biophys., 17:257-278, 1955.

[77]

FitzHugh R. Pulses and physiological states in theoretical models of nerve membrane. Biophys. J., 1:445, 1961.

[78]

FitzHugh R. Mathematical models of excitation and propagation in nerve. in H.P. Schwan, editor, Biological Engineering. McGraw-Hill, New York, (1969).

[79]

Fluhler, Bernham, Loew. Spectra, membrane binding and potentiometric reponses of new charge shift probes. Biochemistry, 24:5749-5755, 1985.

[80]

Fowles G. Introduction to Modern Optics. Dover Publications, 328 pgs, (1989)

[81]

Frazier W T, Kandel E C, Kupfermann I, Waziri R, Coggeshall R E. J. Neurophysiol. 30:1288-1351, 1967.

[82]

Fromherz P, Lambacher A., Spectra of voltage sensitve flourescence of styryl-dye in neuron membrane. Biochimica et Biophysica Acta, 1068:149-156, 1991.

[83]

Foerster P, Muller SC, Hess B. Curvature and propagation velocity of chemical waves. Science, 241:685-687.

[84]

Foerster P, Muller S C, Hess B. Critical size and curvature of wave formation in an excitable chemical medium. Proc. Natl. Acad. Sci. USA., 86:6831-6834, 1989.

[85]

Foerster P, Muller S C, Hess B. Curvature and spiral geometry in aggregation patterns of Dictyostelium discoidium. Development, 109:11-17, 1990.

[86]

Frame L H, Simsom M B. Oscillations of conduction, action potential duration, and refractoriness: A mechanism for spontaneous termination of reentrant tachycardias. Circulation, 78:1277-1287, 1988.

[87]

Frazier D F, Wolf P D, Wharton J M, Tang A S L, Smith W M, Ideker R E. Mechanism for electrical initiation of reentry in normal canine myocaridium. J. Clin. Invest., 83:1039-1047, 1989.

[88]

Gambaudo J M, Glendinning P, Tresser C. The rotation interval as a computable measure of chaos. Phys. Lett., 105A(3):97-100, 1984.

[89]

Gerhardt M, Schuster H. A cellular automaton describing the formation of spatially ordered structures in chemical systems. Physica D., 36:209-221, 1989.

[90]

Gerhardt M, Schuster H, Tyson J J. A cellular automaton model of excitable media including curvature and dispersion. Science, 247:1563-1566 1990.

[91]

Gerhardt M, Schuster H, Tyson J J. cellular automaton modle of excitable media II. Curvature, dispersion, rotating waves and meandering waves. A Phys. D, 46(3):392-415, 1990.

[92]

Gerhardt M, Schuster H, Tyson J J. A cellular automaton modle of excitable media III. Fitting the Belousov-Zhaboinski reaction. Phys. D, 46(3):416-426 1990.

[93]

Gerisch, G. Chapter 6. in Current Topics in Developmental Biology, Vol. 3, eds. Moscona, A. Monroy, A. (Academic Press, New York), 1986.

[94]

Girouard S D, Laurita K R, Rosenbaum D S. Unique properties of cardiac action potentials recorded with voltage-sensitive dyes. Journal of Cardiovascular Electrophysiology, 7(11):1024-1038, 1996.A

[95]

Girouard S D, Pastore J M, Laurita K R, et al. Optical mapping in a new guinea pig model of ventricular tachycardia reveals mechanisms for multiple wavelengths in a single reentrant circuit. Circulation, 93:603-613, 1996.B

[96]

Glass L, Goldberger A, Belair J. Dynamcs of pure parasystole. Am. J. Physiol. 251:H941-H847, 1986.

[97]

Glass L, Mackey M C. From Clocks to Chaos: The Rhythms of Life. Princeton U.P, Princeton NJ, (1988).

[98]

Glass L, Zeng W Z. Complex bifurcations and chaos in simple theoretical models of cardiac oscillations. Mathematical Approaches to Cardiac Arrhythmias, Annals of the New York Academy of Sciences, Vol 591:316-327, 1990.

[99]

Glass L, Josephsom M E. Resetting and annihilation of reentrant abnormally rapid heart beat. Phys. Rev. Lett. 2059-2063, 1995.

[100]

Gonzalez C, Rodriguez M. A flexible perforated microelectrode array probe for action potential recording in nerve and muscle tissues. J. Neurosci. Meth., 72(2):189-195, 1997.

[101]

Gorelova N A, Bures J. Spiral waves of spreading depression in the isolated chicken retina. J. Neurobiol., 14(5):353-363, 1983.

[102]

Graham M D, Kevrekidis I G, Asakura K, Lauterbach J, Krischer K, Rotermund H H, Ertl G. Effects of boundaries on pattern formation: Catalytic oxidation of CO on platinum. Science, 264:80-82, 1994.

[103]

Greenberg J M, Hastings S P. Pattern formation and periodic structures in systems modeled by reaction diffusion equations. Bull. Amer. Math. Soc., 6:84-, 1978.

[104]

Greenberg J M, Hassard B, Hastings S P. Spatial patterns for discrete models of diffusion in excitable cells. SIAM J. Appl. Math., 34:515-523, 1978.

[105]

Grinvald A. Real-time optical mapping of neuronal activity from single growth cones to the intact mammalian brain. Ann. Rev. Neurosci., 8:263-305, 1985.

[106]

Grinvald A, Bonhoeffer T, Malonk D, Shoham D, Bartfield E, Arieli A, Hildesheim R, Ratzlaff E. Optical imaging of archetecture and function in the living brain. in Memory: Organization and Locus of Change. editors Squire L R, Weinberger N M, Lynch G, McGaugh J L Oxford University Press, (1991).

[107]

Grinvald A, Frostig R D, Lieke E, Hildesheim R. Optical imaging of neuronal activity. Physiological Reviews, 68(4):1285-1366, 1988.

[108]

Guevara M R, Glass L, Shrier A. Phase locking, period-doubling bifurcations, and irregular dynamics in periodically stimulated cardiac cells. Science. 214:1350-1353, 1981.

[109]

Guevara M R, Glass L. Phase locking, period doubling bifurcations, and chaos in a mathematical model of a periodically driven oscillator: a theory for the entrainment of biological oscillators and the generation of carediac dysrhythmias. J. Math. Biol., 15:339-349, 1982.

[110]

Hao B.-L. Elementary Symbolic Dynamics and Chaos in Dissipative Systems. (World Scientific, Singapore), (1989).

[111]

Harrison L, Ideker RE, Smith WM, Klein GJ, Kasell J, Wallace AG, Gallagher JJ. The sock electrode array: a tool for determining global epicardial activation during unstable arrhythmias. Pacing Clin. Electrophysiol., 3(5), 531-540, 1980.

[112]

Y. He, R. Chi, N. He. Lyapunov exponents of the circle map in human hearts. Phys. Lett. A., 170:29-32, 1992.

[113]

M. R. Herman in Geometry and Topology, edited by J. Palis, Lecture Notes in Mathematics Vol 597.

[114]

Hodgkin A L, Huxley A F. A quantitative description of membrane currents and its application to conduction and excitation in nerve. J. Physiol., (Lond.) 117:500-544, 1952.

[115]

Hoffman B F, Cranefield P F. Electrophysiology of the heart. New York: McGraw-Hill, 1960.

[116]

Honerkamp J, Mutschler G, Seitz R. Coupling of a slow and a fast oscillator can generate bursting. Bull. Math. Biol., 47:1-21, 1985.

[117]

Horowitz P, Hill W. The art of electronics. Cambridge University Press, (1989).

[118]

Ideker RE, Smith WM, Blanchard SM, Reiser SL, Simpson EV, Wolf PD, Danieley ND. The assumptions of isochronal cardiac mapping. Pacing Clin. Electrophysiol., 12(3):456-478, 1989.

[119]

Ikeda T, Uchida T, Hough D, Lee J J, Fishbein M C, Mandel W J, Chen P S, Karagueuzian H S. Mechanism for spontaneous termination of functional reentry in isolated canine right atrium. Circulation, 94:1962-1973 1996.

[120]

Israel D A, Barrey W, Edell D, Mark R. An array of microelectrodes to stimulate and record from cardiac cells in culture. Am. J. Physiol., 247:H669-H674,1984.

[121]

Israel D A, Edell D J, Mark R G. Time delays in propagation of cardiac action potentials. Am. J. Physiol., 258:H1906-H1917, 1990.

[122]

Ito H, Glass L. Spiral breakup in a new model of discrete excitable media. Phys. Rev. Lett., 66(5):671-674, 1991.

[123]

Ito H. Mesoscopic modeling of wave propagation in excitable media. Physica D., #79:16-40, 1994.

[124]

Jacklet J, Barnes S, Bulloch A, Lukowiak K, Syed N Rhythmic activities of isolated and clustered pacemaker neurons and photoreceptors of Aplasia retina in culture. J. Neurobiol., 31(1):16-28, 1996.

[125]

Jackman W M, Beckman K J, McClelland J H et al. Treatment of supraventricular tachycardia die to atrioventricular nodal reentry, by radiofrequency catheter ablation of slow-pathway conduction. N. Engl. J. Med., 327:313-318, 1992.

[126]

Jalife J, Moe G J. Excitation, conduction and reflection of impulses in isolated bovine and canine cardiac fibers. Circ. Res., 49:233-247, 1981.

[127]

Janse M J, Kleber A G, Capucci A. Electrophysiological basis for arrhythmias caused by acute ischemia. J. Mol. Cell Cardiol., 18:339-355, 1986.

[128]

Janse M J, McGuire M A, Loh P, Thibault B, Hocini M, de Bakker J M. Electrophysiology of the A-V node in relation to A-V nodal reentry. Jpn Heart J., 37(5):785-791, 1996.

[129]

M. H. Jensen, P. Bak, T. Bohr, Complete devils staircase, fractal dimension, and universality of mode-locking structure in the circle map. Phys. Rev. Lett., 50(21):1637-1639, 1983.

[130]

Joyner R W, Veenstra R, Rawling D, Chorro A. Propagation through electrically coupled cells: Effects of a resistive barrier. Biophys. J., 45:1017-1025, 1984.

[131]

Kaneko K. Pattern dynamics and spatiotemporal chaos. Pattern selection, diffusion of defect and pattern competition intermittency. Physica D., 34:1-41, 1989.

[132]

Kaneko K. Overview of coupled map lattices. Chaos, 2(3):279-282, 1992.

[133]

Karma A. Spiral breakup in model equations of action potential propagation in cardiac tissue. Phys. Rev. Lett., 71:1103-1106, 1993.

[134]

Keizer J, Smolen P. Bursting electrical activity in pancreatic B-cells caused by Ca²⁺ and voltage inactivated Ca²⁺ channels. Proc. Natl. Acad. Sci. USA 88:3897-3901, 1991.

[135]

J.P. Keener. Chaotic behaviour in piecewise continuous difference equations. Transactions of the American Mathematical Society, 261(2):589-603, 1980.

[136]

Keener J P, Tyson J J. Spiral waves in the Belusov-Zhabotinskii reaction. Physica D. 21:269-324, 1986.

[137]

Keener J P. Propagation and its failure in coupled systems of discrete excitable cells. SIAM J. Appl. Math., 47:556-572, 1987.

[138]

Keener J P. The dynamics of three-dimensional scroll waves in excitable media. Physica D, 31:269-276, 1988.

[139]

Keener J P. Wave propagation in myocardium. in Theory of Heart: Biomechanics, Biophysics, and Nonlinear Dynamics of Cardiac Function. L Glass, P Hunter and A McCulloch eds. Springer-Verlag New York, (1991).

[140]

Kimura H, Oyamada Y, Ohshika H, Mori M, Oyamada M. Reversible inhibition of gap junctional intercellular communication, synchronous contraction, and synchronism of intracellular Ca2+ fluctuation in cultured neonatal rat cardiac myocytes by heptanol. Exp. Cell Res., 220(2):348-356, 1995.

[141]

Knisley S B, Hill B C. Effects of bipolar point and line stimulation in anisotropic rabbit epicardium: assessment of the critical radius of curvature for longitudinal block. IEEE Trans. Biomed. Eng., 42(10):957-966, 1995.

[142]

Kleber A G. Mechanisms of ventricular arrhythmias: A perspective. Journal of Cardiovascular Pharmacology, 17(Suppl. 6):S1-S8, 1991.

[143]

Kleber A G, Janse M J, Wilms-Schopman F J G, Wilde A A M, Coronel R. Changes in conduction velocity during acute ischemia in ventricular myocardium of the isolated porcine heart. Circulation, 73:189-198, 1986.

[144]

Kleber A G, Riegger C B, Janse M J. Electrical uncoupling and increase of extracellular resistance after induction of of ischemia in isolated, arterially perfused rabbit papillary muscle. Circ. Res., 61:271-279, 1987.

[145]

Fast V G, Kleber A G. Role of wavefront curvature in propagtion of cardiac impulse. Cardiovascular Research, 33:258-271, 1997.

[146]

Koidl B, Tritthart H A, Erkinger S. The effects of ouabain on the electrical and mechanical activities of embryonic chick heart cells in culture. J. Molec. Cell. Cardiol., 12:165-178, 1980.

[147]

Kopell N, Ermentrout G B. Symmetry and phaselocking in chains of weakly coupled oscillators. Commun. Pure. Appl. Math., 39:623-660, 1986.

[148]

Kreisman NR, LaManna JC, Liao SC, Yeh ER, Alcala JR. Light transmittance as an index of cell volume in hippocampal slices: optical differences of interfaced and submerged positions. Brain Res., 25;693(1-2):179-186, 1995.

[149]

Krinsky V I. Problem. Kibern., 20:59-80, 1968.

[150]

Krinsky V I. Biophysics, 16:88-, 1971.

[151]

Kopell N, Howard L N. Target patterns and spiral solutions to reaction-diffusion equations with more than one space dimension. Adv. Appl. Math., 2:417-449, 1981.

[152]

Kukuljan M, Goncalves A, Atwater I. Charybdotoxin-sensitive K-Ca channel is not invovled in glucose-induced electrical activity of pancreatic B-cells. J. Membr. Biol., 119:187-195, 1991.

[153]

Kunysz A, Glass L, A. Overdrive suppression of spontaneously beating chick heart cell aggregates: experiment and theory. Am. J. Physiol., 269:H1153-H1164, 1995.

[154]

Kunysz A, Shrier A, Glass L. Bursting behaviour during fixed delay stimulation of spontaneously beating chick heart cell aggregates. Am. J. Physiol., 42:C331-C346, 1997.

[155]

Kurrer C, Schulten K. Propagation of chemical waves in discrete excitable media: anisotropic and isotropic wave fronts. in Nonlinear Wave Processes in Excitable Media, Holden A V, Markus M, Othmar H G. eds. (Plenum Press, New York), 489-500, 1991.

[156]

Laurita K R, Girouard S D, Rosenbaum D S. Modulation of ventricular repolarization by a premature stimulus. Role of epicardial dispersion of repolarization kinetics demonstrated by optical mapping of the intact guinea pig heart. Circ. Res., 79(3):493-503, 1996.

[157]

Lechleiter J, Girard S, Peralta E, Clapham D. Spiral calcium wave propataion and annihilation in Xenopus laevis oocytes. Science, 252:123-126, 1991.

[158]

Leon L J, Roberge F A, Vinet A. Simulation of two-dimensional anisotropic cardiac reentry: Effects of wavelength on the reentry characteristics. Ann. Biomed. Eng., 22:592-609, 1994.

[159]

Lewis T, Masters A M. Observations upon conduction in the mammalian heart. A-V conduction. Heart, 12:209-269, 1925.

[160]

Lipp P, Niggli E. Modulation of Ca2+ release in cultured neonatal rat cardiac myocytes. Insight from subcellular release patterns revealed by confocal microscopy. Circ. Res., 74(5):979-990, 1994.

[161]

Lipson M J, Naimi S. Multifocal atrial tachycardia: Clinical associations and significance. Circulation, 67:397-407, 1970.

[162]

Lloyd D. Circadian and ultradian clock-controlled rhythms in unicellular microorganisms. Advances in Microbial Physiology, 39:291-338, 1998.

[163]

E. N. Lorenz Noisy periodicity and reverse bifurcations. in Nonlinear Dynamics, Vol 357 of the Annals of the New York Academy of Sciences (1980).

[164]

Momose-Sato Y, Sato K, Hirota A, Kamino K. GABA-Induced intrinsic light-scattering changes associated with voltage-sensitive dye signals in embryonic brain stem slices: coupling of depolarization and cell shrinkage. Neurophysiol., 79(4):2208-2217, 1998.

[165]

Marban E, Kitikaze M, Koretsune Y, Yue D T, Chacks T, Pike M M. Quantification of [Ca²⁺]i in perfused hearts; Evalulation of the 5F-BAPTA/NMR method as applied to the study of ischemia. and reperfusion. Circ. Res., 66:1255-1287, 1990.

[166]

Marcus M, Hess B. Isotropic cellular automaton for modelling excitable media. Science, 347:56-58, 1990.

[167]

McClintock M K. Menstrual synchorony and suppression. Nature, 229:244-246, 1971.

[168]

Martiel J L, Goldbeter A. A model based on receptor desensitization for cyclic AMP signaling in Dictyostelium cells. Biophys. J., 52:807-828, 1987.

[169]

Mears D, Sheppard N F, Atwater I, Rojas E. Magnitude and modulation of pancreatic beta-cell gap junction electrical conductance in situ. J. Membrane Biol., 146:163-176, 1995.

[170]

Meijler F L, Janse M J. Morphology and electrophysiology of the mammalian atrioventricular node. Physiol. Rev., 68:608-647, 1988.

[171]

Meijler F L, Fisch C. Does the atrioventricular node conduct? Br. Heart J., 61:309-315, 1989.

[172]

Mendez C, Moe G K. Demonstration of a dual A-V nodal conduction system in the isolated rabbit heart. Circ. Res., 19:378-393, 1966.

[173]

Mercader M A, Michaels D C, Jalife J. Reentrant activity in the form of spiral waves in mathematical models of the sinoatrial node. in Cariac Electrophysiology: From Cell to Bedside, 2nd Edition. eds. Zipes D P, Jalife J. W. B Saunders, Philadelphia PA, 1995.

[174]

Minamikawa T, Cody S H, Williams D A. In sitiu visualization of spontaneous calcium waves within perfused rat heart by confocal imaging. Am. J. Physiol, 272(Heart. Circ. Physiol. 41):H236-H234, 1997.

[175]

Mines G R, On the dynamic equilibrium of the heart. J. Physiol. (Lond.), 46:349-383, 1913.

[176]

Mironov S L. Theoretical model of slow-wave membrane potential oscillations in molluscan neurons. Neuroscience, 10:899-905, 1983.

[177]

Misgeld U, Swandulla D. J. Basic Clin. Physiol. Pharmacol., 1:57-63, 1990.

[178]

Moe G K, Preston J B, Burlington H. Physiological evidence for a dual A-V transmission system. Circ. Res., 4:357-375, 1956.

[179]

Moe G K, Rheinbolt W C, Abildskov J A A computer model of atrial fibrilation. Am. Heart. J., 67:200-220, 1964).

[180]

Moe G K, Jalife J, Mueller W J, et al. A mathematical model of parasystole and its application to clinical arrhythmias. Circulation, 56:968-979, 1977.

[181]

Morad M, Dillon S, Weiss X. An acousto-optically steered laser scaning system for measurment of action potential spread in intact heart. in Optical Methods in Cell Physiology, eds Weer, Salzberg (1986).

[182]

Morad M, Trautwein W. The effect of duration of the action potential on contraction in the mammalian heart tissue. Pflugers Arch. ges. Physiol., 299:66-82, 1969.

[183]

Muller T H, Misgeld U, Swandulla D. Ionic currents in cultured rat hypothalamic neurons. J. Physiol. (Lond.), 450:341-362, 1991.

[184]

Munk A A, Adjemian R A, Zhao J, Ogbaghebriel A, Shrier A. Electrophysiological properties of morphologically distinct cells isolated from the rabbit atrioventricular node. J. Physiol. 493(3):801-18, 1996.

[185]

Munk A A, Bub G, Petrecca K, Shrier A. Imaging the spread of excitation on the atrioventricular node of the rabbit heart. Noble Symposium. 1996.

[186]

von Neuman J. Theory of self-reproducing automata. University of Illinois Press., Urbana, 1966.

[187]

Nagumo J, Yoshizawa S, Arimoto S. IEEE Trans. Commun. Technol., 12:400-, 1965.

[188]

Nakagawa N, Kuramoto Y. From collective oscillations to collective chaos in a globally coupled oscillator system. Physica D, 75:74-80, 1994.

[189]

Nakamura Y, Tominaga F, Munakata T. Clustering behavior in time-delayed nearest neighbor coupled oscillators. Phys. Rev. E., 49(6):4849-4856, 1994.

[190]

Noma A, Tsuboi A N. Dependence of junctional conductance on proton, calcium and magnesium ions in cardiac paired cells. J. Physiol., 382:193-211, 1987.

[191]

Paes de Carvalho A, De Almeida D F. Spread of activity through the atrioventricular node. Circ. Res., 8:801-809, 1960.

[192]

Palti Y, Ben David G, Lachov E, Mika Y H, Omri G, Schatzberger R. Islets of Langerhans generate wavelike electrical activity modulated by glucose concentration. Diabetes, 45:595-601 1996.

[193]

Petrecca K, Amellal F, Laird D W, Cohen S A, Shrier A. Sodium channel distribution within the rabbit atrioventricular node as analysed by confocal microscopy. Journal of Physiology., 501(2):263-74, 1997.

[194]

Pieper CF, Parsons D, Lawrie GM, Lacy J, Roberts R, Pacifico A Design and implementation of a new computerized system for intraoperative cardiac mapping. J Appl Physiol., 71(4):1529-1539 1991.

[195]

Plant R E, Kim M. On the mechanism underlying bursting in the Aplasia abdominal ganglion R-15 cell. Math. Biosci., 26:357-375, 1975.

[196]

Plant R E, Kim M. Mathematical description of a bursting pacemaker neuron by modification of the Hodgkin-Huxley equations. Biophys. J., 16:227-244, 1976.

[197]

Prystowsky EN. Atrioventricular node reentry: physiology and radiofrequency ablation. Pacing Clin. Electrophysiol., 20(2 Pt 2):552-571, 1997.

[198]

Oberholzer M, Ostreicher M, Christen H, Bruhhmann M. Methods in quantitative image analysis. Histochem. Cell Biol., 105(5):333-355, 1996.

[199]

Orbach H S, Cohen L B. Optical monitoring of activity from many areas of the in vitro and in vivo salamander olfactory bulb: a new method for studying functional organization in the vertebrate central nervous system. J. Neurosci., 3(11):2251-2262, 1983.

[200]

Onimaru H, Arata A, Homma I. Neuronal mechanisms of respiratory rhythm generation. An approach using in vitro preparation. Japanese J. Physiol., 47(5):385-403, 1997.

[201]

Palsson E, Cox E C. Origin and evolution of circular waves and spirals in Dictyostelium discoideum territories. Proc. Natl. Acad. Sci. USA., 93:1151-1155, 1996.

[202]

Panfilov A V, Holden A V. Self-generation of turbulent vortices in a two-dimensional model of cardiac tissue. Phys. Lett. A., 151(1,2):23-26, 1990.

[203]

Panfilov A V, Holden A V. Computer simulation of re-entry sources in myocardium in two and three dimensions. J. Theor. Biol., 161:271-285, 1993.

[204]

Panfilov A V, Hogeweg P. Spiral breakup in a modified FitzHug-Nagumo model. Phys. Lett. A., 176:295-299, 1993.

[205]

Panfilov A V, Vasiev B N. Vortex initiation in a heterogeneous excitable medium. Physica D, 49:107-113, 1991.

[206]

Pertsov A M, Ermakova E A. Mechanism of drift of a spiral wave in an inhomogeneous medium. Biofizika, 33:338-342, 1988.

[207]

Pizarro G, Cleemann L, Morad M. Optical measurement of voltage-dependent Ca2+ influx in frog heart. Proc. Natl. Acad. Sci. USA., 82(6):1864-1868, 1985.

[208]

E. L. C. Pritchett, E. A. McCarthy, K. L. Lee, W. E. Wilkeinson, in Cardiac Electrophysiology: From Cell to Bedside, D. P. Zipes, J. Jalife, Eds (Philadelphia, W. B. Saunders Co., 1990) pp. 703-707.

[209]

Publicover N. Calcium oscillations in freshly dispersed and cultured interstitial cells from canine colon. Am. J. Physiol., 262:C589-C597, 1992.

[210]

Quan W, Rudy Y. Unidirectional block and reentry of cardiac excitation: A model study Circ. Res., 66:367-382, 1990.

[211]

Ratzlaff E H, Grinvald A. A tandem-lens epiflourescence macroscope : hundred-fold brightness advantage for wide-field imaging. J. Neurosci. Meth., 36:127-137, 1991.

[212]

Rhodes F, Thompson C L. Topologies and rotation numbers for families of monotone functions of the circle. J. London Math. Soc. 43(2):156-170, 1991.

[213]

Rhor S, Salzberg B M. Characterization of impulse propagation at the microscopic level across geometrically defined expansions of excitable tissue: multiple site recording of transmembrane voltage (MSORTV) in patterned growth heart cell cultures. J. Gen. Physiol., 104:287-309, 1994.

[214]

Rhor S, Kucera J P, Kleber A G. Slow conduction in cardiac tissue, I: Effects of reduction of excitability versus a reduction of electrical coupling on microconduction. Circ. Res., 83:781-794, 1998.

[215]

Rinzel J, Lee Y S. Dissection of a model for neural parabolic bursting. J. Math. Biol., 25m:653-675, 1987.

[216]

Roberge F A, Vinet A, Victorri B. Reconstruction of propagated electrical activity with a two-dimensional model of anisotropic heart muscle. Circ. Res., 58:461-475, 1986.

[217]

Rosario L M, Barbosa R M, Antunes C M, Silva A M, Abrunhosa A J, Santos R M Bursting electrical activity in pancreatic B-cells: evidence that the channel underlying the burst is sensitive to Ca²⁺ influz throu L-type Ca²⁺ channels. Pfluegers Arch., 424:439-447, 1993.

[218]

Rosenbaum D S, Kaplan D T, Kanai A, Jackson L, Garan H, Cohen R J, Salama G Repolarization inhomogeneities in ventricular myocardium change dynamically with abrupt cycle length shortening. Circulation, 84(3):1333-1345, 1991.

[219]

Rudy Y, Quan W. A model study of the effects of the discrete cellular structure on electrical propagation in cardiac tissue. Circ. Res., 61:815-823, 1987.

[220]

Salama G, Morad M. Merocyanine 540 as an optical probe of transmembrane electrical activity in the heart. Science, 191(4226):485-487, 1976.

[221]

Salama G, Lombardi R, Elson J. Maps of optical action potentials and NADH fluorescence in intact working hearts. Am. J. Physiol., 252(2 Pt 2):H384-394, 1987.

[222]

Salzberg B M, Grinvald A, Cohen L B, Davila H, Ross W N. Optical recording of neuronal activity in an invertebrate central nervous system: simultaneous monitoring of several neurons. J. Neurophysiol., 40:1281-1291, 1977.

[223]

Salzberg B M, Davila H V, Cohen L B. Optical records of impulses in individual neurons of an invertebrate central nervous system. Nature, 246:508-509, 1973.

[224]

Sano T, Takayama N, Shimamoto T. Directional differece in conduction velocity in cadiac ventricular syncytium studied by microelectrodes. Circ. Res., 7:262-267, 1959.

[225]

Satin L, Cook D. Calcium current inactivation in insulin-secreting cells is mediated by calcium influx and membrane depolarization. Pfluegers Arch., 418:417-422, 1989.

[226]

Saxberg B E H, Cohen R J. Cellular automata models of cardiac conduction. in Theory of Heart: Biomechanics, Biophysics, and Nonlinear Dynamics of Cardiac Function. L Glass, P Hunter and A McCulloch eds. Springer-Verlag New York, (1991).

[227]

Schoels W, Gough W B, Restivo M, El-Sherif, N. Circus movement atrial flutter in the canine sterile pericarditis model. Activation patterns during initiation, termination, and sustained reentry in vivo. Circ. Res., 67:35-50, 1990.

[228]

Selverston A I, Moulins M. Oscillatory neural networks. Ann. Rec. Physiol., 47:29-48 1985.

[229]

Sherman A, Rinzel J, Keizer J. Emergence of organized bursting in clusters of pancreatic B-cells by channel sharing. Biophys. J., 54:411-425, 1988.

[230]

Sherman A. Contributions of modeling to understanding stimulus-secretion coupling in pancreatic B-cells. Am. J. Physiol., 271:E362-E372, 1996.

[231]

Shibata J, Bures J. Optimum topographical conditions for reverberating cortical spreading depression in rats. J. Neurobiol., 5:107-118, 1974.

[232]

Smith P A, Ashcroft F M, Rorsman P Simultaneous recordings of glucose dependent electrical activity and ATP-regulated K⁺- currents in isolated mouse pancreatic B-cells. FEBs Lett., 261:187-190, 1990.

[233]

Smolen P, Keizer J. Slow voltage inactivatio of Ca²⁺ currents and bursting mechanisms for the mouse pancreatic B-cell. J. Membr. Biol., 127:9-19, 1992.

[234]

Smolen P, Rinzel J, Sherman A. Why pancreatic islets burst but single B-cells do not: the heterogeneity hypothesis. Biophys., J. 64:1668-1680, 1993.

[235]

Sole R V, Vallis J Spiral waves, chaos and multiple attractors in lattics models of interacting populations. Physics Letters A, 166:123-128, 1992.

[236]

Spach M S, Miller W T, Geselowitz D B, Barr R C, Kootsey J M, Johnson E A. The discontinuous nature of propagation in normal canine cardiac muscle. Circ. Res., 48:39-54, 1981.

[237]

Spach M S, Heidlage J F. The stochastic nature of cardiac propagation at a microscopic level: electrical description of myocardial architecture and its application to conduction. Circ. Res., 76:366-380, 1995.

[238]

Spring K R. Detectors for fluorescent microscopy. Scanning Microsc., 5(1):63-69, 1991.

[239]

Starbin J M, Zilbeiter Y I, Rusnak E M, Starmer C F. Wavelet formation in excitable cardiac tissue: the role of wavefront-obstacle interactions in initiating high-frequency fibrillatory-like arrhythmias. Biophys. J., 70(2):581-594, 1996.

[240]

Strogatz S H, Mirollo R E. Collective synchronization in lattices of non-linear oscillators with randomness. J. Phys. A: Math. Gen., 21:L699-L705, 1988.

[241]

Strogatz S H, Stewart I. Coupled oscillators and biological synchronization. Scientific American., 269(6):102-109, 1993.

[242]

Sung RJ, Lauer MR, Chun H Atrioventricular node reentry: current concepts and new perspectives. Pacing Clin. Electrophysiol., 17(8):1413-1430, 1994.

[243]

Tasaki I, Watanabe A, Sandlin R, Carnay L. Changes in fluorescence, turbidity, and birefringence associated with nerve excitation. Proc. Natl. Acad. Sci. USA., 61:883-888, 1968.

[244]

Tawara S. Das Reizleitungs System des herzens. Jena, Germany: Fischer, 1906.

[245]

Tsay T, Inman R, Wray B, Herman B, Jacobson K. Characterization of low light level video cameras for fluoresence microscopy. in: Optical Microscopy for Biology, Herman B, Jackobson K (eds.), Wiley-Liss, New York, 219-233. (1990).

[246]

Tyson J J, Keener J P. Singular perturbation theory of travelling waves in excitable media (a review). Physica D., 32:327-361, 1988.

[247]

Valaskovic G A, Morrison G H. Quantitative imaging ion microscopy: a short review. Scanning Microsc., 6(2):305-318, 1992.

[248]

Vasalle M. The relationship among cardiac pacemakers: Overdrive suppression. Circ. Res., 41:269-277, 1977.

[249]

Vinet A, Roberge F A.The dynamics of sustained reentry in a ring model of cardiac tissue. Ann. Biomed. Eng., 22:568-591, 1994.

[250]

Wang H X, de Paola R, Norwood W I, Dynamics underlying the patterning of cardiac dysrhythmais. Phys. Rev. Lett. 70(23):3671 - 3674, 1993.

[251]

Wang X -J, Rinzel J. in The Handbook of Brain Theory and Neural Networks ed. Arbib, M. A., (MIT Press, Cambridge MA), 686-691, 1995.

[252]

Weimar J, Tyson J J, Watson L T. Diffusion and wave propataion in cellular automaton models of excitable media. Physica D., 55:309-327, 1992.

[253]

Weiss J, Cabeen W R Jr, Roberts N K. Nodoventricular accessory atrioventricular connection associated with dual atrioventricular pathways: a case report and review of the literature. J. Electrocardiol., 14(2):185-190, 1981.

[254]

Wiener N, Rosenblueth A. The mathematical formulation of the problem of conduction of impulses in a network of connected excitable elements, specifically in cardiac muscle. Arch. Inst. Cardiol. Mex., 16:205-, 1946.

[255]

Wilders R, Jongsma H, Ginneken A. Pacemaker activity of the rabbit sinoatrial node. A comparison of mathematical models. Biophys. J., 60(5):1202-1215, 1991.

[256]

Williams T L, Sigvart K A, Kopell N, Ermentrout G B, Remler M P. Forcing of coupled nonlinear oscillators: studies of intersegmental coordination in the lamprey locomotor central pattern generator. J. Neurophysiol., 64(3):862-871, 1990.

[257]

Willows A O D. Behavioral acts elicited by stimulation of single, identifiable brain cells. Science, 157:570-574, 1967.

[258]

Wilson W A, Wachtel H. Negative resistane characteristics essential for the maintenance of slow oscillations in bursting neurons. Science, 186:932-934, 1974.

[259]

Schaffer P, Ahammer H, Muller W, Koidl B, Windisch H Di-4-ANEPPS causes photodynamic damage to isolated cardiomyocytes. Pflugers Arch., 426(6):548-551, 1994.

[260]

Windisch H, Ahammer H, Schaffer P, Muller W, Platzer D. Optical multisite monitoring of cell excitation phenomena in isolated cardiomyocytes. Pflugers Arch., 430(4):508-518, 1995.

[261]

Winfree A T. An integrated view of the resetting of a circadian clock. J. Theor. Biol., 28:327-352, 1970.

[262]

Winfree A T. Spiral waves of activity. Science, 175:634-636, 1972.

[263]

Winfree A T. Scroll-shaped waves of chemical activity in three dimensions. Science, 181:937-939, 1973.

[264]

Winfree A T, Strogatz S H. Singular filaments organize chemical waves in three dimensions. I. Geometrically simple waves. Physica D, 8:35-49, 1983.

[265]

A. T. Winfree, When Time Breaks Down: The Three-Dimensional Dynamics of Electrochemical Waves and Cardiac Arrhythmias, (Princeton University Press, Princeton, NJ, 1987).

[266]

Winfree A T. Electrical instability in cardiac muscle: Phase singularities and rotors. J. Theor. Biol., 138:353-405, 1989.

[267]

Winfree A T. Varieties of spiral wave behavior: an experimentalists approach to the theory of excitable media. Chaos, 1:303-334, 1991.

[268]

Winfree A T. Alternative stable rotors in an excitable medium. PhysicaD, 49:125-140, 1991.

[269]

Winfree A T. Estimating the ventricular fibrillation threshold. in Theory of Heart: Biomechanics, Biophysics, and Nonlinear Dynamics of Cardiac Function. L Glass, P Hunter and A McCulloch eds. Springer-Verlag New York, (1991).

[270]

Winfree A T. Heart muscle as a reaction-diffusion medium: the roles of electric potential diffusion, activation front curvature, and anisotropy. Int. J. Bifurcation and Chaos., 7(3):487-526, 1997.

[271]

Wit A L, Cranefield P F, Hoffman B F. Slow conduction and reentry in the ventricular conducting system. Circ. Res., 30:304-322, 1972.

[272]

Witkowski F X, Leon L, Penkoske P, Giles W, Spano M, Ditto W, Winfree A. Spatiotemporal evolution of ventricular fibrillation. Nature, 392(6671):78-82, 1998.

[273]

Witkowski F X, Leon L J, Penkoske P A, Clark R B, Spano M L, Ditto W L, Giles W R. A method for visualization of ventricular fibrillation: design of a cooled fiberoptically coupled image intesified CCD data acquisition system incorporating wavelet shrinkage based on adaptive filtering Chaos, 8(1):1054-1500, 1998.

[274]

Glass L, Zeng W Z Complex bifurcations and chaos in simple theoretical models of cardiac oscillations in Mathematical approaches to cardiac arrhythmias, vol 591 of the Annals of the New York Academy of Sciences, 316-327 (1990).

[275]

Zeng W Z, Glass L, Shrier A. Evolution of rhythms during periodic stimulation of embryonic chick heart cell aggregates. Circ. Res., 69:1022-1033, 1991.

[276]

Zykov V S. Cycloidal circulation of spiral waves in excitable medium. Biofizika, 31:862-865, 1986.

[277]

Zykov V S. Analytical evaluation of the dependence of the speed of an excitation wave in a two-dimensional excitable medium on the curvature of it's front. Biophysics, 25:906-911, 1980.

[278]

The technical specification for the AD645 amplifier and the SSM2024 VCA units are available from http:www.analogdevices.com. The technical specification for the MRD500 are available from http:www.motorola.com, and are also listed in the Electrosonic electronics catalogue.

Index (showing section)

Footnotes:

¹ The dispersion relation (Equation (1.4)) can be rewritten to give the dependence on the AV nodal conduction time on recovery time. This relation is known as the recovery curve. Fatigue would shift the dispersion relation downward, decreasing conduction velocity.

²The `Leading circle reentry' model is commonly associated with having no excitable gap, as activations recur at the maximum rate allowed by the tissue. This conceptual model ignores curvature effects (which may not play a major role in defining rotor period in tissue as excitable as the heart) and dispersion effects. Dispersion effects predict that speed is inversely related to period. As there should be a gradient in wavefront speed near the reentry period[270], there should be a partially excitable gap.

³Spontaneous activity occurs at unphysiologically low concentrations of extracellular K⁺[109].

⁴Contraction in the monolayer is visually detectable at magnifications of 20x or greater. It is not observable at low magnification.

⁵The monolayers are homogeneous in the sense that they do not contain connective tissue. They contain heterogeneities in the form of uneven cell spacing.

⁶It is noteworthy that recently a camera has been developed by Fuji that stores only the change in light for each pixel.

⁷Hamamatsu has recently released a photodiode array with amplification circuitry, however the gain of the system is too low for measurements of low light levels and is not customizable.

⁸The rate and resolution depends on what convention is used. For example, a typical video system (using the Radio Electronics Television Manufacturing Association scanning convention) scans in 525 lines per frame at 30 frames per second for a scan rate of 15750 Hz.

⁹Software has previously been developed for generating isochronal maps from electrode arrays[194]

¹⁰ The portable greymap (PGM) format is part of the Extended Portable Bitmap format. Information on this format can be obtained by viewing the manual pages on PGMs on all Unix (and Unix like) systems.

¹¹A conceptually simpler way to think of RMS is that the probability is 68% that the voltage signal will be between -RMS and +RMS.

¹²t = CR, where C is capacitance and R is resistance. Here R=10 gigaohm and t=0.05s

¹³If a high resolution digitization scheme is used, offset and signal can be recorded without loss of data. For example, a 12 bit digitization scheme could potentially record a useful signal with 400 (4096 x 10 %) point resolution.

¹⁴Recently, the SSM2024 has been discontinued. Analog Devices supplies a replacement part with better specifications.

¹⁵machining was done by S. Kekani, Physics Department

¹⁶The system used a 1× Ziess microscope objective on a inverted microscope with a Hammmamatsu intensified high speed imaging system. See [209] for system details.

¹⁷This may not always be true. Michael Guevara (personal communication) claims that he has observed bursting dynamics in (assumed isopotential) aggregates.

¹⁸Distances measured in either the x or y directions from a given grid point are measured parallel to the orientation of the grid, while distances measured in both the x and y directions are measured diagonal to the orientation of the grid.

¹⁹Although a factor of two appears large as a comparison between experimental and theoretical values, it is in good agreement when compared to similar scalings for cardiac muscle, which are up to a factor of six too high.