Next Article in Journal
Do Transitive Preferences Always Result in Indifferent Divisions?
Next Article in Special Issue
Entropy Measures in the Assessment of Heart Rate Variability in Patients with Cardiodepressive Vasovagal Syncope
Previous Article in Journal
Instantaneous 3D EEG Signal Analysis Based on Empirical Mode Decomposition and the Hilbert–Huang Transform Applied to Depth of Anaesthesia
Previous Article in Special Issue
Analyses of Heart Rate, Respiration and Cardiorespiratory Coupling in Patients with Schizophrenia
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Entropy Rate Maps of Complex Excitable Dynamics in Cardiac Monolayers

1
Max Planck Institute for Dynamics and Self-Organization, Am Faßberg 17, 37077 Göttingen, Germany
2
Institute for Nonlinear Dynamics, Georg-August-Universität Göttingen, Am Faßberg 17,37077 Göttingen, Germany
*
Author to whom correspondence should be addressed.
Entropy 2015, 17(3), 950-967; https://doi.org/10.3390/e17030950
Submission received: 16 December 2014 / Revised: 12 February 2015 / Accepted: 13 February 2015 / Published: 26 February 2015
(This article belongs to the Special Issue Entropy and Cardiac Physics)

Abstract

:
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.

1. 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.

2. 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.

2.1. Dataset and Preprocessing

The experimental dataset consists of an image sequence of 512 × 512 with a sampling time of ts =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.

2.1.1. 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].

2.2. 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.

2.3. 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.

2.4. 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:
D φ d r = ± 2 π
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.

2.5. 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 Vx,y(t) containing the information of whether the system at pixel (x, y), x ∈ [1,…, Nx], y ∈ [1, …, Ny] was in an excited state or not at a given time t ∈ [1, …, Nt].
A natural measure of the complexity of the sequence of excitations is then defined as the entropy H x , y ( 1 ) of the binary sequence Vx,y(t):
H x , y ( 1 ) = p 0 log p 0 p 1 log p 1 , p 1 = 1 N t t = 1 N t V x , y ( t ) , p 0 = 1 p 1
where p1 and p0 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:
d x , y ( t ) = i = 0 k 1 V x , y ( t i τ ) 2 i
which leads to 2k possible combinations for each symbol. The entropy can then be computed as:
H x , y ( k ) = i = 0 2 k 1 p i log p i
where pi is the probability for a specific integer within the time series dx,y(t).
Note that while H x , y ( 1 ) is merely a quantification of the number of occurrences of excitations (including their durations), H x , y ( k ) also takes their temporal order into account.

2.6. Entropy Rate Maps

We now consider the differences between entropies H x , y ( k ):
h x , y ( k ) = H x , y ( k ) H x , y ( k 1 )
The entropy rate is then defined as the asymptotic limit h x , y = lim k h x , y ( k ) 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 hx,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 kmax = 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.

3. Simulations of Excitable Media

3.1. 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:
d u d t = k u ( 1 u ) ( u a ) u v + D u
d v d t = ( ε + μ 1 v μ 2 + u ) ( v k u [ u ( a + 1 ) ] )
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.

3.2. 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:
d u d t = 1 ε u ( 1 u ) ( u v + b a ) + D u
d v d t = u 3 v
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.

4. 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.

4.1. 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.

4.1.1. 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.

4.1.2. 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.

4.1.3. 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.

4.2. 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.

4.2.1. 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.

4.2.2. 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.

4.2.3. 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.

4.3. 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.

4.3.1. Dominant Frequency Maps

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.

4.3.2. Phase Singularities

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.

4.3.3. 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.

4.4. 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.

5. 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.

Acknowledgments

The research leading to the results has received funding from the European Community’s Seventh Framework Programme FP7/2007-2013 under Grant Agreement 17 No. HEALTH-F2-2009-241526, EUTrigTreat. We acknowledge support from the German Federal Ministry of Education and Research (BMBF) (project FKZ 031A147, GO-Bio), the German Research Foundation (DFG) (Collaborative 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.

Conflicts of Interest

The authors declare no conflicts of interest.

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 × 104 cells per square centimeter, within 10-mm diameter glass retaining rings. The monolayers were then incubated at 37 °C and 5% CO2 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.1. 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.

References

  1. Cherry, E.M.; Fenton, F.H. Visualization of spiral and scroll waves in simulated and experimental cardiac tissue. New J. Phys. 2008, 10, 125016. [Google Scholar] [CrossRef]
  2. Parlitz, U.; Schlemmer, A.; Luther, S. Synchronization patterns in transient spiral wave dynamics. Phys. Rev. E 2011, 83, 057201. [Google Scholar] [CrossRef]
  3. Gray, R.A.; Pertsov, A.M.; Jalife, J. Spatial and temporal organization during cardiac fibrillation. Nature 1998, 392, 75–78. [Google Scholar]
  4. Borek, B.; Shajahan, T.K.; Gabriels, J.; Hodge, A.; Glass, L.; Shrier, A. Pacemaker interactions induce reentrant wave dynamics in engineered cardiac culture. Chaos 2012, 22, 033132. [Google Scholar]
  5. Shajahan, T.; Krinski, V.; Knyazeva, S.; Luther, S. Eliminating pinned spiral waves in cardiac monolayer by far field pacing. Proceedings of the 2014 8th Conference of the European Study Group on Cardiovascular Oscillations (ESGCO), Trento, Italy, 25–28 May 2014; pp. 151–152.
  6. Press, W.H.; Teukolsky, S.A.; Vetterling, W.T.; Flannery, B.P. Numerical Recipes: The Art of Scientific Computing, 3rd ed; Cambridge University Press: New York, NY, USA, 2007. [Google Scholar]
  7. Pandit, S.V.; Jalife, J. Rotors and the Dynamics of Cardiac Fibrillation. Circ. Res. 2013, 112, 849–862. [Google Scholar]
  8. Herlin, A.; Jacquemet, V. Reconstruction of phase maps from the configuration of phase singularities in two-dimensional manifolds. Phys. Rev. E 2012, 85, 051916. [Google Scholar]
  9. Iyer, A.N.; Gray, R.A. An Experimentalist’s Approach to Accurate Localization of Phase Singularities during Reentry. Ann. Biomed. Eng. 2001, 29, 47–59. [Google Scholar]
  10. Rogers, J. Combined phase singularity and wavefront analysis for optical maps of ventricular fibrillation. IEEE Bio-Med. Eng. 2004, 51, 56–65. [Google Scholar]
  11. Olkkonen, H. Computation of hilbert transform via discrete cosine transform. J. Signal Inf. Process. 2010, 1, 18–23. [Google Scholar]
  12. Schürmann, T.; Grassberger, P. Entropy estimation of symbol sequences. Chaos 1996, 6, 414–427. [Google Scholar]
  13. Faure, P.; Lesne, A. Estimating Kolmogorov Entropy from Recurrence Plots. In Recurrence Quantification Analysis; Webber, C.L., Jr., Marwan, N., Eds.; Understanding Complex Systems; Springer: Berlin/Heidelberg, Germany, 2015; pp. 45–63. [Google Scholar]
  14. pyEntropyRateMap A Python Algorithm for Computing Entropy Rate Maps. Available online: http://www.bmp.ds.mpg.de/software.html accessed on 15 February 2015.
  15. Jung, P. Coherent structure analysis of spatiotemporal chaos. Phys. Rev. E 2000, 61, 2095–2098. [Google Scholar]
  16. Bub, G.; Shrier, A.; Glass, L. Global organization of dynamics in oscillatory heterogeneous excitable media. Phys. Rev. Lett. 2005, 94, 028105. [Google Scholar]
  17. Aliev, R.R.; Panfilov, A.V. A Simple Two-variable Model of Cardiac Excitation. Chaos Solitons Fractals 1996, 7, 293–301. [Google Scholar]
  18. Fenton, F.H.; Cherry, E.M.; Karma, A.; Rappel, W.J. Modeling wave propagation in realistic heart geometries using the phase-field method. Chaos 2005, 15, 013502. [Google Scholar]
  19. Barkley, D. A model for fast computer simulation of waves in excitable media. Physica D 1991, 49, 61–70. [Google Scholar]
  20. Barkley, D. Barkley model. Scholarpedia 2008, 3. [Google Scholar] [CrossRef]
  21. MultiRecorder, a Multiple-Camera Recording and Analysis Software. Available online: http://www.bmp.ds.mpg.de/multirecorder.html accessed on 6 October 2014.
Figure 1. (a) Snapshots of fluorescence intensity; (b) sketch of the monolayer within the conductor generating the global electric field used for stimulation; (c) raw (unfiltered) time series at pixel location (64, 64). The frames are sampled with Δts = 18.29 ms. The cell culture has a diameter of 1 cm.
Figure 1. (a) Snapshots of fluorescence intensity; (b) sketch of the monolayer within the conductor generating the global electric field used for stimulation; (c) raw (unfiltered) time series at pixel location (64, 64). The frames are sampled with Δts = 18.29 ms. The cell culture has a diameter of 1 cm.
Entropy 17 00950f1
Figure 2. (a) Snapshots of the filtered movie (see the supplementary online material) and (b) section of the intensity time series at pixel location (64, 64) after the application of wavelet filtering and normalization (compare Figure 1).
Figure 2. (a) Snapshots of the filtered movie (see the supplementary online material) and (b) section of the intensity time series at pixel location (64, 64) after the application of wavelet filtering and normalization (compare Figure 1).
Entropy 17 00950f2
Figure 3. (a) Snapshots and (b) time series at location (64, 64) extracted from the simulated data using the Aliev–Panfilov model (Equations (6) and (7)). See the text for a description of the dynamics.
Figure 3. (a) Snapshots and (b) time series at location (64, 64) extracted from the simulated data using the Aliev–Panfilov model (Equations (6) and (7)). See the text for a description of the dynamics.
Entropy 17 00950f3
Figure 4. (a) Dominant frequency maps obtained for different windows of the simulated excitable medium using the Aliev–Panfilov model (Equations (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.
Figure 4. (a) Dominant frequency maps obtained for different windows of the simulated excitable medium using the Aliev–Panfilov model (Equations (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.
Entropy 17 00950f4
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)).
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)).
Entropy 17 00950f5
Figure 6. (a) Snapshots and (b) time series at location (64, 64) extracted from the simulated data using the cubic Barkley model (Equations (8) and (9)). See the text for a description of the dynamics.
Figure 6. (a) Snapshots and (b) time series at location (64, 64) extracted from the simulated data using the cubic Barkley model (Equations (8) and (9)). See the text for a description of the dynamics.
Entropy 17 00950f6
Figure 7. (a) Dominant frequency maps and (b) number of phase singularities vs. frame number (time) for the cubic Barkley simulation (Equations (8) and (9)). The two crosses mark the positions where two time series are sampled, which are given in Figure 8b,c.
Figure 7. (a) Dominant frequency maps and (b) number of phase singularities vs. frame number (time) for the cubic Barkley simulation (Equations (8) and (9)). The two crosses mark the positions where two time series are sampled, which are given in Figure 8b,c.
Entropy 17 00950f7
Figure 8. (a) Entropy rate maps for six different non-overlapping windows of a length of 500 frames for simulated data generated with the cubic Barkley simulation (Equations (8) and (9)). Two time series extracted from the fifth window at positions (28, 65) (b) and (45, 42) (c). Both positions are highlighted (white cross) in the corresponding dominant frequency map in Figure 7a and the entropy rate map in (a).
Figure 8. (a) Entropy rate maps for six different non-overlapping windows of a length of 500 frames for simulated data generated with the cubic Barkley simulation (Equations (8) and (9)). Two time series extracted from the fifth window at positions (28, 65) (b) and (45, 42) (c). Both positions are highlighted (white cross) in the corresponding dominant frequency map in Figure 7a and the entropy rate map in (a).
Entropy 17 00950f8
Figure 9. (a) Snapshots and (b) time series at location (64, 64) extracted from the cell culture data. See the text for a description.
Figure 9. (a) Snapshots and (b) time series at location (64, 64) extracted from the cell culture data. See the text for a description.
Entropy 17 00950f9
Figure 10. (a) Dominant frequency maps and (b) number of phase singularities vs. frame number (time) for the cell culture data. The two crosses mark the positions where two time series are sampled in Figure 11b,c. Both crosses are located in frequency regions with a dominant frequency of 1.2 HZ (note the small orange island around the upper left cross).
Figure 10. (a) Dominant frequency maps and (b) number of phase singularities vs. frame number (time) for the cell culture data. The two crosses mark the positions where two time series are sampled in Figure 11b,c. Both crosses are located in frequency regions with a dominant frequency of 1.2 HZ (note the small orange island around the upper left cross).
Entropy 17 00950f10
Figure 11. (a) Entropy rate maps for six different windows of a length of 500 frames without overlap for the cell culture data. Two time series extracted from the third window at positions (60, 90) (b) and (93, 63) (c). Both positions are highlighted (white crosses) in the corresponding dominant frequency map in Figure 10a and in the entropy rate map in subfigure (a).
Figure 11. (a) Entropy rate maps for six different windows of a length of 500 frames without overlap for the cell culture data. Two time series extracted from the third window at positions (60, 90) (b) and (93, 63) (c). Both positions are highlighted (white crosses) in the corresponding dominant frequency map in Figure 10a and in the entropy rate map in subfigure (a).
Entropy 17 00950f11

Share and Cite

MDPI and ACS Style

Schlemmer, A.; Berg, S.; Shajahan, T.K.; Luther, S.; Parlitz, U. Entropy Rate Maps of Complex Excitable Dynamics in Cardiac Monolayers. Entropy 2015, 17, 950-967. https://doi.org/10.3390/e17030950

AMA Style

Schlemmer A, Berg S, Shajahan TK, Luther S, Parlitz U. Entropy Rate Maps of Complex Excitable Dynamics in Cardiac Monolayers. Entropy. 2015; 17(3):950-967. https://doi.org/10.3390/e17030950

Chicago/Turabian Style

Schlemmer, Alexander, Sebastian Berg, T. K. Shajahan, Stefan Luther, and Ulrich Parlitz. 2015. "Entropy Rate Maps of Complex Excitable Dynamics in Cardiac Monolayers" Entropy 17, no. 3: 950-967. https://doi.org/10.3390/e17030950

APA Style

Schlemmer, A., Berg, S., Shajahan, T. K., Luther, S., & Parlitz, U. (2015). Entropy Rate Maps of Complex Excitable Dynamics in Cardiac Monolayers. Entropy, 17(3), 950-967. https://doi.org/10.3390/e17030950

Article Metrics

Back to TopTop