The physics of heart rhythm disorders

2.1. Basic anatomy and activation sequence

The anatomy of the heart is schematically shown in Fig. 1A and is described in great detail in numerous textbooks (see, e.g., [5]). Briefly, it consists of two small chambers, the atria, and two larger ones, the ventricles. The walls of these chambers are comprised of different cells, including cardiac myocytes, a special type of muscle cells that allow these chambers to contract and pump blood. The blood flow in the body is described in detail in the literature and will not be the focus of this review [8]. Suffices to say that deoxygenated blood flows from the body to the right atrium (RA; commonly used abbreviations are listed in Table 1), after which it is pumped into the right ventricle (RV). Contraction of this ventricle pumps the blood to the lungs where it is oxygenated and, subsequently, returned to the left atrium, which transports it to the left ventricle. Left ventricle contraction allows this oxygen-rich blood to flow to the rest of the body. Finally, valves between the different chambers ensures unidirectional blood flow.

Open in a separate window

Figure 1:

Schematic representations of the heart. A. The relevant anatomical parts with arrows indicating direction of blood flow: deoxygenated blood enters the right atrium through the venae cavae and is pumped by the right ventricle to the lungs through the pulmonary artery. Oxygenated blood returns to the heart through the pulmonary veins and is circulated to the rest of the body through the aorta. B. Schematic representation of the electric system of the heart, shown here in blue. 1. Sinoatrial node 2. Atrioventricular node 3. Bundle of His 4. Left bundle branch 5. Left posterior fascicle 6. Left-anterior fascicle 7. Left ventricle 8. Ventricular septum 9. Right ventricle 10. Right bundle branch (By J. Heuser - self made, based upon Image:Heart anterior view coronal section.jpg by Patrick J. Lynch).

Clearly, the contraction of these chambers needs to be carefully orchestrated, as a disruption of the blood flow will result in inadequate oxygen supply to the body. This coordination is made possible by an electrical wave that coherently and periodically sweeps through the tissue of the heart and that, ultimately, results in cell contraction. During the resting phase of the heart, called diastole, cardiac cells are hyperpolarized, with a resting potential of around −85 mV. Just as muscle cells and neurons, however, myocytes can depolarize and, due to inter-cellular coupling, an electrical wave can propagate through the tissue.

This wave propagates through the electrical conduction system, schematically shown in Fig. 1B. It is initiated by a group of cells within the right atrium, called the sino-atrial node [9], which shows oscillatory behavior and spontaneously depolarize in a periodic fashion. This group of cells regulates the frequency of one’s heartbeat, which is a function of the person’s physical and mental state: physical exercise, fear, nervousness, and other conditions all lead to an increase in the frequency of sino-atrial firing while relaxation and calming exercises can lead to a decrease [10].

Once fired, the electrical activation propagates through the neighboring cells, made possible by the electrical coupling of the cells and the excitable nature of the cells, further explained in Section 2.2 [11]. This coupling is mediated by gap junctions, assemblies of channel proteins that allow rapid exchange of ions between cells [12]. The electrical wave sweeps through the right atrium and, immediately following that, the left atrium. The speed of propagation in normal atrial tissue is about 100 cm/s which, together with a typical size of the atria, results in atrial depolarization within less than 90 ms. These useful numbers of the heart, along with others, are listed in Table 2. As a consequence of the cells depolarization, a large amount of calcium (Ca²⁺) rushes into the cell. These calcium ions bind to troponin which, ultimately, leads to cell contraction [13].

The contraction of the atria is slightly earlier than the contraction of the ventricles, which ensures effective pumping of the blood from the small to the large chambers. This is possible since the electrical wave in the atria does not immediately propagate further into the ventricles. In fact, the atria are electrically isolated from the ventricles by a non-conducting barrier, except for a specialized region of cells called the atrioventricular (AV) node (see Fig. 1B). Conduction velocity within the AV node is only 5 cm/s, resulting in a delay of propagation. This delay results in complete atrial depolarization and contraction before the ventricles are activated. Following activation of the AV node, the electrical signal propagates through specialized sets of fibers (bundle of His and the Purkinje fibers) that run from the atria to the endocardial surface of the ventricles and conduct the activation from the atria with minimal time delay. The resulting ventricular depolarization starts at the inner surface and the base of the heart, followed by transmural and base-apex propagation with a conduction velocity around 50 cm/s. The entire heart mass is depolarized within approximately 250 ms, after which the cells slowly return to their resting potential.

2.2. Single cell anatomy and electrophysiology

A single cardiac myocyte is roughly cylindrical in shape, with a length of approximately 120 μm and a diameter of 25 μm. Ventricular cells contain transverse tubules, or T-tubules, narrow invaginations of the sarcolemma that occur at so-called Z-lines and are regularly spaced (~ 1.8μm) [18]. This striation can be clearly seen in the confocal micrograph of a human cardiac myocyte in Fig. 2A, where α-actinin, a protein that anchors actin filaments to the Z-lines, is labeled in green. Many of the proteins responsible for Ca²⁺ handling, and thus contraction, are concentrated within these tubules [13]. It is of note, however, that atrial cells do not have T-tubules [18]. Embedded with the membrane of the myocyte are a large number of ion channels, as schematically shown in Fig. 2B. These channels consist of specialized proteins that are selectively permeable to certain ions. For example, the potassium channel has a conductivity for K⁺ that is much larger than for the similarly sized sodium ion Na⁺. How this is achieved can be addressed by detailed molecular dynamics simulations, reviewed in [19]. The precise working of these channels have been reviewed before and we will only briefly describe their function [20, 21].

Open in a separate window

Figure 2:

Electrophysiology of a myocyte. A. Fluorescent image of a human cardiac myocyte where α-actinin is labeled in green and the nucleus is stained in blue. Scale bar: 10 μm. (From [24]). B. Schematic representation of a ventricular cardiac myocyte and some of its ion channels, ion transporters, signaling cascades, and structure. RyR indicates ryanodine receptors; SR, sarcoplasmic reticulum; Mito, mitochondria; CaMK, calcium-calmodulindependent protein kinase; NCX, Na/Ca exchanger; PLM, phospholemman; PLB, phospholamban; ATP, ATPase; PKA, protein kinase A; β-AR, β-adrenergic receptor; IKr, IKs, IK1, potassium channels; and ICl(Ca), a calcium-activated Cl current. (From [25]) C. Typical action potential shape, along with the most important currents during each phase of the action potential.

At rest, the membrane potential V of cardiac cells is roughly −85 mV, close to the Nernst potential for potassium. At this potential, the net flux of potassium ions is zero and the chemical and electrical forces are in balance [22]. Almost all cardiac cells are excitable. This means that when the membrane potential is raised above a critical value, the potential rapidly increases, followed by a gradual return to equilibrium. In other words, a small stimulus can result in a large response. As a consequence of this excitability, cardiac cells also exhibit a refractory period. That is to say, a cell cannot be excited again unless it has recovered sufficiently. In Section 4 we will see that this refractory period plays an important role in the initiation of arrhythmias. Of course, cardiac cells are not only the only excitable cells in the human body. Neurons are also excitable and display a similar threshold and refractoriness. A significant difference between cardiac cells and neurons, however, is the APD. For neurons, this AP is of the order of 1 ms [23] while in cardiac cells the AP can take on values of 100-300 ms.

In cardiac cells, an increase in the potential results in the opening of sodium channels, which are controlled by the voltage, and sodium ions rush in (I_Na in Fig. 2C). As a result, the potential of the cell increases rapidly, leading to a large and fast upstroke. After several milliseconds, the potential reaches its maximum, after which it remains at a plateau for 50-200 ms during which the influx of calcium ions (blue curve in Fig. 2C) is electrically balanced by the efflux of potassium ions. Once the Ca²⁺ channels close, the potassium currents (I_K in Fig. 2C) return the potential back to its resting state. The ensuing dynamics of the voltage is called the action potential (AP) while the action potential duration (APD) is defined as the time spent by the cell above a certain percentage value of repolarization (often taken to be 90%).

As already mentioned, Ca²⁺ is essential in the coupling between electrical excitation and mechanical contraction, carried out by the motor protein myosin that pulls actin filaments in a ratcheted fashion [13]. Before depolarization, the intracellular (cytosolic) Ca²⁺ concentration is very low, only roughly 100 nM. As a result, the myosin binding site on actin is blocked, preventing contraction. During repolarization, the cytosolic Ca²⁺ concentration increases sharply, and Ca²⁺ binds to the regulatory protein troponin, which exposes the myosin binding site on actin [26]. The subsequent binding of myosin to actin, combined with the energy released from ATP hydrolysis, creates a power stroke, which provides the mechanical force required for muscle contraction [27].

The sharp increase in Ca²⁺ required for contraction is achieved through a positive feedback mechanism called Calcium Induced Calcium Release (CICR) [28]. During depolarization, calcium ions that move into the cell through the voltage gated L-type Ca²⁺ channels (LCCs) enter a small sub-cellular compartment called the dyadic cleft, which separates the membrane (containing the LCCs) and the sacroplasmic reticulum (SR). These SRs are membrane-bound structures that function as Ca²⁺ storage units. Calcium ions within the cleft then bind to and activate the cardiac ryanodine receptors (RyRs) that are embedded in the SR membrane (see Fig. 2B). This activation leads to a large influx of Ca²⁺ into the cytoplasm and a ~ 10-fold increase in cytosolic Ca²⁺ concentration [29, 30]. Also part of the SR are SERCAs, Ca²⁺ pumps that are responsible for removing calcium ions from the cytosol and returning them to the SR. Once the cell has sufficiently repolarized, the activation cycle can start again.

The regulation of the RyR receptor has been subject of intense research, reviewed in, e.g., [31, 32]. Malfunctioning of this receptor has been implicated in heart failure, the condition during which the heart muscle is progressively weakened and unable to pump blood through the body [33]. It has also been connected to the initiation of cardiac arrhythmias, including ventricular and atrial fibrillation. In this case, it is believed that the leakage of the RyR receptor can result in spontaneous depolarization events and aberrant electrical activation, as further discussed in Section 4 [34].

3.1. Single cells

To cast the time-dependence of the membrane potential of a cardiac cell, V, in the form of mathematical equations, one can view the cell membrane as a capacitor, with capacitance C_m per unit area (μF cm⁻²), in parallel with a resistor, representing ionic currents. Since there cannot be any net accumulation of charge on either side of the membrane, the sum of the current through the channels and the capacitive current needs to be zero:

Here, the last term describes all ionic currents (μA cm⁻²) considered in the model and where, by convention, inward currents are negative and outward currents are positive. The dynamics of these currents can be described in terms of ordinary differential equations, as first shown by the seminal work of Hodgkin and Huxley [49]. They studied action potentials in the squid giant axon using voltage clamp experiments, which measured the dynamics of the ion channel conductances. This study revealed that the generic form of a current, say X, is:

where E_X is its reversal potential, g_X represents the conductance per area, and f is a function that depends on gating variables x_i, i = 1, …, N. For example, the potassium current per unit area for the giant squid axon can be accurately described as I_K = g_Kn⁴(V − E_K), where n is a dimensionless gating variable, given by:

In this expression, α_n and β_n are functions of V, chosen to match experimental results.

Ionic channels in myocytes (some of which are shown in Fig. 2A) can be described using Hodgkin-Huxley-like expressions. During the past decades, a large number of cardiac models that formulate these ion channels have been developed with ever increasing complexity and details (see, e.g., the review by Fenton and Cherry [50]). A major contributor to this increase of complexity is the introduction of calcium dynamics, absent in earlier and simplified models. Models have been developed for both ventricular [51, 52, 53] and atrial cells [54, 55, 56] and have been used to give insights into cardiac electrical dynamics, including arrhythmias (see Section 4).

In certain cases, for example for specific mutations of certain channels, the Hodgkin-Huxley formalism may not be sufficiently detailed to describe the electrical properties of the cell. For these cases, it is possible to formulate a Markov-based description of the gating of the channel. In this description, the ion channel is represented by a number of discrete states, corresponding to different configurations of the channel and the model incorporates the transition rates between these states, which can be found using fitting techniques [57]. These Markov formulations have been used to determine how certain genetic defects may be linked to cardiac arrhythmias [58, 59]. Furthermore, this type of model can be used to give insights in the effects of drugs that target channels; we will see an example of this in Section 6 [60]. Several studies have discussed the merits of the Markov vs Hodgkin-Huxley formalism [61, 62].

There are a number of drawbacks to using very detailed models that attempt to incorporate all possible ion channels. First of all, the number of variables is very large, which can make simulations more time consuming and less intuitive. For example, a Markov-based model for the human epicardial ventricular myocyte contains 67 variables [52]. More fundamentally, however, it is challenging to accurately describe all ion channels operating in human cardiac cells. The equations that describe the identified channels are highly nonlinear with many poorly characterized parameters. Some of these channels are not measured in human cells and are, instead, fitted to data from animal cells. Furthermore, a human heart consists of many different cells with different properties and capturing this heterogeneity is challenging at best. Finally, in addition to this patient-specific heterogeneity, there exists considerable patient-to-patient variability. This is the case for healthy individuals but especially for individuals that are suffering from cardiac diseases. Thus, it seems an impossible task to construct a universal electrophysiological cell model.

An alternative approach to developing highly detailed models is to forgo the complexity inherent to the many ion channels and create a model that can replicate some of the essential features of a cardiac cell with only a few variables. Such models should, of course, be able to create realistic action potential, characterized by a fast upstroke (i.e., increase) of the voltage, followed by a gradual return to equilibrium. In other words, these models should have excitable or excitable-like dynamics such that a small stimulus will have a significantly different and smaller response that a larger stimulus. Furthermore, these models should be able to display some of the spatially extended dynamics thought to underlie cardiac arrhythmias (see Section 4). In particular, it should exhibit spiral wave dynamics, which we will discuss in more detail in the next sections. The advantage of these models is that they are much easier to implement, are faster to simulate, and can sometimes be more intuitive. Of course, the price one pays is that determining the role of specific channels is no longer possible. Nevertheless, it was through the use of such models that some of the critical insights into spatio-temporal organization during cardiac arrhythmias were gained [63, 64].

One of the simplest models one can write down for a cardiac cell is the FitzHugh-Nagumo (FHN) model [65, 66]. This model contains only two variables, one of which (v) can be thought of as the membrane potential while the other one (w) represents all currents of the cardiac cell and functions as the recovery variable:

In these equations, the parameter ϵ is taken to be much less than 1, such that v is the fast, or excitation, variable and w is the slow recovery variable. Furthermore, f is typically taken as cubic function while g is a linear function. Although many different forms can be taken, a particular and convenient form corresponds to f(v,w)=v−v33−w and g(v, w) = v + a − bw where a and b are parameters of the model [67].

For such a model, it is possible and instructive to draw a phase diagram in the (v, w) space along with the nullclines, defined as the curves on which v.=w.=0. The intersection of these nullclines then gives the fixed point [68]. Depending on the parameters of the model, this fixed point is either stable or unstable. Fig. 4A shows the nullclines for the convenient choice of f and g and for parameters for which the fixed point, indicated by the yellow dot, is stable [69]. Transiently perturbing the variables away from this fixed point, for example by introducing a non-zero excitation current to either the v or w variable, results in dynamics that depends critically on the amplitude of the perturbation. A small perturbation will result in a small deviation, followed by a return to the fixed point. An example of the dynamics in the (v, w) space following such a perturbation is shown as the red trajectory labeled by 1 in Fig. 4A while the corresponding time trace is plotted in Fig. 4B. This is in sharp contrast to the dynamics following a larger perturbation, represented by the trajectory and timetrace labeled 2 in Fig. 4A&B. In this case, the perturbation results in a large excursion in phase space. Due to the small value of ϵ, the v variable rapidly switches from its perturbed value to the right-most branch of the v nullcline. Once this nullcline is reached, the variables slowly follow it until they reach its maximum value. This is then followed by a rapid transition to the left-most branch of the v nullcline and a slow relaxation to the fixed point. For small values of ϵ, the time of the phase excursion can be estimated by computing the time spent on the v nullcline. For the specific form of the FHN model of Fig. 4, the time spent on the right-most branch of the v nullcline is found by evaluating

where the integral runs from the solution of f (v*, −0.625) = 0 on the right-most branch, to v = 1, corresponding to the maximum of the nullcline. A similar expression can be derived for the time spent on the left-most branch of the nullcline.

Open in a separate window

Figure 4:

The FHN model. A: Nullclines (black) of the FHN model, along with the the streamlines of the phase portrait and the equilibrium point (yellow dot). A small stimulus initiates the red decay trajectory, marked by 1, whereas a large stimulus initiates the red excitation trajectory, marked by 2. Parameter values of the model (see text): ϵ = 0.08, a = 0.7, and b = 0.8, resulting in a fixed point located at v=−1.2 and w=−0.625. B: Time trace of the variable v corresponding to the trajectory #1 and #2 in A.

Even though the FHN model does not posses an actual firing threshold and does not display an all-or-none response [65], its dynamics is qualitatively similar to action potential generation, which is characterized by a rapid upstroke, followed by a slow relaxation to the equilibrium value. Furthermore, it also exhibits the refractoriness seen in cardiac cells: following an excitation as exemplified by trajectory #2 in Fig. 4, a new, large response can only be elicited when the variables have sufficiently recovered. This allows the FHN model to recapitulate generic excitable phenomenon when specific ionic dynamics are not important. However, it lacks several critical elements of the dynamics of cardiac cells and it is certainly not an electrophysiologically accurate model.

A more recent model, specifically developed to mimic cardiac dynamics, was introduced by Fenton and Karma (FK) [70]. In its original version, it contained only three currents (a fast and slow inward, and a slow outward current, which can be viewed as representatives of the Na⁺, Ca²⁺, and K⁺ channels, respectively), along with two gating variables and, of course, the voltage. Importantly, and distinct from the FHN model, the parameters of the original FK model and its extensions can be adjusted to mimic features of more detailed models or clinical data [70, 71, 72, 73]. An example of this is shown in Fig. 5A, where the AP shape of a AF patient undergoing an AF ablation procedure is plotted as a dashed line. The AP shape for this patient, and 4 other patients, was obtained during controlled pacing in a patient using a monophasic action potential (MAP) catheter [73]. This type of catheter is able to faithfully record the membrane potential at a specific location [74, 75]. In this study, the FK parameters were adjusted to fit the clinical AP shapes and the results were excellent (see, e.g., the blue line in Fig. 5A). The same study also showed that the parameters of a detailed model, created by Koivumaki, Korhonen, and Tavi (KKT) and that contains 13 currents [56], can be adjusted to fit the AP shape. The result of fitting the parameters of this “KKT” model is shown in Fig. 5A as a red line. As is the case for the FK model, the KKT model parameters can also be fitted to produce good agreements with clinically obtained AP shapes.

Open in a separate window

Figure 5:

Model fitting. A: Clinical AP morphology from an AF patient compared to model AP morphologies obtained by fitting the model parameters for an AF patient. The average clinical AP morphology is shown as a dashed line while the FK model is shown in blue and the KKT model is in red (From [73]). B: Schematic representation of a cardiac cell activated by stimuli with period T. As a result, the membrane potential increases rapidly, followed by a gradual return to the resting potential. The APD (action potential duration) measures the time at which the potential is above a certain threshold while the DI (diastolic interval) measures the time interval between the end of one and the start of another action potential. C: Clinically determined APD and its polynomial fit for an AF patient (open and closed symbols, respectively) and the results from a fitting procedure for the FK (blue) and KKT (red) model. D: CV restitution curves of several ionic models and of measurements in a guinea pig (solid lines) together with those obtained from the fitted FK model (From [70]). (BR: Beeler Reuter [76], MBR: modified Beeler Reuter with speed-up calcium, MLR-I: Luo-Rudy with speed-up calcium [77], GP: Guinea pig [78]). E: As in C, but now for the clinically determined CV restitution curve.

Aside from the AP shape, the parameters of the FK model can also be adjusted to replicate several key dynamical characteristics of electrophysiological models and of patients. One example is the APD restitution curve, which describes how the APD recovers as a function of the rest period between repolarization and the next excitation [79]. As we will see in the next section, the APD restitution curve determines the stability of a cardiac cell. The curve can be computed by stimulating a cell using a current with period T, as schematically shown in Fig. 5B. Using the time trace of the resulting potential, we can then define the diastolic interval (DI), which is the time spent by the cell below a certain threshold potential value before a new stimulus (Fig. 5B). By varying the period T we can construct the APD restitution curve, which plots the APD as a function of the DI. In general, this curve has a positive slope: as the DI is reduced, the cell has less time to recover from the previous excitation. Consequently, channels do not open or close fully, resulting in a shorter APD. As T is further reduced, the APD continues to decrease until an excitation is no longer possible. We should point out that different qualitative shapes of the APD curve are possible, including multi-branch structures or ones where the slope is negative [80, 81, 82].

An example illustrating that the FK model (and the KKT model) can be fitted to clinically obtained data is shown in Fig. 5C. Here, the APD vs DI (open symbols), together with its high-order polynomial fit (grey closed symbols) are plotted for the AF patient presented in panel A. The fit of the FK model to the polynomial restitution curve for this patient is plotted in Fig. 5C as a blue line, while the fit for the KKT model is plotted as a red line. Both models were found to be able to replicate the restitution curves of all 5 patients [73].

These results demonstrate that simplified models are able to replicate certain cardiac properties, which make them useful tools to investigate cardiac dynamics in general and arrhythmias in particular. Clearly, as stressed before, the FK model, and similar simplified models, cannot describe the relative importance of specific channels. They can, however, simulate voltage dynamics in a more intuitive and efficient way. Both detailed and simplified models have their role in research and perhaps a hybrid approach, in which a detailed model is systematically reduced to a much simpler one, represents an attractive middle ground [83].

3.2. Cable propagation

So far, we have examined the dynamics of a single cell. Let us now consider what happens if we connect cells and form a cable. This 1-dimensional propagation can be captured by the cable equation, investigated in many studies (see, e.g., [84, 85, 22]). To derive this model, we assume that cells with radius a are coupled with a resistance per unit length equal to R_a and have a capacitance C_m. Then, applying Ohm’s law and Kirchhoff’s law, we can derive, in the continuous limit, the cable equation

where D = 1/(4πaR_aC_m) [85]. Thus, the spreading of the potential V through the cable is governed by a reaction-diffusion equation with a diffusion constant D.

Stimulating a cable on one end with a sufficiently strong stimulus will result in a cardiac wave that propagates along the cable to the other end. It can be shown that the conduction velocity (CV) of this propagating wave scales as v∼D [85]. In the case when the cable is paced at regular intervals, similar to the the periodic pacing of a single cell in Fig. 5B, this CV also depends on the period and thus the DI. The property is captured by the CV restitution curve, which plots the CV as a function of the DI. Just as the APD restitution curve, this curve typically has a positive slope: the CV decreases when the pacing frequency is increased. This is shown in Fig. 5D, where the CV restitution curve, obtained by stimulating one end of a 5 cm long cable, is plotted for several electrophysiological models [70]. Mechanistically, the CV restitution arises from the large timescales required to reset the quantities of the cell, most notably calcium. This plot also shows that the FK model parameters can be adjusted to replicate these CV restitution curves. The ability to adjust model parameters to fit CV restitution curves is further demonstrated in Fig. 5E. This plot shows the fits from both the FK (blue line) and KKT (red line) to the clinical CV restitution curve of the AF patient presented in panels A and C.

4.1. Single cell instabilities

4.1.1. Alternans

When a single cell is paced with a constant frequency, the cell can develop an instability, resulting in an APD that is not constant. One particularly well studied example of such an instability is electrical alternans [127, 128, 129, 130, 131, 132, 133, 134]. During alternans, the sequence of APDs of a cardiac cell alternates between long and short and the first reports date back more than a century ago [135]. Alternans became a focus within the cardiac community when the beat-to-beat alternation of the T-wave on the ECG was recognized as a precursor to ventricular arrhythmias and sudden cardiac death [136, 137]. It garnered further attention when it was found that this instability may also play a role in the initiation of atrial fibrillation [138].

The slope of the APD restitution curve (see Fig. 5) has long been recognized as an important determinant for the initiation of electrical alternans [139, 140, 141]. As we recall from Section 2, this curve relates the APD to the preceding DI and in Fig. 7A&B we have plotted a restitution curve with a slope of less than and larger than 1, respectively. For a periodically paced cell, the sum of the APD and DI should be equal be equal to the pacing period T (APD + DI = T) and this line is plotted as a dashed line in Fig. 7A&B. To investigate whether the intersection of this line and the restitution curve (red point in Fig. 7A&B) is stable we can construct a cobweb diagram, a standard graphing technique in dynamical systems [142, 143]. To construct this diagram, we perturb the DI and find the corresponding APD from the restitution curve (blue point in Fig. 7A&B). This APD results in a DI that can be found using the period relationship (dashed line) which, in turn, leads to the next APD. Repeating this procedure, we can see that for restitution curves with a slope < 1 the perturbation decreases, indicating that the point is stable (Fig. 7A). For a slope > 1, however, the perturbation grows and the APD and DI alternate resulting in long-short-long sequences (Fig. 7B). An example of this sequence is shown in Fig. 7C where the FK model is stimulated with a period corresponding to a slope of the restitution curve that is larger than one. As a result, the APD is long for one stimulus and short for the following.

Open in a separate window

Figure 7:

Alternans in cardiac myocytes. A: Cobweb diagram illustrating a stable fixed point, shown in red, on a APD restitution curve with slope < 1. The dynamics of a small perturbation in DI, corresponding to the blue point, can be followed using the restitution curve APD_n+1 = f(DI_n) and the relationship APD_n + DI_n = T (dashed line) and is shown to approach the fixed point, indicating stability. B: As in (C) but now for a restitution curve with a slope larger than 1. In this case, a point on the curve is unstable, as can be seen by following the path of the perturbation. C: Example of alternans in a single cell. Shown here is the response of the FK model (parameter set #4 from Ref. [87]) to periodic pacing with a period of 220 ms. The resulting APD oscillates between long and short.

The instability condition can also be derived algebraically. To determine the stability of a point APD* and DI* on the restitution curve, consider perturbing the diastolic interval by δDI: DI₁ = DI* + δDI. For a linear restitution curve with slope α, the next action potential will have a duration of APD₂ = APD*+αδDI and the subsequent DI is given by DI₂ = DI* − αδDI. Repeating this, we find that the n-th APD is given by

APDn=APD∗+(−1)nαn−1δDI

(15)

As one can see from this result, APD_n will approach APD* as long as α < 1 but will diverge from APD* if α > 1. Thus, points on the restitution curve with slope smaller than one are stable while points with a slope larger than one are unstable.

Aside from electrical alternans, manifested by APs with alternating durations, experiments have also shown alternans in intracellular Ca²⁺ cycling [144]. The fact that Ca²⁺ cycling shows alternans is perhaps not surprising, since there is a directional coupling between voltage and Ca²⁺ handling. As we have discussed in Section 2, Ca²⁺ is released from the SR upon depolarization and pumped back by SERCA pumps during the relaxation phase of the action potential. Thus, voltage regulates the release of Ca²⁺, which means that electrical alternans can automatically result in alternating Ca²⁺ signals. In fact, this coupling is bi-directional since the amount of released Ca²⁺ also affects the LCC (it inactivates in a Ca²⁺ dependent manner) and the sodium-calcium exchanger (positively regulated by intracellular Ca²⁺).

Initially, it was not clear if calcium alternans could occur without electrical alternans. Experiments in which the voltage was clamped and therefore did not alternate have shown, however, that calcium alternans can occur even in the absence of APD alternans [145]. Furthermore, additional experiments have demonstrated that APD and calcium alternans can even be out of phase, with the voltage displaying a short-long-short sequence while the Ca²⁺ concentration shows a long-short-long sequence [146]. An example of calcium alternans in a ventricular myocyte of a rat is shown in Fig. 8A. Here, the Ca²⁺ concentration along a line within the myocyte is plotted as a function of time for three successive beats using a fluorescent Ca²⁺ indicator [147]. Clearly, the Ca²⁺ release is large following the first and third beat, but not following the second beat.

Open in a separate window

Figure 8:

Calcium alternans in experiments and simulations. A: Linescans in a rat ventricular myocyte in response to 3 consecutive voltage pulses, visualizing Ca²⁺ using a fluorescent indicator. From [147]. B: Normalized Ca²⁺ release in a 2D slice of a simulated cell, with green/red indicating minimal/maximal Ca²⁺ release. The pattern shows alternans, with a small release in beat 445 and a large release in beat 446. From [154]. C: A line scan plot of cytoplasmic Ca²⁺ from a simulation of coupled CRU networks, illustrating miniwaves during the release beats. From [148].

The mechanisms responsible for calcium alternans have been actively investigated since its discovery (for a review, see e.g. [148]). To fully understand calcium alternans requires a more detailed look how Ca²⁺ is released from the SR following excitation. It turns out that this release does not necessarily occur uniformly within the entire myocyte. This is because Ca²⁺ is released from so-called calcium release units (CRUs), consisting of a few LCCs and roughly 100 RyRs. A typical myocyte contains approximately 20,000 CRUs that are closely packed in a three-dimensional 3D grid [149]. Calcium release from a single CRU, called a spark, is highly localized both in time (a typical spark lasts roughly 30 ms) and in space ~ 2 μm. A spark is thus the elementary event in cardiac excitation-contraction coupling [150]. Furthermore, calcium sparks are stochastic (since LCCs and RyRs have random opening properties) and refractory: the probability of releasing calcium immediately following a spark is reduced and slowly recovers over ≈400ms [149]. In addition, the CRUs are coupled through the diffusion of Ca²⁺, which can lead to waves of Ca²⁺ within a myocyte [151].

Computational models of coupled CRUs have demonstrated calcium alternans, as shown in Fig. 8B&C. Fig. 8B shows the calcium release of a 2D slice of a simulated cell and, as in the experimental result, the Ca²⁺ release alternates between large and small values. Furthermore, as shown from the line scan in Fig. 8B taken from a simulation of coupled CRUs, a wave of Ca²⁺ is triggered and propagate through the myocyte, consistent with experimental results. To further gain insight into calcium alternans, it is possible to derive a recurrence equation that relates the number of sparks in one beat to the number of sparks in the preceding beat [148]. For this, it is assumed that a CRU can fire (i.e., release Ca²⁺) either through the direct activation by LCCs, with probability α, or in response to the firing of a neighboring CRU, with probability γ. The latter can arise from an increase in local Ca²⁺ due to diffusion. If we assume that CRUs recover with a probability β, that there are N₀ total CRUs present, and that the fraction of recruited CRUs through diffusion is f, we obtain for the number of sparks in beat i + 1, N_i+1, the following expression:

Ni+1=α(N0−βNi)+(1−α)(N0−βNi)f

(16)

In this expression, the first term represents firing from CRUs induced by the opening of LCCs while the second term represents firing due to diffusion from neighboring CRUs [152, 153]. For a well-mixed system, where the firing probability is uniform in space, the fraction f can be shown to be given by

f=1−(1−αγ[1−βNk∕N0])M

(17)

where M is the number of nearest neighbors of a CRU. The iterative map can be analyzed to determine its steady state and its stability. This analysis reveals that a period-doubling instability, which represents alternans, can only occur for non-zero values of γ and β and when α deviates from one. In other words, alternans results from recruitment, refractoriness, and the stochastic firing of the CRUs.

4.2. Wave block

Alternans as described above is a single cell instability. However, most cardiac arrhythmias, including AF and VF, are spatio-temporal events and need to be investigated in spatially extended domains. One critical ingredient in the initiation of spatially extended arrhythmias is wave block, which occurs when the activation front is no longer able to propagate further and stalls. This can be either because the region is non-conductive due to, for example, injury or because it is temporarily refractory due to its specific electrophysiological properties. In either case, the wave stalls and, in a 1D cable, activation would terminate. In 2D or 3D, however, a different scenario can develop in which the wave is blocked in a particular region but continues to propagate around it. Under certain conditions, this can result in a wave that curls around this region but “reenters” it after it has sufficiently recovered. This reentry is believed to be of paramount importance in the initiation of arrhythmias.

Wave block and reentry can be nicely illustrated by a conceptual “pinwheel” experiment, proposed by Winfree [155, 156] and schematically shown in Fig. 9A. In this figure, an initial stimulus has produced a planar activation wave, indicated by the green line and propagating from left to right. Behind the activation front the tissue is slowly recovering, illustrated by the graded gray scale, and is in the refractory state. After some time, a second stimulus is applied and the resulting activation pattern will depend critically on the timing of this second stimulus. Given immediately following the activation front, the entire stimulus will be firmly in refractive tissue and will not result in any additional activation. If the extra stimulus is provided a long time after the initial wave front has passed, the tissue will no longer be refractory and a new circular wave front will propagate radially outwards. For stimuli between these two extremes, something more interesting will happen, as indicated in the figure: the stimulus will result in propagation in the retrograde direction while it will be blocked in the opposite direction, as indicated by the green and red symbols, respectively. As the retrograde front progresses, the refractory tissue will further recover, eventually becoming excitable. As a result, the wave front can now curl around and reenters the previously inexcitable tissue, as depicted by the two curved green arrows in Fig. 9A.

Open in a separate window

Figure 9:

Initiation of spiral waves using a premature stimulus. A: Schematic drawing of a stimulus (white circle) that follows a planar wave that propagated from left to right. Proper timing of this stimulus results in wave block to the right, indicated by the red symbol, but not to the left, as indicated by the green arrows. As the newly created wave propagates in the retrograde direction, tissue in the right-half of the domain continues to recover, eventually leading to reentry and spiral wave formation. B-E: Snapshots of a simulation showing the formation of a figure-of-eight reentry pattern as the result of a properly timed extra stimulus. Panel B shows the premature stimulus along with the planar wave propagating from left to right. Panel C demonstrates that only propagation in the retrograde direction is possible while panel D illustrates the ensuing reentry into the recovered tissue. Panel E shows the figure-of-eight reentry pattern. Snapshots were created using the FK model and numbers correspond to time in ms.

The result of this properly timed extra-stimulus are thus two reentry waves. Given enough space, these reentry waves can continue to curl, leading to two counter-rotating spiral waves and a so-called figure-of-eight reentry pattern. We should point out that these spiral wave solutions are not only found in simulations of cardiac tissue but are generic solutions of spatially extended excitable systems. They are found in a variety of other systems, included in aggregation patterns of cells of the social amoeba Dictyostelium discoideum cells [157], in chicken retinas [158], in the oxidation of surface catalytic systems [159], and in chemical Belousov-Zhabotinsky systems [160, 161].

An example of the creation of spiral waves using the above pinwheel scenario is shown in Fig. 9B-​-E.E. Here, we used the FK model [70] and first created a planar wave that propagated from left to right in a 5x5 cm² computational domain. We then excited a region behind the activation front by adding a constant current to the equation of the voltage for a duration of 1ms (panel B). As can be seen from the snapshots, the part of the stimulus directed towards the planar front failed to propagate since it falls within the refractory period. The back part of the stimulus, however, was able to propagate into freshly depolarized tissue (panel C) and reentered into the previously refractory tissue (panel D). In this simulation, the amount of space was sufficient for the two reentry waves to curl backwards (panel D), resulting in two spiral wave and a figure-of-eight pattern (panel E).

As mentioned above, the timing of the additional stimulus is critical in the formation of a reentry pattern: a stimulus that is given too early will not result in any excitation while a stimulus that is given too late will result in centrifugal propagation. This means that there is a specific time interval ΔT, termed the vulnerable window, which will lead to reentry [162, 163]. The size of this vulnerable window was estimated by Winfree as the ratio of size of the circular stimulus and the conduction velocity of the planar front, v_plane [156]. This can be intuitively understood by realizing that a very fast moving wave should have a small vulnerable window. In addition, a stimulus with a large spatial extent will more easily cover both the refractive and non-refractive part of the tissue and will thus have a larger vulnerable window. A more precise result, using simple geometric arguments, was obtained by Yang et al.: ΔT=(R−Rc)∕vplane [164]. In this expression, R_c is the critical stimulus size, below which propagation is not possible.

The generation of spiral wave reentry using a properly timed premature stimulus has been demonstrated in the laboratory as well [165, 166, 167]. The patterns observed in the experiments, however, are typically more complex than shown in Fig. 9. One reason for this is that cardiac tissue is anisotropic while the simulation in Fig. 9 was carried in an isotropic domain. Furthermore, simulations using the bidomain model, introduced in Section 3 and which takes into account both intracellular and extracellular resistivities, showed that more intricate patterns are possible [168, 169]. Some of these patterns were subsequently verified in an experimental study [170].

4.2.1. Alternans and reentry

Let us come back to alternans, characterized as an APD with alternating long and short duration. In Section 4.1.1 we introduced alternans as a single cell instability and it is not immediately clear how an instability at the single cell level can result in dangerous arrhythmias, as reported in the clinical literature [136, 137]. After all, coupling identical cells in a 1D cable or a 2D sheet will result in so-called spatially concordant alternans. In this dynamical state, not considered to be pro-arrhythmic, the entire spatial domain shows alternans with all elements in phase. In other words, during concordant alternans the alternating APD sequence is uniform in space with each point exhibiting the same short-long-short sequence. The pro-arrhythmic nature of single-cell alternans is related to a spatio-temporal instability that can develop when we couple these cells. This additional instability was observed in experiments in rings of cardiac tissue, which showed that the propagation of a pulse can become unstable if the circulation time around the ring was reduced [171]. This instability was also observed in computer simulations of a ring of cardiac tissue and in optical mapping experiments of the activation activity in guinea pig hearts [172, 129]. These studies showed that spatially extended systems can exhibit spatially discordant alternans, characterized by spatially distinct regions in which the alternans is out of phase. In other words, the alternans is no longer spatially uniform and regions that show small-large-small APD are adjacent to regions that show large-small-large APD. Furthermore, these studies also showed that discordant alternans creates large APD gradients, which can result in unidirectional propagation block. This wave block can then result in reentry, which can initiate arrhythmias, including life-threatening ventricular fibrillation. For a review on spatially discordant alternans, see e.g. [173, 174].

One obvious way of creating spatially discordant alternans is to have heterogeneous tissue with different properties in different regions. As shown by Watanabe et al., however, homogeneous tissue that display concordant alternans can also create spatially discordant alternans [175]. Specifically, they showed that such discordant alternans can be initiated in one of two distinct scenarios. In the first, a cable is paced on one end and a extra stimulus is given at the other end. For the right pacing rate and timing of the premature stimulus, the interaction between the activation wave from the pacing site and from the extra stimulus can create two regions in which the alternans exhibits different phases. This scenario, however, requires fine tuning of the pacing rate and timing of the premature stimulus and may not be of clinical relevance.

The second scenario involves CV restitution, already encountered in Section 3 (see Fig. 5). This scenario is illustrated in a computer simulation shown in Fig. 10A. The cable is paced at one end (top) with fixed rate and the first stimulus, shown as the left-most sequence and labeled A0, is applied after the tissue was at rest for a prolong period of time. As a result, the APD is long but the following DI is short. Due to CV restitution, the second stimulus (A1) has a smaller CV than A0, resulting in an increase in APD, and thus DI, as the pulse travels down the cable. The third stimulus follows a long DI and has thus a large CV. This, together with the increase in APD along the cable from the previous beat, results in decreasing DI and APD as the pulse travels down the cable. The next pulse, DI and APD increases again from top to bottom of the cable. This sequence of different speeds in successive beats is repeated, resulting in a spatial gradient of the APD. For steep enough CV restitution curves, eventually a node is formed, which separates a spatial region that exhibits small-large-small APD with a region that shows large-small-large APD. This is shown in the remainder of this panel, which plots the 13th and higher stimulus: both the top and bottom show alternating APDs that are out-of-phase. In 2D, this scenario results in out-of-phase regions that are separated by a nodal line, as illustrated in Fig. 10B. This nodal line separates the two out-of-phase regions and points on this line do not exhibit alternans.

Open in a separate window

Figure 10:

Alternans in spatially extended domains. A: Initiation of discordant alternans in an 8 cm long computational cable using the Beeler-Reuter model (vertical axis denotes space and horizontal axis denotes time). All stimuli were applied at the top end of the cable, using constant 310 ms intervals. After multiple stimuli, the top and bottom of the cable are alternating out of phase. From [175]. B: Spatially discordant APD alternans in simulated 2D cardiac tissue using the Luo-Rudy model and a pacing cycle of 220 ms. The top panel displays the action potential at sites a and b, showing a long-short pattern. The lower panel shows that the spatial APD distribution, color-coded with blue for short and red for long, for two consecutive beats. The white line represents the nodal line, which separates the out of phase domains and where alternans is absent. C: Sequence of events following a premature ectopic beat (asterisk) that can result in reentry. The beat is initially blocked when it attempts to propagate into the region with long APD. It can reenter this region, however, after it has propagated along the nodal line, resulting in a figure-eight re-entry. B&C from [173].

The spatiotemporal dynamics of alternans of 1D and 2D paced tissue was mathematically analyzed by Echebarria and Karma [176, 177]. They derived an equation for small amplitude alternans and showed that a linear instability could lead to nodes in paced tissue that are either stationary or traveling. Furthermore, they showed how the wavelength of the pattern is a function of the electrotonic coupling and the CV restitution curve. Finally, the fact that discordant alternans can spontaneously arise in homogeneous tissue was also shown experimentally using optical mapping in paced rabbit hearts [178].

How is alternans in spatially extended system pro-arrhythmic? The answer lies in conduction block and reentry. Envision an extra stimulus in the spatial region of short APD, as indicated by the asterisk in the blue part of Fig. 10C. Within this spatial region, the resulting activation wave is able to propagate. Once it crosses the nodal line, however, it encounters the region of long APD. For an appropriately timed extra stimulus, this encounter will result in propagation that is blocked along the black line in the figure. However, as the wave is blocked, it continues to propagate in the direction parallel to the line of block. Eventually, the tissue on the long APD region repolarizes and activation can reenter. This will then potentially initiate spiral wave formation.

4.3. Tissue heterogeneities and wave breakup

Key ingredient in the above scenario is the large gradient of APD at the nodal line. This gradient, created even though the tissue is homogeneous, is purely due to the dynamics of the tissue and results in wave block and failure of propagation. It is also possible, however, to create wave block and subsequent reentry in the presence of tissue heterogeneities. These heterogeneities can be pathological, due to disease, but also occur naturally in the heart. For example, the ventricular shows a significant transmural gradient in the APD, even in non-failing hearts [179]. Such heterogeneities can result in the local failure of propagation and in reentry.

The effect of heterogeneities on wave propagation has been studied extensively in both 2D and 3D computer models [180, 181, 182, 183, 184, 185]. An example of a heterogeneity-induced spiral wave is shown in the simulation snapshots of Fig. 11A [180]. Here, the FHN model is simulated in a rectangular domain with non-periodic boundaries. This domain contains a heterogeneity bounded by the red lines within which the refractory period is larger. Planar waves enter the domain from the left hand side. As a result of the increase in refractoriness in the heterogeneous region, the second planar wave is blocked when it encounters the heterogeneity. However, since the second wave has a lower velocity than the first one, eventually it is able to enter the heterogeneous region, resulting in reentry and a spiral wave.

Open in a separate window

Figure 11:

Examples of wave breakup in models with heterogeneities. A: Snapshots showing successive times in a simulation of the FHN in a domain that contains which a heterogeneous zone of refractoriness, indicated by the red lines. Due the heterogeneity, the wave is originally blocked but then reenters the heterogeneous zone, after which a spiral wave is formed. After [180]. B: Reentry by interaction of an activation wave with a circular and fibrotic region in a computational model using the modified Beeler-Reuter model. Total area 77 cm². From [186]. C: Initiation of reentry in a patient-specific model shown in 4 snapshots, 60ms apart. Reentry is shown using white arrows and is due to fibrosis in that region. From [187].

Heterogeneities can also form in diseased or aging hearts, either acutely or over time. An example of the former is present in acute ischemia, during which part of the heart is deprived from oxygen. This lack of oxygen affects several ion channels and can result in significantly reduced conduction velocities and wave block [188]. An example of the latter can be found in the fibrotic remodeling of cardiac tissue. Fibroblasts are the most numerous cell type in the heart and play a critical role in maintaining the mechanical structure of the heart [189]. Due to aging or injury, these fibroblasts can be activated and differentiate into myofibroblasts, resulting in the deposition of excessive amounts of collagen fibers and fibrosis [190, 191, 192, 193, 194]. These collagen fibers can then act as conduction barriers between adjacent myocyte strands, altering the conduction properties of the tissue. Furthermore, fibrotic remodeling also affects the gap junction coupling between myocytes, resulting in a decrease in the longitudinal and transverse conduction [195, 104].

To incorporate the effects of fibrotic remodeling in computational models is challenging and is an active field of research, as reviewed in [196]. These models need to provide a description of the coupling between myofibroblasts and myocytes, the alteration of gap junction coupling, and the inclusion of non-conductive collagen fibers. The latter is particularly challenging since fibrosis can manifest itself not only in compact patterns, characterized by a relatively large area of non-conductive tissue, but also in diffuse patches, with short collagen strands scattered among myocardial fibers, or patchy patterns, characterized by small regions of fibrotic tissue [197]. These microstructure patterns require resolving propagation on small length scales, and it is still unclear if continuous modeling can appropriately mimic fibrosis [198]. For this reason, several studies are employing discrete approaches in which tissue is represented by discrete elements [199, 200, 201]. Furthermore, hybrid models, which combine continuous and discrete modeling approaches, are also being developed [202].

Regardless of the computational approach, modeling has shown that the inclusion of fibrosis can result in wave breakup and reentry [203, 204, 205, 206, 207, 186, 208]. Especially diffuse or patchy fibrosis is pro-arrhythmic. An example from one of these simulations is shown in Fig. 11B, where a circular patch of fibrosis, modeled as a region with a high fraction of non-conducting links, is embedded in a 2D computational sheet [186]. A propagating planar front breaks inside the fibrotic region and reenters the homogeneous region, resulting in reentry and continuous activation of the tissue.

Clinical studies using advanced magnetic resonance imaging techniques have shown that fibrotic regions are linked to ventricular arrhythmia [209]. In addition, imaging studies have demonstrated that patients with more severe fibrosis have a larger probability of recurrent arrhythmia after AF ablation [210]. Recently, the results of these imaging studies were incorporated into patient-specific computer models of atrial wave propagation [206, 104, 207, 187, 106, 211]. An example of such a study is shown in Fig. 11C, where four snapshots of propagation patterns after the initiation of AF are plotted [187]. In this simulation, which uses a patient-specific geometry and fibrosis distribution, premature stimulation was able to generate reentry at a spatially conserved location. Importantly, the locations of these reentry events were highly correlated with zones of increased fibrosis, demonstrating that fibrotic remodeling can play a significant role in the generation of reentry and arrhythmia.

4.4. Spiral wave breakup

A single or a pair of spiral waves, as in the figure-of-eight pattern shown in Fig. 9, is likely not fatal to patients even if they form in the ventricles. It may result in an elevated heart rate (tachycardia) since the spiral wave will have a period that is shorter than the typical sino-atrial node period [212]. The activation of the heart, however, will still be mostly coherent, resulting in reduced but non-fatal pumping ability. The activation patterns during fibrillation, on the other hand, are less coherent than those created by a single of a pair of spiral waves and, most likely, involve multiple and possibly migrating spiral waves [213, 214]. This would mean that the initiation of fibrillation would encompass not only the creation of a spiral wave but also the break up of this spiral wave [127, 215].

Moe and co-workers, using a cellular automaton model, were the first to hypothesize that the breakup of spiral waves could underlie cardiac fibrillation [63]. Their computational model showed that a single spiral wave could break up in multiple wavelets in a stochastic fashion. These wavelets migrated through the domain and continued to be annihilated and born, resulting in a self-sustained computational arrhythmia.

Since this seminal work, many computational studies have demonstrated spiral wave breakup [216]. These simulations include algorithms that treat the medium as a continuum and cast the equations for the voltage in terms of the partial differential equations described in Section 3 (Eq. 8 and Eq. 10). Furthermore, the mechanisms behind spiral wave breakup have also been thoroughly investigated. Karma, for example, demonstrated that alternans can destabilize spiral waves [217, 128]. Specifically, he showed that alternans results in oscillations of the APD of the spiral wavefront. Once these oscillations become large enough, they lead to conduction block and to spiral wave break at that location.

Alternans is not the only mechanism that can create spiral wave breakup in homogeneous tissue. In an exhaustive numerical study, Fenton et al. revealed a number of potential mechanisms for spiral wave breakup [87]. For example, aside from alternans, which will result in APD oscillations and wave break far from the tip, they showed that a steep APD restitution can also result in breakup close to the tip. In this case, the CVs of the front and back velocity of the wave are no longer identical, resulting in wave block and breakup. Another mechanism involves a Doppler shift in the wave frequency due to tip meandering. The Doppler effect arises since the tip, which can be viewed as a wave source, is moving so that the frequency is higher in the direction in which the spiral is moving than in the opposite direction. This higher frequency can then lead to wave block and breakup, even in the absence of a steep APD restitution curve (i.e., a curve with a slope smaller than one), as illustrated in Fig. 12. A final example comes from so-called supernormal CV restitution curves, observed in, for example, the His-Purkinje system [218]. In these curves, the CV increases rather then decreases as the DI becomes smaller. Due to this increase in the CV, there is a crowding of spiral waves, which can result in breakup. As with the Doppler effect, this breakup can occur for APD restitution curves with slopes that are smaller than one.

Open in a separate window

Figure 12:

Spiral wave breakup. A-D: Snapshot of a simulation of the FK model in which a spiral wave breaks up due to the Doppler shift that is induced by the meandering tip. This shift results in wave block (B) and spiral wave breakup (C-D). The voltage is displayed using a color map with orange/black corresponding to polarized/depolarized tissue, respectively. From [87].

4.5. Spiral tips and their trajectories

It is clear from our discussion so far, that spiral waves play an important role in arrhythmias. How can we quantify the motion of these spiral waves? This is most commonly achieved by determining the location of the spiral wave tip. This tip can be thought of as the point where the activation front, separating the freshly excited tissue from recovered tissue, and the de-activation front, indicating where the tissue is recovered, intersect. To determine the location of this tip, it is common to construct a phase maps, which shows the progression from depolarization to repolarization [213]. The spiral tip then corresponds to a singularity, which is defined as the point around which the integral of the gradient of the phase ϕ is not equal to zero. Thus, the tip can be determined by computing the line integral around a closed path s→ [213, 219, 220]:

If this closed path contains a single phase singularity, then PS = ±1, depending on the chirality of the spiral wave. Several different techniques can be used to compute these phase maps and the most common one employs applying the Hilbert transform to the time trace of the voltage V(t):

where P is the Cauchy principal value of the integral. This operation, which is equivalent to taking the convolution of V(t) with 1/πt, results in a signal that is phase shifted by π/2. The analytical signal, defined as V(t) + iH(V(t)) = A(t)e^iϕ(t), can then by used to construct the phase, which can be taken to be between −π and +π.

An example of this is shown in Fig. 13A where we have plotted a snapshot of the voltage in a simulation of a single spiral wave. The voltage in this simulation was normalized to be between 0 and 1 and the isocontours of V and H(V) (chosen to be 0.5 for both) are shown in red and green, respectively. The phase singularity, corresponding to the spiral wave tip, is the intersection between these two lines and is indicated as a yellow dot.

Open in a separate window

Figure 13:

Spiral tip dynamics A: Snapshot of a simulation of a single spiral with the voltage shown using a gray scale with white (black) corresponding to depolarized (repolarized) tissue. The intersection of the red and green lines is marked by a yellow dot and represents a phase singularity and can be identified as the spiral wave tip. B: Example of spiral wave breakup in a simulation, with voltage shown using a color scale, activation fronts shown as green line, recovery lines shown as black lines, tips of clockwise rotating spiral waves as black dots and tips of counter clockwise rotating spiral waves as white dots. C&D: Snapshots of counterclockwise rotating spiral waves for two different parameter sets of the FK model. In these snapshots, the voltage is color-coded with red corresponding to high and low corresponding to low values, and the tip trajectory is shown in white. E: Comparison of tip trajectories obtained using the full models (upper row) and the single particle model (lower row). The left and middle column correspond to the trajectories from the FK model shown in C&D while the right column corresponds to a tip trajectory from the KKT model. Scale bars are 0.5 cm. From [221].

It is trivial to extend this computation to a snapshot of a simulation with spiral wave breakup. In these simulations, spiral waves can rotate either counterclockwise or clockwise. This is illustrated in 13B, which shows a simulation of spiral wave breakup in a domain with non-conducting boundaries. In this panel, we have marked the tips of all clockwise rotating spirals as a black dots and the tips of all counterclockwise rotating spirals as a white dots. Furthermore, we have also plotted the activation and repolarization fronts as black and green lines, respectively. For topological reasons both fronts need to start and end at either a phase singularity of a non-conducting boundary. Therefore, for computational domains with periodic boundary conditions, spiral tips always come as pairs, corresponding to a counterclockwise and a clockwise rotating spiral.

Computing the spiral tip as a function of time results in a tip trajectory. This trajectory can take on a wide variety of shapes, depending on the parameters of the model (see, e.g., [222, 223]). For example, the trajectory can be perfectly circular or can meander through space, resulting in a “flower-like” pattern. These two possibilities are presented in Fig. 13C and ​andDD where a change in one parameter of the FK model dramatically alters the trajectory [221]. More complicated trajectories, including ones in which the path contains large linear segments, are also possible.

Even though these trajectories are the results of non-linear PDEs, previous studies have shown that the dynamics of spiral tips can be analyzed in terms of ordinary differential equations near bifurcation points [224, 225] and that simple equations can describe circular tip trajectories in the presence of periodic modulations [226]. Furthermore, the tip trajectories be faithfully reproduced by a single particle (SP) model, which describes the x and y coordinate of the spiral wave tip and, thus, ignores the spatial extent of the spiral wave arms. This is particularly the case for trajectories with dynamics that contain only a few frequencies. In this case, this model describes a single particle subject to periodic forcing terms with frequencies ω₁ and ω₂ and amplitudes F₁ and F₂:

where the phase ϕ determines whether the pattern is inward (ϕ = 0) or outward (ϕ = π) and where the third term represents a damping that prevents the model from drifting and allows it to converge to a steady state quickly. This SP model can be trivially solved and the parameters can be immediately determined from the tip trajectory computed using the spatially extended electrophysiological models [221]. Specifically, the angular frequencies are obtained from the power spectrum of the trajectory, with ω₁ corresponding to the overall period of the spiral and ω₂, corresponding to the angular frequency of the petals in the flower patterns. Furthermore, F₁ and F₂, are determined by the overall size of the trajectory, and the size of the petals.

As shown in Fig. 13E, the SP model can faithfully capture the spiral wave trajectory of full, spatially extended electrophysiological models. In this panel, the trajectory of the full model is plotted in blue while the trajectory of the corresponding single particle model is plotted in red. The left and middle panel show the trajectories of the FK model corresponding to the snapshot in Fig. 13C&D. Both the simple circular trajectory as well as the meandering trajectory (middle column) are reproduced by the SP model. The SP model can also reproduce tip trajectories obtained using more complicated electrophysiological models [221]. An example of this is shown in the right-most trajectory in Fig. 13E, which represents the trajectory in the detailed KKT model [56], already encountered in Fig. 5, and its single particle counterpart (red). Even though this detailed model contains 13 currents and more than 40 parameters, its trajectory can still be accurately described by a single particle subject to two periodic forces. For more complicated trajectories, it should be possible to extend the number of frequencies appearing in the force terms.

Recently, it has been shown that this SP model can also be used to replicate more complex tip trajectories found in spatially extended models in the presence of heterogeneities. These complex trajectories were found in simulations of a single spiral wave in the presence of a pair of equal-sized heterogeneities, in which the conductance was decreased. Specifically, when the conductance inside these heterogeneities was decreased, trajectories in both the simple FK model and the detailed KKT model were found to display chaotic tip trajectories. Simulations revealed that the size of the heterogeneities and their spacing critically determined the tip dynamics [221]. This is quantified in the phase diagram shown in Fig. 14A, which also plots three trajectories for fixed heterogeneity size and different spacing. The phase diagram shows that for both small and large spacings, the spiral tip dynamics is regular: for small spacings, the spiral wave is anchored to both heterogeneities and orbits around them (blue region, Fig. 14A) while for large spacings, the spiral wave rotates around one of the two heterogeneities (purple region, Fig. 14A). For intermediate spacings, however, the trajectory becomes irregular and the spiral tip alternately circulates around one or the other heterogeneity. By computing the Lyapunov exponent of the time series of the tip [227], it was shown that these trajectories were chaotic. These dynamics can also be captured by the SP model when the x and y force equations (Eqns. 21) are supplemented by the terms −dU(x,y)dx and −dU(x,y)dy, respectively. Here, the potential U(x, y), shown in Fig. 14B, represents the two heterogeneities that can be described by two Gaussian wells with circular symmetry and with local minima at the locations of the heterogeneities (c_x1, c_y1) and (c_x2; c_y2):

where g parametrizes the depth of the wells. This extended SP model was able to qualitatively reproduce the results of the full model, as can seen in the phase diagram presented in Fig. 14C [221]. As in the full model, for small spacing, the trajectory migrates around both heterogeneities and for large spacing the particle migrates around one of the two heterogeneities. For intermediate spacings, however, the spiral tip trajectory becomes chaotic with a positive Lyapunov exponent (Fig. 14C).

Open in a separate window

Figure 14:

Complex spiral tip dynamics. A: Phase diagram of the spiral wave dynamics in the FK model in the presence of two identical heterogeneities. Representative tip trajectories, corresponding to the three dots, are shown as black lines, with the location of the heterogeneities indicated as red crosses. The tip dynamics in the yellow region has a positive Lyapunov exponent (λ = 7.1) and is thus chaotic. B: Schematic representation of a particle moving in a potential landscape with two wells, representing tissue heterogeneities. C: Phase diagram of the single particle model in the presence of the potential landscape shown in panel B. The phase diagram is qualitatively similar to the one of the full model, presented in panel A. D-E: Snapshots of a simulation of the FK model in the presence of a circular heterogeneity with a radius of 0.275 cm, plotted by the dashed white line. Snapshot D represents the start of the simulation interval while snapshot E represents its end. The trajectory of the spiral tip (cyan dot) is plotted in red and shows that the spiral wave is temporarily trapped by the heterogeneity. F: Distance r from the tip to the center of the heterogeneity as a function of time. Trapped intervals are indicated by the red bars. G: Trajectory of an intermittently trapped particle, moving on a potential surface with a single well. H: Distance r of the particle to the minimum of the potential well as a function of time. Trapped intervals during are indicated by the red bars. Panel A-C from [221]; panel D-H from [228].

Another example where the single particle model can qualitatively reproduce results from the spatially extended cardiac model equations comes from the intermittent trapping of spiral waves [228]. This trapping was found to occur when a circular heterogeneity with depressed excitability was introduced into a computational domain. This is illustrated in 14D-​-E,E, where the tip trajectory of a spiral wave in the FK model is shown in red. The spiral wave migrates through the computational domain, encounters the heterogeneity and is trapped for several rotations, after which is dislodges and migrates further through the domain. This intermittent trapping can also be seen in Fig. 14F, which plots the distance r from the tip to the center of the heterogeneity as a function of time. For this time interval, the particle was trapped three times, as indicated by the red bars. To replicate these dynamics, a single well potential was added to the SP model (Fig. 14G). As in the full model, the SP model also exhibited trajectories that were intermittently trapped (Fig. 14H). Furthermore, the phase diagram in size vs. heterogeneity strength space was qualitatively similar, indicating that the particle model can be used to simulate spiral tip dynamics [228].

5.1. Required spatial and temporal resolution

What is the minimum required spatial and temporal resolution to capture atrial and ventricular fibrillation in a clinical setting? The ability to determine accurate intracardiac maps of AF from noisy signals depends on the level of contamination, as was recently discussed [247]. For signals that are noise-free, on the other hand, the resolution requirements can be estimated from a general scale analysis [248, 249]. The minimum temporal resolution required to map human VF and AF can be estimated from tissue characteristics. The mapping must be able to distinguish activation on neighboring recording sites, and can be found by dividing the spatial resolution by the conduction velocity of the activation fronts. Typical conduction velocities of cardiac fronts are in the range of 50-150 cm/s [250]. Thus, for a spatial resolution of 5 mm applicable to the basket of Fig. 16D, this results in a required temporal resolution in the range of 3-10 ms. This temporal resolution is easily achievable using clinical mapping systems.

Open in a separate window

Figure 16:

Experimental and clinical mapping of spiral wave patterns. A: Voltage sensitive dye recordings of reentry in in a cultured monolayer of neonatal rat cardiac cells. The colors indicate the normalized voltage level, with blue being the resting state (0%) and red the peak action potential (100%). From [260]. One full rotation of the reentry is illustrated, with frames 20 ms apart. B: Snapshot of the phase map on a rabbit heart surface during sustained fibrillation, obtained using a voltage sensitive dye. The phase map reveals two counter-rotating spiral waves, indicated by the white and black dots. From [213]. C: Activation maps during sustained atrial fibrillation in the RA of a human explanted heart, recorded using a voltage sensitive dye. The activation is simultaneously measured on the endo- and epicardial surface and reveals a reentry pathway on both surfaces. From [261]. D: Isochronal maps of the LA, separated by 237 days and obtained using a basket electrode, demonstrate temporal conservation of a spiral wave in human AF. Time in ms. Targeted ablation at the spiral tip eliminated AF. From [262]. E: Mechanical filaments (red) inside pig ventricles (grey shaded) between endocardial (endo) and epicardial (epi) heart surface during ventricular fibrillation imaged using 4D ultrasound. From [263].

The required spatial resolution depends on the length scale of the mapped event. We will consider two different activation patterns: spiral waves and focal sources. For spiral waves, we can identify three spatial scales. The first scale is the wavelength, λ, defined as the spatial separation between two successive arms of the spiral. If this spiral is rotating in a rigid fashion, we can estimate that to accurately capture this spiral we will need approximately 4 spatial grid points between successive wavefronts, resulting in Δx_req ≈ λ/4. Second, a migrating spiral wave has a length scale given by R_tip, typically smaller than the wavelength. Finally, there is also a minimum required field of view size, reflecting the scale of the coherent domain of the spiral wave, R_rotor. Choosing fields of view that are smaller than R_rotor can miss the spiral tip and rotational activity.

These different length scales have been determined experimentally. The spatial scale of the spiral wave tip, for example, ranges from R_tip = 1 cm to R_tip = 3 cm, requiring a resolution of 0.5 cm to theoretically map the trajectory of spiral tips [251, 252]. The coherent domain of a spiral wave can be estimated to have an area of > 5 cm², corresponding to R_rotor ≈ 2.5 cm [253]. Therefore, the field of view of the mapping technique should be at least 2.5x2.5 cm. The spiral wavelength can be estimated as the product of the minimum conduction velocity and shortest refractory period [254]. In human AF, the former is approximately 40 cm/s [250] while the latter is roughly 110 ms [255], resulting in a minimum wavelength of λ ≈ 4 − 5 cm. Taken together, this means that the smallest spiral waves require a spatial resolution of approximately 1 cm although a larger resolution may be sufficient if the spiral waves are larger. For human VF, the conduction velocity varies within the ventricular wall but is within the same range as AF velocities [256], and the expected spatial resolution requirements are therefore similar to the ones of human AF.

5.2. Optical mapping

Optical mapping is by many considered the gold standard of cardiac imaging. In this technique, voltage sensitive dyes are infused into the heart after which the heart is illuminated with light. The photo-emitting properties of the dyes change as a function of the membrane potential, thus allowing the visualization of cardiac electrical activation and recovery at spatial and temporal resolutions that can easily satisfy the requirements discussed in Sec. 5.1. [257, 213, 258, 259]. These features make optical mapping ideal to map complex activation patterns, including fibrillation.

An example of an in vitro experiment using optical mapping is shown in Fig. 16A [260]. This example shows a series of images of a monolayer of cardiac cells where the fluorescence level was directly converted into a color map. From this map, it is possible to determine the spatio-temporal dynamics of tissue activation, which in this case reveals a clear counter-rotating spiral wave. Importantly, in this experiment, the monolayer was prepared such that it was relatively homogeneous and the observed reentry was therefore functional and not dependent on tissue heterogeneities. A second example of optical mapping is shown in Fig. 16B, where the activation on the surface of a rabbit heart during an arrhythmia is visualized using an optical dye [213]. In this case, the fluorescent data is converted into a phase map, which demonstrates the presence of a spiral wave. A final example shows the activation pattern in an explanted human heart (Fig. 16C) [261]. In this case, the activity was measured on both surfaces of the atrium, revealing a rotating reentry pattern.

5.3. Electrode mapping

Currently, optical mapping techniques are not available in humans due to the toxicity of dyes and other technical problems [264]. It is possible, however, to obtain spatio-temporal information of the electrical activity within the atria and ventricles using an array of contact electrodes. Some of these arrays contain many electrodes with spacing as little as ≈ 1 mm but can only map a relatively small area of the tissue [265, 266]. Other arrays cover a large part of the cardiac chamber but with lower spatial resolution [267, 262]. An example of the latter involves advancing basket electrodes into the atria or ventricles of patients and recording signals from the contact electrodes [214, 267, 262, 268]. This basket consists of 8 flexible splines, each containing 8 electrodes, that are inserted into the atria within a sheath. After removing the sheath, the splines extend and make contact with the atrial wall, resulting a spatial resolution of approximately 5 mm and a coverage of approximately 70% of the atrium, as schematically shown in Fig. 16D. Recordings can be made in both chambers simultaneously or sequentially, after which the electrograms can be analyzed to determine the local activation times. These times can then be used to construct activation or phase maps, visualizing the spatio-temporal dynamics of AF or VF. Importantly, these large electrode baskets meet the spatial resolution requirements discussed above.

The difficulty in interpreting the signals from the contact electrodes, however, is to extract the correct activation time. Often, the electrograms are quite noisy due to, e.g., poor contact with the wall or far field effects. These activation times are required to create activation maps, which visualizes the spatio-temporal dynamics of the cardiac rhythm. Alternatively, they can be used to construct a phase map, This is accomplished by using the activation times to assign a phase as a function of time for each electrode location. Once these phase maps have been constructed, they can be used to compute the location of the singularities, corresponding to the tips of spiral waves (see Section 4 and Eq. 18).

The 64-electrode baskets have demonstrated the existence of spiral waves in both the atria and ventricles of patients in AF and VF [214, 267, 262, 268]. Moreover, for AF patients, it was shown that the tips of these spiral waves are often localized in a conserved location, representing a spiral wave that is spatially stable. An example of spatially conserved activations maps obtained using electrode mapping in a patient is presented in Fig. 16D. These maps were obtained 237 days apart and during AF and revealed a counter-clockwise rotating spiral wave in roughly the same location. Since no procedure was carried out during this interval, these maps suggest that in this patient the spiral wave activity during AF is conserved. In other AF patients, however, tip trajectories were shown to be meandering with a significant spatial extent [262].

We should note that these observations are not completely without controversy. Several reports, using high resolution electrode arrays, have failed to observe consistent rotational activity [266, 269, 265]. There are a number of possible explanations for these disparate results. First, these reports typically record activity patterns on the epicardial wall, which may be more complex than endocardial patterns [270]. Second, the spatial extent of these arrays are small, making it more difficult to capture large scale rotational activity. Third, the spatio-temporal patterns are analyzed using activation times at the site of each electrode. As discussed before, the electrical recordings are noisy and assigning activation times is challenging. Incorrect identification of the activation times results in flawed activation maps, which may hinder the identification of rotational patterns.

6.1. Drug therapy

Drugs that target ion channels can be used to treat patients with cardiac arrhythmias. In the case of AF, for example, patients can be given drugs that either aim to restore and maintain sinus rhythm or that attempt to control the ventricular rate [277, 278]. Unfortunately, the success rate of these pharmacological interventions is relatively low. The problem is that it is challenging to predict how changing the properties of a specific ion channel may affect the spatio-temporal dynamics of heart tissue. What might be antiarrhythmic at the single cell level could be pro-arrhythmic at the organ level. A well known example of this comes from the Cardiac Arrhythmia Suppression Trial (CAST) trial, which was conducted in the 1980s [279]. This trial tested whether so-called class I drugs, drugs that block the fast sodium channel, would reduce the occurrence of VF. Single cell experiments had shown that these drugs were able to suppress spontaneous electrical activity and the hypothesis was that fewer extra stimuli would reduce the risk of fibrillation [280]. Unfortunately, the trial actually showed the opposite and had to be stopped prematurely [279, 281]. It was not until 2011 that a modeling study pinpointed the reason why these drugs were pro-arrhythmic [60]. This study showed that even though these sodium channel blockers may reduce spontaneous electrical activity they can also increase the size of the vulnerable window and, thus, lead to unidirectional block, as discussed in Section 4. This unidirectional block can then result in retrograde conduction and the initiation of reentry and, possibly, fibrillation.

6.2. Termination of fibrillation using shocks

VF is lethal within minutes and the only chance of survival for a patient is defibrillation using a defibrillator. This device creates a large electrical field across to the heart, after which the heart is hopefully be able to restart its normal sinus rhythm [282, 283]. Defibrillation can be applied using an external defibrillator or a implantable cardioverter defibrillator (ICD). The latter, used in patients that have a high risk of developing VF, can automatically detect a ventricular arrhythmia. Repeated shocks are, however, associated with an increase in mortality [284] and to limit the amount of shocks, ICDs apply an antitachycardia pacing (ATP) protocol to terminate detected VTs. In this case, the device delivers a train of stimuli to the heart and the idea is that the activation waves generated by the stimuli interfere with the tachycardia rhythm and take control of the tissue [285, 286]. Several different ATP algorithms are currently in use [287] and a recent study used personalized computational modeling to determine the efficacy of two of these algorithms [288].

The precise mechanism by which defibrillation works has long been a topic of discussion but it is clear that it has to affect the membrane potential far away from the electrodes of the defibrillating device. Applying an electric field in homogeneous tissue, however, would result in changes in the membrane potential that decay exponentially with distance x from the electrodes: V ~ exp(−x/λ) [289]. Here, λ is the space constant of cardiac tissue, which is on the order of 1 mm. Thus, the applied electric field decays quickly, will be essentially negligible several millimeters away from the electrodes, and will not be able to affect the membrane potential in most of the tissue.

Over the last few decades, insights from modeling studies have provided plausible defibrillation scenarios. Crucial in these scenarios is the fact that cardiac tissue is not homogeneous but consists of cardiac cells that are coupled through gap junctions that have a higher resistance than the rest of the cell. The resulting heterogeneities, on the scale of the cell, can lead to oscillatory variations in the membrane potential that do not decay away from the defibrillator electrodes [290, 291, 292, 293]. This scenario is shown in Fig. 17A, where the membrane potential is computed in a cable comprised of cells with a high junctional resistance. Due to this junctional resistance, the membrane potential is not smoothly varying but, instead, has a characteristic sawtooth shape. For a sufficiently large shock, these oscillations can excite quiescent cells, which leads to additional propagating activation waves, even far away from the stimulus site. These waves can then interact with the irregular activation responsible for VF, leading to termination of all propagation.

Open in a separate window

Figure 17:

Low-energy defibrillation. A: Sawtooth potential in a simulation of a strand of 30 cells and an intracellular junction resistance of 50 MΩ (From [290]). B: Upper panel: Unipolar anodal stimulation at the center of a thin computational sheet with equal intra- and extracellular anisotropy ratios. In this case, the changes in V_m have identical sign, with a shape that is given by the the fiber orientation (taken here to be along the horizontal direction, as indicated by the arrow). Lower panel: Identical stimulus but now with unequal intra- and extracellular anisotropy ratio. Changes in V_m are both positive and negative, with tear-drop shaped regions of depolarization (called virtual cathodes). (From [122].) C: Example of termination of spiral waves by wave emission from heterogeneities in a simulation using the Luo-Rudy model. Electric field pulses initiated waves from the obstacles, resulting in the termination of the rotating waves. (From [303].) D: Snapshots of the electrical activation of a dog’s atrium, visualized using optical mapping, during low-energy defibrillation (From [304]). The upper row visualizes AF while the lower row corresponds to the activation maps after each low energy pulse (LEAP). The last panel shows the removal of fibrillation, resulting in tissue that is at rest.

The amount of energy delivered by a single shock from an external and an implantable defibrillator is approximately 150-360 J and 10-30 J, respectively. This large amount of energy is above the pain threshold for humans and can have deleterious effects [294]. Therefore, it would be advantageous to implement low-energy defibrillation techniques that are below the pain threshold and that are not harmful to the patient (for a recent review, see [295]). One of these techniques rely on the concept of so-called virtual electrodes [296, 297, 298, 299, 283]. These virtual electrodes arise from the difference in anisotropy between the intracellular domain (the cardiac cells) and the extracellular domain (the space surrounding the cells).

The appearance of virtual electrodes can be shown by simulating the bidomain equations, which describe the potential in the intracellular and in the extracellular domain (see Eqn. 14 in Section 3). The effect of the anisotropy ratio is illustrated in Fig. 17B. The upper panel shows the activation pattern following anodal stimulation in a piece of tissue with equal anisotropy ratio. This stimulation results in an ellipsoidal hyperpolarized region indicated by AN (anode) that will not result in tissue activation. Note that the observed pattern would be identical to the one obtained using the monodomain equation formulation (see Section 3).

For unequal anisotropy ratios, however, the bidomain equations will give a different polarization pattern for the same stimulus since the potential in the intra- and extracellular space will evolve differently. In particular, close to the stimulus, it is now possible to have regions that are both depolarized and hyperpolarized. This is shown in the lower panel of Fig. 17B: the anodal stimulus results not only in a hyperpolarized region (again labeled as AN) but also in two depolarized regions. These regions are called virtual electrodes (and in this case virtual cathodes and labeled as VC since the regions are depolarized) and can, for strong enough stimulation, result in tissue activation even in the case of an anodal stimulus. These virtual electrode patterns can occur during both onset and termination of the stimulation and for different electrode configurations and have been verified in multiple optical mapping experiments [300, 301, 302].

The fact that virtual electrodes can hyperpolarize or depolarize regions near a stimulation electrode can play a significant role in the success or failure of a defibrillating shock [299, 305]. Furthermore, the creation of virtual electrodes near obstacles can be used to lower the required energy for defibrillation. Specifically, simulation studies showed that the application of an electric field can induce regions of hyperpolarization and depolarization near an anatomical heterogeneity where the conductivity changes in a discontinuous manner. These studies also showed that the resulting activation sites were able to unpin a rotating wave anchored to this heterogeneity [306, 307]. Subsequent experimental work demonstrated that the required energy to unpin is significantly smaller than a defibrillating shock [308].

Pumir and Krinsky provided an analytical approach that showed how the depolarization or hyperpolarization near an obstacle depends on the applied electric field E₀ [306]. They linearized the monodomain equation around the resting potential, V_r, to obtain

where e = V_m − V_r, V_m is the membrane potential, and λ is the already encountered space constant. For a circular obstacle with radius R, it is possible to find a solution using the appropriate boundary conditions: e(r) → 0 when ∣r∣ → ∞ and ∇(e + E₀ · r) · n = 0 on the boundary of the obstacle. The latter condition means that the current across the obstacle’s boundary is zero, with n being its normal. The solution of this problem, expressed in polar coordinates (r, θ), is written in terms of Bessel functions:

Thus, the maximum value of e is achieved at the boundary of the obstacle and increases for increasing values of R.

The scenario above focused on the appearance of a virtual electrode from a single obstacle. If the domain contains multiple discontinuities in conductivity, then a virtual electrode will appear at each of these discontinuities. For a sufficiently strong electric field, these virtual electrodes depolarize the tissue above the excitation threshold and can act as sites of activation. In that case, the additional activation waves can excite large portions of the domain, potentially extinguishing any existing propagating waves. Indeed, computational studies showed that the combination of small scale heterogeneities and an low energy electric field can effectively defibrillate cardiac tissue [303]. An example of this is shown in Fig. 17C, where spiral wave breakup was created in a square computational domain that contains multiple small and large non-conducting obstacles (panel 1). Electric field pulses resulted in the initiation of waves from the heterogeneities (panel 2 and 3) and the interaction with these waves and the original spiral waves led to the termination of any activation (panel 4).

Experimental studies, both in vitro and in vivo have tested the idea of low-energy defibrillation based on heterogeneity-induced virtual electrodes [309, 304]. These studies found that it is possible to terminate atrial and ventricular fibrillation using a limited number of shocks with an energy that is approximately 10% of the energy of a standard defibrillating shock. An example of the in vivo work is presented in Fig. 17D, which shows sequential snapshots of the optical mapping of AF in a dog before (upper row) and during (lower row) the delivery of low energy electric pulses. The pulses were able to create additional activation sites as is evident from the increase in size of the excited area in the lower row. After the last pulse, all activation terminated and the tissue became quiescent. Although a number of structural properties can act as the heterogeneities required for the generation of additional activation sites, including abrupt changes in the fiber orientation and fibrosis, the study by Luther et al. suggested that blood vessels are sufficient to explain the appearance of the additional activation sites [304].

6.3. Ablation treatments

It has become increasingly common to use ablation to treat a variety of cardiac arrhythmias. In this procedure, cardiac tissue is destroyed and rendered non-conductive, typically using either radiofrequency ablation, which burns cells, or cryoablation, which freezes cells. Ablation is used to treat relatively “simple arrhythmias”, including supraventricular tachycardia and ventricular tachycardia but also the most complex heart rhythm disorder, fibrillation.

6.3.2. Ablation of AF: targeted ablation

One such alternative strategy uses ablation to target the sites of rotational or focal activity [214]. Using the electrode basket described in Section 5, it is possible to pinpoint the location of the tips of spiral waves (see Fig. 16D) [267]. These tip locations can then be targeted using localized ablation, rendering the area of the spiral wave tip non-excitable. This technique was shown to be able to abruptly terminate AF in some patients [214, 267, 314]. Furthermore, this procedure also reported better long term outcomes than PVI ablation, although outcomes vary be center [214, 315]. The observed termination was perhaps surprising since it was thought that the creation of a non-excitable region may anchor spiral waves, resulting in waves that indefinitely rotate around the ablation spot. Using computational analysis, several possible mechanisms have been proposed to explain the success of this targeted ablation [316]. For example, it is possible that the ablation lesion connects the tip location to a non-excitable boundary, thus creating a line of block that prevents the spiral wave to continue its rotation. It is also possible that ablation removes a particular reentry path within the atrium, as shown in Fig. 18A-​-D.D. In this case, an isthmus between two nonconducting regions, indicated by the hatches in the figure, provides a conduction pathway and ablating this pathway will remove the reentry.

Open in a separate window

Figure 18:

Modeling ablation of fibrillation. A-F: Examples of ablation mechanisms that lead to acute termination of atrial fibrillation (Adapted from [316]). In (AD), a slow-conducting zone is sandwiched between nonconducting regions (hatches). The resulting figure-of-eight reentry can be abruptly terminated by ablating the slow-conducting isthmus. E: Simulation of termination in a high excitability region, indicated by the dashed circle, that is surrounded by low excitability tissue. The excitability is a linear function of the radial direction. Before ablation (left panel), a wave is attached to a circular obstacle. After ablation, the wavefront finds itself in a region of low excitability, detaches from the obstacle, and migrates towards the edge of the computational domain, resulting in the termination of reentry (right panel). F: Simulation of termination in a region of low excitability. Within the dashed circle, the tissue excitability is a linear function of the radius and is lowest at the obstacle. Before ablation, a counter clockwise rotating activation front is attached to a nonconducting zone (left panel). After ablation (middle and left panels), the wavefront encounters the wave back, resulting in block (double bars in right panel). The wavefront detaches from the ablation lesion and terminates. Scale bars, 1 cm. (G-H) Targeted ablation based on fibrotic tissue regions (From [211]). (G) Posterior (top) and anterior (bottom) views of a patient-specific atrial mode as reconstructed from segmented LGE-MRI scans, including the distribution of fibrotic tissue shown in green. (H) Sites of ablation delivery shown as red dots in the left atrium of the patient in (G). The dashed ellipse indicates the locations ablated based on the locations of persistent reentry circuits, identified by computational analysis.

In addition, spatial heterogeneities in excitability in combination with ablation lesions, can play a role in the termination of reentry. In Fig. 18E&F we show two examples of this scenario where the excitability parameter in the Fenton-Karma model (τ_d) is varied in a radial fashion. In the first, a spiral wave is stably attached to a nonconducting zone and the excitability is higher close to this zone than in the surrounding tissue. Ablating the heterogeneous high-excitability domain results in a nonconducting zone, detachment of the spiral wave, and consequent migration and collision with anatomic obstacles, which are modeled as no-flux boundaries, and termination. In the second scenario, illustrated in Fig. 18F, a spiral wave is located in a region of low excitability surrounded by higher excitability. Initially, the activation front is anchored at a small, circular nonconducting obstacle. Creating an ablation lesion of larger radius results in the activation front propagating in tissue with higher excitability and correspondingly higher wave speed, allowing the activation front and back to meet, leading to wave block (double bars in the right panel of Fig. 18F) and termination.

It is also possible to determine the effect of localized ablation using patient-specific atrial geometries that include areas of fibrosis, as was carried out in a recent computational study [211]. These fibrotic areas were determined using late gadolinium enhancement magnetic resonance imaging (LGE-MRI) [210]. This technique can visualize gadolinium-based contrast agents within the extracellular space of the myocardium, although its reproducibility and accuracy are still under investigation [317]. An example of the results of this imaging technique in a patient is shown in Fig. 18G, where the fibrotic regions are plotted in green. This computational study then applied virtual ablation targeted to these sites of fibrosis and continued to add additional lesions until AF was no longer inducible [211]. The lesions for this particular patient are shown in Fig. 18H as red dots and connect the region of fibrosis with a non-conducting boundary. As a proof of concept, this procedure was prospectively tested in a clinical setting with a limited number of patients, showing promising results [211].

6.4. Control of fibrillation

A number of alternative options have been proposed to treat fibrillation. Most of these options are based on theoretical and computational studies and have not yet been implemented in clinics. One proposed treatment option attempts to use control theory to prevent alternans from occurring. As discussed in Section 4, alternans manifests itself in APDs or Ca²⁺ signals that alternate between long and short and can also materialize at a spatial scale. Focusing on electrical alternans, we can write the restitution curve (see Fig. 5) as APD_n+1 = f(APD_n, T), where APD_i is the action potential duration of the i-th beat and where T is the pacing period. For restitution curves with a slope larger than one, we saw that the fixed point of this recurrence map is unstable. Employing a feedback control scheme, however, that changes the pacing period in a judicious way, makes it possible to drive the system to its stable fixed point. For example, replacing the pacing period at each cycle to T−α2(APDn−APDn−1), results in the following 2D map

The fixed point of this system, APD* = Y* is stable as long as

where A=∂f∂APD and B=∂f∂T [336]. This and other control schemes have been shown to be successful in computational and small scale in vitro systems although adding the dynamics of Ca²⁺ may make it more difficult to achieve adequate control [336, 337, 338, 339, 340, 341, 342].

The heart is, of course, a spatially extended system and the control of spatio-temporal instabilities arising from alternans is more challenging [343]. A computational study investigated the break up of a spiral wave in the presence of equally spaced control points [344]. It showed that applying a feedback current at these discrete points during the repolarizing phase of the action potential could suppress alternans and, therefore, prevent spiral wave break up. It also determined that the success of the control scheme critically depends on two factors: the amplification time of the instability of spiral wave break up and the electrode spacing. For successful control, the amplification time needs to be larger than the time required for control and the the electrode spacing must be comparable or smaller than the spiral wavelength. These factors, in addition to other practical issues, make it difficult to apply this control scheme to in vivo settings. Interestingly, a recent study showed that a proprietary closed-loop system was able to significantly suppress the amplitude of alternans voltage in swine hearts [345]. Furthermore, this study showed a trend toward reduced arrhythmia susceptibility in the presence of control. It would be interesting to further understand the mechanisms that underlie these results.

This work was supported by National Institutes of Health R01 HL122384, RO1 HL149134, and HL083359. I would like to thank Charlotte Rappel, Timothy Tyree and Michael Reiss for a careful reading of the manuscript.