 Technical advance
 Open Access
 Published:
Dynamical density delay maps: simple, new method for visualising the behaviour of complex systems
BMC Medical Informatics and Decision Making volume 14, Article number: 6 (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.
Background
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.
A Poincaré plot (map) of a time series is the scatter plot of x(n) versus x(n + 1) or equivalently, x(n1) versus x(n), hence the general term delay map ^{b}[1–7]. In Figure 1 (top), we show the RR interval time series from a healthy subject during sleeping hours, along with its histogram and the standard (“black and white”) Poincaré plot. As described further below, this dataset was selected because it exhibits a common type of nonstationarity characterized by relatively abrupt changes of “state”. The “cometlike” shape of the Poincaré plot of the RR interval time series (Figure 1, top) indicates a statistical relationship between consecutive data points where a given RR interval is likely followed or preceded by one of similar duration. In fact, the spread of the data points in the direction perpendicular to the diagonal of the Poincaré plot is a measure of the time series’ local variance (“short term variability”).
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).
This example illustrates how the standard Poincaré plot resolves one limitation of the histogram, namely that the latter does not provide information about correlations among the data points. However, standard (“black and white”) Poincaré plots have a salient limitation of their own; namely, they do not capture information about the density of the data points. This limitation can be overcome by modifying the standard Poincaré plot to include the relative frequency of pairs of consecutive data points by constructing a 2D histogram^{c} of RR(n), RR(n + 1) (Figure 2).
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
The D3Map approach consists of five simple steps:

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.
In this specific case, (Figure 3, colourised 2D contour plot), there are two states (termed S1 and S2), centered around RR intervals of 0.65 sec and 1 sec (corresponding to periods of faster and slower heart rates of about 92 beats/min and 60 beats/min, respectively). These two states likely correspond to different sleep stages [10]: the faster rates to periods of rapid eye movement (REM) sleep (or transient awakenings) and the slower rates to deeper sleep. A movie of the 2D contour plots (Animation 1, Sleep2D.avi) providing information about the nature (rapid or gradual) of the transitions between the two states is included in the supplemental material. From a physiologic point of view, the “flipping” between these two states represents a highly nonstationary behavior, reminiscent of bistability[11].
Figure 4 shows “snapshots” of the D3Map animation included in the online supplementary material (Animation 2, Sleep3D.avi) for four different data segments. These snapshots demonstrate system transitions between different states over time. The D3Map animations show that the states, themselves, are not static. Instead they appear to consist of multiple substates, manifesting as local undulations in the 3D surface, which continuously emerge and disappear. This finding is not apparent from visual inspection of the raw interbeat interval time series.
To further explore the utility of this visualisation method, we also show the analysis of cardiac interbeat interval time series from two patients with chronic atrial fibrillation, a common cardiac arrhythmia that may be associated with stroke or heart failure. The RR dynamics during atrial fibrillation are presumed to be random. In fact, clinicians use the expression “irregularly irregular” to describe the unpredictable timings of the heartbeat during this arrhythmia. This assumption is supported by the circular shape of the contour map in Figure 5 from the first of these subjects (dataset a1nn from http://www.physionet.org/challenge/chaos/), which is typical of time series of uncorrelated values. The D3Map movie (snapshot shown in Figure 6) shows that in this case, the dynamics is relatively stationary.
In contrast, for the second patient with chronic atrial fibrillation (dataset a2nn from http://www.physionet.org/challenge/chaos/), the contour map (Figure 7) shows four loci of high density consistent with a more structured pattern of ventricular response, which are not readily discerned from the time series. The structure centered around (0.4, 0.4) sec in the contour map indicates that the RR intervals fluctuate around 0.4 sec for at least three consecutive beats. Similarly, the structure centered around (0.8, 0.8) sec indicates that the RR intervals fluctuate around 0.8 sec for at least three consecutive beats. The changes from 0.4 to 0.8 and from 0.8 to 0.4, respectively, generate the structures centered around (0.4, 0.8) and (0.8, 0.4) sec. The D3Map (snapshot shown in Figure 8) animation confirms the dynamical nature of these fluctuation patterns.
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.
References
 1.
Kantz H, Schreiber T: Nonlinear Time Series Analysis. 1997, Cambridge: Cambridge University Press
 2.
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.
 3.
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.
 4.
Piskorski J, Guzik P: Filtering Poincaré plots. CMST: Comput Meth Sci Tech. 2005, 11: 3948.
 5.
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.
 6.
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.
 7.
Piskorski J, Guzik P: Asymmetric properties of longterm and total heart rate variability. Med Biol Eng Comput. 2011, 49: 12891297. 10.1007/s115170110834z.
 8.
Eilers PH, Goeman JJ: Enhancing scatterplots with smoothed densities. Bioinformatics. 2004, 20: 623628. 10.1093/bioinformatics/btg454.
 9.
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.
 10.
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.
 11.
Gammaitoni L, Hanggi P, Jung P, Marchesoni F: Stochastic resonance. Rev Mod Phys. 1998, 70: 223287. 10.1103/RevModPhys.70.223.
 12.
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.
 13.
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.
Prepublication history
The prepublication history for this paper can be accessed here:http://www.biomedcentral.com/14726947/14/6/prepub
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).
Author information
Affiliations
Corresponding author
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
The conceptualization and initial implementation of the dynamical density delay map method were developed by AB, MDC and ALG. All four authors contributed to the refinement of the method, analysis of the findings and to the composition and editing of the manuscript. All authors have read and approved this manuscript.
Electronic supplementary material
(dscatter2.m).
Additional file 10: Matlab function that creates a scatter plot coloured by density. This is a modified version of the original dscatter function http://www.mathworks.com/matlabcentral/fileexchange/8430flowcytometrydatareaderandvisualization/content/dscatter.m. All Matlab scripts listed above use this function. (M 7 KB)
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
About this article
Cite this article
Burykin, A., Costa, M.D., Citi, L. et al. Dynamical density delay maps: simple, new method for visualising the behaviour of complex systems. BMC Med Inform Decis Mak 14, 6 (2014). https://doi.org/10.1186/14726947146
Received:
Accepted:
Published:
Keywords
 Atrial fibrillation
 Delay map
 Heart rate variability
 Nonlinear dynamics
 Poincaré plot
 Sleep
 Time series
 Visualisation