Front. Physiol. Frontiers in Physiology Front. Physiol. 1664-042X Frontiers Media S.A. 10.3389/fphys.2017.01139 Physiology Methods Multifractal Desynchronization of the Cardiac Excitable Cell Network During Atrial Fibrillation. I. Multifractal Analysis of Clinical Data Attuel Guillaume 1 Gerasimova-Chechkina Evgeniya 2 Argoul Francoise 3 Yahia Hussein 1 Arneodo Alain 3 * 1Geometry and Statistics in Acquisition Data, Centre de Recherche INRIA, Talence, France 2Laboratory of Physical Foundation of Strength, Institute of Continuous Media Mechanics UB RAS, Perm, Russia 3Laboratoire Ondes et Matières d'Aquitaine, Université de Bordeaux, Centre National de la Recherche Scientifique, UMR 5798, Talence, France

Edited by: Zbigniew R. Struzik, The University of Tokyo, Japan

Reviewed by: Alok Ranjan Nayak, International Institute of Information Technology, India; Alessandro Giuliani, Istituto Superiore di Sanità, Italy

*Correspondence: Alain Arneodo alain.arneodo@u-bordeaux.fr

This article was submitted to Fractal Physiology, a section of the journal Frontiers in Physiology

 26 03 2018 2017 8 1139 05 10 2017 24 12 2017 Copyright © 2018 Attuel, Gerasimova-Chechkina, Argoul, Yahia and Arneodo. 2018 Attuel, Gerasimova-Chechkina, Argoul, Yahia and Arneodo

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Atrial fibrillation (AF) is a cardiac arrhythmia characterized by rapid and irregular atrial electrical activity with a high clinical impact on stroke incidence. Best available therapeutic strategies combine pharmacological and surgical means. But when successful, they do not always prevent long-term relapses. Initial success becomes all the more tricky to achieve as the arrhythmia maintains itself and the pathology evolves into sustained or chronic AF. This raises the open crucial issue of deciphering the mechanisms that govern the onset of AF as well as its perpetuation. In this study, we develop a wavelet-based multi-scale strategy to analyze the electrical activity of human hearts recorded by catheter electrodes, positioned in the coronary sinus (CS), during episodes of AF. We compute the so-called multifractal spectra using two variants of the wavelet transform modulus maxima method, the moment (partition function) method and the magnitude cumulant method. Application of these methods to long time series recorded in a patient with chronic AF provides quantitative evidence of the multifractal intermittent nature of the electric energy of passing cardiac impulses at low frequencies, i.e., for times (≳0.5 s) longer than the mean interbeat (≃ 10−1 s). We also report the results of a two-point magnitude correlation analysis which infers the absence of a multiplicative time-scale structure underlying multifractal scaling. The electric energy dynamics looks like a “multifractal white noise” with quadratic (log-normal) multifractal spectra. These observations challenge concepts of functional reentrant circuits in mechanistic theories of AF, still leaving open the role of the autonomic nervous system (ANS). A transition is indeed observed in the computed multifractal spectra which group according to two distinct areas, consistently with the anatomical substrate binding to the CS, namely the left atrial posterior wall, and the ligament of Marshall which is innervated by the ANS. In a companion paper (II. Modeling), we propose a mathematical model of a denervated heart where the kinetics of gap junction conductance alone induces a desynchronization of the myocardial excitable cells, accounting for the multifractal spectra found experimentally in the left atrial posterior wall area.

 mechanisms of atrial fibrillation heart electrical activity multifractal analysis wavelet transform modulus maxima method two-point magnitude correlation analysis multifractal noise 2014-1R60212-00003295 14.W01.17.2674-MK Secrétariat Général pour les Affaires Régionales, Etat en Région Aquitaine10.13039/501100002733 Russian Science Foundation10.13039/501100006769 1. Introduction

AF is an arrhythmia originating in the rapid and irregular electrical activity of the atria (the heart's two upper chambers) that causes their pump function to fail, increasing up to fivefold the risk of embolic stroke (Wolf et al., 1978, 1991; Attuel et al., 1986; Lip and Lane, 2015). Isolated short episodes of tachy-arrhythmias may be normal, but if they become more frequent and last longer than 1 min then paroxysmal AF is declared. This condition alone requires treatment of the atrial electrophysiological substrate, all the more so as AF often coexists with and predisposes to heart failure, with increased morbidity and mortality (Middlekauff et al., 1991; Stevenson and Stevenson, 1999; Wang et al., 2003). Management of AF by drug therapy aims at controlling either the ventricular rate, for instance by use of negative chronotropes such as beta-adrenergic blockers, or the rhythm, by use of anti-arrhythmic drugs interfering with ionic membrane currents of the excitable cells to prolong the action potential duration (APD), in combination with anticoagulants, but it does not lead to a cure (Roy et al., 2008; Al-Khatib et al., 2014). Since the work of Haïssaguerre et al. (1998), radio-frequency ablation of the pulmonary veins in the left atria has been developed for patients with paroxysmal AF as drug efficacy was found to be poor or even to become pro-arrhythmic (Echt et al., 1991; The CAST II investigators, 1992). This type of intervention seeks to punctually destroy sources of abnormal focal electrical activity susceptible to trigger the arrhythmia. Alternative strategies that have been developed lately include compartmenting the atrium in order to block possible reentrant circuits, or even directly targeting areas of abnormally fractionated activity (Nademanee et al., 2004; Camm et al., 2010). Despite the association of various strategies, clinical efficiency remains disappointing (Ganesan et al., 2013; Takigawa et al., 2014). The rate of AF recurrences after an initial ablation procedure treating paroxysmal AF increases with time (Bertaglia et al., 2010), necessitating multiple redos, and most patients suffering persistent AF are resistant to treatment (Camm et al., 2010; Verma et al., 2015; Wynn et al., 2016).

The prevailing electrophysiological concepts describing tachy-arrhythmias are more than a century old. They involve abnormal automaticity and conduction (Janse, 1997). Initiation and maintenance are thought to arise from a vulnerable substrate prone to the emergence of multiple self-perpetuating reentry circuits, also called “multiple wavelets” (Moe and Abildskov, 1959; Moe et al., 1964). Reentries may be driven structurally, for instance because of locally high fibrous tissue content which badly conducts, or functionally because of high spatial dispersion of decreased refractoriness and APD (Misier et al., 1992). The latter is coined the leading circle concept with the clinically more relevant notion of a critical “wavelength” (in fact the length) of the cardiac impulse (Allessie et al., 1977; Smeets et al., 1986; Rensma et al., 1988; Attuel et al., 1989). The related concept of vulnerability was originally introduced to uncover a physiological substrate evolving from normality to pathology. It was found in vulnerable patients that high rate frequency would invariably lead to functional disorder as cardiac cells would no longer properly adapt their refractoriness (Attuel et al., 1982). Mathematical models have managed to exhibit likewise phenomena, with the generation of breaking spiral waves in various conditions (Ito and Glass, 1991; Karma, 1993). The triggering role of abnormal ectopic activity of the pulmonary veins has been demonstrated on patients with paroxysmal AF resistant to drug therapy (Haïssaguerre et al., 1998), but its origin still remains poorly understood. This region is highly innervated with sympathetic and parasympathetic stimulation from the ANS (Tan et al., 2007; Ulphani et al., 2007; Arora, 2012). In particular, Coumel et al. (Coumel et al., 1978; Coumel, 1994) have revealed the pathophysiological role of the vagal tone on a vulnerable substrate. It is frequently observed that rapid tachycardia of ectopic origin transits to AF. This is known to result from electrical remodeling. As described for the first time by Allessie et al. (Allessie, 1998), remodeling is a transient and reversible process by which the impulse properties such as its refractory period are altered during the course of the arrhythmia, promoting its perpetuation: “AF begets AF” (Wijffels et al., 1995). Under substantial beating rate increase, cells may undergo remodeling to overcome the toxicity of their excessive intercellular calcium loading, by a rapid down regulation (a few minutes) of their L-type calcium membrane current (Yue et al., 1997; Nattel et al., 2008). Moreover, other ionic channel functions are also modified such as the potassium channel function, inducing a change in the conduction properties including the conduction velocity (Iwasaki et al., 2011; Nattel and Harada, 2014). The intercellular coupling at the gap junction level shows also alterations of their connexin expression and dispersion (Severs et al., 2008). For more details, one may consult van Marion et al. (2015) and Zipes et al. (2017).

In this study, we delve into the complexity of voltage signals recovered with bipolar electrodes in the CS during chronic AF. Attempts to assess visually the spatio-temporal complexity of voltage signals using maps of unipolar electrodes revealed various complex patterns of activity not all compatible with reentries (Konings et al., 1994, 1997). We use here two declinations of a wavelet-based multi-scale method, the moment (partition function) method and the magnitude cumulant method (Muzy et al., 1994; Arneodo et al., 2008), as originally introduced in the field of fully developed turbulence (Muzy et al., 1991). This methodology has been extensively applied in different domains of fundamental and applied sciences, including geophysics (Venugopal et al., 2006), econophysics (Muzy et al., 2001), biology (Arneodo et al., 2011) and medicine (Gerasimova et al., 2014; Gerasimova-Chechkina et al., 2016). In the context of cardiac physiology, this methodology was shown to be valuable in assessing congestive heart failure from the monitoring of sinus heart rate variability (Ivanov et al., 1999, 2001; Goldberger et al., 2002). The wavelet-based multifractal analysis of the electric energy of passing cardiac impulses during AF reported in the paper provides unprecedented experimental estimates of the multifractal spectra in different heart areas. The reported results show that the electric energy dynamics looks like a log-normal “multifractal white noise” with no underlying multiplicative time-scale structure. Long time series recordings, at various locations throughout the whole atria of many patients, are mandatory to obtain and analyze with the aim at identifying how many multifractal scalings exist. For instance, what is the effect on multifractal scaling of altered cell network topologies, especially when fibrosis is present? Or how is such scaling correlated to the disease temporal evolution, with or without drug therapy or surgical intervention, when the substrate is remodeled? In a companion paper (Attuel et al., submitted), we start tackling the problem at source, in the simplest (1D) cell network topology, by exploring the possibility that the substrate function is modulated by the kinetics of conduction. A simple reversible mechanism of short term remodeling under rapid pacing is demonstrated, by which ionic overload acts locally (dynamical feedback) on the kinetics of gap junction conductance. The whole process may propagate and pervade the myocardium via electronic currents. No influence of the ANS is included and no structural inhomogeneities are taken into account. Then the complete network of excitable cells becomes desynchronized, with induced dispersion of remodeled refractoriness and APD, and abnormal automaticity. Contrary to existing mathematical models based on circuit reentries, a spatio-temporal multifractal intermittent dynamics emerges similar to the one found in the CS next to the left atrial posterior wall area, opening a new avenue toward the understanding of AF mechanisms of perpetuation.

 2. Methods of analysis

The wavelet transform (WT) is a mathematical microscope (Arneodo et al., 1988, 1995b, 2008; Muzy et al., 1991, 1994) that is well suited for the analysis of complex non-stationary time-series such as those found in genomics (Nicolay et al., 2007; Arneodo et al., 2011; Audit et al., 2013) and physiological systems (Ivanov et al., 1999, 2001; Goldberger et al., 2002; Ciuciu et al., 2012; Chudácek et al., 2014; Gerasimova et al., 2014; Richard et al., 2015). Thanks to its ability to be blind to non-stationary low-frequency (polynomial) trends in the analyzed signal E(t) (or E(x)), it has been early recognized as well adapted to reveal the hierarchy that governs the temporal (or spatial) distribution of singularities of multifractal signals including singular measures and functions (Arneodo et al., 1988; Holschneider, 1988; Jaffard, 1989; Muzy et al., 1991, 1994; Mallat and Hwang, 1992). It is therefore implemented in robust methods capturing the self-similar intricate fractal structures hidden in signals that exhibit a typical “1/f noise” scaling as seen in the Fourier spectral density (Mandelbrot, 1982, 1998; West and Shlesinger, 1988).

2.1. The wavelet transform microscope: a singularity scanner

The WT is a time-scale decomposition method which consists in expanding signals in terms of wavelets constructed from a single function, the “analyzing wavelet” ψ, by means of translations and dilations (Grossmann and Morlet, 1984; Daubechies, 1992; Meyer, 1992; Mallat, 1998). The WT of a real-valued function E is defined as:

Tψ[E](t0,a)=1a∫−∞+∞E(t)ψ(t−t0a)dt,

where t0 is a time parameter and a (> 0) a scale parameter (inverse of frequency). By choosing a wavelet ψ which has its first nψ moments null [∫tmψ(t)dt=0,0≤m<nψ], it can be proven that the behavior of Tψ[E](t0, a) as a function of the scale a, as a → 0+, characterizes the local behavior of E(t) (Arneodo et al., 1988, 1995b, 2008; Jaffard, 1989; Muzy et al., 1991, 1994; Mallat and Hwang, 1992):

Tψ[E](t0,a)~ah(t0), a→0+,

provided nψ > h(t0), where h(t0) is the point-wise Hölder exponent that characterizes the maximum regularity of the signal E at point t0. If n < h(t0) ≤ n + 1, the (n − 1)th derivative of E(t) is regular and its nth derivative is singular at t0. Thus the larger h(t0), the smoother the function, the faster the power-law decrease of Tψ[E] when a → 0+. For h(t0) = 0, E(t) is discontinuous and bounded at t0 and the wavelet transform no longer depends on a. For discontinuous “noise” signals, h(t0) < 0 and Tψ[E](t0, a) increases when a → 0+. For instance, h(t0) = −1 corresponds to a delta distribution at t = t0, while if almost everywhere h(t)=-12, the single exponent H=12 is characteristic of a “white” noise (Muzy et al., 1994). To resolve all the cusp singularities present in a function, the analyzing wavelet must be chosen to have enough vanishing moments to resolve the singularities with Hölder exponent hmax, namely nψ ≥ hmax. Since hmax is not known a priori, the most appropriate way to correctly estimate all singularities is to analyze the given function with analyzing wavelets of increasing order nψ until a robust estimate of the so-called spectrum of singularities is obtained (Bacry et al., 1993; Muzy et al., 1994; Arneodo et al., 1995b) (see section 2.2). In the present study, we use the successive derivatives of a Gaussian function g(N)(t)=dNdtN(e-t2/2) as analyzing wavelets with nψ = N (Muzy et al., 1994; Arneodo et al., 1995b) (Figure S1).

 2.2. A wavelet-based canonical multifractal formalism: the wavelet transform modulus maxima method

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 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). From the deep analogy that links the multifractal formalism to thermodynamics (Bohr and Tél, 1988; Arneodo et al., 1995b), the WTMM method provides a canonical description (Arneodo et al., 1995b) of the distribution of point-wise Hölder exponents for finite a. The WTMM method allows therefore some control and mastering of finite-size effects and statistical convergence issues. As a → 0+, the thermodynamic limit is reached formally guaranteeing an equivalence with micro-canonical approaches such as found in Turiel et al. (2008). To account for the possible presence of oscillating (chirps) (Arneodo et al., 1995a), a grand-canonical multifractal formalism has also been developed for signals involving cusp and oscillating singularities (Arneodo et al., 1997b). The canonical 1D WTMM method has been generalized in 2D for the multifractal analysis of rough surfaces (Arneodo et al., 2000, 2003; Decoster et al., 2000; Roux et al., 2000) and for the analysis of 3D scalar and vector fields (Kestener and Arneodo, 2003, 2004; Arneodo et al., 2008) with successful applications in astrophysics (Khalil et al., 2006; Kestener et al., 2010; McAteer et al., 2010), geophysics (Venugopal et al., 2006; Roux et al., 2009), surface science (Roland et al., 2009), image processing (Mallat, 1998; Arneodo et al., 2003, 2008; Antoine et al., 2008), cellular biology (Khalil et al., 2007; Snow et al., 2008; Goody et al., 2010; Grant et al., 2010; Martinez-Torres et al., 2014, 2015) and medicine (Kestener et al., 2001; Arneodo et al., 2003; Khalil et al., 2009; Batchelder et al., 2014; Gerasimova-Chechkina et al., 2016; Marin et al., 2017). Note that alternative approaches to the WTMM method have been developed using discrete wavelet bases, including the recent use of wavelet leaders (Jaffard et al., 2007; Wendt et al., 2007).

2.2.1. The method of moments

In 1D, the WTMM method (Muzy et al., 1991, 1993, 1994; Bacry et al., 1993; Arneodo et al., 1995b) consists in computing the WT skeleton defined, at each fixed scale a, by the local maxima L(a) of the WT modulus |Tψ[E](t, a)|. These WTMM are disposed on curves connected across scales called maxima lines (see Figure 2C). Along these maxima lines l, Mallat and Hwang (1992) have shown that Equation (2) also applies for the WTMM that behave as |Tψ[E](t,a)|∽ah(t), where h(t) is the Hölder exponent characterizing the singularity of the signal E at time t. The canonical multifractal formalism (Muzy et al., 1994; Arneodo et al., 1995b) characterizes the relative contributions of each Hölder exponent value via the estimate of the singularity spectrum D(h) defined as the fractal (Hausdorff) dimension of the set of points t where h(t) = h. This spectrum can be obtained by investigating the scaling behavior of partition functions defined in terms of WTMM (and which correspond to the moments of the WTMM probability distribution function):

Z(q,a)=∑l∈L(a)|Tψ[E](t,a)|q~aτ(q), a→0+,

where q ∈ ℝ, and L(a) is the set of all maxima lines l that satisfy: l∈L(a), if ∀a′ ≤ a, ∃(t, a′) ∈ l. In the framework of the analogy with thermodynamics (Bohr and Tél, 1988; Arneodo et al., 1995b), q and τ(q) play respectively the role of an inverse temperature and a free energy. The main result of the canonical wavelet-based multifractal formalism is that in place of energy and entropy (i.e., the variables conjugated to q and τ), we have h, the Hölder exponent, and D(h), the singularity spectrum. This means that the singularity spectrum of E(t) is a convex function that can be calculated from the Legendre transform of the partition function scaling exponents τ(q) (Bacry et al., 1993; Muzy et al., 1993, 1994; Arneodo et al., 1995b):

D(h)=minq[qh−τ(q)].

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 (4):

h(q,a)=∂∂qln(Z(q,a))=∑l∈L(a)ln(|Tψ[E](t,a)|)·Wψ[E](q,l,a),

and

D(q,a)=q∂∂qln(Z(q,a))−ln(Z(q,a))           =∑l∈L(a)Wψ[E](q,l,a)·ln(Wψ[E](q,l,a)),

where Wψ[E](q,l,a)=|Tψ[E](t,a)|q/Z(q,a) corresponds to the Bolzmann weight in the analogy that connects the multifractal formalism to thermodynamics (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. For further mathematical developments on the 1D WTMM method, we refer the reader to Bacry et al. (1993) and Jaffard (1997a,b).

 2.2.2. The method of magnitude cumulants

With the previous method of moments, to compute the entire τ(q) curve, we need to perform linear regression fits of ln Z(q, a) vs. ln a (Equation 3) for a wide range of q values and then to proceed to a polynomial fit of the τ(q) data prior to the Legendre transform (Equation 4) to get the D(h) singularity spectrum. An alternative method based on magnitude cumulants has been introduced by Delour et al. (2001) to minimize the number of linear regression fits (as few as 3) while still adequately inferring and accurately estimating the nonlinear behavior of the τ(q) spectrum. This method is based on the following reasoning. The computation of the partition function Z(q, a) amounts to computing the following arithmetic mean of the WTMM to the power q:

〈|Ta|q〉=1NaZ(q,a),

where we simplified notations Ta ≡ Tψ[E](·, a), and where Na is the number of maxima lines at scale a, which scales as ∽a-Df, where Df = −τ(0) is the fractal dimension of the support set of the singularities in the signal E(t). From Equations (3) and (7), we get the expansion

[τ(q)+Df]lna~ln{〈eqln|Ta|〉},                        ~∑n=1∞Cn(a)qnn!,

where Cn(a) are the cumulants of the magnitude ln |Ta|. Then from the behavior of these cumulants:

C1(a)≡〈ln|Ta|〉~c1ln(a),C2(a)≡〈ln2|Ta|〉−〈ln|Ta|〉2~−c2lna,C3(a)≡〈ln3|Ta|〉−3〈ln2|Ta|〉+2〈ln|Ta|〉3~c3lna,

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

τ(q)=−Dfq00!+∑n=1∞[Cn(a)lna]qnn!,      =−c0+c1q−c2q2/2!+c3q3/3!⋯

where the coefficients cn > 0 are estimated as the slope of Cn(a) vs. ln (n = 1, 2, 3, , …), and c0 = Df.

The implication of the above developments is that we can estimate τ(q) from the polynomial expansion of Equation (10), where the coefficients are obtained from the log-log linear regressions of the cumulants of the magnitude Cn(a) vs. ln(a) (Equation 9) (Delour et al., 2001). A quadratic log-normal τ(q) approximation would need only three such linear regressions, cn = 0, ∀n > 2.

 2.3. Monofractal vs. multifractal functions

Homogeneous monofractal signals (distributions) are signals with singularities of unique Hölder exponent H. Their τ(q) spectrum is a linear function of q with slope c1 = H (Equation 9). Monofractal scaling indeed means that the shape of the probability distribution function (pdf) of rescaled wavelet coefficients does not change across scales as expressed by the following relationship between the WTMM pdfs Pa(T) and Pa′(T) at scale a and a′ > a respectively (Arneodo et al., 2002, 2011):

Pa(T)=(a′a)−HPa′((a′a)−HT).

A nonlinear τ(q) is the signature of multifractal signals with Hölder exponent h(t) fluctuating over time t (Muzy et al., 1991, 1994; Bacry et al., 1993; Arneodo et al., 1995b, 2002, 2008). In this study, we fit the τ(q) data by the so-called log-normal quadratic approximation τ(q)=-c0+c1q-c2q2/2. The corresponding singularity spectrum has a quadratic single humped shape:

D(h)=c0−(h−c1)2/2c2,

where c0 = −τ(0) = Df is the fractal dimension of the support of singularities of E(t), c1 is the value of h that maximizes D(h), and the intermittency coefficient c2 (Delour et al., 2001) characterizes the width of the D(h) spectrum as an indication of a change in WTMM coefficient statistics across scales. If h(t) fluctuates according to a pdf ρ(h), then (Castaing et al., 1990, 1993; Arneodo et al., 1997a, 1998c, 1999):

Pa(T)=∫ρ(h)(a′a)−hPa′((a′a)−hT) dh,

meaning that the pdf at scale a can be expressed as a weighted sum of dilated pdfs at larger scales a ′ >a. Let us point out that the monofractal situation (Equation 11) is recovered when assuming that ρ(h) = δ(h − H) in Equation (13).

Note that τ(2) = c0 + 2c1 − 2c2, also called the correlation dimension, is related to the “1/f” scaling exponent of the Fourier spectral density (Muzy et al., 1994; Mandelbrot, 1998):

|E^(f)|2~f−β,  with τ(2)=β−2. 2.4. In quest of an underlying time-scale multiplicative structure: the two-points magnitude correlation method

Multiplicative cascade processes (Arneodo et al., 1998b) are paradigmatic mechanisms generating multifractal distributions, with as historical examples the Kolmogorov–Obukhov log-normal energy cascade model of fully developed turbulence (Kolmogorov, 1962; Oboukhov, 1962; Mandelbrot, 1974), and the Multifractal Random Walk (MRW) model recently introduced to account for the intermittency observed in financial time series (Muzy et al., 2000; Bacry et al., 2001). But, if multifractal scaling implies some evolution of the WTMM statistics across scales, it does however not require any correlation of the wavelet coefficients across scales. In addition to the above one-point WTMM statistics, it is thus useful to study the two-point correlation function of the logs of the WTMM coefficients ln |Ta(t)|, which determines the way the correlation structure of the Hölder exponents h (or singularities) changes with scale (Arneodo et al., 1998a,b). Defining the two-point magnitude correlation function C(a, Δt) as:

C(a,Δt)=〈(ln|Ta(t)|−〈ln|Ta(t)|〉)·(ln|Ta(t+Δt)|                 −〈ln|Ta(t)|〉)〉,

and seeing how this correlation changes as a function of Δt at scale a, provides information about the time-scale structure that underlies the multifractal properties of the considered signal. As demonstrated by Arneodo et al. (1998a,b) for random multiplicative cascades on wavelet dyadic trees (see also Meneveau and Sreenivasan, 1991):

C(a,Δt)~−c2 lnΔt,      Δt>a,

where the proportionality coefficient c2 is the intermittency coefficient defined in Equation(9) (Note that C(a, Δt = 0) ≡ C2(a) ∽ − c2ln a). Thus, by computing C(a, Δt) from Equation (15) 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).

 3. Description of data 3.1. Study design and population

The experimental data are hospital-based. We have analyzed data recorded in the atria of 8 patients with persistent or chronic AF, chosen without any prior explicit exclusion criteria. These patients were enrolled to undergo radio frequency ablation between 2010 and 2012, in the international cardiac electrophysiology service of public hospital CHU Haut-Lévêque in Pessac, France. All patients gave written informed consent to the investigation of data from the intervention. A protocol for clinic research was approved by the institutional Clinical Research and Ethics Committee. For this specific investigation of the data, the authors accessed fully anonymized and de-identified data. As representative of the results obtained with our set of patients, we report in this manuscript the results of a detailed wavelet-based multifractal analysis of five long time series specially recorded in one of the 8 patients with chronic AF.

 3.2. Electric potential recording

A steerable decapolar catheter, equipped with five 1 mm distant pairs of electrodes, each pair separated by 5 mm (Xtrem, Sorin Medical ©), was positioned in the CS as recommended in (Figures 1A,B). The distal leads of the catheter tip are positioned in a region near the left pulmonary veins, while the proximal leads lie closer to the right atrium orifice of the CS. This catheter was immobile and probed the electrical activity of those areas in the left atrium. Monitoring typically lasted the whole intervention which could take hours. Our file consists of 5 simultaneous recordings at points Pt1 to Pt5, from distal to proximal positioning along the vein (Figure 1A), each lasting 422 s with sampling time 10−3 s, a few minutes before the first ablation procedure started. The potential difference Δϕ(t) between each of the two electrodes in each pair was recorded, with the convention of distal minus proximal (Figure 1C). The normal to rapid frequency in sinus rhythm varies in the range 1Hz ≲ f ≲ 3 Hz, whereas during AF it is typically in the range 3Hz ≲ f ≲ 10Hz (Figure 1D). But the most obvious observation is that on-site recordings during AF contrast with the ones during sinus rhythm as the former seem to fluctuate randomly at even higher frequencies. Physiologically, a natural high frequency cut-off is somewhere in between 100Hz ≲ fc ≲ 1, 000 Hz, which corresponds to the shortest characteristic time scale in a cardiac cell cycle, that is depolarization. Furthermore, AF is considered as the most irregular cardiac arrhythmia (Konings et al., 1994, 1997).

Electric potential recording. (A) Sketch of the positioning of the pairs of electrodes (Pt1–Pt5) along the catheter in the coronary sinus (CS); left posterior view (LPV) of the heart left ventricle (LV), left inferior pulmonary vein (LIPV), left atrium (LA), and inferior vena cava (IVC). (B) A Left anterior fluoroscopy view (LAV) of the atria showing the pairs of electrodes in the CS and one other catheter. (C) 6 s portions of Δϕ(t) recorded at the points Pt1, Pt2, Pt3, Pt4, and Pt5, from top to bottom. (D) Corresponding Fourier power spectra computed over the whole 422 s time series.

 3.3. Local impulse energy

Each electrode averages the electric potential over its surface. During depolarization, ions flow through the cell membrane channels and the gap junction channels, inducing a rapid change of the electric potential (∽10−2s). Then, repolarization is a much smoother event. Bipolar electrodes are separated by a few millimeters which is typically the length scale of a depolarizing front in the atria with a conduction velocity c ∽ 10−1 m/s and a refractory period RP ∽ 10−1 s (Figures S2A,B), thus defining the so-called “wavelength” scale c × RP (Smeets et al., 1986; Rensma et al., 1988). Any spatio-temporal variation of the vulnerable substrate function happens over larger time scales. Thus the conduction velocity can be considered constant over such small scales and the bipolar electric potential difference is therefore “frozen.” It follows that the bipolar electric potential difference is locally advected with velocity c:

∂Δϕ(t)∂t=−c∇→Δϕ(t).

Under this assumption, the evaluation of the local electric energy of a cardiac impulse is straight-forward:

ε2ℰ2(t)=ε2c2(∂∂tΔϕ(t))2,

where E is the electric field magnitude and ε is the dielectric bulk permitivity of the (inter-) cellular medium. Energy will thus peak when the impulse travels between the two electrodes. It fully incorporates ionic flux through membrane channels and electrotonic currents, specifically those taking place at the gap junctions. Because we have no means to assess the conduction velocity c, we will use in this study the following definition of the energy:

E(t)=(∂Δϕ(t)∂t)2,

after dropping the term in front of the r.h.s. of Equation (18), i.e., ε/2c2 which remains constant at first order as long as the conduction velocity c does not fluctuate too much. To practically derive E(t) from the recorded Δϕ(t), we used an order 4 finite difference scheme on a oversampled (Δt0=10-4s) cubic-spline fitting of the data (Figures S2C,D). This is needed to estimate peaks in the energy within temporal windows as narrow as ∽10−3 s (Figures S2E,F). We have checked that the scaling properties displayed by E(t) in the low frequency range of interest here (0.08Hz ≲ f ≲ 2Hz) are not affected by this discretization scheme.

Let us note that Equation (19) is not without reminding the 1D surrogate dissipation approximation Edis(t)∽(∂υ/∂t)2, where υ is the longitudial velocity, in fully developed (homogeneous and isotropic) turbulence under the Taylor hypothesis of “frozen” turbulence (Meneveau and Sreenivasan, 1991; Frisch, 1995).

 3.4. Software and documentation

The numerical procedure to perform the WTMM analysis of 1D signals can be downloaded at

http://perso.ens-lyon.fr/benjamin.audit/LastWave

LastWave is open source software written in C. We recommend interested users to read the LastWave C-Application Programming Interface documentation and to contact the corresponding author to be directed to the part of the code of most relevance to them.

 4. Results 4.1. One-point multifractal analysis of local impulse energy data

We first present an exhaustive step-by-step analysis of one of the impulse energy time-series recorded at the electrode Pt2 as an illustration of the intricacies involved in the methodology and as a demonstration that without a priori knowledge about the signal, a reliable multifractal analysis requires an iterative process between diagnosis and estimation until robustness is achieved. In Figure S3 is illustrated how the WT microscope is able to filter out the nonstationarities (polynomial trends) in E(t) (Figure S3A) when using analyzing wavelets g(N) (Figure S1) of increasing order (Figures S3B–D).

4.1.1. Multifractal analysis of the impulse energy data with the WTMM method of moments

When applying the WTMM method to the impulse energy time-series recorded at Pt2 (Figure 2A), we revealed that the partition functions Z(q, a) (Equation 3) obtained from the WT computed with the analyzing wavelet g(3) (Figure 2B) and its skeleton (Figure 2C), display nice scaling properties for q = −1 to 5 over a range of time-scales larger than the mean interbeat ∽0.5 s (Figure 1D). We strictly limited this range to (0.6, 10 s) for linear regression fit estimates in a logarithmic representation (Figure 3A). The τ(q) so-obtained is well approximated by a quadratic spectrum with parameters [c0, c1, c2] = [1.01, −0.34, 0.053] and cn = 0 for n > 2 (Equation 10) (Figure 4A). This signature of multifractality with a support of singularities of fractal dimension Df ≈ 1, and an intermittency coefficient c2 = 0.053 ± 0.010 (Table 1) is confirmed when respectively plotting h(q, a)/ln2 (Equation 5) and D(q, a)/ln2 (Equation 6) vs. log2 a Figures 3B,C. From the estimate of the slopes h(q) and D(q), we get the single humped D(h) spectrum shown in Figure 4B which is well approximated by the quadratic spectrum defined in Equation (12) with the above parameter values obtained from a polynomial fitting of the τ(q) data. Interestingly, when comparing the results obtained with the analyzing wavelet g(3) with those obtained with g(1) and g(2) in Figures 3, 4, we notice that except some slight differences observed when using the first-order analyzing wavelet g(1), the multifractal spectra obtained with the second-order wavelet g(2) and the third-order wavelet g(3) almost superimpose (Figure 4, Table 1) with singularities of Hölder exponent h ≤ 0, characteristic of a multifractal “noise” signal.

Wavelet transform of local impulse energy time-series. (A) A 100 s portion of E(t) (Equation 19) recorded at the electrode Pt2. (B) Time-scale WT representation of E(t) with the analyzing wavelet g(3) (Figure S1). The modulus of the WT is coded, independently at each scale a, using 256 colors from black (|Tg(3)(t,a)|=0) to red (maxt|Tg(3)(t,a)|). (C) WT skeleton defined by the maxima lines. The scale a = Δt/Δt0, where Δt0=10-4 s. In (B) the white horizontal dotted lines delimit the range of time scales (29 ≤ a ≤ 213) used to perform linear regression fit estimates of the τ(q) and D(h) multifractal spectra.

Multifractal analysis of local impulse energy time-series recorded at the electrode Pt2 with the WTMM method. (A) log2 Z(q, a) vs. log2 a (Equation 3). (B) h(q, a)/ln 2 vs. log2 a (Equation 5). (C) D(q, a)/ln 2 vs. log2 a (Equation 6). The computation were performed for different values from q = −1 to 5 with the analyzing wavelet g(1) (▿), g(2) (□) and g(3) (◦) (Figure S1). The vertical dashed lines delimit the range of scales (29 ≤ a ≤ 213) used for the linear regression estimate of τ(q), h(q) and D(q) in Figure 4.

Multifractal spectra of local impulse energy time-series (Pt2) obtained with the WTMM method of moments. (A) τ(q) vs. q estimated by linear regression fit of log2 Z(q, a) vs. log2 a (Figure 3A). (B) D(h) vs. h obtained from linear regression fits of h(q, a) (Figure 3B) and D(q, a) (Figure 3C) vs. log2 a. The symbols correspond to the analyzing wavelets g(1) (▿), g(2) (□) and g(3) (◦) (Figure S1). The curves correspond to quadratic spectra (Equations 10 and 12) with parameters [c0, c1, c2] = [1.01, −0.33, 0.055] (⋯ , g(1)), [1.00, −0.35, 0.049] (−−−, g(2)), [1.01, −0.34, 0.053] (—, g(3)) (see Table 1).

Results of the WTMM multifractal analysis of the local impulse energy time-series recorded along the CS vein at electrodes Pt1, Pt2, Pt3, and Pt5.

g(1) g(2) g(3)
Point1
c0 0.961 ± 0.001 0.995 ± 0.001 1.009 ± 0.002
c1 −0.351 ± 0.003 −0.298 ± 0.005 −0.281 ± 0.007
c1* −0.353 ± 0.017 −0.297 ± 0.011 −0.274 ± 0.011
c2 0.048 ± 0.006 0.063 ± 0.011 0.064 ± 0.014
c2* 0.050 ± 0.032 0.082 ± 0.018 0.096 ± 0.020
Point2
c0 1.011 ± 0.002 0.998 ± 0.001 1.011 ± 0.001
c1 −0.333 ± 0.007 −0.348 ± 0.003 −0.337 ± 0.005
c1* −0.307 ± 0.012 −0.335 ± 0.011 −0.331 ± 0.011
c2 0.055 ± 0.013 0.049 ± 0.007 0.053 ± 0.010
c2* 0.076 ± 0.030 0.031 ± 0.027 0.047 ± 0.028
Point3
c0 1.023 ± 0.003 1.005 ± 0.002 1.021 ± 0.003
c1 −0.472 ± 0.014 −0.496 ± 0.008 −0.481 ± 0.011
c1* −0.445 ± 0.011 −0.480 ± 0.008 −0.464 ± 0.007
c2 0.082 ± 0.028 0.091 ± 0.015 0.098 ± 0.022
c2* 0.164 ± 0.031 0.103 ± 0.023 0.164 ± 0.020
Point5
c0 1.044 ± 0.003 1.017 ± 0.003 1.029 ± 0.002
c1 −0.320 ± 0.014 −0.365 ± 0.011 −0.383 ± 0.009
c1* −0.256 ± 0.015 −0.335 ± 0.018 −0.365 ± 0.025
c2 0.176 ± 0.028 0.167 ± 0.022 0.152 ± 0.018
c2* 0.175 ± 0.043 0.106 ± 0.029 0.114 ± 0.031

c0, c1, c2 are the coefficients of the polynomial expansion of τ(q) (Equation 10) obtained with the WTMM method of moments when using the analyzing wavelets g(1), g(2), and g(3) respectively (Figure S1). c1* and c2* are the corresponding coefficients obtained with the magnitude cumulant method.

This multifractal diagnosis is confirmed in Figure 5 where the WTMM pdfs obtained at different scales with g(3) (Figure 5A) are shown to collapse on each other when using the propagative equation of the statistics across scales (Equation 13) with the quadratic τ(q) spectrum estimated just above (Figure 5B) (Castaing et al., 1990, 1993; Arneodo et al., 1997a, 1998c, 1999; Venugopal et al., 2006).

Demonstration of the WTMM pdf rescaling via the propagative equation across scales (Equation 13). (A) Original pdfs of the logs of WTMM coefficients (magnitude coefficients) of the local impulse energy time series recorded at the electrode Pt2, computed with the analyzing wavelet g(3), at scales a (= Δt/Δt0, where Δt0=10-4 s) = 29(*), 210(•), 211(▴) and 212(■). (B) Rescaled pdfs using the multifractal quadratic estimate of the τ(q) spectrum (Equation 10) with parameters [c0, c1, c2] = [1.01, −0.34, 0.053] (see Table 1).

 4.1.2. Multifractal analysis of the impulse energy data with the method of magnitude cumulants

After the WTMM partition function approach, we turn our attention to the alternate magnitude cumulant analysis methodology. The first-, second- and third-order cumulants were computed using Equation (9) and are plotted vs. the logarithm of the scale in Figure 6. As expected C1(a), C2(a) and C3(a) display consistent scaling behavior over the same range of scales 29 ≤ a ≤ 213 (a = Δt/Δt0, where the oversampling time is Δt0=10-4 s) and this for the three analyzing wavelets g(1), g(2) and g(3). The results obtained for C3(a) (Figure 6C) confirm that with the limited statistical sample at our disposal (422 s long time series), there is no way to conclude about the possible departure from a log-normal quadratic τ(q) spectrum (c3 ≡ 0). Nicely, the quadratic τ(q) spectrum obtained with g(3) with parameters c1*=-0.33±0.01 and c2*=0.047±0.028 is found in good agreement with the one previously estimated with the method of moments, confirming the multifractal diagnosis of the local impulse energy at low frequencies. Let us point out that, as reported in Table 1, the τ(q) spectrum obtained with g(1) is again slightly different from the ones obtained with g(2) and g(3) which turn out to be indistinguishable. This is the numerical demonstration that a robust estimate of the multifractal spectra is achieved when using the third-order analyzing wavelet g(3).

Magnitude cumulant analysis of local impulse energy time-series (Pt2). (A) C1(a)/ln 2 vs. log2 a. (B) C2(a)/ln 2 vs. log2 a. (C) C3(a)/ln 2 vs. log2 a. The computation of the Cn(a) (Equation 9) was performed with the analyzing wavelets g(1) (▿), g(2) (□) and g(3) (◦) (Figure S1). The vertical dashed lines delimit the range of scales (29 ≤ a ≤ 213) used for the linear regression estimate of coefficients c1*, c2* and c3* of τ(q) (Equation 10) reported in Table 1.

 4.1.3. Summary of one-point multifractal analysis of the impulse energy data

Similar analysis was performed on the four other local impulse energy time-series recorded at electrodes Pt1 (Figures S4, S8, S11), Pt3 (Figures S5, S9, S12), Pt4 (Figures S6, S13), and Pt5 (Figures S7, S10, S14). As already noticeable on the power spectrum in Figure 1D, the electric potential Δϕ(t) recorded at electrode Pt4 presents some subharmonic oscillatory component at frequency ∽2.5 Hz that dramatically spoils the scaling behavior previously obtained at Pt2 over the range of time-scales (0.6, 10 s) (Figures S6, S13). We can speculate on the signature of the nearby mitral valve influence. This explains that Pt4 will be singled out in the rest of our study. Figure 7 displays the τ(q) (Figure 7A) and D(h) (Figure 7B) spectra of the local impulse energy E(t) obtained with the WTMM method of moments when using the analyzing wavelet g(3) (the corresponding spectra obtained with analyzing wavelets g(1) and g(2) are shown in Figures S15, S16 respectively). The spectra obtained for the time-series recorded at electrode Pt1 situated, as Pt2, in the ligament of Marshall anatomic area innervated by the ANS (Figure 1A), are very similar to the ones observed at Pt2. Both τ(q) and D(h) spectra are well approximated by quadratic spectra (Equations 9 and 10, respectively) with parameters [c0, c1, c2] = [1.01, −0.28, 0.064] (Figure 7) with a definitely positive finite intermittency coefficient c2 = 0.064 ± 0.014 (Table 1). Interestingly, time series recorded at electrodes Pt3 and Pt5 in a different anatomical area next to the left atrial posterior wall (Figure 1A) both show rather similar τ(q) and D(h) multifractal spectra but significantly different from the ones obtained at electrodes Pt1 and Pt2 (Figure 7). Again these spectra are found nearly quadratic with parameters [c0, c1, c2] = [1.02, −0.48, 0.098] for Pt3, and [1.03, −0.38, 0.152] for Pt5 (Table 1). Local impulse Local impulse energy time series at Pt3 and Pt5 show higher intermittency with larger c2 values whereas they display weaker long-range correlations c1 = < h > ∽ − 0.45 (i.e., closer to c1 = −0.5 a value characteristic of uncorrelated white noise), instead of the value c1∽ − 0.3 for Pt1 and Pt2 characteristic of positive long-range correlations. As reported in Table 1, this regionalization of the multifractal properties of the local impulse energy is quantitatively confirmed when using the magnitude cumulant method. It is further corroborated when reproducing this multifractal analysis for other patients with paroxysmal, persistent or chronic AF (Figure 8A) and for a patient at different periods of time preceding ablation procedure (Figure 8B) as an indication of stationarity.

Multifractal spectra of local impulse energy time-series recorded along the CS vein. (A) τ(q) vs. q estimated by linear regression fit of log2 Z(q, a) vs. log2 a. (B) D(h) vs. h obtained from linear regression fits of h(q, a) and D(q, a) vs. ln2 a. The analyzing wavelets is g(3). The colored symbols correspond to the electrodes Pt1 (black), Pt2 (red), Pt3 (blue) and Pt5 (green). The curves correspond to quadratic spectra (Equations 10 and 12) with parameters [c0, c1, c2] = [1.01, −0.28, 0.064] (black, Pt1), [1.01, −0.34, 0.053] (red, Pt2), [1.02, −0.48, 0.098] (blue, Pt3), and [1.03, −0.38, 0.152] (green, Pt5) (see Table 1).

τ(q) spectra of local impulse energy time-series recorded along the CS vein at the electrodes Pt2 (red), Pt3 (blue) and Pt5 (green). The curves represent quadratic polynomial fit of the data (Equation 10). (A) The symbols correspond to the reference Patient 1 (chronic AF, ▿) and to Patients 2 (chronic AF, ◦), 3 (paroxysmal AF, □) and 4 (persistent AF, △). (B) The symbols correspond to the reference Patient 1 (▿) and to three different time-series for Patient 4 (◦, □, △) recorded at different periods of time preceding ablation procedure.

 4.2. 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 recorded at the positions Pt1, Pt2, Pt3, and Pt5 along the CS vein are shown in Figure 9. C(a, Δt) (Equation 15) computed with the analyzing wavelet g(3) is represented vs. Δt for two scales a = 29 and 210 in the scaling range. Strikingly for all four time series, for Δt ≳ a, C(a, Δt) drops to zero as a clear indication that the magnitudes are uncorrelated. As a reference, we put in each panel of Figure 9, a dashed straight line of slope −c2 as predicted by Equation (16) for multifractal signals exhibiting a cascading multiplicative structure along a time-scale tree (Arneodo et al., 1998b). The slow decay predicted by the “multiplicative” log-normal model with intermittency coefficient c2 is definitely not observed. Thus, local impulse energy time-series look much more like what has been called log-normal “mutlifractal white noise” in pioneering works to distinguish “uncorrelated” and “multiplicative” log-normal models (Arneodo et al., 1998a). A similar absence of magnitude correlation is observed when reproducing this two-point magnitude analysis with the analyzing wavelets g(1) (Figure S17) and g(2) (Figure S18).

Two-point magnitude analysis of local impulse energy time-series recorded along the CS vein. Two-point correlation function C(a, Δt)/C(a, 0) vs. ln (Δt) (Equation 15) for local impulse energy E(t) computed with the analyzing wavelet g(3). The two curves correspond to scales a = 29 (black) and 210 (gray) within the scaling range. (A) Pt1, (B) Pt2, (C) Pt3, and (D) Pt5.

 5. Discussion

To summarize, we showed that the wavelet-based multifractal analysis of long time series of the local impulse energy recorded in the CS of a patient with chronic AF was able to reveal and quantify the intermittent nature of these signals at low frequency (f ≲ 2 Hz). To our knowledge, our study is the first to report on the observation and quantification of such multifractal dynamics of the endocavitary electrical activity during chronic AF which is found more complex than previously suspected. On the basis of the analysis of the time-series recorded at 4 catheter electrodes out of 5 positioned in the CS, two main observations can be made: (i) the local impulse energy displays different multifractal properties in the left atrial wall area than in the ligament of Marshall area consistently with different anatomical substrate conditions, and (ii) while recorded along the CS vein, the local impulse energy does not exhibit long-range dependence associated with an underlying multiplicative cascade, or in other words the multifractal distribution of the singularities inferred by the two-point magnitude analysis does not display any correlation across scales just like a log-normal “multifractal white noise” (Arneodo et al., 1998a). The nature of this study was exploratory, with a data set limited to a few patients, and with a few time series rather long for clinical practice (422 s) but not so long regarding the range of time scales [0.6, 10 s] where scaling was observed. This is the reason for the different complementary analyses employed in this paper including the WTMM method of moments, the WTMM method of magnitude cumulants, and the two-points magnitude cumulant method, using analyzing wavelets of different orders, until reliable estimates were obtained. Of course our results deserve to be confirmed over a large set of patients at different stages of AF development and to be explored in different areas of the atria. The goal would be to exhibit precise determinants of the diseased substrate using multifractal scaling analysis. But this preliminary analysis definitely challenges current knowledge in physical, physiological and clinical fundamentals of AF arrhythmia. Specifically, it challenges the mechanistic approach of AF based on circuit reentries.

The absence of an underlying cascading process is not such a surprise since underlying the multifractal properties displayed by the local impulse energy at low frequencies (f ≲ 2 Hz), there is no clear 3D “fragmentation” (Mandelbrot, 1982) process inducing some cascading of energy from large to small time scales and also no obvious 2D “aggregation, coalescence or growth” (Vicsek, 1989) process bringing energy from small to large time scales. What are the physical and physiological mechanisms that drive the multifractal nature of local impulse energy and give rise to the observed differences according to area is still an open question. Nonetheless, these results already undermine the commonly accepted concepts revolving around circuit reentries, and a fortiori spiral waves, as being basic mechanisms for the onset and perpetuation of AF. The mechanistic “wavelength” criterion indeed conveys the idea that random spatio-temporal dispersion of refractoriness, or more generally of functional properties, leads to random mixing of circuit reentries. The “wavelength” scale adjusts naturally to the typical scale λ of dispersion when it exists c × RP ≲ λ, as would be the case for Gaussian statistics of dispersion. In that case, the statistics of the local impulse energy remains Gaussian throughout scales. On the contrary, to fit our new observations we have seen that the statistics is not Gaussian and evolves across scales through a log-normal propagation law which accounts for the intermittency observed over the range of a few beat cycles (~0.6 s) to several tens (~10 s) (and possibly more), therefore spanning the whole atria. Although the ligament of Marshall area is highly innervated (Tan et al., 2007; Ulphani et al., 2007; Arora, 2012), it is quite unlikely that modulations by the ANS, that affects heart rate, play a significant role in the intermittent dynamics since the documented three peak frequencies at ∽0.4, ∽0.15, and ∽0.04 Hz (Akselrod et al., 1981) do not show up in our analysis. Furthermore, we have found at least two areas with different multifractal regimes. Thus, our findings raise new challenging questions calling for ongoing efforts to develop physiological heart tissue models that account for the low frequency intermittent nature of local impulse energy. Recent studies in animal models suggest the protective role of connexin gene transfer to prevent sustained AF (Bikou et al., 2011; Igarashi et al., 2012). In this spirit, in a companion modeling paper (Attuel et al., submitted), we propose a model of gap junction conduction remodeling in a denervated heart that accounts for the observed intermittent dynamics over large time scales, as resulting from incoherent random back scatterings leading to the desynchronization of the network of cardiac excitable cells.

 Author contributions

Conception and design: GA, HY, AA. Development and methodology: EG-C, FA, AA. Analysis and interpretation of data: GA, EG-C, FA, HY, AA. Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): GA, HY. Writing, review, and/or revision of the manuscript: GA, EG-C, FA, HY, AA. Administrative, technical or material support (i.e., requiring and organizing data, constructing databases): GA, EG-C.

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.

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.

 Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2017.01139/full#supplementary-material

 References Akselrod S. Gordon D. Ubel F. A. Shannon D. C. Barger A. C. Cohen R. J. (1981). Power spectrum analysis of heart rate fluctuations : a quantitative probe of beat-to-beat cardiovascular control. Science 213, 220–222. 10.1126/science.6166045 Al-Khatib S. M. LaPointe N. M. A. Chatterjee R. Crowley M. J. Dupre M. E. Kong D. F. . (2014). Rate-and rhythm-control therapies in patients with atrial fibrillation: a systematic review. Ann. Int. Med. 160, 760–773. 10.7326/M13-146724887617 Allessie M. A. (1998). Atrial electrophysiologic remodeling: another vicious circle? J. Cardiovasc. Electrophysiol. 9, 1378–1393. 10.1111/j.1540-8167.1998.tb00114.x9869538 Allessie M. A. Bonke F. I. Schopman F. J. (1977). Circus movement in rabbit atrial muscle as a mechanism of tachycardia. Circ. Res. 41, 9–18. 10.1161/01.RES.41.1.9862147 Antoine J. P. Murenzi R. Vandergheynst P. Ali S. T. (2008). Two-Dimensional Wavelets and their Relatives. Cambridge: Cambridge University Press. Arneodo A. Audit B. Decoster N. Muzy J.-F. 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 Bunde A. Kropp J. Schellnhuber H. J. (Berlin: Springer Verlag), 26–102. Arneodo A. Audit B. Kestener P. Roux S. G. (2008). Wavelet-based multifractal analysis. Scholarpedia 3:4103. 10.4249/scholarpedia.4103 Arneodo A. Bacry E. Jaffard S. Muzy J.-F. (1997a). Experimental analysis of self-similarity and random cascade processes: application to fully developed turbulence data. J. Phys. II 7, 363–370. 10.1051/jp2:1997130 Arneodo A. Bacry E. Jaffard S. Muzy J.-F. (1997b). Oscillating singularities on Cantor sets: a grand-canonical multifractal formalism. J. Stat. Phys. 87, 179–209. 10.1007/BF02181485 Arneodo A. Bacry E. Manneville S. Muzy J.-F. (1998a). Analysis of random cascades using space-scale correlation functions. Phys. Rev. Lett. 80, 708–711. 10.1103/PhysRevLett.80.708 Arneodo A. Bacry E. Muzy J.-F. (1995a). Oscillating singularities in locally self-similar functions. Phys. Rev. Lett. 74, 4823–4826. 10.1103/PhysRevLett.74.482310058608 Arneodo A. Bacry E. Muzy J.-F. (1995b). The thermodynamics of fractals revisited with wavelets. Physica A 213, 232–275. 10.1016/0378-4371(94)00163-N Arneodo A. Bacry E. Muzy J.-F. (1998b). Random cascades on wavelet dyadic trees. J. Math. Phys. 39, 4142–4164. 10.1063/1.532489 Arneodo A. Decoster N. Kestener P. Roux S. G. (2003). A wavelet-based method for multifractal image analysis: from theoretical concepts to experimental applications. Adv. Imaging Electr. Phys. 126, 1–92. 10.1016/S1076-5670(03)80014-9 Arneodo A. Decoster N. Roux S. G. (2000). A wavelet-based method for multifractal image analysis. I. Methodology and test applications on isotropic and anisotropic random rough surfaces. Eur. Phys. J. 15, 567–600. 10.1007/s100510051161 Arneodo A. Grasseau G. Holschneider M. (1988). Wavelet transform of multifractals. Phys. Rev. Lett. 61, 2281–2284. 10.1103/PhysRevLett.61.228110039072 Arneodo A. Grasseau G. Kostelich E. J. (1987). Fractal dimensions and f(α) spectrum of the Hénon attractor. Phys. Lett. A 124, 426–432. 10.1016/0375-9601(87)90546-9 Arneodo A. Manneville S. Muzy J.-F. (1998c). Towards log-normal statistics in high Reynolds number turbulence. Eu. Phys. J. B 1, 129–140. 10.1007/s100510050162 Arneodo A. Manneville S. Muzy J.-F. Roux S. G. (1999). Revealing a lognormal cascading process in turbulent velocity statistics with wavelet analysis. Philos. Trans. R. S. Lond. A 357, 2415–2438. 10.1098/rsta.1999.0440 Arneodo A. Muzy J.-F. Sornette D. (1998d). Direct” causal cascade in the stock market. Eur. Phys. J. B 2, 277–282. 10.1007/s100510050250 Arneodo A. Vaillant C. Audit B. Argoul F. d'Aubenton-Carafa Y. Thermes C. (2011). Multi-scale coding of genomic information: from DNA sequence to genome structure and function. Phys. Rep. 498, 45–188. 10.1016/j.physrep.2010.10.001 Arora R. (2012). Recent insights into the role of the autonomic nervous system in the creation of substrate for atrial fibrillation. Circ. Arrhythm. Electrophysiol. 5, 859–859. 10.1161/CIRCEP.112.97227322895601 Attuel P. Childers R. Cauchemez B. Poveda J. Mugica J. Coumel P. (1982). Failure in the rate adaptation of the atrial refractory period: its relationship to vulnerability. Int. J. Cardiol. 2, 179–197. 10.1016/0167-5273(82)90032-86185445 Attuel P. Pellerin D. Gaston J. (1989). Latent atrial vulnerability: new means of electrophysiologic investigations in paroxysmal atrial arrhythmias, in The Atrium in Health and Disease, eds Attuel P. Coumel P. Janse M. (Mount Kisco: Mount KiscoFutura Publishing Co., Inc.), 159–200. Attuel P. Rancurel G. Delgatte B. Golcher E. Ghazoullieres P. Friocourt P. . (1986). Importance of atrial electrophysiology in the work-up of cerebral ischemic attacks. PACE 9, 1121–1126. 10.1111/j.1540-8159.1986.tb06680.x2432519 Audit B. Bacry E. Muzy J.-F. Arneodo A. (2002). Wavelet-based estimators of scaling behavior. IEEE Trans. Info. Theory 48, 2938–2954. 10.1109/TIT.2002.802631 Audit B. Baker A. Chen C.-L. Rappailles A. Guilbaud G. Julienne H. . (2013). Multiscale analysis of genome-wide replication timing profiles using a wavelet-based signal-processing algorithm. Nat. Protoc. 8, 98–110. 10.1038/nprot.2012.14523237832 Bacry E. Delour J. Muzy J. F. (2001). Multifractal random walk. Phys. Rev. E 64:026103. 10.1103/PhysRevE.64.02610311497647 Bacry E. Muzy J.-F. Arneodo A. (1993). Singularity spectrum of fractal signals from wavelet analysis: exact results. J. Stat. Phys. 70, 635–674. 10.1007/BF01053588 Batchelder K. A. Tanenbaum A. B. Albert S. Guimond L. Kestener P. Arneodo A. . (2014). Wavelet-based 3D reconstruction of microcalcification clusters from two mammographic views: new evidence that fractal tumors are malignant and euclidean tumors are benign. PLoS ONE 9:e107580. 10.1371/journal.pone.010758025222610 Bertaglia E. Tondo C. De Simone A. Zoppo F. Mantica M. Turco P. . (2010). Does catheter ablation cure atrial fibrillation? Single-procedure outcome of drug-refractory atrial fibrillation ablation: a 6-year multicentre experience. Europace 12, 181–187. 10.1093/europace/eup34919887458 Bikou O. Thomas D. Trappe K. Lugenbiel P. Kelemen K. Koch M. . (2011). Connexin 43 gene therapy prevents persistent atrial fibrillation in a porcine model. Cardiovasc. Res. 92, 218–225. 10.1093/cvr/cvr20921799069 Bohr T. 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. Camm A. J. Kirchhof P. Lip G. Y. Schotten U. Savelieva I. Ernst S. . (2010). Guidelines for the management of atrial fibrillation. Eur. Heart J. 31, 2369–2429. 10.1093/eurheartj/ehq27820802247 Castaing B. Gagne Y. Hopfinger P. J. (1990). Velocity probability density functions of high Reynolds number turbulence. Physica D 46, 177–200. 10.1016/0167-2789(90)90035-N Castaing B. Gagne Y. Marchand M. J. (1993). Log-similarity for turbulent flows? Physica D 68, 387–400. 10.1016/0167-2789(93)90132-K Chudácek V. Andén J. Mallat S. Abry P. Doret M. (2014). Scattering transform for intrapartum fetal heart rate variability fractal analysis: a case-control study. IEEE Trans. Biomed. Eng. 61, 1100–1108. 10.1109/TBME.2013.229432424658235 Ciuciu P. Varoquaux G. Abry P. Sadaghiani S. Kleinschmidt A. (2012). Scale-free and multifractal time dynamics of fMRI signals during rest and task. Front. Physiol. 3:186. 10.3389/fphys.2012.0018622715328 Coumel P. (1994). Paroxysmal atrial fibrillation: a disorder of autonomic tone? Eur. Heart J. 15, 9–16. 10.1093/eurheartj/15.suppl_A.98070496 Coumel P. Attuel P. Lavallée J. Flammang D. Leclercq J. Slama R. (1978). The atrial arrhythmia syndrome of vagal origin. Arch. Mal. Coeur. Vaiss. 71, 645–656. 28709 Daubechies I. (1992). Ten Lectures on Wavelets. CBMS-NSF Regional Conference Series in Applied Mathematics, vol. 61. Philadelphia, PA: SIAM. Decoster N. Roux S. G. Arneodo A. (2000). A wavelet-based method for multifractal image analysis. II. Applications to synthetic multifractal rough surfaces. Eur. Phys. J. 15, 739–764. 10.1007/s100510051179 Delour J. Muzy J.-F. Arneodo A. (2001). Intermittency of 1D velocity spatial profiles in turbulence: a magnitude cumulant analysis. Eur. Phys. J. B 23, 243–248. 10.1007/s100510170074 Echt D. S. Liebson P. R. Mitchell L. B. Peters R. W. Obias-Manno D. Barker A. H. . (1991). Mortality and morbidity in patients receiving encainide, flecainide, or placebo. The Cardiac Arrhythmia Suppression Trial. N. Engl. J. Med. 324, 781–788. 10.1056/NEJM199103213241201 Frisch U. (1995). Turbulence: The Legacy of Kolmogorov. New York, NY: Cambridge University Press. Ganesan A. N. Shipp N. J. Brooks A. G. Kuklik P. Lau D. H. Lim H. S. . (2013). Long-term outcomes of catheter ablation of atrial fibrillation: a systematic review and meta-analysis. J. Am. Heart Assoc. 2:e004549. 10.1161/JAHA.112.00454923537812 Gerasimova E. Audit B. Roux S. G. Khalil A. Gileva O. Argoul F. . (2014). Wavelet-based multifractal analysis of dynamic infrared thermograms to assist in early breast cancer diagnosis. Front. Physiol. 5:176. 10.3389/fphys.2014.0017624860510 Gerasimova-Chechkina E. Toner B. Marin Z. Audit B. Roux S. G. Argoul F. . (2016). Comparative multifractal analysis of dynamic infrared thermograms and X-ray mammograms enlightens changes in the environment of malignant tumors. Front. Physiol. 7:336. 10.3389/fphys.2016.0033627555823 Goldberger A. L. Amaral L. A. N. Hausdorff J. M. Ivanov P. C. Peng C.-K. Stanley H. E. (2002). Fractal dynamics in physiology: alterations with disease and aging. Proc. Natl. Acad. Sci. U.S.A. 99, 2466–2472. 10.1073/pnas.01257949911875196 Goody M. F. Kelly M. W. Lessard K. N. Khalil A. Henry C. A. (2010). Nrk2b-mediated NAD+ production regulates cell adhesion and is required for muscle morphogenesis in vivo: Nrk2b and NAD+ in muscle morphogenesis. Dev. Biol. 344, 809–826. 10.1016/j.ydbio.2010.05.51320566368 Grant J. Verrill C. Coustham V. Arneodo A. Palladino F. Monier K. . (2010). Perinuclear distribution of heterochromatin in developing C. elegans embryos. Chrom. Res. 18, 873–885. 10.1007/s10577-010-9175-221116703 Grossmann A. Morlet J. (1984). Decomposition of Hardy functions into square integrable wavelets of constant shape. SIAM J. Math. Anal. 15, 723–736. 10.1137/0515056 Haïssaguerre M. Jaïs P. Shah D. C. Takahashi A. Hocini M. Quiniou G. . (1998). Spontaneous initiation of atrial fibrillation by ectopic beats originating in the pulmonary veins. N. Engl. J. Med. 339, 659–666. 10.1056/NEJM1998090333910039725923 Holschneider M. (1988). On the wavelet transformation of fractal objects. J. Stat. Phys. 50, 963–993. 10.1007/BF01019149 Igarashi T. Finet J. E. Takeuchi A. Fujino Y. Strom M. Greener I. D. . (2012). Connexin gene transfer preserves conduction velocity and prevents atrial fibrillation. Circulation 125, 216–225. 10.1161/CIRCULATIONAHA.111.05327222158756 Ito H. Glass L. (1991). Spiral breakup in a new model of discrete excitable media. Phys. Rev. Lett. 66, 671–674. 10.1103/PhysRevLett.66.67110043869 Ivanov P. C. Amaral L. A. N. Goldberger A. L. Havlin S. Rosenblum M. G. Stanley H. E. . (2001). From 1/f noise to multifractal cascades in heartbeat dynamics. Chaos 11, 641–652. 10.1063/1.139563112779503 Ivanov P. C. Amaral L. A. N. Goldberger A. L. Havlin S. Rosenblum M. G. Struzik Z. R. . (1999). Multifractality in human heartbeat dynamics. Nature 399, 461–465. 10.1038/2092410365957 Iwasaki Y.-K. Nishida K. Kato T. Nattel S. (2011). Atrial fibrillation pathophysiology: implications for management. Circulation 124, 2264–2274. 10.1161/CIRCULATIONAHA.111.01989322083148 Jaffard S. (1989). Hölder exponents at given points and wavelet coefficients. C. R. Acad. Sci. Paris Ser. I 308, 79–81. Jaffard S. (1997a). Multifractal formalism for functions part I: results valid for all functions. SIAM J. Math. Anal. 28, 944–970. 10.1137/S0036141095282991 Jaffard S. (1997b). Multifractal formalism for functions part II: self-similar functions. SIAM J. Math. Anal. 28, 971–998. 10.1137/S0036141095283005 Jaffard S. Lashermes B. Abry P. (2007). Wavelet leaders in multifractal analysis, in Wavelet Analysis and Applications, eds Qian T. Vai M. I. Xu Y. (Basel: Birkhäuser Basel), 201–246. 10.1007/978-3-7643-7778-6_17 Janse M. J. (1997). Why does atrial fibrillation occur? Eur. Heart J. 18, C12–C18. 10.1093/eurheartj/18.suppl_C.129152670 Karma A. (1993). Spiral breakup in model equations of action potential propagation in cardiac tissue. Phy. Rev. Lett. 71, 1103–1106. 10.1103/PhysRevLett.71.110310055449 Kestener P. Arneodo A. (2003). Three-dimensional wavelet-based multifractal method: the need for revisiting the multifractal description of turbulence dissipation data. Phys. Rev. Lett. 91:194501. 10.1103/PhysRevLett.91.19450114611583 Kestener P. Arneodo A. (2004). Generalizing the wavelet-based multifractal formalism to random vector fields: application to three-dimensional turbulence velocity and vorticity data. Phys. Rev. Lett. 93:044501. 10.1103/PhysRevLett.93.04450115323764 Kestener P. Conlon P. A. Khalil A. Fennell L. McAteer R. Gallagher P. T. . (2010). Characterizing complexity in solar magnetogram data using a wavelet-based segmentation method. Astrophys. J. 717:995. 10.1088/0004-637X/717/2/995 Kestener P. Lina J.-M. Saint-Jean P. Arneodo A. (2001). Wavelet-based multifractal formalism to assist in diagnosis in digitized mammograms. Image Anal. Stereol. 20, 169–174. 10.5566/ias.v20.p169-174 Khalil A. Aponte C. Zhang R. Davisson T. Dickey I. Engelman D. . (2009). Image analysis of soft-tissue in-growth and attachment into highly porous alumina ceramic foam metals. Med. Eng. Phys. 31, 775–783. 10.1016/j.medengphy.2009.02.00719297233 Khalil A. Grant J. L. Caddle L. B. Atzema E. Mills K. D. Arneodo A. (2007). Chromosome territories have a highly nonspherical morphology and nonrandom positioning. Chrom. Res. 15, 899–916. 10.1007/s10577-007-1172-817926137 Khalil A. Joncas G. Nekka F. Kestener P. Arnéodo A. (2006). Morphological analysis of HI features. II. Wavelet-based multifractal formalism. Astrophys. J. Suppl. Ser. 165, 512–550. 10.1086/505144 Kolmogorov A. N. (1962). A refinement of previous hypotheses concerning the local structure of turbulence in a viscous incompressible fluid at high Reynolds number. J. Fluid Mech. 13, 82–85. 10.1017/S0022112062000518 Konings K. T. Kirchhof C. J. Smeets J. R. Wellens H. J. Penn O. C. Allessie M. A. (1994). High-densiy mapping of electrically induced atrial fibrillation in humans. Circulation 89, 1665–1680. 10.1161/01.CIR.89.4.16658149534 Konings K. T. Smeets J. L. Penn O. C. Wellens H. J. Allessie M. A. (1997). Configuration of unipolar atrial electrograms during electrically induced atrial fibrillation in humans. Circulation 95, 1231–1241. 10.1161/01.CIR.95.5.12319054854 Lip G. Y. Lane D. A. (2015). Stroke prevention in atrial fibrillation: a systematic review. JAMA 313, 1950–1962. 10.1001/jama.2015.436925988464 Mallat S. (1998). A Wavelet Tour of Signal Processing. New York, NY: Academic Press. Mallat S. Hwang W. L. (1992). Singularity detection and processing with wavelets. IEEE Trans. Inf. Theory 38, 617–643. 10.1109/18.119727 Mandelbrot B. B. (1974). Intermittent turbulence in self similar cascades; divergence of high moments and dimension of the carrier. J. Fluid Mech. 62, 331–358. 10.1017/S0022112074000711 Mandelbrot B. (1982). The Fractal Geometry of Nature. San Francisco, CA: Freeman. Mandelbrot B. (1998). Multifractal Noise and 1/f Noise : Wild Self-Affinity in Physics. New York, NY: Springer. Marin Z. Batchelder K. A. Toner B. C. Guimond L. Gerasimova-Chechkina E. Harrow A. R. . (2017). Mammographic evidence of microenvironment changes in tumorous breasts. Med. Phys. 44, 1324–1336. 10.1002/mp.1212028112408 Martinez-Torres C. Berguiga L. Streppa L. Boyer-Provera E. Schaeffer L. Elezgaray J. . (2014). Diffraction phase microscopy: retrieving phase contours on living cells with a wavelet-based space-scale analysis. J. Biomed. Opt. 19:036007. 10.1117/1.JBO.19.3.03600724615643 Martinez-Torres C. Laperrousaz B. Berguiga L. Boyer-Provera E. Elezgaray J. Nicolini F. E. . (2015). Deciphering the internal complexity of living cells with quantitative phase microscopy: a multiscale approach. J. Biomed. Opt. 20:096005. 10.1117/1.JBO.20.9.09600526334978 McAteer R. T. J. Kestener P. Arneodo A. Khalil A. (2010). Automated detection of coronal loops using a wavelet transform modulus maxima method. Solar Phys. 262, 387–397. 10.1007/s11207-010-9530-7 Meneveau C. Sreenivasan K. (1991). The multifractal nature of turbulent energy dissipation. J. Fluid Mech. 224, 429–484. 10.1017/S0022112091001830 Meyer Y. (1992). Wavelets and Operators. Cambridge: Cambridge University Press. Middlekauff H. R. Stevenson W. G. Stevenson L. W. (1991). Prognostic significance of atrial fibrillation in advanced heart failure. Circulation 84, 40–48. 10.1161/01.CIR.84.1.402060110 Misier A. R. Opthof T. van Hemel N. M. Defauw J. J. de Bakker J. M. Janse M. J. . (1992). Increased dispersion of “refractoriness” in patients with idiopathic paroxysmal atrial fibrillation. J. Am. Coll. Cardiol. 19, 1531–1535. 10.1016/0735-1097(92)90614-S1593049 Moe G. K. Abildskov J. A. (1959). Atrial fibrillation as a self-sustaining arrhythmia independent of focal discharge. Am. Heart J. 58, 59–70. 10.1016/0002-8703(59)90274-113661062 Moe G. K. Rheinboldt W. C. Abildskov J. A. (1964). A computer model of atrial fibrillation. Am. Heart J. 67, 200–220. 10.1016/0002-8703(64)90371-014118488 Muzy J.-F. Bacry E. Arneodo A. (1991). Wavelets and multifractal formalism for singular signals: application to turbulence data. Phys. Rev. Lett. 67, 3515–3518. 10.1103/PhysRevLett.67.351510044755 Muzy J.-F. Bacry E. 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. 10.1103/PhysRevE.47.8759960082 Muzy J.-F. Bacry E. Arneodo A. (1994). The multifractal formalism revisited with wavelets. Int. J. Bifurc. Chaos 4, 245–302. 10.1142/S0218127494000204 Muzy J. F. Delour J. Bacry E. (2000). Modelling fluctuations of financial time series: from cascade process to stochastic volatility model. Eur. Phys. J. B 17, 537–548. 10.1007/s100510070131 Muzy J.-F. Sornette D. Delour J. Arneodo A. (2001). Multifractal returns and hierarchical portfolio theory. Quant. Finance 1, 131–148. 10.1080/713665541 Nademanee K. McKenzie J. Kosar E. Schwab M. Sunsaneewitayakul B. Vasavakul T. . (2004). A new approach for catheter ablation of atrial fibrillation: mapping of the electrophysiologic substrate. J. Am. Coll. Cardiol. 43, 2044–2053. 10.1016/j.jacc.2003.12.05415172410 Nattel S. Burstein B. Dobrev D. (2008). Atrial remodeling and atrial fibrillation mechanisms and implications. Circ. Arrhythm. Electrophysiol. 62, 62–73. 10.1161/CIRCEP.107.754564 Nattel S. Harada M. (2014). Atrial remodeling and atrial fibrillation: recent advances and translational perspectives. J. Am. Coll. Cardiol. 63, 2335–2345. 10.1016/j.jacc.2014.02.55524613319 Nicolay S. Brodie of Brodie E. B. Touchon M. Audit B. d'Aubenton Carafa Y. Thermes C. . (2007). Bifractality of human DNA strand-asymmetry profiles results from transcription. Phys. Rev. E 75:032902. 10.1103/PhysRevE.75.03290217500744 Oboukhov A. M. (1962). Some specific features of atmospheric tubulence. J. Fluid Mech. 13, 77–81. 10.1017/S0022112062000506 Rensma P. L. Allessie M. A. Lammers W. J. Bonke F. I. Schalij M. J. (1988). Length of excitation wave and susceptibility to reentrant atrial arrhythmias in normal conscious dogs. Circ. Res. 62, 395–410. 10.1161/01.RES.62.2.3953338122 Richard C. D. Tanenbaum A. Audit B. Arneodo A. Khalil A. 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. 10.1016/j.jneumeth.2014.12.01625549550 Roland T. Khalil A. Tanenbaum A. Berguiga L. Delichère P. Bonneviot L. . (2009). Revisiting the physical processes of vapodeposited thin gold films on chemically modified glass by atomic force and surface plasmon microscopies. Surf. Sci. 603, 3307–3320. 10.1016/j.susc.2009.09.021 Roux S. G. Arneodo A. Decoster N. (2000). A wavelet-based method for multifractal image analysis. III. Applications to high-resolution satellite images of cloud structure. Eur. Phys. J. 15, 765–786. 10.1007/s100510051180 Roux S. G. Venugopal V. Fienberg K. Arneodo A. Foufoula-Georgiou E. (2009). Evidence for inherent nonlinearity in temporal rainfall. Adv. Water Resourc. 32, 41–48. 10.1016/j.advwatres.2008.09.007 Roy D. Talajic M. Nattel S. Wyse D. G. Dorian P. Lee K. L. . (2008). Rhythm control versus rate control for atrial fibrillation and heart failure. N. Engl. J. Med. 358, 2667–2677. 10.1056/NEJMoa070878918565859 Severs N. J. Bruce A. F. Dupont E. Rothery S. (2008). Remodelling of gap junctions and connexin expression in diseased myocardium. Cardiovasc. Res. 80, 9–19. 10.1093/cvr/cvn13318519446 Smeets J. L. Allessie M. A. Lammers W. J. Bonke F. I. Hollen J. (1986). The wavelength of the cardiac impulse and reentrant arrhythmias in isolated rabbit atrium. The role of heart rate, autonomic transmitters, temperature, and potassium. Circ. Res. 58, 96–108. 10.1161/01.RES.58.1.963943157 Snow C. J. Goody M. Kelly M. W. Oster E. C. Jones R. Khalil A. . (2008). Time-lapse analysis and mathematical characterization elucidate novel mechanisms underlying muscle morphogenesis. PLoS Genet. 4:e1000219. 10.1371/journal.pgen.100021918833302 Stevenson W. G. Stevenson L. W. (1999). Atrial fibrillation in heart failure. N. Eng. J. Med. 341, 910–911. 10.1056/NEJM19990916341120910486424 Takigawa M. Takahashi A. Kuwahara T. Okubo K. Takahashi Y. Watari Y. . (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. 10.1161/CIRCEP.113.00047124610740 Tan A. Y. Chen P.-S. Chen L. S. Fishbein M. C. (2007). Autonomic nerves in pulmonary veins. Heart Rythm 4, S57–S60. 10.1016/j.hrthm.2006.12.01117336886 The CAST II investigators (1992). The cardiac arrhythmia suppression trial. N. Engl. J. Med. 327, 227–233. Turiel A. Yahia H. Pérez-Vicente C. J. (2008). Microcanonical multifractal formalism – a geometrical approach to multifractal systems: part I. Singularity analysis. J. Phys. A: Math. Theor. 41:015501. 10.1088/1751-8113/41/1/015501 Ulphani J. S. Arora R. Cain J. H. Villuendas R. Shen S. Gordon D. . (2007). The ligament of Marshall as a parasympathetic conduit. Am. J. Physiol. Heart Circ. Physiol. 293, H1629–H1635. 10.1152/ajpheart.00139.2007 van Marion D. M. Lanters E. A. Wiersma M. Allessie M. A. Brundel B. B. de Groot N. M. (2015). Diagnosis and therapy of atrial fibrillation: the past, the present and the future. J. Atr. Fibrillation 8:1216. 10.4022/jafib.121627957185 Venugopal V. Roux S. G. Foufoula-Georgiou E. Arneodo A. (2006). Revisiting multifractality of high-resolution temporal rainfall using a wavelet-based formalism. Water Resour. Res. 42:W06D14. 10.1029/2005WR004489 Verma A. Jiang C.-Y. Betts T. R. Chen J. Deisenhofer I. Mantovan R. . (2015). Approaches to catheter ablation for persistent atrial fibrillation. N. Engl. J. Med. 372, 1812–1822. 10.1056/NEJMoa140828825946280 Vicsek T. (1989). Fractal Growth Phenomena. Singapore: World Scientific. Wang T. J. Larson M. G. Levy D. Vasan R. S. Leip E. P. Wolf P. A. . (2003). Temporal relations of atrial fibrillation and congestive heart failure and their joint influence on mortality. Circulation 107, 2920–2925. 10.1161/01.CIR.0000072767.89944.6E12771006 Wendt H. Abry P. Jaffard S. (2007). Bootstrap for empirical multifractal analysis. IEEE Signal Process. Mag. 24, 38–48. 10.1109/MSP.2007.4286563 West B. J. Shlesinger M. F. (1988). The noise in natural phenomena. American Scientist 78, 40–45. Wijffels M. C. Kirchhof C. J. Dorland R. Allessie M. A. (1995). Atrial fibrillation begets atrial fibrillation. A study in awake chronically instrumented goats. Circulation 92, 1954–1968. 10.1161/01.CIR.92.7.19547671380 Wolf P. A. Abbott R. D. Kannel W. B. (1991). Atrial fibrillation as an independent risk factor for stroke: the Framingham Study. Stroke 22, 983–988. 10.1161/01.STR.22.8.9831866765 Wolf P. A. Dawber T. R. Thomas H. E. Kannel W. B. (1978). Epidemiologic assessment of chronic atrial fibrillation and risk of stroke: the Framingham Study. Neurology 28, 973–973. 10.1212/WNL.28.10.973570666 Wynn G. J. El-Kadri M. Haq I. Das M. Modi S. Snowdon R. . (2016). Long-term outcomes after ablation of persistent atrial fibrillation: an observational study over 6 years. Open Heart 3:e000394. 10.1136/openhrt-2015-00039427547426 Yue L. Feng J. Gaspo R. Li G.-R. Wang Z. Nattel S. (1997). Ionic remodeling underlying action potential changes in a canine model of atrial fibrillation. Circ. Res. 81, 512–525. 10.1161/01.RES.81.4.5129314832 Zipes D. P. Jalife J. Stevenson W. G. (2017). Cardiac Electrophysiology: From Cell to Bedside. Philadelphia, PA: Elsevier Saunders.

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 No 2014-1R60212-00003295), by the President of Russian for the young scientists (Grant No 14.W01.17.2674-MK), and by the Metchnikov Program (EG-C visit to LOMA).