Multifractal Desynchronization of the Cardiac Excitable Cell Network During Atrial Fibrillation. II. Modeling

In a companion paper (I. Multifractal analysis of clinical data), we used a wavelet-based multiscale analysis to reveal and quantify the multifractal intermittent nature of the cardiac impulse energy in the low frequency range ≲ 2Hz during atrial fibrillation (AF). It demarcated two distinct areas within the coronary sinus (CS) with regionally stable multifractal spectra likely corresponding to different anatomical substrates. The electrical activity also showed no sign of the kind of temporal correlations typical of cascading processes across scales, thereby indicating that the multifractal scaling is carried by variations in the large amplitude oscillations of the recorded bipolar electric potential. In the present study, to account for these observations, we explore the role of the kinetics of gap junction channels (GJCs), in dynamically creating a new kind of imbalance between depolarizing and repolarizing currents. We propose a one-dimensional (1D) spatial model of a denervated myocardium, where the coupling of cardiac cells fails to synchronize the network of cardiac cells because of abnormal transjunctional capacitive charging of GJCs. We show that this non-ohmic nonlinear conduction 1D modeling accounts quantitatively well for the “multifractal random noise” dynamics of the electrical activity experimentally recorded in the left atrial posterior wall area. We further demonstrate that the multifractal properties of the numerical impulse energy are robust to changes in the model parameters.

1. Introduction

Atrial fibrillation (AF) is the most common sustained tachyarrhythmia encountered in clinical practice (Nattel and Harada, 2014). It is sometimes not diagnosed until the occurrence of a severe complication such as embolic stroke. Often associated with heart disease, clinical investigations do not always uncover any preexisting cardiovascular comorbidity (idiopathic or lone AF). Physiologically, current understanding of the onset and perpetuation of most tachyarrhythmias including AF presumes the involvement of circuit reentries. This scenario was established historically from the observation of reciprocating rhythms initiated in the atrio-ventricular node via the fast and slow pathways of impulse conduction. Atrial flutter and regular tachycardias were thus inferred to be rooted in circling conduction pathways, as going around anatomical obstacles or scars for example. Tachyarrhythmias may also spontaneously evolve toward more irregular arrhythmias such as AF. Experiments show that AF can be induced by in situ injection of toxins like aconitine as well as by ectopic stimulation, i.e., under extreme conditions enforcing local functional changes of the excitable conducting substrate. AF may then persist independently of the inciting protocol (Macfarlane et al., 2011; Zipes et al., 2017). These observations paved the way to the concept of multiple circuit reentries, not necessarily linked to the anatomy but to a vulnerable atrial substrate because of functional dispersion in space and time (such as non uniform dispersion of refractoriness) (Moe and Abildskov, 1959; Moe et al., 1964; Allessie et al., 1977; Attuel et al., 1989; Winfree, 1989). But clinically, the question remained whether abnormal conducting pathways and ectopic triggers do stabilize AF. In that respect, intervention procedures were developed either to surgically create an electrical maze in the atria or, in a less invasive and safer way, to isolate abnormal ectopic activity as found in the pulmonary veins areas by radio frequency ablation. Both procedures led to high clinical success rates in stopping paroxysmal AF (Cox et al., 1991; Ha¨ıssaguerre et al., 1998). Unfortunately, the diverse procedures instigated since then remain suboptimal because the risk of relapse increases with time, and the disease often evolves toward a chronic state (Ganesan et al., 2013; Takigawa et al., 2014).

Cardiomyocytes belong to the family of excitable cells which are ubiquitous in animals and plants (Hille, 2001). They are distinguishable from non-excitable cells by their ability to reach an electrically depolarized state of their extra-cellular phospolipid bi-layer membranes. Action potentials (APs) correspond to cycle events in which the membrane reaches a depolarized state before relaxing back to the polarized rest state. In the wake of the seminal work by Hodgkin and Huxley on the giant squid axon AP (Hodgkin and Huxley, 1952a,b), a cardiac impulse is similarly described by the nonlinear coupling between a diffusing activator, the electric potential across the excitable cell membrane, and a non-diffusing inhibitor, the overall ion currents permeating through the membrane (Noble, 1962, 1965; Beeler and Reuter, 1977; Plonsey and Barr, 2007; Fenton and Cherry, 2008; Macfarlane et al., 2011). This nonlinearity underlies the fact that the AP amplitudes depend very little on the intensity of the exciting perturbation, provided they are suprathreshold. Various transmembrane proteins selectively allow some solutes to permeate. Leaking (potassium) channels are balanced by (sodium-potassium) pump exchangers forcing the cell membrane into a negatively polarized state, which compensates for the hypertonic activity of internally sequestered vital substances (Tosteson and Hoffman, 1960; Armstrong, 2003). Excitable cells take advantage of this situation to generate electrical signals. Their plasma membrane incorporates a large number of ion channels, sensitive to various other species such as e.g., calcium. They are proteins forming pores that greatly facilitate ion transport down electrochemical gradients. Ion channels act as voltage dependent gates and their reaction rates reflect the height of the free energy barrier separating the open and closed conformation states (Hille, 2001). The membrane depolarizes in a few milliseconds to a near Nernst-Planck resting equilibrium, as for instance triggered by a supra-threshold electrical stimulus. In the heart, in addition, each cardiomyocyte cycle lasts a definite amount of time, typically a few hundreds milliseconds, incorporating a refractory period during which re-excitation is impossible.

Electric pulses travel by ohmic conduction within a single cell membrane, whereas the transport of these pulses across adjacent cardiac cells is ensured by gap junction channels (GJCs). These GJCs connect adjacent cardiac cells along a preferred longitudinal direction, each via tight assemblies of hundreds to thousands of pores, themselves formed by two hemichannels of six bound proteins called connexins. Also found in nearly every connected animal tissues, and on the contrary to ion channels, GJCs are wide and open at rest, thereby prominently contributing to homeostasis, chemical messaging, and electrical permeability throughout the whole network of connected cells (Weidmann, 1966; Severs, 1990; Kumar and Gilula, 1996; Harris, 2001; Evans and Martin, 2002). The reduction of connexin genetic expression was shown to be source of conduction inhomogeneities increasing susceptibility to arrhythmias in ventricles and atria, in various animal models and in humans with diverse heart conditions, potentially causing sudden cardiac death (van der Velden et al., 2000; Dupont et al., 2001; van der Velden and Jongsma, 2002; Ausma et al., 2003; Danik et al., 2004; Severs et al., 2008). Recent studies also suggest that missense mutations in connexin encoding genes predispose to AF (Gollob et al., 2006), whereas connexin gene transfer plays a protective role in preventing sustained AF (Bikou et al., 2011; Igarashi et al., 2012).

The phenomenology of irregular arrhythmias is classically interpreted as emanating from the chaotic dynamics of excitable reaction-diffusion systems in which ohmic conduction is assumed. Typical routes to deterministic chaos with specific rhythms were theoretically identified in simple models based on cell cycle phase resetting (Glass and Mackey, 1979, 1987; Guevara et al., 1981; Guevara and Glass, 1982). In particular, period doubling bifurcations in the dynamics of excitable pulses propagating on a ring (Courtemanche et al., 1993) were put forward as a mechanism prior to the onset of AF. More complete ion channel models such as the historical Beeler-Reuter model (Beeler and Reuter, 1977), were shown to succeed in generating generic spatio-temporal patterns (Jensen et al., 1984; Chialvo et al., 1990; Fenton and Cherry, 2008). Spiral waves in 2D spatial dimensions (theoretically scroll waves in 3D) were observed experimentally during ventricular tachycardia as well as before the onset of fibrillation (Pertsov et al., 1993; Gray et al., 1995; Garfinkel et al., 1997; Zipes et al., 2017). These spiral waves can be seen as focal wave trains swirling periodically around the analog of the “leading circle” at their core (Allessie et al., 1977; Winfree, 1989). More recently, cellular automata on a two-dimensional square lattice was used to demonstrate that phase-dependent spiral attenuation could reproduce wave propagation in excitable media of myocardial cells (de la Casa et al., 2007a,b). Within this framework, fibrillation was interpreted as the break-up of unstable spiral waves spoiling rhythmic regularity (Gerhardt et al., 1990; Ito and Glass, 1991; Bär and Eiswirth, 1993; Karma, 1993; Panfilov and Hogeweg, 1993; Bär and Brusch, 2004; Zipes et al., 2017). Other models based on nonlinear stochastic feedback mechanisms were also proposed to explain the regulation of cardiac dynamics (Ivanov et al., 1998), suggesting that it could be intrinsically random (physiological noise). To our knowledge, it has never been reported in cardiac models the existence of a multifractal intermittent dynamics of the cardiac impulse energy as observed experimentally over large time scales in our companion paper I (Attuel et al., 2017). In this study, we elaborate on a tentative interpretation of the observed intermittent dynamics during AF as the signature of synaptic plasticity. Typical individual GJC transition times between open and closed states were shown to be much longer than those of membrane polarization but compare well with membrane recovery time (≳ 100ms) (Spray et al., 1984; Neyton and Trautmann, 1986; Wang et al., 1992; Bukauskas and Verselis, 2004; Desplantez et al., 2007). Moreover, slow gating modulations have been evidenced due to cytoplasmic protons (low pH) and free calcium (Spray et al., 1984; Burt and Sray, 1988; Kumar and Gilula, 1996; Harris, 2001; Bukauskas and Verselis, 2004; Perrachia, 2004; Swietach et al., 2013). Thus, a perturbation of the GJC opening and closing due to electric charge may induce some time lag or advance in the activation of the cell, slowing down or boosting the propagation of the AP, even impeding or reversing it (after reexcitation), resulting in a local departure from ohmic conduction law.

Here, we propose a mathematical model of cardiac cell excitability which includes their dynamical coupling by GJC kinetics. As (i) gap junctions electrically bind cardiac cells preferentially along their elongated direction (Severs, 1990; Evans and Martin, 2002)), and (ii) in the left atrial lateral wall area, the CS has a thin surrounding muscular structure traversed by myocardial strands (Ho et al., 2012), we simply consider a one-dimensional (1D) spatial model to describe the transport of AP along and across myocardial cells via the temporal interplay of voltage-gated channels and GJCs. We show that (if probably not minimal) this 1D model robustly accounts for the intermittent modulation of cardiac pulse trains experimentally observed in the clinical data recorded in the left atrial posterior wall area of the CS (The study reported in Companion paper I (Attuel et al., 2017) was carried out with the recommendations of the International Cardiac Electrophysiological Service of public hospital CHU Haut-Lévêque, Pessac, France. The protocol for clinic research was approved by the Institutional Clinical Research and Ethics Committee: CPP (Comité de Protection des Personnes) and AFSSaPS (Agence Française de Sécurité Sanitaire des Produits de Santé). All patients involved gave written informed consent to the investigation of data. For this specific investigation of the data, the authors accessed fully anonymized and de-identified data.

2. Model and Numerical Data

To reproduce the spatio-temporal multifractal intermittent dynamics of voltage signals collected from the CS of patients with chronic AF (companion paper I Attuel et al., 2017), we propose the following system of four nonlinearly coupled partial differential equations (PDEs) (Attuel et al., 2016):

$\begin{array}{ll}\{\begin{array}{c}{c}_m\frac{\partial }{\partial t}{U}_m=\text{F}(U_m,w_m)-c_g^{-1}\frac{\partial }{\partial x}(g{\rho }_g)+\kappa \frac{{\partial }^2}{\partial x^2}{U}_m,\hfill \\ \frac{\partial }{\partial t}{w}_m=\text{G}(U_m,w_m),\hfill \\ \frac{\partial }{\partial t}g={\omega }^2{\rho }_g-{\nu }_{1}g,\hfill \\ \frac{\partial }{\partial t}{\rho }_g=-g\frac{\partial }{\partial x}{U}_m-{\nu }_2{\rho }_g,\hfill \end{array}\hfill & (1)\hfill \end{array}$

where U_m(x, t) is the membrane electric potential drop, w_m(x, t) the ionic channel gating variable, g(x, t) the gap junction conductance deviation from normal, and ρ_g(x, t) the gap junction capacitive charge density. w_m(x, t) is generally a vector of state variables describing the generation of the APs that is related to the intracellular concentration variations of different ions (Na⁺, K⁺, Ca²⁺, Cl⁻). As explained in section 2.1, the kinetics of the voltage gated channels will be described by the simplified FitzHugh-Nagumo model (FitzHugh, 1962; Nagumo et al., 1962; Izhikevich and FitzHugh, 2006). Let us note that, as compared to previous modeling attempts, our model (Equation 1) lies on the assumption that the gap junction conductance is not static, and that the GJC gating is driven by both the local transmembrane field produced by the charging of the gap junction and its conductance (nonlinear coupling term). The first two equations correspond to the standard mono-domain model for cardiac AP conduction, in which an additional term for the GJC current has been introduced, responsible for an imbalance between depolarizing and repolarizing membrane currents. The last two equations describe the GJC conduction and capacitive charging.

2.1. Standard Modeling of Cardiac AP Conduction

Let us consider an idealized elongated fiber of excitable cells (Figure 1), along which a traveling depolarization front (AP upstroke) is classically modeled by a 1D cable equation, assuming Kirchhoff's law of conservation of currents (Plonsey and Barr, 2007; Niebur, 2008; Macfarlane et al., 2011):

$\begin{array}{ll}\{\begin{array}{c}{c}_m\frac{\partial }{\partial t}{U}_m=\text{F}(U_m,w_m)+\kappa \frac{{\partial }^2}{\partial x^2}{U}_m,\\ \frac{\partial }{\partial t}{w}_m=\text{G}(U_m,w_m),\hfill \end{array}& (2)\end{array}$

where U_m (in V) is the electric potential across the cell membrane, c_m (in F/m) is the fiber's insulating membrane capacitance per unit length, κ = σS (in Ω⁻¹ × m) is an inverse resistance per unit length with σ the mono-domain conductivity¹, and S (in m²) is the fiber cross-section. The coupling of the membrane electric potential U_m with inhibiting membrane currents responsible for repolarization by ion specific voltage-gated channels (Na⁺, K⁺, Ca²⁺, Cl⁻) is represented by the variable w_m (in A/m) in the nonlinear function F. The function G symbolizes the kinetics of the repolarizing voltage gated ion channels. One example of a boundary condition consists in imposing an external current stimulus I_ext (in A), possibly time dependent, at the x = 0 extremity of the cell array:

$\begin{array}{ll}\kappa \frac{\partial }{\partial x}{U}_m(x=0)=-I_{ext}.& (3)\end{array}$

A popular model for the nonlinear function F(U_m, w_m) in Equation 2 is the so-called FitzHugh-Nagumo (FHN) model (FitzHugh, 1962; Nagumo et al., 1962; Izhikevich and FitzHugh, 2006), which was constructed from a damped van der Pol oscillator model (named the Bonhoeffer-van der Pol model by Fitzhugh). Indeed, the FHN model was introduced as a nerve model to simplify the Hodgkin-Huxley model (Hodgkin and Huxley, 1952a,b) and to facilitate analytic calculations:

$\begin{array}{ll}\text{F}(U_m,w_m)=\mu U_m(A-U_m)(U_m-B)-{\gamma }^2{w}_m,& (4)\end{array}$

where μ, A, B and γ are real parameters. Note that c_m/(ABμ) ≈ τ_d corresponds to a depolarization time scale, A and B are coefficients establishing the magnitude of the saturating “plateau” potential while the sign of A tells whether the portion of fiber is excitable (A > 0) or self-sustained oscillating (A < 0), and γ² > 0 is the coupling rate between U_m and w_m. The coupling to membrane currents is described by the following form of the function G in Equation (2):

$\begin{array}{ll}\text{G}(U_m,w_m)={\alpha }^2{U}_m-{\nu }_0{w}_m+\lambda ,& (5)\end{array}$

where α² > 0 (real positive) is an inverse inductance per unit length, λ is a leaking current parameter that will be put to zero for simplicity, and ${\nu }_0^{-1}$ a relaxation time scale for voltage gated channel inactivation. Since Equation (2) is a RLC circuit analog, it is straightforward to show that the saturating plateau potential lasts for a typical refractory period RP ≈ (αγ)⁻¹, similarly to the damped van der Pol oscillator (Van der Pol and Van der Mark, 1928; Takashi, 2007).

FIGURE 1

Figure 1. (A) Microscopic view of cardiac muscle cells (longitudinal section) stained with hematoxylin (nuclei) and eosin (cytoplasm). The myocardial cells can be recognized by their elongated aspect (50–100 μm long, with sections ~ 10μm), a longitudinal striated organization, multiple branchings and connections at their extremities via intercalated disks (IDs). (B) Schematic GJCs at the onset of depolarization. Depolarizing (resp. polarized) cells are stained in red (resp. blue). Surface charges are depicted with ± symbols, on the inner and outer side of the depolarizing (gray) and polarized (black) membranes and on the dysfunctional GJCs (red) with capacitive charge loading Q_g such that 〈Q_g 〉 = ∫ ρ_g dx. Voltage-gated channel ionic flows are marked with green vertical arrows. Normal GJCs: red resistor symbol, dysfunctional GJCs: gray capacitor symbol. Total currents I (blue horizontal arrow) circulate in opposite directions inside and outside the cell, splitting into the normally flowing ohmic current I_Ω (red), and residual closed GJCs current I_g (dashed red) building up charge. Typical membrane and intercalated disk dimensions are l_m ≈ l_g ∽ 5 nm.

It takes a time τ_d ∽ 1 to 10 ms to depolarize one cell membrane. Actually, τ_d acts as a cut-off time scale for continuous models as Equation (2). For a linear capacitance $c_m\lesssim 10^{-6}$ F/m, and a linear conductivity κ ≲ 10⁻¹⁰ Ω⁻¹ × m, an upstroke (ion channels releasing the electric energy stored by conversion of chemical bond energy) is estimated to travel at a typical velocity $c\backsim \sqrt{\kappa {(c_m{\tau }_d)}^{-1}}\gtrsim 0.1$ m/s. This is consistent with the conduction velocity measured in the atria under various conditions c ∽ 0.1 to 1 m/s (Zipes et al., 2017). The upstroke spans a typical membrane length of order l_d ≈ cτ_d ∽ 0.1 to 10 mm (Figure 1B), which acts as a spatial cut-off scale in such continuous models. Note that for very slow conduction velocity and rapid upstroke, the upstroke length scale can decrease down to the longitudinal size of one cell. Mathematically, a single boundary stimulus gives rise to a unique traveling upstroke: U_m(x, t) = U_m(x − ct). Immediately following depolarization, other voltage gated ion channels start contributing in a way as to set the membrane potential back to its original negative value. The resting polarized state is then restored back in a typical time of order τ_r ≳ 100 ms in the heart (Hille, 2001; Plonsey and Barr, 2007; Macfarlane et al., 2011). Time history dependencies of channel activation define the absolute refractory period RP of duration ≳ 100 ms in the heart.

A traveling pulse solution of the partial differential (Equation 2) indeed corresponds to a homoclinic cycle in the FHN model biasymptotic to the resting electric potential stable fixed point (Beuter et al., 2003; Izhikevich, 2007; Fenton and Cherry, 2008; Guckenheimer and Kuehn, 2009). Homoclinic cycles are triggered locally by the advancing foot of the depolarizing front or by an external stimulus. To reproduce the observed variety of physiological properties of an excitable myocardium, nonlinear modifications or extra membrane currents were added to the coupling functions F and G as for instance in Bär and Eiswirth (1993), Karma (1993), Panfilov and Hogeweg (1993), Fenton and Karma (1998), Fenton et al. (2002), Fenton and Cherry (2008), and Zipes et al. (2017). However, in their FHN type, Equations (2), (4), and (5) are quite adapted to simulate the homoclinic orbits underlying AP cycles and cardiac pulse trains. The dynamical complexity of the original FHN model has attracted a lot of mathematical and numerical attention (Glass and Mackey, 1987; Izhikevich, 2007). Under periodic stimulation, period doubling bifurcations were shown to precede a transition to deterministic chaos (Nolasco and Dahlen, 1968; Glass and Mackey, 1979, 1987; Guevara et al., 1981; Chialvo et al., 1990; Fenton and Cherry, 2008). Interestingly, along the line of the diffusely coupled Fitzhugh-Nagumo (Equations 2–5), the generation of higher-dimensional hyperchaotic spatio-temporal dynamics was suggested as a possible explanation of the dynamic transition to fibrillatory states in cardiac tissue (Baier and Müller, 2004).

2.2. Inter-cellular Ion Conduction and GJC Dynamics

The cardiac cell-cell contacts, where electrical signal conduction occurs, are found at intercalated discs located mostly at the narrow end of elongated cardiomyocytes. On their lateral side, the cardiomyocytes are ensheathed by cell-matrix contacts (Figure 1) with weaker electrical coupling. This organization favors a synchronized unidirectional propagation of electrical signals through serial strands of cardiomyocytes. When this synchronization is compromised, life-threatening arrhythmias can develop, and ultimately can become an obstacle for the regeneration of damaged heart tissues. The cardiac intercalated discs (IDs) must therefore be resilient to both mechanical and electrical disturbances to ensure a fast and reproducible propagation of the electrical signal that initiates contraction throughout the heart. ID includes three main structures: (i) the desmosome and (ii) the adherens junction (AJ) that provide the mechanical strength and continuity of the cell-cell contact, and that are both connected to the cytoskeleton, and (iii) the gap junction (GJ) which couples the cells electrically and metabolically. Other proteins which are not directly involved in the cell-cell contact also reside in the ID, including ion channels. The close contact and communication between cardiac myocytes is therefore essential for proper heart functioning as a syncytium, for both electrical and mechanical signal conductions. The GJCs (type I) which are found in cardiac muscles are organized in hexagonal arrays (connexins), with 8–9 center-to-center spacing, and overall thickness ~ 15-16 nm. Despite their negligible size, as compared to the myocyte length, these junctions were soon recognized as discontinuous conduction zones of the myocardium (Spach, 1983), and they were suspected to offer a low resistance to the transfer of electrical signals and more surprisingly to enhance this transfer nonlinearly (Cole et al., 1988).

Several numerical and theoretical attempts have tried to reconcile the difference of scales (from membrane and GJCs thickness to cell length) in a continuous or quasi-continuous treatment of cardiac fiber excitability (Diaz et al., 1983; Joyner et al., 1984; Cole et al., 1988; Keener, 1990; Bub et al., 2005; Pumir et al., 2005; Hand and Griffith, 2010; Hand and Peskin, 2010; Lin and Keener, 2010, 2013). At the heart tissue level, it appears at first legitimate to discard the GJC discontinuities, except maybe for radical instances of static or permanent alterations of GJCs due e.g., to connexin expression depletion (Cole et al., 1988) or gap junction blocker like α-glycerrhetinic acid (Bub et al., 2005) over large areas. As the longitudinal electric field generated by repeatedly passing APs drives the GJCs, we propose in this study to consider the spatio-temporal dynamics of the GJC conductances especially at cardiac intercalated discs where most GJCs reside. Indeed, we do not explicitly account for the discreteness of the GJCs spatially, considering instead an average behavior, but we introduce a nonlinear coupling term in the conduction current mono-domain cable (Equation 2) to mark the distinct propagation and polarization dynamics of GJCs and channel-gated ion channels.

Transjunctional voltage gated inactivation and recovery rates are comparable to repolarization rates (Neyton and Trautmann, 1986; Wang et al., 1992; Desplantez et al., 2007). This is a clue that we regard as pivotal for the behavior of ion transport along the direction of propagation at the locations of GJCs, especially during rapidly beating electrical activity (AP). As compared to previous AP propagation models involving GJCs (Burt and Sray, 1988; Kumar and Gilula, 1996; Harris, 2001; Bukauskas and Verselis, 2004; Danik et al., 2004; Perrachia, 2004), we put more attention on the dynamical implication of GJC (slow) kinetics in response to positive charge accumulation (capacitive charging). By considering that the capacitive charging of closed GJCs is proportional to the voltage drop across the interstitial space, we show that the GJC slow kinetics can enhance the effect thereby amplifying fluctuations around a stationary conductance. The whole AP propagation can therefore be strongly affected because small initial fluctuations of the GJC conductance may be amplified to non negligible values corresponding to a sort of nonlinear resonance of the GJC dynamics.

Averaging over typical depolarization time and length scales, corresponding to several thousands of intercalated discs, allows to derive a continuous mathematical model of inter-cellular ion condition and GJC dynamics. When combined with the well controlled FHN mono-domain mode, we get the following 1D system of nonlinearly coupled PDEs:

$\{\begin{array}{ll}{c}_m\frac{\partial }{\partial t}{U}_m=\mu U_m(A-U_m)(U_m-B)\hfill & \hfill \\ -{\gamma }^2{w}_m-c_g^{-1}\frac{\partial }{\partial x}(g{\rho }_g)+\kappa \frac{{\partial }^2}{\partial x^2}{U}_m,\hfill & \left(6\text{a}\right)\hfill \\ \frac{\partial }{\partial t}{w}_m={\alpha }^2{U}_m-{\nu }_0{w}_m+\lambda ,\hfill & \left(6\text{b}\right)\hfill \\ \frac{\partial }{\partial t}g=+{\omega }^2{\rho }_g-{\nu }_{1}g,\hfill & \left(6\text{c}\right)\hfill \\ \frac{\partial }{\partial t}{\rho }_g=-g\frac{\partial }{\partial x}{U}_m-{\nu }_2{\rho }_g,\hfill & \left(6\text{d}\right)\hfill \end{array}$

with the same boundary conditions as defined in Equation (3). The new dynamical variables are: (i) the GJC conductance fluctuations g (around its mean value accounted for in κ) (in S ≡ Ω⁻¹), which can be interpreted as the non zero gradient part $g\equiv \frac{\partial }{\partial x}\kappa \ne 0$, and (ii) the average capacitive charge linear density ρ_g. The new parameters are c_g: an average linear capacitance analog to c_m but for the intercalated disks, ω²: a reaction rate for the zeroth-order kinetic Equation (6c), and ν₁ a relaxation rate (ν₂ → 0). Note that the zeroth-order (independent of g) reaction rate is justified because only a small fraction of abnormally closing GJCs have to grow from zero initially, while a majority of others remain conducting. At the same time, capacitive charges act on all nearby hexamers. In addition, current conservation modifies Equation (2) so as to include a current term coming from the capacitive charging of GJCs. Precisely, the local ohmic current flowing through open GJCs experiences losses due to closed GJCs. Thus, adding the GJC current as $I_g=c_g^{-1}g{\rho }_g$, to the unpertubed conduction current $I_{\Omega }=-\kappa \frac{\partial }{\partial x}{U}_m$, gives a total longitudinal current I = I_Ω + I_g (Figure 1B). From charge conservation hypothesis, taking the divergence of I yields Equation (6a). Equation (6d) is obtained by considering that the net capacitive charging current density $\frac{\partial }{\partial t}{\rho }_g$ is non zero as soon as the GJCs start closing massively. This current density is written as the product of the GJC conductance and the local electric field produced by the AP wave ($-\frac{\partial }{\partial x}{U}_m$). Finally, ν₁ accounts for the GJC relaxation rate while ν₂ (with ν₂ ≪ ν₁) accounts for other charge leakages. We refer the reader to Table 1 for more details on units and numerical values.

TABLE 1

Table 1. Parameters used in the simulations of Equations (6)–(8).

The boundary conditions are defined such that the value of A is chosen opposite (A → −A) for x = 0, making this fiber end self-sustained oscillating, and unconstrained for x = L with null electric field $-\frac{\partial }{\partial x}{U}_m(x=L)=0$, while no external current is added I_ext = 0. Thus, we set :

$\begin{array}{l}{c}_m\frac{\partial }{\partial t}{U}_m=\mu U_m(-A-U_m)(U_m-B)-{\gamma }^2{w}_m\\ -c_g^{-1}g{\rho }_g\mathrm{\text{ for }}x=0,& (7)\end{array}$

and

$\begin{array}{ll}\frac{\partial }{\partial t}{\rho }_g=-{\nu }_2{\rho }_g\mathrm{\text{ for }}x=L.& (8)\end{array}$

2.3. Numerical Scheme

We use standard finite difference techniques to integrate numerically the system of partial differential Equation (6) with boundary conditions defined in Equations (7) and (8). The linear Laplacian operator is handled with the standard “Crank-Nicholson” scheme using tri-diagonalization. But, one notable unusual trait, as compared to other cardiac models, is that by construction here, the spatial operator is no longer elliptic because of the divergence term. This calls for a special treatment of the gap junction fluctuating quantities.

One way to find out a stable scheme is by considering the following heuristic. For the sake of exposition simplicity, the mean membrane potential is assumed adiabatically constant $\frac{\partial }{\partial t}{U}_m=0$. Equation (6d) possesses an important mirror antisymmetry upon changing U_m → −U_m and x → −x, which distinguishes flow directions of GJ charging or discharging, for a given fluctuating conductance. Consistently, Equations (6b, c) are invariant under the change ρ_g → −ρ_g and g → −g. Hence, the fluctuating vector $W=\left[\begin{array}{l}{\rho }_g\hfill \\ g\hfill \end{array}\right]$ obeys the following temporal evolution equation:

$\begin{array}{ll}\frac{\partial }{\partial t}W=\mathrm{\text{M}}W,& (9)\end{array}$

where the matrix local evolution operator is

$\begin{array}{ll}\mathrm{\text{M}}=\left[\begin{array}{cc}-{\nu }_2& -\frac{\partial }{\partial x}{U}_m\\ {\omega }^2& -{\nu }_1\end{array}\right],& (10)\end{array}$

which shows (ν₂ → 0) that the quantity $\theta =\omega \sqrt{\left|\frac{\partial }{\partial x}{U}_m\right|}$ can act, for some fluctuation, as a typical growth rate (as in the case of a right moving depolarizing front for which the gradient is negative) or as a typical pulsation frequency otherwise. To handle this particular linear, (temporal) evolution, the so-called “leap-frog” method is usually well adapted and was found numerically stable.

When fluctuations can be assumed as happening over a local spatial distance δℓ such that a positive variation from zero is $W\simeq -\delta \ell \frac{\partial }{\partial x}W$, one can rewrite dimensionally the temporal evolution operator (Equation 10) into a spatial transfer operator: $\mathrm{\text{M}}W\to -{\mathrm{\text{M}}}^{\prime }\frac{\partial }{\partial x}W-\nu W$ where

$\begin{array}{ll}{\mathrm{\text{M}}}^{\prime }=\left[\begin{array}{cc}0& \delta \ell \frac{\partial }{\partial x}{U}_m\\ -\delta \ell {\omega }^2& 0\end{array}\right],& (11)\end{array}$

and $v=\left[\begin{array}{l}{v}_2\hfill \\ v_1\hfill \end{array}\right]$. Thus, an advection velocity appears of magnitude v = δℓθ, directed in one or the other way depending on the sign of the fluctuating variables. If not further driven, these fluctuating variables are damped over a spatial distance (along a characteristic) Δx = ±vΔt, in a time $\Delta t~{\nu }_1^{-1}$. A standard “upwind” scheme takes good care of the numerical stability of advection, with sufficient precision, but favors a particular direction of propagation for positive fluctuations, here in agreement with the mirror antisymmetry. This also explains why we choose to place on the outer left side of our 1D system a rapidly beating cell (Equation 7) mimicking an ectopic source, whereas the dynamics is monitored only on its right, at different spatial positions. The outer right boundary is chosen non-excitable in the present study, which mimics any kind of connection to non-conducting tissues. Numerical integration values of time step δt = 0.4 ms for all simulations in Table 1 (except Simul #5 for which δt = 0.25 ms), and spatial step δx = 0.3 mm are chosen to abide by the Courant-Friedrichs-Lewy condition vδt/δx < 1, for which an upper bound estimate of the conduction velocity is c ≃ 0.1 m/s.

2.4. Software and Documentation

The numerical integration code can be downloaded directly from the following link: https://geostat.bordeaux.inria.fr/images/fwd-matlab-code.zip or through Geostat team software page: https://geostat.bordeaux.inria.fr/index.php/downloads.html (last item).

2.5. From the Formal Model Variables to Experimentally Recorded Potential Values

We simulated the external electric potential recorded by a bipolar electrode Δ_b Φ by use of the dipole layer approximation. Within the mono-domain framework, the electric potential felt just outside a cardiac fiber of thickness e ∽ 2 mm is evaluated as (Plonsey and Barr, 2007; Macfarlane et al., 2011):

$\begin{array}{l}{\Delta }_b\Phi (x,t)={\Pi }_p\star i_d(x,t)={\displaystyle \int }{\Pi }_p(x-x^{\prime };e)\\ \left(\kappa \frac{{\partial }^2}{\partial x^2}{U}_m(x^{\prime },t)-i_g(x^{\prime },t)\right)\mathrm{\text{d}}{x}^{\prime },& (12)\end{array}$

where Π_p is the convolution function for the bipolar probe, $i_d(x,t)=\kappa \frac{{\partial }^2}{\partial x^2}{U}_m(x,t)-i_g(x,t)$ the total dipole layer current divergence (or flux across a fiber section), and $i_g(x,t)=\frac{\partial }{\partial x}{I}_g(x,t)$ the new contribution coming from GJC losses. Since we are interested in variations over distances greater than the depolarization length scale, a good approximation is to compute ${\Pi }_p(x)=\frac{\partial }{\partial x}\Pi \left(x,e\right)$, where $\Pi (x,e)=2\pi \left(\frac{x}{\left|x\right|}-\frac{x}{\sqrt{x^2+e^2}}\right)\approx \int |r{|}^{-3}\overrightarrow{r}·\overrightarrow{\mathrm{\text{d}}S}$ is the solid angle spanning the fiber section at a distance $\overrightarrow{r}=(x,e)$ of the depolarization front from the probe. Figure 2 shows typical pseudo bipolar potential time series numerically simulated with our 1D system of PDEs (Equation 6) with, as boundary condition at x = 0 (Equation 7), an automatically beating source of frequency αγ ~ 5 Hz so as to match the cardiac pulse trains observed experimentally. Tuning the newly introduced parameters ω² and ν₁ in Equation (6), we have found quite easily paths leading from a phase of coherent propagation of AP pulses to a phase of quite incoherent and intermittent electrical activity (Figure 2) that strongly reminds the very irregular behavior of electric potential time series recorded during AF (see for comparison Figure 1 in our companion paper I Attuel et al., 2017). Besides the obvious interest of analyzing the succession of bifurcations and transition events encountered along these paths in parameter space, we will focus in this paper on a comparative study of the complex and highly intermittent modulation of cardiac pulse trains simulated numerically with our model of cardiac AP conduction and GJC dynamics and the one observed experimentally in the coronary sinus during episodes of AF (Attuel et al., 2017).

FIGURE 2

Figure 2. Numerically simulated pseudo bipolar potential. (A) 6 s time-series Δ_b Φ (Equation 12) are generated at positions x (in δx = 0.3 mm units) = 8, 18, 75 and 134 (from top to bottom) with our model defined in Equations (6)–(8) with parameter values defined in Simul #2 (Table 1), L = 150 (in δx = 0.3 mm units). (B) Corresponding Fourier transform power spectra. No attempt was made to reproduce the high and low pass filtering used in real bipolar catheter acquisitions. The natural point source frequency at the boundary is ∽ α γ ≈ 5 Hz. Rarefaction of pulses occurs as one moves away from the source which is a hint at some randomly back-scattered pulses that collide and annihilate up-following pulses without reexciting new pulses.

3. Methods of Analysis

3.1. Local Impulse Energy

With the same convention as in the companion paper I (Attuel et al., 2017), the local 1D impulse electric energy traveling with scarcely any alteration at velocity c over a depolarization time period τ_d, is evaluated from Equation (12) as:

$\begin{array}{ll}E(x,t)={\left(\frac{\partial \Delta {\varphi }_b(x,t)}{\partial t}\right)}^2,& (13)\end{array}$

after dropping some constant prefactors. To practically derive E(x, t) from the numerically simulated Δ_b ϕ (x, t), we used the same order 4 finite difference scheme as in the companion paper I (Attuel et al., 2017).

3.2. Zooming on the Local Impulse Energy With the Wavelet Transform Microscope

In Figure 3 is shown in a comparative time-scale decomposition of two 100 s long local (x fixed) impulse energy time-series, the first one was experimentally recorded at the electrode Pt3 in a patient with chronic AF (companion paper I, Attuel et al., 2017) and the second ones was numerically simulated with Equations (6)–(8) with model parameters defined in Simul #2 (Table 1). The continuous wavelet transform (WT) consists in expanding signals in terms of wavelet constructed from a single function, the “analyzing wavelet” ψ, by means of translations and dilations (Grossmann and Morlet, 1984; Daubechies, 1992; Meyer, 1992; Mallat, 1998):

$\begin{array}{ll}{T}_{\psi }\left[E\right](t_0,a)=\frac{1}{a}{\displaystyle {\int }_{-\infty }^{+\infty }}E(t)\psi (\frac{t-t_0}{a})dt,& (14)\end{array}$

where t₀ is a time parameter and a (> 0) a scale parameter (inverse of frequency). Interestingly, by choosing a wavelet ψ which has its first n_ψ moments null $\left[\int t^m\psi (t)dt=0,0\le m<n_{\psi }\right]$, it can be proven (Arneodo et al., 1988, 1995b, 2008; Jaffard, 1989; Muzy et al., 1991, 1994; Mallat and Hwang, 1992) that:

$\begin{array}{ll}{T}_{\psi }\left[E\right](t_0,a)\backsim a^{h(t_0)},a\to 0^{+},& (15)\end{array}$

provided n_ψ > h(t₀), where h(t₀) is the point-wise Hölder exponent that characterizes the maximum regularity of the signal E at point t₀. As experienced in the companion paper I (Attuel et al., 2017) for experimental local impulse energy time-series, we will use in this work the third derivative of a Gaussian function as analyzing wavelet with n_ψ = 3 (Muzy et al., 1994; Arneodo et al., 1995b) (Figure S1 in the companion paper I):

$\begin{array}{ll}{g}^{(3)}(t)=\frac{d^3}{d{t}^3}\left(e^{-t^2/2}\right).& (16)\end{array}$

Interestingly, to track cusp singularities (for oscillating singularities like chirps see Arneodo et al., 1995a, 1997), the WT skeleton defined by the so-called maxima lines of local WT modulus maxima (WTMM) (Figures 3C,C'), was proved to be of practical use since along these maxima lines, Equation (15) was shown to apply (Mallat and Hwang, 1992).

FIGURE 3

Figure 3. Wavelet transform of local impulse energy time-series. (A) A 100 s portion of E(t) recorded at the electrode Pt3 (Companion paper I Attuel et al., 2017). (B) Time-scale WT representation of E(t) with the third-order analyzing wavelet g⁽³⁾ (Equation 16). The modulus of the WT is coded, independently at each scale a, using 256 colors from black ($|T_{g^{(3)}}(t,a)|=0$) to red (${max}_t|T_{g^{(3)}}(t,a)|$). (C) WT skeleton defined by the maxima lines of a 10 s portion of E(t). The scale a = αΔt/δt, where α is an analyzing wavelet dependent constant (α = 8.6 10⁻² for g⁽³⁾ with the lastwave software), and δt = 0.4 ms. (A′–C′) same as (A–C) for a numerical time series E(x, t) generated at position x = 75 (in δx = 0.3 mm units) with our model defined in Equations (6)–(8) with the set of parameter values used in Simul #2 (Table 1), and a total system length L = 150 (in δx = 0.3 mm units). In (B′) the white horizontal dashed-dotted lines delimit the range of time scales (2^8.5 ≤ a ≤ 2^13.5) used to perform linear regression fit estimates of the τ(q) and D(h) multifractal spectra.

3.3. A Wavelet-Based Multifractal Formalism

The wavelet transform modulus maxima (WTMM) method (Muzy et al., 1991, 1994; Bacry et al., 1993; Arneodo et al., 1995b, 2008) was originally developed to generalize box-counting techniques (Arneodo et al., 1987) and to remedy the limitations of the structure functions method (Muzy et al., 1993) in performing multifractal analysis of one-dimensional (1D) velocity signals in fully-developed turbulence. It has proved very efficient to estimate scaling exponents and multifractal spectra (Muzy et al., 1994; Audit et al., 2002; Arneodo et al., 2008; Ivanov et al., 2009). The WTMM method has been extensively applied in different domains of fundamental and applied science, including the analysis of complex time-series found in genomics (Nicolay et al., 2007; Arneodo et al., 2011; Audit et al., 2013) and physiological systems (Ivanov et al., 1999, 2001, 2009; Nunes Amaral et al., 2001; Goldberger et al., 2002; Ciuciu et al., 2012; Chudácek et al., 2014; Gerasimova et al., 2014; Richard et al., 2015; Gerasimova-Chechkina et al., 2016). The WTMM method and its wavelet leaders discrete generalization (Jaffard et al., 2007; Wendt et al., 2007) have already been applied to cardiac signals but mainly to characterize the fluctuations of inter-beat intervals (Ivanov et al., 1999, 2001; West et al., 2004; Leonarduzzi et al., 2010; Wendt et al., 2014; Gadhoumi et al., 2018). As for the analysis of experimental local impulse energy time-series in the companion paper I (Attuel et al., 2017), we will use here two declinations of the WTMM method, the moment (partition function) method and the magnitude cumulant method (Muzy et al., 1994; Arneodo et al., 2008).

The method of moments consists in investigating the scaling behavior of partition functions defined in term of WTMM:

$\begin{array}{ll}Z(q,a)={\displaystyle \sum _{l\in {L}(a)}}|T_{\psi }\left[E\right](t,a){|}^q\backsim a^{\tau \left(q\right)},a\to 0^{+},& (17)\end{array}$

where q ∈ ℝ, ${L}(a)$ is the set of the maxima lines l that defines the WT skeleton and the exponents q and τ(q) play, respectively, the role of an inverse temperature and a free energy in the analogy that links the multifractal formalism and thermodynamics (Bohr and Tél, 1988; Arneodo et al., 1995b). Then, from the Legendre transform of τ(q):

$\begin{array}{ll}D(h)={\displaystyle \underset{q}{min}}\left[qh-\tau (q)\right],& (18)\end{array}$

we get as equivalent of entropy, the so-called D(h) singularity spectrum defined as the fractal (Hausdorff) dimension of the set of points t where the Hölder exponent (equivalent of energy) h(t) = h (Bacry et al., 1993; Muzy et al., 1993, 1994; Arneodo et al., 1995b; Jaffard, 1997a,b).

In practice, to avoid instabilities in performing the Legendre transform, we instead compute the following expectation values (Muzy et al., 1994; Arneodo et al., 1995b), analogous to the fundamental thermodynamic relations, by inversion of Equation (18):

$\begin{array}{ll}h(q,a)=\frac{\partial }{\partial q}ln(Z(q,a))={\displaystyle \sum _{l\in {L}(a)}}ln(|T_{\psi }\left[E\right](l,a)|)·W_{\psi }\left[E\right](q,l,a),& (19)\end{array}$

and

$\begin{array}{l}D(q,a)=q\frac{\partial }{\partial q}ln(Z(q,a))-ln(Z(q,a))\\ ={\displaystyle \sum _{l\in {L}(a)}}{W}_{\psi }\left[E\right](q,l,a)·ln(W_{\psi }\left[E\right](q,l,a)),& (20)\end{array}$

where $W_{\psi }\left[E\right](q,l,a)=|T_{\psi }\left[E\right](l,a){|}^q/Z(q,a)$ corresponds to the Bolzmann weight (Arneodo et al., 1995b). Then, from the slopes of h(q, a) and D(q, a) vs. ln a, we get h(q) and D(q), and therefore the D(h) singularity spectrum as a curve parametrized by q.

An alternative method that can be used as a double check of the predictions of the method of moments is the so-called method of magnitude cumulants (Delour et al., 2001). This method consists in computing the cumulants C_n(a) of the magnitude ln |T_a|. Then from the behavior of the cumulants:

$\begin{array}{ll}\begin{array}{c}{C}_1(a)\equiv \langle ln\left|T_a\right|\rangle \backsim c_{1}ln(a),\hfill \\ C_2(a)\equiv \langle {ln}^2\left|T_a\right|\rangle -{\langle ln\left|T_a\right|\rangle }^2\backsim -c_{2}lna,\hfill \\ C_3(a)\equiv \langle {ln}^3\left|T_a\right|\rangle -3\langle {ln}^2\left|T_a\right|\rangle +2{\langle ln\left|T_a\right|\rangle }^3\backsim c_{3}lna,\hfill \\ ···\hfill \end{array}\hfill & (21)\hfill \end{array}$

we get the following expansion formula for τ(q):

$\begin{array}{ll}\tau (q)-c_0+c_{1}q-c_2{q}^2/2!+c_3{q}^3/3!···& (22)\end{array}$

where the coefficients c_n > 0 are estimated as the slope of C_n(a) vs. ln a (n = 1, 2, 3,, ··· ), and c₀ = D_F as the fractal dimension of the support of singularities of E(t).

Multifractal analysis allows us to distinguish monofractal signals of unique Hölder exponent H. Their τ(q) spectrum is a linear function of q with slope c₁ = H, all the other c_i = 0, i ≥ 2. The corresponding D(h) singularity spectrum reduces to a single point D(h = c₁) = D_F. In contrast, a nonlinear τ(q) is the signature of multifractal signals with Hölder exponent h(t) fluctuating over time (Muzy et al., 1991, 1994; Bacry et al., 1993; Arneodo et al., 1995b, 2002, 2008). As in the companion paper I (Attuel et al., 2017), in this study we will fit the numerical data by so-called log-normal quadratic approximation

$\begin{array}{ll}\tau (q)=-c_0+c_{1}q-c_2{q}^2/2,& (23)\end{array}$

leading to a quadratic single hump shaped D(h) singularity spectrum

$\begin{array}{ll}D(h)=c_0-{(h-c_1)}^2/2c_2,& (24)\end{array}$

where c₀ = −τ(0) = D_f is the fractal dimension of the support of singularities of E(t), c₁ is the value of h that maximizes D(h), and the intermittency coefficient c₂ characterizes the width of the D(h) spectrum (Delour et al., 2001).

3.4. Beyond One-Point Statistics: The Two-Point Magnitude Correlation Method

Many studies have misleadingly extrapolated a multifractal diagnosis to the existence of an underlying multiplicative cascade process. To address this issue, we indeed need to investigate two-point statistics. The two-point magnitude correlation method amounts to investigate how the two-point magnitude correlation function (Arneodo et al., 1998a)

$\begin{array}{ll}\begin{array}{l}C(a,\Delta t)=\langle \left(\ln |T_a(t)|-\langle \ln |T_a(t)|\rangle \right)·\hfill \\ \left(\ln |T_a(t+\Delta t)|-\langle \ln |T_a(t)|\rangle \right)\rangle \hfill \end{array}& (25)\end{array}$

changes as a function of Δt at scale a. As demonstrated by Arneodo et al. (1998a,b) for random multiplicative cascades on wavelet dyadic trees (see also Meneveau and Sreenivasan, 1991):

$\begin{array}{ll}C(a,\Delta t)\backsim -c_{2}ln\Delta t,\Delta t>a,& (26)\end{array}$

where the proportionality coefficient c₂ is the intermittency coefficient defined in Equations (21) and (22) [Note that C(a, Δt = 0) ≡ C₂(a) ∽ − c₂ln a]. Thus, by computing C(a, Δt) from Equation (25) and plotting it as a function of ln Δt, inferences can be made about long-range dependence and consistency with a multiplicative cascading process (Arneodo et al., 1998a,b). Applications of the two-point magnitude correlation method have already provided insight into a wide variety of problems, e.g., the validation of the log-normal cascade phenomenology of fully developed turbulence (Arneodo et al., 1998a,c, 1999) and of high resolution temporal rainfall (Venugopal et al., 2006; Roux et al., 2009), and the demonstration of the existence of a causal cascade of information from large to small scales in financial time series (Arneodo et al., 1998d; Muzy et al., 2001). In the companion paper I (Attuel et al., 2017), we have applied this method to experimental local impulse energy time-series during episodes of AF. Surprisingly, this study has revealed the absence of an underlying multiplicative time-scale structure that will be used in this work as a multifractal random noise numerical test of the pertinence of our cardiac excitable network modeling (section 2) during AF.

4. Results

4.1. One-Point Multi-Fractal Analysis of Local Impulse Energy Numerical Data

4.1.1. Multifractal Analysis of the Simulated Impulse Energy Time-Series With the WTMM Method of Moments

When applying the WTMM method to the impulse energy time-series E(x, t) obtained by numerically integrating Equation 6 with periodically driven boundary condition (Equation 7) at x = 0 (Figure 3A'), we confirmed that the partition function Z(q, a) (Equation 17) obtained from the WT computed with the analysing wavelet g⁽³⁾ (Figure 3B') and its skeleton (Figure 3C'), displays scaling properties for q = −1 up to ≲ 5 over a range of time-scales larger than the typical interbeat ~ 0.2 s (Figure 2B). We strictly limited this range to (1.7, 54 s) for linear regression fit estimates in a logarithmic representation (Figure 4A). The τ(q) spectrum so-obtained at different spatial positions x = 8, 18, 75 and 134 for a 1D system of total length L = 150 (x and L expressed in δx = 0.3 mm units) is robustly well approximated by a quadratic spectrum (Equation 23) with parameter values that do not change much when moving away from the periodic beating source at x = 0 (Table 2). The support of singularities of E(x, t) has a fractal dimension c₀ = D_F ~ 1, independently of x. The singularities of Hölder h = c₁ ~ 0.45−0.48 that maximizes D(h) are also consistently observed all along our 1D spatial system. Interestingly, the intermittency coefficient c₂ is significantly different from zero (the hallmark of multifractal signals) and increases from values c₂ ~ 0.03 close to the source up to values c₂ ~ 0.1 far away from the source. Overall the values of c₀, c₁ and c₂ parameters reported in Table 2 are very similar to the ones obtained in Table 1 of the companion paper I (Attuel et al., 2017) for electrodes Pt3 and Pt5 located in the left atrial posterior wall. Our numerical simulations thus reveals some increased intermittency of E(x, t), when moving away from the periodic beating source, rather rapidly converging to a τ(q) spectrum in quantitative agreement with the one observed experimentally (Figure 5A). This is confirmed when, respectively, plotting h(q, a)/ln 2 (Equation 19) and D(q, a)/ln 2 (Equation 20) vs. log₂a (Figures 4B,C). Despite some deterioration of scaling for large q ≳ 4 values, from the estimate of the slopes h(q) and D(q), we get the single humped D(h) spectra shown in Figure 5B and which clearly widen when moving away from the source to ultimately match the D(h) spectra observed experimentally (Figure 5B). As a check of the reliability of the so-obtained D(h) spectra, they are found quite well approximated by the quadratic spectra defined in Equation 24 with the parameter values obtained from a polynomial fitting of the τ(q) data (Table 2, Figure 5).

FIGURE 4

Figure 4. Multifractal analysis of local impulse energy time-series E(x, t) generated with our model defined in Equations (6)–(8) with the parameter values defined in Simul #2 (Table 1), L = 150 (in δx = 0.3 mm units). (A) log₂Z(q, a) vs. log₂a (Equation 17). (B) h(q, a)/ln 2 vs log₂a (Equation 19). (C) D(q, a)/ln 2 vs. log₂a (Equation 20). The computation were performed with the WTMM method (Paper I, Methods of Analysis Attuel et al., 2017) for different values from q = −1 to 5 with the analyzing wavelet g⁽³⁾ (Equation 16). The vertical dashed lines delimit the range of scales (2^8.5 ≤ a ≤ 2^13.5) used for the linear regression estimate of τ(q), h(q) and D(q) in Figure 5. The symbols correspond to the time-series E(x, t) computed at the spatial positions x = 8 (◦), 18 (□), 75 (▽) and 134 (△) (in δx = 0.3 mm units).

TABLE 2

Table 2. Results of the WTMM multifractal analysis of local impulse energy time-series E(x, t) numerically simulated with our model defined in Equations (6)–(8) with the set of parameter values defined in Simul #2 (Table 1), L = 150 (in δx = 0.3 mm units).

FIGURE 5

Figure 5. Multifractal spectra of local impulse energy time-series E(x, t) generated with our model defined in Equations (6)–(8) with the set of parameter values defined in Simul #2 (Table 1), L = 150 (in δx = 0.3 mm units). (A) τ(q) vs. q estimated by linear regression fit of log₂Z(q, a) vs. log₂a (Figure 4A). (B) D(h) vs. h obtained from linear regression fits of h(q, a) (Figure 4B) and D(q, a) (Figure 4C) vs. log₂a. The symbols have the same meaning as in Figure 4. The curves correspond to quadratic spectra Equations (23) and (24) with parameters [c₀, c₁, c₂] reported in Table 2 for time-series E(x, t) computed at the spatial positions x = 8 (◦), 18 (□), 75 (▽) and 134 (△) (in δx = 0.3 mm units). For comparison are reported the spectra previously obtained for the experimental time-series recorded at the electrodes Pt3 (blue ▼) and Pt5 (green ▼) in the left atrial posterior wall (Companion paper I Attuel et al., 2017).

4.1.2. Multifractal Analysis of the Impulse Energy Data With the Method of Magnitude Cumulants

As previously performed in our companion paper I (Attuel et al., 2017) for the analysis of experimental impulse energy data, we have reproduced our multifractal analysis of numerical time-series using the alternate magnitude cumulant methodology. The first-, second- and third-order cumulants where computed using Equation (21) and are plotted vs. the logarithm of the scale in Figure 6. As expected, C₁(a), C₂(a) and C₃(a) display consistent scaling behavior over the same range of scales (2^8.5 ≤ a ≤ 2^13.5). The results obtained for C₃(a) (Figure 6C) confirm that with limited 422 s long time series as recorded in experiments, there is no way to conclude about the possible departure from a log-normal quadratic τ(q) spectrum (c₃ ≡ 0). Interestingly, the quadratic τ(q) spectra with parameter values $c_1^{*}$ and $c_2^{*}$ listed in Table 2, are found in good agreement with the ones previously estimated with the method of moments. This not only strengthens the multifractal diagnosis of the local impulse energy at low frequencies but it further confirms that farther from the source, larger the intermittency coefficient c₂, and closer to the experimental multifractal spectra obtained in the companion paper I (Attuel et al., 2017).

FIGURE 6

Figure 6. Magnitude cumulant analysis of local impulse energy time-series E(x, t) generated with our model defined in Equations (6)–(8) with the set of parameter values defined in Simul #2 (Table 1), L = 150 (in δx = 0.3 mm units). (A) C₁(a)/ln2 vs. log₂a. (B) C₂(a)/ln2 vs. log₂a. (C) C₃(a)/ln2 vs. log₂a. The computation of the C_n(a) (Equation 21) was performed with the third-order analyzing wavelet g⁽³⁾ (Equation 16). The vertical dashed lines delimit the range of scales (2^8.5 ≤ a ≤ 2^13.5) used for the linear regression estimate of coefficients $c_1^{*}$, $c_2^{*}$ and $c_3^{*}$ of τ(q) (Equation 22) reported in Table 1. The symbols have the same meaning as in Figures 4, 5.

4.2. Robustness of the Multifractal Properties of Local Impulse Energy Under Model Parameter Changes

4.2.1. 1D Spatial System Length L

As a first test of the robustness of the computed multifractal properties of local impulse energy time series, we have performed additional simulations of our 1D PDE system (Equations 6–8) for the same parameter values as before but changing the total length (in δx = 0.3 mm units) of our cellular array L = 210 (Simul #1), 150 (Simul #2), 90 (Simul #3) and 30 (Simul #4). As shown in Figure 7 and Table 3, when using both the methods of moments and of magnitude cumulants to compute the τ(q) and D(h) spectra of E(x = L/2, t) at the midpoint of the array, we recover qualitatively similar multifractal spectra with c₀ = D_F = 1, c₁ ~ 0.45–0.50 and an intermittency coefficient c₂ that increases when increasing L. This is nothing but a confirmation that E(x, t) becomes more and more intermittent when moving the spatial position X = L/2 away from the periodically beating source at x = 0. For 1D arrays as small as L = 30, c₂ ~ 0.015 corresponding to a very narrow D(h) spectrum (Figure 7B). When increasing L, this D(h) spectrum widens progressively to become comparable to the ones obtained experimentally (c₂ ~ 0.1) in Table 1 of the companion paper I (Attuel et al., 2017) for the electrodes Pt3 and Pt5 located in the left atrial posterior wall.

FIGURE 7

Figure 7. Multifractal spectra of local impulse energy time-series E(x, t) generated with our model defined in Equations (6)–(8) with the sets of parameter values defined in Table 1. (A) τ(q) vs. q. (B) D(h) vs. h. These spectra were computed at the same relative spatial position x = L/2 for different lengths L = 30 (Simul #4, ◦, ··· ), 90 (Simul #3, □, - - - - - - 150(Simul#2, ▽, _____), and 210 (Simul #1, △, ) (in δx = 0.3 mm units). The curves correspond to quadratic spectra (Equations 23 and 24) with parameters [c₀, c₁, c₂] reported in Table 3. For comparison are reported the spectra previously obtained from the experimental time-series recorded at the electrodes Pt3 (blue ▼) and Pt5 (green ▼) in the left atrial posterior wall (Companion paper I Attuel et al., 2017).

TABLE 3

Table 3. Results of the WTMM multifractal analysis of local impulse energy time-series E(x = L/2, t) numerically simulated with our model defined in Equations (6)–(8) with the set of parameter values defined in Table 1 for different lengths: L = 30 (Simul #4), L = 90 (Simul #3), L = 150 (Simul #2), and L = 210 (Simul #1) (in δx = 0.3 mm units).

4.2.2. Fiber Conductivity κ

When changing the fiber conductivity parameter κ that accounts for the diffusive coupling of the cells in our 1D cell array (Equation 6a), not much modification of the τ(q) and D(h) spectra is observed (Figure 8). Some weak narrowing of D(h) is obtained when increasing κ corresponding to a small but systematic decrease of the intermittency coefficient from c₂ ~ 0.07 for κ = 0.01 mm.Ω⁻¹ to c₂ ~ 0.04 for κ = 0.08 mm.Ω⁻¹ (Table 4). This is an indication that when strengthening the inter-cell conduction coupling, keeping all the other model parameters fixed, one somehow reduces the multifractal (intermittent) desynchronization of our 1D excitable cell network.

TABLE 4

Table 4. Results of the WTMM multifractal analysis of local impulse energy time-series E(x = L/2, t) numerically simulated with our model defined in Equations (6)–(8) for different values of the conductivity parameter κ defined in Simul #2, Simul #6 and Simul #7 (Table 1), L = 150 (in δx = 0.3 mm units).

FIGURE 8

Figure 8. Multifractal spectra of local impulse energy time-series E(x, t) generated with our model defined in Equations (6)–(8) with the sets of parameter values defined in Simul #6 (Table 1) and L = 150 (in δx = 0.3 mm units). (A) τ(q) vs. q. (B) D(h) vs. h. The spectra were computed at the same relative spatial position x = L/2 = 75 for different conductivities of the fiber: Simul #2 (▽, _____), Simul #6 (+, - - - - -) and Simul #7 (*, ··· ). The curves correspond to quadratic spectra (Equations 23 and 24) with parameters [c₀, c₁, c₂] reported in Table 4.

4.2.3. New Set of Parameters

We have also reproduced the one-point multifractal analysis reported in section 4.1.1 (Simul #2) to local impulse energy numerical time-series generated with the parameter set defined in Simul #5 (Table 1) and the same 1D system spatial size L = 150 (in δx = 0.3 mm units). Many parameters were changed, including κ, μ, ν₀, ν₁, γ², α² and ω² (Table 1). As shown in Figure 9, the τ(q) and D(h) spectra obtained at different spatial positions x = 8, 18, 75, and 134 are quite similar to the ones previously obtained in Figure 5. They are robustly approximated by quadratic spectra (Equations 23 and 24, respectively) with comparable c₀, c₁, and c₂ parameter values (Table 5).

TABLE 5

Table 5. Results of the WTMM multifractal analysis of local impulse energy time-series E(x, t) numerically simulated with our model defined in Equations (6)–(8) with the parameters defined in Simul #6 (Table 1), L = 150 (in δx = 0.3 mm units).

FIGURE 9

Figure 9. Multifractal spectra of local impulse energy time-series E(x, t) generated with our model defined in Equations (6)–(8) with the set of parameter values defined in Simul #5 (Table 1) and L = 150 (in δx = 0.3 mm units). (A) τ(q) vs. q. (B) D(h) vs. h. The spectra were computed with the WTMM method with the third-order analyzing wavelet g⁽³⁾ (Equation 16). The symbols correspond to the spatial positions x = 8 (•, ··· ), 18 (■, - - - - -), 75 (▼, ___________), and 139 (▲, ) (in δx = 0.3 mm units). The curves correspond to quadratic spectra (Equations 23 and 24) with parameters [c₀, c₁, c₂] reported in Table 5. For comparison are reported in open symbols (◦, □, ▽ △), the corresponding spectra previously obtained with the set of parameter values defined in Simul #2 (Table 1) and L = 150 in Figure 5.

4.3. Two-Point Magnitude Analysis of Local Impulse Energy Data

The results of the two-point magnitude correlation analysis of the local impulse energy time series numerically generated with our 1D PDE system (Equations 6–8) with the set of parameter values defined in Simul #2 are shown in Figure 10. C(a, Δt)/C(a, 0) (Equation 25) computed with the third-order analyzing wavelet g⁽³⁾ (Equation 16) is represented vs. ln (Δt) for two scales a = 2⁹ and 2¹⁰ in the scaling range. Strikingly, for the four numerical time-series corresponding to different spatial positions x = 8 (Figure 10A), 18 (Figure 10B), 75 (Figure 10C), and 134 (Figure 10D), for time-lag Δt ≳ a, C(a, Δ)/C(a, 0) drops to zero as a clear indication that the magnitudes are uncorrelated. As a reference, we put in each panel in Figure 10, a dashed straight line of slope −c₂ as predicted by Equation (26) for multifractal signals exhibiting a cascading multiplicative structure along a time-scale tree (Arneodo et al., 1998a,b). The slow decay predicted by the “multiplicative” log-normal model with intermittency coefficient c₂ is definitely not observed. Thus, the numerical local impulse energy time series look much more like what has been called log-normal multifractal random noise in pioneering works to distinguish “uncorrelated” and “multiplicative” log-normal models (Arneodo et al., 1998a). Importantly, a similar absence of magnitude correlation was observed with the experimental time-series recorded at the electrodes located in the left atrial posterior wall in our companion paper I (Attuel et al., 2017). This is a strong evidence that our cardiac excitable cell network model indeed accounts for both one-point and two-point statistics of local impulse energy time-series during AF.

FIGURE 10

Figure 10. Two-point magnitude analysis of local impulse energy time-series E(x, t) generated with our model defined in Equations (6)–(8) with the sets of parameter values defined in Simul #2 (Table 1), L = 150 (in δx = 0.3 mm units). Two-point correlation function C(a, Δt)/C(a, 0) vs. ln (Δt) (Equation 25) computed with the third-order analyzing wavelet g⁽³⁾ (Equation 16). The two curves correspond to scales a = 2⁹ (black) and 2¹⁰ (gray) within the scaling range. The different panels correspond to different spatial positions x = 8 (A), 18 (B), 75 (C), and 134 (D) (in δx = 0.3 mm units).

5. Conclusions

To summarize, we proposed a model of cardiac excitable cell network which accounts for the transport of AP along and across myocardial cells via the spatio-temporal interplay of voltage-gated and gap junction channels kinetics. We demonstrated that this model robustly reproduces the multifractal intermittent nature of the cardiac impulse energy experimentally recorded in the left atrial posterior wall area over times (≳ 0.5 s) longer that the mean interbeat (≃ 10⁻¹ s) during AF (companion paper I, Attuel et al., 2017). In particular, this model gives full account of the experimental observation of the absence of a multiplicative time-scale structure underlying multifractal scaling. To our knowledge, our combined experimental and numerical studies are the first to report the observation, quantification and modeling of such multifractal dynamics which is found more complex than previously suspected. Preliminary exploration of the model bifurcation diagram suggests that it shares with other models a good reproducibility of the spectrum of rhythmic and AP disorders, such as early after depolarizations (EADs) or salvos of premature beats, found experimentally prior to the onset of AF. This stems here from the membrane current imbalance between depolarization and repolarization, originating in the capacitive currents building up at the GJs. In the model, the nonlinearity of the GJC temporal response was proposed to be due to a nonlinear coupling of the local electric field with the GJC charging during AP propagation. However, the nature of our studies was exploratory, with a data set limited to a few patients, and although it was performed on time-series rather long for clinical practice (422 s), they were not so long regarding the range of time scales [0.6, 10 s] where scaling was observed. The relevance of our modeling would definitely benefit of the analysis of new data over a large set of patients at different stages of AF development and to be explored in different areas of the atria. Recording electric potential time series during AF concomitantly to non intrusive Cardiovascular Magnetic Resonance (CMR) imaging could help further the assessment of atrial remodeling features such as increased expression of intercellular gap junction and conduction velocity shortening, in addition to sinus node dysfunction. Even more instructive, the comparative analysis of different types of rhythms as atrial flutter, AV junctional rhythm and various other annotated rhythms (Gadhoumi et al., 2018) would allow us to evaluate the practicality of multifractal cardiac impulse energy in the discrimination of AF from other rhythms. Also, by combining our wavelet-based multifractal analysis (low frequency) to a more classical dynamical system analysis including bifurcation diagrams, Lyapunov exponent computations, we should be in position to improve and refine our physiological heart tissue modeling and to open new perspectives toward the understanding of mechanisms of AF perpetuation.

Author Contributions

GA, FA, HY, and AA: conception and design; GA, FA, and AA: development and methodology; GA, EG-C, FA, and AA: analysis and interpretation of data; GA, EG-C, FA, HY, and AA: writing, review, and/or revision of the manuscript; GA, EG-C, and FA: administrative, technical or material support (i.e., requiring and organizing data, constructing databases).

Funding

This work was partially supported by Contrat Conseil Région Aquitaine CAVERNOM Cardiac Arrhythmia Complexity and Variability by means of Robust Nonlinear Methods (Grant N^o 2014-1R60212-00003295), by the President of Russian for the young scientists (Grant N^o 14.W01.17.2674-MK), and by the Metchnikov Program (EG-C visit to LOMA).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The data used in this study have been provided by IHU Liryc, Hôpital Xavier Arnozan, Avenue du Haut Lévêque, 33604 Pessac cedex.

Footnotes

1. ^${\sigma }^{-1}={\sigma }_i^{-1}+{\sigma }_e^{-1}$, where σ_i and σ_e are the internal and external medium conductivities in bi-domain models.

References

Allessie, M. A., Bonke, F. I., and Schopman, F. J. (1977). Circus movement in rabbit atrial muscle as a mechanism of tachycardia. Circ. Res. 41, 9–18. doi: 10.1161/01.RES.41.1.9

PubMed Abstract | CrossRef Full Text | Google Scholar

Armstrong, C. M. (2003). The Na/K pump, Cl ion, and osmotic stabilization of cells. Proc. Natl. Acad. Sci. U.S.A. 100, 6257–6262. doi: 10.1073/pnas.0931278100

PubMed Abstract | CrossRef Full Text | Google Scholar

Arneodo, A., Audit, B., Decoster, N., Muzy, J.-F., and Vaillant, C. (2002). “A wavelet based multifractal formalism: application to DNA sequences, satellite images of the cloud structure and stock market data,” in The Science of Disasters: Climate Disruptions, Heart Attacks, and Market Crashes, eds A. Bunde, J. Kropp, and H. J. Schellnhuber (Berlin: Springer Verlag), 26–102.

Google Scholar

Arneodo, A., Audit, B., Kestener, P., and Roux, S. G. (2008). Wavelet-based multifractal analysis. Scholarpedia 3:4103. doi: 10.4249/scholarpedia.4103

CrossRef Full Text | Google Scholar

Arneodo, A., Bacry, E., Jaffard, S., and Muzy, J.-F. (1997). Oscillating singularities on Cantor sets: a grand-canonical multifractal formalism. J. Stat. Phys. 87, 179–209. doi: 10.1007/BF02181485

CrossRef Full Text | Google Scholar

Arneodo, A., Bacry, E., Manneville, S., and Muzy, J.-F. (1998a). Analysis of random cascades using space-scale correlation functions. Phys. Rev. Lett. 80, 708–711. doi: 10.1103/PhysRevLett.80.708

CrossRef Full Text | Google Scholar

Arneodo, A., Bacry, E., and Muzy, J.-F. (1995a). Oscillating singularities in locally self-similar functions. Phys. Rev. Lett. 74, 4823–4826. doi: 10.1103/PhysRevLett.74.4823

PubMed Abstract | CrossRef Full Text | Google Scholar

Arneodo, A., Bacry, E., and Muzy, J.-F. (1995b). The thermodynamics of fractals revisited with wavelets. Physica A 213, 232–275. doi: 10.1016/0378-4371(94)00163-N

CrossRef Full Text | Google Scholar

Arneodo, A., Bacry, E., and Muzy, J.-F. (1998b). Random cascades on wavelet dyadic trees. J. Math. Phys. 39, 4142–4164. doi: 10.1063/1.532489

CrossRef Full Text | Google Scholar

Arneodo, A., Grasseau, G., and Holschneider, M. (1988). Wavelet transform of multifractals. Phys. Rev. Lett. 61, 2281–2284. doi: 10.1103/PhysRevLett.61.2281

PubMed Abstract | CrossRef Full Text | Google Scholar

Arneodo, A., Grasseau, G., and Kostelich, E. J. (1987). Fractal dimensions and f(α) spectrum of the Hénon attractor. Phys. Lett. A 124, 426–432. doi: 10.1016/0375-9601(87)90546-9

CrossRef Full Text | Google Scholar

Arneodo, A., Manneville, S., and Muzy, J.-F. (1998c). Towards log-normal statistics in high Reynolds number turbulence. Eur. Phys. J. B 1, 129–140. doi: 10.1007/s100510050162

CrossRef Full Text | Google Scholar

Arneodo, A., Manneville, S., Muzy, J.-F., and Roux, S. G. (1999). Revealing a lognormal cascading process in turbulent velocity statistics with wavelet analysis. Philos. Trans. R. S. London Ser. A 357, 2415–2438. doi: 10.1098/rsta.1999.0440

CrossRef Full Text | Google Scholar

Arneodo, A., Muzy, J.-F., and Sornette, D. (1998d). “Direct” causal cascade in the stock market. Eur. Phys. J. B 2, 277–282. doi: 10.1007/s100510050250

CrossRef Full Text | Google Scholar

Arneodo, A., Vaillant, C., Audit, B., Argoul, F., d'Aubenton-Carafa, Y., and Thermes, C. (2011). Multi-scale coding of genomic information: from DNA sequence to genome structure and function. Phys. Rep. 498, 45–188. doi: 10.1016/j.physrep.2010.10.001

CrossRef Full Text | Google Scholar

Attuel, G., Gerasimova-Chechkina, E., Argoul, F., Yahia, H., and Arneodo, A. (2017). Multifractal desynchronization of the cardiac excitable cell network during atrial fibrillation. I. Multifractal analysis of clinical data. Front. Physiol. 8:01139. doi: 10.3389/fphys.2017.01139

PubMed Abstract | CrossRef Full Text | Google Scholar

Attuel, G., Pont, O., Xu, B., and Yahia, H. (2016). “Sudden cardiac death and turbulence,” in The Foundations of Chaos Revisited: From Poincaré to Recent Advancements, ed C. Skiadas (Switzerland: Springer), 235–248.

Google Scholar

Attuel, P., Pellerin, D., and Gaston, J. (1989). “Latent atrial vulnerability: new means of electrophysiologic investigations in paroxysmal atrial arrhythmias,” in The Atrium in Health and Disease, eds P. Attuel, P. Coumel, and M. Janse (Mont Kisco, NY: Mount KiscoFutura Publishing Co., Inc), 159–200.

Google Scholar

Audit, B., Bacry, E., Muzy, J.-F., and Arneodo, A. (2002). Wavelet-based estimators of scaling behavior. IEEE Trans. Info. Theory 48, 2938–2954. doi: 10.1109/TIT.2002.802631

CrossRef Full Text | Google Scholar

Audit, B., Baker, A., Chen, C.-L., Rappailles, A., Guilbaud, G., Julienne, H., et al. (2013). Multiscale analysis of genome-wide replication timing profiles using a wavelet-based signal-processing algorithm. Nat. Protoc. 8, 98–110. doi: 10.1038/nprot.2012.145

PubMed Abstract | CrossRef Full Text | Google Scholar

Ausma, J., van der Velden, H. M. W., Lenders, M.-H., van Ankeren, E. P., Jongsma, H. J., Ramaekers, F. C. S., et al. (2003). Reverse structural and gap-junctional remodeling after prolonged atrial fibrillation in the goat. Circulation 107, 2051–2058. doi: 10.1161/01.CIR.0000062689.04037.3F

PubMed Abstract | CrossRef Full Text | Google Scholar

Bacry, E., Muzy, J.-F., and Arneodo, A. (1993). Singularity spectrum of fractal signals from wavelet analysis: exact results. J. Stat. Phys. 70, 635–674. doi: 10.1007/BF01053588

CrossRef Full Text | Google Scholar

Baier, G., and Müller, M. (2004). Excitable chaos in diffusively coupled FitzHugh-Nagumo equations. Rev. Mex. Fis. 50, 422–426.

Google Scholar

Bär, M., and Brusch, L. (2004). Breakup of spiral waves caused by radial dynamics: eckhaus and finite wavenumber instabilities. New J. Phys. 6:5. doi: 10.1088/1367-2630/6/1/005

CrossRef Full Text | Google Scholar

Bär, M., and Eiswirth, M. (1993). Turbulence due to spiral breakup in a continuous excitable medium. Phys. Rev. E 48:R1635. doi: 10.1103/PhysRevE.48.R1635

PubMed Abstract | CrossRef Full Text | Google Scholar

Beeler, G. W., and Reuter, H. (1977). Reconstruction of the action potential of ventricular myocardial fibres. J. Physiol. 268, 177–210. doi: 10.1113/jphysiol.1977.sp011853

PubMed Abstract | CrossRef Full Text | Google Scholar

Beuter, A., Glass, L., Mackey, M. C., and Titcombe, M. S. (2003). Nonlinear Dynamics in Physiology and Medicine. New York, NY: Springer Science+Business Media.

Google Scholar

Bikou, O., Thomas, D., Trappe, K., Lugenbiel, P., Kelemen, K., Koch, M., et al. (2011). Connexin 43 gene therapy prevents persistent atrial fibrillation in a porcine model. Cardiovasc. Res. 92, 218–225. doi: 10.1093/cvr/cvr209

PubMed Abstract | CrossRef Full Text | Google Scholar

Bohr, T., and Tél, T. (1988). “The thermodynamics of fractals,” in Directions in Chaos, Vol. 2, ed B.-L. Hao (Singapore: World Scientific Publishing Co), 194–237.

Google Scholar

Bub, G., Shrier, A., and Glass, L. (2005). Global organization of dynamics in oscillatory heterogeneous excitable media. Phys. Rev. Lett. 94:028105. doi: 10.1103/PhysRevLett.94.028105

PubMed Abstract | CrossRef Full Text | Google Scholar

Bukauskas, F. F., and Verselis, V. K. (2004). Gap junction channel gating. Biochim. Biophys. Acta 1662, 42–60. doi: 10.1016/j.bbamem.2004.01.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Burt, J. M., and Sray, D. C. (1988). Single-channel events and gating behavior of the cardiac gap junction channel. Proc. Natl. Acad. Sci. U.S.A. 85, 3431–3434. doi: 10.1073/pnas.85.10.3431

PubMed Abstract | CrossRef Full Text | Google Scholar

Chialvo, D. R., Gilmour, R. J., and Jalife, J. (1990). Low dimensional chaos in cardiac tissue. Nature 343, 653–657. doi: 10.1038/343653a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Chudácek, V., Andén, J., Mallat, S., Abry, P., and Doret, M. (2014). Scattering transform for intrapartum fetal heart rate variability fractal analysis: a case-control study. IEEE Trans. Biomed. Eng. 61, 1100–1108. doi: 10.1109/TBME.2013.2294324

PubMed Abstract | CrossRef Full Text | Google Scholar

Ciuciu, P., Varoquaux, G., Abry, P., Sadaghiani, S., and Kleinschmidt, A. (2012). Scale-free and multifractal time dynamics of fMRI signals during rest and task. Front. Physiol. 3:186. doi: 10.3389/fphys.2012.00186

PubMed Abstract | CrossRef Full Text | Google Scholar

Cole, W. C., Picone, J. B., and Sperelakis, N. (1988). Gap junction uncoupling and discontinuous propagation in the heart. Biophys. J. 53, 809–818. doi: 10.1016/S0006-3495(88)83160-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Courtemanche, M., Glass, L., and Keener, J. P. (1993). Instabilities of a propagating pulse in a ring of excitable media. Phys. Rev. Lett. 70, 2182–2185. doi: 10.1103/PhysRevLett.70.2182

PubMed Abstract | CrossRef Full Text | Google Scholar

Cox, J. L., Schuessler, R. B., D'Agostino, H. J. J., Stone, C. M., Chang, B. C., Cain, M. E., et al. (1991). The surgical treatment of atrial fibrillation. III. Development of a definitive surgical procedure. J. Thorac. Cardiovasc. Surg. 101, 569–583.

PubMed Abstract | Google Scholar

Danik, S. B., Liu, F., Zhang, J., Suk, H. J., Morley, G. E., Fishman, G. I., et al. (2004). Modulation of cardiac gap junction expression and arrhythmic susceptibility. Circ. Res. 95, 1035–1041. doi: 10.1161/01.RES.0000148664.33695.2a

PubMed Abstract | CrossRef Full Text | Google Scholar

Daubechies, I. (1992). Ten Lectures on Wavelets. CBMS-NSF Reg. Conf. Ser. Appl. Math. Vol 61. Philadelphia, PA: SIAM.

Google Scholar

de la Casa, M. A., de la Rubia, F. J., and Ivanov, P. C. (2007a). Patterns of phase-dependent spiral wave attenuation in excitable media. Phys. Rev. E 75:051923. doi: 10.1103/PhysRevE.75.051923

PubMed Abstract | CrossRef Full Text | Google Scholar

de la Casa, M. A., de la Rubia, F. J., and Ivanov, P. C. (2007b). Patterns of spiral wave attenuation by low-frequency periodic planar fronts. Chaos 17:015109. doi: 10.1063/1.2404640

PubMed Abstract | CrossRef Full Text | Google Scholar

Delour, J., Muzy, J.-F., and Arneodo, A. (2001). Intermittency of 1D velocity spatial profiles in turbulence: a magnitude cumulant analysis. Eur. Phys. J. B 23, 243–248. doi: 10.1007/s100510170074

CrossRef Full Text | Google Scholar

Desplantez, T., Dupont, E., Severs, N. J., and Weingart, R. (2007). Gap junction channels and cardiac impulse propagation. J. Membr. Biol. 218, 13–28. doi: 10.1007/s00232-007-9046-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Diaz, P. J., Rudy, Y., and Plonsey, R. (1983). Intercalated discs as a cause for discontinuous propagation in cardiac muscle: a theoretical simulation. Ann. Biomed. Eng. 11, 177–189. doi: 10.1007/BF02363285

PubMed Abstract | CrossRef Full Text | Google Scholar

Dupont, E., Matsushita, T., Kaba, R. A., Vozzi, C., Coppen, S. R., Khan, N., et al. (2001). Altered connexin expression in human congestive heart failure. J. Mol. Cell. Cardiol. 33, 353–371. doi: 10.1006/jmcc.2000.1308

PubMed Abstract | CrossRef Full Text | Google Scholar

Evans, W. H., and Martin, P. E. M. (2002). Gap junctions: structure and function (review). Mol. Membr. Biol. 19, 121–136. doi: 10.1080/09687680210139839

PubMed Abstract | CrossRef Full Text | Google Scholar

Fenton, F. H., and Cherry, E. M. (2008). Models of cardiac cell. Scholarpedia 3:1868. doi: 10.4249/scholarpedia.1868

CrossRef Full Text | Google Scholar

Fenton, F. H., Cherry, E. M., Hasting, H. M., and Evans, S. J. (2002). Multiple mechanisms of spiral wave breakup in a model of cardiac electrical activity. Chaos 12, 852–892. doi: 10.1063/1.1504242

PubMed Abstract | CrossRef Full Text | Google Scholar

Fenton, F. H., and Karma, A. (1998). Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: filament instability and fibrillation. Chaos 8, 20–47. doi: 10.1063/1.166311

CrossRef Full Text | Google Scholar

FitzHugh, R. (1962). Mathematical Models of Excitation and Propagation in Nerve. New York, NY: McGraw-Hill, Biological Engineering.

Google Scholar

Gadhoumi, K., Do, D., Badilini, F., Pelter, M. M., and Hu, X. (2018). Wavelet leader multifractal analysis of heart rate variability in atrial fibrillation. J. Electrocardiol. 51, S83–S87 doi: 10.1016/j.jelectrocard.2018.08.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Ganesan, A. N., Shipp, N. J., Brooks, A. G., Kuklik, P., Lau, D. H., Lim, H. S., et al. (2013). Long-term outcomes of catheter ablation of atrial fibrillation: a systematic review and meta-analysis. J. Am. Heart Assoc. 2:e004549. doi: 10.1161/JAHA.112.004549

PubMed Abstract | CrossRef Full Text | Google Scholar

Garfinkel, A., Chen, P.-S., Walter, D. O., Karagueuzian, H. S., Kogan, B., Evans, S. J., et al. (1997). Quasiperiodicity and chaos in cardiac fibrillation. J. Clin. Invest. 99, 305–314. doi: 10.1172/JCI119159

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerasimova, E., Audit, B., Roux, S. G., Khalil, A., Gileva, O., Argoul, F., et al. (2014). Wavelet-based multifractal analysis of dynamic infrared thermograms to assist in early breast cancer diagnosis. Front. Physiol. 5:176. doi: 10.3389/fphys.2014.00176

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerasimova-Chechkina, E., Toner, B., Marin, Z., Audit, B., Roux, S. G., Argoul, F., et al. (2016). Comparative multifractal analysis of dynamic infrared thermograms and X-ray mammograms enlightens changes in the environment of malignant tumors. Front. Physiol. 7:336. doi: 10.3389/fphys.2016.00336

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerhardt, M., Schuster, M., and Tyson, J. J. (1990). A cellular automaton model of excitable media including curvature and dispersion. Science 247, 1563–1566. doi: 10.1126/science.2321017

CrossRef Full Text | Google Scholar

Glass, L., and Mackey, M. C. (1979). A simple model for phase locking of biological oscillations. J. Math. Biol. 6, 207–223. doi: 10.1007/BF02547797

CrossRef Full Text | Google Scholar

Glass, L., and Mackey, M. C. (1987). From Clocks to Chaos: The Rhythms of Life. Princeton, NJ: Princeton University Press.

Google Scholar

Goldberger, A. L., Amaral, L. A. N., Hausdorff, J. M., Ivanov, P. C., Peng, C.-K., and Stanley, H. E. (2002). Fractal dynamics in physiology: alterations with disease and aging. Proc. Natl. Acad. Sci. U.S.A. 99, 2466–2472. doi: 10.1073/pnas.012579499

PubMed Abstract | CrossRef Full Text | Google Scholar

Gollob, M. H., Jones, D. L., Krahn, A. D., Danis, L., Gong, X.-Q., Shao, Q., et al. (2006). Somatic mutations in the connexin 40 gene (gja5) in atrial fibrillation. N. Eng. J. Med. 354, 2677–2688. doi: 10.1056/NEJMoa052800

PubMed Abstract | CrossRef Full Text | Google Scholar

Gray, R. A., Jalife, J., Panfilov, A. V., Baxter, W. T., Cabo, C., Davidenko, J. M., et al. (1995). Mechanisms of cardiac fibrillation. Science 270, 1222–1225. doi: 10.1126/science.270.5239.1222

PubMed Abstract | CrossRef Full Text | Google Scholar

Grossmann, A., and Morlet, J. (1984). Decomposition of Hardy functions into square integrable wavelets of constant shape. SIAM J. Math. Anal. 15, 723–736. doi: 10.1137/0515056

CrossRef Full Text | Google Scholar

Guckenheimer, J., and Kuehn, C. (2009). Homoclinic orbits of the Fitzhugh Nagumo equation: the singular-limit. Discrete Contin. Dyn. Syst. 2, 851–872. doi: 10.3934/dcdss.2009.2.851

CrossRef Full Text | Google Scholar

Guevara, M. R., and Glass, L. (1982). Phase locking, period-doubling bifurcations, and chaos in a mathematical model of a periodically driven oscillator: a theory for the entrainement of biological oscillators and the generation of cardiac dysrhythmias. J. Math. Biol. 14, 1–23. doi: 10.1007/BF02154750

CrossRef Full Text | Google Scholar

Guevara, M. R., Glass, L., and Shrier, A. (1981). Phase locking, period-doubling bifurcations, and irregular dynamics in periodically stimulated cardiac cells. Science 214, 1350–1353. doi: 10.1126/science.7313693

PubMed Abstract | CrossRef Full Text | Google Scholar

Haïssaguerre, M., Jaïs, P., Shah, D. C., Takahashi, A., Hocini, M., Quiniou, G., et al. (1998). Spontaneous initiation of atrial fibrillation by ectopic beats originating in the pulmonary veins. N. Engl. J. Med. 339, 659–666. doi: 10.1056/NEJM199809033391003

PubMed Abstract | CrossRef Full Text | Google Scholar

Hand, P. E., and Griffith, B. E. (2010). Adaptive multiscale model for simulating cardiac conduction. Proc. Natl. Acad. Sci. U.S.A. 107, 14603–14608. doi: 10.1073/pnas.1008443107

PubMed Abstract | CrossRef Full Text | Google Scholar

Hand, P. E., and Peskin, C. S. (2010). Homogenization of an electrophysiological model for a strand of cardiac myocytes with gap-junctional and electric-field coupling. Bull. Math. Biol. 72, 1408–1424. doi: 10.1007/s11538-009-9499-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Harris, A. L. (2001). Emerging issues of connexin channels: biophysics fills the gap. Q. Rev. Biophys. 34, 325–472. doi: 10.1017/S0033583501003705

PubMed Abstract | CrossRef Full Text | Google Scholar

Hille, B. (2001). Ion Channels of Excitable Membranes, 3rd Edn. Sunderland: Sinauer Associates, Inc.

Google Scholar

Ho, S. Y., Cabrera, J. A., and Sanchez-Quintana, D. (20012). Left atrial anatomy revisited. Circ. Arrhythm. Electrophysiol. 5, 220–228. doi: 10.1161/CIRCEP.111.962720

PubMed Abstract | CrossRef Full Text | Google Scholar

Hodgkin, A. L., and Huxley, A. F. (1952a). Currents carried by sodium and potassium ions through the membrane of the giant axon of Loligo. J. Physiol. 116, 449–472. doi: 10.1113/jphysiol.1952.sp004717

PubMed Abstract | CrossRef Full Text | Google Scholar

Hodgkin, A. L., and Huxley, A. F. (1952b). A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117, 500–544. doi: 10.1113/jphysiol.1952.sp004764

PubMed Abstract | CrossRef Full Text | Google Scholar

Igarashi, T., Finet, J. E., Takeuchi, A., Fujino, Y., Strom, M., Greener, I. D., et al. (2012). Connexin gene transfer preserves conduction velocity and prevents atrial fibrillation. Circulation 125, 216–225. doi: 10.1161/CIRCULATIONAHA.111.053272

PubMed Abstract | CrossRef Full Text | Google Scholar

Ito, H., and Glass, L. (1991). Spiral breakup in a new model of discrete excitable media. Phys. Rev. Lett. 66, 671–674. doi: 10.1103/PhysRevLett.66.671

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanov, P. C., Amaral, L. A. N., Goldberger, A. L., Havlin, S., Rosenblum, M. G., Stanley, H. E., et al. (2001). From 1/f noise to multifractal cascades in heartbeat dynamics. Chaos 11, 641–652. doi: 10.1063/1.1395631

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanov, P. C., Amaral, L. A. N., Goldberger, A. L., Havlin, S., Rosenblum, M. G., Struzik, Z. R., et al. (1999). Multifractality in human heartbeat dynamics. Nature 399, 461–465. doi: 10.1038/20924

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanov, P. C., Amaral, L. A. N., Goldberger, A. L., and Stanley, H. E. (1998). Stochastic feedback and the regulation of biological rhythms. Europhys. Lett. 43, 363–368. doi: 10.1209/epl/i1998-00366-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanov, P. C., Ma, Q. D. Y., Bartsch, R. P., Hausdorff, J. M., Nunes Amaral, L. A., Schulte-Frohlinde, V., et al. (2009). Levels of complexity in scale-invariant neural signals. Phys. Rev. E 79:041920. doi: 10.1103/PhysRevE.79.041920

PubMed Abstract | CrossRef Full Text | Google Scholar

Izhikevich, E. M. (2007). Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. Cambridge: MIT Press.

Izhikevich, E. M., and FitzHugh, R. (2006). FitzHugh Nagumo model. Scholarpedia 1:1349. doi: 10.4249/scholarpedia.1349

CrossRef Full Text | Google Scholar

Jaffard, S. (1989). Hölder exponents at given points and wavelet coefficients. C. R. Acad. Sci. Paris Ser. I 308, 79–81.

Google Scholar

Jaffard, S. (1997a). Multifractal formalism for functions part I: results valid for all functions. SIAM J. Math. Anal. 28, 944–970. doi: 10.1137/S0036141095282991

CrossRef Full Text | Google Scholar

Jaffard, S. (1997b). Multifractal formalism for functions part II: self-similar functions. SIAM J. Math. Anal. 28, 971–998. doi: 10.1137/S0036141095283005

CrossRef Full Text | Google Scholar

Jaffard, S., Lashermes, B., and Abry, P. (2007). “Wavelet leaders in multifractal analysis,” in Wavelet Analysis and Applications, eds T. Qian, M. I. Vai, and Y. Xu (Basel: Birkhäuser Basel), 201–246.

Google Scholar

Jensen, J. H., Christiansen, P. L., and Scott, A. C. (1984). Chaos in the Beeler-Reuter system for the action potential of ventricular myocardial fibres. Physica D 13, 269–277. doi: 10.1016/0167-2789(84)90283-5

CrossRef Full Text | Google Scholar

Joyner, R. W., Veenstra, R., Rawling, D., and Chorro, A. (1984). Propagation through electrically coupled cells. Effects of a resistive barrier. Biophys. J. 45, 1017–1025. doi: 10.1016/S0006-3495(84)84247-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Karma, A. (1993). Spiral breakup in model equations of action potential propagation in cardiac tissue. Phys. Rev. Lett. 71, 1103–1106. doi: 10.1103/PhysRevLett.71.1103

PubMed Abstract | CrossRef Full Text | Google Scholar

Keener, J. P. (1990). The effects of gap junctions on propagation in myocardium: a modified cable theory. Ann. N. Y. Acad. Sci. 591, 257–277. doi: 10.1111/j.1749-6632.1990.tb15094.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, N. M., and Gilula, N. B. (1996). The gap junction communication channel. Cell 84, 381–388. doi: 10.1016/S0092-8674(00)81282-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Leonarduzzi, R. F., Schlotthauer, G., and Torres, M. E. (2010). “Wavelet leader based multifractal analysis of heart rate variability during myocardial ischaemia,” in Annual International Conference of the IEEE Engineering in Medicine and Biology (Buenos Aires), 110–113.

Google Scholar

Lin, J., and Keener, J. P. (2010). Ephaptic coupling in cardiac myocytes. IEEE Trans. Biomed. Circuits Syst. 60, 576–582.

PubMed Abstract | Google Scholar

Lin, J., and Keener, J. P. (2013). Modeling electrical activity of myocardial cells incorporating the effects of ephaptic coupling. Proc. Natl. Acad. Sci. U.S.A. 107, 20935–20940. doi: 10.1073/pnas.1010154107

PubMed Abstract | CrossRef Full Text | Google Scholar

Macfarlane, P. W., van Oosterom, A., Kligfield, P., Pahlm, O., Janse, M., and Camm, J. (eds.) (2011). Comprehensive Electrocardiography. London: Springer.

Mallat, S. (1998). A Wavelet Tour of Signal Processing. New York, NY: Academic Press.

Google Scholar

Mallat, S., and Hwang, W. L. (1992). Singularity detection and processing with wavelets. IEEE Trans. Inf. Theory 38, 617–643. doi: 10.1109/18.119727

CrossRef Full Text | Google Scholar

Meneveau, C., and Sreenivasan, K. (1991). The multifractal nature of turbulent energy dissipation. J. Fluid Mech. 224, 429–484. doi: 10.1017/S0022112091001830

CrossRef Full Text | Google Scholar

Meyer, Y. (1992). Wavelets and Operators. Cambridge: Cambridge University Press.

Moe, G. K., and Abildskov, J. A. (1959). Atrial fibrillation as a self-sustaining arrhythmia independent of focal discharge. Am. Heart J. 58, 59–70. doi: 10.1016/0002-8703(59)90274-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Moe, G. K., Rheinboldt, W. C., and Abildskov, J. A. (1964). A computer model of atrial fibrillation. Am. Heart J. 67, 200–220. doi: 10.1016/0002-8703(64)90371-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Muzy, J.-F., Bacry, E., and Arneodo, A. (1991). Wavelets and multifractal formalism for singular signals: application to turbulence data. Phys. Rev. Lett. 67, 3515–3518. doi: 10.1103/PhysRevLett.67.3515

PubMed Abstract | CrossRef Full Text | Google Scholar

Muzy, J.-F., Bacry, E., and Arneodo, A. (1993). Multifractal formalism for fractal signals: the structure-function approach versus the wavelet-transform modulus-maxima methods. Phys. Rev. E 47, 875–884. doi: 10.1103/PhysRevE.47.875

PubMed Abstract | CrossRef Full Text | Google Scholar

Muzy, J.-F., Bacry, E., and Arneodo, A. (1994). The multifractal formalism revisited with wavelets. Int. J. Bifurc. Chaos 4, 245–302. doi: 10.1142/S0218127494000204

CrossRef Full Text | Google Scholar

Muzy, J.-F., Sornette, D., Delour, J., and Arneodo, A. (2001). Multifractal returns and hierarchical portfolio theory. Quant. Finance 1, 131–148. doi: 10.1080/713665541

CrossRef Full Text | Google Scholar

Nagumo, J., Arimoto, S., and Yoshizawa, S. (1962). An active pulse transmission line simulating nerve axon. Proceedings of the IRE 50, 2061–2070. doi: 10.1109/JRPROC.1962.288235

CrossRef Full Text | Google Scholar

Nattel, S., and Harada, M. (2014). Atrial remodeling and atrial fibrillation: recent advances and translational perspectives. J. Am. Coll. Cardiol. 63, 2335–2345. doi: 10.1016/j.jacc.2014.02.555

PubMed Abstract | CrossRef Full Text | Google Scholar

Neyton, J., and Trautmann, A. (1986). Physiological modulation of gap junction permeability. J. Exp. Biol. 124, 93–114.

PubMed Abstract | Google Scholar

Nicolay, S., Brodie of Brodie, E. B., Touchon, M., Audit, B., d'Aubenton Carafa, Y., Thermes, C., et al. (2007). Bifractality of human DNA strand-asymmetry profiles results from transcription. Phys. Rev. E 75:032902. doi: 10.1103/PhysRevE.75.032902

PubMed Abstract | CrossRef Full Text | Google Scholar

Niebur, E. (2008). Neuronal cable theory. Scholarpedia 3:2674. doi: 10.4249/scholarpedia.2674

CrossRef Full Text | Google Scholar

Noble, D. (1962). A modification of the Hodgkin-Huxley equations applicable to Purkinje fibre action and pacemaker potentials. J. Physiol. 160, 317–352. doi: 10.1113/jphysiol.1962.sp006849

CrossRef Full Text | Google Scholar

Noble, D. (1965). Electrical properties of cardiac muscle attributable to inward-going (anomalous) rectification. J. Cell Comp. Physiol. 66, 127–136. doi: 10.1002/jcp.1030660520

CrossRef Full Text | Google Scholar

Nolasco, J. B., and Dahlen, R. W. (1968). Graphic method for the study of alternation in cardiac action potential. J. Appl. Physiol. 5, 191–196. doi: 10.1152/jappl.1968.25.2.191

CrossRef Full Text | Google Scholar

Nunes Amaral, L. A., Ivanov, P. C., Aoyagi, N., Hidaka, I., Tomono, S., Goldberger, A. L., et al. (2001). Behavioral-independent features of complex heartbeat dynamics. Phys. Rev. Lett. 86, 6026–6029. doi: 10.1103/PhysRevLett.86.6026

PubMed Abstract | CrossRef Full Text | Google Scholar

Panfilov, A., and Hogeweg, P. (1993). Spiral breakup in a modified FitzHugh-Nagumo model. Phys. Lett. A 176, 295–299. doi: 10.1016/0375-9601(93)90921-L

CrossRef Full Text | Google Scholar

Perrachia, C. (2004). Chemical gating of gap junction channels Roles of calcium, pH and calmodulin. Biochim. Biophys. Acta 1662, 61–80. doi: 10.1016/j.bbamem.2003.10.020

CrossRef Full Text | Google Scholar

Pertsov, A. M., Davidenko, J. M., Salomonsz, R., Baxter, W. T., and Jalife, J. (1993). Spiral waves of excitation underlie reentrant activity in isolated cardiac muscle. Circ. Res. 72, 631–650. doi: 10.1161/01.RES.72.3.631

PubMed Abstract | CrossRef Full Text | Google Scholar

Plonsey, R., and Barr, R. C. (2007). Bioelectricity: A Quantitative Approach. 3rd Edn. London: Springer.

Google Scholar

Pumir, A., Arutunyan, A., Krinsky, V., and Sarvazyan, N. (2005). Genesis of ectopic waves: role of coupling, automaticity, and heterogeneity. Biophys. J. 89, 2332–2349. doi: 10.1529/biophysj.105.061820

PubMed Abstract | CrossRef Full Text | Google Scholar

Richard, C. D., Tanenbaum, A., Audit, B., Arneodo, A., Khalil, A., and Frankel, W. N. (2015). Swdreader: a wavelet-based algorithm using spectral phase to characterize spike-wave morphological variation in genetic models of absence epilepsy. J. Neurosci. Methods 242, 127–140. doi: 10.1016/j.jneumeth.2014.12.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Roux, S. G., Venugopal, V., Fienberg, K., Arneodo, A., and Foufoula-Georgiou, E. (2009). Evidence for inherent nonlinearity in temporal rainfall. Adv. Water Resourc. 32, 41–48. doi: 10.1016/j.advwatres.2008.09.007

CrossRef Full Text | Google Scholar

Severs, N. J. (1990). The cardiac gap junction and intercalated disc. Int. J. Cardiol. 26, 137–173. doi: 10.1016/0167-5273(90)90030-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Severs, N. J., Bruce, A. F., Dupont, E., and Rothery, S. (2008). Remodelling of gap junctions and connexin expression in diseased myocardium. Cardiovasc. Res. 80, 9–19. doi: 10.1093/cvr/cvn133

PubMed Abstract | CrossRef Full Text | Google Scholar

Spach, M. S. (1983). The discontinuous nature of electrical propagation in cardiac muscle. Ann. Biomed. Eng. 2, 209–261. doi: 10.1007/BF02363287

CrossRef Full Text | Google Scholar

Spray, D. C., White, R. L., Campos De Carvalho, A., Harris, A. L., and Bennet, M. V. L. (1984). Gating of gap junction channels. Biophys. J. 45, 219–230. doi: 10.1016/S0006-3495(84)84150-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Swietach, P., Youm, J.-B., Saegusa, N., Leem, C.-H., Spitzer, K. W., and Vaughan-Jones, R. D. (2013). Coupled Ca²⁺/H⁺ transport by cytoplasmic buffers regulates local Ca²⁺ and H⁺ ion signaling. Proc. Natl. Acad. Sci. U.S.A., 110, E2064–E2073. doi: 10.1073/pnas.1222433110

PubMed Abstract | CrossRef Full Text | Google Scholar

Takashi, K. (2007). Van der Pol oscillator. Scholarpedia 2:2202. doi: 10.4249/scholarpedia.2202

CrossRef Full Text | Google Scholar

Takigawa, M., Takahashi, A., Kuwahara, T., Okubo, K., Takahashi, Y., Watari, Y., et al. (2014). Long-term follow-up after catheter ablation of paroxysmal atrial fibrillation: the incidence of recurrence and progression of atrial fibrillation. Circ. Arrhythm. Electrophysiol. 7, 267–273. doi: 10.1161/CIRCEP.113.000471

PubMed Abstract | CrossRef Full Text | Google Scholar

Tosteson, D. C., and Hoffman, J. F. (1960). Regulation of cell volume by active cation transport in high and low potassium sheep red cells. J. Gen. Physiol. 44, 169–194. doi: 10.1085/jgp.44.1.169

PubMed Abstract | CrossRef Full Text | Google Scholar

Van der Pol, B., and Van der Mark, J. (1928). The heartbeat considered as a relaxation oscillation and an electrical model of the heart. Philos. Mag. 7, 763–775. doi: 10.1080/14786441108564652

CrossRef Full Text | Google Scholar

van der Velden, H. M. W., Ausma, J., Rook, M. B., Hellemons, A. J. C. G. M., van Veen, T. A. A. B., Allessie, M. A., et al. (2000). Gap junctional remodeling in relation to stabilization of atrial fibrillation in the goat. Cardiovasc. Res. 46, 476–486. doi: 10.1016/S0008-6363(00)00026-2

PubMed Abstract | CrossRef Full Text | Google Scholar

van der Velden, H. M. W., and Jongsma, H. J. (2002). Cardiac gap junctions and connexins: their role in atrial fibrillation and potential as therapeutic targets. Cardiovasc. Res. 54, 270–279. doi: 10.1016/S0008-6363(01)00557-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Venugopal, V., Roux, S. G., Foufoula-Georgiou, E., and Arneodo, A. (2006). Revisiting multifractality of high-resolution temporal rainfall using a wavelet-based formalism. Water Resour. Res. 42:W06D14. doi: 10.1029/2005WR004489

CrossRef Full Text | Google Scholar

Wang, H.-Z., Jian, J. L., Lemanski, F. L., and Veenstra, R. D. (1992). Gating of mammalian cardiac gap junction channels by transjunctional voltage. Biophys. J. 63, 139–151. doi: 10.1016/S0006-3495(92)81573-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Weidmann, S. (1966). The diffusion of radiopotassium across intercalated disks of mammalian cardiac muscle. J. Physiol. 187, 323–342. doi: 10.1113/jphysiol.1966.sp008092

PubMed Abstract | CrossRef Full Text | Google Scholar

Wendt, H., Abry, P., and Jaffard, S. (2007). Bootstrap for empirical multifractal analysis. IEEE Signal Process. Mag. 24, 38–48. doi: 10.1109/MSP.2007.4286563

CrossRef Full Text | Google Scholar

Wendt, H., Kiyono, K., Abry, P., Hayano, J., Watanabe, E., and Yamamoto, Y. (2014). Multiscale wavelet p-leader based heart rate variability analysis for survival probability assessment in CHF patients. Proceedings of the International Conference of the IEEE Engineering in Medicine and Biology Society (Piscataway, NJ), 2809–2812. doi: 10.1109/EMBC.2014.6944207

PubMed Abstract | CrossRef Full Text | Google Scholar

West, B. J., Scafetta, N., Cooke, W. H., and Balocchi, R. (2004). Influence of progressive central hypovolemia on Hölder exponent distributions of cardiac interbeat intervals. Ann. Biomed. Eng. 32, 1077–1087. doi: 10.1114/B:ABME.0000036644.69559.ad

PubMed Abstract | CrossRef Full Text | Google Scholar

Winfree, A. T. (1989). Electrical instability in cardiac muscle: phase singularities and rotors. J. Theor. Biol. 138, 353–405. doi: 10.1016/S0022-5193(89)80200-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Zipes, D. P., Jalife, J., and Stevenson, W. G. (eds.) (2017). Cardiac Electrophysiology: From Cell to Bedside. Philadelphia, PA: Elsevier Saunders.

Google Scholar