## Abstract

Cardiac and uterine muscle cells and tissue can be either autorhythmic or excitable. These behaviours exchange stability at bifurcations produced by changes in parameters, which if spatially localized can produce an ectopic pacemaking focus. The effects of these parameters on cell dynamics have been identified and quantified using continuation algorithms and by numerical solutions of virtual cells. The ability of a compact pacemaker to drive the surrounding excitable tissues depends on both the size of the pacemaker and the strength of electrotonic coupling between cells within, between, and outside the pacemaking region.

We investigate an ectopic pacemaker surrounded by normal excitable tissue. Cell–cell coupling is simulated by the diffusion coefficient for voltage. For uniformly coupled tissues, the behaviour of the hybrid tissue can take one of the three forms: (i) the surrounding tissue electrotonically suppresses the pacemaker; (ii) depressed rate oscillatory activity in the pacemaker but no propagation; and (iii) pacemaker driving propagations into the excitable region.

However, real tissues are heterogeneous with spatial changes in cell–cell coupling. In the gravid uterus during early pregnancy, cells are weakly coupled, with the cell–cell coupling increasing during late pregnancy, allowing synchronous contractions during labour. These effects are investigated for a caricature uterine tissue by allowing both excitability and diffusion coefficient to vary stochastically with space, and for cardiac tissues by spatial gradients in the diffusion coefficient.

## 1. Introduction

An excitable cell has a stable resting potential, and only generates action potentials (APs) in response to an input or under abnormal conditions. Pacemaking cells are autorhythmic in the absence of an input and in physiological conditions, and can respond to inputs and changes in conditions by a change in the rate or pattern of their periodicity. The qualitative differences between these two types of cell behaviours—excitability or autorhythmicity—result from the different parameters and kinetics of their excitation equations, all of which have the same form. Ionic current flow through the membrane charges or discharges the membrane capacitance, and for a unit area of cell membrane(1.1)where *C*_{m} is the specific membrane capacitance, *V* is membrane potential and the membrane current density, which is the sum of the current densities for the different ion-selective channels, pumps and exchangers. Differences in APs from different tissues, and different parts of the same tissue, result from differences in the expression of the different ion transport and sequestration proteins.

A bifurcation occurs in a dynamical system, such as equation (1.1), as the stability of a solution changes under variation of parameters, and is observable as a qualitative change in behaviour. The bifurcation is located at the parameter value at which a change in the stability of a solution occurs. An example is the emergence of periodic solutions at the point where an equilibrium solution loses its stability. For the Hodgkin–Huxley equations (Hodgkin & Huxley 1952), a subcritical Hopf bifurcation, into small amplitude, unstable periodic solutions, separates the stable resting excitable state from autorhythmicity, induced by an increase in extracellular calcium (Ca^{2+}) concentration or a reduction in membrane maximal potassium (K^{+}) conductance (Holden & Yoda 1981). These unstable solutions connect to large amplitude, stable periodic solutions that correspond to a periodic discharge of AP.

Cardiac sinoatrial pacemaker cells are autorhythmic (Zhang *et al*. 2000), while ventricular cells are excitable (Kleber & Rudy 2004), and normally require repetitive current flows from neighbouring excited cells to generate repetitive APs. However, under abnormal conditions, such as a reduction in the time-independent potassium current conductance *g*_{K1} (Miake *et al*. 2002, 2003) or intracellular sodium concentration ([Na^{+}]_{i}) overload (Harrison *et al*. 1992), periodic activity can be induced in ventricular cells. The difference in behaviour between excitable and autorhythmic cells may be considered dynamically as a bifurcation occurring in parameter space, rather than mechanistically in terms of the ionic currents that contribute towards the pacemaker potential. In §3, we characterize the bifurcations into periodicity produced in virtual ventricular cells, brought about by: (i) block of the pump current , as seen during ischemia (Carmeliet 1999) and digitalis intoxication (Kass *et al*. 1978; Gheorghiade *et al*. 2004) for example, and (ii) downregulation of the time-independent inward rectifying current , which helps to maintain the resting membrane potential in ventricular myocytes and suppresses autorhythmicity. Downregulation of has been proposed as a possible means of clinically inducing pacemaker activity in the ventricles by gene transfer, as an alternative to ventricular pacing by implanted electronic pacemaker (Miake *et al*. 2002). For autorhythmic activity in a compact cluster of pacemaking cells—such as the sinoatrial node or an ectopic focus—to drive surrounding tissues, the cluster needs to be large enough to generate sufficient current and appropriately coupled to the surrounding tissue. In §4, we examine effects of the width of a boundary zone between the autorhythmic and excitable tissues, and the coupling between them. For gravid uterine muscle, there are changes in both the excitability of cells and the intercellular coupling between cells during late pregnancy. In §5, we examine these effects on the synchronization of electrical activity in simplified tissue models.

## 2. Numerical methods

### (a) Cardiac cells and tissues

The rate of change of *V* for a single ventricular myocyte is modelled using equation (1.1). was provided by the Luo–Rudy dynamic model (LRd00; Luo & Rudy 1994) or by the human ventricular model (HVM) of ten Tusscher *et al*. (2004). Changes in ionic concentrations are also modelled:(2.1)where [B] is the concentration of ion B, is the sum of the currents carrying ion B, is capacitive membrane area, is the volume of the compartment whose concentration is being updated, is the valence of ion B and *F* is Faraday's constant. Source codes written in C/C++ can be found at http://www.cwru.edu/med/CBRTC/LRdOnline/LRdModel.htm for LRd00, and at http://www-binf.bio.uu.nl/khwjtuss/HVM for HVM. For both models, we solved equations of the form (1.1) and (2.1) using a simple forward Euler method with a time-step of Δ*t*=0.01 ms.

The bifurcations between different qualitative types of behaviour of nonlinear dynamical systems can be studied in state space as one or more parameters in a dynamical system are varied smoothly. These can be mapped in parameter space by path-following or continuation algorithms that track bifurcations as lines in two-dimensional parameter space (Ermentrout 2002), e.g. auto86 (Doedel & Kernévez 1986; see http://sourceforge.net/projects/auto2000/). These packages reduce the identification, location and tracking of equilibrium and periodic solution branches and other bifurcation problems to the continuation of implicitly defined curves, using predictor–corrector methods. For simplicity, consider autonomous dynamical systems and one-parameter bifurcation diagrams that contain equilibria, periodic solutions and Hopf and homoclinic bifurcations. Once a point on a particular branch is located, a continuation package follows the locus of either equilibrium or periodic states on the branch. While following an equilibrium branch, tracking the real parts (and signs) of the eigenvalues gives the stability of the branch that is exchanged between stable and unstable as the real part passes through zero. If a single complex conjugate pair of eigenvalues crosses the imaginary axis, small amplitude (stable or unstable) oscillations emerge at a Hopf bifurcation. Once the bifurcation point is found, continuation can be switched from the equilibrium to the periodic branch. The stability of periodic oscillations is followed using Floquet multipliers.

An autonomous dynamical system with one parameter of interest may be written as(2.2)where **x** is an *n*-dimensional vector of real numbers, and *μ* is a parameter. **x** may denote the vector of variables in the system of ordinary differential equations (ODEs) defining a cardiac cell model. *μ* is the parameter whose variation causes a qualitative change in the behaviour of the model, such as a membrane maximal conductance. Equilibria are found from the solution of(2.3)The Jacobian of the above vector is defined as an matrix, . If all the eigenvalues of **J** have negative real parts, then the steady state is stable. If an eigenvalue has a positive real part, then the steady state is unstable, at least along a sub-manifold. If the equilibrium is a saddle point (i.e. eigenvalues with positive and negative real parts coexist), and the stable and unstable manifolds overlap, then it is associated with homoclinic orbits. Such orbits start and end at the saddle point and are characterized by infinite periods. While computing these orbits, we have to follow the periodic branches emerging from Hopf bifurcation points in one-parameter bifurcation diagrams. Homoclinic bifurcations occur when homoclinic orbits separate phase space regions of periodic orbits and non-periodic orbits. Although homoclinic bifurcations are global and cannot be detected by local continuation, they can be studied numerically by the methods given in Champneys *et al*. (1996).

Since excitation equations are constructed from a series of voltage and patch clamp experiments under different conditions, although they may successfully reproduce APs, they need not be electrochemically neutral, and so may not have physiologically meaningful equilibrium solutions. Further, they may not be posed in the form of dynamical systems: for example, in LRd00 the intracellular Ca^{2+} is triggered a delay after the peak of the AP. Continuation algorithms may only be applicable to reduced parts of the full excitation system, and so may have to be combined with numerical integration of the full system.

For studying propagation in cardiac tissues, we model one-dimensional tissue strands using the equation(2.4)where ∇ is a spatial gradient operator and *D* a diffusion coefficient. No-flux boundary conditions were imposed at the ends of the strands. A time-step of Δ*t*=0.02 ms and a space step of Δ*x*=0.2 mm were used in one-dimensional simulations.

### (b) Uterine tissue

We simulate excitation in two-dimensional uterine tissue using a modified form of the FitzHugh–Nagumo model (FitzHugh 1961; Nagumo *et al*. 1962):(2.5)(2.6)where *u* is the excitation variable, *v* is the recovery variable, *I* is an input current controlling excitation, *x* and *y* denote space, and the constants *k*=8, *u*_{a}=0.1, *ϵ*=0.02 and *G*=5 are membrane kinetic parameters. Behaviour of the tissue was examined as *I* and *D* were varied over the range 0–0.5. The two-dimensional tissue was simulated as a 300×400 node Cartesian grid with a space step of Δ*x*=Δ*y*=0.4 space units, and was solved using a simple forward Euler step method with a time-step of Δ*t*=0.01 time units.

## 3. Bifurcations in ventricular cells

### (a) I_{NaK} block and intracellular Ca^{2+} oscillations

Block of was examined in the endocardial LRd00 model. Kass *et al*. (1978) proposed that block can induce ectopic pacemaking activity by increasing [Na^{+}]_{i} which in turn causes [Ca^{2+}]_{i} overload through the action of the exchanger. This can result in [Ca^{2+}]_{i} oscillations as Ca^{2+} is spontaneously released from, and is subsequently sequestered back into, the sarcoplasmic reticulum. These [Ca^{2+}]_{i} oscillations can interact with membrane currents, which in turn affect *V* and induce pacemaker activity. As the use of continuation algorithms is not feasible for LRd00 (§2*a*), numerical solutions were obtained in order to identify bifurcations that lead to autorhythmicity. We examined the emergence of [Ca^{2+}]_{i} oscillations during block and the subsequent [Na^{+}]_{i} overload by clamping *V* and all Na^{+} and K^{+} concentrations, and so reducing the full LRd00 system to those equations describing Ca^{2+} handling. Behaviour of the Ca^{2+} handling equations was determined by numerical integration over a period of 120 s, where [Ca^{2+}]_{i} either settled to a stationary stable state or oscillated. With *V* clamped at −90 mV and as the primary bifurcation parameter [Na^{+}]_{i} is increased to approximately 16.1 mM, large period [Ca^{2+}]_{i} oscillations emerge (figure 1*a*). The period of the oscillations decreases rapidly as [Na^{+}]_{i} is further increased (figure 1*b*), indicative of a homoclinic bifurcation rather than the Hopf bifurcation identified by Varghese & Winslow (1993) in the calcium subsystem of the DiFrancesco–Noble Purkinje fibre model (DiFrancesco & Noble 1985). Changes to the amplitude and period of the [Ca^{2+}]_{i} oscillations also occur with changes in the secondary bifurcation parameter *V* (figure 1*c*). With an increase in clamped *V* from −90 to −70 mV (and with [Na^{+}]_{i} clamped at 20 mM), increases in diastolic and systolic [Ca^{2+}]_{i} are observed along with a decrease in the period of the oscillations. The effects of changing both the primary and secondary bifurcation parameters on the [Ca^{2+}]_{i} subsystem dynamics are shown in figure 1*d*.

In the LRd00 cell model with and included, and parameters modified to induce [Ca^{2+}]_{i} oscillations—constant (Varghese & Winslow 1994), block, , [CSQN]_{th}=7.0 mM and the time constants of activation and inactivation of increased to 5 ms—the oscillations seen in the Ca^{2+} handling system interact with the exchanger current and the non-specific Ca^{2+}-activated current to cause pacemaker activity. With [Na^{+}]_{i} set to 20 mM, autorhythmic APs seen during the first 20 s of integration (figure 2*a*) become stable *V* oscillations between −13 and −21 mV. The corresponding [Ca^{2+}]_{i} oscillations (figure 2*b*) that drive these *V* oscillations take approximately 40 s to become stable. This transient behaviour is shown in the phase portraits of figure 2*c*,*d*; figure 2*c* shows the full 120 s of integration including the transient behaviour, figure 2*d* shows the period 60–120 s, when *V* and [Ca^{2+}]_{i} oscillations have become stable.

### (b) *I*_{K1} downregulation

Both LRd00 and the epicardial HVM were used to identify bifurcations and autorhythmicity due to downregulation of . Clinically, mutations in a subunit of the channel are implicated in long QT syndrome 7, which has been associated with a loss of function (Tristani-Firouzi *et al*. 2002), and experimentally, downregulation of has been shown to induce autorhythmicity in ventricular cells (Miake *et al*. 2002, 2003). Downregulation of was modelled by multiplying by a fractional term *x*, where represents downregulation.

In both systems, autorhythmic APs were seen when was completely blocked (fractional ); the membrane potential remained stable when was not suppressed (fractional ). Between these two extremes, a bifurcation occurred. In LRd00 (figure 3*a*), the bifurcation points of different ventricular cell types, although different, are very close (fractional ). However, for a given value of fractional below the bifurcation point, oscillations in M cells have the fastest rate, followed by endocardial then epicardial cells; channel densities determine the location of the bifurcation in parameter space as well as the rate of the oscillating system. As they have the fastest rate, any M cells within an ectopic focus will drive the pacemaking. HVM showed similar behaviour but with smaller bifurcation values in both the full system (figure 3*b*, fractional ) and a reduced system with no ionic concentration dynamics (figure 3*c*, fractional ). The different bifurcation values of the two HVM systems shows that, in addition to channel densities, ionic concentrations play a part in determining the location of the bifurcation and the emergent behaviour. For downregulation to be an effective therapy, a ‘reasonable’ downregulation needs to be effective: comparing these results to guinea-pig experiments (Miake *et al*. 2002, 2003) and computational models (Silva & Rudy 2003), a greater reduction of is required for pacemaker activity to be induced in human ventricular cells.

All systems showed a rapid decrease in the period of the oscillations immediately after the bifurcation. This numerical evidence suggests homoclinic bifurcations. If pacemaking is to be induced by progressively developing downregulation of , the induced pacemaker would start with an unreasonably long period; this would be suppressed by the electrical driving of the residual normal pacemaking and is potentially arrhythmogenic.

## 4. Pacemaker driving of surrounding tissue

### (a) Effects of a border zone between pacemaking and excitable tissues

The effects of pacemaker size on the ability of that pacemaker to drive surrounding tissue have been examined for various cardiac tissues (e.g. Fozzard & Schoenberg 1972; Lindemans & Vandergon 1978; Winslow *et al*. 1993; Joyner *et al*. 2000; Wilders *et al*. 2000; Carey *et al*. 2001). These studies assumed that the border between the two tissues was distinct, and that parameter values inducing autorhythmicity changed in a stepwise fashion between the two tissues. We include a border zone between the pacemaking tissue and the normal excitable tissue by adding an area between the two tissues where parameter values and initial conditions change smoothly in space. Using block of to induce pacemaking (§3*a*) located at one end of a 15 mm LRd00 endocardial one-dimensional strand, the parameters are changed either linearly or sigmoidally between the two tissues. For sigmoidal changes, parameters are scaled by , where *a* is a constant and *x* denotes space within the border zone. Figure 4 shows the effects of these different spatial changes of parameter values within the border zone on the liminal length of pacemaker tissue required to drive propagation. In all cases, the absolute changes in parameter values between the pacemaking and excitable tissues are the same: only the spatial rate of change of parameter values in the border zone varies. When the spatial change is linear, liminal pacemaker length decreases from 3.4 to 3.0 mm as the border zone is increased in size from 0 to 3.0 mm. Further increases in border zone size do not cause a further reduction in liminal pacemaker length. With sigmoidal changes of parameter values within the border zone where *a*=0.5, 1 or 2, the liminal pacemaker length is gradually decreased as border zone length increases. When the border zone length is larger than 11.6, 9.0 and 7.6 mm for *a*=0.5, 1 and 2, respectively, the border zone itself can act as a pacemaker and so no central pacemaker tissue is required for propagation. The differences between the effects of linear and the different sigmoidal changes within the border zone suggest that the rate of change of the parameter values with respect to space is more relevant than the absolute change in parameter value. When the spatial rate of change of parameter values in the border zone near the pacemaking tissue is low (e.g. with a sigmoidal change where *a*=2), more border zone tissue is raised above the threshold at which it becomes autorhythmic. If the rate of change is high (e.g. with a linear change), less border zone tissue is above the autorhythmic threshold, therefore increasing the length of pacemaker tissue required to drive propagation.

### (b) Effects of spatial changes in the diffusion coefficient

Propagation through excitable media, such as cardiac and uterine muscle, depends on the diffusion coefficient, *D*, that is used to model current flow through intercellular gap junctions. The examples of propagation along a one-dimensional strand presented in §4*a* have all been in uniformly coupled tissues, i.e. *D* is spatially homogeneous. However, the coupling in real tissues is heterogeneous, e.g. cell–cell uncoupling during ischemia (Carmeliet 1999) or transmural differences (Yan *et al*. 1998; Gima & Rudy 2002). We altered intercellular coupling by changing *D* within pacemaking tissue and within a border zone separating the pacemaker from the excitable tissue. Figure 5*a* shows the effects of changing the length of the pacemaking tissue (tissue with block, §3*a*) and the diffusion coefficient within the pacemaker, *D*_{P}, located at one end of a 15 mm endocardial LRd00 strand with no border zone. Three types of behaviour are observed. At low pacemaker lengths, the electrotonic influence of the excitable tissue suppresses the pacemaker activity, resulting in sub-threshold *V* oscillations within the pacemaking tissue. For values of *D*_{P} approximately larger than 0.018 mm^{2} ms^{−1} and for pacemaker lengths larger than the liminal length for that particular tissue, the voltage oscillations in the pacemaker tissue reach threshold, resulting in APs that propagate into the excitable tissue and along the strand. For values of *D*_{P} approximately smaller than 0.018 mm^{2} ms^{−1}, however, propagation is not achieved. The electrotonic influence of neighbouring tissues is minimal and the tissue behaves as a series of uncoupled individual cells rather than a functional syncytium. Here, there are APs within the pacemaker but, due to the low coupling, no propagation into the excitable tissue. These three types of behaviour are illustrated in the space–time plots of figure 5*b–d*. Figure 5*b* shows a pacemaking tissue of length 2 mm with *D*_{P}=0.05 mm^{2} ms^{−1} in which pacemaker activity is suppressed. In figure 5*c*, the pacemaker length (4 mm) is supra-liminal for the diffusion coefficient *D*_{P}=0.05 mm^{2} ms^{−1} and so propagation occurs along the strand. In figure 5*d*, *D*_{P} is reduced to 0.01 mm^{2} ms^{−1} with the pacemaker length at 4 mm: periodic APs are observed in the pacemaker but they do not propagate into the excitable tissue. Qualitatively similar results are found when either a linearly changing or sigmoidally changing (*a*=1) 2 mm border zone is introduced between the pacemaker and excitable tissues (results not shown). Similar behaviour was observed by Wilders *et al*. (2000) in a two-dimensional sheet of ventricular cells where a single sinoatrial node model acted as the pacemaker, with the coupling currents associated with this node scaled in order to simulate variations in the size of the pacemaker.

## 5. Synchronization and bursting

During early pregnancy, local contractions of the uterine muscle occur, and as pregnancy progresses, the amplitude of contractions increases (Wray 1993). In advanced pregnancy, the contractions involve most of the uterine tissue. There is an increase in activity and in synchronization of activity during late pregnancy: the increase in synchronization correlates with an increase in the density of gap junctions (Garfield *et al*. 1995). In terms of dynamics, there is an increase in activity, from a predominantly resting to a predominantly periodic state, and an increase in coupling. We caricature these patterns of activity by a two-dimensional excitable medium, in which activity and coupling were distributed randomly with FitzHugh–Nagumo kinetics (equations (2.5) and (2.6)). As there are no anatomically defined pacemaker sites in the uterus, we allowed the two parameters *I* and *D* to vary randomly throughout the medium. At each gridpoint in the medium, *I* and *D* had a fixed value, selected from a Gaussian distribution around a given mean , in the range from 0 to 0.5, with standard deviation of 0.15. The parameter *I* is a bifurcation parameter for the FitzHugh–Nagumo ordinary differential system, and for our parameter set a Hopf bifurcation into periodic, AP-like oscillations occurs at *I* close to 0.25. The equilibrium solution regains its stability at *I* close to 2.2, and so, within the range used, the parameter *I* specifies whether a node is excitable or autorhythmic and, if autorhythmic, its rate. The diffusion coefficient *D* was used to simulate the cell–cell coupling. The qualitative behaviour of the system is presented in parameter space of and , where qualitatively different regions were identified, as illustrated in figure 6. We constructed the parameter space with either only *I* or *D*, and with both *I* and *D* randomly distributed over the two-dimensional model. Four types of behaviour were classified from space–time movies:

Activity is spatially sparse, leading to a low level of spatially integrated activity. At any one active node, the size of the AP is reduced, as a result of electrotonic spread to surrounding tissue.

Spatial clustered activity, each cluster acting as a localized focus emitting periodic waves. Activity at any one node is burst-like, alternating between different rates of repetitive APs, or between repetitive activity and silence.

Moderate level of almost spatially uniform excitation: at any one node there is a maintained depolarization, or small amplitude oscillations around a maintained depolarization.

Predominantly spatially coherent synchronous discharge of APs.

For low values of parameters and *D* or , activity was spatially localized and asynchronous. As the value of was increased, the activity increased, in frequency and integrated amplitude. With an increase in the value of *D* or , there was an increase in spatial coherence, i.e. larger space scales and temporal synchronization of activity was observed.

## 6. Conclusions

The differences between normal pacemaking and excitable cells can be understood in terms of the stability of their equilibrium states: an excitable tissue, such as normal ventricular cells, has a stable resting potential that corresponds to a stable equilibrium, while a pacemaking cell has a stable periodic state. These states are separated in parameter space by bifurcations. For ODE models of the electrical activity of single cells, if they are electrochemically neutral with physiologically sensible equilibrium states, bifurcations can be mapped in parameter space using continuation algorithms. If the cell model is not electrochemically neutral, and so does not have a stable resting state, or if the cell model is not formulated as a differential system, bifurcations can still be mapped by a combination of numerical integrations and the application of continuation algorithms to reduced systems.

Automaticity can be induced in ventricular cell models, by reducing a repolarizing conductance, or by blocking an ion-exchange process. In both cases, automaticity is produced by a homoclinic bifurcation—the emergence of very long-period oscillations—rather than by a Hopf bifurcation. If such automaticity is to act as an ectopic pacemaker, driving the surrounding tissue, the effectiveness and minimal size of the pacemaking tissue depends on the coupling within and between the pacemaker and the excitable tissue.

In virtual uterine tissue, spatial heterogeneity is seen to assist transitions between quiescence, bursting, synchrony and quenched systems. A spatial heterogeneity in intercellular coupling on its own does not permit large-scale synchronous behaviour within the parameter range investigated. Simultaneous heterogeneity in excitation and intercellular coupling, however, increases the possibility of bursting.

## Editors' note

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

## Acknowledgements

This work is supported by grants from the MRC (G0000315), the EPSRC (GR/S43498/01, GR/R92592/01, GR/R40838/01) and the EU Network of Excellence ‘BioSim’ (005137). A.P.B. is supported by an MRC bio-informatics (computational biology) priority area research studentship (G74/63), W.C.T. by a BHF research studentship (FS/03/075/15914).

## Footnotes

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

- © 2006 The Royal Society