Dynamical density delay maps: simple, new method for visualising the behaviour of complex systems
 Anton Burykin^{1},
 Madalena D Costa^{1, 2},
 Luca Citi^{1, 3} and
 Ary L Goldberger^{1, 2}Email author
DOI: 10.1186/14726947146
© Burykin et al.; licensee BioMed Central Ltd. 2014
Received: 23 September 2013
Accepted: 6 January 2014
Published: 18 January 2014
Abstract
Background
Physiologic signals, such as cardiac interbeat intervals, exhibit complex fluctuations. However, capturing important dynamical properties, including nonstationarities may not be feasible from conventional time series graphical representations.
Methods
We introduce a simpletoimplement visualisation method, termed dynamical density delay mapping (“D3Map” technique) that provides an animated representation of a system’s dynamics. The method is based on a generalization of conventional twodimensional (2D) Poincaré plots, which are scatter plots where each data point, x(n), in a time series is plotted against the adjacent one, x(n + 1). First, we divide the original time series, x(n) (n = 1,…, N), into a sequence of segments (windows). Next, for each segment, a threedimensional (3D) Poincaré surface plot of x(n), x(n + 1), h[x(n),x(n + 1)] is generated, in which the third dimension, h, represents the relative frequency of occurrence of each (x(n),x(n + 1)) point. This 3D Poincaré surface is then chromatised by mapping the relative frequency h values onto a colour scheme. We also generate a colourised 2D contour plot from each time series segment using the same colourmap scheme as for the 3D Poincaré surface. Finally, the original time series graph, the colourised 3D Poincaré surface plot, and its projection as a colourised 2D contour map for each segment, are animated to create the full “D3Map.”
Results
We first exemplify the D3Map method using the cardiac interbeat interval time series from a healthy subject during sleeping hours. The animations uncover complex dynamical changes, such as transitions between states, and the relative amount of time the system spends in each state. We also illustrate the utility of the method in detecting hidden temporal patterns in the heart rate dynamics of a patient with atrial fibrillation. The videos, as well as the source code, are made publicly available.
Conclusions
Animations based on density delay maps provide a new way of visualising dynamical properties of complex systems not apparent in time series graphs or standard Poincaré plot representations. Trainees in a variety of fields may find the animations useful as illustrations of fundamental but challenging concepts, such as nonstationarity and multistability. For investigators, the method may facilitate data exploration.
Keywords
Atrial fibrillation Delay map Heart rate variability Nonlinear dynamics Poincaré plot Sleep Time series VisualisationBackground
Physiologic and physical systems often generate highly complex output signals. An inherent feature of these signals is nonstationarity, defined as changes in their statistical properties (mean, variance, and higher moments) over time [1]. Furthermore, nonstationarities are important because they may contain information of both basic and translational interest. However, conventional representations, such as time series graphs, may not fully capture these timevarying properties. This challenge motivates the development of alternative ways to visualise the dynamics of complex time series.
We present a simple, new visualisation approach, based on the concept of delay maps, which highlights nonstationarities and related features. A classical delay map is the Poincaré plot^{a}[1], widely used in biomedicine by investigators probing heart rate variability [2]. Our dynamical density delay map method, termed D3Map, extends this concept to generate animated, colourised two and threedimensional representations.
This paper is intended to serve both as a theoretical introduction and a practical tutorial so that readers can recreate all the figures and animations using our data or their own. Therefore, all videos, data and source code are made publicly available (http://reylab.bidmc.harvard.edu/.D3Map/).
Although Poincaré plots are among the simplest representation of a system’s phase space [1], they may provide important information about the dynamics. For illustration purposes here we use heartbeat time series, i.e., a sequence of consecutive ventricular interbeat (QRS) intervals (termed RR) intervals obtained from an electrocardiographic (ECG) recording [2]. However, our method is applicable to any time series of sufficient length.
To highlight the added information that Poincaré plots provide, we shuffled (randomised) the order of the RR intervals. The histogram for this shuffled time series is unchanged from that of the original signal. However, the Poincaré plot is dramatically altered by the breakdown of correlations caused by the randomisation procedure (Figure 1, bottom).
However, even this socalled “3D Poincaré plot” omits information about the time evolution of a system’s dynamics. For example, consider a time series from a dynamical system with two states, A and B. The 3D Poincaré plot constructed from the entire time series will be the same whether the system spends the first half of the time in state A and the second half in state B, or whether the system changes states on a much shorter time scale. Thus, the 3D Poincaré plot for the entire time series does not capture the time sequence of state changes. This limitation is particularly important when probing the dynamics of physiologic systems, which are typically nonstationary. Here, we introduce a visualization technique, termed D3Map, based on a moving window analysis of the data.
Methods
 I.
Divide the initial time series into a sequence of either overlapping or nonoverlapping segments. Here, we analyse overlapping segments selected using a 500 sec moving window, shifted 10 sec at a time.
 II.
For each data segment, construct a normalised 2D histogram of RR(n) vs RR(n + 1) and its contour map. Let m be the number of counts in the highest bin of the histogram. The normalisation procedure consists of dividing the number of counts in each bin of the histogram by m. Here, we also opted for smoothing the histogram, using the dscatter^{d} Matlab function with the default parameter, which implements the method described in [8].
 III.
Colourise both the normalized 2D histogram and its contour map. In this implementation, we selected the jet colourmap,^{e} where blue represents the lowest and dark redbrown the highest density values, respectively. The colourised, normalised and smoothed 2D histogram of RR(n), RR(n + 1) is called the 3D surface Poincaré plot.
 IV.
In a single figure, present the 3D surface Poincaré plot, its colourised contour map^{f} and the graph of the original data segment used to generate them.
 V.
Sequentially present the figures with the analysis of each data segment as an animation [5]. Here, we used free video editing software VirtualDub (http://www.virtualdub.org) to create AVI movies and animated GIF images.^{g}
The cardiac interbeat interval time series used in the examples were from openaccess, deidentified data freely available at the NIHsponsored PhysioNet website (http://www.physionet.org) [9]. The links to the specific PhysioNet data sets used are indicated in parenthesis in the text below.
Results and discussion
The D3Map visualisation method was first applied to the cardiac interbeat interval time series from the healthy subject whose data (Additional file 1) are shown in Figures 1 and 2. All figures presented here for this subject were generated from the same data set (nsr001 at http://physionet.org/physiobank/database/nsr2db/), which we chose because it exhibits a type of nonstationarity characterised by relatively abrupt state transitions, which are commonly observed in the output of “freerunning” physiologic signals.
This interesting and anomalous pattern of heartbeat timings during AF illustrated in Figures 7 and 8 has been noted by other investigators [12]. However, its electrophysiologic mechanism (possibly involving dual AV nodal pathways or His bundle alternans) and its bedside (diagnostic and prognostic) implications have not been fully addressed. (Clinical details are not available in the two cases presented here). The D3Map visualization tool may help foster basic and clinical work on different dynamical subsets of AF, which in current practice are all grouped together.
Conclusions
We present a new animation method for visualising the temporal evolution of complex signal dynamics. For trainees, such simpletoconstruct plots (termed “D3Maps”) may be used to illustrate technical concepts, such as nonstationarity and multistability, which are difficult to grasp, but essential to understanding the dynamics of physiologic control systems. In addition, the animations may reveal unexpected patterns in data structure, which make the D3Maps useful in exploratory research, facilitating hypothesis generation and development and testing of mathematical/physiologic models.
Animations
All animations are available at http://reylab.bidmc.harvard.edu/.D3Map/.
Animation 1 (Sleep2D.avi) Movie of the 2D contour Poincare plots. Each frame shows both the colourised 2D contour plot and the 500 sec interbeat interval time series segment used to generate it. Movie frames were generated using "RR.dat" dataset (Additional file 1) and "Sleep2D.m" Matlab script (Additional file 2).
Animation 2 (Sleep3D.avi) D3Map animation. Each frame shows a 500 sec interbeat interval time series segment (bottom) and the corresponding colourised 3D surface Poincaré (top) and colourised 2D contour (middle) maps. Movie frames were generated using "RR.dat" dataset (Additional file 1) and "Sleep3D.m" Matlab script (Additional file 3). Note: The duration of the Sleep2D.avi and Sleep3D.avi movies is 3 min 20 sec, representing approximately 6 hours of heart rate dynamics obtained during sleep. A 500 sec moving window shifted 10 sec at a time was used to select the data segments displayed in each frame.
Animation 3 (Typical_AF_2D.avi) Movie of the 2D contour Poincare plots for a patient with a typical pattern of atrial fibrillation. Each frame shows both the colourised 2D contour plot and the 1200 sec interbeat interval time series segment used to generate it. Movie frames were generated using "Typical_AF.dat" dataset (Additional file 4) and "Typical_AF_2D.m" Matlab script (Additional file 5).
Animation 4 (Typical_AF_3D.avi) D3Map animation derived from the analysis of the RR interval time series of a patient with a typical pattern of atrial fibrillation. Each frame shows a 1200 sec interbeat interval time series segment (bottom) and the corresponding colourised 3D surface Poincaré (top) and colourised 2D contour (middle) maps. Movie frames were generated using "Typical_AF.dat" dataset (Additional file 4) and "Typical_AF_3D.m" Matlab script (Additional file 6).
Animation 5 (Atypical_AF_3D.avi) Movie of the 2D contour Poincare plots for a patient with an atypical pattern of atrial fibrillation. Each frame shows both the colourised 2D contour plot and the 1200 sec interbeat interval time series segment used to generate it. Movie frames were generated using "Atypical_AF.dat" dataset (Additional file 7) and "Atypical_AF_2D.m" Matlab script (Additional file 8).
Animation 6 (Atypical_AF_3D.avi) D3Map animation derived from the analysis of the RR interval time series of a patient with an atypical pattern of atrial fibrillation. Each frame shows a 1200 sec interbeat interval time series segment (bottom) and the corresponding colourised 3D surface Poincaré (top) and colourised 2D contour (middle) maps. Movie frames were generated using "Atypical_AF.dat" dataset (Additional file 7) and "Atypical_AF_3D.m" Matlab script (Additional file 9).
Note: The duration of the Typical_AF_2D.avi and Typical_AF_3D.avi movies is 56 sec, representing approximately 24 h of heart rate dynamics. The duration of the Atypical_AF_3D.avi, and Atypical_AF_3D.avi movies is 47 sec, representing approximately 20 h of heart rate dynamics. A 1200 sec moving window shifted 300 sec at a time was used to select the data segments displayed in each frame of the movies. All Matlab scripts listed above use dscatter2 Matlab function (Additional file 10) to calculate smoothed Poincare map.
Endnotes
^{a}The Poincaré plot is sometimes referred to as a Lorenz plot, a delay map, or a return map.
^{b}Matlab code for creating Poincaré plot from the sequence of RR intervals is listed below:
data = load('RR.dat’); % load the Additional file 1 with the RR time series
RR = data(:,2);
RRx = RR(1:end1);
RRy = RR(2:end);
^{c}Due to the normalisation, the height of the 2D histogram is not exactly equivalent to the probability density, but represents a relative occupancy of the phase space mapped onto the [0:1] interval.
^{d}Available at http://www.mathworks.com/matlabcentral/fileexchange/8430flowcytometrydatareaderandvisualisation.
^{e}A colourmap is defined as a specific chromatic sequence, which is linearly mapped onto the [0:1] interval, so any number between 0 and 1 corresponds to a particular colour.
^{f}Since we cannot see the “back” of the 3D surface Poincaré plot, we included its 2D contour Poincaré plot in each figure.
^{g}We opted to create animations using standalone specialised software, rather than Matlab see Appendix in ref. [13]. Specialised multimedia editing software has additional options for postprocessing, such as compression algorithms, annotation editing (e.g., adding extra frames with a movie title, explanatory text or subtitles), as well as audio (e.g., narrations or music).
Abbreviations
 D3Map:

Dynamical density delay map
 HR:

Heart rate.
Declarations
Acknowledgments
We gratefully acknowledge support from the Wyss Institute for Biologically Inspired Engineering, the G. Harold and Leila Y. Mathers Charitable Foundation, the James S. McDonnell Foundation and the National Institutes of Health (grants R00AG030677 and R01GM104987).
Authors’ Affiliations
References
 Kantz H, Schreiber T: Nonlinear Time Series Analysis. 1997, Cambridge: Cambridge University PressGoogle Scholar
 Camm AJ, Malik M, Bigger JT, Breithardt G, Cerutti S, et al: Heart rate variability. Standards of measurement, physiological interpretation, and clinical use. Eur Heart J. 1996, 17: 354381. 10.1093/oxfordjournals.eurheartj.a014868.View ArticleGoogle Scholar
 Brennan M, Palaniswami M, Kamen P: Do existing measures of Poincare plot geometry reflect nonlinear features of heart rate variability?. IEEE Trans Biomed Eng. 2001, 48: 13421347. 10.1109/10.959330.View ArticlePubMedGoogle Scholar
 Piskorski J, Guzik P: Filtering Poincaré plots. CMST: Comput Meth Sci Tech. 2005, 11: 3948.View ArticleGoogle Scholar
 Piskorski J, Guzik P: Dynamic decomposition of Poincaré plots for multivariate analysis and visualization of simultaneously recorded physiological time series. CMST: Comp Meth Sci Tech. 2010, 16: 181186.View ArticleGoogle Scholar
 Piskorski J, Guzik P: Geometry of the Poincare plot of RR intervals and its asymmetry in healthy adults. Physiol Meas. 2007, 28: 287300. 10.1088/09673334/28/3/005.View ArticlePubMedGoogle Scholar
 Piskorski J, Guzik P: Asymmetric properties of longterm and total heart rate variability. Med Biol Eng Comput. 2011, 49: 12891297. 10.1007/s115170110834z.View ArticlePubMedPubMed CentralGoogle Scholar
 Eilers PH, Goeman JJ: Enhancing scatterplots with smoothed densities. Bioinformatics. 2004, 20: 623628. 10.1093/bioinformatics/btg454.View ArticlePubMedGoogle Scholar
 Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PC, Mark RG, Mietus JE, Moody GB, Peng CK, Stanley HE: PhysioBank, PhysioToolkit, and PhysioNet: components of a new Research Resource for Complex Physiologic Signals. Circulation. 2000, 101: e215e220. 10.1161/01.CIR.101.23.e215.View ArticlePubMedGoogle Scholar
 Tsuboi K, Deguchi A, Hagiwara H: Relationship between heart rate variability using Lorenz plot and sleep level. Conf Proc IEEE Eng Med Biol Soc. 2010, 2010: 52945297.PubMedGoogle Scholar
 Gammaitoni L, Hanggi P, Jung P, Marchesoni F: Stochastic resonance. Rev Mod Phys. 1998, 70: 223287. 10.1103/RevModPhys.70.223.View ArticleGoogle Scholar
 Climent AM, de la Salud GM, Husser D, Castells F, Millet J, Bollmann A: Poincare surface profiles of RR intervals: a novel noninvasive method for the evaluation of preferential AV nodal conduction during atrial fibrillation. Biomed Eng IEEE Trans Biomed Eng. 2009, 56: 433442.View ArticlePubMedGoogle Scholar
 Burykin A, Peck T, Krejci V, Vannucci A, Kangrga I, Buchman TG: Toward optimal display of physiologic status in critical care: I. Recreating bedside displays from archived physiologic data. J Crit Care. 2011, 26 (105): e101e109.Google Scholar
 The prepublication history for this paper can be accessed here:http://www.biomedcentral.com/14726947/14/6/prepub
Prepublication history
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.