Entropy Rate Maps of Complex Excitable Dynamics in Cardiac Monolayers

The characterization of spatiotemporal complexity remains a challenging task. This holds in particular for the analysis of data from fluorescence imaging (optical mapping), which allows for the measurement of membrane potential and intracellular calcium at high spatial and temporal resolutions and, therefore, allows for an investigation of cardiac dynamics. Dominant frequency maps and the analysis of phase singularities are frequently used for this type of excitable media. These methods address some important aspects of cardiac dynamics; however, they only consider very specific properties of excitable media. To extend the scope of the analysis, we present a measure based on entropy rates for determining spatiotemporal complexity patterns of excitable media. Simulated data generated by the Aliev–Panfilov model and the cubic Barkley model are used to validate this method. Then, we apply it to optical mapping data from monolayers of cardiac cells from chicken embryos and compare our findings with dominant frequency maps and the analysis of phase singularities. The studies indicate that entropy rate maps provide additional information about local complexity, the origins of wave breakup and the development of patterns governing unstable wave propagation. Entropy 2015, 17 951


Introduction
Quantifying spatiotemporal complexity is a frequently occurring task for image sequences obtained in optical mapping experiments in excitable media, such as cardiac excitation waves [1]. For example, electrical excitation waves in cardiac tissue provide important information about structural features and distinct dynamical states related to normal heart beat or arrhythmias. Suitable spatiotemporal complexity measures are necessary for improving classification and forecasting techniques and may also reveal information about the underlying structure of the observed medium (although this inference may be prone to artifacts and has to be carried out with care [2]). In the context of excitable cardiac dynamics, there are some well-established analysis procedures employing dominant frequency maps and analysis of phase singularities [3]. These concepts primarily focus on periodic dynamics. Arrhythmias, however, correspond to dynamically complex states, and therefore, alternative approaches for describing the spatiotemporal dynamics are required. In this article, we develop a method based on entropy rates of binary symbolic strings with the aim of estimating the spatiotemporal complexity in data obtained in optical mapping experiments using monolayers of cardiomyocytes obtained from embryonic chicken. We validate the entropy rate approach using synthetic data generated data with the Aliev-Panfilov model and the cubic Barkley model. Then, we apply it to experimental data from embryonic chicken cell cultures and compare the results with dominant frequency maps and the analysis of phase singularities. All movies that we refer to in this paper are available as supplementary online material.

Experimental Setup
Cultured cardiac cell monolayers are frequently used to investigate the propagation of excitation waves. Despite structural and functional differences from native tissue, the mechanisms underlying the normal function of the heart, as well as arrhythmic behavior can be studied in a simplified two-dimensional multicellular substrate. The dataset discussed in the following section has been chosen because it displays different types of dynamics often visible in the dynamics of cardiac tissue, and therefore, seems to be appropriate for method development, benchmarking and comparison.
The experimental protocol of the cell culture experiments has been carried out as described previously in [4,5]. For convenience, the corresponding texts have been slightly adapted and are given in the Appendix A.1 and A.2.

Dataset and Preprocessing
The experimental dataset consists of an image sequence of 512 × 512 with a sampling time of t s = 18.29 ms. During preprocessing, these frames are 4 × 4 binned, resulting in 128 × 128 images with a spatial resolution of ∆x = ∆y = 81.3 µm per pixel. Figure 1a shows six snapshots of the cell culture movie (see the supplementary online material). The first part of the movie data shows a slightly moving, rotating spiral centered in the upper half of the field of view. Starting from Frame 595 (approximately 11 s), a series of monophasic electric field pulses of amplitude 2 V/cm and a duration of 20 ms were applied at a frequency of 0.8 Hz from two electrodes, as shown in Figure 1b, leading to wave emission within the medium and, thus, establishing a more complex type of dynamics. The remainder of the movie, which has been analyzed here up to Frame 3000, shows different regions of interacting spiral waves and wave breakup. Visual inspection suggests increasing stability of the medium as soon as the series of pulses ends at Frame 2739 (approximately 50 s). Figure 1c shows a section of a typical unfiltered fluorescence intensity time series at pixel (64, 64). The raw video data has a poor signal-to-noise ratio, which makes detection of the action potential very difficult. As the explanatory power of the following entropy calculations will be dependent on a precise detection of the excitation peaks, we need to preprocess the data and apply filters that decrease the noise.

Signal Processing
Gaussian smoothing can be applied by convolving the data with a Gaussian kernel. This can be done temporally, spatially or spatiotemporally with the respective 1D-, 2D-or 3D-kernels. In this way, any pixel value in the dataset is replaced by the weighted average of its neighboring pixels. Depending on the standard deviation used for the Gaussian kernel, different filtering strengths can be implemented. Although this method is commonly used for imaging data, it has the disadvantage of blurring out the shape of the action potential. We therefore selected wavelet filtering, which is briefly presented in the next paragraph, in order to filter our video data.
The wavelet transform [6] provides a straightforward method of filtering, as it decomposes the data signal into specific representations at different scales. This decomposition is carried out using a convolution of the signal with window functions called wavelets. Wavelets are localized in time or space, can be moved and scaled and form an orthogonal (and in some cases, orthonormal) basis. The method of wavelet filtering consists of wavelet transforming the original data, setting the components at undesired (here, fine) scales to zero and performing an inverse transformation. A particular strength of the wavelet transform is that the shape of the wavelet can be tuned to the shape of the signal to optimize the filtering performance. Similar to the Fourier transform, there exists a fast algorithm for computing the wavelet transform, called fast wavelet transform. The wavelet filtering method can be applied to arbitrarily many dimensions. We use one-and two-dimensional transforms here separately in order to be able to use different wavelets in each dimension. The mother wavelets were chosen to resemble the shape of the action potential. Daubechies wavelets of order four are employed for spatial filtering and a Daubechies wavelet of order three for temporal filtering [6].

Normalization
The amplitude of the fluorescence intensity signal can vary with respect to its spatial position in the medium, which imposes difficulties for peak detection. This inhomogeneity is due to the inhomogeneity of the medium. In order to cope with this difficulty, our data is normalized after the filtering procedure by subtracting the temporal mean and by dividing by the standard deviation. As this would amplify very noisy regions, we threshold the signal, setting all time series to zero where the (temporal) standard deviation does not exceed a threshold value. This has the advantage that the data are masked automatically, so that only the region of interest, containing the actual medium, remains visible. To obtain a smooth border, the regions that are of no interest are masked out afterwards using a circular mask. Figure 2 shows the result of our preprocessing procedure, revealing clearly visible waves and action potentials.  Figure 1).

Dominant Frequency Maps
Dominant frequency maps are a widely used tool to investigate the spatial structure of (temporal) frequency components within an excitable medium [7]. At each pixel of the movie (see the supplementary online material), the power spectrum is computed from the corresponding time series after subtracting the mean. The frequency of the Fourier component with the highest amplitude is calculated from the spectrum. This frequency value is plotted according to its spatial position, which leads to a color-coded image relating the dominant frequency components to their spatial locations. This procedure can be applied separately for specific temporal windows, which allow investigating the temporal evolution of dominant frequency maps.

Phase Singularity Analysis
A prominent feature visible in many excitable media is spiral wave dynamics [7]. A straight forward approach is therefore to look at the dynamics of spiral cores, which includes movement, pairwise appearance and annihilation [8]. Spiral core centers manifest as singularities in the phase map, and many efficient algorithms have been developed to identify and track them in their temporal evolution (e.g., [9,10]).
Phase maps, which display the instantaneous phase at each spatial location, can be efficiently computed from preprocessed signals using the Hilbert transform [11]. Knowing the phase as a function of space, we find the spiral cores by integrating the path integral of the gradient of the phase ∇ϕ over closed paths ∂D. The integral vanishes if the region D either contains no phase singularity or cores of an equal number of clockwise and counter clockwise spirals. If D contains only a single phase singularity, we obtain: where the sign indicates clockwise or counter clockwise rotation of the spiral wave [9]. Trajectories of phase singularities are reconstructed using nearest-neighbor tracking. Phase singularities that are only visible for less than five frames are considered artifacts and are removed.

Entropy of Binary Sequences
One of the most prominent features of excitable media is the presence of an excited state that is very clearly distinguishable from the resting state. As visible in the snapshots of the time series shown in Figure 2, these states can be easily differentiated in the normalized videos using a thresholding procedure for the individual time series. Thresholding is carried out by setting the signal to one whenever the signal amplitude is above the threshold and to zero elsewhere. We use a fixed threshold of 0.6 for all systems investigated in this study. This step leads to a binary video V x,y (t) containing the information of whether the system at pixel (x, y), x ∈ [1, . . . , N x ], y ∈ [1, . . . , N y ] was in an excited state or not at a given time t ∈ [1, . . . , N t ].
A natural measure of the complexity of the sequence of excitations is then defined as the entropy H (1) x,y of the binary sequence V x,y (t): where p 1 and p 0 denote the relative frequencies of pixel (x, y) being excited or not, respectively. This definition can be extended to symbol strings of length k using delay embedded binary sequences with delay τ. Each binary sequence can then be represented by an integer: which leads to 2 k possible combinations for each symbol. The entropy can then be computed as: where p i is the probability for a specific integer within the time series d x,y (t).
Note that while H x,y is merely a quantification of the number of occurrences of excitations (including their durations), H (k) x,y also takes their temporal order into account.

Entropy Rate Maps
We now consider the differences between entropies H (k) x,y : The entropy rate is then defined as the asymptotic limit h x,y = lim k→∞ h (k) x,y and quantifies the amount of entropy per symbol [12]. Because of its relation to the Kolmogorov-Sinai entropy [13], this quantity has an intuitive interpretation: the Kolmogorov-Sinai entropy takes a value of zero for periodic systems, a finite value larger than zero for chaotic systems and goes to infinity when there are no correlations in phase space, as is expected for very noisy systems. To obtain a spatially-resolved estimate for the entropy rate, we create maps of h x,y in the results section and call them entropy rate maps [14]. Since our data are of finite length, there is an important tradeoff to take into account: for larger k, the entropy rate estimates will become better, but the estimates of the individual symbol probabilities will become worse. For all our computations presented in the results section, we chose a maximum k max = 14. Furthermore, we found an embedding delay of τ = 4 frames to be a good value to capture all relevant time scales of our systems. As previously mentioned in Section 2.1.1, a binarization threshold of 0.6 has been used for all videos.
It is pointed out that the method presented here does not imply an explicit quantification of both spatial and temporal aspects of the dynamics, as the quantification is computed for each time series at each pixel individually. Presentation in the form of a map allows for an investigation of the spatial distribution of the dynamical properties within the extended medium, which could also be analyzed using, e.g., pattern recognition techniques.
Another measure for spatiotemporal complexity for excitable media has been proposed in [15], and its application to data obtained from embryonic chicken heart cell cultures has been demonstrated in [16]. This method uses the size distribution of spatiotemporal clusters of wavefronts to compute an entropy. However, it relies on a more fine-grained wavefront structure, which is not present in our case.

Aliev-Panfilov Model
For the first type of simulated data, the Aliev-Panfilov model [17] was used. This model is a simple two-variable model for the canine myocardium, defined by: The simulations were performed using the parameter values ε = 0.01, µ 2 = 0.3, k = 8 and D = 1.6 with a grid constant of ∆x = ∆y =1. The PDEs (Equation (7)) were solved with Euler stepping using a time step of 0.1s. After six time steps, a frame was generated and stored (i.e., the time interval between the frames that are used for analysis is 0.6s). No-flux boundary conditions were implemented using the phase field method [18]. The circular simulation region has a radius of 59 pixels.
To simulate a case with different dynamics in different regions, a linear ramp between Pixels 32 and 96 in the vertical direction was introduced at Frame 600, where the parameters a and µ 1 were varied between a = [0.06, 0.11] and µ 1 = [0.09, 0.11], where a = µ 1 = 0.11 leads to stable, non-chaotic dynamics. This linear ramp establishes a smooth transition between a stable region and a region with continuously increasing wave-breakup, which allows for an investigation of the detectability of spatially distributed features within the data. At Frame 1800 of the simulation, we set a = 0.02 and µ 1 = 0.09 in the whole domain to achieve another complex type of dynamics. This latter parameter set has been introduced to analyze the response of the analysis methods to another type of complexity, although this specific parameter set is considered artificial and not related to a situation found in real cardiac tissue.

Cubic Barkley Model
The Barkley model [19] is a simple model for a reaction diffusion system. As a second type of simulated data, we use a variation of it, where a cubic term is introduced in the second variable to achieve spatio-temporal chaos [20]. The model equations are: For the simulations, we used the parameter values a = 0.75, b = 0.0006 and D = 8. The parameter ε was varied between ε = [0.06, 0.08], where 0.6 gives stable spirals. The simulation was made using Euler stepping with a time step of 0.01s, and the time interval between frames is 0.7s. The circular simulation region has a radius of 59 pixels. Again a linear ramp was introduced from top to bottom on Frame 600 between Pixels 48 and 80 in the vertical direction. On Frame 1800, the parameter ε was set to ε = 0.08 in the whole medium.

Results and Discussion
In this section, we present results obtained for three types of analyses for three different datasets: • Simulated data using the Aliev-Panfilov model introduced in Section 3.1.
• Simulated data using the cubic Barkley model introduced in Section 3.2.
• The cell culture dataset described in Section 2.1.
For each dataset, a dominant frequency analysis, an analysis of the development of the number of phase singularities and an analysis of entropy rate maps have been carried out.

Aliev-Panfilov Model
The Aliev-Panfilov model simulation shows a stable spiral wave visible in the upper half of the movie (see the supplementary online material) until Frame 600. After the activation of the linear parameter variation (Frame 600), the waves continuously break up in the lower half of the video. From Frame 1800 onwards, an interaction of two different types of wave patterns leads to a very complex situation: a few very fast waves move through pipes formed by very slow waves and create reentrant dynamics. Note that, while the parameter ranges up to Frame 1800 were designed to reflect the dynamics observed in real cardiac tissue, this parameter range is solely used for the investigation of very complex dynamics and does not correspond to a realistic type of heart tissue. Figure 3 shows snapshots and a typical time series for this simulated dataset.   (6) and (7)). See the text for a description of the dynamics.

Dominant Frequency Maps
The dominant frequency maps, displayed in Figure 4a, allow for a quick identification of regions of interest with respect to the dominant frequency component in the underlying time series. The single spiral in the first window, as well as the wave breakup regions in Windows 2-4 seem to match the observed dynamics very well. A reduced frequency is observed in the wave breakup regions, due to the stretched duration of the action potential. This behavior is also visible in the time series shown in Figure 3b. In the last two windows, the dominant frequencies decrease even more and show a more spatially inhomogeneous pattern.   (6) and (7)). All windows have a length of 500 frames and are moved successively in steps of 500 frames. (b) The number of phase singularities vs. frame number (time) for simulated data using the Aliev-Panfilov model.

Phase Singularities
The result of the phase singularity statistics is displayed in Figure 4b. It clearly identifies the stable spiral, which is visible during the first 600 frames. The second phase with spiral breakup and non-periodic dynamics shows an increased number of phase singularities varying between one and seven. The complex wave pattern visible from Frame 1800 onwards manifests itself in an even higher variation of phase singularities reaching from 15 to no spirals at all.

Entropy Rate Maps
As described in Section 2.6, entropy rate maps can be thought of as an estimate of the local complexity of the dynamics of the medium. Values close to zero correspond to periodic signals, as the entropy production converges to zero. Chaotic signals correspond to high entropy production. Figure 5 shows the resulting entropy rate maps for the Aliev-Panfilov data and confirms the intuitive understanding of the entropy rate to some extent: the first three windows clearly differentiate the stable from the more complex regions, showing fine grained variations of local complexity. Note that in the first window, high entropy rate values are visible close to the lower boundary. These are due to the transient behavior due to initial conditions.
As the fourth window already contains the even-more complex behavior starting from Frame 1800, these distinct patterns are already visible in the corresponding entropy rate map, together with a clearly visible symmetry breakup.
The last two windows highlight the differences between the fast waves with very complex and unpredictable motions and the slow waves, which lead to a much lower entropy production within the analyzed time frames.  Figure 5. The entropy rate maps for six different non-overlapping windows of a length of 500 frames for simulated data generated with the Aliev-Panfilov model (Equations (6) and (7)).

Cubic Barkley Model
The video obtained from the cubic Barkley model simulations initially shows single spiral activity in the upper half of the video, resulting in plane waves traveling over the lower half of the video. Snapshots and a representative time series are shown in Figure 6. Starting from approximately Frame 800, the first visible changes of the introduced inhomogeneity are visible in the lower half of the simulated dish. This is followed by increasingly complex behavior due to wave breakup. Around Frame 2000, the simulated medium seems to exhibit even more complex behavior, showing many spiral waves.  (8) and (9)). See the text for a description of the dynamics.   (9)). The two crosses mark the positions where two time series are sampled, which are given in Figure 8b,c.

Dominant Frequency Maps
The dominant frequency maps show initially a quite homogeneous frequency distribution with a small area containing a frequency change in the second window. From the third window on, the number of different dominant frequency components increases.

Phase Singularities
As shown in Figure 7, the number of phase singularities stays constant at one, as expected, until approximately Frame 880. From that frame on, the number of phase singularities continuously increases. From Frame 2500 on, there seems to be a stabilization with a fluctuating number of 7 to 17 singularities, until the end of the dataset.

Entropy Rate Maps
The entropy rate maps, which are visible in Figure 8, show an initially overall stable picture with a more stable region close to the spiral core and a bit less stable region close to the lower and right borders of the simulated disk. In the second temporal window, the first signs of wave breakup are visible in the lower left part of the medium, which resembles the optical impression from the video very well. Windows 3 and 4 show the increasingly stronger wave breakup effect originating from the lower left part of the dataset. In Windows 5 and 6, a separation of small stable regions and larger unstable regions is visible, where the stable regions seem to correspond to meandering spiral cores and periodic wave emission sites.

Cell Culture Data
The dataset has been described in Section 2.1. Figure 9 shows snapshots of the cell culture dataset at different frames together with a complete time series extracted from the video data at a central location. The third (upper right) snapshot already shows activity due to the pulses applied (see Section 2.1).
The following snapshots show unordered wave dynamics involving several meandering and interacting spirals and continuous wave emission (due to pulsing), mainly from the right parts of the disk.   Figure 10a shows dominant frequency maps, which have been generated for non-overlapping windows of size 500 frames for the whole video. The first window shows a homogenous situation with a dominant frequency of 0.55 Hz. The small spot at (52, 99) with a frequency of 1.1 Hz marks the core of the spiral. Windows 2 to 5 show similar diagrams with frequencies between 0.5 Hz and 1.3 Hz. The higher frequencies are mainly located in the right and lower parts. In the last window the large region at frequency 1.2 Hz seems to disintegrate into smaller areas. Figure 10b shows the number of phase singularities over time for the cell culture data. It is easily visible that during the initial one-spiral phase, only one phase singularity is detected, succeeded by a steep increase of the number of singularities as soon as the pulsing starts. In the remainder of the data, the number of singularities shows similar behavior, varying between 6 and 31.

Entropy Rate Maps
The first window of the entropy rate maps visible in Figure 11 includes only one spiral; therefore, the map shows homogeneously low values. The second window, which contains the overlap between the stable spiral and the beginning of the pulses, allows one to distinguish between stable areas and more complex areas. The stable area to the lower right becomes even better separable from the more complex regions in the third window. The more complex regions show inhomogeneous patterns and seem to move over time, while the border to the stable region seems to stay almost constant. The final window of the video, which contains only a few pulses, results in the spreading of the complex inhomogeneous entropy rate patterns over the whole disk.

Comparison of Entropy Rate Maps and Dominant Frequency Maps
In order to obtain a better understanding of the different features of excitable media, which are highlighted in dominant frequency maps and entropy rate maps, respectively, we are now going to discuss two pairs of time series extracted from two specific windows of a length of 500 frames from our datasets.
The first pair for comparison has been selected from two points within the fifth window of the cubic Barkley simulation, where the dominant frequency clearly suggests regions of different dynamic behavior (Figures 7a and 8). While the shift in the dominant frequency from 0.17 Hz at (28, 65) to 0.20 Hz at (45, 42) defines the main structures within this dominant frequency map, the decrease in the entropy rate from 0.12 nat to 0.10 nat seems to be of less importance within the corresponding entropy rate map. This confirms the intuitive understanding of the entropy rate, which contains information about periodicity, regardless of the underlying frequency.
The second pair for comparison has been extracted from the third window of the cell culture data at points (60, 90) ( Figure 11b) and (93, 63) ( Figure 11c). In both cases, the dominant frequency equals 1.2 Hz, while the entropy rate decreases from 0.14 nat at (60, 90) to 0.09 nat at (93, 63). This behavior corresponds to the difference between the large blue area, which is considered more periodic in the entropy map, and the more inhomogeneous area, which is considered more complex. The reason for the low variability in the blue area is that the effect of the external pacing is strongest in this region. The time series in Figure 11b displays a higher variation in amplitude, while the visual impression of Figure 11c is a more periodic one. On the other hand, this difference is not highlighted by the dominant frequency map in Figure 10a.

Conclusions
Entropy rate maps provide a robust and efficient way to compute a spatially resolved measure of the local complexity in optical mapping and simulated data of cardiac cell culture. Specifically, it allows one to differentiate regions with spiral wave dynamics or other periodic wave propagation from regions with more complex/chaotic behavior. As was seen in simulations using the cubic Barkley model (Equations (8) and (9)), the origins of wave breakup can be spotted, and their spreading over the tissue can be investigated in detail. An important application might be the screening of the dynamics of large amounts of data, where a quick overview of the local complexity is often desired. Further investigations of this method are needed in order to address open questions, such as the optimal choice of the threshold parameter, the effect of different types of normalization and the quantification of the complexity in the spatial domain.
Research Centres SFB 1002 Project C3 and SFB 937 Project A18) and the German Center for Cardiovascular Research (DZHK e.V.).

Author Contributions
A. Schlemmer implemented algorithms and applied them to the experimental and numerical data. T.K. Shajahan did the cell culture experiments and performed the measurements. S. Berg created the simulations and produced the numerical data. U. Parlitz and S. Luther suggested the topic and wrote the manuscript, along with A. Schlemmer and S. Berg. All authors have read and approved the final manuscript.

A. Appendix
The following sections give some brief information about the experimental techniques used. More details can be found in [4,5].

A.1. Preparation of Cardiac Monolayers
Cardiac monolayers were prepared using ventricular myocytes from eight-day old chicken embryos, as described previously [4]. Ventricles were isolated, dissociated in trypsin (Sigma Aldrich, St. Louis, MO, USA), resuspended in growth medium 818A and then plated on fibronectin coated cover glass at confluent densities, 36 × 10 4 cells per square centimeter, within 10-mm diameter glass retaining rings. The monolayers were then incubated at 37 • C and 5% CO 2 for 48 h in the growth medium 818A. To map calcium activity, each of these was loaded with Ca-green I dye (2 µg/2 mL) in Hanks' Balanced Salt Solution (HBSS) containing 5 µL pluronic acid (Invitrogen, Waltham, MA, USA) for 30 min. They were then washed with HBSS and transferred to the imaging chamber, maintained at 33 • C.

A.2. Optical Mapping
Calcium activity in the monolayer was optically mapped using an Olympus MVX10 microscope and an EMCCD (Electron Multiplying Charge Coupled Device) camera (photometrics cascade II). About a square centimeter of the monolayer could be imaged with a spatial resolution of 256 × 256 (2 × 2 camera binning) pixels at a frequency of 55 frames per second, which is sufficient to capture electrical waves in this monolayer. Since the voltage signals from the monolayer do not have sufficient signal strength, we map intracellular calcium activity using the dye, calcium green (Invitrogen). Custom written software in Python and matplotlib are used for live imaging and data analysis [21].
Before further analysis, the data were again binned using a 2 × 2 binning, since the spatial resolution is much higher than needed for our analysis.

Conflicts of Interest
The authors declare no conflicts of interest.