Causal Unit of Rotors in a Cardiac System

The heart exhibits complex systems behaviors during atrial fibrillation (AF), where the macroscopic collective behavior of the heart causes the microscopic behavior. However, the relationship between the downward causation and scale is nonlinear. We describe rotors in multiple spatiotemporal scales by generating a renormalization group from a numerical model of cardiac excitation, and evaluate the causal architecture of the system by quantifying causal emergence. Causal emergence is an information-theoretic metric that quantifies emergence or reduction between microscopic and macroscopic behaviors of a system by evaluating effective information at each spatiotemporal scale. We find that there is a spatiotemporal scale at which effective information peaks in the cardiac system with rotors. There is a positive correlation between the number of rotors and causal emergence up to the scale of peak causation. In conclusion, one can coarse-grain the cardiac system with rotors to identify a macroscopic scale at which the causal power reaches the maximum. This scale of peak causation should correspond to that of the AF driver, where networks of cardiomyocytes serve as the causal units. Those causal units, if identified, can be reasonable therapeutic targets of clinical intervention to cure AF.


I. INTRODUCTION
The heart is a multi-scale complex system consisting of five billion autonomous cardiomyocytes, each of which is a nonlinear dynamical system. However, it exhibits relatively simple system behaviors under physiologic conditions. During the regular heart rhythm, the macroscopic behavior of the heart is reducible to microscopic causal behavior and interactions of the cell population in the sinoatrial node ("supervenience"). In contrast, complex system behaviors emerge under pathologic conditions where the macroscopic collective behavior of the heart causes the microscopic behavior ("supersedence"). For example, as soon as the heart undergoes an order-disorder phase transition into fibrillation (Ashikaga and Asgari-Targhi, 2017), it controls the behaviors of individual cardiomyocytes to maintain itself.
This downward causation is clinically observable in a phenomenon of atrial fibrillation (AF) called "AF begets AF", where a longer duration of pacing-maintained AF results in a longer maintenance of AF after cessation of pacing (Wijffels et al., 1995).
The downward causation from macroscopic to microscopic behaviors of the cardiac system is quantifiable as inter-scale information flow that can be used as a surrogate for the mechanism that maintains AF ("AF driver"). In our previous work (Ashikaga and James, 2017b), we demonstate that transfer entropy accurately quantifies the upward and downward information flow between microscopic and macroscopic descriptions of the cardiac system with one of the potential AF drivers, a rotor, the rotation center of spiral waves (Narayan et al., 2012;Haissaguerre et al., 2014;Mandapati et al., 2000). We have also found that the downward information flow significantly decreases as the description of the system becomes more macroscopic. This subtle but important finding indicates that the relationship between the downward causation and scale is nonlinear. It is possible that, as the system is coarse-grained, it reaches a macroscopic scale at which the causal power peaks. Further coarse-graining removes the fine details of the causal architecture.
We hypothesize that the cardiac system with rotors has a scale where causal power reaches the maximum. To test the hypothesis, we describe rotors in multiple spatiotemporal scales by generating a renormalization group from a numerical model of cardiac excitation, and evaluate the causal architecture of the system by quantifying causal emergence (CE). CE is an information-theoretic metric that quantifies emergence or reduction between microscopic and macroscopic behaviors of a system by evaluating effective information (Hoel, Albantakis, and Tononi, 2013) at each spatiotemporal scale. Effective information (EI) is a quantity that captures causal interactions of a system between its unconstrained repertoire of possible cause and a specific state of possible effect (Tononi and Sporns, 2003).

II. METHODS
We perform the simulation and the data analysis using Matlab R2016b (Mathworks, Inc.).

A. Model of spiral waves
We use a monodomain reaction-diffusion model that was originally derived by Fitzhugh (FitzHugh, 1961) and Nagumo (Nagumo, Arimoto, and Yoshizawa, 1962) as a simplification of the biophysically based Hodgkin-Huxley equations describing current carrying properties of nerve membranes (Hodgkin and Huxley, 1952), which was later modified by Rogers and Mc-Culloch (Rogers and McCulloch, 1994) to represent cardiac action potential. This model accurately reproduces several important properties of cardiac systems, including slowed conduction velocity, unidirectional block owing to wavefront curvature, and spiral waves.
Here, v is the transmembrane potential with a finite action potential duration (APD), r is the recovery variable, and I ex is the external current (Pertsov et al., 1993). D is the diffusion tensor, which is a diagonal matrix whose diagonal and off-diagonal elements are 1 mm 2 /msec and 0 mm 2 /msec, respectively, to represent a 2-D isotropic system (Rogers and McCulloch, 1994). We solved the model equations using a finite difference method for spatial derivatives and explicit Euler integration for time derivatives assuming Neumann boundary conditions. We generate 1,000 sets of a 2-D 120 × 120 isotropic lattice of components (= 11.9 cm × 11.9 cm) by inducing spiral waves with 40 random sequential point stimulations in 40 random components of the lattice (Supporting Movie 1 ) (Ashikaga and James, 2017a).
In each component, we computed the time series for 10 seconds excluding the stimulation period with a time step of 0.063 msec, which was subsequently downsampled at a sampling frequency of 400 Hz.
We then defined the instantaneous phase φ(t) of the time series s(t) in each component via construction of the analytic signal ξ(t), which is a complex function of time.
Here the function s H (t) is the Hilbert transform of s(t) where p.v. indicates that the integral is taken in the sense of the Cauchy principal value.
We defined the rotor of the spiral wave as a phase singularity (Winfree, 1987), where the phase is undefined because all phase values converge. The phase singularity can be localized through calculation of the topological charge n t (Goryachev and Kapral, 1996;Mermin, 1979).
where φ( r) is the local phase, and the line integral is taken over the path l on a closed curve c surrounding the singularity (Bray and Wikswo, 2002).
In this study |n t | was used to quantify the average number of rotors over the entire time series.

B. Renormalization group
We generate a renormalization group of the system by a series of spatial and temporal transformation including coarse-graining and rescaling of the original microscopic description of the system. For each component, the time series of cardiac excitation is descretized to 1 when excited (during the APD at 90% repolarization, or APD 90 ) or 0 when resting ( Fig.1A) (Ashikaga et al., 2015). Then we coarse-grain the system spatially and temporally by decimation by a factor of 2 ( Figure 1B). Spatial decimation transforms a n × n lattice into a n 2 × n 2 lattice by extracting the top left component of each 2 × 2 block (Supporting Movie 2 ). Temporal decimation downsamples the binary time series of each component by a factor of 2. Using a combination of iterative coarse-graining in spatial and temporal axes we create a renormalization group of a total of 36 spatiotemporal scales of the system. The renormalization group includes spatial scales 1 (30×30 lattice), 2 (15×15 lattice), 3 (8×8 lattice), 4 (4×4 lattice), 5 (2×2 lattice), and 6 (1×1 lattice) ( Figure 1C), and temporal scales 1 (400 Hz), 2 (200 Hz), 3 (100 Hz), 4 (50 Hz), 5 (25 Hz) and 6 (12 Hz) (Fig.1D).

C. Effective information
We treat each component on the lattice as a time-series process X. Entropy H of each time-series process X is where p(x) denotes the probability density function of the time series generated by X.
Effective information quantifies the information generated when the system enters a specific state of possible effect Y out of its unconstrained probability distribution of possible cause X (Tononi and Sporns, 2003).
where X has a uniform probability distribution so that it provides the maximum entropy H(X) max . I(X; Y ) is mutual information, p(x, y) and H(X, Y ) denote the joint probability density function and the joint entropy of X and Y , respectively. Mutual information is originally a measure of statistical dependence to quantify how much information is shared between a source and a destination (Shannon, 1948). In this context, however, mutual information is applied between two time series of a system that is first perturbed into all possible states with equal probability and then observed as a sepcific state. Because of the system perturbations, mutual information here is a causal measure, and thus effective information of the system is a state-independent informational measure of a systems causal architecture (Hoel, Albantakis, and Tononi, 2013).
One can describe a n × n lattice at time t as a binary string of length n × n. Therefore, the unconstrained repertoire of all possible causes X at time t 0 consists of 2 n 2 possible states with equal probability 1/2 n 2 at each time point. We define the bin number b (b < 2 n 2 ) to calculate the probability distribution of X and Y , and we use b = 2 10 = 1, 024 in this study.
Analytically, because X has a uniform probability distribution, the probability that X falls in one of the b bins at each time point is 1/b. Therefore, entropy of X is equal to the maximum entropy ( Figure 2A).
Numerically, X can be defined as a vector of uniformly distributed random numbers between 1 and 2 n 2 -1 for a time series of finite duration. Due to the discretization effect, the probability is non-uniform. Entropy is close to but not identical to the maximum entropy ( Fig.2A).
Similarly, Y can be defined as a vector of decimal numbers between 1 and 2 n 2 -1, each of which represents a specific state of the system with rotors (Fig.2B).
Causal emergence is a difference in effective information between scales.
where m and n are different scales of the system description from the renormalization group.
When scale m is more macroscopic than scale n(m > n), a positive CE indicates that the macroscopic behavior is emergence (supersedence), whereas a negative CE indicates that the macroscopic behavior is reduction (supervenience) (Hoel, Albantakis, and Tononi, 2013).
In this study we quantify causal emergence with respect to the most microscopic system description with spatial scale = temporal scale = 1.

A. Evaluation of variance of effective information to quantify rotor dynamics
First we evaluate the variance of effective information in rotor dynamics. We repeat 1,000 numerical computations of X and Y in a representative spiral wave data set to calculate entropy H(X), H(Y ), H(X, Y ), then calculate EI(X → Y ). Numerically, H(X) is not uniquely determined due to the discretization effect, but the variance is small (Fig.3). spatial coarse-graining has minimal impact on the probability ditribution of H(X) from scales 1 through 4, but H(X) steeply falls in scales 5 and 6. In contrast, temporal coarse-graining gradually shifts the distribution of H(X) to the left. H(Y ) is uniquely determined because the probability of all bins is uniformly 1/b (shown in blue), and thus entropy is equal to the maximum entropy at log 2 b=10 bits. In contrast, numerically, the probability is non-uniform due to the discretization effect (shown in red). Entropy is 9.829 bits, which close to but not identical to the maximum entropy. B. Probability distribution of a specific state of possible effect Y . The probability is non-uniform. Entropy is 2.289 bits in this case. C. Bivariate probability distribution of cause X and effect Y . Joint entropy is 10.220 bits in this case. Effective information from case X to effect Y is equal to mutual information between X and Y , thus is calculated as 1.898 bits.
it represents a specific state of the system rgardless of the spatiotemporal scale (Fig.4).
In this case, spatial coarse-graining clearly increases the distribution of H(Y ) to the right, which peaks at scale 4 and decreases at scales 5 and 6. Similarly, temporal coarse-graining increases the distribution of H(Y ) to the right, which peaks at scale 4 and decreases at scales 5 and 6. The relationship between the spatiotemporal coarse-graining and the probability distribution of joint entropy H(X, Y ) is similar to that of H(X) (Fig.5), and the variance remains small. Effective infromation EI(X → Y ) peaks at spatial scale of 4 and temporal scale 5, and the variance of EI(X → Y ) remains small (Fig.6). This result indicates that, despite the discretization effect, numerical computation of EI(X → Y ) is robust with a high reproducibility, and thus EI(X → Y ) can be used to quantify the information of rotor dynamics at each spatiotemporal scale.

B. Evaluation of effective information in aggregate data sets
Next, we quantify effective information of the renormalization group of a total of 36 spatiotemporal scales of the system in aggregate data of 1,000 sets (Fig.7). Overall, effective information increases as the scale increases from microscopic to macroscopic descriptions of the system. However, it reaches the global maximum at spatial scale = temporal scale = 4, beyond which effective information decreases (Fig.7). The difference in effective information between scales is larger in spatial coarse-graining (Fig.7B) than that of temporal coarse-graining (Fig.7C), indicating that the impact of spatial coarse-graining on effective information is higher than that of temporal coarse-graining. is not uniquely determined due to the discretization effect, but the variance is small.
Each subplot represents the probability distribution of H(X, Y ). The columns represent the spatial scales (1 through 6) and the rows represent the temporal scales (1 through 6).

C. Relationship between the number of rotors and causal emergence
Lastly, we evaluate the relationship between the number of rotors and causal emergence in individual data sets. The number of rotors ranges from 0 to 7, with a median of 3 ( Fig.8). For system descriptions at spatial scale≤4 and temporal scale≤4, causal emergence is positive for all the data sets except a few where a rotor prematurely disappears on its own (number of rotors≤1, red dots in Fig.9). There is a significant positive correlation between the number of rotors and causal emergence, indicating that causal emergence consistently increases as the number of rotors increases. For system descriptions at spatial scale≥5, causal emergence is negative for all the data sets, and there is a significant negative correlation between the number of rotors and causal emergence. This indicates that the macroscopic behavior at those scales are reducible to the microscopic behavior. For system descriptions at spatial scale=1 and temporal scale≥5, causal emergence scatters in positive and negative values. This indicates that the causal relationship at those scales is inconsistent. There is still a significant negative correlation between the number of rotors and causal emergence, but the correlation coefficients are small. For system descriptions at spatial scale=2,3,4 and temporal scale≥5, causal emergence is almost always positive and there is a significant positive correlation between the number of rotors and causal emergence. This result indicates that temporal coarse-graining has a smaller impact than spatial coarse-graining on the causal architecture.

A. Main findings
First, we find that numerical computation of effective information in the cardiac system with rotors is robust with high reproducibility despite the discretization effect. Therefore effective information can be used to quantify the information of rotor dynamics at each spatiotemporal scale.
Next, we find that there is a spatiotemporal scale at which effective information peaks in the cardiac system with rotors. This suggests the presence of causal units at this scale consisting of networks of components that serve as the AF driver.
Lastly, we find that there is a positive correlation between the number of rotors and causal emergence up to the scale of peak causation. This indicates that the causal relationship between the macroscopic and the microscopic behaviors reverses beyond that scale.

B. Causal architecture of rotors
AF currently impacts the lives of 33 million patients worldwide (Chugh et al., 2013;Rahman, Kwan, and Benjamin, 2014). Importantly, AF is associated with a five-fold increased  risk of thromboembolism such as stroke (Wolf, Abbott, and Kannel, 1991), and accounts for 15% of strokes overall (Investigators et al., 1994). In addition, AF increases the risk of cognitive impairment (Kalantarian et al., 2013;Thacker et al., 2013) independent of clinical stroke. AF is also associated with a 2-fold increased risk of dementia (Ott et al., 1997), and more than 10% of AF patients develop dementia over 5 years (Miyasaka et al., 2007).
Furthermore, AF is a powerful risk factor of myocardial infarction (Soliman et al., 2014) and death (Benjamin et al., 1998). Although the mechanism that initiates AF is ascribed to focal triggers primarily from the pulmonary veins (Haissaguerre et al., 1998), the AF driver remains unknown.
In this study we describe rotors in multiple spatiotemporal scales by generating a renormalization group of the cardiac system and evaluate the causal architecture of the system by quantifying causal emergence. Causal emergence was originally developed in neuroscience but is applicable to any multi-scale systems (Hoel, 2017). Our analysis using causal emergence confirms that rotors are emergent behaviors of the heart, that is, macroscopic collective behaviors that cause microscopic behaviors. This indicates that a multi-scale, complex systems approach is an appropriate direction of investigation to understand the AF driver, rather than the reductionistic approach to understanding the AF driver by describing microscopic behaviors of the system with near-infinite degrees of freedom.
In this particular cardiac system, we find that effective information peaks at spatial scale = temporal scale = 4. The spatial scale of 4 divides the original 11.9cm × 11.9cm lattice of cardiac system into 4×4 units, each of which occupies 3cm×3cm in size. This 2-D unit contains 5-7.5×10 5 cells with dimensions of human atrial cardiomyocytes (e.g., 120 µm length and 10-15µm diameter (Nygren et al., 1998)). It is likely that different cardiac systems and patients have different sizes of causal unit of rotors just as the "critical mass" needed to sustain fibrillation (McWilliam, 1887;Garrey, 1914) is different for different patients ("effective size" (Panfilov, 2006)). However, the important point is that one can coarsegrain the cardiac system with rotors to identify a macroscopic scale at which the causal power reaches the maximum. This scale of peak causation should correspond to that of the AF driver, where networks of cardiomyocytes serve as the causal units.

C. Limitations
We used a modified Fitzhugh-Nagumo model, which is a relatively simple model of excitable media, with a homogeneous and isotropic lattice. It is possible that our findings may not directly be extrapolated to a more realistic cardiac system with tissue heterogeneity and anisotropy. However, the information-theoretic metrics in this study are independent of any specific trajectory of each rotor. Therefore, our approach is applicable to any other cardiac system.

D. Conclusions
One can coarse-grain the cardiac system with rotors to identify a macroscopic scale at which the causal power reaches the maximum. This scale of peak causation should correspond to that of the AF driver, where networks of cardiomyocytes serve as the causal units.
Those causal units, if identified, can be reasonable therapeutic targets of clinical intervention to cure AF.

FUNDING
This work was supported by the Fondation Leducq Transatlantic Network of Excellence.

SUPPLEMENTARY MATERIAL
Supporting Movie 1. Random sequential point stimulations. We induce spiral waves by introducing 40 random sequential point stimulations in 40 random components of the lattice.
In this example, random sequential point stimulations induce five spiral waves.
Supporting Movie 2. Renormalzation group. The movie shows a renormalization group of the cardiac system with two spiral waves by a series of transformation including coarsegraining and length rescaling (scale 1 through 6). For each component, the time series of cardiac excitation is descretized to 1 (black) when excited (during the APD at 90% repolarization, or APD 90 ) or 0 (white) when resting.