Recurrence quantification analysis across sleep stages

Biomedical Signal Processing and Control 20 (2015) 107–116

Contents lists available at ScienceDirect

Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc

Recurrence quantification analysis across sleep stages Jerome Rolink a,∗ , Martin Kutz a , Pedro Fonseca b,c , Xi Long b,c , Berno Misgeld a , Steffen Leonhardt a a b c

Philips Chair for Medical Information Technology (MedIT), RWTH Aachen University, Germany Department of Electrical Engineering, Eindhoven University of Technology, The Netherlands Personal Health Group, Philips Group Innovation Research, The Netherlands

a r t i c l e

i n f o

Article history: Received 3 September 2014 Received in revised form 26 January 2015 Accepted 9 April 2015 Available online 16 May 2015 Keywords: Recurrence quantification analysis Sleep stages Feature extraction Cardio-respiratory features

a b s t r a c t In this work we employ a nonlinear data analysis method called recurrence quantification analysis (RQA) to analyze differences between sleep stages and wake using cardio-respiratory signals, only. The data were recorded during full-night polysomnographies of 313 healthy subjects in nine different sleep laboratories. The raw signals are first normalized to common time bases and ranges. Thirteen different RQA and cross-RQA features derived from ECG, respiratory effort, heart rate and their combinations are additionally reconditioned with windowed standard deviation filters and ZSCORE normalization procedures leading to a total feature count of 195. The discriminative power between Wake, NREM and REM of each feature is evaluated using the Cohen’s kappa coefficient. Besides kappa performance, sensitivity, specificity, accuracy and inter-correlations of the best 20 features with high discriminative power is also analyzed. The best kappa values for each class versus the other classes are 0.24, 0.12 and 0.31 for NREM, REM and Wake, respectively. Significance is tested with ANOVA F-test (mostly p < 0.001). The results are compared to known cardio-respiratory features for sleep analysis. We conclude that many RQA features are suited to discriminate between Wake and Sleep, whereas the differentiation between REM and the other classes remains in the midrange. © 2015 Elsevier Ltd. All rights reserved.

1. Introduction In order to diagnose sleep or sleep related disorders, nowadays, polysomnography (PSG) is performed to record physiological data over a whole night in a sleep laboratory. Its key component is the measurement of cerebral activity with electroencephalography (EEG). Several electrodes have to be applied to the patients scalp at defined positions which can only be done by trained persons. In addition, electrocardiogram (ECG), respiratory activity, muscular activation, i.e. electromyogram (EMG), body temperature and eye movements, i.e. electrooculogram (EOG), are monitored. The amount of different parameters already evokes a high complexity of such recordings. With increasing age, chronic illnesses and sleep disorders occur more often. The demographic change will not allow each subject with potential sleep disorder to be analyzed in a sleep laboratory as this is too cost- and time-intensive. Alternatives

∗ Corresponding author at: Philips Chair for Medical Information Technology (MedIT), RWTH Aachen University, Pauwelsstrasse 20, 52074 Aachen, Germany. Tel.: +49 241 8023211; fax: +49 241 8082442. E-mail address: [email protected] (J. Rolink). http://dx.doi.org/10.1016/j.bspc.2015.04.006 1746-8094/© 2015 Elsevier Ltd. All rights reserved.

need to be established to explore sleep in more convenient environments, for example at home. Several studies show that it is possible to perform sleep staging to a certain extent using cardio-respiratory signals and body movements [4,6,25,30,31,41]. For example frequency and time-frequency analysis are applied on cardiac data to determine sleep or sleep-related disorders, such as sleep disordered breathing or sleep apnea [17,23]. The measurement of cardiac and respiratory signals can be transferred to the home environment more easily than EEG measurements, e.g. by using wearable sensor suites [21] or integrating the sensors into the bed mattress or pillow [9]. New methods and algorithms have to be developed to extract relevant data and multi-modal (coupled) signals using less sensors. Park et al. and Mita for example show improved algorithms to derive respiration from ECG measurements [36,38]. The aim of the research in this domain is to get an accurate estimate of sleep stages without the need of EOG, EMG or even EEG recordings. Recording in a familiar environment helps reducing the so-called first night effect biasing the regular sleep behavior. It also allows long-term measurements that are more representative than only one single night. Sleep is a very complex physiological process, not yet completely understood, in which the parasympathetic nervous system

108

J. Rolink et al. / Biomedical Signal Processing and Control 20 (2015) 107–116

(NS) predominantly controls the sleep depth by regulating the body temperature, heart rate and breathing activity. However, the sympathetic NS acts against it, so that the body is always alert to external stimuli [42]. It is obvious that sleep is essential and vital as sleep deprivation results in serious behavioral and psychological disorders that finally lead to obesity, cardiovascular morbidity, traffic accidents and death [1]. In this work we analyze cardio-respiratory signals across wake and sleep stages at night. We employ recurrence quantification analysis (RQA), which is known to be a powerful tool to study complex (physiological) systems. RQA has already successfully been applied to EEG data [44] to distinguish between sleep stages. Smietanoswki et al. determined several RQA features using recurrence plots on heart rate (HR) signals to derive dynamics of the heart rate variability (HRV) in healthy subjects and patients with obstructive sleep apnea (OSA) [43]. They found out that the change of recurrences between healthy and OSA subjects is different, arising from a more complex heart rate dynamic for the rapid eye movement (REM) stage. RQA is also an adequate method to study the nonlinear dynamic properties of QT and RR intervals during acute myocardial ischemia [39]. Terrill et al. apply RQA features, amongst others, solely on respiratory effort (RE) signals to classify 30 s epochs of infant sleep into Wake, REM and NREM stages [46]. They show that most of the computed RQA features are top-rated by a feature selection algorithm and contribute to better performance. Another promising approach for sleep staging is the exploration of phase synchrograms by means of RQA, as proposed by Nguyen et al. [37]. That study was also performed on data acquired from healthy infants and additionally on simulated data to show the improvements of the method compared to conventional methods to analyze synchrograms. Cardio-respiratory phase synchronization (position and number of heartbeats within respiratory cycles), which is often inspected with synchrograms, significantly changes with sleep stages, as discovered by Bartsch et al. [2]. However, to the knowledge of the authors, no work was published employing RQA on both, exclusively raw cardiac and respiratory data, for adult sleep stage investigation. Therefore, in this work, we concentrate on extracting thirteen of the most common RQA features from ECG, RE, HR signals and the two combinations of {ECG+RE} and {RE+HR} to determine the most discriminative features between different sleep stages and wake. This paper first gives a general overview of RQA as a quantifying description of recurrence plots (RPs), succeeded by a list of extracted RQA features. Then, it describes the employed data set and signal preprocessing steps for RQA, followed by a summary of the RQA embedding parameters. Finally, the ability of 195 computed cardio-respiratory RQA features derived from 313 whole night recordings of different subjects – representing more than 2300 h of data – to distinguish between NREM, REM and Wake stages is presented and discussed.

2. Materials and methods 2.1. Recurrence quantification analysis Recurrence quantification analysis (RQA) is a method for describing recurring states in the phase space of dynamic systems. The concept of recurrence was introduced by Henri Poincaré in 1890: the trajectory of a classical system will return infinitely many times to a limited region in phase space [15,40]. Often, only one time series xn (t) of a complex system x in d dimensions (x ∈ RD ) is observable. For this case, Takens proposed in 1981 a time delay embedding technique to reconstruct the phase space of the original system. Embedding m ≥ 2 · D + 1 (correlation dimension D) sampling points

of the discrete time series with delay  for each time ti , then xˆ is a reconstruction of the original trajectory: xˆ (ti ) =

m 

xn · (ti + (k − 1) · ) · ek ,

(1)

k=1

with ek as the unit vector of the dimension k. A diffeomorphism exists, so the embedding preserves the properties of the original attractor [34,45]. To visualize recurrences in multidimensional phase space, Eckmann et al. presented recurrence plots (RPs) in 1987 [11]. An RP is  The entry Ri,j is an illustration of the binary recurrence matrix R. set to 1 (black pixel in the plot) if two states at different times are similar, i.e. the points xˆ (ti ) and xˆ (tj ) in phase space fall inside a ball of radius ε: Ri,j (ε) = (ε − xˆ (ti ) − xˆ (tj ))

i, j = 1, . . ., N,

(2)

where  is the Heaviside step function and || · || is a norm, e.g. the Euclidean. RPs can be interpreted as a generalized form of autocorrelation. The main diagonal line represents the identity. Diagonal lines in general indicate deterministic behavior, regular patterns of parallel lines reveal a periodicity. Isolated points are typical for chaotic systems, whereas a rectangular patch results from an accumulation of subsequent states in a limited region of phase space during laminar sections of a process. For cross-recurrences of two time series xn and yn , embedded in the same phase space, (2) is modified for cross-recurrence plots (CRP) in an analog way [34,48,50]: CRi,j (ε) = (ε − xˆ (ti ) − yˆ (tj )).

(3)

Thus, it is possible to study interrelations between the two signals, e.g. synchronization, time shift or distortion. Unlike RPs, CRPs are usually not symmetric and therefore bowed lines may appear. Examples of RPs and CRPs applied to ECG, respiratory effort and heart rate samples are shown in Fig. 1. Their creation and interpretation will be discussed in the upcoming sections. In order to quantify the two dimensional plots for computerized analysis of complex dynamic systems, Zbilut and Webber offered in the 1990s a few features, that measure the percentage of recurring states REC and the amount of diagonal lines: DET describes the degree of determinism and ENT is the Shannon entropy concerning deterministic sections [47,49]. Additionally, Marwan et al. demonstrated in 2002, that the information contained in vertical structures (laminarity LAM, trapping time TT and maximum length Vmax ) is useful to detect chaos-chaos transitions. These RQA features were successfully applied to heart-rate variability data to predict life-threatening cardiac arrhythmia [35]. Besides, recurrence times and their entropy describe the periodicity [14,27]. In recent years, features from graph theory were added, which quantify the density of recurring states in phase space [3,10,33]. The Cross Recurrence Plot Toolbox for MATLAB® [32], which is used in this work, provides a way to compute windowed RQA. Especially for long data sets containing several hundreds of thousands data points and to achieve an acceptable resolution it is necessary to analyze small sub-RPs stepping through the signal. A total of thirteen different RQA features is computed. They are REC, DET, ENT, L, Lmax , LAM, TT, T(1), T(2), Vmax , Trec , clust and trans. Detailed descriptions and feature interpretations are specified in the following Section 2.2. On the one hand, particularly for physiological systems of complex dynamics an important advantage of RQA is that no assumptions on the time series are required, like linearity in Fourier transforms [50]. On the other hand, it is challenging to find optimum embedding parameters, because no strict rules exist. For an adequate embedding dimension to reconstruct the original

J. Rolink et al. / Biomedical Signal Processing and Control 20 (2015) 107–116

109

Fig. 1. Recurrence and cross-recurrence plots for 30 s Wake (W), REM (R) and non-REM (N) samples (from left to right) using ECG, RE and HR signals (vertically aligned time series from top to bottom), taken from the same subject and normalized as described in Section 2.4. The horizontally aligned time series all represent the same RE sample of the analyzed sample. The global Euclidean distance for each plot is color coded from dark to light.

trajectory we choose the lowest m for which the number of false nearest neighbors drops below a threshold [22]. The time delay  is set to meet the recommended first minimum of the mutual information for all subjects [13]. For the choice of threshold ε, Marwan et al. [34] lists several “rules of thumb”. To cover a broad range of the dynamics, we select a subject-dependent threshold so that REC (percentage of entries in the RP above a certain threshold ε, which is a measure for auto-/cross-correlation, see Section 2.2) is about 10% for five randomly selected sections of each analyzed signal. In the next sections we describe the RQA features, the data set and the preprocessing of the cardio-respiratory data. Furthermore, we present how the RQA parameters are determined.

as the histogram of the lengths l; lmin is the minimal length. The RP of a continuous sinus, for example, contains only diagonal lines: the trajectory is predictable. • Average (L) and maximum length (Lmax ) of diagonal lines

N

L=

l=l

Nmin

l=lmin

P(l)

,

(7)

N

l ) : Lmax = max({li }i=1

(8)

Corresponds to the mean and maximum time, for which the process is deterministic. Nl is the number of diagonal lines. • Entropy (ENT) ENT = −

2.2. RQA features

l · P(l)

N 

p(l) · ln (p(l)) ,

(9)

l=lmin

This section aims at describing and giving possible interpretations of all thirteen RQA features computed using the Cross-Recurrence Plot (CRP) for MATLAB® [32,34]. • Recurrence rate (REC)

(4)

i,j=1

Percentage of black pixels in a recurrence plot of size N × N. It is a measure for the auto- (RP) or cross-correlation (CRP). • Determinism (DET) l=l

Nmin

l · P(l)

R i,j=1 i,j

:

(5)

Ratio of diagonal line structures and recurrence point density with

i,j=1

v=v

N min

v · P(v)

v=1 v · P(v)

l−1 

(1 − Ri−1,j−1 ) · (1 − Ri+l,j+l )

Ri+k,j+k

k=0

P(v) =

(10)

N 

v−1 

(1 − Ri,j ) · (1 − Ri,j+v )

Ri,j+k

(11)

k=0

in analogy to DET: The likelihood that the system will maintain in certain states, e.g. in the “quiet” section between two R-peaks in an ECG. • Trapping Time (TT) and longest vertical line (Vmax )

N

TT = (6)

,

with vertical line length v, minimum vertical line length vmin and with

i,j=1

N

P(l) =

degree of complexity. For uncorrelated noise

with no diagonals as well as for deterministic systems with one diagonal, the entropy tends to zero. • Laminarity LAM LAM =

N

N 

P(l) : Nl

N

1  REC = 2 Ri,j : N

DET =

where p(l) =

v=v

N min

v · P(v)

v=vmin P(v)

v ). Vmax = max({vl }N l=1

,

(12) (13)

110

J. Rolink et al. / Biomedical Signal Processing and Control 20 (2015) 107–116

Corresponds to the mean and maximum time, in which laminar states will maintain. • Recurrence times of 1st and 2nd type, T(1) and T(2) (cf. [14]): (1)

Tj

= {ik − ik−1 | Rik ,j = 1},

(14)

(2) Tj

= {ik − ik−1 | Rik ,j = 1 ∧ ik − ik−1 > 1}, k ∈ N

(15)

are the mean times between consecutive points within the ε-ball (1) (2) around xj with (Tj ) and without counting sojourn points (Tj ), i.e. considering only the first point on each trajectory that enters the ball. The recurrence times of 1st and 2nd type, T(1) and T(2), (1) (2) are the averaged times of Tj and Tj for all j = 1, . . ., N. • Recurrence period density entropy Trec (cf. [27]):

 1 P(t) · ln (P(t)) , · ln (Tmax ) Tmax

Trec = −

(16)

t=1

T

max R(k) is the recurrence time probabilwhere P(T ) = R(T )/ k=1 ity density, R(T) the histogram of recurrence times and Tmax the maximum recurrence time (within the meaning of the first time a trajectory re-enters the ε-ball of a state xj ). The entropy provides information about the aperiodicity. • Clustering coefficient clust and transitivity trans [3,10,33]

clust =

N 

N

R j,k=1 i,j

i=1

N

trans =

· Rj,k · Rk,i

REC i R i,j,k=1 i,j

N

· Rj,k · Rk,i

R i,j,k=1 i,j

· Rk,i

,

,

(17)

(18)

N

R , are methods of graph theory, counting where REC i = j=1 i,j triangles and triples. clust specifies the mean amount of very close related states among the recurrent ones and trans characterizes how strong recurrent states are connected. 2.3. Data set One part of the data set used in this study comprises full-night PSG and simultaneously recorded actigraphy (Actiwatch, Philips Respironics) of 15 healthy subjects with a Pittsburgh Sleep Quality Index (PSQI) [5] of less than 6 and no evidence of respiratory or sleep diseases. The subjects had an average age of 31.0 (±10.4) year. Nine subjects were monitored in the Sleep Health Center, Boston, USA (Alice 5 PSG, Philips Respironics) and six subjects were monitored in the Philips Experience Lab, Eindhoven, the Netherlands (Vitaport 3 PSG, TEMEC). Sleep staging was performed by a sleep technician according to the guidelines of the American Academy of Sleep Medicine (AASM) [18], distinguishing between Wake, REM, N1, N2 and N3. The trials were approved by the ethics committee of the local organizations and all subjects signed a consent form. The other part of the data set consists of two-night recordings of 164 subjects (298 recordings in total). It was acquired in the context of the EU SIESTA project between 1997 and 2000 in seven different sleep laboratories all following the same measurement protocol [24]. The database was restricted to subjects considered to be healthy (e.g. no shift workers, no depression, usual bedtime before midnight), with a PSQI of 5 at most. Sleep stages were scored by trained sleep technicians in six classes (Wake, REM, S1, S2, S3, S4) according to the Rechtschaffen & Kales (R&K) guidelines [19]. In this work we employ AASM annotations as they are more recent than the R&K and were derived from the R&K guidelines [18]. Note, that we only analyze feature differences in the classes Wake, NREM and REM, where N1–N3 and S1–S4 are merged into a common class NREM. In total, RQA is performed on 313 healthy subjects

representing more than 2300 h of sleep data (average record duration of 7.39 h per night). Thirteen different RQA features are derived from the following three modalities and two signal combinations: 1) 2) 3) 4) 5)

RECG : RQA of the electrocardiogram (ECG) RRE : RQA of thoracic respiratory effort (RE) RHR : RQA of heart rate (HR), derived from ECG XECG,RE : cross-RQA of ECG and RE XRE,HR : cross-RQA of RE and HR

In addition to the unnormalized data a ZSCORE normalization (remove mean value and scale by overall standard deviation) and a standard deviation filter with a sliding window of 30 s are applied to each subject in order to remove subject-dependent properties. The latter works like a moving average filter, but instead of taking the mean value it computes the standard deviation of the window (windowed standard deviation filter). The subject-bysubject differences are either caused by the general physiology or the employed measurement technique itself. For example, the amplitude of the RE signal is affected by the positioning of the sensor whereas the baseline and the variations of the HR are very individual. The following symbols are used to indicate the transformation on the RQA feature ⇒ : no normalization, : ZSCORE normalization, ◦: windowed standard deviation. Afterwards, the resulting features are interpolated to the ground truth annotation time stamps, resulting in 30 s epoch time series aligned with the reference. This leads to a total count of 195 RQA features1 that are analyzed concerning their discriminative power between sleep stages and wake. The total numbers of 30 s epochs of each analyzed class are NNREM =178,667 (64.39%), NREM =45,695 (16.47%) and NWake =53,129 (19.15%). 2.4. Signal preprocessing To ensure that all signals are in a valid range for RQA, they need to be normalized and, in the case of cross-RQA, resampled to a common time base. The employed procedures are described in the following section. Examples of three different sleep stages – Wake, REM and NREM – after normalization are shown in Fig. 1. The structures of the (cross-)recurrence plots are different for each sleep stage making it possible to distinguish them. 2.4.1. ECG Raw ECG acquisition was performed with different sampling rates of 100, 200, 256 or 500 Hz, depending on the sleep laboratory. Since most of these sampling rates are too high for efficient RQA computations and need to be equal for all subjects the ECG was re-sampled to a common sampling rate of 150 Hz. In addition, baseline drift has been removed with a high-pass filter (0.8 Hz cutoff frequency) and an amplitude min/max-normalization to a range of ±10[a . u .]. Most of the ECG signal, except the QRS complex, is then located within a ±1[a . u .] range. 2.4.2. Respiratory effort The raw thoracic effort signal of the PSG has a sample rate of 10 Hz, which is appropriate for RQA computations. Due to body motion and position changes of the subject, the signal is likely to change over time. RQA is working best, when the signals are normalized to similar amplitude ranges for different subjects. By default, each window within the RQA procedure of the CRP Toolbox is normalized to the standard deviation [32]. Not only valid respiration cycles, but also motion artifacts causing very high amplitudes

1 (3 modalities+2 signal combinations)·(1 original feature data + 2 normalizations)· 13 features = 195.

J. Rolink et al. / Biomedical Signal Processing and Control 20 (2015) 107–116

111

would then be normalized in the same way. This makes it impossible to distinguish between valid and invalid respiration segments. Therefore the windowed normalization is deactivated and an alternative procedure is developed. We determine the fundamental frequency of the respiratory signal using the Hilbert transform and use an adaptive comb-like bandpass [8] to extract the real respiration signal envelope changing over time. Afterwards, the raw respiratory data is scaled with twice the envelope amplitude leading to a harmonized range of about ±0.5 [a.u.] of regular breathing. Motion artifacts are still present since they are generally much higher in amplitude than the respiratory signal and have a much higher frequency component that is not considered during the normalization procedure. An example of a normalized signal segment with corresponding RPs is presented in Fig. 2. At around 50 s, the signal amplitude suddenly drops, although a breathing component is visible. By applying RQA to the raw signal segment it is possible to detect this artifact but it is impossible to extract information of the respiration itself. The main concern is the correct choice of ε. A too low/high ε destroys the RP structure, which is the basis for the computation of most RQA features (see Section 2.2). The choice of ε of a normalized RP is easier since the artifact is almost annihilated in the RP even if it is still visible in the time domain.

the well-known Hamilton–Tompkins algorithm [16]. The resulting differences in time intervals are used to determine the instantaneous HR. To obtain a signal with equidistant samples, the HR is up-sampled by linear interpolation to a sample rate of 10 Hz, equal to the RE sampling rate. Generally, absolute HR and HRV are varying from subject to subject. To overcome the inter-individual differences the up-sampled HR is normalized by dividing the median HR over the whole night and subtracting 1 to obtain a zero-centered HR signal. The variation of the HR is about 10% of the median value, so we finally scaled the HR by 10:

2.4.3. Heart rate The instantaneous HR is determined by localizing R-peaks in the raw ECG. The detection is using an improved QRS detector based on

2.4.5. RQA parameters To obtain useful information from the signals, it is necessary to choose adequate RQA parameters as mentioned before: dimension

ˆ [a.u.] = HR





HR [bpm] −1 median(HR) [bpm]

· 10 [a.u.] .

(19)

2.4.4. Signal adaptation for cross-RQA The cross-RQA requires two signals with the same sample rate. Two different cross-RQA, between ECG-RE and RE-HR (see combinations 4 and 5 in listing in Section 2.3), are computed in this work. To avoid elimination of signal components that might be important for RQA, we up-sample the signal with the lower sample rate and equalize both, i.e. RE to match the 150 Hz of the ECG in the first case. HR and RE are already at 10 Hz and therefore do not need to be adapted.

Fig. 2. Respiratory effort (RE) signal normalization example (top: high-pass filtered/normalized data, bottom: RP of the segment). Left: high-pass filtered RE is shown with the detected amplitude envelope determined with an adaptive comb-filter. Right: normalized RE signal by scaling raw signal with amplitude envelope. The sudden drop at 50 s is clearly visible in the RP and makes it impossible to extract RQA features correctly. Color coding of the RP ranges from dark to light, respectively the minimum ε and maximum ε of the corresponding segment.

112

J. Rolink et al. / Biomedical Signal Processing and Control 20 (2015) 107–116

Table 1 RQA parameters for each modality (see Section 2.3). m is the dimension,  the time delay,  the detection threshold, wsize the window size and wstep the step size. [smp] indicates number of samples at the sampling frequency Fs of the modality. Modality

Fs [Hz]

m

 [smp]

 (mean±SD)

RECG RRE RHR XECG,RE XRE,HR

150 10 10 150 10

4 4 4 4 4

9 10 20 2 10

0.355 0.480 0.377 1.702 0.925

m, time delay , threshold , window size wsize and window step size wstep (see Table 1). • Time delay : It can be estimated by computing the mutual information (MI). The optimal value for  is the occurrence of the first minimum of the MI. Afterwards, the embedding parameter m can be determined. For periodic signals  is in the order of T/4, with T the period length of the signal. In the case of RECG , the period length TECG is in the order of the length of a PQRST complex in the ECG (i.e. TECG ≈ 200–300 ms). Since respiration and HR changes are slow physiological processes. Their  values are higher than for ECG, i.e. average TRE ≈ 4 s and THR ≈ 8 s are used. When computing MI, it becomes apparent that for cross-RQA the values for  should be smaller than or equal to the smallest  of a single modality to achieve a good compromise for embedding both signals in the same phase space [34]. • Dimension m: On the one hand, high values for dimension are recommended to include the dynamics [50] and to prevent false recurrences [35]. On the other hand, each additional dimension increases the computational complexity. After applying the false nearest neighbors algorithm with the selected  to estimate the embedding parameter m, we find that m = 4 is sufficient for all signals and subjects. • Window size wsize : Since generating an RP for a continuous fullnight signal exceeds the limits of memory and CPU resources on most personal computers, RQA was done on sliding windows of size wsize . To ensure that enough data points are taken into account for the computations the windows are of different sizes for each modality: for ECG it is chosen to contain about five heart beats, for RE at least one respiratory cycle and for HR about thirty consecutive heart beats. For the crossrecurrence XECG,RE and XRE,HR it is also ensured that at least five consecutive heart beats and one respiratory cycle are within the processing window. This suggests that the time stamps of the windowed analysis have to be shifted backwards by wsize − (m − 1) ·  samples to match the original time stamps [32]. The corresponding window sizes for each modality are listed in Table 1. • Progression step of the sliding window wstep : For ECG, RE and HR signals wstep is chosen to guarantee an overlap of one heart beat for RECG /RHR and one respiratory cycle for RRE . Adequate wstep values for XECG,RE and XRE,HR were derived from RECG and RRE in the same way, respectively. • ε: As mentioned above, ε is determined subject-by-subject before the extraction of all RQA features. The aim is to obtain an ε so that the RQA feature REC, which is computed on five randomly chosen sections for each RECG , RRE , RHR , XECG,RE and XRE,HR , is at about 10% [34]. It ensures a wide dynamic range of the RQA features. 3. Results and discussion The mean (±SD) ε for all 313 subjects is shown in Table 1. The ε values are stable for RECG , RRE , RHR and XRE,HR , showing that the applied normalization procedure on the raw signals is valid. In

± ± ± ± ±

0.267 0.224 0.207 2.039 0.299

wsize [smp]

wstep [smp]

750 150 300 750 150

525 58 20 630 58

contrast, it is not possible to find similar ε values for XECG,RE that vary from subject to subject. One explanation is the difference in dynamics of the underlying modalities. Where RE is a slowly changing signal the ECG has high dynamic properties. This leads to very high ε values during inspiration and/or expiration, see top right XECG,RE CRPs of Fig. 1. The threshold ε therefore has to be increased in some subjects to obtain the desired 10% recurrence rate. The selected rule of thumb described by Marwan et al. [34] to determine ε might not be very well suited in this case. To evaluate the discriminative power of each computed feature between the classes Wake, NREM and REM we need a scoring that indicates the separability of the classes. We chose a metric originally developed to quantify the agreement between two independent observers, the Cohen’s Kappa coefficient of agreement : p0 = pc = =

TP + TN Nt

(20)

(TP + FP) · (TP + FN) + (FP + TN) · (FN + TN) Nt2 p0 − pc , 1 − pc

(21) (22)

using the True Positive (TP), True Negative (TN), False Positive (FP), False Negative (FN) and Total Number (Nt ) values of the confusion matrix at a certain decision threshold th where the positive class is either NREM, REM or Wake. By varying this parameter from the minimum to maximum feature value, not only performance curves (e.g. Receiver Operating Characteristic (ROC) or Precision Recall (PR) curves) are obtained but also values for accuracy (acc), sensitivity (sens) or specificity (spec) of the decision. The optimal th is then set at the maximal  value. The  metric is directly interpretable as the proportion of joint judgments in which there is agreement, after chance agreement is excluded [7] regardless of the possible imbalance between the classes. Due to the highly unbalanced class data of this work (see Section 2.3) the  score is much more adequate than the absolute standardized mean difference (ASMD) or the Mahalanobis distance for feature discrimination analysis. Landis and Koch describe a subjective interpretation of this coefficient, ranging from 0 = poor to 1 = (almost) perfect agreement [26]. Afterwards, the  values of the features are sorted in a top 20 and top 10 list for each class (NREM, REM and Wake) and for each modality (RECG , RRE , RHR , XECG,RE and XRE,HR ), respectively (see Table 2). In total 277,491 epochs of 30 s have been analyzed, representing 2312 h of data. NREM class covers 178,667 epochs (≈64% of the whole analyzed data) whereas REM stage makes up the smallest fraction with 45,695 epochs (≈16%). The percentage of Wake epochs is about 19% (53,129 epochs). Note, that Wake vs. all is the same as performing a Wake vs. Sleep analysis. The best features to discriminate the NREM class from others are obtained with the RECG and RRE RQA modalities (see Table 2). The ZSCORE () or windowed standard deviation (◦) normalization

J. Rolink et al. / Biomedical Signal Processing and Control 20 (2015) 107–116

113

Table 2  values for top ten features by modality for classes ‘NREM’/‘REM’/‘Wake’ vs. all. Overall top 20 features for each class are marked in bold text. Underlined features indicate that they are present within at least two classes of the top ten for each modality and at least within a top 20 list. Number of 30 s epochs: NNREM =178,667 (64.39%), NREM =45,695 (16.47%) and NWake =53,129 (19.15%). #

RECG

RRE

RHR

XECG,RE

XRE,HR

NREM vs. all 1 2 3 4 5 6 7 8 9 10

◦T (1) (0.18) ◦LAM (0.17)  Vmax (0.17) REC (0.17) T (1) (0.17)  TT (0.16)  LAM (0.16) ◦DET (0.16) REC (0.16)  ENTR (0.16)

◦clust (0.24) ◦T (2) (0.24) ◦ LAM (0.22) ◦trans (0.21) ◦Trec (0.21) ◦ T(1) (0.21) ◦ENTR (0.20) ◦ DET (0.19) T (2) (0.19) REC (0.18)

T (1) (0.16) L (0.16)  T(1) (0.15) ENTR (0.13)  L (0.13)  trans (0.13) trans (0.12) DET (0.11) ENTR (0.11)  Trec (0.10)

◦ENTR (0.16) ◦ T(1) (0.15) ◦clust (0.14) ◦ REC (0.14) ◦DET (0.12) ◦LAM (0.12) ◦ T(2) (0.11) ◦ Lmax (0.11) ◦ Vmax (0.10)  T(2) (0.10)

◦TT (0.14) ◦ T(2) (0.13) ◦ ENTR (0.13) ◦ Trec (0.13) ◦ L (0.13) ◦ T(1) (0.13)  Lmax (0.13)  Vmax (0.13)  Vmax (0.12)  Lmax (0.12)

REM vs. all 1 2 3 4 5 6 7 8 9 10

◦ clust (0.07)  clust (0.05) TT (0.05)  L (0.05)  clust (0.05) ◦ Vmax (0.03)  Lmax (0.03) ◦DET (0.03) ◦LAM (0.03) REC (0.03)

 L (0.13)  Lmax (0.13) ◦clust (0.12) T (1) (0.12) T (1) (0.11)  DET (0.11)  LAM (0.11)  ENTR (0.10) REC (0.10)  trans (0.10)

 clust (0.11) ENTR (0.10)  TT (0.10) trans (0.09) L (0.09)  clust (0.09) T (1) (0.09) DET (0.08) ENTR (0.08)  TT (0.08)

◦clust (0.09)  T(2) (0.05) ◦ T(1) (0.04)  clust (0.04) ◦ REC (0.03) ◦ENTR (0.03) ◦ TT (0.03) ◦ Lmax (0.03)  Trec (0.03) ◦ T(2) (0.03)

 T(2) (0.07)  ENTR (0.07)  L (0.07)  TT (0.07)  DET (0.06)  clust (0.06)  LAM (0.06)  Lmax (0.06)  L (0.06)  ENTR (0.06)

Wake vs. all 1 2 3 4 5 6 7 8 9 10

◦LAM (0.31) ◦DET (0.30) ◦T (1) (0.28) ◦ trans (0.25) ◦ ENTR (0.25) REC (0.24) T (1) (0.23)  Vmax (0.22) REC (0.22) TT (0.22)

◦T (2) (0.23) ◦Trec (0.22) ◦ Vmax (0.22) T (2) (0.21) ◦ TT (0.19) T (1) (0.19) ◦ENTR (0.19)  T(2) (0.19) ◦trans (0.18) T (1) (0.18)

 T(1) (0.12) L (0.11) T (1) (0.10)  L (0.10)  Trec (0.10)  Trec (0.10) ◦ L (0.09)  T(2) (0.09)  T(2) (0.09)  Lmax (0.08)

◦LAM (0.22) ◦DET (0.22) ◦ENTR (0.19) ◦ T(1) (0.18) ◦ REC (0.17)  T(2) (0.14)  TT (0.14)  DET (0.13)  LAM (0.13) ◦ T(2) (0.13)

 Vmax (0.13)  Lmax (0.13) ◦ REC (0.13)  REC (0.12) ◦ Trec (0.12) ◦ Lmax (0.11) ◦ ENTR (0.11)  L (0.11) ◦ T(2) (0.11) ◦ T(1) (0.11)

: no normalization, : ZSCORE normalization, ◦: windowed standard deviation

is most frequently represented within the top 20 features. During NREM sleep, respiration and heart activity are very regular compared to other classes. This physiological behavior is expressed by RQA features, such as DET, LAM, T(1), T(2) and ENTR, which can be found within the top features for this class (max() = 0.24). The amount of REM data and the derived feature values are mostly comparable to Wake (see boxplots in Fig. 4). This leads to overall low  values (max() = 0.13) for all modalities in this class. Nevertheless, the best RQA modalities for REM are RRE and RHR . Within the top 20 the ZSCORE normalization procedure () is most represented for good  values. In contrast, the  values, especially for the RECG modality, are much higher for Wake discrimination than for the other two classes (max() = 0.31). The windowed standard deviation normalization (◦) is represented most frequently among the best features. During Wake the regularity seen in respiration and heart activity is not present anymore. Also, more movements occur causing artifacts [12] in the signals that are misinterpreted by the RQA feature extraction algorithms. Note, the cross-RQA modalities XECG,RE and XRE,HR only have moderate  values over all classes, meaning that RQA features derived of coupled signals do not contribute in distinguishing between the Wake, NREM and REM. To reduce the total feature amount for the following analyses, we decide to select the features according to those criteria: 1) the feature is at least represented once in the top 20 of one class for all modalities

2) the feature is at least represented two times in the top ten (for each modality) 3) the  score is greater than 0.18 (slight to fair agreement [26]). This ensures that the selected features discriminate best between at least two or even all classes and we obtain reduced subsets of relevant features for each modality (in total 20 different features): 1) 2) 3) 4) 5) 6) 7) 8) 9) 10) 11) 12) 13) 14) 15) 16) 17) 18)

RECG  REC RECG  TT RECG  T(1) RECG  REC RECG ◦ DET RECG ◦ LAM RECG ◦ T(1) RRE  T(1) RRE  REC RRE  T(1) RRE  T(2) RRE ◦ ENTR RRE ◦ T(2) RRE ◦ Trec RRE ◦ clust RRE ◦ trans RHR  ENTR XECG,RE ◦ DET

114

J. Rolink et al. / Biomedical Signal Processing and Control 20 (2015) 107–116

Fig. 3. Performance measures , sens, spec and acc for the classes NREM (N), REM (R) and Wake (W), respectively versus the remaining classes for features a–d that are known from literature (see Table 3) and for the final RQA feature subset (1–20) as described in Section 3.

Table 3 Features known from literature applied on the same data set described in Section 2.3. #

Feature

Description

Reference

a

RE v9

Redmond et al. [41]

b

RRmean

c d

ECG˛ REDTW

RE frequency variation over 9 epochs Normalized mean ECG RR-interval ECG ˛ coefficient of the DFA RE dynamic time warp distance

Redmond et al. [41] Kantelhardt et al. [20] Long et al. [28]

19) XECG,RE ◦ ENTR 20) XECG,RE ◦ LAM with : no normalization, : ZSCORE normalization, ◦: windowed standard deviation. Note that the numbering of this list is only used for the plotting of the individual feature performances (see Fig. 3) and feature inter-correlation (see Fig. 4) and does not give any ranking or interpretation of the features at this point. The features have been grouped by RQA modality for better overview and are also underlined in Table 2. We compute the performance values , sens, spec and acc of this subset and of four additional existing features known to work well for sleep staging (see description in Table 3 [20,28,41]). Note, that the data basis is the same as for the RQA derived features. The existing features are labeled with a, b, c and d. The performance results of all the prior features are shown in

Fig. 3. Feature a (RE v9 ,  = 0.37, sens = 0.62, spec = 0.79, acc = 0.68) and feature c (ECG˛ ,  = 0.34, sens = 0.68, spec = 0.68, acc = 0.68) demonstrate remarkable performance characteristics in discriminating between NREM and the other classes. Feature d (REDTW ,  = 0.31, sens = 0.32, spec = 0.94, acc = 0.82) has good performance identifying Wake. The performance during REM sleep is moderate with the best feature a (RE v9 ,  = 0.16, sens = 0.83, spec = 0.49, acc = 0.54). Nevertheless, the RQA features – the best are shown below the existing features in Fig. 3 – are very effective in discriminating the Wake class, too. Especially, the ECG related feature 6 (RECG ◦ LAM,  = 0.31, sens = 0.36, spec = 0.92, acc = 0.81) is outperforming the literature RE feature d in terms of sensitivity at almost equal specificity and accuracy. Regarding the definite class imbalance between Wake and the other classes, each increase in sensitivity has a much stronger impact on  than the remaining performance metrics. It also shows that ECG features are at least as powerful as RE features for Sleep/Wake discrimination. For REM detection, the described features remain in the midrange. The best feature for this class is number 15 (RRE ◦ clust,  = 0.12, sens = 0.58, spec = 0.61, acc = 0.61), having higher spec and acc than feature a at the cost of lower sens. However, the ability to distinguish NREM is outperformed by almost all non-RQA features. Nevertheless, the best is feature 15 again (RRE ◦ clust,  = 0.24, sens = 0.62, spec = 0.64, acc = 0.63). This feature also has remarkable sens for Wake:  = 0.13, sens = 0.65, spec = 0.59, acc = 0.58. To get a visual impression of the data distribution of the 20 selected features, boxplots for the classes NREM (N), REM (R) and Wake (W) are shown in Fig. 4 at the top. Each box illustrates the median (blue line) and quartiles of the data (upper and lower bound of the box). As the  values already demonstrate in Fig. 3, the ability to discriminate between Sleep and Wake is clearly visible in many features. It is also noticeable that REM is often very similar either to NREM or Wake. The significances of the class differences for each feature are tested with an ANOVA F-test computed with MATLAB® (The Mathworks Inc., Massachusetts, USA). Their corresponding p-values are shown with red symbols above the boxplots of each feature (* : p < 0.001, † : p < 0.01, ‡ : p < 0.05, ♦ : p < 0.1, ◦ : p ≥ 0.1). Some of the RQA features have zero-inflated feature distributions, which means that the distribution is peaking at one edge, e.g. 0, and has a very long tail, e.g. towards high values. Examples can be found looking at features 5, 6, 18 or 20 (RECG ◦ DET, RECG ◦ LAM, XECG,RE ◦ DET or XECG,RE ◦ LAM) of the extracted subset. These features might have much higher discriminative power in combination with other features [29] than the  value suggests. Inter-feature correlation is an important point to consider as well. Most classifiers need uncorrelated data to prevent over-fitting. In addition to the boxplots, the Pearson correlation of all features is computed and illustrated in the lower part of Fig. 4. The diagonal line means perfect correlation (=1) since this is the auto-correlation. The lighter the color coding the higher the feature correlation. Some of the extracted features are strongly correlated. Features 5 and 6 show one example, i.e. RECG ◦ DET and RECG ◦ LAM. Generally, features derived from the same modality have stronger correlations, which can be seen by the light colored regions for features 5–6, 8, 10, 11 and 18–20. The contrast to other modalities are evident in feature 2 (RECG  TT) and feature 17 (RHR  ENTR). Most of the features presented in Fig. 4 have significance values for the class discrimination of p < 0.001 or p < 0.01. Only features 7, 18 and 20 (RECG ◦ T(1), XECG,RE ◦ DET and XECG,RE ◦ LAM) are not very significant in discriminating between NREM and REM, as well features 2 and 17 (RECG  TT and RHR  ENTR) between REM and Wake discrimination.

J. Rolink et al. / Biomedical Signal Processing and Control 20 (2015) 107–116

115

Fig. 4. Top: boxplots with median and quartiles for NREM, REM and Wake class of the selected subset consisting of 20 features. Significance of the class differences was tested using an ANOVA F-test (* : p < 0.001, † : p < 0.01, ‡ : p < 0.05, ♦ : p < 0.1, ◦ : p ≥ 0.1 in red above the classes for each feature). Bottom: Inter-feature Pearson correlation matrix, grouped by RQA modalities (RECG , RRE , RHR , XECG,RE , XRE,HR ). : no normalization, : ZSCORE normalization, ◦: windowed standard deviation

4. Conclusions In this work we apply RQA on ECG, RE, HR and the two combinations of {ECG+RE} and {RE+HR} to derive features for sleep stage detection from the results. In total we extract 195 features from 313 full-night PSG records. A couple of these are suited to distinguish between NREM, REM and/or Wake classes. The significance of class discrimination is p < 0.001 for almost all of the selected features. Especially for Wake/Sleep distinction, we obtain  values that indicate very acceptable discrimination. The introduced normalization procedure for the RE signal aims at reducing movement artifacts and regularizing the amplitude over time and different subjects. Another approach to detect Wake might be to analyze the differences between the normalized and

non-normalized signals, since movements are likely to occur during the Wake stages. Also during REM the raw RE signal contains many artifacts compared to NREM [12]. Since we concentrate on the RQA features in this work, no further investigations were done on this matter. To conclude, although not all selected RQA features correlate they can be used together to enhance the differentiation between the classes. Note, that the features are not meant to replace existing features from literature. They can rather be used in addition to other cardio-respiratory features with classifiers to optimize sleep staging tasks since RQA aims at extracting the dynamic, often very non-linear, characteristics of data without the need for temporal or spectral a-priori knowledge or assumptions (see Section 2.1). Using two modalities to distinguish Sleep and Wake gives the opportunity

116

J. Rolink et al. / Biomedical Signal Processing and Control 20 (2015) 107–116

to either enhance redundancy in case of corrupted measurement signals or to increase detection performance using additional information. Acknowledgements This work was funded by Philips Research. The authors would like to thank all reviewers for their perceptive comments. References [1] S. Banks, D.F. Dinges, Behavioral and physiological consequences of sleep restriction, J. Clin. Sleep Med. 3 (5) (2007) 519–528. [2] R.P. Bartsch, A.Y. Schumann, J.W. Kantelhardt, T. Penzel, P.C. Ivanov, Phase transitions in physiologic coupling, Proc. Natl. Acad. Sci. U. S. A. 109 (26) (2012) 10,181–10,186, http://dx.doi.org/10.1073/pnas.1204568109 [3] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, D. Hwang, Complex networks: structure and dynamics, Phys. Rep. 424 (4) (2006) 175–308, http://dx.doi.org/ 10.1016/j.physrep.2005.10.009 ´ J. Salinger, S. Nevsímalová, Spectral analysis of [4] P. Busek, J. Vanková, J. Opavsky, the heart rate variability in sleep, Physiol. Res. 54 (4) (2005) 369–376. [5] D. Buysse, C. Reynolds, T. Monk, S. Berman, The Pittsburgh Sleep Quality Index: a new instrument for psychiatric practice and research, Psychiatry Res. 28 (1989) 193–213. [6] C. Canisius, T. Ploch, T. Penzel, A. Jerrentrup, K. Kesper, Classifying sleep stages using the heart rate spectrum – comparison and feasibility in healthy subjects and sleep apnea patients, in: Proceedings of Biosignal, Berlin, Germany, 2010, pp. 16–19. [7] J. Cohen, A coefficient of agreement for nominal scales, Educ. Psychol. Meas. 20 (1) (1960) 37–46, http://dx.doi.org/10.1177/001316446002000104 [8] D. Cyrill, J. McNames, M. Aboy, Adaptive comb filters for quasiperiodic physiologic signals, in: Proceedings of the 25th Annual International Conference of the IEEE, vol. 3, IEEE, Cancun, Mexico, 2003, pp. 2439–2442. [9] S. Devot, A.M. Bianchi, E. Naujokat, M. Mendez, A. Brauers, S. Cerutti, Sleep monitoring through a textile recording system, IEEE Eng. Med. Biol. Soc. 2007 (2007) 2560–2563, http://dx.doi.org/10.1109/IEMBS.2007.4352851 [10] R. Donner, J. Heitzig, J. Donges, Y. Zou, N. Marwan, J. Kurths, The geometry of chaotic dynamics – a complex network perspective, Eur. Phys. J. 84 (4) (2011) 653–672. [11] J.P. Eckmann, S.O. Kamphorst, D. Ruelle, Recurrence plots of dynamical systems, Europhys. Lett. 973 (4) (1987) 973–977. [12] J. Foussier, X. Long, P. Fonseca, B. Misgeld, S. Leonhardt, On the relationship of arousals and artifacts in respiratory effort signals, in: Y.T. Zhang (Ed.), in: IFMBE Proceedings, ICHI 2013, vol. 42, Springer International Publishing, Vilamoura, Portugal, 2014, pp. 31–34, http://dx.doi.org/10.1007/978-3-319-03005-0 [13] A. Fraser, H. Swinney, Independent coordinates for strange attractors from mutual information, Phys. Rev. A 33 (2) (1986) 1134–1140. [14] J. Gao, Recurrence time statistics for chaotic systems and their applications, Phys. Rev. Lett. 83 (16) (1999) 3178–3181, http://dx.doi.org/10.1103/ PhysRevLett.83.3178 [15] P. Gujrati, Poincaré Recurrence, Zermelo’s Second Law Paradox, and Probabilistic Origin in Statistical Mechanics, 2008, pp. 1–5, arXiv:0803.0983. [16] P.S. Hamilton, W.J. Tompkins, Quantitative investigation of QRS detection rules using the MIT/BIH arrhythmia database, IEEE Trans. Biomed. Eng. 33 (12) (1986) 1157–1165. [17] M.F. Hilton, R.A. Bates, K.R. Godfrey, M.J. Chappell, R.M. Cayton, Evaluation of frequency and time-frequency spectral analysis of heart rate variability as a diagnostic marker of the sleep apnoea syndrome, Med. Biol. Eng. Comput. 37 (6) (1999) 760–769. [18] C. Iber, S. Ancoli-Israel, A.L. Chesson, S.F. Quan, The AASM Manual for the Scoring of Sleep and Associated Events: Rules, Terminology and Technical Specifications, American Academy of Sleep Medicine (AASM), Westchester, IL, 2007. [19] A. Kales, A. Rechtschaffen, in: A. Rechtschaffen, A. Kales (Eds.), A manual of standardized terminology, techniques and scoring system for sleep stages of human subjects, U. S. National Institute of Neurological Diseases and Blindness, Neurological Information Network, Bethesda, MD, 1968, no 204. [20] J. Kantelhardt, Y. Ashkenazy, P. Ivanov, A. Bunde, S. Havlin, T. Penzel, J.H. Peter, H. Stanley, Characterization of sleep stages by correlations in the magnitude and sign of heartbeat increments, Phys. Rev. E 65 (5) (2002) 1–6, http://dx.doi. org/10.1103/PhysRevE.65.051908 [21] W. Karlen, C. Mattiussi, D. Floreano, Sleep and wake classification with ECG and respiratory effort signals, IEEE Trans. Biomed. Circuits Syst. 3 (2) (2009) 71–78, http://dx.doi.org/10.1109/TBCAS.2008.2008817 [22] M. Kennel, R. Brown, H. Abarbanel, Determining embedding dimension for phase-space reconstruction using a geometrical construction, Phys. Rev. A 45 (6) (1992) 3403–3411. [23] K. Kesper, S. Canisius, T. Penzel, T. Ploch, W. Cassel, ECG signal analysis for the assessment of sleep-disordered breathing and sleep pattern, Med. Biol. Eng. Comput. 50 (2) (2012) 135–144, http://dx.doi.org/10.1007/s11517-0110853-9

[24] G. Klösch, B. Kemp, T. Penzel, A. Schlögl, P. Rappelsberger, E. Trenker, G. Gruber, J. Zeitlhofer, B. Saletu, W. Herrmann, S.L. Himanen, D. Kunz, M.J. Barbanoj, J. Röschke, A. Värri, G. Dorffner, The SIESTA project polygraphic and clinical database, IEEE Eng. Med. Biol. Mag. 20 (3) (2001) 51–57, http://dx.doi.org/10. 1109/51.932725 [25] Y. Kurihara, K. Watanabe, Sleep-stage decision algorithm by using heartbeat and body-movement signals, IEEE Trans. Syst. Man Cybern. 42 (6) (2012) 1450–1459. [26] J.R. Landis, G.G. Koch, The measurement of observer agreement for categorical data, Biometrics 33 (1) (1977) 159–174. [27] M.A. Little, P.E. McSharry, S.J. Roberts, D.A.E. Costello, I.M. Moroz, Exploiting nonlinear recurrence and fractal scaling properties for voice disorder detection, Biomed. Eng. Online 6 (23) (2007), http://dx.doi.org/10.1186/1475-925X-623 [28] X. Long, P. Fonseca, J. Foussier, R. Haakma, R. Aarts, Using dynamic time warping for sleep and wake discrimination, in: IEEE Engineering in Medicine and Biology Society, vol. 25, Hong Kong/Shenzhen, China, 2012, pp. 886–889. [29] X. Long, P. Fonseca, J. Foussier, R. Haakma, R. Aarts, Sleep and wake classification with actigraphy and respiratory effort using dynamic warping, IEEE J. Biomed. Health Inform. (2013) 1–12, http://dx.doi.org/10.1109/JBHI.2013.2284610 [30] X. Long, P. Fonseca, R. Haakma, R.M. Aarts, J. Foussier, Time-frequency analysis of heart rate variability for sleep and wake classification, in: IEEE Conference on Bioinformatics and Bioengineering (BIBE), November, IEEE, Larnaca, Cyprus, 2012, pp. 85–90. [31] X. Long, J. Foussier, P. Fonseca, R. Haakma, R.M. Aarts, Analyzing respiratory effort amplitude for automated sleep stage classification, Biomed. Signal Process. Control 14 (2014) 197–205, http://dx.doi.org/10.1016/j.bspc.2014.08.001 [32] N. Marwan, Cross Recurrence Plot Toolbox for Matlab, Version 5.15, Release 28.6, 2010, URL: www.recurrence-plot.tk [33] N. Marwan, J.F. Donges, Y. Zou, R.V. Donner, J. Kurths, Complex network approach for recurrence analysis of time series, Phys. Lett. A 373 (46) (2009) 4246–4254, http://dx.doi.org/10.1016/j.physleta.2009.09.042 [34] N. Marwan, M.C. Romano, M. Thiel, J. Kurths, Recurrence plots for the analysis of complex systems, Phys. Rep. 438 (5-6) (2007) 237–329, http://dx.doi.org/10. 1016/j.physrep.2006.11.001 [35] N. Marwan, N. Wessel, U. Meyerfeldt, A. Schirdewan, J. Kurths, Recurrence-plotbased measures of complexity and their application to heart-rate-variability data, Phys. Rev. E 66 (2) (2002) 026,702, http://dx.doi.org/10.1103/PhysRevE. 66.026702 [36] M. Mita, Algorithm for the classification of multi-modulating signals on the electrocardiogram, Med. Biol. Eng. Comput. 45 (3) (2007) 241–250, http://dx. doi.org/10.1007/s11517-006-0130-5 [37] C.D. Nguyen, S.J. Wilson, S. Crozier, Automated quantification of the synchrogram by recurrence plot analysis, IEEE Trans. Biomed. Eng. 59 (4) (2012) 946–955, http://dx.doi.org/10.1109/TBME.2011.2179937 [38] S.B. Park, Y.S. Noh, S.J. Park, H.R. Yoon, An improved algorithm for respiration signal extraction from electrocardiogram measured by conductive textile electrodes using instantaneous frequency estimation, Med. Biol. Eng. Comput. 46 (2) (2008) 147–158, http://dx.doi.org/10.1007/s11517-007-0302-y [39] Y. Peng, Z. Sun, Characterization of QT and RR interval series during acute myocardial ischemia by means of recurrence quantification analysis, Med. Biol. Eng. Comput. 49 (1) (2011) 25–31, http://dx.doi.org/10.1007/s11517-0100671-5 [40] H. Poincaré, Sur le problème des trois corps et les équations de la dynamique, Acta Math. 13 (1) (1890) A3–A270. [41] S.J. Redmond, P. de Chazal, C. O’Brien, S. Ryan, W.T. McNicholas, C. Heneghan, Sleep staging using cardiorespiratory signals, Somnologie 11 (4) (2007) 245–256, http://dx.doi.org/10.1007/s11818-007-0314-8 [42] J. Siegel, Why we sleep, Sci. Am. 289 (5) (2003) 92–97. [43] M. Smietanowski, W. Szelenberger, A. Trzebski, Nonlinear dynamics of the cardiovascular parameters in sleep and sleep apnea, J. Physiol. Pharmacol. 57 (11) (2006) 55–68. [44] I.H. Song, D.S. Lee, S.I. Kim, Recurrence quantification analysis of sleep electoencephalogram in sleep apnea syndrome in humans, Neurosci. Lett. 366 (2) (2004) 148–153, http://dx.doi.org/10.1016/j.neulet.2004.05.025 [45] F. Takens, Detecting strange attractors in turbulence, in: Dynamical Systems and Turbulence, Warwick 1980, Springer, 1981, pp. 366–381. [46] P.I. Terrill, S.J. Wilson, S. Suresh, D.M. Cooper, C. Dakin, Application of recurrence quantification analysis to automatically estimate infant sleep states using a single channel of respiratory data, Med. Biol. Eng. Comput. 50 (8) (2012) 851–865, http://dx.doi.org/10.1007/s11517-012-0918-4 [47] C. Webber, J. Zbilut, Dynamical assessment of physiological systems and states using recurrence plot strategies, J. Appl. Physiol. 76 (1994) 965–973. [48] J. Zbilut, A. Giuliani, C. Webber, Detecting deterministic signals in exceptionally noisy environments using cross-recurrence quantification, Phys. Lett. A 246 (1998) 122–128. [49] J. Zbilut, C. Webber, Embeddings and delays as derived from quantification of recurrence plots, Phys. Lett. A 171 (1992) 199–203. [50] J.P. Zbilut, C.L. Webber, Recurrence quantification analysis Wiley Encyclopedia of Biomedical Engineering, vol. 2, John Wiley & Sons, Inc., 2006, http://dx.doi. org/10.1002/9780471740360.ebs1355.